Complex Dynamical Systems,

Sponsored by the U.S. National Science Foundation, a workshop on the boundary element method (BEM) was held on the campus of the University of Akron during September 1–3, 2010 (NSF, 2010, “Workshop on the Emerging Applications and Future Directions of the Boundary Element Method,” University of Akron, Ohio, September 1–3). This paper was prepared after this workshop by the organizers and participants based on the presentations and discussions at the workshop. The paper aims to review the major research achievements in the last decade, the current status, and the future directions of the BEM in the next decade. The review starts with a brief introduction to the BEM. Then, new developments in Green's functions, symmetric Galerkin formulations, boundary meshfree methods, and variationally based BEM formulations are reviewed. Next, fast solution methods for efficiently solving the BEM systems of equations, namely, the fast multipole method, the pre-corrected fast Fourier transformation method, and the adaptive cross approximation method are presented. Emerging applications of the BEM in solving microelectromechanical systems, composites, functionally graded materials, fracture mechanics, acoustic, elastic and electromagnetic waves, time-domain problems, and coupled methods are reviewed. Finally, future directions of the BEM as envisioned by the authors for the next five to ten years are discussed. This paper is intended for students, researchers, and engineers who are new in BEM research and wish to have an overview of the field. Technical details of the BEM and related approaches discussed in the review can be found in the Reference section with more than 400 papers cited in this review.


Introduction
The boundary element method (BEM) is a numerical method for solving boundary-value or initial-value problems formulated in boundary integral equations (BIEs).In some literature, it is also called the boundary integral equation method because of this relationship.The most important advantage of the BEM as compared with other domain-based numerical methods is the reduction of the dimensions of the problems to be solved, from 3D to 2D, or from 2D to 1D, leading to much easier mesh generation as compared with domain-based methods.The second advantage is the accuracy the BEM offers, due to the nature of integrals used in the formulations.The third advantage of the BEM is its capability of accurately modeling infinite domain problems without introducing additional, often cumbersome, conditions at infinity.Table 1 is a comparison of the BEM with the finite element method (FEM), which may change in the future based on new developments in either one or both methods.
1.1 A Brief History.The direct BIE formulations (in which the unknown functions have clear physical meanings) and their modern numerical solutions with boundary elements originated more than 40 years ago during the 1960s.The 2D potential problem was first formulated in terms of a direct BIE and solved numerically by Jaswon, et al. [1][2][3].This work was later extended to the 2D elastostatic case by Rizzo in the early 1960s as his dissertation work at the University of Illinois, which was later published as a journal article in 1967 [4].Following these early works, extensive research efforts were carried out in the BIE formulations of many problems in applied mechanics and in the numerical solutions during the 1960s and 1970s [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20].The name boundary element method appeared in the middle of the 1970s in an attempt to make an analogy with the FEM [21][22][23].
Some personal accounts of the historical developments of the BIE and BEM in various places of the U.S., Europe, and China were given by Rizzo [24], Cruse [25], Watson [26], Shippy [27], Mukherjee [28], Telles [29], and Yao and Du [30] in a special issue of the Electronic Journal of Boundary Elements in 2003 [31].A comprehensive review of the heritage and early history of the BIE and BEM was given by Cheng and Cheng in 2005 [32].This review provides vivid descriptions of the long history and major contributions to the mathematical foundations of the integral methods by the pioneers from the 18th to the first half of the 20th centuries and the early days of the BEM until the late of 1970s.Another comprehensive review of the BEM research up to the early of 1980s was given by Tanaka in 1983 [33].A special review was given by Tanaka, Sladek and Sladek in 1994 [34] on the research efforts in the late 1980s and early 1990s regarding the various approaches to dealing with the singular and hypersingular integrals in the BIEs.Ten representative books on the BEM, published during the last three decades, can be found in Refs.[35][36][37][38][39][40][41][42][43][44] which can be consulted to further study the details of the BEM.
1.2 A Review of the BEM Formulation.We use the 3D potential theory problem as an example to show the key steps and main results in the BEM formulation.For a domain V with boundary S (Fig. 1), we have the following boundary-value problem: where / is the potential field.The boundary conditions (BCs) are /ðxÞ ¼ /ðxÞ; x 2 S / ðDirichlet BCÞ (2) qðxÞ @/ @n ðxÞ ¼ qðxÞ; x 2 S q ðNeumann BCÞ (3) in which the overbar indicates the prescribed value, S / [ S q ¼ S, and n is the outward normal of the boundary S (Fig. 1).
The fundamental solution (also known as full-space Green's function) Gðx; yÞ for 3D potential problems is given by where r is the distance between source point x and field point y.
The normal derivative of G is Fðx; yÞ @Gðx; yÞ @nðyÞ ¼ À 1 4pr 2 r; k n k ðyÞ (5) with r; k ¼ @r=@y k ¼ ðy k À x k Þ=r.The fundamental solution Gðx; yÞ represents the response (potential) at y due to a unit source at x (See Sec.2.1 for the definition of a Green's function).
If we apply the second Green's identity with / and Gðx; yÞ in domain V with boundary S, we obtain the following representation integral: /ðxÞ ¼ ð S Gðx; yÞqðyÞ À Fðx; yÞ/ðyÞ ½ dSðyÞ; x 2 V (6) Equation ( 6) is an expression of the solution / inside the domain V for Eq.(1).Once the boundary values of both / and q are known on the entire boundary S, Eq. ( 6) can be applied to calculate / everywhere in V, if needed.To solve the unknown boundary values of / and q on S, we let x tend to the boundary S to obtain the following BIE: cðxÞ/ðxÞ ¼ ð S Gðx; yÞqðyÞ À Fðx; yÞ/ðyÞ ½ dSðyÞ; x 2 S (7) in which cðxÞ is a coefficient and cðxÞ ¼ 1=2 if S is smooth around x. Equation ( 7) is called a conventional or singular BIE for solving potential problems.
If we take the normal derivative of Eq. ( 6) and let the source point x go to boundary S, we obtain the so-called hypersingular BIE:   where the two new kernels for 3D are Kðx; yÞ @Gðx; yÞ @nðxÞ ¼ 1 4pr 2 r; k n k ðxÞ (9) Hðx; yÞ @Fðx; yÞ @nðxÞ ¼ 1 4pr 3 n k ðxÞn k ðyÞ À 3r; k n k ðxÞr; l n l ðyÞ ½ (10) Hypersingular BIEs are needed in the BEM when it is applied to model cracks, thin shapes and exterior acoustic wave problems, as discussed in other sections.The Fðx; yÞ and Kðx; yÞ kernels are strongly singular.Therefore, the second integral in BIE (Eq.( 7)) and the first integral in BIE (Eq.( 8)) are Cauchy-principal value (CPV) integrals in general.The Hðx; yÞ kernel is hypersingular, and the second integral in BIE (Eq.( 8)) is a Hadamard-finite part (HFP) integral.There are different approaches [34,[45][46][47] to evaluate or avoid the singular and hypersingular integrals in the BIE formulations which are no longer the barriers in the implementation of the BEM.
The BIE as shown in Eq. ( 8) is called a direct BIE formulation in which the density functions / and q have clear physical meanings.The so-called indirect BIE formulations are also available in the literature, in which only one integral is applied and the density function has no clear physical meanings.In some cases, the indirect BIE formulation is advantageous because of the reduction of the integrals with which to work.More details about the indirect BIE formulations can be found in many textbooks on the BEM.
Complex variable BIE formulations in the BEM are also popular for dealing with various 2D problems (see, e.g., Refs.[48][49][50][51]; [44] (Chapter 4 on 2D elasticity)), because of the efficiencies in using the complex functions.However, this approach has limited applicability as to solving 3D problems.
As an example to see how to discretize the BIEs, assume that we use N constant boundary elements (surface patches for 3D, line segments for 2D) on boundary S. We can obtain the following discretized equation of BIE (7):  q 1 q 2 . . . with ð DSj F i dS; for i; j ¼ 1; 2; :::; N in which G i and F i are the kernels with the source point x placed at node i, DS j represents element j, and d ij is the Kronecker d symbol.
In the conventional BEM approach, we form a standard linear system of equations by applying the boundary condition at each node and switching the columns in the two matrices in Eq. ( 11) to obtain: where A is the coefficient matrix, k is the unknown vector, and b is the known vector.The disadvantages of the conventional BEM formulations are: (1) The construction of matrix A requires OðN 2 Þ operations using Eq. ( 12), (2) The size of the memory required for storing A is also OðN 2 Þ because A is in general a nonsymmetrical and dense matrix, and (3) The solution of the system in Eq. ( 13) using direct solvers such as Gauss elimination requires OðN 3 Þ operations.Thus, for several decades, the conventional BEM approach of solving Eq. ( 13) directly has been limited to small-scale models with, at most, a few thousands of elements, for example, on a desktop PC running a 32-bit operating system.In Secs. 2 and 3, new BIE formulations and new fast solution methods aimed to address these shortcomings of the conventional BEM approach are discussed.
The rest of this review paper is organized as follows: In Sec. 2, new developments in the last decade in the theoretical formulations related to the BIE and BEM are discussed.These developments are: new development in the Green's functions which are the key ingredients in the BEM, the symmetric Galerkin BEM formulation, the BIE related boundary meshfree methods, and the variationally based BEM formulations.In Sec. 3, three fast solution methods that emerged in the last decade for solving the BEM systems of equations are reviewed.These fast methods are the fast multipole method, the pre-corrected fast Fourier transformation method, and the adaptive cross approximation method.These fast solution methods are the key factors in the last decade that have contributed to revitalizing the BEM.In Sec. 4, modern applications of the BEM in solving largescale, critical research and industrial problems using the advanced BEM are presented.The areas of applications include: modeling microelectromechanical systems, composite materials, functionally graded materials, fracture mechanics, acoustic, elastic and electromagnetic waves, time-domain problems, and coupling of the BEM with other methods.In Sec. 5, we discuss the future directions of the BEM research, development and education for the next five to ten years, based on our own opinions.More than 400 references that highlight the history and advances of the BEM research in the last decade or so are listed at the end of this paper for readers to further study the related topics.
2 New Formulations 2.1 Green's Functions for Anisotropic and Layered Materials.Green's function (GF) is named after George Green, who wrote the famous essay on the application of mathematical analysis to the theories of electricity and magnetism [52].In his essay, Green coined the term "potential" to denote the results obtained by adding the masses of all the particles of a system, each divided by its distance from a given point.The general properties of the potential function were subsequently developed and applied to electricity and magnetism.The formula connecting the surface and volume integrals, now known as the Green's theorem, was also introduced in that work.A Green's function, also called singular function or fundamental solution (a full-space GF), represents the solution in a given system due to a point source.For example, the GF denoted by Gðx; yÞ for potential theory problems governed by Eq. (1) and in an infinite domain satisfies r 2 Gðx; yÞ þ dðx; yÞ ¼ 0; x; y 2 R 2 =R 3 (14) where r 2 ¼ @ 2 ð Þ=@y i @y i , and R 2 and R 3 indicate the full 2D and 3D spaces, respectively.The Dirac d function dðx; yÞ in Eq. ( 14) represents a unit source (e.g., heat source) at the source point x, and Gðx; yÞ represents the response (e.g., temperature) at the field point y due to the source.The solution of Eq. ( 14) in 3D is given in Eq. (4).GFs play a key role in the BEM.Without the GFs, the boundary nature of the BEM will be hard to obtain.On the other hand, due to the use of the GF, the BEM is often regarded as a semi-analytical method and thus more accurate, because the BEM solution is built upon the analytical result provided by the GF.More information about the GFs can be found in Refs.[53,54].The GF concept is now extensively used in the solution of partial differential equations.The recently re-emerged method of fundamental solutions also makes direct use of the GFs [55,56].In the Eshelby-based micromechanics theory, the GFs further plays an important role [57].
The GF also played an important role in the development of various advanced material systems in the last decade.Anisotropy and multiphase coupling are important features in advanced smart material systems.The static GFs in 3D transversely isotropic piezoelectric infinite, half, and bimaterial spaces were derived by Ding's group at Zhejiang University [58,59] and Dunn's group at University of Colorado at Boulder [60] where the potential function method was employed.Furthermore, under the assumption of transverse isotropy, Ding's group has further derived the corresponding magnetoelectroelastic GFs in the infinite space, half space, and bimaterial space [61].
The corresponding 2D GFs in general anisotropy (including also the piezoelectric and magnetoelectroelastic coupled systems) can be best represented by the Stroh formalism based on the complex variable method for infinite, half, and bimaterial planes [62][63][64].In general, GFs of the degenerated case where repeated eigenvalues appear cannot be directly reduced from these expressions.However, those of slightly perturbed cases (i.e., the quasiisotropic) can give very accurate solutions [65].On the other hand, Yin [66] employed the Lekhnitskii's formalism and derived the Green's functions for infinite, half, and bimaterial planes where both general anisotropy (with distinguished eigenvalues) and quasi-isotropic case (with repeated eigenvalues in different ranks) were cast into unified expressions.
Three-dimensional GFs in an anisotropic and elastic full space were derived using various mathematical methods, including eigenvalues/vectors, Fourier transform, Radon transform, among others.A comprehensive review on this topic can be found in Lavagnino [67].Very recently, Buroni et al. [68] derived the set of solutions which contain any kind of mathematical degeneracy of the eigenvalues.The corresponding half-space and bimaterial GFs were obtained in Ref. [69].Similar GFs were also derived for the half-space case with general surface conditions [70] and for the bimaterial case with general (or imperfect) interface conditions [71].Simple line-integral expressions for the interfacial GFs and trimaterial GFs can also be obtained.
Three-dimensional GFs in the corresponding anisotropic piezoelectric and even magnetoelectroelastic full, half, and bimaterial space can be derived by following a similar procedure.Actually, to derive these GFs, one only needs to extend the problem dimension from three to four for the piezoelectric case and to five for the magnetoelectroelastic case.For instance, for the 3D magnetoelectroelastic case, the GFs can be found in Refs.[72,73].Besides the point "force" GFs, some point "dislocation" GFs were also derived, which have been shown to be very efficient in dealing with fracture problems in piezoelectric and magnetoelectroelastic solids [74,75].
Layered structures are very common in advanced composite structures, and therefore, derivation of the GFs in these systems is also important.Typical methods to find the GFs is to first apply the double Fourier transform or Hankel transform in the horizontal layer plane and solve the problem in the transformed domain via the propagator matrix method.The physical domain solution is obtained by applying the inverse transform, which usually requires numerical integration.Besides the GFs for elastic isotropic, transversely isotropic, thermoelastic, and poroelastic layered spaces, GFs for the piezoelectric layered half space and even piezoelectric functionally graded half space were also derived [76,77].Similar approaches have also been developed to derive the GFs in general anisotropic and layered elastic spaces [78].
Significant efforts have been conducted by different authors to derive dynamic GFs in anisotropic solids and further extending these ideas to the piezoelectric and magnetoelectroelastic cases.Denda et al. [79] derived the 2D time-harmonic GFs by means of the Radon transform, and implemented them for the eigenvalue analysis of piezoelectric solids.Following the same approach, Wang and Zhang [80] presented their time-domain, Laplace's domain, and 3D counterparts of the GFs.Wu and Chen [81] later derived explicit dynamic GFs for the 2D case.With respect to the magnetoelectroelastic case, Ren and Liu [82] presented a timeharmonic dynamic fundamental solution for transversely isotropic magnetoelectroelastic media under anti-plane deformation.Chen et al. [83] derived the explicit expressions for the dynamic potentials of an inclusion embedded in a "quasi-plane" magnetoelectroelastic medium of transverse isotropy and for the dynamic GFs of such a medium both in the time domain and in the frequency domain.Rojas-Dı ´az et al. [84] later derived the 2D and 3D timeharmonic GFs by means of Radon transform.
For the development of GFs in other fields, readers may consult a special issue of the journal Engineering Analysis with Boundary Elements published in 2005 [78] where some recent research on GFs and related BEM issues are discussed.

Symmetric
Galerkin BEM.There are two basic procedures that are generally used to reduce the continuous boundary integral Eqs. ( 7) and ( 8) to finite systems.The simpler procedure is collocation, wherein the BIEs are explicitly enforced at a finite set of points.In its simplest form, these collocation points are chosen to be the nodes used to discretize the boundary.If the boundary potential and flux are interpolated from their values at these N points, then the pointwise enforcement of Eqs. ( 7) and ( 8) provides the N equations needed to solve for the unknown values.Collocation necessarily leads to non-symmetric matrices system.
In contrast to collocation, the Galerkin approach does not require that the integral equations be satisfied at any point.Instead, the equations are enforced in a weighted average sense as follows: where w k are the chosen weight functions.The needed N equations can be generated by an appropriate choice of N weights.The strict definition of Galerkin is that the weight function w k is composed of the shape functions that are non-zero at the node x, the shape functions being the local basis functions used to interpolate the boundary functions.For a linear interpolation, the Galerkin weight functions are illustrated in Fig. 2.
In mathematical terminology, collocation is a strong solution, the equations being satisfied at the specified points, whereas Galerkin is a weak solution.The requirement that the equations hold in an integrated sense has a geometric interpretation: The approximate Galerkin solution is obtained by projecting the exact solution onto the subspace consisting of all functions which are a linear combination of the shape functions.The Galerkin solution is therefore the linear combination that is the closest to the exact solution.In general, the Galerkin method is more accurate than the collocation method, and also provides a more elegant treatment of boundary corners.However, the primary advantage of the Galerkin method is that the treatment of hypersingular integrals is actually much simpler than with collocation.The Galerkin form of the hypersingular integral exists if the interpolation of the boundary potential is continuous C 0 .For collocation, a smoother interpolation (differentiable C 1 ) is required for the existence of the integral.Standard elements, e.g., linear or quadratic, are continuous across element boundaries, and it is quite a bit more complicated to implement a differentiable interpolation.
2.2.1 Symmetric Galerkin Formulation.With Galerkin, x and y are treated equally, and thus it is possible to produce a symmetric coefficient matrix, the symmetric Galerkin boundary element method (SGBEM).This algorithm is based upon the symmetry properties of the Green's function, Gðx; yÞ ¼ Gðy; xÞ @Gðx; yÞ @nðyÞ ¼ À @Gðx; yÞ @nðxÞ ¼ @Gðy; xÞ @nðxÞ @ 2 Gðx; yÞ @nðxÞ@nðyÞ ¼ @ 2 Gðy; xÞ @nðxÞ@nðyÞ (18) It follows that if the potential is specified everywhere on the boundary, then the Galerkin form of the potential equation will yield a symmetric matrix; a similar statement holds for flux boundary data and the flux equation.For a general mixed boundary value problem, symmetry results if each equation is applied on its appropriate surface.That is, the potential equation is applied on the Dirichlet portion of the boundary and the flux equation is applied on the Neumann portion of the boundary.After discretization, the set of Eq. ( 17) can be written in blockmatrix form as f ½ / ¼ h ½ q, and in block-matrix these equations become The symmetry of the coefficient matrix, 21 , now follows from the properties of the kernel functions.The advantages of Galerkin, of course, come at a price: although, as shown in Fig. 2, the weight functions have local support, the additional boundary integration with respect to x is nevertheless computationally expensive.These costs can be reduced somewhat by exploiting symmetry or by fast solution methods, but Galerkin BEM will nevertheless require more computation time than collocation.

Brief
History.Almost all of the early BEM computations employed a collocation approximation.In 1977, Bui [85] produced a Galerkin formulation for fracture analysis, and in 1979, Sirtori [86] proposed the symmetric Galerkin formulation for linear elasticity.It was recognized almost after a decade that collocation BEM failed to retain many of the nice properties of finite elements: sign-definite symmetric matrix operators, a variational formulation, and proofs of convergence and stability.In addition, a collocation implementation for fracture analysis had to resort to higher order C 1 elements, e.g., Overhauser or Hermite elements to adequately deal with the hypersingular kernel, and these elements proved difficult to handle, especially in three dimensions.These negative aspects of collocation fueled the interest in the symmetric Galerkin boundary element formulation, and subsequent development was primarily carried out in Europe.

Recent Applications.
One of the important application areas of symmetric Galerkin is fracture.An implementation for fracture analysis in plane orthotropic elasticity was reported by Gray and Paulino [103].A weak form integral equation for 3D fracture analysis was proposed by Li, Mear, and Xiao [104].Frangi [105,106] developed a similar crack analysis based on integration by parts, and applied it to crack propagation and fatigue crack growth.Salvadori [107] presented a crack formulation for cohesive interface problems, and Sutradhar and Paulino [108] presented a symmetric Galerkin formulation to evaluate T-stress and stress intensity factors (SIFs) using the interaction integral method.Xu, Lie and Cen [109] presented a 2D crack propagation analysis using quasi-higher order elements.Phan et al. [110][111][112][113][114] presented SIF analysis with frictional contact sliding at discontinuity and junctions, SIF calculations for crack-inclusion interaction problems, and a SGBEM based quasi-static crack-growth prediction tool to investigate crack-particle interactions.Galerkin methods have found a broad range of engineering applications.For instance, an SGBEM formulation for 2D steady and incompressible flow was developed by Capuani et al. [115].A Stokes problem with general boundary condition with slip condition was reported by Reidinger and Steinbach [116].Other works using SGBEM include; for example, steady state harmonic solution of the Navier equation by Perez-Gavilan and Aliabadi [117], interface and multi-zone problems by Gray and Paulino [118], lower bound shakedown analysis by Zhang et al. [119], dynamic soilstructure interaction by Lehman and Antes [120], problems of finite elasticity with hyperelastic and incompressible materials by Polizotto [121], dynamic frequency domain viscoelastic problems subjected to steady-state time-harmonic loads by Perez-Gavilan and Aliabadi [122,123], analysis of Kirchhoff elastic plates by Frangi and Bonnet [124].

Meshfree Methods
Based on BIEs.Boundary-based meshfree methods are first introduced in this section.This is followed by some details of the boundary node method (BNM) and the extended boundary node method (EBNM).
The primary idea in meshfree methods is to divorce the traditional coupling between spatial discretization (meshing) and interpolation of the unknown solution, as is commonly practiced in the FEM or the BEM.Instead, a "diffuse" approximation is used to represent the unknown functions; cells, with a very flexible structure (e.g., any cell can be arbitrarily subdivided without affecting its neighbors), are used for integration.A domain-based meshfree method was first proposed by Nayroles et al. [125]; this idea was improved and expanded upon by Belytschko et al. [126] who proposed the element-free Galerkin method (EFG).The EFG couples the moving least squares (MLS) approximation scheme with the Galerkin weak form to obtain a domain-based meshfree method.The BNM, first proposed by Mukherjee and Mukherjee [127], on the other hand, is a combination of the MLS approximation scheme and the standard BIE, thus producing a boundarybased meshfree method.Thus, the BNM retains the meshless attribute of the EFG and the dimensionality advantage of the BEM.As a consequence, the BNM only requires the specification of points on the 2D bounding surface of a 3D body (including crack faces in fracture mechanics problems), together with unstructured surface cells, thereby practically eliminating the meshing problem.In contrast, the FEM needs volume meshing, the BEM needs surface meshing, and the EFG needs points throughout the domain of a body.Some publications related to the BNM are [42,[128][129][130][131][132][133].Other examples of boundary-based meshfree methods are the local BIE (LBIE) approach [134,135], the boundary only radial basis function method (BRBFM) [136], the boundary cloud method (BCM) [137,138], the boundary point interpolation method (BPIM) [139], the hybrid boundary node method (HBNM) [140], the boundary face method (BFM) [141], the boundary element-free method [142], and the Galerkin boundary node method (GBNM) [143][144][145][146][147].
The BCM [137,138] is very similar to the BNM in that scattered boundary points are used for constructing approximating functions and these approximants are used with the appropriate BIEs for the problem.However, a key attractive feature of these papers is that, unlike the BNM where boundary curvilinear coordinates must be employed, the usual Cartesian coordinates can be used in the BCM.Use of Cartesian, rather than curvilinear boundary coordinates, is certainly preferable, especially for 3D problems.One disadvantage of the variable basis approach [138] (as well as the standard BNM [127]), on the other hand, is that these methods do not properly model possible discontinuities in the normal derivative of the potential function across edges and corners.Telukunta and Mukherjee [148] have tried to combine the advan-tages of the variable basis approach [138], together with allowing possible discontinuities in the normal derivative of the potential function, across edges and corners, in an approach called the extended boundary node method (EBNM).A detailed formulation for the EBNM for 2D potential theory, together with numerical results for selected problems, appear in [148]; while the EBNM for 3D potential theory is the subject of [149].A brief description of the EBNM for 3D potential theory is given below.Details are available in Ref. [149].
It is assumed that the bounding surface S of a solid body occupying the region V is the union of piecewise smooth segments called panels.The BNM employs a diffuse approximation in which the value of a variable at a boundary point is defined in terms of its values at neighboring boundary points within its domain of dependence (DOD).Correspondingly, a boundary node affects points within its range of influence (ROI).These regions are shown in Fig. 3(a).
Let /ðxÞ be the sought after harmonic function in 3D and qðxÞ be its outward normal derivative.The first step is to write: Appropriate selection of the approximation functions pðxÞ and qðxÞ (each is a vector of length m) is of crucial importance, and is discussed in detail in Ref. [149].
The coefficients a i and b i are obtained by minimizing weighted discrete L 2 norms.Finally, one gets where Here, Eq. ( 22) relates the nodal approximations of / and q to their real nodal values.The matrices A, B, C, D are defined in Ref. [149].
A moving least squares (MLS) approach is adopted in the BNM and in the EBNM.Now, one has variable weight functions within each cloud, i.e., w I ðx; x I Þ.The basic idea behind the choice of a weight function is that its value should decrease with distance In the standard BNM, the ROI of a node near an edge, e.g., node 4, is truncated at the edges of a panel.In the EBNM, the ROI can reach over to neighboring panels and contain edges and=or corners -see, e.g., node 5 (b) Gaussian weight function defined on the ROI of a node (from Ref. [149] where d ¼ gðx; x I Þ is the minimum distance measured on the surface S (i.e., the geodesic) between x and the collocation node I with co-ordinates x I ; and the quantities d I determine the extent of the ROI (the compact support) of node I.
The EBNM combines the advantages of the variable basis approach presented in Ref. [138] together with allowing discontinuities in q.A basis (i.e., the functions in p or q in Eq. ( 21)) for a cloud must satisfy two competing requirements -it must be broad enough to include all cases, yet it must be narrow enough such that the matrices A and C in Eq. ( 23) are nonsingular [149].
The next important step is to use the MLS approximations in Eq. ( 22) for the functions / and q in the regularized BIE for 3D potential theory.This regularized BIE is: The boundary S is partitioned into N c cells S k and MLS approximations (Eq.( 22)) are used in Eq. ( 25).The result is: where M I ðxÞ and M I ðyÞ are the contributions from the I th node to the collocation (source) point x and field point y, respectively.Also, n y nodes lie in the DOI of the field point y and n x nodes lie in the DOI of the source point x.When x and y belong to the same cell, the cell is treated as a singular cell and special techniques are needed to numerically evaluate the weakly singular integrals that involve the kernel Gðx; yÞ.This method is described in Ref. [42] for 3D problems.Regular Gaussian integration is used for numerical evaluation of nonsingular integrals.The final assembled equations from the EBNM are of the form: Equation ( 27) is combined with Eq. ( 22) and the given boundary values in order to solve for the unknown quantities on the boundary S.
It is well known that surface meshing (that suffices for BEM applications in linear problems) is much easier than (3D) domain meshing (as is typically required for the FEM).One might then wonder about the need for a (pseudo) mesh-free version of the BEM.It is important to point out, however, that certain situations, such as problems with moving boundaries, or those in which optimal shape design and/or adaptive meshing is carried out, typically require multiple remeshings during the solution process.This process can be painful even for the BEM.The BNM has a very flexible cell structure; the cells, used just for integration, only need to cover the surface of a body and not overlap.No other topological restrictions are imposed.Therefore, this method offers an attractive alternative to the standard BEM, especially for this class of problems.As an illustration of the power of the BNM, the initial and final cell distribution on a cube, from an ONE-step adaptive (standard) BNM calculation for 3D potential theory [131], is reproduced as Fig. 4 in this paper.This is, in fact, a Dirichlet problem with the prescribed potential on the cube surface given by The advantage of a flexible cell structure is apparent in Fig. 4(b) in which the cube surface has 1764 cells of different sizes.

2.4
Variationally Based BEM Formulations.In this section, we review the conventional, collocation BEM from a point of view that is not common in the literature, but which helps to understand some key concepts and also indicates how the BEM relates to some variationally based developments [150].This section also shows the connection of the BEM with the FEM, that is not highlighted in the literature.
We start with the derivation of the BEM.For brevity, we restrict to elastostatics, out of which the formulation of potential problems may be obtained as a special case.Let an elastic body be subjected to body forces b i in domain V and traction forces t i on part S r of the boundary S. Displacements u i are known on the complementary part S u of S. We look for an approximation of the stress field that best satisfies equilibrium in V: and r ji n j ¼ t i on S r , with a corresponding displacement field that best corresponds to u i ¼ u i on S u , where n j is the outward unit normal to S. Index notation is used here.

From Variational to Consistent Weighted-Residual
Formulations of the BEM.The present problem might be formulated in the following framework of the strong form of the principle of stationary (minimum) total potential energy [151], only assuming that r ij is a symmetric tensor that satisfies a priori the constitutive equation r ij ¼ C ijkl u k;l : for a variation du i of u i , already extending the boundary integral from S r to S, since du i ¼ 0 on S u .The weak form of Eq. ( 30) is the basis of the displacement finite element method [152] and this is explicitly assumed by the superscript ð Þ d attached to the stress field in Eq. (30).For brevity, the body force b i is assumed to be zero (see Ref. [150] for the general case).For the BEM, the above problem may be formulated in terms of weighted residuals, in a non-variational, less restrictive framework than Eq. ( 30), using as weighting functions a field of variations of fundamental solutions: The sequential presentation of the two latter equations helps to clarify the main difference between the FEM and the BEM: besides integrability issues that are inherent to any integral statement, the premise for Eq. ( 30) is that du i ¼ 0 on S u , whereas the premise for Eq. ( 31) is that du Ã i be the fundamental solution of Eq. ( 29).
The BEM is derived from Eq. (31) for variations of the following series of fundamental solutions r Ã ij and u Ã i , given in terms of force parameters p Ã m , where u r is , for s ¼ 1…n r , are n r rigid-body displacements that are multiplied by, in principle, arbitrary constants C sm .In the above equations, m characterizes both location and direction of application of p Ã m .Then, r Ã ijm and u Ã im are functions -with global support -of the coordinates and directions of p Ã m referred to by m (the source point), as well as of the coordinates and directions referred to by i (the field point), at which the effects of p Ã m are measured.For simplicity of notation, the coordinates of source and field points are not explicitly represented.
Substituting for dr Ã ijm and du Ã im in Eq. ( 31) according to Eq. (32), we obtain, for arbitrary dp Ã m , the modified expression of the Somigliana's identity, which is used to evaluate displacements u m (and, subsequently, stresses) at a domain point m for prescribed forces t i , and boundary displacements u i .The last term vanishes only if t i are in equilibrium, which is not necessarily true when we are dealing with approximations.Then, the results are in principle influenced by some arbitrary constants C sm [153].Equation (33) may also be used to evaluate displacements u i and traction forces t i as the problem's unknowns along S r and S u when the collocation point is placed on the boundary.Using boundary elements and placing the collocation point at each node, we arrive at the basic equation of the conventional, collocation boundary element method, given in matrix form as where im dS is a flexibility-type matrix.The error term e in Eq. ( 34) corresponds to residuals whose magnitude depends on the amount of rigid-body displacements that are implicit in the fundamental solution, Eq. ( 32), as well as on how refined is the boundary mesh.This vector of residuals is usually disregarded in the implementations shown in the literature [153].A consistent numerical formulation must take this term explicitly into account and end up with matrices that are independent from C sm instead of just disregarding it.Otherwise, uncontrolled ill-conditioning related to the matrix G may appear in the equation system of Eq. ( 13).This specific issue has already been the subject of a thorough theoretical investigation [150,154] and shall not be addressed in the present review.For simplicity, it is henceforth assumed that, for a sufficiently fine mesh, e ¼ 0.

The Hybrid BEM.
A variational formulation was proposed in the year 1987 as an attempt to preserve the mechanical consistency of Eq. ( 30) and the boundary-only features of the BEM [151].It turned out to be a generalization of the Hellinger-Reissner potential [152] and was called hybrid BEM in reference to Pian's developments on finite elements [155].We require that the potential be stationary for boundary displacements u d i interpolated along S and such that u d i ¼ u i on S u , and for domain stresses given by Eq. (32).The superscripts ð Þ d and ð Þ s characterize that two independent fields -of displacements and stresses -are used as trial functions.The displacement field corresponding to r s ij is given by Eq. ( 32).After substituting for u d i , r s ij and u s i in Eq. ( 35) according to the element discretization, we arrive at the following sets of matrix equations, for arbitrary variations dd and dp Ã : In these equations, H is the same matrix of the BEM, introduced in Eq. ( 34), which performs a kinematic transformation in the second set of equations and, in its preceding transpose form, an equilibrium transformation.The nodal force vector p ¼ p m ½ Ð S t i u im dS is equivalent, in terms of virtual work, to the applied forces t i on S. Any terms related to C sm are naturally void.
In Eq. ( 36), F Ã is a symmetric flexibility matrix that transforms point forces p Ã into equivalent nodal displacements d Ã , as referred to the sets of nodal quantities used to interpolate the domain stress field (fundamental solutions).The definition of F Ã is given in the following, also repeating the definition of H in order to clarify some key concepts [150,151]: The matrix F Ã shares the singularity features of the matrices H (singularity of r Ã jim ) and G (improper integral due to u Ã im ), except 031001-8 / Vol.64, MAY 2011 Transactions of the ASME when m and n refer to the same nodal point.The evaluation of the singular integral is carried out in terms of a finite part plus a discontinuous term, as indicated in the second and third lines of Eq. (37).The impossibility of obtaining the coefficients about the main diagonal of F Ã in Eq. ( 37) is evident, as the coefficients of the matrix U Ã U Ã mn are the displacement fundamental solution evaluated at point n for a source point applied at m: If m ¼ n, this would correspond to a meaningless value at a point that lies outside the domain of interest V.The response of U Ã U Ã mn for m ¼ n is actually not infinite, but rather undefined in the frame of an integral statement.The evaluation of such terms, for a finite domain, must succeed in a linear algebra framework, as briefly explained in the following.
Let W ¼ W ns ½ 2 R n d Ân r be a matrix whose columns form an orthogonal basis of the nodal displacements d of Eq. ( 36) related to the rigid-body displacements of a finite domain.Then, as by construction W ¼ NðHÞ, we may define a matrix V ¼ NðH T Þ and check that, for linear algebra consistency of Eq. ( 36), for any set of force parameters p Ã .It results that, for a finite domain V, which is a linear algebra means of evaluating the coefficients of F Ã for m and n referring to the same nodal point, which could not be obtained from Eq. ( 37) [150,151,156].Solving for p Ã in Eq. ( 36), we obtain the stiffness relation with an inverse matrix that comes from the theory of generalized inverses [150,151] The same developments outlined above apply to an unbounded domain that is complementary to V. The resulting equations help to understand several topological issues related to Eq. ( 39), for a general non-convex, multiply connected domain.Such relations have eventually led to the simplified hybrid BEM [157], which is in general as accurate as the hybrid BEM or the CBEM, although computationally less expensive [158], since the time-consuming evaluation of F Ã in Eq. ( 37) can be circumvented.
In either hybrid or simplified hybrid BEM, once d, p and p Ã are solved, for adequate boundary conditions, the results at internal points can be directly obtained via Eq.(32), thus avoiding the time-consuming evaluation of the Somigliana's identity of Eq. ( 33) [150,151].

2.4.3
The Hybrid Displacement BEM.Motivated by the hybrid BEM, Figueiredo and Brebbia proposed an alternative variational formulation, which was called hybrid displacement BEM [159].It consists of the application of the Hu potential [152], given in the following using the notation already introduced for the BEM and hybrid BEM: also with fundamental solutions as domain interpolation functions, but in terms of displacements.The compatibility between domain and boundary displacements, u s i and u d i , is weakly enforced by using the Lagrange multiplier t i .After numerical discretization, Eq. ( 41) leads to All vectors and matrices above have already been introduced, except for L, which is expressed as L 'm ¼ Ð S u im t i' dS [150].We may obtain from the above equations the stiffness relation where G ðÀ1Þ must be expressed in terms of generalized inverses, in a consistent formulation, as G is in general rectangular and possibly not full rank [150].In the latter formulation, the undefined coefficients of F Ã are obtained by means of a spectral property related to G, which takes the place of Eq. ( 39), although equivalent to it.The same topological concerns made with respect to the hybrid BEM affect the hybrid displacement BEM [150].

Fast Solution Methods
In addition to the new mathematical formulations for the BEM, as reviewed in the previous section, a few new fast solution methods for solving the BEM systems of equations also emerged in the last decade.With the developments of these fast solution methods, the BEM can now be applied to solve large-scale and practical engineering problems with the number of boundary elements on the order of a few millions on desktop computers.In this section, we will review three fast solution methods, namely, the fast multipole method (FMM), the pre-corrected fast Fourier transformation (pFFT) method, and the adaptive cross approximation (ACA) method.The first two methods were mainly developed in the U.S., while the third one was originated in Europe.Application examples of the new formulations in the BEM and these fast solution methods will be presented in Sec. 4.

Fast Multipole Method.
In the mid-1980s, Rokhlin and Greengard [170][171][172] pioneered the innovative fast multipole method that can be used to accelerate the solutions of the BEM by several orders, promising to reduce the CPU time in the FMMaccelerated BEM to O(N).We call the fast multipole accelerated BEM the fast multipole BEM (or FMBEM).Some of the early work on the fast multipole BEM in applied mechanics can be found in Refs.[173][174][175][176][177]. A comprehensive review paper on the fast multipole BIE and BEM up to 2002 can be found in Ref. [178], a tutorial paper in Ref. [179], and a textbook in Ref. [44].
The main idea of the fast multipole BEM is to apply iterative solvers (such as GMRES [180]) to solve Eq. ( 13) and use the FMM to accelerate the matrix-vector multiplication (Ak) in each iteration, without forming the entire matrix A explicitly.Direct integrations are still needed when the elements are close to the source point, whereas fast multipole expansions are used for elements that are far away from the source point.A main reason for the reduction in operations in the FMM is due to the fact that the Green's functions in the BIEs can be expanded in the following form: where y c is an expansion point near y.This can be achieved by use of various forms of expansions, including, but not limited to, Taylor series expansions.By using an expansion as in Eq. ( 44), we can write the original integral, such as the one with the G kernel in BIE (Eq.( 7)), as Applied Mechanics Reviews MAY 2011, Vol.64 / 031001-9 where S c is a subset of S away from x.In the conventional BEM, the integral is computed with the expression on the left-hand side of Eq. ( 45) directly.Any changes in the location of the source point x will require reevaluation of the entire integral.In the fast multipole BEM, when the source point x is far away from S c , the original integral is computed with the expression on the righthand side of Eq. ( 45), in which the new integrals need to be evaluated only once, independent of the locations of the source point x.That is, the direct relation between x and y is cut off by use of the expansion and introduction of the new "middle" point y c .Additional expansions and translations, as well as a hierarchical tree structure of the elements, are introduced in the fast multipole BEM to further reduce the computational costs.
As an example, we discuss the expansions in the FMM for 2D potential problems and the main procedures and algorithms in the fast multipole BEM.Consider the following integral with the G kernel in BIE (Eq.( 7)): in which S c is a subset of boundary S and away from the source point x, and Gðx; yÞ , we can write where is the fundamental solution in complex notation and Ref g indicates the real part of the function.Thus the integral in ( 46) is equivalent to the real part of the following integral: We first expand the kernel function to separate the source point z 0 and field point z.To do this, we introduce an expansion point z c close to the field point z, that is, z À z c j j( z 0 À z c j j .Applying the Taylor series expansion, we obtain where the two auxiliary functions and O 0 ðzÞ ¼ À logðzÞ.The integral in ( 49) is now evaluated by the following multipole expansion: where ÞqðzÞdSðzÞ; k ¼ 0; 1; 2; ::: are called moments about z c , which are independent of the collocation point z 0 and need to be computed only once.If the expansion point z c is moved to a new location z c 0 , we compute the moment at the new location by the moment-to-moment (M2M) translation Next, we introduce the so-called local expansion about the source point z 0 .Suppose z L is a point close to the source point z 0 , that is, From the multipole expansion in Eq. ( 51), we have the following local expansion: where the coefficients L l ðz L Þ are given by the following momentto-local (M2L) translation: If the point is moved from z L to z L 0 , we have the following localto-local (L2L) translation: Multipole expansion moments for the integrals with the F kernel, as well as those in the hypersingular BIE, can be defined similarly.The M2M, M2L, and L2L translations remain the same for these integrals.
The algorithms in the fast multipole BEM can be described as follows: Step 1. Discretize the boundary S in the same way as in the conventional BEM.For example, we apply constant elements to discretize boundary S of a 2D domain as shown in Fig. 5.
Step 2. Determine a tree structure of the elements.For a 2D domain, consider a square that covers the entire boundary S and call this square the cell of level 0. Next, start dividing this parent cell into four equal child cells of level 1. Continue dividing in this way cells that contain elements.Stop dividing a cell if the number of elements in that cell is fewer than a pre-specified number (this number is 1 in the example in Fig. 5).A cell having no child cells is called a leaf Fig. 5 A hierarchical cell structure covering all the boundary elements (shaded cells in Fig. 5).A quad-tree structure of the cells covering all the elements is thus formed after this procedure is completed (Fig. 6).
Step 3. Upward pass.Compute the moments on all cells at all levels l ! 2 and trace the tree structure upward (Fig. 6).For a leaf, Eq. ( 52) is applied directly (with S c being the set of the elements in the leaf and z c the centroid of the leaf).For a parent cell, calculate the moment by summing the moments on its four child cells using the M2M translation (Eq.( 53)), in which z c 0 is the centroid of the parent cell and z c is the centroid of a child cell.
Step 4. Downward pass.In the downward pass, we compute the local expansion coefficients on all cells starting from level 2 and tracing the tree structure downward to all the leaves (Fig. 6).The local expansion associated with a cell C is the sum of the contributions from the cells in the interaction list of cell C and from all the far cells (see Refs. [44,179] for definitions).The former is calculated by use of the M2L translation, Eq. ( 55), with moments associated with cells in the interaction list, and the latter is calculated by use of the L2L translation, Eq. ( 56), for the parent cell of C with the expansion point being shifted from the centroid of C's parent cell to that of C. Step 5. Evaluation of the integrals.Compute integrals from elements in leaf cell C and its adjacent cells directly as in the conventional BEM.Compute the integrals from all other cells (in the interaction list of C and far cells) using the local expansion (Eq.( 54)).Step 6. Iterations of the solution.The iterative solver updates the unknown solution vector k in the system Ak ¼ b and continues at Step 3 to evaluate the next matrix and vector multiplication (Ak) until the solution k converges within a given tolerance.
More details of the fast multipole BEM can be found in Refs.[44,178,179].The fast multipole algorithm discussed in this section is the original algorithm [172], which is efficient for models in which the elements are about the same size and distributed uniformly in a bulky domain.For BEM models with non-uniform element distributions and especially with large elements adjacent to smaller elements, the so-called adaptive FMMs are more efficient, in which the definitions of the adjacent cells and cells in the interaction list are further refined.Discussions on the adaptive algorithms can be found in Refs.[181][182][183][184][185].

Pre-Corrected
Fast Fourier Transformation Method.The pre-corrected-FFT accelerated technique was proposed in 1997 by Phillips and White [186,187] in an effort to develop a kernel independent acceleration technique to rapidly solve integral equations associated with Laplace and Helmholtz prob-lems.During the next decade, this technique has been improved, extended and applied together with the BEM to obtain solutions of various large-scale engineering problems governed by Laplace equation [188,189], Stokes equation [190,191], linear elasticity [192], coupled electrostatic and elasticity [193], Poisson's equation [194] and even quasi-linear equations [195].In recent years, extensions to dynamic problems have also been conducted [196,197], resulting in a more powerful and versatile method.A distinct advantage of the pre-corrected FFT technique, in addition to its kernel independent nature, is its simplicity.Unlike the FMM, this method requires only one Cartesian grid and there is no need to construct a hierarchical structure.As such, the implementation is straightforward.The complexity of the method is OðNÞ þ OðN g log N g Þ, where N is the total number of surface elements and N g is the total number of grid points.For a discussion of the complexity, readers are referred to Ref. [187].
Similar to most acceleration schemes, the pre-corrected-FFT technique achieves acceleration also by computing the far-field interaction, that is, Ak, where the source point is located far away from the field element, approximately via efficient means.Consider, for example, the evaluation of potential at point x due to charges q i at points y i distributed randomly inside a cell.The cell size is much smaller than the distances between x and y i .Instead of clustering all the discrete charges q i into one lumped charge Q and evaluating the potential at x due to the lumped charge Q as is done in the FMM, the pre-corrected FFT scheme computes the potential at x by first replacing each charge q i with point charges laying on a uniform grid; for example, the cell vertices.The grid representation allows the fast Fourier transformation to be used to efficiently perform potential computation on the grid.This grid potential is then interpolated onto point x to obtain the potential at x.In the section that follows, we use a 2D potential problem to illustrate the steps of the pFFT scheme and discuss the algorithms.
As mentioned previously, in the pFFT scheme, the far-field interaction is calculated approximately with the aid of a regular Cartesian grid.As such, the first step in this scheme is to construct the grid.This is done by first constructing a parallelepiped and superimposing it onto a 3D problem domain with its boundary being discretized into n surface elements.This parallelepiped is then subdivided into an array of small cubes so that each small cube contains only a few surface elements.The vertices of the cubes form the grid which can be further refined if needed.In addition, these cubes are also used to partition the near and the far fields.In most implementation schemes, the near field of the elements within a cube includes only the nearest neighbors of that cube.This, however, is not rigid.The optimum number of cubes should be determined based on the balance between accuracy and efficiency.For a fixed surface discretization, more cubes mean a larger far field and hence a higher efficiency, but a lower accuracy because the direct evaluation region would be smaller.Figure 7(a) shows a 2D view of a parallelepiped Fig. 6 A hierarchical quad-tree structure for the 2D boundary element mesh superimposed onto a discretized sphere.The parallelepiped has been divided into a 3 Â 3 Â 3 array of small cubes.Each cube has been further discretized to form the Cartesian grid indicated by black dots.The grid is extended outside the parallelepiped in order to ensure that the number of grid points is an even number.Note that the surface elements and the pFFT cubes can intersect with each other.There is no need to maintain any consistency between the surface elements and the cubes.
To evaluate the far-field interaction, for example, to compute the potential uðxÞ at the centroid of the purple element induced by the charge qðyÞ distributed on the yellow element as shown in Fig. 7(b), three steps are involved.First, the element charge is projected onto the surrounding grid points q located on or within the cube that encloses the yellow element.The potential at each grid point u l surrounding the evaluation element, that is, the purple element, is then computed using the fast Fourier transformation technique.Finally, the potential uðxÞ is obtained by interpolating the grid potentials u l back to the element.Mathematically this procedure is equivalent to approximating uðxÞ ¼ Ð S Gðx; yÞqðyÞdSðyÞ by uðxÞ , where W l is an interpolation operator that interpolates the grid potential u l to the potential at x and P is the projection operator that projects the charge at point y to charges at grid points y surrounding y.In the implementation of this scheme, the projection is done for all elements before the FFT is performed; hence the grid charges represent the distributed charges of the entire problem domain.The resulting grid potentials thus contain the contributions from the neighboring grid charges.To ensure accuracy, this part of the potential must be subtracted from the total potential and the contributions from the neighboring elements, that is, the near-field interactions, which are computed by directly evaluating uðxÞ ¼ Ð S Gðx; yÞqðyÞdSðyÞ, should then be added back.This step is the pre-correction step.
Both projection and interpolation can be performed using polynomial interpolation techniques and the complexity of this process is OðNÞ.The convolution step, which is conducted via FFT, requires about OðN g log N g Þ operations, where N g is the total number of grid points.For most problems, N g % N and hence the complexity of the pFFT scheme is OðN log NÞ.The direct calculation also requires OðNÞ operations.
The main steps in the pre-corrected-FFT acceleration technique are summarized as follows: (1) Construction and superposition of the 3D Cartesian grid and the discretized problem domain (Fig. 7(a)).( 2) Determination of the near and far fields for each element using the cubes of the parallelepiped; for example, for the yellow element inside cube S, its near field includes all elements inside cube S and its nearest neighboring cubes, as illustrated by the gray region in Fig. 7(b).
(3) Projection of the element charge onto the surrounding grid points based on polynomial interpolation, i.e., to compute q ¼ Ð S P qðyÞdSðyÞ ð Þ .(4) Calculation of the grid-to-grid interaction using the fast Fourier transformation, i.e., to compute u l ¼ P Gðx l ; y Þ Á q .(5) Interpolation of the grid potentials back to elements, i.e., to compute uðxÞ % P l W l Á u l .(6) Subtraction of the near-field interaction.(7) Computation of the near-field interaction using direct calculation.(8) Summation of the near-field and far-field contributions.
For a detailed description of this technique, readers are referred to Refs.[187,193].
3.3 Adaptive Cross Approximation Method.The adaptive cross approximation (ACA) follows a strategy different from the above discussed techniques to reduce the complexity of the BEM with respect to operations as well as to storage.However, the principal idea is the same and the same properties of the underlying kernel functions are used.The above discussed methodologies allow for use of the matrix-vector multiplication with almost linear complexity.However, the only approach that allows for all matrix operations (matrix-vector, matrix-matrix product, matrixmatrix addition, matrix inversion, LU-decomposition, etc.) with almost linear complexity are the so-called hierarchical matrices (H-matrices) introduced by Hackbusch [198].They can be understood an as algebraic structure reflecting a geometrically motivated partitioning into sub-blocks.Each sub-block is classified to be either admissible or not (similar to the clustering in Figs. 5  and 6).This block structure points out the fact that H-matrix arithmetic is easily parallelizable [199].
After having concluded the setup of an H-matrix, admissible blocks have to be approximated.The FMM deals with the analytical decomposition of integral kernels; hence, the procedure becomes problem dependent.The problem-independent classes are the so-called algebraic approximation methods.The singular value decomposition (SVD) leads to the optimal approximation, however, with OðN 3 Þ complexity.Less expensive algorithms are the Mosaic skeleton method developed in Ref. [200] and the successively developed ACA.It has been applied by Bebendorf [201] to the approximation of BEM matrices for the first time.The outstanding feature of ACA compared to SVD is that it requires only the evaluation of some original matrix entries and the approximation is still almost optimal.The main idea is the approximation of a matrix with a small rank k compared to t and s.This low rank representation can be found whenever the generating kernel function G(x,y) in the computational domain of A is asymptotically smooth.As shown in Refs.[202] and [203] all kernel functions G(x,y) of elliptic operators with constant coefficients and x 6 ¼ y have this property.Only in the case x ¼ y they become singular and are not smooth.
Due to the kernel independent idea, ACA can be used in a black-box like manner.Its coding and adaptation to existing codes is straight forward.The algorithm is robust and it is based on a stopping criterion depending on a prescribed approximation accuracy e. Applications in elasticity may be found in Refs.[204,205], for crack problems in [206,207], and for electromagnetism in Ref. [208], to cite a few.

Hierarchical
Matrices.Due to the previously described kernel properties of elliptic operators, the necessity to separate the near-from the far-field becomes evident.Low-rank approximations of the type (Eq.( 57)) can be obtained only for well-separated computational domains x 6 ¼ y.Thus, H-matrices [202,203] are used.Their setup is based on the following idea: The index sets I Transactions of the ASME and J of row and column degrees of freedom are permuted in such a way that those who are far away from each other do also obtain indices with a large offset.First, by means of a distance based hierarchical subdivision of I and J cluster trees T I and T J are created.In each step of this procedure, a new level of son clusters is inserted into the cluster trees.A son cluster is not further subdivided and is called to be a leaf if its size reaches a prescribed minimal size b min .Basically, two approaches can be distinguished.First, the subdivision based on bounding boxes splits the domain into axis-parallel boxes which contain the son clusters.Second, the subdivision based on principal component analysis splits the domain into well-balanced son clusters leading to a minimal cluster tree depth.Details for both approaches can be found in Ref. [203].
Now, the H-matrix structure is defined by the block cluster tree T IÂJ :¼ T I Â T J .Its setup is performed by means of the following admissibility criterion: with the clusters t & T I , s & T J , and the admissibility parameter 0 < g < 1.The diameter of the clusters t and s and their distance is computed as usual, that is, Each cluster t and s is associated with its computational domain x t and y s on C. The supports of the corresponding degrees of freedom of row i and column j are denoted by x i and y j , i.e., If Eq. ( 58) is fulfilled, a block b ¼ t Â s is admissible.If the condition in Eq. ( 58) is not fulfilled, the admissibility is recursively checked for their son clusters, until either Eq. ( 58) holds or both clusters t and s become leafs.In the latter case, block b is not admissible.Admissible blocks have well separated computational domains x t and y s and the algorithm presented in the next section is used to approximate them.Not admissible blocks must be evaluated in the standard fashion.An example for the hypersingular operator with 3279 entries is presented in Fig. 8.
In this example, 1171 clusters have been used.The red colored blocks are not admissible, whereas the green ones are admissible.The numbers represent the final low rank of the sub-blocks.They are computed with ACA, which is discussed in the next section.

Approximation of the Admissible Blocks.
A remark to the notation in this section: ðAÞ ij denotes the ij-th entry, whereas ðAÞ i and ðAÞ j are the ith row vector and jth column vector, respectively, for matrix A. The idea of ACA is to split up a matrix A 2 C tÂs into A ¼ S k þ R k where S k denotes the rank k approximation of A and R k the residuum to be minimized.Starting from a first pivot c 1 ¼ ðR 0 Þ À1 ij has to be found, where i and j are the row and column indices of the actual (0-th in this case) approximation step.Hints for the right choice of the initial pivot can be found in Ref. [202].In each ongoing step , the scaled outer product of the pivot row and column is subtracted from R and added to S with the ith row vector and jth column vector defined as The residuum R is minimized and the rank of the approximant S is increased step by step.The pivot c þ1 is chosen to be the largest entry in modulus of either the row ðR Þ i or column ðR Þ j .Finally, the approximation stops if the following criterion holds: Note, the entire matrix A will never be generated.Therefore, special care has to be taken in order to find the pivot such that the algorithm converges to the prescribed accuracy e [202,209].The numbers in Fig. 8 are the obtained ranks.As already mentioned above, a crucial point is the right choice of the pivot value.This becomes even more important in case of vector valued problems like, e.g., elasticity.In this case each entry is of matrix type A ij 2 C 3Â3 and the pivot might be defined as In no norm k k p it is guaranteed that a proper pivot entry can be found.For example, if ðR Þ ij contains one very small entry compared to the remaining entries of ðR Þ i the pivot row U is scaled up by c and Eq. ( 63) does not hold.If ðR Þ ij contains a zero entry the pivot c is not even defined.This becomes evident if, e.g., the double layer operator, Eq. ( 5) is evaluated on some plane C k fx; y 2 R 3 : y k À x k ¼ 0g that lies perpendicular to the coordinate axis k.A work around is the partitioning of the system matrix with respect to each degree of freedom (for details, see Ref. [210]).the FEM package ABAQUS for mechanical analysis together with a BEM code FastCap [212] for electric field analysis.

Emerging Applications
A generic MEMS structure, consisting of two parallel conducting beams, is shown in Fig. 9.The electric potential /, in the region V exterior to the two beams, is governed by Laplace's equation (Eq.( 1)) and is prescribed on the boundary S of the two beams (it equals / 1 on the surface of beam 1 and À/ 1 on the surface of beam 2).
The corresponding BIE in 2D space is of the form (see, e.g., Ref. [213]) where r ¼ eð@/=@nÞ ¼ eq is the charge density at a point on a beam surface and e is the dielectric constant of the exterior medium.The total charge Q on the surface S is Given the potential /, and the total charge Q, Eqs. ( 65) and ( 66) can be solved for the charge density r on the total beam surface S.This charge density is then used to determine the traction h at any point on the beam surface according to The mechanical problem of deformation of the beams, subjected to traction h, is generally solved by the FEM.One can employ an iterative procedure, as has been done by many authors, or a total Lagrangian approach, first proposed in Ref. [214].

Inertial Sensors.
Another successful application of the BEM in the modeling and design of MEMS/NEMS is air damping analysis of inertial sensors, like the currently widely used accelerometers and gyroscopes, and other types of micro resonators.Micro resonators are a class of MEMS devices with applications ranging from various physical, chemical and biological sensors to high-frequency filters, oscillators in wireless communication systems.In a typical micro resonator, a mechanical structure is driven into its resonance by an integrated micro actuator.The vibration of the structure is then measured either electrically or optically and its resonant frequency and the quality factor, defined as the ratio of the total energy to the dissipated energy within one cycle, can be extracted.In most applications, the sensing mechanism relies on the fact that the change in the resonant frequency is related to the physical quantities to be measured.The resolution of these sensors thus depends strongly on air damping, a dominant dissipation source for devices operating in open air and even in a low vacuum.
Air damping is determined by gas transport and its interaction with the structure.A key parameter that influences the transport is the Knudsen number (Kn), defined as the ratio of the mean free path of gas molecules to the characteristic length of the fluid field.Based on the value of Kn, gas can be classified roughly into four regimes, the continuum regime (Kn < 0.01), the slip regime (0.01 < Kn < 0.1), the transition regime (0.1 < Kn < 10) and the free-molecule regime ((Kn > 10).In the free-molecule regime, gas transport is dominated by the collisions between molecules and structures, while in the continuum and slip regimes it is dominated by the collisions between gas molecules.This fact results in different governing equations and hence requires different modeling approaches for different gas regimes.
Continuum and slip regimes.For resonators oscillating at relatively low frequencies, the incompressible Stokes equation describes accurately gas behavior in both continuum and slip regimes.In the absence of body force, it reads where u and P denote the amplitudes of the velocity and pressure of gas, q is the density of gas, l is the coefficient of viscosity, and x is the oscillation frequency.The equivalent integral equation of Eq. ( 68) can be formulated as In Eq. ( 69), U and T are the free-space Green's functions of oscillatory Stokes flows [231] and t are tractions exerted by the structure on the fluid.The boundary conditions in the continuum regime follow the classical no-slip conditions, that is, the tangential gas velocity equals the wall velocity at the boundary.In the slip regime, gas rarefaction causes slip between the tangential velocities of the flow and that of the solid surface (the wall).A phenomenological slip model is shown in Eq. (70).In a non-dimensional form, the slip velocity at an isothermal and locally flat wall is given as where U g ; U w are normalized tangential velocities of the flow and the wall at the interface, respectively, @Ug=@n is the normal derivative of the normalized tangential velocity at the interface, r v is the tangential momentum accommodation coefficient, and b is a high-order slip coefficient [232].An integral equation for the normal derivative of the velocity is thus required to solve Stokes equation with slip boundary conditions.The first implementations of Eq. ( 69) have been proposed in the quasi static limit where inertia forces can be neglected [190,[233][234][235].The effect due to the inertial force is addressed in Ref. [190].In the quasi-static case, Eq. ( 69) reduces to a classical Dirichelet problem of incompressible elasticity but turns out to be severely ill-conditioned.This issue has been resolved in Refs.[233][234][235] with a Burton-Miller like approach coupling the velocity and the traction BIEs.
The extension to the moderate-frequency regime and to the slip regime has been addressed in Refs.[191,236] and acceleration techniques such as the pre-corrected FFT technique [190,191,237,238], Transactions of the ASME the fast multipole method [233][234][235][236]239] and wavelets [240] have been employed to accelerate the computation.
Transition and free molecule regimes.The situation is, however, less developed and understood in the transition regime, based on the formal classification given above.Even though a boundaryvolume integral equation approach has been formulated theoretically in Ref. [241], due to the complicated collision integral the BEM is not suitable and other techniques of deterministic and statistical nature are generally preferred.On the contrary, the collisionless (or free-molecule) flow lends itself again to the development of robust and competitive BIE approaches.The free-molecule flow represents the limiting case where the Knudsen number tends to infinity and collisions between molecules can be neglected.In inertial MEMS, this applies typically at pressures in the range of a few mbars and below, as is often encountered in package sealed sensors (e.g., gyroscopes, magnetometers and resonators).
Typical techniques for the simulation of gas dissipation in this regime belong to the class of test particle Monte Carlo (TPMC) methods [241] where the trajectories of large ensembles of molecules are analyzed with suitable models for gas-surface interaction in order to compute averages of the quantities of interest.Compared to deterministic techniques, the computational cost of TPMC is much higher due to its statistical nature.In Refs.[242,243] a deterministic technique based on integral equations has been put forward and has been shown to be indeed very competitive with TPMC methods for the typical working conditions of inertial MEMS.

Modeling of Fiber-Reinforced Composites.
In this section, we review the BEM formulation for multi-domain elasticity problems [44], that can be applied to model fiber-reinforced composite materials, functionally graded materials, and other inclusion problems in elasticity.
Consider a 2D or 3D elastic domain V 0 with an outer boundary S 0 and embedded with n elastic inclusions V a with interface S a , where a ¼ 1; 2; :::; n (Fig. 11).Note: In this section, n is used to indicate the number of inclusions.For the matrix domain V 0 , we have the following BIE: Hypersingular or traction BIE can also be applied in the matrix as well as in the inclusion domains.In fact, the dual BIE formulation (a linear combination of the conventional BIE and hypersingular BIE) is preferred for modeling inclusion problems in which thin shapes often exist and could present difficulties for the conventional BIE formulation when it is applied alone.
Assume that the inclusions are perfectly bonded to the matrix, that is, there are no gaps or cracks and no interphase regions, we have the following interface conditions: for a ¼ 1; 2; :::; n.That is, the displacements are continuous and the tractions are in equilibrium at the interfaces.Based on these assumptions, we can write the discretized form of the multidomain BIEs as follows, which gives faster convergence with iterative solvers [44]:  . . . in which u 0 and t 0 are the displacement and traction vector on the outer boundary S 0 , u i and t i are the displacement and traction vector on the interface S i (viewed from the matrix domain), A ij and B ij are the coefficient submatrices from the matrix domain, and A f i and B f i are the coefficient submatrices from inclusion i.Studies of 2D composite models with periodic conditions can be found in Refs.[244,245].

Inclusion S n
The fast multipole BEM and other fast solution methods can be applied to solve the BEM systems of equations shown in Eq. ( 74) for the inclusion problems.The multipole expansions and related translations for 2D and 3D elasticity BIEs can be found in Refs.[44,[246][247][248][249][250].Some other related work on the fast multipole BEM for general 3D elasticity problems can be found in Ref. [251], for 3D inclusion problems in Refs.[175,252], and for crack problems in Refs.[176,253,254].In the following, we present two examples using the BEM to show its potentials in modeling composite materials.
Figure 12 shows 2D BEM models of fiber-reinforced composites used to extract the effective properties in the transverse directions.The fibers are distributed randomly and with an decreasing density in one direction.Thus, the composites may be considered as functionally graded materials.The plane-strain models are used to extract the moduli of the composites in the transverse directions.Two cases are studies, one with fibers in a circular shape and the other with fibers in a star shape.The computed effective Young's moduli for the two composite models are reported in Ref. [44], which agree with the estimates given in Ref. [255].
Figure 13 shows a 3D representative volume element (RVE) model with 5832 fibers and 10,532,592 DOFs for a composite.The volume fraction of the fibers is 3.85%, the fibers are distributed randomly.Theoretically, the moduli of the materials should be independent of the sizes of the RVE models.However, when the RVE sizes are small (with only a few fibers), there are significant changes in the estimated moduli as we increase the model sizes, which suggests that the models may not be representative of the composites [249].The estimated moduli approach constant with the increase of the fibers in the RVE as revealed in the results [249].

Modeling of Functionally Graded Materials.
In general, a functionally graded material (FGM) is a special type of composite in which the constituent volume fraction varies gradually leading to a non-uniform microstructure with continuously graded macroproperties, e.g., thermal conductivity, density, and specific heat.For instance, for problems governed by potential theory, e.g., steady-state heat transfer, the thermal conductivity is a function of position, i.e., k :k(x).Approaches to treat problems of potential theory in non-homogeneous media include the Green's function approach, domain integral evaluation, and variable transformation approach.In the context of the BEM, problems of potential theory in non-homogeneous media have been previously studied by Cheng [256], Ang et al. [257], Shaw and Makris [258], Shaw [259], Harrouni et al. [260], Divo and Kassab [261], and recently by Gray et al. [262] and Dumont et al. [165].Most of these works are focused on obtaining the Green's function.In the Green's function approach, the Green's function has to be derived and a boundary only formulation can be obtained.
Several approaches have been developed to evaluate domain integrals associated with boundary element formulations including approximate particular solution methods, dual reciprocity methods, and multiple reciprocity methods.The particular solution methods and dual reciprocity methods can be considered more or less to be equivalent in nature.These methods have been widely used on the axiom that the domain integral in the boundary integral formulation is somehow removed.In these methods the inhomogeneous term of the governing differential equation is approximated by a simple function such as (1 þ r) or radial basis functions (RBFs).The mathematical properties and the convergence rates of the RBF approximations have been studied extensively.In these techniques, the boundary-only nature of the BEM is compromised.
A transformation approach, called the "simple BEM," for potential theory problems in non-homogeneous media where nonhomogeneous problems are transformed into known problems in homogeneous media can be used [263,264].The method leads to a pure boundary-only formulation.This idea has been successfully implemented in three dimensions for steady state, transient heat conduction problems and crack problems where the material property varies in one, two, and three dimensions.There is a class of material variations, which can transform the problem equation to a Laplace or standard/modified Helmholtz equation.The idea of the "simple BEM" consists of transforming problems in nonhomogeneous media to known problems in homogeneous media such that existing codes for homogeneous media can be reused with simple modifications.By means of this approach, which consists of simple changes in the boundary conditions of existing homogeneous heat conduction computer codes, the solutions for Transactions of the ASME nonhomogeneous media with quadratic, exponential, and trigonometric material variations can be readily obtained.
For problems described by potential theory, the governing differential equation for a potential function / defined on a region X bounded by a surface C, with an outward normal n, can be written as r Á ðkðx; y; zÞr/Þ ¼ 0 (75) where k(x, y, z) is a position-dependent material function and the dot represents the inner product.Equation ( 75) is the field equation for a wide range of problems in physics and engineering such as heat transfer, fluid flow motion, flow in porous media, electrostatics, and magnetostatics.The boundary conditions of the problem can be of the Dirichlet or Neumann type: respectively, with C ¼ C 1 þ C 2 for a well-posed problem.

Variable transformation approach:
By defining the variable one rewrites Eq. ( 75) as or, alternatively where By setting k 0 ðx; y; zÞ ¼ 0; þb 2 ; and À b 2 , three classical homogeneous equations namely, the Laplace, standard Helmholtz and the modified Helmholtz can be obtained, respectively.From these cases, a family of variations of k(x, y, z) can be generated, as given in Refs.[263,264].From an engineering point of view (for applications such as FGMs), material variation in one coordinate is of practical importance, which include thermal barrier coatings, bone and dental implants, piezoelectric and thermoelectric devices, graded cementitious materials, and optical application with graded refractive indexes.

Boundary conditions:
In order to solve the boundary value problem based on the modified variable v, the boundary conditions of the original problem have to be incorporated in the modified boundary value problem.Thus, for the modified problem, the Dirichlet and the Neumann boundary conditions given by Eq. ( 76) change to the following: Notice that the Dirichlet boundary condition of the original problem is affected by the factor ffiffi ffi k p .Moreover, the Neumann boundary condition of the original problem changes to a mixed boundary condition (Robin boundary condition).This later modification is the only significant change on the transformed boundary value problem.
A numerical example of an FGM rotor with eight mounting holes having an eight-fold symmetry is presented (Fig. 14).Owing to the symmetry, only one-eighth of the rotor is analyzed.The top view of the rotor, the analysis region, and the geometry of the region are illustrated in Fig. 14(a).The grading direction for the rotor is parallel to its line of symmetry, which is taken as the zaxis.The thermal conductivity for the rotor varies according to the following expression: The BEM mesh employs 1584 elements and 3492 nodes.A schematic for the thermal boundary conditions and the BEM mesh employed is shown in Fig. 14(b).Here, the solution of the problem is verified using the commercial software ABAQUS.The FEM mesh consists of 7600 20-node (quadratic) brick elements Fig. 13 A 3D RVE with 5832 fibers and with the total DOFs 5 10,532,592 [249] and 35,514 nodes.Figure 14(c) shows the comparison of the BEM and FEM results for the temperature around the hole.A contour plot of the temperature distribution using the BEM is shown in Fig. 14(d).All the results of the BEM and FEM solutions are in very good agreement.

Modeling of Fracture Mechanics Problems.
Fracture mechanics is one of the areas where the BEM has been clearly shown to be a powerful and effective numerical method when compared to other computational techniques.Among its advantages one may cite that: (1) discretization of only the boundary is required; thus simplifying preprocessing and remeshing in crackgrowth analysis, (2) it shows improved accuracy in stress concentration problems, since there are no approximations imposed on the stress solution at the interior domain points, (3) fracture parameters (stress intensity factors (SIF), energy release rates, etc.) can be accurately determined from the computed nodal data in a straightforward manner, and (4) modeling of problems involving infinite and semi-infinite domains (like in wave propagation phenomena) is simple and accurate, because the radiation conditions at infinity are automatically satisfied.
However, the direct application of the displacement BIE to cracks where the two surfaces coincide leads to a mathematically degenerate problem.This has been a sticking issue in the BEM for modeling fracture problems (see, e.g., Ref. [36]).Different approaches have been adopted in the BEM to tackle fracture problems, such as the displacement discontinuity approach [265,266] and the use of hypersingular BIEs [267][268][269].
Other alternative approaches also exist.The first approach is based on the adoption of a Green's function derived for an infinite medium containing a crack with the same geometry and boundary conditions (generally traction-free) as in the problem to be solved.In this way, there is no need to mesh the crack, since the fundamental solution already accounts for it.The other approach divides the domain into several regions via the definition of artificial interfaces that contain the cracks; thus it links the cracks to the boundary of the resulting subdomains.Equilibrium and compatibility conditions are further enforced on the nodes of the newly introduced interfaces to solve the problem.
These two approaches have limited applicability: ad-hoc fundamental solutions are needed in the first case, while the artificial boundaries required in the second approach are not uniquelydefined and they also increase the number of unknowns.Therefore, none of them is suitable for large-scale BEM models involving a great number of cracks.Following the works of Hong and Chen [270], Krishnasamy et al. [267][268][269], and Portela et al. [271], an alternative approach is further proposed to overcome these difficulties, that is, the dual or hypersingular BIE approach that results in a single region formulation by applying the displacement BIE on one of the crack faces and the traction BIE (that follows from differentiation of the displacement BIE and substitution into the constitutive law) on the other.
A review of the major contributions of the BEM in fracture mechanics between 1972 and 1997 may be found in Ref. [272], with applications covering from linear elastic fracture mechanics to SIF computation, dynamics, composite materials and interface cracks, thermoelastic problems, crack identification techniques, and nonlinear problems.
Later efforts have focused on applying the BEM to fracture of multifield materials with coupled piezoelectric and/or elastomagnetic Fig. 14 Analysis of an FGM rotor part using the BEM [263] 031001-18 / Vol.64, MAY 2011 Transactions of the ASME properties [273][274][275][276][277][278], crack growth analysis using the BEM alone or together with FEM [279][280][281], or addressing fracture phenomena in composites or at different scales [282][283][284][285].However, when dealing with large-scale fracture problems, the BEM demands great computer memory storage and significant CPU time unless algorithms to accelerate the solution come into play.The development of fast BEM approaches in the last 15 years oriented towards applied mechanics problems has revived the potential applications of the BEM to fracture mechanics.These techniques include the fast multipole method [44,170,172,178], panel clustering method [286], mosaic-skeleton approximation [287] and the methods based on hierarchical matrices and adaptive cross approximation [198,202].All these techniques aim at reducing the computational complexity of the matrix-vector multiplication, which is the core operation in iterative solvers for linear systems.FMM and panel clustering address the problem from an analytical point of view and require the knowledge of some kernel expansion prior to integration.Meanwhile, mosaic-skeleton approximations and hierarchical matrices are based on purely algebraic tools for the approximation of boundary element matrices which are thus suitable for applications where analytical closed-form expressions of the kernels are not available, as in anisotropic fracture applications.
Nishimura et al. [176] applied the original FMM by Rokhlin [170] to solve crack problems for 3D Laplace's equation.Yoshida et al. [288] further improved the procedure by implementing Greengard and Rokhlin's diagonal form [289], leading to a CPU time reduction of about 30% as compared to the original FMM when a dense array of cracks is distributed in the domain.In addition, Yoshida et al. [253] developed a Galerkin formulation that took advantage of the improved accuracy of the Galerkin method.
Helsing [290] presented a 2D complex variable formulation based on an integral equation of Fredholm second kind together with the FMM.The potential of the resulting approach is illustrated by accurately computing stress fields in a mechanically loaded material containing 10,000 randomly oriented cracks.Wang et al. [291] improved the rate of convergence of fast multipole BEM by considering preconditioners based on sparse approximate inverse type.Convergence within 10 steps is achieved for a 2D elastic infinite body containing 16,384 cracks with 1,572,864 DOF.Liu [247] developed a new FMM with the dual BEM for 2D multi-domain elasticity based on a complex variable representation of the kernels.The resulting formulation was more compact and efficient than previous FMM approaches, as shown in tests involving a large number of crack-like inclusions.
For the application of the ACA method in the modeling of crack problems, Kolk et al. [292] investigated 3D fatigue crack growth using a dual BEM.In order to accelerate the solution, the resulting collocation matrix is compressed by a low-rank approximation performed by the ACA method.The efficiency of the formulation is illustrated with an industrial application in which a corner crack located at the dedendum of a gear propagates.This approach was later revisited and improved by minimizing the computation time of each incremental loop needed to simulate the crack growth process [293].
More recently, Benedetti et al. [294] developed a fast solver based on the use of hierarchical matrices for the representation of the collocation matrix and implemented it for dual BEM analysis of 3D anisotropic crack problems.The low-rank blocks were computed by the ACA.Numerical tests demonstrated that the solver was rather insensitive to the degree of anisotropy of different materials, which results in significant savings in memory storage and assembly and solution time.An order of magnitude improvement has been observed for a system of order 10 4 .Benedetti and Aliabadi [295] further extended the approach to the dynamic analysis of 3D cracked isotropic bodies in the Laplace transform domain.Based on this approach, Benedetti et al. [296] developed a model to analyze 3D damaged solids with adhesively bonded piezoelectric patches used as strain sensors.The cracked structure was modeled using a fast dual BEM with a hierarchical ACA solver, which permitted the accurate analysis of crack parameters.
Figure 15 shows an example of large-scale fracture analysis using the fast multipole BEM.In this case, arrays of 12 Â 12 Â 12 penny-shaped cracks in an elastic domain were solved using the BEM.A total of 1,285,632 DOFs is used in this BEM model.The crack opening displacements were computed and plotted in the figure [297].
4.5 Acoustic Wave Problems.The BEM has been applied to solve acoustic and other wave propagation problems for several decades, due to its features of boundary discretizations and accurate handling of radiation conditions at the infinity for exterior problems.The developments of fast methodologies in the last decade for solving large-scale acoustic problems are perhaps the most important advances in the BEM that have made this method unmatched by other methods in modeling such problems.
The fast multipole method developed by Rokhlin and Greengard [170][171][172] has been extended to solving the Helmholtz equation governing the acoustic wave problems for quite some time (see papers in Refs.[183,184,[298][299][300][301][302][303][304][305][306][307][308][309], a review in Ref. [178], and Chapter 6 of Ref. [44]).Large-scale acoustic BEM models with a total number of DOFs around several millions have been attempted by the fast multipole BEM on desktop PCs.For example, acoustic BEM models of entire airplanes [44,309,310], wind turbines [184], and a submarine [44,310] have been solved for noise predictions.Coupled structural acoustic problems have also been solved successfully using the fast multipole BEM with the FEM for analyzing ship structures [311], and, most interestingly, for computing physically based and realistic sound in computer animations [312,313].
The adaptive cross approximation suited for acoustics has not found so much attention like the FMM until recently.There are mathematical studies concerning ACA applied in acoustics which show that increasing the frequency decreases the efficiency [314].The counterpart in the FMM is that for higher frequencies, more terms in the series expansion have to be used.One application of ACA in acoustics in time domain can be found in [210].An impressive example in determining the noise level in an aircraft can be found in Ref. [315].Preliminary comparison studies on the performance of the fast multipole BEM and ACA BEM in modeling acoustic problems are provided in Refs.[315,316].In the following, we first review the BIE formulation for solving acoustic wave problems and then show a couple of applications.
We consider the linear time-harmonic acoustic wave problems in domain E, which can be an interior domain, or an exterior domain outside a structure with boundary S. The acoustic pressure / is governed by the Helmholtz equation and boundary conditions given on S, such as pressure or velocity conditions.For exterior (infinite domain) acoustic wave problems; in addition to the boundary conditions on S, the field (radiated or scattered wave) at infinity must also satisfy the Sommerfeld radiation condition.
Boundary integral equations have been applied in solving acoustic wave problems for decades [317][318][319][320].In the direct BIE formulation, the representation integral of the solution / is (for the derivation, see Refs.[320][321][322] and Chapter 6 of Ref. [44]) where x 2 E, q ¼ @/=@n, / I ðxÞ is a possible incident wave for an exterior problem, x the circular frequency, and the two kernels for 3D are given by Gðx; y; Fðx; y; xÞ ¼ @Gðx; y; xÞ @nðyÞ where k is the wavenumber.Let the source point x approach the boundary.We obtain the following conventional boundary integral equation (CBIE) [320][321][322]: where cðxÞ ¼ 1/2, if S is smooth around x.This CBIE can be employed to solve for the unknown / and q on S. The integral with the G kernel is a weakly singular integral, while the one with the F kernel is a strongly singular (CPV) integral, as in the potential case.A regularized or weakly-singular form of CBIE (Eq.( 86)) exists and can be derived by use of the identities for the static kernels [44,46,47,323].
It is well known that the CBIE has a major defect for exterior problems, that is, it has non-unique solutions at a set of fictitious eigenfrequencies associated with the resonate frequencies of the corresponding interior problems [318].This difficulty is referred to as the fictitious eigenfrequency difficulty.A remedy to this problem is to use the normal derivative BIE in conjunction with the CBIE.Taking the derivative of integral representation (Eq.( 83)) with respect to the normal at a point x on S and letting x approach S, we obtain the following hypersingular boundary integral equation (HBIE): In HBIE (Eq.( 87)), the integral with the kernel K is a strongly singular integral, while the one with the H kernel is a hypersingular (HFP) integral.HBIE (Eq.( 87)) can also be written in a regularized or weakly singular form [321,322].For exterior acoustic wave problems, a dual BIE (or composite BIE [321]) formulation using a linear combination of CBIE (Eq.( 86)) and HBIE (Eq.( 87)) can be written as where b is the coupling constant.This formulation is called Burton-Miller formulation [318] for acoustic wave problems and has been shown to yield unique solutions at all frequencies, if b is a complex number (which, for example, can be chosen as b ).As the first example, a submarine model (Fig. 16) is presented which is solved by use of the fast multipole BEM (FastBEM AcousticsV R , V.2.2.0).This is an interesting example in solving large-scale underwater acoustic problems, which has been a challenging task for other domain-based methods.The Skipjack submarine is modeled, which has a length of 76.8 m and a total of 250,220 boundary elements are used in the discretization, with a typical element size equal to 0.14 m.Velocity BCs are applied to the propeller and an incident plane wave in the direction (1, 0, À1) is specified.The model is solved at a non-dimensional wavenumber ka ¼ 384 (frequency f ¼ 1233 Hz).The computed sound pressure on the surface of the model is shown in Fig. 16.The BEM model was solved in 85 min on a DellV R PC with IntelV R Duo Core CPU and with the tolerance for convergence set at10 À4 .
Next, we show an example using the ACA BEM for solving large-scale 3D acoustic problems.Figure 17 shows a BEM model of five wind turbines discretized with 557,470 boundary elements.The wind turbines have the same height of 29 m and are spaced .This large BEM model (with more than 500,000 complex-valued DOFs) was solved in 210 min using the ACA BEM on a Dell PC with an IntelV R Duo Core CPU.For comparison, the same large-scale acoustic model was solved on the same PC in 108 min using the fast multipole BEM [184].At present, the fast multipole BEM is faster than the ACA BEM for models with more than about 100,000 DOFs and the fast multipole BEM also uses smaller memory for storage.However, the ACA BEM has the kernel-independent advantage and is easier to implement.
4.6 Elastic Wave Problems.Boundary element methods in elastodynamics have a long history.Actually, they date back to 1930s or 1940s where Kupradze [324] started his theoretical work on singular integral equations for elastodynamics.Kupradze's [324] contributions are not restricted to theoretical works.Indeed, his numerical method known as Kupradze's functional equation method (collocate direct BIEs in the exterior of the domain under consideration) is not just well-known, but also has been rediscovered by many authors!It is noteworthy that Cruse and Rizzo, who are usually considered to be the founders of BEM (called boundary integral equation methods by these authors), investigated the use of the BEM in elastodynamics [6] soon after their first work in elastostatics.It is probable that these pioneers were aware of the potential of the BEM in elastodynamics since it can deal with exterior problems easily.Indeed, the BEM in wave problems gives highly accurate results which are not polluted by reflections from artificial boundaries, while many other numerical methods require special techniques to avoid such reflections.Subsequent developments of the BEM in elastodynamics were quite rich with many applications in areas such as nondestructive testing or evaluation, earthquake engineering and others, as one can see in the review article by Kobayashi [325].
As in other branches of the BEM, however, the O(N 2 ) complexity of the algorithm turned out to be a serious problem.With this problem, ordinary BEMs were found to be inefficient in solving large problems in spite of the boundary-only nature of the approaches.This problem, however, is now almost solved with the development of the so called fast BEMs.Indeed, the fast methodologies are undoubtedly among the most important developments in BEMs in elastodynamics in the last decade, which will be the subject of the following subsection.This will be followed by another subsection devoted to stability issues and inverse problems.

Fast BEMs in Elastodynamics
Frequency domain.Fast methods in elastodynamics are usually obtained from certain fast methods for Helmholtz equation in frequency domain or from those for the scalar wave equation in time domain.In frequency domain; for example, one notes that the fundamental solution C ij for elastodynamics is written as yÞ þ e ipr e jqr @ @x p @ @y q G kT ðx; yÞ !
where G k and k L;T are the fundamental solution of Helmholtz equation and the wave numbers of P and S waves given by and q, ðk; lÞ and x are the density, Lame ´'s constants, and frequency, respectively.This shows that one may formulate fast methods for elastodynamics simply by differentiating those for the Helmholtz equation.
The first fast method used with the BEM in elastodynamics seems to be the fast multipole method.Indeed, Chew's group [326] presented a 2D FMM formulation extending the highfrequency FMM (diagonal form) in Helmholtz proposed by Rokhlin [298].They also showed an example of the scattering from a rough interface between two elastic materials.Fukui [327] proposed a low-frequency FMM in 2D with which they solved scattering problems for many holes.Fujiwara [328] also proposed a low-frequency FMM in 2D and solved scattering problems for many cavities and cracks.
In 3D, Fujiwara [329] presented a diagonal form FMM, but he solved low-frequency problems related to earthquakes.Yoshida et al. [330,331] considered low-frequency crack problems in 3D, which they solved with a low-frequency FMM utilizing a decomposition of the kernel in Eq. ( 91).Yoshida et al. [332] also proposed a diagonal form version of the elastodynamic crack analysis.Unfortunately, their error control scheme was immature.As a matter of fact, the use of the diagonal form is effective only when the wavenumber times the cell (box) size used for FMM is above a certain threshold.In problems where the smaller cells do not satisfy this criterion, one has to switch to the low-frequency FMM in order to maintain the accuracy.Namely, one may have to use both diagonal forms and low-frequency FMM in one problem.This type of FMM formulation is called a wideband FMM [309,332].
Subsequent developments in elastodynamic FMMs were rather slow until 2007 where Chaillat et al. [334,335] reconsidered diagonal forms in elastodynamics.Sanz et al. [336] utilized several techniques including SPAI preconditioners.Tong and Chew [337] applied techniques developed for Maxwell's equations to the S wave part in their diagonal form formulation.Chaillat et al. [338] considered half-space and multi-domain problems for seismological applications.Isakari et al. [339] presented a low-frequency FMM for periodic inclusion problems in elastodynamics which may be used in applications related to phononic crystals.These Fig. 17 Sound field from five wind turbines evaluated using the ACA and fast multipole BEM authors also considered pre-conditioners based on Calderon's formulae [340].
FMMs are not the only acceleration methodology useful in elastodynamic BEMs.Indeed, one may mention the use of pre-corrected FFT approach by Yan et al. [197], and the H-matrix (ACA) approach in crack problems by Benedetti and Aliabadi [295].
4.6.1.2Time domain.As we have seen, papers on fast methods in elastodynamics in frequency domain are scarce in spite of the importance of the subject in applications.Still more scarce are the time-domain counterparts.Takahashi et al. [341] extended the so called PWTD (plane wave time domain) approach developed by Ergin et al. [342] for the wave equation to elastodynamics in 3D.Roughly speaking, this approach is the time-domain counterpart of the diagonal forms in the frequency domain with the Hankel functions in the translation operators replaced by Bessel functions.This approach has been parallelized and extended to the anisotropic case by Otani et al. [343].Figure 18 shows a typical example from this investigation, i.e., the response of 1152 spherical cavities in an infinite isotropic elastic solid subject to a plane P wave.Yoshikawa et al. [344] applied this approach to a simple inverse problem of crack determination.
Another fast approach in time domain has been proposed by Fukui's group [345], in which these authors accelerated a convolution quadrature formulation in 2D elastodynamics with a lowfrequency FMM.

Other Developments in Elastodynamic BEM
Stability issues.Stability in the time domain approach still seems to be an issue, as can be seen from the continuing publications on this subject.See, for example, the work by Soares and Mansur [346] and a very recent paper by Panagiotopoulos and Manolis [347].The former approach works as a filtering to the recent-in-time convolution operation and is related to similar efforts in electromagnetics (e.g., Rynne and Smith [348]).We also note that mathematical proofs of the stability of some variational approaches are available (e.g., the work by Ha-Duong et al. [349]), although collocation methods are of more interests to engineers.
Convolution quadrature method (CQM) is a new approach for the time-domain BEM which improves the stability of the method.It also extends the applicability of the BEM in the time domain because CQM needs only the Laplace transforms of the fundamental solutions with respect to time, which may be easier to obtain than their time-domain counterparts.We here mention Kielhorn and Schanz [350] in addition to Fukui's contribution [345].See the section on time domain approaches for more references.
Inverse problems.Applications to problems related to ultrasonic non-destructive evaluation have been among major motivations of developing BEMs for elastodynamics.The past decade has seen many interesting applications of BEMs to inverse problems of determining defects in structures, of which we cite just a few in the rest of this section.
Solution techniques of inverse problems of determining unknown defects in an elastic structure by BEMs can be classified into two categories, namely, those which require successive solutions of direct problems and those which do not.Methods of the first category typically determine the unknown defects by minimizing cost functions which represent the difference between the experimental data and the computed value of certain physical quantities.Approaches of this type are considered suitable for determining the detail of the defects, although they may not necessarily be effective for obtaining global information such as the location or the number of defects.For methods of this category, one needs effective ways of evaluating the sensitivities of the cost functions with respect to parameters which characterize the defects.Adjoint variable (field) methods are known to be effective for such a purpose particularly when one has many parameters to identify the defects.For example, Bonnet et al. [351] considered shape sensitivities for cavities and cracks using the time-domain BEM.Bonnet and Guzina [352] discussed similar problems for inclusions in the frequency domain.
As examples of inverse techniques of the second category, we mention topological derivatives and linear sampling methods.Topological derivatives are essentially the rates of changes of a certain cost function when an infinitesimal cavity is produced at a given position in a structure.This concept turned out to be useful as a "void indicator" which can be used in order to obtain a rough idea about the defects, as we can see in Bonnet and Guzina [353].In the near-field linear sampling method discussed in Fata and Guzina [354], these authors consider an elastic half space with unknown defects.To determine the defects, they excite a part of the plane boundary and measure the response on another part of the plane boundary.They can delineate the defects by examining the blow-up of solutions of certain integral equations whose kernel functions are determined from the experimental data.An adjoint of this formulation is easier to use when one has only a small number of excitations (sources), as has been discussed in Fata and Guzina [355].

Electromagnetic
Wave and Periodic Problems.Electromagnetics (EM) is undoubtedly one of academic/industrial areas where the BEM is used widely.This is partly because there are many good applications where the BEM is quite effective.For example, propagation of electromagnetic waves often takes place in an open air, which is a loss-less homogeneous exterior domain.The BEM is considered to be one of ideal tools for solving such problems and Maxwell's equations describe the phenomenon perfectly.It is therefore not surprising that the BEM and related methods are accepted enthusiastically in the EM community.Indeed, BEM software for EM applications including those accelerated with FMM are now available commercially or even freely.
There are many techniques common to the BEM in EM and in mechanics, and the mechanical community has a lot to learn from developments in EM community.It is therefore deemed worth the efforts for those from mechanics to have some knowledge on the BEM for Maxwell's equations.The first part of this section is thus intended as a quick introduction to the BEM for Maxwell's equations for those who are familiar with the BEM in mechanics.The interested reader is referred to Chew et al. [356] for further details.We then discuss more recent topics on the periodic FMM for Maxwell's equations, which have many interesting applications such as photonic crystals and metamaterials.Finally we review where E and H are the electric and magnetic fields, x is the frequency (with e Àixt time dependence), and e and l are the dielectric constant and the magnetic permeability for the material occupying D, respectively.We may combine these equations to have In addition, we have an incident wave denoted by (E inc , H inc ) in D and scattered waves E sca ¼ E À E inc and H sca ¼ H À H inc which satisfy radiation conditions.The fundamental solutions C ip for Maxwell's equations in Eq. ( 92) have the following form: Note that C ij has stronger singularity than its counterparts in mechanics.Also noteworthy is the fact that the curl of C ij is less singular than C ij and that the double curl of C ij is proportional to C ij modulo Dirac delta.We have the following integral representations for E and H: In these formulae, j and m are the surface electric and magnetic current vectors defined by jðyÞ ¼ nðyÞ Â HðyÞ; mðyÞ ¼ EðyÞ Â nðyÞ where n stands for the outward unit normal vector on the surface of the domain D, and div S indicates the surface divergence.From Eqs. ( 93) and ( 94), we have where m inc and j inc are defined in an obvious manner.Equations ( 95) and (96) where j 0 is the basis function used for the discretization for j.
Galerkin approaches based on Eq. ( 97) are usually called method of moments (MoM) in the EM community.We note that the basis functions for j have to satisfy the continuity of normal components across the edges of elements, or "divergence conforming," for the use of the integration by parts to be justifiable.RWG functions [357] are among the most popular choices of divergence conforming basis functions used in MoM.These functions are vector valued and have degrees of freedom on edges of elements.Use of RWG functions for both of j and m, however, may not be a good choice for problems with dielectric scatterers since it may lead to ill-conditioned systems.Use of dual [358] (or Buffa-Christiansen [359]) basis function is a possible solution to this problem.

Fast
Methods for Maxwell's Equations.Fast methods for equations of the Helmholtz type seem to be investigated most vigorously in the EM community because of many important applications such as radars, etc. Examples of very large problems solved with fast solvers of integral equations are found in papers such as Taboada et al. [360] (150 Â 10 6 unknowns), Ergu ¨l and Gu ¨rel [361] (more than 200 Â 10 6 unknowns), etc.Not surprisingly, there are overwhelming amount of publications devoted to the subject of fast integral equation solvers which make it impossible to cover even just the main developments in a review article of this size.The reader may find it more informative to refer to a book [356] or to introductory or review articles by specialists [362][363][364][365]. 4.7.3Periodic Problems in Maxwell's Equations.As a recent progress in computational electromagnetics, we can mention periodic FMM.Periodic boundary value problems for Maxwell's equations are increasingly important these days because there are new micro and nano devices which utilize periodic structures, such as photonic crystals and metamaterials.Photonic crystals are a class of periodic structures of dielectric materials having a periodicity comparable to the wavelength of light [366].Photonic crystals are known for their peculiar properties such as frequency selectivity, etc., and are said to guide and store the light quite freely.Metamaterials are a class of composite materials with subwavelength periodic metallic structures.They exhibit curious behaviors such as negative apparent refractive indices for a certain frequency range.Applications of metamaterials are being sought in optical fibers, in the so called "super lens" [367] with extremely high resolution, etc.
BEMs for periodic boundary value problems are obtained easily as one replaces fundamental solutions by periodic Green's functions in the BIEs.Periodic Green's functions for wave problems usually take the form of lattice sums of fundamental solutions.This means that a solution of a periodic boundary value problem can be viewed as an infinite repetition of the solution in the unit cell.This structure of the solution matches well with the tree structure of FMM.Indeed, one may identify the unit cell with the level 0 box of FMM to formulate a periodic version of FMM.The effects from far replicas of the unit cell are then computed as an infinite sum of M2L formulae between the unit cell and its replicas.An attempt of this type for 3D doubly periodic problems in Maxwell's equations is found in Otani and Nishimura [368] where they considered the case of a cubic unit cell.This approach has been extended to the case of a parallelepiped unit cell in [369] and to a tall unit cell case in [370].Figure 19 shows an example from Ref. [370] where layers of glass fibers in water are considered.This is supposed to be a model of the skin of worms, which is known to look blue if viewed from a certain angle: an example of the so-called structural color.
Figure 19(b) shows the energy transmittance when the wave length of the incident wave is 475 nm (blue).This result shows that one will observe a strong reflection of blue light when the incident angle is 30, which is in agreement with experiments.Periodic problems are characterized by anomalous behaviors of solutions related to the existence of guided modes.The performance of the periodic FMM for Maxwell's equations near anomalies are examined in Ref. [371].
It is also possible to formulate FMMs for periodic boundary value problems without using periodic Green's functions or lattice sums.Barnett and Greengard [372] presented well-conditioned integral equations with layer potentials distributed not just on the surface of the scatterers but on the boundary of the unit cell as well.Their efforts of this type in 2D Helmholtz equation are found in Ref. [372] for doubly periodic problems and in Ref. [373] for singly periodic problems.They also proposed a novel approach to compute lattice sums accurately [372].4.7.4Periodic Problems in Mechanics.Periodic problems are found in mechanics also.Actually, the seminal paper by Greengard and Rokhlin [171] already includes discussions on periodic boundary conditions in many-body Coulomb systems.Also, periodic boundary value problems appear in microscopic analyses in the homogenization theory which is used in estimating macroscopic moduli of composite materials.Houzaki et al. [374] and Otani and Nishimura [375] considered periodic boundary value problems for rigid inclusion problems in 2D antiplane problems and 3D problems of elastostatics.We note that static periodic problems are delicate mathematically because the periodic Green's functions do not exist.However, their numerical treatment is easier than the corresponding dynamic problems once the mathematical problem is resolved.
In regard to dynamic problems, we mention Otani and Nishimura [376] for singly periodic crack problems in 2D for Helmholtz equation and Isakari et al. [339,377] for elastodynamic scattering problems in 3D by doubly periodic inclusions.
Different to the integral equation in Eq. ( 86) is the convolution integral in time and the fundamental solutions are those in Eqs.(84) and (85) transformed to time domain.Fortunately, for most physical problems, the singular behavior of the time-dependent fundamental solutions with respect to space is the same as in the respective static counterparts.A second hypersingular BIE as counterpart to Eq. ( 87) can be as well deduced by applying the gradient operator on the first integral equation.Details from a mathematical point of view can be found in [378], where an engineering representation may be found in [41,379].
To obtain a boundary element formulation, the domain as well as the time (T ¼ MDt) is discretized and the shape functions in space and time have to be chosen.For the acoustic example this is with the nodal values in space and time / j i and q j i .Note, the summation over i is the spatial counter, whereas the summation over j is the counter in time.This approach is different from timedependent (standard) finite element approaches where the time discretization is applied on the PDE.Such approaches applied in BEM would result in domain integrals.The last step is the integration in space and time which yields the following matrix formula: In Eq. ( 100), k m and b m have the same meaning as in Eq. ( 13), however, with changing values in each time step.The convolution integral produces the summation over the past time steps.Hence, a standard implementation has a complexity of O(N 2 M) for M An efficient implementation to reduce this effort can be found in Ref. [380].
An alternative approach to the above sketched procedure is the transformation to Laplace or Fourier domain and to solve the elliptic problem.This avoids the high storage but requires a numerical inverse transformation with all its difficulties (see, e.g., Refs.[381,382]).Furthermore, at each frequency the equation system must be solved in contrast to Eq. ( 100) where only one equation solution is necessary.An overview on dynamic BEM formulation until the year 1996 can be found in Refs.[383,384].
4.8.1 Analytical Integration in Time.The integration in time is mostly performed analytically.This is due to the structure of several time domain fundamental solutions used up to now.In acoustics the fundamental solution is a Dirac distribution in time and, hence, the time integral is reduced to the function at the retarded time t À r=c.In elastodynamics, the integration can as well be done analytically, but in this case first the shape functions (Eq.( 99)) have to be chosen (mostly linear for the displacements and constant for the tractions) and then the integration is performed within each time step.This procedure has been published initially by Mansur [385].In three dimensions, due to Huygens principle, there exists a cut-off time at which the system matrices in Eq. ( 100) becomes zero.However, this is not true in two dimensions and for any problem with damping.4.8.2Convolution Quadrature Method.An alternative methodology which also extends the time domain formulations to problems where the time domain fundamental solution is not known is the so-called convolution quadrature method (CQM)-based BEM.In this formulation, the convolution integral is approximated numerically by the formula ð t 0 ð suppðu p j Þ Gðx; y; t À sÞ q j ðsÞ u p j ðyÞ dSðyÞds % X n k¼0 x nÀk j ð Ĝ; DtÞ q j ðkDtÞ with The CQM has been developed by Lubich [386] and is applied to elastodynamics in Ref. [387].In Eq. ( 101), R is a parameter to be chosen (often R ¼ 10 À 5 2M ) and c(z) denotes an A-and L-stable multistep method.The second-order backward differential formula is a good choice [388] and A-and L-stable Runge-Kutta methods can be used [389].The essential advantage of this formula is that both the time dependent fundamental solution and its Laplace transform Ĝðx; y; sÞ can be used.For many timedependent problems, especially if damping is included, only the transformed fundamental solutions are available.The applications to visco-and poroelasticity can be found in Refs.[390,391].Beside this advantage, the stability of the algorithm is improved compared to the analytical integration.Finally, the same recursion formula (Eq.( 100)) is obtained where the matrix entries are now the integration weights x n j which present the time-integrated fundamental solutions.
As an example, the computed sound pressure level on the boundary of an amphitheatre using the time-domain BEM are depicted in Fig. 20 at the times t ¼ 0.00625 s, t ¼ 0.046875 s, and t ¼ 0.078125 s.The amphitheatre is discretized with 3474 linear triangles with 1761 nodes and a time step size of Dt ¼ 0.003125 s is used.Sound reflecting boundary conditions are considered and the sound propagates in the air with a speed of c ¼ 346m=s.The sound field is excited by a prescribed sound flux of q ¼ 1 N=m 2 in the middle of the stage.Looking on the lowest sound pressure level, it is observed that there are high negative values.Recalling the definition of the sound pressure level, those values denote sound below the auditory threshold.At all areas which cannot be reached by the sound wave at the observation time such low values are obvious.Also, after some time nothing more can be heard.

Numerical
Instabilities and Weak Formulations.The classical space-collocation approach based on Eq. ( 98) suffers, in some cases, from the numerical instability of the time-marching scheme which has been put into evidence, e.g., in Refs.[392][393][394][395].This instability is often called intermittent owing to the often unpredictable way in which it appears upon variation of the time step adopted in the analysis.
It has been shown that weak formulations can be very beneficial with respect to this issue.A first proposal in the time domain was put forward and implemented by Ha Duong [396] for acoustics; it couples the weak version -both in time and space -of the integral representation Eq. ( 98) with the weak gradient BIE which is obtained from Eq. ( 98) by applying the gradient operator with respect to x.Recently, Aimi et al. [397,398] have proposed an alternative direct space-time Galerkin BEM which leads to unconditionally stable schemes.Instead of exploiting the gradient BIE, this latter formulation makes use of the velocity BIE obtained by applying to Eq. ( 98) the time derivative with respect to t. 4.8.4Fast Methodologies.One of the first approaches for a fast evaluation of the retarded potentials related to scattering problems was published in Ergin et al. [342].The main idea is a plane wave expansion of the respective integral kernel.This idea has been extended to elastodynamics in Refs.[341,343].Another approach has been proposed by Hackbusch [399] for acoustics using the panel clustering technique to sparsify each matrix per time step.An improvement of this approach can be found in Ref. [400].For the heat equation, an approach based on a Chebyshev interpolation of the kernel has been proposed by Tausch [401].An FMM approach in viscoelasticity has been proposed in Zhu et al. [402].Beside the above-mentioned approaches, the reformulation of the CQM by Banjai and Sauter [403] allows the application of fast techniques developed for elliptic problems.In Ref. [404] this reformulated CQM is discussed as well for collocation approaches.In combination with ACA, a fast elastodynamic BEM has been proposed by Messner and Schanz [210].The CQM in combination with an FMM has been studied by Saitoh et al. [405].
Finally, we mention related but different use of fast BEM in time-domain problems in fluid mechanics in which the governing equations are integro-differential equations.Solution methods for such equations require evaluations of integrals, which are accelerated with fast methods.Examples of approaches of this type are found in vortex methods (e.g., Ref. [406]) and in vesicular flow problems (e.g., Ref. [407]), to mention just a few.
4.9 Coupling of the BEM with the FEM.Even if the finite element method has a dominant status in the field of computational methods in engineering, mostly because of its great flexibility and wide range of applicability, approaches based on integral equations are clearly superior for certain classes of problems.In fact, the FEM and the BEM have relative benefits and limitations depending on the problems to be solved.Hence, there are many applications where a coupling between the two approaches is in principle very attractive.For example, different parts of a structure can be modeled independently by the FEM and the BEM in order to exploit the advantages of both techniques.4.9.1 Methodologies.The first investigations on this idea dated back to 1977 with the pioneering work of Zienkiewicz, Kelly and Bettes [408] based on the collocation non-symmetric BEM.A symmetric coupling of the BEM and FEM was proposed later by Costabel and Stephan [409] followed by the variationally based coupling procedure between the FEM and an indirect version of the Galerkin BEM developed by Polizzotto and Zito [410].In many papers, the BEM subdomains, see for instance the work by Haas and Kuhn [411] and by Ganguly, Layton and Balakrishna [412], or the FEM subdomains, e.g., the work by Ganguly, Layton and Balakrishna [413], are handled as equivalent finite or boundary macro-elements, respectively.In recent years, however, the idea of adapting one of the two techniques to the other in a sort of master-slave technique with macroelements has been gradually abandoned for a more balanced approach.
The need to maintain the different subdomains as independent as possible has become a crucial point for different reasons.First, this aims at preserving the structure of the system matrices which are sparse positive definite for the FEM and full, possibly nonsymmetric and not definite for the BEM.Second, the computational cost and memory requirements of the BEM part must be imperatively reduced by means of algorithms such as fast multipole methods, panel clustering and ACA; this entails that the system matrix is actually never available.Third, different subdomains might require different discretizations (both in time and space) resulting in non-conforming meshes and interpolations.This last point is closely related to the enforcement of the matching conditions (e.g., continuity of displacements and tractions in solid mechanics) at the interface between the subdomains.When displacement continuity is enforced in strong form, the BEs at the interface and the sides of the FEs adjacent to the interface should coincide, together with their nodes.On the other hand, if the kinematic continuity condition is enforced in weak form, the matching of the two discretizations at the interface is no longer necessary, with an evident flexibility gained in the meshing process.This permits the integration of independently discretized models possibly coming from different sources.
A non-matching node coupling of SGBEM and FEM is developed by Ganguly, Layton and Balakrishna [413] in the context of a BE-based approach.In fracture mechanics, a coupling approach which preserves the independence of the FE and BE meshes is the so-called hybrid surface-integral-finite-element technique, see the works by Keat, Annigeri and Cleary [414], Han and Atluri [415] and Forth and Staroselsky [280].Springhetti et al. [416] address a SGBEM-FEM coupling procedure in elastostatics characterized by a weak enforcement of the both matching conditions at the interface.
Two different families of approaches seem to guarantee, at best, the complete independence of the various subdomains: iterative (or alternating) techniques and Lagrange multiplier techniques.The former type of approach could be in principle applied to any kind of coupling.The individual subdomains are treated independently by either method based on an initial guess of the interface unknowns.Then, the newly computed interface variables are synchronized and exploiting these updated values in another subdomain solution yields enhanced results.The procedure is repeated until convergence is achieved.A comprehensive overview of such methods is given by von Estorff and Hagen [417].Within the iteration procedure, a relaxation operator is applied to the interface boundary conditions in order to enable and speed up convergence.In this sense, the iterative coupling approaches are better called interface relaxation FEM-BEM.Although the iterative coupling is very attractive from the point view of software design, the convergence commonly depends on the relaxation parameters which are rather empirical (see Refs. [418][419][420][421].
In the Lagrange multiplier approach matrix entries are never mixed from the different subdomains and the interface conditions are posed in a weighted sense such that nonconforming interfaces can be handled.The largest family in this category, represented by the FETI-BETI approaches [422,423], can be formulated as follows.At the subdomain level, Dirichlet-to-Neumann maps are realized by means of each discretization method.In elastostatics, for instance, this amounts to generating a Steklov Poincare ´operator expressing surface tractions in terms of surface displacements t ¼ Su.The realization of this operator for FEM and BEM is discussed at length, e.g., in Refs.[422][423][424].Then, tractions along the interface are treated as Lagrange multipliers k and the equilibrium and continuity conditions are enforced in a weak manner along the interfaces in terms of u and k.Interface displacements and Lagrange multipliers are the unknowns of the reduced problem which can be solved by means, e.g., of a conjugate gradient method [425].The overall approach can be interpreted as a preconditioned conjugate gradient solver [422].
As far as evolving problems are concerned, one additional reason for avoiding the solution of the global coupled system of equations is that the use of different time steps for the subdomains becomes possible.This is an important advantage, especially in the BEM, where the range of applicable time step lengths, resulting in a stable and accurate solution, is limited.In the staggered solution approach, the equations for each subdomain are solved once at each time step.However, the stability and accuracy of both the BEM itself and the staggered coupling approach impose requirements on the choice of the time step durations, which may be contradictory.Therefore, it is desirable to use a coupling procedure that is stable and accurate for a wide range of time steps.This can be obtained by introducing corrective iterations into the staggered algorithm [421].Within the iteration procedure, a relaxation operator may be applied to the interface boundary conditions in order to speed up convergence.An overview of a variety of 031001-26 / Vol.64, MAY 2011 Transactions of the ASME such techniques involving FEM-BEM as well as BEM-BEM coupling, is given by Refs.[417][418][419].
4.9.2Applications.A rich literature on applications of coupled FEM-BEM approaches has flourished in recent years combining and exploiting at best the potentials of both the FEM and BEM techniques.The FEM is the standard method of choice for the analysis of complex structural parts and, more in general, of structural mechanics in the presence of shells, plates and beams.It allows for an easy handling of general constitutive laws (encompassing nonlinearities and anisotropy) and is widely used in the industry.However, the BEM is clearly superior in the presence of domains of infinite extension featuring linear isotropic constitutive behaviors, evolving domains, strong gradients and singularities.When both situations occur in a single analysis, then a coupled approach turns out to be very promising.
A typical application with these requirements is represented by 3D linear fracture mechanics, especially in the presence of propagating fractures in complex structures.The FEM is inefficient in the modeling of singular stress fields (or regions with strong gradients) and in the presence of evolving domains.Such difficulties have been elegantly circumvented by a coupling procedure (with the symmetric Galerkin BEM in particular), e.g., by Frangi and Novati [426], Forth and Staroselsky [280], Aour, et al. [427], and Lucht [428].
In elastoplastic analyses with a limited spread of plastic deformations, the FEM is utilized where the plastic material behavior is expected to develop.The remaining bounded/unbounded linear elastic regions are best approximated by the BEM [409].A central aspect of coupling approaches is that they require the user to predefine and manually localize the FEM and BEM sub-domains prior to the analysis.Recently, Elleithy [429], Elleithy and Grzhibovskis [430] focused on the development of an adaptive FEM--BEM coupling algorithm in which they provide fast and helpful estimation of the FEM and BEM sub-domains.
Soil structure interaction in statics and dynamics has stimulated intensive research in coupled approaches where the BEM is employed for the infinite (or semi-infinite) and possibly layered soil, while the FEM is adopted for the finite structure and structural elements like piles [431][432][433][434][435][436].Recently, applications in vibration analysis of transport systems have been addressed by similar techniques [437][438][439].
In magnetostatics industrial applications, the BEM-FEM coupling is often applied to the analysis of finite metallic parts of a mechanism (like a relay) excited by coils.The surrounding air is considered as an isotropic and linear medium addressed by the BEM with total or reduced scalar potential approaches.The nonlinear constitutive behavior of the metallic parts and their high surface-to-volume ratio call, on the contrary, for the application of the FEM.This approach also allows performing coupled analyses encompassing large displacements and rotations of the pivoting lamina.Investigations along these guidelines have been presented, e.g., by Kuhn [440], Balac and Caloz [441], Kuhn and Steinbach [442], Frangi et al. [443,444], Salgado and Selgas [445], Pusch and Ostrowski [446].The simulation of eddy currents can also be conveniently addressed with similar approaches [447].
BEM-FEM techniques have a dominant status in the simulation of elasto-acoustic coupling.While the FEM is the preferred tool in the field of structural vibrations, using the BEM for the unbounded linear domain has the additional benefit that the Sommerfeld radiation condition for exterior domains is inherently fulfilled.Advanced applications have been presented by Chen et al. [448], Czygan and von Estorff [449], Gaul and Wenzel [450], Langer and Antes [451], Fischer and Gaul [425], and Soares [452].
Although the above categories have attracted most of the research effort, fluid structure interaction is also emerging as promising.For example, Wang et al. [453] focused on the simulation of cell motion in a biological fluid flow.The membrane of a moving cell is represented by a thin shell composed of incompressible neo-Hookean elastic materials and the movement of liquids around the membrane are approximated as incompressible Newtonian flows with low Reynolds numbers.Seghir et al. [454] present a numerical model coupling boundary and finite elements suitable for dynamic dam-reservoir interaction employing a standard finite element idealization of the dam structure and a symmetric boundary element formulation for the unbounded reservoir domain.Other important applications include diffuse optical tomography by Elisee et al. [455] and modeling of photonic crystals for semiconductor laser beams by Jerez-Hanckes et al. [456].

Future Directions
In this closing section, we discuss some future directions of the BEM and related research topics.The views expressed herein are not meant to be exclusive.They are merely those of the authors of this review article, aiming to shed some light on the research efforts to be made in the coming years.All these efforts should further improve the formulations, implementations, solutions, and applications of the BEM to make it an even more effective and efficient numerical method in the fields of computational science and engineering.
5.1 Multi-scale, Multi-physics, and Nonlinear Problems.Multi-scale and multi-physics applications of the BEM should be a productive area of research for the BEM in the coming years.Many of the multi-physics and multi-scale problems involve complicated geometries and large-scale numerical models, such as modeling of nanomaterials [250], biomaterials, and MEMS, that can be solved efficiently now by the BEM.Any new development of the BEM to solve realistic models in these areas will further demonstrate the unique capabilities of the BEM and benefit industry.
The authors believe that the BEM is ready to solve some nonlinear problems on a larger scale as compared with the conventional BEM used in the past.Because of the new formulations (symmetric Galerkin, various preconditioning techniques, etc.) which offer more stable solutions and fast solution methods (FMM, pFFT, ACA, etc.) with iterative solvers for large-scale problems, the BEM and the related methods should be able to handle effectively nonlinear problems with elastic-plastic materials [35], large deformation [39,457,458], and moving-boundary problems [38] such as in hydrodynamics [459].
5.2 More Efficient Fast Solution Methods.The BEM accelerated with the FMM, pFFT, or ACA can now solve largescale BEM models of several million DOFs very efficiently, if they work.In some cases, the convergence of the FMM, pFFT, or ACA with iterative solvers may not be achieved within reasonable time limits.Therefore, developing better pre-conditioners for the BEM systems of equations has been, and will continue to be, a focus for improving the fast multipole, pFFT, and ACA BEMs.Wavelet methods [460,461] can also be used to compress the BEM matrices to yield better conditioned systems of equations for solving large-scale BEM models.To make kernel-independent fast multipole BEM, the black box, or generalized fast multipole method (bbFMM) [462] can be further explored.

Applications in Urgent Engineering
Fields.The BEM is uniquely positioned to address some of the urgent problems in biomedical engineering, alternative energy and environmental engineering.For example, the computer image-based BEM can be applied to study thermal and electrical responses of bones or brains [44,463], which are difficult to model using other methods due to the complicated domains.The fluid-structure interaction problem of red-blood cells in Stokes flow can also be attempted by the BEM [44,464].This is currently an active research topic.
We believe that the BEM can play a significant role in the development of the next generation, more efficient and environmentfriendly energy harvesting systems in the next decade to tackle urgent energy-related problems.For example, the BEM can be applied to predict thermal, electrical, fluid flow, and other mechanical responses, fatigue life and noise levels of arrays of solar panels and wind turbines, which are often in open spaces and in large scales.With full-scale and large-scale models of solar panels and wind turbines, we will be able to study their interactions with the environment and evaluate their impacts on human beings and other species.
In other areas of environmental engineering, the BEM can be applied to predict noise along a highway stretch caused by cars or in a city block caused by high-speed trains or explosives, to track the distribution and propagation of pollutants over a city, or to evaluate manmade noise in the ocean and the negative impacts on marine lives.

Software Development.
The transition of research codes to commercial software for the BEM has been lagging for the last three decades as compared with those for other numerical methods.This is not easy for researchers who are not familiar with the marketplace.However, for the BEM to be useful and to be relevant in computational science and engineering, researchers on the BEM have to be aware of this issue.BEM researchers should solve large-scale, realistic, industry relevant problems with their developed BEM approaches and codes, moving beyond the academic problems they are used to solve.They should interact more often with engineers to be aware of the type of the problems that they can solve with the BEM and thus contribute to industry research and development efforts.
In fact, there are many advantages of the BEM regarding the development of a software package, as compared with the FEM.The data structure for a BEM model is simpler and requires much less storage space due to the use of boundary elements.The required geometric (CAD) model of a problem can also be less stringent.For example, only the surface data are required in the BEM, and that surface data need not to be perfect (or "air tight").Meshing surfaces for the BEM are straightforward and there are many codes for meshing available in the public domain.Therefore, optimizing and integrating a BEM research code with other software components to develop a BEM software package is relatively straightforward.
There have been some efforts by researchers in the BEM to develop commercial BEM software for industrial applications.Among the commercial BEM software available, we mention the packages GPBEST by BEST Corporation, BEASY TM by Computational Mechanics Inc., COMET BEAT by Comet Technology Corporation, Virtual.Lab by LMS International, and FastBEM AcousticsV R by Advanced CAE Research.A list of free educational computer codes developed at several universities for solving various problems using the BEM can also be found at the Boundary Element Resources Network (BENET) [465].BEM source codes for education can also be found in textbooks, such as the code for Galerkin BEM in Ref. [43] and the code for fast multipole BEM in Ref. [44].
More coordination among the BEM researchers and collaborations with industries are urgently needed in order to develop the next generation BEM software that will be unique, efficient, and user-friendly, and to add more valuable computational tools in the computational toolbox for engineers.Establishing an open-source platform for the BEM, where researchers and developers can contribute and share their BEM source codes and work together toward an integrated BEM software system, may be a first step towards achieving these goals.
5.5 Education.Education of a new generation of students who are interested in conducting research and development of the BEM will be crucial for further advances in the BEM.Educators need to consider more innovative ways to teach the subjects of the BEM that have sometimes been abstract and theoretical in the past.Developing new course materials on the BEM, that include the BIE theory, element formulations, programming, and realworld applications, is needed.Integration of the BEM with other computational methods in a course on computational science and engineering may be a better way to introduce the BEM to students.
There have been several conference series dedicated to the BEM in the last decade, such as the International Conference on Boundary Elements and other Mesh Reduction Methods (BEM/ MRM), International Conference on Boundary Element Techniques (BETeq), International Association for Boundary Element Methods (IABEM) Conference, and numerous minisymposia at the World Congress on Computational Mechanics and US National Congress on Computational Mechanics.In addition, two workshops on the fast multipole BEM were conducted, one in Los Angeles in 2006 and one in Kyoto in 2007.There were also two GF related workshops earlier in the 1990s.The first workshop was on "Green's functions and boundary element analysis for modeling the mechanical behavior of advanced materials" which was held in 1994 in Boulder, Colorado and sponsored by the National Institute of Standards and Technology (NIST).In 1998, another workshop titled "Library of Green's Functions and Its Industrial Applications" was held at Gaithersburg, Maryland.This workshop, also sponsored by NIST, was to demonstrate the idea and application of a library of discretized GFs to industrial problems.There were two GF related educational projects supported by NSF: One was the GREEN Project (Green's Functions Research & Education Enhancement Network), the other was an extension of the GREEN project "NSF-NSDL GREEN Project: A Digital Library Partnership of Academia, Government, and Industry." Nationwide and worldwide annual or biannual conferences, workshops, and symposia on the BEM and related methods should continue so that more students and young researchers can have a chance to learn the BEM, develop interests in the BEM and other computational methods, and eventually help to advance BEM research and development in the future.

Fig. 1 A
Fig. 1 A 3D domain V with boundary S

Fig. 2
Fig. 2 Illustration of the Galerkin weight functions for 2D BEM

Fig. 3
Fig. 3 Domain of dependence and range of influence.(a) The nodes 1, 2 and 3 lie within the DOD of the evaluation point E. The ROIs of nodes 1, 2, 3, 4 and 5 are shown as gray regions.In the standard BNM, the ROI of a node near an edge, e.g., node 4, is truncated at the edges of a panel.In the EBNM, the ROI can reach over to neighboring panels and contain edges and=or corners -see, e.g., node 5 (b) Gaussian weight function defined on the ROI of a node (from Ref.[149]).

Fig. 7 (
Fig. 7 (a) 2D view of a parallelepiped superimposed onto a discretized sphere.The black dots indicate grid points and the cubes are represented by the solid lines; (b) Illustration of the four steps in the pFFT scheme; the gray region denotes the near-field region of the yellow element.

4. 1
Modeling of MEMS/NEMS 4.1.1Electrostatically Driven MEMS/NEMS Devices.Numerical solutions of electrically actuated microelectromechanical systems (MEMS) have been carried out for nearly 20 years by using the BEM to model the exterior electric field and the finite element method or the BEM to model deformation of the structure.The commercial software package MEMCAD [211], for example, uses

Fig. 8 H
Fig. 8 H-matrix with 3279 entries and rank numbers of ACA

Fig. 10 SEM
Fig. 10 SEM pictures of (a) a laterally oscillating beam resonator, and (b) a biaxial accelerometer

Fig. 16 A
Fig. 16 A BEM model of the Skipjack submarine impinged upon by an incident wave

Fig. 18 BEM
Fig. 18 BEM model of 1152 spherical cavities in an infinite elastic solid [343].The spatial DOF is 1,105,290 and the number of time steps is 200.The CPU time is 10 h 47 min and the memory requirement is 152.8GB.

4. 8
Time-Domain Problems.The BEM in the time domain is based on the representation formula of the underlying timedependent partial differential equation.A BIE in time domain can be deduced either with partial integration or based on reciprocal theorems.In the case of acoustic problems, the acoustic pressure / x; t ð Þ can be written as cðxÞ/ x; t ð y; t À sÞ qðy; sÞ À Fðx; y; t À sÞ /ðy; sÞ ½ Â dSðyÞds 8ðx; tÞ 2 S Â 0; T ½

Fig. 20
Fig. 20 Sound pressure Level (dB) at the boundary of the amphitheatre
yÞ are the two kernel functions of inclusion a.They are functions of the shear modulus, Poisson's ratio, and outward normal for inclusion a.
1 2 u i ðxÞ ¼ ð S U ij ðx; yÞt j ðyÞ À T ij ðx; yÞu j ðyÞ Â Ã dSðyÞ; 8x 2 S (71) where u i and t i are the displacement and traction, respectively; S ¼ [ a¼0 n S a is the total boundary of domain V 0 (assuming that S is smooth around x), and U ij ðx; yÞ and T ij ðx; yÞ are the two kernel functions (Kelvin's solution).For each inclusion, the BIE can be written as 1 2 u ðaÞ i ðxÞ ¼ ð Sa U ðaÞ ij ðx;yÞt ðaÞ j ðyÞ ÀT ðaÞ ij ðx;yÞu ðaÞ j ðyÞ h i dSðyÞ; 8x 2 S a (72) for a ¼ 1; 2; :::; n, in which u ðaÞ ij ðx; yÞ and T ðaÞ ij ðx; 64, MAY 2011 Transactions of the ASME Downloaded 21 Apr 2012 to 131.175.24.250.Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm periodic boundary value problems in mechanics and the BEM to solve them.4.7.1 BEM and MoM for Maxwell's Equations.We consider an exterior domain D 2 R 3 in which Maxwell's equations are satisfied: are called EFIE (electric field integral equation) and MFIE (magnetic field integral equation).Conceptually, EFIE (MFIE) corresponds to the displacement (traction) BIE in mechanics or vice versa.Note, however, that EFIE and MFIE have essentially the same structure, which is in contrast to BIEs in mechanics in which hypersingularity appears only in traction BIEs.Boundary value problems for Maxwell's equations are solved as one solves either EFIE or MFIE for unknown boundary densities j and/or m.These density functions have two independent components in 3D problems since they are tangential to the boundary, while their mechanical counterparts have three independent components.A typical boundary condition is m ¼ 0 on S for perfect electric conducting scatterers.As in mechanics, solutions of EFIE or MFIE may not be unique for particular values of x called irregular frequencies.In that case, one may consider a