Marc Bonnet

Giulio Maier

Castrenze Polizzotto

Symmetric Galerkin boundary element method

Keywords:

published or not. The documents may come

INTRODUCTION

The boundary integral equation (BIE) approaches referred to in this article as 'traditional' stem from the mathematical work on integral equations (more specifically on Somigliana's integral representation formulas established in 1886), as formulations of linear boundary value problems alternative to those in terms of partial differential equations. Traditional boundary element methods (BEMs) developed from these BIE formulations in the late 60s and in the 70s, through modelling of boundary variables and, mostly, through collocation as their algebraized version for approximate numerical solutions. In the last two decades, thus originated BEMs have been the subject of a very considerable effort in computational mechanics, gaining a still growing popularity, competing with finite element methods (FEMs) in some areas and spreading in engineering applications, also by means of commercial computer codes.

Despite their undeniable success, besides remarkable and now well understood advantages in some kinds of problems, traditional BEMs in the above sense are known to exhibit certain unpleasant Transmitted by Associate Editor Xxx Yyy features, among which most pertinent here is the lack of symmetry in some matrix operators which play a key role in various theoretical developments and analysis procedures. These nonsymmetric operators typically include: the matrix which in linear problems relates the boundary unknowns to the vector containing the data; the matrix which transforms the vector collecting the modelled plastic strains into the consequent selfstress vector; the stiffness matrix which relates kinematic boundary variables to static ones, with subsequent well-expected difficulties in FE-BE coupling.

Symmetry can be understood as a synonym for harmony. In fact, ethymologically, σνµµετ ρια meant commensuration and was referred to musical notes in Pythagoras' canons. However, lack of symmetry in traditional BEMs entails undesirable effects not only of aesthetical nature. It is accompanied also by lack of sign-definiteness (or semi-definiteness) and the negative consequences are far-reaching and have both theoretical and computational significance. In fact, typically, e.g., continuous solutions of the BIEs (or of BIE-based formulations for nonlinear problems) and their approximations as discrete solutions of BE models fail to be characterized by variational properties; criteria for convergence and algorithmic stability in time-stepping solution techniques, known in FE elastic and inelastic analyses, cannot be proved for their traditional BE counterparts; as a third example, important concepts of structural plasticity such as shakedown and bounding theorems do not carry over to systems discretized in the traditional BEM fashion.

It was said that 'advances in engineering and mathematical sciences usually spring from a feeling of dissatisfaction, much as the oyster responds to the grains of sand'. The symmetric Galerkin boundary element methods (SGBEM) of concern in this review article can be regarded as a response to grains of sand contained in traditional BEM and pointed out above. A number of responses to those grains , i.e. remedies to the above undesirable features of traditional BIE-BEMs, can be found in the recent literature in terms of BIE-BE formulations endowed with symmetries. Only the SGBEM will be focused herein. However, it represents the symmetrization provision most frequently dealt with in the present literature and, in the writers' opinion, is the most mature and promising one as for future prospects from both the theoretical and computational standpoint.

The SGBEM can be singled out by means of a few of its distinct peculiarities as follows: its continuous formulation rests on BIEs generated with recourse to both single-layer and double-layer sources, in such a way that the integral operator turns out to be symmetric with respect to some suitably constructed bilinear form; the BE-discrete model (BEs and domain cells simultaneously in some problems) is constructed by weighted residuals with shape and weight functions correlated in the sense of Galerkin and so that the governing variables turn out to be 'generalized variables' according to a concept originally introduced by W. Prager in structural plasticity, i.e. such that the dot products and their energy meanings are preserved in the transition from fields of (local) variables to vectors of (global, nodal) variables in the discrete model.

In the early literature, the symmetric formulation for BIEs as understood here was first proposed for linear elastic analysis in a 1979 paper by Sirtori Since the late eighties the literature on SGBEM, to be surveyed in this paper, has grown and still grows at an accelerated pace. Two recent comprehensive books on BEM in general devote chapters to the fundamentals of the SGBEM: the former by Kane appeared in 1994, Ref. [START_REF] Kane | Boundary element analysis in engineering continuum mechanics[END_REF], the latter by Bonnet in 1995, Ref. [START_REF] Bonnet | Equations intégrales et éléments de frontière[END_REF], its English translation being expected in 1997, Ref. [START_REF] Bonnet | Boundary integral equations methods[END_REF]. Some aspects of the SGBEM have been dealt with in review papers either ad hoc, ref. [START_REF] Maier | Symmetric Galerkin boundary element method in plasticity and gradientplasticity[END_REF], or covering subareas of BEMs, such as the recent AMR article by Beskos on dynamics, Ref. [START_REF] Beskos | Boundary element methods in dynamic analysis, part. II (1986-1996[END_REF].

The present article is believed to provide the first comprehensive systematic survey of the SGBEM. The research results achieved up to the end of 1996, and available in the general literature, will be reviewed in twelve sections to follow. The symmetric continuum formulations of BIEs and their Galerkin discretizations (with generalized variables in Prager's sense) by BEs are outlined for the categories of problems so far considered in the SGBEM context, namely: linear elastostatic (and potential) problems in Sec. 2; elastic-plastic analysis in Sec. 3; fracture mechanics (linear elastic and quasi-brittle) in Sec. 4; time-dependent linear problems in Sec.5 (primarily linear elastodynamics, briefly also transient heat conduction and its analogues, viscoelasticity, poroelasticity and acoustics). Section 6 is devoted to such diverse subjects as limit analysis and analysis of Kirchhoff plates by the SGBEM. Section 7 presents energy approaches and variational theorems with reference primarily to linear elasticity and elastoplasticity. The integrations are dealt with in Sec. 8 , with special attention to hypersingular integrands which represent a crucial aspect in SGBEM and has been a braking factor in its development. Other mathematical and computational issues are presented in Sec. 9, specifically: convergence, algorithmic stability, error estimates, papproaches, self-adaptive solution techniques.

Sensitivity analysis has recently become an effective tool for solving inverse problems such as parameter identification and structural shape optimization. Therefore it is briefly discussed in Sec. 10 in view of the remarkable potentialities of the SGBEM in this area.

Sec. 11 concerns coupling between SGBEM and FEM, since this may represent an attractive novelty in domain decomposition, which is another fashionable issue in today's computational mechanics.

Section 12 surveys the still inadequate but fast growing computer implementations of the SGBEM and attempts some comparative assessments of computational merits. Some of the diverse, alternative or related, symmetric formulations that appeared in the literature are briefly surveyed in Sec. [START_REF] Banerjee | The boundary element method in engineering[END_REF], with reference to representative contributions. Section 14 gathers closing remarks on the consolidated results, limitations, potentialities and research needs of the SGBEM, as emerging from the present state-of-the-art study. In assembling the reference list, completeness in terms of meaningful contributions has been pursued together with objectivity of selection criteria; however, undue omissions are regrettably inevitable even in a restricted subject like SGBEM and the authors apologize for them.

NOTATION. The formalism and nomenclature of matrix algebra are adopted in this paper, except where special purposes make indicial tensor notation more convenient. Matrices and vectors are represented by bold-face characters; exponent T marks transpose; vector inequalities apply componentwise. The most frequently used symbols are defined here: x, ξ are coordinate vectors in a Cartesian orthonormal reference system, for the field point and source point (or load point), respectively; u, p, b are displacements, tractions and body forces; σ, ε, ϑ are vectors gathering the independent (allowing for tensor symmetry) stresses, strains and imposed or inelastic strains (with the "engineering" definition of strains). The open domain where a problem is formulated is denoted by Ω; its boundary, assumed to be piecewise smooth, by Γ, the outward unit normal to Γ by n at x and by ν at ξ. The symbol G is used for Green's functions defined on the homogeneous unbounded space Ω∞. Other symbols are defined in the text where they are employed for the first time.

LINEAR ELASTIC PROBLEMS

Formulation of the boundary integral equations

All BEMs are deeply rooted in the mathematical theory of linear elasticity, which provided, since its golden age in the 19th century, central concepts such as effect superposition, influence functions and reciprocity relationships, and basic ingredients such as Kelvin's fundamental solution (1848), Somigliana's identity (1886) and Gebbia's kernels (1891). Therefore, quite naturally, we will refer below to linear elastic analysis in order to start discussing the SGBEM, formulating it and elucidating its distinction from, and links with, the traditional BEM. Naturally, numerous papers on SGBEM contain variants of its formulation and developments in elastostatics; a representative sample and source of information may be Refs. [START_REF] Bonnet | Boundary integral equations methods[END_REF][START_REF] Kane | Boundary element analysis in engineering continuum mechanics[END_REF][START_REF] Polizzotto | An energy approach to the boundary element method; Part I: Elastic solids[END_REF][START_REF] Sirtori | A Galerkin symmetric boundary element method in elasticity: formulation and implementation[END_REF] Consider a homogeneous elastic body which occupies the (closed) domain Ω = Ω ∪ Γ resulting from the union of the two-or three-dimensional open domain Ω and its boundary Γ. The boundary Γ is assumed to be smooth, i.e. endowed everywhere with a unique outward normal denoted by n or ν, depending on the symbol x or ξ, respectively, used for the point cartesian coordinates. Two elastic states, say the actual one (unstarred symbols) under given external actions and some fictitious one (starred symbols) can be involved in the following Betti's equation:

Γ (p T u ⋆ -u T p ⋆ ) dΓ ξ + Ω (b T u ⋆ -u T b ⋆ ) dΩ ξ = Ω (σ T ϑ ⋆ -ϑ T σ ⋆ ) dΩ ξ (1)
For brevity, no domain external action (b = 0, ϑ = 0) will be considered in what follows. If the fictitious state is identified as the Kelvin state associated to a point x ∈ Ω, i.e. if one sets ϑ ⋆ = 0 and models b ⋆ by a Dirac distribution ∆(ξx), then Eq. ( 1) becomes Somigliana's integral representation formula for displacements, which, 'taken to the boundary' (i.e. with x ∈ Γ) provides the usual starting point for the traditional BEM. Instead, let the (starred) fictitious state be identified as part of the response of the homogeneous elastic space Ω∞ embedding Ω, to distributions over Γ of static discontinuities f ⋆ (or single-layer sources) and of kinematic discontinuities d ⋆ (or double-layer sources). In order to set a convenient ground for the subsequent rigorous treatment of the singular integrals (Sec. 8), it would be suitable to conceive here the boundary Γ as the limit configuration (say Γ + ) of a moving surface Γ ⋆ outside Ω (cf. e.g. Ref. [START_REF] Sirtori | A Galerkin symmetric boundary element method in elasticity: formulation and implementation[END_REF]). Thus, having set b = 0 and ϑ = 0, we can write for any ξ ∈ Ω (not on Γ):

u ⋆ (ξ) = Γ Guu(ξ, x)f ⋆ (x) dΓx + Γ Gup(ξ, x)d ⋆ (x) dΓx (2) p ⋆ (ξ) = Γ Gpu(ξ, x)f ⋆ (x) dΓx + Γ Gpp(ξ, x)d ⋆ (x) dΓx (3)
where the source density functions are interpreted as traction and displacement jumps, respectively, across Γ:

f ⋆ (x) = p ⋆ (x -) -p ⋆ (x + ) (4)
d ⋆ (x) = u ⋆ (x + ) -u ⋆ (x -) (5)
The meaning of the new symbols are as follows: superscriptsand + denote points (ξ -, ξ + ; x -, x + ) or sets of points (Γ -, Γ + ) belonging to Ω (outward normal: n -= n) and to the exterior domain Ωe = Ω∞ -(Ω ∪ Γ) (outward normal: n + = -n), respectively; infinitely close to Γ; matrices G gather Green's functions for Ω∞, the former subscript specifying the effect (u for displacement, p for traction at ξ -with normal ν -), the latter the source that causes it, namely: u for static discontinuity concentrated as unit force, in the load point (here denoted by x), like in Kelvin fundamental solution; p for kinematic (displacement) discontinuity concentrated in x with normal n with respect to Γ and with unit integral over Γ, like in Gebbia solution, Refs. 29, 63, 87. The two-point influence functions contained in matrices G hk exhibit the following properties, denoting by r the distance defined as Euclidean norm (r = |x -ξ|):

(a) Singularities for r → 0; specifically, in three-dimensional problems:

Guu = O(1/r), Gup = O(1/r 2 ) Gpu = O(1/r 2 ), Gpp = O(1/r 3 ) (6)
and in two-dimensional problems:

Guu = O(ln r), Gup = O(1/r) Gpu = O(1/r), Gpp = O(1/r 2 ) (7)
(b) Reciprocity relationships flowing from Betti's theorem for r = 0:

G hk (ξ, x) = G T kh (x, ξ), (h, k = u, p) (8)
(c) Sign definiteness of quadratic forms expressing strain energies (see Sec. 7):

1 2 Γ Γ f T (x)Guu(x, ξ)f (ξ) dΓ ξ dΓx > 0, ∀f ≡ 0 (9) - 1 2 Γ Γ d T (x)Gpp(x, ξ)d(ξ) dΓ ξ dΓx > 0, ∀d ≡ 0 (10)
The latter inequality should be relaxed to ≥ (semidefiniteness) if d may represent a relative rigid-body motion between Ω and Ωe.

(d) Links are worth noticing between kernels concerning different sources: e.g. Gebbia's fundamental solution can be obtained from Kelvin's by taking first derivatives of it with respect to load-point coordinates ξ.

Basically, the elastic SGBIE formulation results from inserting u ⋆ (ξ), p ⋆ (ξ) (ξ ∈ Γ), expressed as integral representations, into Betti's equation [START_REF] Aimi | Numerical integration schemes for evaluation of (hyper)singular integrals in 2D BEM[END_REF]. However, the expressions (2-3) of u ⋆ , p ⋆ cannot be directly used to this purpose because of the kernel singularities [START_REF] Antes | The boundary integral approach to static and dynamic contact problems[END_REF][START_REF] Arnold | On the asymptotic convergence of collocation methods[END_REF]. Recalling classical results of elastic potential theory given

Γ Guu(x, ξ)p(ξ) -- Γ Gup(x, ξ)u(ξ) dΓ ξ dΓx + Γ d ⋆T (x) - 1 2 p(ξ) + Γ Gpu(x, ξ)p(ξ) -Gpp(x, ξ)u(ξ) dΓ ξ dΓx = 0 (13)
Let the given external actions p, ū act on two complementary disjoint parts of the boundary, say Γp and Γu respectively. Since equation (13) holds for any source distribution f ⋆ and d ⋆ , the two expressions in brackets must vanish separately, thus providing an overdetermined redundant set of BIEs. However, let us set f ⋆ = 0 on Γp and d ⋆ = 0 on Γu and denote by -ḡ u and -ḡ p the sums of the integrals containing the boundary data p and ū (and the additional four domain integrals in square brackets in Eq. ( 13) containing domain data b and θ, if these had not been ignored for brevity). Thus, partitioning Γ into Γp and Γu, the following two BIEs are generated from the Betti equation [START_REF] Banerjee | The boundary element method in engineering[END_REF] in view of the arbitrariness of f ⋆ (x + ) on Γu and d ⋆ (x + ) on Γp:

Γu Guu(x, ξ)p(ξ) dΓ ξ - Γp Gup(x, ξ)u(ξ) dΓ ξ = ḡu (x), (x ∈ Γu) (14)
and

Γu Gpu(x, ξ)p(ξ) dΓ ξ - Γp Gpp(x, ξ)u(ξ) dΓ ξ = ḡp (x), (x ∈ Γp) (15) with ḡu (x) = 1 2 ū(x) - Γp Guu(x, ξ)p(ξ) dΓ ξ + - Γu Gup(x, ξ)ū(ξ) dΓ ξ (x ∈ Γu) (16
)

ḡp (x) = 1 2 p(x) -- Γp Gpu(x, ξ)p(ξ) dΓ ξ + = Γu Gpp(x, ξ)ū(ξ) dΓ ξ (x ∈ Γp) (17)
As functions of x ∈ Ω∞, ḡu and ḡp are discontinuous across Γ (namely:

ḡu (x + ) -ḡu (x -) = ū(x), ḡp (x + ) -ḡp (x -) = p(x)).
The singularity issue can also be dealt with using a regularization approach (Sec. 8), which avoids Cauchy principal values or finite part integrals.

Let matrices Ψ collect shape functions defined over the whole boundary Γ, identically vanishing outside the relevant nodal 'support' and endowed with suitable continuity properties (dictated by integration requirements, see Sec. 8). The discretizations of both the actual boundary fields and the source distributions, respectively, materialize as follows:

p(ξ) = Ψp(ξ)P , u(ξ) = Ψu(ξ)U f ⋆ (x) = Ψ ⋆ p (x)F , d ⋆ (x) = Ψ ⋆ u (x)D (18)
Let the governing vectors in (14-15) be partitioned according to the subdivision of the boundary Γ into its free Γp and constrained Γu portions, i.e. (with h = u, p):

P T = [ P T , X T p ], U T = [X T u , Ū T ], Ψ h = [ u Ψ h , p Ψ h ] F T = [F T u , F T p ], D T = [D T u , D T p ], Ψ ⋆ h = [ u Ψ ⋆ h , p Ψ ⋆ h ] (19)
The matrices resulting from the double integrations implied by Eq. ( 13) when Eqs. ( 18) and ( 19) are substituted into it, can be represented in the form:

Ĝij hk = Γ i Γ j i ′ Ψ ⋆T h ′ (x)G T hk (x, ξ) j ′ Ψ k ′ (ξ) dΓ ξ dΓx ( 20
)
where subscripts h, k (= u, p) refer to the nature of the source and of the actual field; h ′ , k ′ (= p, u) to the variable which is workconjugate to h, k; i, j (= u, p) refer to the boundary portion Γu or Γp; i ′ , j ′ (= p, u) to its complement (Γp or Γu).

Through modelling, Eqs. [START_REF] Belytschko | A variationally coupled finite element-boundary-element method[END_REF], and partitioning, Eqs. [START_REF] Belytschko | Singular integration in variationally coupled FE-BE method[END_REF], and account taken of Eq. ( 20), Eq. ( 13) becomes:

F T p [ Ĝpu uu Xp - Ĝpp up Xu -B f p ] + F T u [ Ĝuu uu Xp - Ĝup up Xu -B f u ] + D T p [ Ĝpu pu Xp - Ĝpp pp Xu + B d p ] + D T u [ Ĝuu pu X p - Ĝup pp Xu + B d u ] = 0 (21)
where vectors B contain boundary data, Eqs. [START_REF] Becache | A variational boundary integral equation method for an elastodynamic antiplane crack[END_REF][START_REF] Becache | Regularization in 3D for anisotropic elastodynamic crack and obstacle problems[END_REF], and would contain also domain data b, θ if any.

In the discretized Betti equation ( 21) the vectors premultiplying the expressions in brackets are arbitrary and the diverse ways in which they are used to generate linear equations (as many equations as boundary unknowns contained in X T = [X T p , X T u ]) characterize various BEMs as briefly discussed in Sec. 2.2.

The SGBEM of concern herein arises when the following provisions are taken in Eq. ( 21): (i) as in the classical Galerkin's weighted residual approach, the source interpolations acting as weight functions are chosen equal to the interpolation functions adopted in order to model the actual fields within each pair of static (p, f ⋆ ) and of kinematic (u, d ⋆ ) variables, namely:

Ψ ⋆ i = Ψi, (i = p, u) (22)
(ii) setting F u = 0 and Dp = 0, the selected equations emerging from ( 21) are:

Ĝuu uu - Ĝup up - Ĝpu pu Ĝpp pp X p X u = B f u B d p (23)
or, more compactly, with self-evident meaning of the new symbols:

AX = B, with A T = A (24)
The symmetry of the coefficient matrix A, Eq. ( 23), can be proven by taking into account the Galerkin assumption Eq. ( 22) and the reciprocity properties (8) of the kernels in the double integration formulae, Eq. ( 20), which must be suitably interpreted and implemented as discussed in Sec. 8.

Miscellaneous issues and remarks

It is worth noticing that the discretized Betti equation ( 21) provides an unified basis from which, besides SGBEM, other BEMs proposed in the literature naturally descend with different choices of the source distributions and interpolation functions. Three of these alternative choices and consequent BE approaches are mentioned below.

(A) If corresponding shape and weight functions are equal according to Eq. ( 22) but kinematic discontinuities (double-layer) sources are renounced by setting d ⋆ = 0, then a nonsymmetric Galerkin BEM is formulated as proposed in Refs. 128, 129. The loss of symmetry, which occurs unless Γp = ∅, is partly compensated for by the avoidance of hypersingular kernels implied by the displacement jumps as sources. More recently, spectral discretizations in connection with this nonsymmetric GBEM have been proposed in Ref. 115.

(B) If the weight functions Ψ ⋆ are assumed as Dirac distributions associated to boundary points (so that the Galerkin approach, Eq. ( 22), is abandoned in favor of collocation at nodes) and only kinematic sources are employed (i.e. f ⋆ p = 0), then a nonsymmetric BEM emerges from Eq. ( 21), with attractive features in linear elastic fracture mechanics. This is basically the method developed by Crouch and Starfield, though mostly in the variant of an indirect formulation, Ref. [START_REF] Crouch | Boundary element methods in solid mechanics[END_REF] (a) The two BIEs in Eqs. [START_REF] Banerjee | Boundary element methods in engineering sciences[END_REF][START_REF] Banerjee | Advanced boundary element analysis of two-and three-dimensional problems of elastoplasticity[END_REF] can be rewritten in a compact (operatorial) fashion:

Ly = g (25)
where y stands for both the unknown fields on the boundary Γ = Γu ∪ Γp, g for the data functions on the r.h.s. of Eqs. (14-15) and L denotes the (linear, integral) operator which transforms the former into the latter. Making use of the reciprocity properties of the kernel G in (14-15), it can be easily proven, e.g. Refs. [START_REF] Bonnet | Boundary integral equations methods[END_REF][START_REF] Polizzotto | An energy approach to the boundary element method; Part I: Elastic solids[END_REF][START_REF] Sirtori | A Galerkin symmetric boundary element method in elasticity: formulation and implementation[END_REF], that for any pair of boundary fields y and z:

Γ z T Ly dΓ = Γ y T Lz dΓ ∀y, z (26)
Equation ( 26) expresses the symmetry (or self-adjointness) of the operator L with respect to the bilinear form defined there over Γ. An important consequence of this symmetry (attainable by a straightforward and traditional path of reasoning, e.g. Ref. [START_REF] Tonti | On the variational formulation for linear initial value-problems[END_REF]) is a variational characterization (i.e. stationarity as a sufficient and necessary condition) of the boundary solution, say ŷ, namely:

δ (1) Φ(y) = 0 (27) with Φ = 1 2 Γ y T Ly dΓ - Γ g T y dΓ
Because of the sign properties ( 9), [START_REF] Balakrishna | Efficient analytical integration of symmetric Galerkin boundary integrals over curved elements: thermal conduction formulation[END_REF], the functional Φ turns out to exhibit a saddle point at solution ŷ, with minimization and maximization with respect to p over Γu and u over Γp separately, Refs. 135, 154. The energy approach to, and interpretation of, the SGBEM in elasticity and elastoplasticity developed in Sec. 7 will cover the above variational properties as a special case and will put them in broader perspectives.

(b) So far, the provisions necessary to confer a meaning to the singular integrals have been tacitly assumed for granted and ignored. However, the presence of singularities stronger than those involved in the traditional BEM and in most of its variants represents a crucial feature of the SGBEM (and perhaps the very reason of its delayed developments). Therefore, Sec. 8 and parts of Secs. 9 and 12 of this review will primarily be devoted to various meaningful singularityrelated issues.

(c) In the consolidated literature on traditional BEMs, a distinction is often made between direct and indirect approaches. [START_REF] Babuska | Survey lectures on the mathematical foundations of the finite element method[END_REF][START_REF] Bamberger | Formulation variationnelle pour le calcul de la diffraction d'une onde acoustique par une surface rigide[END_REF] in view of the formal simplification entailed by the scalar nature of the unknown fields. In fact, the essential mathematical features (ellipticity in primis) are the same as in linear elastostatics and this makes it redundant to deal with such situations herein explicitly.

ELASTIC-PLASTIC QUASI-STATIC ANALYSIS

Integral equations and their space discretization

Clearly, the formulation of BIEs for elastoplastic analysis can be obtained as a generalized version of that for elastic analysis. However, instead of extending the path of reasoning based on Betti's theorem of Sec. 2, an alternative and simpler path will be followed below for the sake of approach diversity.

In the presence of material nonlinearities the SGBEM is believed by the writers to exhibit special advantages over traditional BEMs because a number of potentially useful results can be transferred from continua to discrete models. Therefore, in this Section, a rather detailed conspectus is presented of results which concern quasi-static elastoplastic analysis (see Refs. 42, 80, 81, 96, 103, 104, 106, 108-110, 112, 126, 133, 140, 142, 143), to be supplemented by those aspects of plasticity briefly dealt with in Secs. 4.3,5.1,6.1.

Let Ωp ⊂ Ω be the subdomain where plastic yielding is reasonably expected. Ignoring plasticity for a while, assume that the homogeneous elastic body considered be embedded in a suitable homogeneous solid (here, like in Sec. 2, the unbounded space Ω∞), of which Green's influence functions G are analytically known. Besides the static and kinematic sources on the boundary Γ (f ⋆ and d ⋆ ) employed in Sec. 2,imposed strains ϑ ⋆ are assumed in Ωp as proposed in Ref. [START_REF] Maier | A Galerkin approach to boundary element elastoplastic analysis[END_REF]. The following effects on Ω∞ due to all these sources are expressed by superposition: (i) displacements, to identify with data ū in the actual body, in points x of Ω neighbouring the constrained boundary Γu; (ii) tractions, to identify with data p in points x of Ω infinitely close to the free boundary Γp; (iii) stresses in Ωp.

The linearity of Ω∞ permits to superpose effects by means of influence functions. These functions, index h running over the above three kinds of effects (h = u, p, σ), can be identified with fundamental solutions or Green's functions of Ω∞, and precisely with: Kelvin's kernels G hu for static sources f ⋆ ; Gebbia's kernels G hp for displacement jump sources d ⋆ ; Bui's kernels G hσ for strain sources ϑ ⋆ . We name the Green's functions of the third kind after Bui in view of the substantial correction (addition of the 'convective term') provided in Ref. [START_REF] Bui | Some remarks about the formulation of threedimensional thermoelastoplastic problems by integral equations[END_REF] to the kernel Gσσ as it was widely used, until 1978, without that term in elastoplastic analysis by traditional BEM.

The actual elastic-plastic state in the solid considered is recovered by means of three provisions: (i) interpreting ϑ ⋆ as unknown plastic strains ε p ; (ii) making explicit the circumstance that the exterior domain Ω∞-Ω is undeformed and unstressed, i.e. setting u ⋆ (ξ + ) = 0 and p ⋆ (ξ + ) = 0 in the expressions of the boundary sources interpreted as discontinuities accross Γ ; (iii) entering the boundary data, i.e. u ⋆ (ξ -) = ū on Γu and p ⋆ (ξ -) = p on Γp and the boundary unknowns u ⋆ (ξ -) = u on Γp and p ⋆ (ξ -) = p on Γu.

Thus the integral expressions originally written as superpositions of the effects u and p on Γ -and σ in Ω, due to sources f ⋆ , d ⋆ on Γ and ϑ ⋆ = ε p in Ω, yield the three integral equations:

Γu Guu(x, ξ)p(ξ) dΓ ξ - Γp Gup(x, ξ)u(ξ) dΓ ξ + Ωp Guσ(x, ξ)ε p (ξ) dΩ ξ = ḡu (x), (x ∈ Γu) (28)
- Γu Gpu(x, ξ)p(ξ) dΓ ξ + Γp Gpp(x, ξ)u(ξ) dΓ ξ - Ωp Gpσ(x, ξ)ε p (ξ) dΩ ξ = ḡp (x), (x ∈ Γp) (29) Γu Gσu(x, ξ)p(ξ) dΓ ξ - Γp Gσp(x, ξ)u(ξ) dΓ ξ + Ωp Gσσ(x, ξ)ε p (ξ) dΩ ξ = ḡσ (x) + σ(x). (x ∈ Ω) (30)
Like in Eqs. [START_REF] Banerjee | Boundary element methods in engineering sciences[END_REF][START_REF] Banerjee | Advanced boundary element analysis of two-and three-dimensional problems of elastoplasticity[END_REF], ḡh (h = u, p, σ) denote the resulting terms gathering data on the boundary and also on the domain (body forces and thermal strains were ignored in Sec. 2 for brevity). Among the nine kernels that show up in Eqs. [START_REF] Bonnet | Regularized direct and indirect symmetric variational BIE formulations for three-dimensional elasticity[END_REF][START_REF] Bonnet | Boundary integral equations methods[END_REF][START_REF] Bonnet | Sensitivity approach for internal crack or void perturbation using boundary integral equation and adjoint approach[END_REF], four have been employed for elastic analysis and commented upon in Sec. 2.

As for the kinds of properties considered there, the new five kernels can be easily seen to exhibit the following features:

(α) Those concerning stresses as effects, i.e. G σh (h = u, p, σ) are singular like those concerning traction effects, i.e. G ph (h = u, p, σ); Bui kernels have the same singularities as Gebbia kernels, (β) The reciprocity relationships (8) hold for all kernels, i.e. for (h, k = u, p, σ), (γ) The quadratic form associated, for any domain, to Gσσ, is negative semi-definite, like that associated to Gpp over any closed surface, in view of the common meaning of the sign-inverse of the elastic strain energy due to (generally non-compatible) strains imposed in Ω∞ as external actions.

If the left-hand sides of Eqs. [START_REF] Bonnet | Regularized direct and indirect symmetric variational BIE formulations for three-dimensional elasticity[END_REF][START_REF] Bonnet | Boundary integral equations methods[END_REF][START_REF] Bonnet | Sensitivity approach for internal crack or void perturbation using boundary integral equation and adjoint approach[END_REF] together are interpreted as consisting of an integral operator L acting on fields {p T , u T , ε pT } T , then, using properties (β) and (γ), L can be shown to exhibit symmetry with respect to a bilinear form defined over Γ and Ωp. Clearly, remark (a) of Sec. 2.2 does not apply, as Eq. ( 30) contains the unknown field σ(x).

Symmetry is preserved when the discretization is performed in accordance with two provisions: (i) field modelling [START_REF] Belytschko | A variationally coupled finite element-boundary-element method[END_REF] extended to Ωp by ε p (ξ) = Ψ θ (ξ)Θ p , denoting by Θ p the vector of 'generalized' plastic strains; (ii) Galerkin weighted residual statement according to Eqs. [START_REF] Beskos | Boundary element methods in dynamic analysis, part. II (1986-1996[END_REF] with (h = u, p, σ), carried out for Eqs. [START_REF] Bonnet | Regularized direct and indirect symmetric variational BIE formulations for three-dimensional elasticity[END_REF][START_REF] Bonnet | Boundary integral equations methods[END_REF][START_REF] Bonnet | Sensitivity approach for internal crack or void perturbation using boundary integral equation and adjoint approach[END_REF]. In fact, the double integrations like in Eq. ( 20) with h, k and i, j running over the augmented index set (u, p, σ), if suitably executed in view of the integrand singularities (cfr. Sec. 8), generate symmetric algebraic linear equations. These can be written in the following compact form, which augments Eq. ( 24) to elastic-plastic analysis purposes:

AX + CΘ p = BΓ (31)
C T X + ĜσσΘ p = BΩ + Σ (32)
where:

A T = A, ĜT σσ = Ĝσσ, Σ = Ω Ψ T θ (ξ)σ(ξ) dΩ ξ (33)
Equation ( 31) condenses the Galerkin-discretized BIEs [START_REF] Bonnet | Regularized direct and indirect symmetric variational BIE formulations for three-dimensional elasticity[END_REF] and [START_REF] Bonnet | Boundary integral equations methods[END_REF].

Let its solution with respect to the vector X of boundary unknowns be substituted into Eq. [START_REF] Bui | An Integral Equation Method for Solving the Problem of a Plane Crack of Arbitrary Shape[END_REF], which is the Galerkin discretized version of integral equation [START_REF] Bonnet | Sensitivity approach for internal crack or void perturbation using boundary integral equation and adjoint approach[END_REF]. This move is implicitly or explicitly recurrent in the traditional elastoplastic analysis by BEM, e.g. Refs. 15, 31, 49, 158 (it is not so in computational plasticity by FEMs, where it would require the inversion of the whole elastic stiffness matrix). In the present SGBEM the equation resulting from these trivial manipulations,namely

ZΘ + Σ e = Σ (34)
turns out to exhibit the following special features:

Z T = Z ; - 1 2 Θ pT ZΘ p ≥ 0 ∀Θ p (35)
Clearly, the interpretation of Eq. ( 34) is equal to that of its counterpart from traditional BEMs, i.e.: actual stresses are the sum of the linear elastic stress response Σ e to external actions and selfequilibrated (in a suitable approximate way, Ref. 136) stresses due to plastic strains through the influence matrix Z. The properties (35) (i.e. symmetry of Z and the meaning of strain energy for the quadratic form associated to it) are in the SGBEM the same as those of kernel Gσσ for the continuum Ω∞ in local variables, whereas they are not valid in conventional BEMs. As a consequence of Eqs. [START_REF] Carini | A variational formulation of viscoelastic analysis by boundary integral equations[END_REF], the eigenvalues of Z are real, those that vanish correspond to stressless systems of plastic strains and their number can be regarded as a measure of the ability of the cell-discretization to accommodate mechanisms in plastic collapse and shakedown analysis of perfectly plastic solids and structures (see Ref. [START_REF] Maier | Symmetric Galerkin boundary element method in plasticity and gradientplasticity[END_REF] and Sec. 6.1).

Constitutive laws in local and generalized variables

The main peculiarities of SGBEM in inelastic analysis emerge even in the narrow context of associative, single-yield mode, perfect (nonhardening) plasticity, to which therefore the considerations in this Subsection will be restricted for brevity. Let the elastoplastic material model be described in the following classical fashion:

σ = E(ε -ε p ) ϕ = ∂ϕ ∂σ T (σ)σ -y ≤ 0, ; ε p = ∂ϕ ∂σ (σ) λ λ ≥ 0, ϕ λ = 0 (36)
Here, ϕ denotes the yield function, y a yield limit (meant to be a material constant), λ the plastic multiplier. As usual, the 'equivalent stress' y + ϕ is assumed to be positively homogeneous of order one in the stresses. The elasticity law (36a) postulates strain additivity; Eqs. (36b) express the yield criterion and the normality rule; Eqs. (36c) the so-called 'consistency' or 'loading-unloading' rule. The space discretization in the boundary and domain integral equations and their Galerkin approximate enforcement (Sec. 3.1) has entailed a link between the plastic strain model

ε p (x) = Ψ θ (x)Θ p (37)
and the stress-governing parameters Σ over Ωp, defined by Eq. (333).

It is highly desirable that vectors Θ and Σ gather 'generalized variables' in Prager sense (W. Prager and M. Save introduced this notion for structural plasticity in the late fifties), namely that:

Σ T Θ p = Ωp σ T (x)ε p (x) dΩ, ∀Σ, Θ p (38)
The identity [START_REF] Carini | Analytical integrations for two-dimensional elastic analysis by the symmetric Galerkin boundary element method[END_REF] entails the preservation of the scalar product (and its energy meaning) of conjugate variables occurring in pairs (kinematic and static) in passing from the local quantities to the global ones which govern the discrete model. This preservation is necessary (even if not sufficient) in order to endow the discrete model with the essential features of the original continuum solid or structure. In particular, Eq. ( 38), if combined with the symmetries of the boundary integral operator (Sec. 2.1), has far-reaching consequences in terms of attractive peculiarities of the SGBEM, as discussed in Sec. 7. Condition [START_REF] Carini | Analytical integrations for two-dimensional elastic analysis by the symmetric Galerkin boundary element method[END_REF] is complied with if the shape functions modelling stresses are suitably derived from those chosen for strains, i.e. if:

σ(x) = Ψσ(x)Σ, with Ψσ = Ψ θ Ωp Ψ T θ (x)Ψ θ (x) dΩ -1 (39
) Then Θ p in turn becomes weighted averages of ε p , shape functions Ψσ acting as weights, in full similarity to Eq. [START_REF] Bui | Some remarks about the formulation of threedimensional thermoelastoplastic problems by integral equations[END_REF]. Clearly, the same shape functions will be used for total as well as plastic strains, i.e. ε(x) = Ψ θ (x)E.

A fully analogous modelling procedure will be applied to the other pair of conjugate variable fields λ and ϕ by introducing interpolations Ψ λ and consequent, through a dependence like Eq. [START_REF] Cen | Bifurcations and instabilities in fracture of cohesive-softening structures: a boundary element analysis[END_REF], interpolations Ψϕ, with vectors Λ and Φ gathering Prager-generalized plastic multipliers and yield functions, respectively. The same procedure would hold for possible internal variables occurring in pairs in cases of hardening (instead of ideal) plasticity.

After the above preliminaries, straightforward manipulations lead from the material model in local variables to constitutive laws in weighted-average variables, generalized in Prager's sense, namely:

Σ = K(Θ -Θ p ), Φ = ∂Φ ∂Σ T (Σ)Σ -Y ≤ 0 Θp = ∂Φ T ∂Σ (Σ) Λ, Λ ≥ 0, Φ T Λ = 0 (40)
where:

K = Ωp Ψ T θ (x)EΨ θ (x) dΩ Y = Ωp Ψ T λ (x)y(x) dΩ
K and Y having the meanings of cell-wise stiffness matrix (symmetric, positive-definite, block-diagonal and intrinsic or 'natural', i.e. unaffected by rigid-body motions) and yield limit vector for all cells, respectively.

The relations [START_REF] Cen | On Boundary Element Inelastic Analysis in the Presence of the Softening[END_REF] can be interpreted as a description of the elastic-plastic global behaviour of the potentially yielding subdomain Ωp, cell-by-cell in a generally decoupled format. It is worth noting that the cell-wise constitution (40) reflects both the space modelling and the local material model; of the latter it possesses all essential features (perfect plasticity and Drucker's stability, i.e. normality and convexity) as it can be easily proven; clearly, this does not rule out local violations of the material model [START_REF] Carini | Variational formulation of the boundary element method in transient heat conduction[END_REF], (see Refs. [START_REF] Comi | Generalized variable finite element modelling and extremum theorems in stepwise holonomic elastoplasticity with internal variables[END_REF][START_REF] Maier | A Galerkin symmetric boundary-element method in plasticity: formulation and implementation[END_REF][START_REF] Maier | A Galerkin approach to boundary element elastoplastic analysis[END_REF].

Supplemented by displacement modelling (missing in quasistatic analysis by all BEMs), the above kind of relationships can be generated also in the FEM context, and in fact variants of them have been dealt with there under various labels such as 'consistent' FEs, multifield models (e.g. Ref. [START_REF] Comi | Generalized variable finite element modelling and extremum theorems in stepwise holonomic elastoplasticity with internal variables[END_REF], discretizations in 'natural' variables).

The space-discrete finite-step problem

The incremental elasto-plastic constitutive laws [START_REF] Cen | On Boundary Element Inelastic Analysis in the Presence of the Softening[END_REF] for all cells in Ωp have to be integrated in time t. Among various popular time-integration schemes, we choose here the implicit backwarddifference method, according to which the (nonholonomic, pathdependent) rate relations [START_REF] Cen | On Boundary Element Inelastic Analysis in the Presence of the Softening[END_REF] generate the following (stepwise holonomic) finite-step relations:

∆Σ = K(∆Θ -∆Θ p ), Φ = ∂Φ ∂Σ T (Σ)( Σ + ∆Σ) -Y ≤ 0 ∆Θ p = ∂Φ T ∂Σ (Σ)∆Λ, ∆Λ ≥ 0, Φ T ∆Λ = 0 (41)
where all variables (like Σ) at the starting instant t are known, the unknowns are the finite increments denoted by ∆ over the time interval ∆t, and the gradients of generalized yield functions Φ are meant at the final instant t + ∆t (i.e their argument is Σ = Σ + ∆Σ).

The elastoplastic cell constitution [START_REF] Cole | A numerical boundary integral equation method for elastodynamics I[END_REF] for Ωp must be associated to the linear Eq. ( 34) resulting from the elasticity integral equations, rewritten here for increments:

Z∆Θ p + ∆Σ e = ∆Σ (42)
Equations ( 41) and (42) together constitute the relation set that governs the SGBEM-discretized nonlinear boundary-value problem over the finite time step ∆t, the input increments of external actions being captured in the linear elastic stress response ∆Σ e . For the numerical solution to this nonlinear problem various procedures are available in the FE context of computational inelasticity and can be transferred to the BE context. The so-called 'modified' Newton-Raphson iterative method seems to be the most directly transferable and the most popular so far in the BE literature (see e.g. Refs. [START_REF] Banerjee | Advanced boundary element analysis of two-and three-dimensional problems of elastoplasticity[END_REF][START_REF] Brebbia | Boundary element techniques[END_REF][START_REF] Cruse | Elasto-plastic BIE analysis of cracked plates and related problems. Part. 1: formulation; Part2: numerical results[END_REF][START_REF] Telles | Implicit procedures for the solution of elastoplastic problems by the boundary element method[END_REF]. In fact, each iteration (say the r-th) basically consists of two phases: (i) 'prediction', namely a linear, global computation of ∆Θr with elastic stiffness for known ∆Θ p r-1 (provided by previous iteration or assumed at initialization); (ii) 'correction', namely a nonlinear local (cellwise) solution of Eqs. [START_REF] Cole | A numerical boundary integral equation method for elastodynamics I[END_REF] for given ∆Θr. The use of the original elastic stiffnesses in the predictive phase is more natural in BEMs, but by no means mandatory: classical Newton-Raphson procedures with consistent tangent matrices can be adapted at the price of some manipulations, both in traditional BEMs (e.g. Refs. [START_REF] Poon | Numerical implementation of an CTO-based implicit approach for the BEM solution of usual and sensitivity problems in elasto-plasticity[END_REF][START_REF] Telles | Implicit procedures for the solution of elastoplastic problems by the boundary element method[END_REF] and in SGBEM, Ref. [START_REF] Maier | Dynamic elastoplasticity analysis by Galerkin discretization of symmetric integral equations with elastostatic kernels[END_REF].

The finite-step problem of inelastic analysis represents one of the topics where the SGBEM exhibits clear advantages over traditional BEMs. In fact, using the symmetry of the matrix operators A and Z in the elasticity equations, combined with the domain fields modelling by Prager-generalized variables, permitted to achieve, among others, the following results, Ref. 42, on the finite-step problem [START_REF] Cen | On Boundary Element Inelastic Analysis in the Presence of the Softening[END_REF][START_REF] Cole | A numerical boundary integral equation method for elastodynamics I[END_REF] (actually in a more general context allowing for hardening and multi-yield modes): (a) the solution, if any, is characterized by an extremum property; (b) the modified Newton-Raphson procedure converges to the solution, if any, and the objective function contemplated in theorem (a) monotonically decreases along the iteration sequence; (c) algorithmic stability, i.e. the contractivity of a 'natural norm' (in Simo's sense, Ref. 152) along the time-step sequence is ensured unconditionally, i.e. for any time amplitude ∆t.

Like in classical plasticity of stable materials, only lack of hardening may jeopardize solution existence and uniqueness in the incremental boundary value problem. The above conclusions, and others of similar type expected in the same line of thought, allow an insight into, and the controllability of, inelastic analyses. They privilege the SGBEM since they are not available, neither are likely to be, in the traditional BEMs.

Material instabilities

In traditional BEMs only classical, Druckerian elastoplasticity appears to have been envisaged as constitutive model. Unstable material behaviours (in the sense of negative second-order work for some deformation disturbances) are technically important and form at present a subject of intensive research. Violations of Drucker's stability postulate such as softening (negative hardening) and lack of normality, may cause material instability which in turn has farreaching mathematical and computational consequences in elasticplastic analyses. In the BE area, loss of ellipticity and mesh dependence has been pointed out, Refs. 40, 103. Various provisions intended to restore objectivity and, in particular, to remove hypersensitivity to space discretizations, have been proposed in recent years, mostly with reference to FEMs: link through fracture energy between mesh length and constitutive softening; adaptive remeshing (Zienkiewicz et al.); fictitious time-dependence; polar Cosserat media; nonlocal continua; gradient plasticity.

Nowadays the last one of the above provisions appears to materialize a good compromise between regularization effectiveness and computational economy. It is also the only one investigated in the BE area so far, specifically in the context of SGBEM only, Ref. [START_REF] Maier | Symmetric Galerkin boundary element method in plasticity and gradientplasticity[END_REF]. A typical constitutive law of nonlocal, gradient plasticity is obtained if the yield criterion of Eq. (36b1) is generalized to:

ϕ = ∂ϕ ∂σ T (σ)σ -y -hλ + c∇ 2 λ ≤ 0 (43)
The two additional terms are: a linear softening term (with negative hardening constant: h < 0); a diffusive term containing the Laplacian ∇ 2 of the cumulative (time-integrated) plastic multiplier λ and a material constant c. The presence of a partial differential operator in the plastic constitutive laws requires suitable conditions to be assumed on the boundary ∂Ωp of the subdomain Ωp where plastic yielding cannot a priori be ruled out. The (nonlinear) boundary conditions, of dubious mechanical interpretation, proposed by Aifantis and Mühlhaus involves the space gradient ∇ and read:

(n T ∇λ)∆λ = 0 on ∂Ωp (44)
A variational saddle-point theorem, proposed for FEM in Ref. [START_REF] Comi | A generalized variable formulation for gradient dependent softening plasticity[END_REF], permits one to regard Eq. ( 44) as a natural boundary condition and to achieve a generalized, gradient-plasticity version of the cell constitutive laws without the explicit presence of that boundary condition, but with the addition in the yield criterion, Eq. (41a), of the hardening-softening terms (H + C)Λ where matrices H and C depend on the material constants h < 0 and c, respectively, and on the cell mesh as well (Ref. [START_REF] Maier | A Galerkin symmetric boundary-element method in plasticity: formulation and implementation[END_REF]. Thus the SGBEM of the preceding section was extended to softening associative plasticity with consequent regularization (in the sense of mesh independence also in the presence of softening), however at the price of a coupling among the cells in the correction phase of the solution procedure (coupling due to the C 0 continuity requirement generated on the plastic multiplier field by the Laplace operator).

FRACTURE MECHANICS

Linear elastic fracture analysis

Integral equation methods are frequently applied to the solution of fracture mechanics problems, either in a (semi)-analytical fashion, for simple geometrical configurations, or as a boundary element numerical solver. The main advantages provided by BEMs in this area are the much increased ease of geometrical modelling, especially in three-dimensional situations and for the simulation of crack propagation, and the high accuracy attainable in the evaluation of stress intensity factors. Moreover, the BIE/BEM approach is almost unavoidable for the numerical solution of dynamical problems like scattering of waves by cracks, fault modelling in seismology, and other related situations which frequently involve very large or unbounded domains.

Ref. [START_REF] Cruse | Boundary Element Analysis in Computational Fracture Mechanics[END_REF] provides an effective introduction to the traditional BEM in fracture mechanics. Due to the well-known degeneracy of the Somigliana identity when collocated at points of the crack faces, hypersingular traction BIE (HBIE) and their space discretization by collocation are often considered, and a substantial amount of research is undertaken in this area. Implementations of HBIE must be based on either regularization methods or direct evaluation of finite-part integrals. Either way, the 'density function' (e.g. the crack opening displacement (COD), i.e. the displacement jump accross the crack surface) must be C 1,αcontinuous at the collocation point. This requirement has severe consequences on the available choices of shape functions and collocation points, especially for 3D problems. The so-called 'dual' BEM for cracks combines displacement and traction collocation BIEs, and thus faces similar implementational difficulties.

On the other hand, it seems that only scattered efforts have been directed towards the application of SGBIE formulations to crack problems. Compared to HBIE formulations and collocation HBEM, the main advantage of the SGBEM is that it needs only C 0,αcontinuity for the kinematic unknown on the crack surface, thereby allowing the use of standard C 0 interpolations. Moreover, it will be pointed out in Sec. 10 that the symmetric character of Galerkin BIE formulations is valuable for sensitivity analysis; this fact leads to interesting techniques, briefly described in Sec. 10, for the evaluation of the energy release rate considered as the domain derivative of the potential energy at equilibrium.

Focussing now on a SGBIE formulation for linear elastic fracture mechanics, let Ω denote an elastic solid containing an internal crack idealized by an open surface Γ d , and denote by Γ the external boundary. In addition to the usual boundary data on the external portions

Gup(ξ, x)w ⋆ (x) dΓx (46) p ⋆ (ξ) = Γ [Gpu(ξ, x)f ⋆ (x) + Gpp(ξ, x)d ⋆ (x)] dΓx + Γ d Gpp(ξ, x)w ⋆ (x) dΓx (47)
The right-hand side in Betti theorem, Eq. ( 45), uses the fact (known from potential theory) that

u ⋆ (ξ + ) -u ⋆ (ξ -) = w ⋆ (ξ) (ξ ∈ S)
with u ⋆ defined by [START_REF] Costabel | Coupling of finite and boundary elements methods for elastoplastic interface problem[END_REF]. Then, the considerations developed in Sec. 2 can be followed again with obvious modifications caused by the presence of the new crack surface Γ d and relevant kinematic unknown w. Modelling the latter using shape functions Ψw(x) = Ψw(x) defined over Γ d , one eventually arrives at the symmetric linear system of equations:

    Ĝuu uu - Ĝup up Ĝuc up - Ĝpu pu Ĝpp pp - Ĝpc pp Ĝcu pu - Ĝcp pp Ĝcc pp          X p Xu Xw      =      B f u B d p B w c      (48)
using notation [START_REF] Belytschko | A variationally coupled FE-BE method for transient problem[END_REF] but with superscripts i, j = u, p, c referring to the surfaces Γu, Γp, Γ d , respectively. The new right-hand side B w c will be made explicit in Sec. 8.3. SGBIE formulations for the special case of a crack (or a set of cracks) isolated in an infinite elastic medium, easier to conceive (except for the regularization issue, which is left aside until Sec. 8), have been studied since about twenty years, see e.g. Ref. [START_REF] Bui | An Integral Equation Method for Solving the Problem of a Plane Crack of Arbitrary Shape[END_REF]. In fact, the system (48) becomes simply:

Ĝcc pp X w = - Γ d Ψ T w (x)p(x) dΓx (49)
Similar SGBIE formulations have been obtained for scattering of elastodynamic waves by isolated cracks, see e.g. Ref. 12. It is interesting to note that direct and indirect Galerkin BIE formulations coincide for problems involving cracks in infinite media, the primary unknown being in both cases the COD w; moreover, the strain energy W stored in the infinite body is given (see Eq. ( 80) of Sec. 7) by

W = - 1 2 X T w G cc pp Xw
The interpolations functions Ψw(x) for the COD, usually taken of the same type as the Ψu(x), need only to be continuous over Γ d . Near the crack front, one may use any of the special interpolation functions that have been developed in conventional BEM to represent the square-root behaviour of the COD, e.g. quarter-node elements. The SGBEM system of equations ( 48) can be used to perform usual linear fracture mechanics analyses. The implementation of standard post-treatments like J-integral computation or stress intensity factors evaluation using COD extrapolation is straightforward.

Numerical experiments in linear elastic fracture mechanics analysis were first achieved for cracks in infinite media, like in e.g. Refs. 72, 169. However, 2D numerical experiments for bounded isotropic domains with internal cracks are reported in Ref. 68, 171, based respectively upon direct and indirect Galerkin BIE formulations. An implementation for fracture analysis in plane orthotropic elasticity is reported in Ref. [START_REF] Gray | Symmetric Galerkin boundary integral fracture analysis for plane orthotropic elasticity[END_REF]. Finally, Ref. 95 presents a SGBEM implementation for three-dimensional linear fracture analysis for bounded elastic bodies, based on a formulation very close to Eq. ( 48).

Quasi-brittle fracture analysis

The computer simulation of fracture processes in concrete and concrete-like materials (often collectively referred to as 'quasibrittle' and including many geomaterials and several kinds of ceramics) is mostly based on the 'cohesive crack model' stemming from pioneering works of Barenblatt (1962) and Dugdale (1967): displacement discontinuity loci are allowed and endowed with an interface rigid-plastic softening law while linear elasticity is assumed everywhere else. This model rules out stress singularities at the tips and entails a 'process zone' behind the tip. Such now popular idealization naturally suggests to conceive the analysis as a coalescence of two substantially diverse ingredients (a) an experimentally corroborated constitutive nonlinear dependance of the tractions p from the displacement jump w (and sometimes also of its previous time history) in all points of the locus, say Γ d , where nonzero w are reasonably expected; (b) a linear relationship between the same variables through the body or structure deprived of Γ d (i.e. of any kinematic discontinuity) and subjected to the given external actions in a linear elastic regime. In compact symbols, the above two ingredients at time instant t read, respectively:

p(x, t) = f (w(x, τ ), 0 ≤ τ ≤ t), x ∈ Γ d (50) p(x, t) = p E (x, t) + Γ d Z d (x, ξ).w(ξ) dΓ ξ (51)
The latter equation exhibits the following noteworthy features: (i) it reflects the overall geometry of the domain and the elasticity of the material in Ω; (ii) all loads are allowed for through the (linear elastic) traction p E response to them across Γ d only; (iii) the influence function matrix Z, or Green's kernel of the considered body, is symmetric due to the Betti's reciprocity identity, and negative semidefinite due to its energy meaning (definite if Ω is simply connected with respect to Γ d ):

Z d (x, ξ) = Z T d (ξ, x) - 1 2 Γ d Γ d w T (x)Z d (x, ξ).w(ξ) dΓx dΓ ξ ≥ 0 (52)
The confinement of all nonlinearity to the locus Γ d with lesser dimensionality with respect to Ω and the consequent formulation (50) of the analysis problem, clearly advocate a BIE-BE approach for its numerical solution. The above listed features, especially those in (iii), provide a strong motivation for solving and analyzing quasi-brittle fracture problems by the present SGBEM, because this method leads to approximate discrete formulations which preserve all those properties, while the traditional BEM does not. Such privilege of the SGBEM rests not only on expected computational gains (e.g. in terms of accuracy at equal number of unknowns), but especially on the achieved possibility of physical insight and evaluations. The importance of an enrichment of the theory underlying the numerical computations becomes quite apparent in quasi-brittle fracture mechanics if the constitutive law of Γ d is focused. In general, the relationship symbolically represented by Eq. ( 50) is nonlinear, nonholonomic, softening, sometimes nonassociative and multidissipative (i.e. with a multiplicity of yield modes) and can often be described by means of internal variable rigid-plastic models.

The following kinds of analysis problems can be singled out, depending on the engineering situation and on the consequent specialization of the relationship [START_REF] Davì | A hybrid displacement variational formulation of BEM for elastostatics[END_REF].

(A) Problems in rates (or in infinitesimal increments) allowing for the irreversible, nonholonomic nature of the quasi-brittle fracture process. Then the locus Γ d reduces to the current process zone (where currently w = 0 and p = 0). An important purpose of rate formulations is to check overall stability (characterized by the positiveness, in any possible virtual kinematic disturbance, of the overall secondorder work done by an external agency which preserves equilibrium while promoting the kinematic disturbance), see Ref. [START_REF] Maier | Symmetric Galerkin boundary element method for quasi-brittle fracture and frictional contact problems[END_REF]. Another purpose is to capture the onset of possible bifurcation (in the sense of path-equilibrium branching) and to compute the whole multiplicity of rate solutions as beginning of other alternative fracturing processses, see e.g. Ref 39. Finally, in the very special but frequent case of mode I only (i.e. because of symmetry) along an a priori known path with piecewise-linear decay of traction p for increasing opening displacement w, the rate solution can be amplified (by solving a trivial linear programming problem) up to the activation of a new yielding mode, thus reducing the time-stepping 'exact' analysis to a sequence of steps, each one of which simply consisting of a nonlinear rate solution and its linear amplification.

(B) Problems in total variables, i.e. based on a holonomic formulation. For instance, in simulating three-point bending tests (e.g. to parameter identification purposes), the crack propagation path along the symmetry axis and a monotonic increase of the displacement jump w ('regularly progressive yielding' in the plasticity jargon) can be reasonably conjectured. In cases like this, where irreversibility manifestations can a priori be ruled out and the discontinuity locus Γ d assumed, the analysis can be formulated as a non-linear elastic problem (in the spirit of the 'deformation theory' of plasticity), cf. Ref 24.

(C) Formulations of the cohesive model which are in a sense intermediate between (A) and (B) naturally occur when the time evolution of the fracture process has to be simulated without prior knowledge and with full nonlinearity of the interface law. Then, like in computational plasticity (cf. Sec. 3.3), a stepwise-holonomic time-integration scheme has to be derived, including in each finite step the updating of internal variables to allow for irreversibility and a direction search algorithm for the (process zone) tip advancement.

Focussing now on the linear-elastic background described by the integral equation [START_REF] Davì | A symmetric and positive definite variational BEM for 2-D free vibration analysis[END_REF], its discretization according to the SGBEM can be carried out by the same path of reasoning followed in Sec. 3 for inelastic analysis, basically with the only difference that all constitutive nonlinearities envisaged are confined to the kinematic discontinuity locus Γ d with the same dimensionality as the boundary Γ (instead of a subdomain Ωp). This difference implies that double-layer (displacement jump) sources, to be identified with the unknowns w, are adopted within the domain Ω on Γ d (instead of concentrated strain sources θ to be identified with plastic strains in Ωp). Otherwise the conceptual and operative sequence is the same as in Sec 3, namely: embedding in Ω∞; superposition of effects (p and u on Γ -and p on Γ d ) in Ω∞ through its Green functions Gij with i, j = u, p (no longer also σ); introduction of boundary data and unknowns on Γ -and, as for Γ + , attribution of a stressless undeformed state to Ω∞ -Ω, thus leading to BIEs analogous to Eqs. ( 28)-( 30), formally with subscript p instead of σ, Γ d instead of Ωp and w instead of ε; modelling the unknown fields and enforcement of the integral equation in the Galerkin weighted-residual fashion; double integrations on elements along Γ and Γ d to generate matrices of coefficients and vectors of data, like in Eqs. ( 31)- [START_REF] Bui | Some remarks about the formulation of threedimensional thermoelastoplastic problems by integral equations[END_REF]; condensation of the boundary unknowns, cf. Eqs. [START_REF] Bui | On the variational boundary integral equations in elastodynamics with the use of conjugate functions[END_REF]. This sequence leads to a discrete counterpart to Eq. ( 51) in the form

P = P E + Z d W (53)
The noteworthy circumstances in this SGBE formulation is that matrix Z d is endowed with the same essential features of its continuum counterpart, Eqs. ( 52), and that P and W governing the unknowns on Γ d are vectors of generalized variables (tractions and displacement jumps, respectively), i.e. they are conjugate in a virtual work sense. Once again, it is worth stressing that the above features, on one hand, would be missing in a traditional BEM and, on the other hand, are especially beneficial in quasi-brittle fracture, where material instability (softening) and its structural consequences are the peculiar ingredient and the main origin of computational difficulties.

In fact, let us consider as an illustration rate problems (A), and only them for brevity. Borrowing again the formalism of incremental 'nonstandard generalized' plasticity, a broad class of interface relationships for the discontinuity locus Γ d (now reduced to the instantaneous process zone) can be represented in terms of rates as follows, after the transition from local to generalized variables consistently with the Galerkin discretization which led to the elasticity Eq. ( 53):

Ẇ = ∂ ΦT ∂P Λ, Ṡ = - ∂ ΦT ∂Q Λ, Q = ∂ 2 Π ∂S∂S T Ṡ (54) Φ = ∂Φ ∂P T Ṗ + ∂Φ ∂Q T Q ≤ 0, Λ ≥ 0, ΦT Λ = 0 (55)
Here Φ and Φ are vectors of yield functions and plastic potentials, respectively; S and Q denote kinematic and static internal variables, respectively, and Π their potential. All derivatives are evaluated in the current situation, at which the incremental process starts. The combination of Eqs. [START_REF] Diligenti | Finite part integrals: their occurrence and computation[END_REF][START_REF] Dostal | Symmetric FE-BE coupling with iteration on the interface[END_REF] with Eq. ( 53) in rates leads to a linear complementarity problem (LCP) of the type (Ref. 102):

-Φ = Ḃ + M Λ ≥ 0, Λ ≥ 0, ΦT Λ = 0 (56)
having set

M = ∂Φ ∂Q T ∂ 2 Π ∂S∂S T ∂ ΦT ∂Q - ∂Φ ∂P T Z d ∂ ΦT ∂P , Ḃ = - ∂Φ ∂P T Ṗ E
Among several remarks and developments on the LCP (56), those mentioned below are especially pertinent to SGBEM. (a) Matrix M is symmetric if Φ = Φ, i.e. with associative cohesive crack models; it is not definite nor semidefinite in general, because the former addend on the right-hand side is not so for softening behaviour (though the latter addend is positive definite or at least semidefinite).

(b) In view of (a) the solution set can be empty or contain a finite number of elements, the latter event having the mechanical meaning of a branching of the system evolution (see e.g. Ref. 39). To compute the whole multiplicity of solutions of a LCP like ( 56) is here an important task but also still a challenge to applied mathematics. In fact, enumerative methods, which guarantee to provide all solution (with finite termination), entail a computing time exponentially growing with the problem size; on the other hand, efficient, asymptotic iterative Newton-type methods lead to a solution depending on the initialization, but hardly lead to all of them (see Refs. [START_REF] Bolzon | Some aspects of quasibrittle fracture analysis as a linear complementarity problem[END_REF][START_REF] Bolzon | On multiplicity of solutions in quasi-brittle fracture computations[END_REF].

(c) In various ways the LCP ( 56) is amenable to generally nonconvex quadratic programming, hence to optimization procedures alternative to those mentioned in (b).

(d) The quadratic form associated to matrix M , if Φ = Φ, can be shown to provide with its copositiveness a necessary and sufficient condition for overall stability, with its positive semidefiniteness a sufficient one for it (Ref. [START_REF] Bolzon | Some aspects of quasibrittle fracture analysis as a linear complementarity problem[END_REF].

(e) The interface law drastically simplifies for fracturing processes in the case of cohesive crack model for mode I with linear softening. In the holonomic (B) and stepwise-holonomic (C) formulations outlined earlier, the LCP structure may be preserved in the case of piecewise linear models or may give way to more general nonlinear complementarity problems (NLCP), remarks (a) to (c) still basically holding with suitable adjustments.

As a conclusive remark, it can be stressed that the SGBEM is ideally suited to quasi-brittle fracture analysis. In fact, SGBEM not only leads to formulations which economically exploit the linear background of the problems and the confinement of nonlinearity to the discontinuity locus Γ d , but also faithfully reflects, through the properties of Z d , the underlying mechanics of the problem and, hence, gives rise to beneficial consequences, as those listed (a) to (e) above, which can hardly be proved in the context of traditional BEMs.

Further remarks on fracture mechanics

A variety of issues in mechanics of materials and structures are closely related to the developments outlined in what preceeds and equally susceptible to advantageous applications of the SGBEM. They will be briefly mentioned below.

Delamination phenomena in inhomogeneous structural components can be regarded as a special decohesion process with a priori known locus Γ d of potential displacement discontinuities and, hence, is implicitly covered by the remarks of Sec. 4.2 on quasi-brittle fracture, provided geometric effects are negligible (e.g. simulations of local buckling in laminates would require extensions of SGBEM not available so far).

Unilateral contact, without or with friction or with more sophisticated interface models, represent a subject per se, extensively dealt with in the BE literature, see e.g. Antes and Panagiotopoulos, Ref. [START_REF] Antes | The boundary integral approach to static and dynamic contact problems[END_REF]. In some important engineering situations, unilateral contact is related to fracture: e.g. it must be allowed for in integrity assessments of concrete dams in the presence of cracks (and possibly of joints) under variable loads (e.g. earthquakes, seasonal thermal cycles, reservoir level fluctuations). SGBEM in contact was investigated specifically in Ref. [START_REF] Polizzotto | Variational boundary-integral-equation approach to unilateral contact problems in elesticity[END_REF] and in association to quasi-brittle fracture in Ref. [START_REF] Maier | Symmetric Galerkin boundary element method for quasi-brittle fracture and frictional contact problems[END_REF]. Methodologically, from the present standpoint unilateral contact may be considered as basically covered by the discussion of Sec. 4.2.

Elastic-plastic fracture mechanics encompasses at present several unsettled questions in physical interpretations, mathematical models and computing procedures. Traditional BEMs and their recent developments applied to ductile (large yield zone) fracture, have been extensively discussed by Leitão, Ref. [START_REF] Leitao | Boundary elements in nonlinear fracture mechanics[END_REF]. As for the SGBEM, the following remarks may supplement this monograph.

(a) Elastic-plastic fracture processes in metal structures are fre-quently simulated in industrial environments by geometric approaches (tip advancement criteria based on crack-tip opening angle (CTOA) at the onset of and during crack propagation, respectively) or energy approaches (based on J-integral or other path-independent integrals). Since these approaches still imply discrete modelling of fractures in association with classical constitutive elastic-plastic laws, stable in the sense of Drucker, the comparisons between SGBEM and traditional BEMs pointed out in Secs 3.1 to 3.3 (plasticity) and 4.1, 4.2 (brittle and quasi-brittle fracture) can be transferred to ductile fracture analysis by the aforementioned approaches.

(b) The present trend in elastic-plastic fracture mechanics is towards the use of micromechanically corroborated material models (such as Gurson model with porosity internal variables), capable to describe in a diffused fashion the damage process up to material separation. These models exhibit softening and, hence, require provisions apt to remove subjectivities like mesh dependence. A particular one of such remedies, i.e. gradient plasticity, was discussed with reference to SGBEM in Sec. 3.4. It is worth noting that the cohesive models for discrete crack simulations (Sec. 4.2) do not exhibit the above subjectivity, since dissipative damage process are 'squeezed' into the discontinuity locus Γ d .

(c) As emphasized e.g. in Ref. [START_REF] Leitao | Boundary elements in nonlinear fracture mechanics[END_REF], in elastic-plastic fracture by traditional BEMs, recourse is made to the so-called 'dual' formulation (in the sense that displacements are Somigliana-represented on one side of the potential crack and tractions on the other, in order to avoid both ill-conditioning and zonewise multidomain formulations). It is worth noting that in SGBEM the issue does not arise despite single-domain formulation, as could be observed in Secs. 4.1-4.2.

(d) Finally, ductile fracture processes usually imply yielding zones large with respect to the crack extension and local geometric features (e.g. thickness), but small with respect to a typical length of the overall geometry. This is the case of a fracture event in a pressure vessel or a pipeline, the economical simulation of which may materialise into three-dimensional BE modelling around the propagating cracks and shell FE modelling elsewhere in the structure. Thus, BE-FE coupling becomes a particularly favorable option, for which SGBEM has a special appeal, as shown in Sec. 11 of this survey.

TIME-DEPENDENT PROBLEMS

Linear elastic dynamics

SGBEM formulations of linear elastodynamic analysis problems for three-dimensional solids were formulated in Ref. [START_REF] Polizzotto | A variational approach to boundary element methods[END_REF]. The method therein employed is similar to that of elastostatics, but the Wheeler-Sternberg formulas, Refs. 41, 58, were used instead of Somigliana's. With a notation similar to that previously employed in elastostatics, the Wheeler-Sternberg formulas can be written:

u p σ = Γu G hu ⋆ p dΓ - Γp G hp ⋆ u dΓ -ḡh , (h = u, p, σ) (57)
where the asterisks denote convolution integrals, i.e., for instance,

G hu ⋆ p = t 0 G hu (x, ξ, t -τ )p(ξ, τ ) dτ. ( 58
)
The two-point time-dependent matrix-valued functions G hk = G hk (x, ξ, tτ ), (h, k = u, p, σ), collect Stokes fundamental solutions. Their explicit expressions are presented in the Appendix of Ref. 100. Interpreted as influence or Green's functions, the kernels G hk give in x ∈ Ω∞ at time t the effect specified by the first subscript of G hk (i.e. displacements for h = u, tractions on the surface element of normal n(x) for h = p and stresses for h = σ) and caused by a unit discontinuity source applied at ξ = x in Ω∞ at time τ < t, this source being specified by the second subscript of G hk through a conjugation rule (i.e. unit impulse of force for k = u, relative displacement layer on the surface element with normal n(ξ) for k = p and distortion for k = σ). The symmetry properties, Eq. ( 8), still hold: G T hk (x, ξ, t) = G kh (ξ, x, t) ∀x, ξ = x, ∀t ≥ 0 and for all h, k = u, p, σ. It can be assumed that τ ≥ 0 and, hence, G hk (x, ξ, 0) = Ġhk (x, ξ, 0) = 0 ∀x, ξ = x. The terms -ḡ h (h = u, p, σ) account for the given load history and initial conditions, i.e.:

-ḡ h = Γp G hu ⋆ p dΓ - Γu G hp ⋆ ū dΓ + Ω G hu ⋆ b dΩ + Ω G hσ ⋆ θ dΩ + Ω ρ Ġhp u0 dΩ + Ω ρG hu u0 dΩ (h = u, p, σ) (59)
where ρ is the mass density and the overbars denote given functions of space coordinates and time.

Equations ( 57) provide the displacements u(x, t), tractions p(x, t) = σ(x, t) • n(x) and the stresses σ(x, t) at point x ∈ Ω∞ as superpositions of the responses ḡh to loads and initial conditions (u0, u0), as well as to the unknown force layer p on Γu and relative displacement layer -u on Γp. By their very nature, they satisfy the motion equations in Ω but not the boundary conditions. The latter conditions, enforced in a symmetric Galerkin mode, read:

Γu Ψ T u (u -ū) dΓ = 0, Γp Ψ T p (p -p) dΓ = 0 ( 60
)
where Ψu and Ψp are shape function matrices modeling in space the boundary unknowns u and p, namely:

u(x, t) = Ψu(x)Xu(t) on Γu p(x, t) = Ψp(x)X p(t) on Γp (61)
Substituting from (57) for u and p into Eq. ( 60) and using Eq. ( 61), one obtains:

Auu ⋆ X p -Aup ⋆ Xu = Bu(t) Apu ⋆ Xp -App ⋆ Xu = Bp(t) (62)
having set:

A hk (t) = Γ h Γ k Ψ T h (x)G hk (x, ξ, t)Ψ k (ξ) dΓ dΓ = A T kh (t), (h, k = u, p) (63)
B h (t) = Γ h Ψ T h (x) ḡh (x, t) dΓ, (h = u, p) (64)
where:

ḡu = ū + ḡu on Γu, ḡp = p + ḡp on Γp (65)
Equations ( 62) are the space-discretized equations of motion and constitute a symmetric system of Volterra integral equations of the first kind in time, with unknowns Xp(t) and -Xu(t), these unknowns being the time histories of the nodal values of the boundary tractions Xp(t) and of the sign-reversed displacements -Xu(t).

Their integration can be achieved using a time-stepping procedure, see e.g. Ref. [START_REF] Karabalis | A simplified 3-D time-domain BEM for dynamic soil-structure interaction problems[END_REF].

An alternative is as follows. Let Eq. ( 62) be rewritten in the compact form A ⋆ X = B(t), ∀t ≥ 0 (66) where

A = Auu Aup Apu App X = Xp -Xu B = Bu Bp
and let the time-primitive Ã(t) of A(t) be introduced, i.e.

A(t) = Ȧ(t), Ã(0) = 0.
Note that à can be derived from the primitive fundamental solutions Ghk as A is derived from the G hk = Ġhk ( Ghk = 0 for t = 0), Ref. 141. Thus, through integration by parts, Eq. ( 66) can be reformulated as

à ⋆ Ẋ = B(t) -Ã(t)X(0). (67)
Denoting by t0 = 0, t1, t2, . . . a subdivision of the time axis into equal intervals of length ∆t, using the notation X (n) = X(tn) and writing Eq. ( 67) at tn with initial conditions at tn-1, one arrives at the equation

tn t n-1 Ã(tn -τ ) Ẋ(τ ) dτ = ∆B (n) -Ã(∆t)X (n-1) (68)
where

∆B (n) = B (n) -B (n-1)
. Modelling X as e.g. piecewiselinear in time, i.e. Ẋ (τ ) = Ẋ(n) = ∆X (n) /∆t in the n-th interval and noting that tnτ = ∆t -(τtn-1) = ∆tτ ′ , one obtains from Eq. ( 68):

J ∆X (n) = ∆B (n) -Ã(∆t)X (n-1) (69)
where the matrix J is symmetric and step-independent, i.e.

J = 1 ∆t ∆t 0 Ã(∆t -τ ′ ) dτ ′ .
The principal submatrices of J , i.e. J uu and J pp, can be shown to be, respectively, positive and negative definite. X (0) can be computed taking into account the initial conditions (u0, u0). Applying Eq. ( 69) sequentially for n = 1, 2, . . ., the vectors X (1) , X (2) , . . . collecting the station values of the boundary unknowns that solve the motion equations ( 67) can be computed. The numerical implementation of the above integration procedure is still lacking. In Ref. 141, the Galerkin boundary conditions [START_REF] Frangi | Symmetric BE method in twodimensional elasticity: evaluation of double integrals for curved elements[END_REF] were derived making use of a variational procedure based on the Hu-Washizu principle for dynamics (cf. Sec. 7). The same equations ( 60) can be obtained by an alternative variational approach to the dynamic problem proposed in Refs. [START_REF] Maier | A variational formulation of the boundary element method in elastodynamics[END_REF] and 100 (see Sec. 7).

Elastic-plastic dynamics

The SGBEM for elastodynamics outlined in Sec. 5.1 was extended to elastic-plastic dynamics (Refs. 125 and 127) through a procedure quite similar to that employed in statics. The key concept consists of replacing the imposed strains θ of Eq. ( 59) by the sum θ + ε p , where ε p (x, t) denotes the unknown plastic strains, but adding the new integral with ε p as a third integral before the data vectors ḡh in Eq. [START_REF] Dumont | The hybrid boundary element method[END_REF]. With this obvious changes, Eq. ( 57) reads:

   u p σ    = Γu G hu ⋆ p dΓ - Γp G hp ⋆ u dΓ + Ω G hσ ⋆ ε p dΩ -ḡh , (h = u, p, σ) (70)
where the data vectors ḡh are still given by Eqs. [START_REF] Frangi | A Galerkin symmetric and direct BIE method for Kirchhoff elastic plates: formulation and implementation[END_REF]. These equations describe the elastic response of Ω∞ to the given loads and initial conditions, as well as to the plastic strains ε p (x, t) treated as imposed strains. The Galerkin boundary conditions (60) still hold. A domain discretization is required for the field ε p and can be achieved in exactly the same way as in statics (see Sec. 3), that is, subdividing the domain in appropriate cell elements and writing

ε p (x, t) = Ψ ϑ (x)Θ p (t) in Ω
where Ψ ϑ is the relevant shape function matrix and Θ p (t) collects node values of plastic strains. Thus, following the same path as in elastodynamics, one arrives at an equation similar to [START_REF] Giroire | A new system of boundary integral equations for plates with free edges[END_REF], i.e.

A ⋆ X + C ⋆ Θ p = B(t), ∀t ≥ 0 (71)
where:

C = Cu Cp , C h (t) = Γ h Ω Ψ T h G hσ Ψ ϑ dΩ dΓ (h = u, p). (72)
The motion equation ( 71) must be supplemented by the equation

Σ = C T ⋆ X + Ĝσσ ⋆ Θ p -BΩ, ∀t ≥ 0 ( 73
)
where, similarly to the quantities defined in Secs. 2 and 3 but now with the dependence on time:

Σ(t) = Ω Ψ T θ σ dΩ, BΩ(t) = Ω Ψ T θ ḡ dΩ, (74) Ĝσσ
(t) = Ω Ω Ψ T θ GσσΨ θ dΩdΩ. (75)
The generalized stresses Σ are conjugate (or dual) to Θp in the sense of Prager that the product Σ T Θp equals the total continuum plastic dissipation, see Sec. 3.2. Equation ( 73) is obtained by substituting into Eq. ( 74)1 the expression (57) of σ. Equations ( 71) and ( 73) are Volterra equations of the first kind which must be solved, through a typical step-by-step procedure of numerical plasticity, with the aid of a set of additional relations which, like in quasi-static plasticity (cf. Eqs. 40, Sec. 3), account for the plastic material model suitably approximated for the cell aggregate, see Ref. [START_REF] Panzeca | Dynamic plasticity analysis by boundary-interior element[END_REF].

The SGBEM in the above formulation for elastic-plastic dynamics was not numerically tested so far, to the writers' knowledge. Probably the main obstacle to overcome in computer implementations and practical applications rests in the generation, through Eqs. ( 63), ( 72) and (74)2, of matrices A (which intervenes also in linear elastic dynamics, Sec. 5.1), C, Ĝσσ in Eqs. ( 71) and [START_REF] Guiggiani | A general algorithm for multidimensional Cauchy principal value integrals in the boundary element method[END_REF]. In fact, as a remarkable and computationally crucial difference from the quasi-static time-independent context, the entries of these matrices are functions of time, not numerical coefficients. Such difficulty can be attenuated, as proposed in Refs. 99 and 100, if matrices of numbers are generated by performing, once and for all in a coordinated fashion, the four integrations required (two in space and two in time over the step or, in elasticity, over the time interval of interest).

As an alternative to the preceding formulations, a SGBEM for dynamic inelastic analysis has recently been formulated and implemented (Ref. [START_REF] Maier | Dynamic elastoplasticity analysis by Galerkin discretization of symmetric integral equations with elastostatic kernels[END_REF]) with recourse to static fundamental solutions. This approach avoids the aforementioned difficulty. In fact, its starting point merely consists in adding to each l.h. side of Eqs. ( 28)-( 30) a domain integral containing inertial forces -ρü(ξ). To compensate for the presence of these unknown body forces in Ω, displacements u(x) in Ω have to be represented as effects as well, giving rise to a further integral equation. The unpleasant consequence is that the domain integrals with unknowns to model are not only the usual ones with plastic strains but also these reflecting the inertia effects. The price seems to be worth paying since, after space-modelling, the construction (by double integrations) of the coefficient matrices is carried out like in the quasi-static regime, Sec. 3.1. Subsequently, any integration scheme of inelastic dynamics may be employed to formulate the finite-step problem basically like in Secs. 3.2 and 3.3. In the presence of constitutive instabilities, Sec. 3.4, the stabilising effects of inertia should be properly allowed for and exploited like in FEM.

Numerical tests and comparisons with solutions generated by a commercial FE code (ABAQUS) and by traditional BEM from the literature, partly presented in Ref. [START_REF] Maier | Dynamic elastoplasticity analysis by Galerkin discretization of symmetric integral equations with elastostatic kernels[END_REF], turn out to be quite encouraging.

The use of static kernels in dynamic problems has already been developed in traditional BEM, sometimes with the provision of transforming the additional inertial domain integral into boundary integrals by Nardini

Other time-dependent problems

Like in other branches of computational mechanics, the developments in the SGBEM of concern herein tend to primarily focus on elasticity and plasticity, partly because of the centrality of these areas in structural mechanics and engineering, where the motivations of those developments and the background of most researchers are located. Outside the core areas dealt with in what precedes, some attention has been paid so far to the SGBIE-BE approaches in the following categories of linear evolutive analysis problems in chronological time (in contrast to event-ordering time variable of quasi-static inelastic analysis). (a) Transient heat conduction, implicitly covering through analogy other diffusion phenomena described by strictly identical mathematical models, such as Darcy filtration in porous media.

(b) Linear viscoelasticity based on the classical constitutive law expressed by a convolutive integral understood in the sense of Stieltjes in order to cover input variables (stress or strain tensor) discontinuous in time.

(c) Poroelasticity in Biot sense of an elastic solid skeleton the porosity of which is fully saturated by a liquid. By analogy, other twophase coupled problems can be regarded as implicitly covered, such as linear thermoelasticity.

(d) Linear acoustics in media with uniform (and in some cases piecewise uniform) constitutive properties. The mathematical model is usually the scalar or vector wave equation, and other physical kinds of wave propagation, e.g. electromagnetic waves, are covered. Most applications of computational acoustics deal with unbounded media.

Traditional BEM have been extensively studied in the last few decades for the numerical solution of the above initial-boundary value problems arising in a number of diverse technologies. Surveys of the abundant relevant literature can be found, e.g., in the books Refs. 13, 29, 31. As for the present SGBEM the few results available so far can briefly be outlined as follows.

The studies so far available on SGBIE-BEM for problem categories (a), (b) and (c) (Refs. 35, 36, 153 respectively) have in common a strategy consisting of the following phases: (i) time-dependent Green's functions (fundamental solutions) were taken from the literature or ad hoc constructed for the homogeneous space Ω∞ (of a linear conductive, viscoelastic or poroelastic medium, respectively) considering the effects of not only single-layer but also double-layer sources; (ii) using both these kinds of sources on suitably chosen complementary parts of the boundary of the body Ω, BIEs are generated such that they can be interpreted as a (linear) integral operator (in space and convolutive in time) which transforms the unknown fields on Γ over a given time interval T , say y(ξ, τ ), into known data-capturing fields, say g(x, t), x ∈ Γ, 0 ≤ t ≤ T ; (iii) the integral operator L is shown to be self-adjoint with respect to a bilinear form, convolutive in time, defined over Γ × T ; (iv) the above bilinear form and symmetry property of L provide a quadratic functional (convolutive in time) of y, the stationarity of which characterizes the solution over Γ × T ; (v) field modelling (under suitable continuity constraints dictated by the kernel singularities in space) and its substitution into the above functional, through the variational property (iv) and by four integrations (two in space and two convolutive in time), leads to a linear algebraic equation system, which governs the modelled approximate solution over Γ × T and is endowed with a symmetric matrix of coefficients.

The same strategy outlined in what precedes will be described with some formal details in Sec. 7 with reference to linear elastodynamics and was pointed out for elastostatics in Sec. 2.2, comment (a). Therefore, only the supplementary remarks which follow appear to be appropriate here.

(A) The line of thought centered on the generation of a symmetric integral operator over Γ × T and on the consequent variational theorem, represents a noteworthy unifying framework (alternative to the Galerkin weighted-residual approach) to formulations of SGBEM in a variety of problems, especially viable in those where energy concepts are not as natural and perceptible as in elasticity and plasticity. The methodology has origin in classical works (by Gurtin, Ref. 74, Tonti, Ref. 162 and others) on linear mechanical initial-boundary value problems formulated in terms of partial differential equations.

(B) A variational characterization of the boundary solution over a time interval can be regarded as a central ingredient in the development of SGBEM for time-dependent problems. Further progress can be achieved in this direction by generating a variational characterization in terms of saddle point with separation of variables, similar to that easily proved in elastostatics (Secs. 2.2 and 7; Refs. [START_REF] Polizzotto | An energy approach to the boundary element method; Part I: Elastic solids[END_REF][START_REF] Sirtori | A Galerkin symmetric boundary element method in elasticity: formulation and implementation[END_REF], namely with a minimum with respect to the static variables and a maximum with respect to the kinematic ones. For transient heat conduction and its analogues, problems of kind (a), this result (emerging in elasticity as a straightforward corollary of the variational properties (A)) has been established in Ref. [START_REF] Carini | Symmetric boundary integral formulations of transient heat conduction: saddle-point theorems for BE analysis and BE-FE coupling[END_REF]. In fact, the boundary solution over the unbounded time interval 0 ≤ t ≤ ∞ has been shown to be characterized by a minimum with respect to the time-histories of fluxes on the Dirichlet boundary and by a maximum with respect to the time-histories of temperatures on the Neumann boundary, when the functional involved is constructed based on a special weighted bilinear form in the Laplace transforms of the unknown boundary fields according to Rafalski's (1969) 'orthogonal projection method' for linear partial differential equations.

(C) In the problems listed at the beginning of this Section, the doublelayer, or 'kinematic', sources needed for the SGBEM consist of discontinuities, concentrated (i.e. Dirac-modelled) in space and time, of the following quantities, respectively: (a) temperatures, (b) displacements, (c) displacement and pressure, (d) acoustic pressure or velocity potential (while the more familiar single-layer sources concern: heat flux, tractions, tractions and fluid flux, normal velocity, respectively). Explicit expressions (concerning 2D and/or 3D 'free space' Ω∞) could be found in the earlier literature only for some of these kernels. Other ad hoc expressions, for isotropic media, have (D) For the classes (a), (b) and (c) of time-dependent problems, the SGBEM has been formulated and preliminarly discussed from some standpoints of computational interest (e.g.: marching solution procedures, Ref. 36, BE-FE coupling, Ref. 37). However, engineeringoriented implementations and comparative numerical tests of the SGBEM (see Sec. 12) are still missing in these areas, and so are, to the writers' knowledge, investigations of relevant computational theoretical issues such as algorithmic stability typically important in time-dependent problems.

(E) In linear acoustics and wave propagation the SGBEM has attracted since several years a considerable attention and applicationoriented effort research effort, e.g. Refs. 76, 77, even if unfortunately with limited crossfertilization among neighbouring areas.

OTHER ANALYSIS PROBLEMS

Limit and shakedown analysis, deformation bounds

Solving a limit analysis problem requires consideration of either equilibrated stress states with the static approach, or compatible plastic strain fields with the kinematic approach, whereas the elastic properties of the material behavior turn out to be irrelevant, i.e. limit analysis exhibits the nature of a rigid-plastic theory. Therefore, the boundary integral equation method (BIEM) and the relevant discretized forms as the BEM seem to be unsuitable to solve limit analysis problems, since the BIEM is deeply rooted in the material elastic behaviour through the fundamental solutions and essentially operates by superposing elastic states, i.e. complete elastic solutions to suitable auxiliary problems. As pointed out in Refs 110,[START_REF] Maier | On Shakedown Analysis by Boundary Elements[END_REF], this apparent drawback can be overcome by expressing any equilibrated stress field as the superposition of two (elastic) stress fields, one of which is the BIEM response to the pertinent loads and the other is a self-stress field obtained as the BIEM response to an arbitrary field of initial plastic strains; analogously, a compatible plastic strain field is one whose BIEM stress response identically vanishes. From the same standpoint, the BIEM appears to be applicable to shakedown analysis problems, which also require consideration of analogous stress and strain fields. In fact, limit analysis can be regarded as a special case of shakedown analysis obtainable by considering a proportional loading instead of an amplification of a domain of variable repeated loads.

Another difficulty, pointed out in Refs 110, 111, 135, 136, arises when the BIE is approximated by the conventional direct BEM, which does not preserve self-adjointness and leads to nonsymmetric solving equation systems. This implies that the BEM stress response to any initial plastic strain field is not necessarily a self-stress field and that, as a consequence, the identically vanishing of the stress response does not guarantee that the related initial strain field is compatible. For these reasons, the conventional direct BEM is not well suited to limit and shakedown analysis, whereas the SGBEM is. It has been proven in Refs 135,[START_REF] Polizzotto | An energy approach to the boundary element method; Part II: Elastic-plastic solids[END_REF] that the SGBEM stress response to any initial plastic strain field is a self-stress field and that any initial plastic strain field whose SGBEM stress response vanishes identically is compatible; in other words, the SGBEM does provide a means to construct self-stresses and compatible initial strains, with a degree of redundancy depending on the mesh, just like the FEM.

The conventional direct BEM was employed in Ref 110 to address limit analysis and in Ref 111 to address shakedown analysis in perfect plasticity. The validity of nonsymmetric formulations rests on the assumption that the mentioned deficiencies of the conven-tional direct BEM can be interpreted as modelling errors and thus can be reduced to within acceptable limits by adopting suitable discretizations. The safety factor formulated by SGBEM or traditional BEM exhibits the same mathematical programming setting (in particular, linear programming if the yield surface is piecewise linear) as in the FEM, but with the following peculiarities: in the static approach, the self-stress variables required by the static (lower bound) theorem are expressed in terms of initial plastic strain variables; in the kinematic approach, the compatible initial plastic strain variables required by the kinematic (upper bound) theorem are replaced by stress-free strain variables.

A discussion of the shakedown analysis problems for perfect plasticity within the framework of BIEM is given in Ref 124 where the classic Melan's and Koiter's theorems, Ref 114, are reproposed with the pertinent BIEM language, and their duality relationship assessed. The SGBEM is used in association with a consistent domain discretization cells aimed at the interpolation of the plastic flow laws by means of the maximum plastic work theorem (see also Refs 97,[START_REF] Polizzotto | A consistent formulation of the BEM within elastoplasticity[END_REF][START_REF] Polizzotto | An energy approach to the boundary element method; Part I: Elastic solids[END_REF][START_REF] Polizzotto | An energy approach to the boundary element method; Part II: Elastic-plastic solids[END_REF].

A related topic where BIE methods can be usefully applied is concerned with the deformation bounding techniques, Ref. [START_REF] Polizzotto | A unified treatment of shakedown theory and related bounding techniques[END_REF]. These are analytical/numerical procedures which are able to provide useful information about the actual plastic deformations produced in a body subjected to some load history, without performing a cumbersome step-by-step analysis. Deformation bounds may be a necessary complement of shakedown analysis in consideration that this analysis gives no information regarding the plastic deformation occurred in the transient elastic-plastic response preceding elastic adaptation. An early attempt to formulate bounding techniques by BIE was presented in Ref 105, where symmetry is shown to be a pre-requisite for any successful formulation of this type, and procedures of 'forced symmetrization' were proposed in order to make the conventional direct BEM applicable to this aim (Ref. [START_REF] Maier | On elastoplastic analysis by boundary elements[END_REF]). The SGBEM was employed in Ref 132 to show that the so-called 'pertubation method" previously devised in Ref 131 can be given a BIE format. A general bounding principle was there formulated which is able, in principle, to provide bounds to several types of plastic deformations, e.g. plastic strains, plastic dissipation, generalized residual displacements at the boundary, and the like. No numerical applications were presented in the papers quoted above.

Elastic analysis of plates

Nonsymmetric boundary integral formulations for plate bending, based on the use of the fundamental solution for biharmonic problems, have been studied by various authors, see the book Ref. [START_REF]Boundary Element Analysis of Plates and Shells[END_REF] or the relevant chapters in Refs. 13, 31. They consist of two coupled boundary integral equations, one strongly singular and the other hypersingular (i.e. containing singularities of order 1/r and 1/r 2 , respectively), associated with the representation formula of the flexural displacement and its normal derivative, respectively, at a point of the boundary. The traditional BE approach rests, as usual, on collocating those integral equations at a finite number of boundary points.

On the other hand, not much effort has been directed to SGBIE formulations. Tottenham, Ref. 163, outlines how a symmetric formulation can be obtained by weighting in a Galerkin sense suitable equations containing different kind of sources: forces, moments, bicouples and tri-couples. Singular integrations are next performed analytically under restrictive hypotheses on geometry and field modelling. More recently, Galerkin BIE formulations for plates with free edges have been derived in Refs. 66, 119, the latter allowing for the presence of corners.

In a recent work, Ref. 59, a direct symmetric boundary integral formulation has been sought in the form of the stationarity condi-tions of the Lagrangian functional obtained by incorporating kinematic boundary conditions, in the form of constraint terms, into the potential energy associated with bending when first-order variations of the unknown bending displacement solve the homogeneous elastic equilibrium equation. In deriving such a formulation one has to deal with very high potential singularities on the boundary, of order up to 1/r 4 ; this fact may explain the past lack of research effort in this area. The evaluation of these terms is tackled using integrations by parts, thanks to the fact that the most singular kernels are shown to be derivatives of other, less singular, kernels with respect to the arc length along the boundary, in an approach similar in principle to that developed in Sec. 8.1. This result in a SGBIE formulation which is direct, i.e. in terms of mechanical unknowns on the boundary (bending displacement, normal slope, normal moment and Kirchhoff shear), accommodates general boundary conditions, and include unknown jumps of twisting moment at corners. For the case of free edges, this formulation and that of Ref. 66 are identical. The singularities involved are at most logarithmic, like in two-dimensional elasticity after regularization, so that the numerical quadrature techniques briefly discussed in Sec. 12 are applicable. This formulation has been implemented with either straight or circular boundary elements, and successfully tested against several analytical solutions. Finally, a BE expression of the stiffness matrix of an unloaded plate has been derived following a method similar to that presented in Sec. 11 for 2-D and 3-D elasticity.

Only very few practical applications of SGBEM have yet been published. One of the initial motivations for the aforementioned work was the study of cracks and delamination, using energy concepts like energy release rate. Like in classical linear fracture mechanics, Secs. 4 and 10, a domain derivative approach to the computation of energy release rate in cracked or delaminated plates crucially benefits from the symmetry of the SGBIE formulation. In a first attempt in this direction, Ref. 61, encouraging numerical results have been achieved.

VARIATIONAL AND ENERGY APPROACHES

The energy principles of mechanics can be used through suitable variational procedures to characterize the system's response to given external actions. These principles, together with the related variational procedures, have a remarkable role in the framework of the discretization methods in order to give firm basis to the discretization operations; that is, in order to derive, from a given continuum a discrete system with qualitatively the same essential features. In this way, the discrete solution is not only an approximate solution of the original problem, but also the true solution for the discretized mechanical system. Such an approach to the continuum solution problem has several beneficial consequences from both the numerical and theoretical points of view. Namely, on one hand, error and convergence analyses are in general easier to accomplish and convergence of the numerical solution to the exact one is likely to occur even with a greater convergence rate; on the other hand, the self-adjointness of the original problem is preserved and the discretized equation set can be employed as an effective analytical basis for theoretical investigations within the framework of discrete mechanics.

The Finite Element Method (FEM) is potentially susceptible to possess the above mentioned requisites of the energy-based discretization methods, but the same cannot be stated for the conventional BEM. The first attempt to formulate an energy-based BEM can apparently be attributed to Reissner, who proposed in Ref. [START_REF] Reissner | On some variational theorems in elasticity[END_REF] All these formulations, often referred to as hybrid BEMs (HBEM), are characterized by symmetry and require single integrations over the boundary. They are thus quite interesting to numerical solution purposes, but have not been shown to possess the desirable requisites of a genuine energy-based discretization method.

As mentioned in Sec. 2.1, kernels Guu and Gpp can be associated with strain energies. Specifically, consider the displacement field u generated in the elastic space Ω∞ = I R 3 by a single-layer source f :

u(x) = Γ Guu(x, ξ)f (ξ) dΓ ξ ( 76
)
where Γ is any closed oriented surface of I R 3 . Then, the strain energy W∞ for the whole Ω∞ is given by:

W∞ = 1 2 Γ (p T (x -)u(x -) -p T (x + )u(x + )) dΓx (77)
having taken into account that the displacement field (76) (i) solves the local elastic equilibrium equation and (ii) verifies radiation conditions at infinity. Hence, using Eq. ( 4) and the continuity across Γ of the field (76), one can write:

W∞ = 1 2 Γ Γ f T (x)Guu(x, ξ)f (ξ) dΓ ξ dΓx (78)
Similarly, considering the displacement field u generated in Ω∞ = I R 3 by double-layer sources d:

u(x) = Γ Gup(x, ξ)d(ξ) dΓ ξ (79)
Eq. ( 5) and computing p(x ± ) from ( 79), W∞ still given by Eq. ( 77) becomes:

W∞ = - 1 2 Γ Γ d T (x)Gpp(x, ξ)d(ξ) dΓ ξ dΓx (80)
Hence, the interpretation of kernels Guu and Gpp in terms of (positive) strain energies clearly entails the sign-definiteness of the relevant quadratic forms, Eqs. ( 9), [START_REF] Balakrishna | Efficient analytical integration of symmetric Galerkin boundary integrals over curved elements: thermal conduction formulation[END_REF]. The kernels Guu, Gpp appear also in the expression of the strain energy stored in a bounded elastic solid (see Sec. 11). In Sec. 2, SGBEM elastostatic formulations have been discussed from a weighted-residual viewpoint. The same final formulation can be obtained from variational principles, as explained in what follows. Hu-Washizu and Hellinger-Reissner principles of linear elasticity were used in Refs. 131-137 to derive energy-based BEM formulations which happen to coincide with the symmetric Galerkin ones, see Refs. 112, 155. With such formulations, the BE discretized model turns out to be well defined, so that the related discrete equation set has a unique solution in terms of generalized boundary displacements and tractions. For given generalized boundary data (either diplacements or tractions at every boundary node), there is a unique boundary response in terms of generalized variables (either tractions, or displacements). This set of generalized boundary variables is representative of a class of continuum solutions that are equivalent to one another, in the sense that the boundary displacements and tractions of every solution have different distributions, but result into the same (unique) generalized boundary values, respectively. Following similar lines of thought, continuous SGBIE formulations for three-dimensional elastostatics (Ref. [START_REF] Bonnet | Regularized direct and indirect symmetric variational BIE formulations for three-dimensional elasticity[END_REF]) and Kirchhoff plate theory (Ref. 59) are shown to express the stationarity conditions for an augmented elastic potential energy in which the kinematical boundary data appear as constraints.

The same BEM formulation can also be derived using another variational principle, that is, the so-called min-max principle proposed in Refs. 133, 138. For an elastic body subjected to volume forces b in Ω, imposed displacements ū on Γu and tractions p on Γp, let us consider the functional

Π[p, u] = 1 2 Γu Γu p T (x)Guu(x, ξ)p(ξ) dΓx dΓ ξ - Γu Γp p T (x)Gup(x, ξ)u(ξ) dΓx dΓ ξ + 1 2 Γp Γp u T (x)Gpp(x, ξ)u(ξ) dΓx dΓ ξ - Γu p T ḡu dΓx + Γp u T ḡp dΓx ( 81
)
where ḡu (x ∈ Γu) and ḡp (x ∈ Γp) are functions of the data, i.e.:

ḡu (x) = 1 2 ū(x) - Γp Guu(x, ξ)p(ξ) dΓ ξ - Ω Guu(x, ξ) b(ξ) dΩ ξ + - Γu Gup(x, ξ)ū(ξ) dΓ ξ ḡp (x) = 1 2 p(x) -- Γp Gpu(x, ξ)p(ξ) dΓ ξ - Ω Gpu(x, ξ) b(ξ) dΩ ξ + = Γu Gpp(x, ξ)ū(ξ) dΓ ξ
It can be proved that the solution in a suitable function set to the saddle-point problem:

min p max u Π[p, u] (82)
provides the tractions p = p|Γ u and the displacements u = u|Γ p pertaining to the solution of the boundary-value problem, and that conversely the latter solution solves problem [START_REF] Holzer | A p-version of the symmetric BEM[END_REF]. The proof rests on two requisites of the functional ( 81): (i) the symmetry properties (8) of the kernels Grs; and (ii) the positive and negative definiteness of, respectively, the first and the third double integral of [START_REF] Holzer | A p-version of the symmetric BEM[END_REF]. Modelling p as p h = Ψp(x)Xp on Γu and u as u h = Ψu(x)Xu on Γp, one has

Π[p h , u h ] = Π(P , U ) ≡ 1 2 P T Ĝpp pp P -P T Ĝpu pu U + 1 2 T Ĝuu uu U -P T Rp + U T Ru ( 83
)
where Ĝij rs are defined by Eq. ( 20) and

Ru = Γp Ψ T u ḡp dΓ Rp = Γu Ψ T p ḡu
Π(P , U ) (84)
is the discrete counterpart of [START_REF] Holzer | A p-version of the symmetric BEM[END_REF]. The relevant Kuhn-Tucker conditions read:

W = 1 2 Γ Γ p T (x)Guu(x, ξ)p(ξ) dΓx dΓ ξ - 1 2 Γ Γ u T (x)Gpp(x, ξ)u(ξ) dΓx dΓ ξ (86)
Then, the following augmented potential energy L is introduced:

L =W - Γ pT (x)u(x) dΓx + Γ Γ λ T (x) Gup(x, ξ)u(ξ) -Guu(x, ξ)p(ξ) dΓx dΓ ξ
where λ(ξ) is the Lagrange multiplier field associated with the boundary compatibility constraint (i.e. the fact that u, p are the traces on Γ of an elastostatic state u, σ in Ω). Then, by expressing the stationarity of L, one finds that (i) λ = p and (ii) the unknown parts p|Γ u and u|Γ p solve the now usual SGBIE formulation. The above interpretation is an a posteriori one in the sense that the weighted residual SGBIEs are used (with particular choices of the test functions) to establish Eq. ( 86), see Sec. 11. Energy based variational formulations were also achieved for elastodynamics in Ref. [START_REF] Polizzotto | A variational approach to boundary element methods[END_REF]. The Hu-Washizu principle extended to dynamics is there employed together with the classical Stokes fundamental solutions. Through semidiscretization by boundary elements, the original initial-boundary value problem is transformed into a discrete set of Volterra time-integral equations of the first kind, which can be numerically solved through a suitable time-stepping algorithm. The potentialities of this formulation, though already extended to plasticity in Ref. [START_REF] Panzeca | A BEM formulation of the dynamic elastic-plastic structural problem via variational principles[END_REF]

y(x, t) = p(x|Γ u , t) u(x|Γ p , t) , ḡ(x, t) = ḡu (x|Γ u , t) ḡp (x|Γ p , t) (87)
where ḡu and ḡp depend on the assigned loading and initial conditions like in Eq. ( 59) but with θ = 0, that is:

ḡu = ū -φu , ḡp = p -φp (88)
where

φh = Ω G hu b dΩ + Γp G hu p dΓx - Γu G hp ū dΓx -ρ Ω G hu | (t=0) u0 + Ġhu | (t=0) u0 dΓ ξ (h = u, p)
denoting by u0 and u0 the given initial displacements and velocities.

The above mentioned boundary operator is defined as follows:

Ly ≡              t 0 Γu Guu(x|Γ u , ξ; t -τ )p(ξ, τ ) dΓ ξ dτ - t 0 Γp Gup(x|Γ u , ξ; t -τ )u(ξ, τ ) dΓ ξ dτ t 0 Γu Gpu(x|Γ p , ξ; t -τ )p(ξ, τ ) dΓ ξ dτ - t 0 Γp Gpp(x|Γ p , ξ; t -τ )u(ξ, τ ) dΓ ξ dτ              . ( 89
)
This operator turns out to be symmetric (or self-adjoint) both in space and time with respect to the time-convolutive bilinear form:

L y, y ′ = T 0 Γ y ′ T (x, T -t) L y (x, t) dΓx dt, ( 90
)
such that the equality

L y, y ′ = L y ′ , y (91)
holds for arbitrary y and y ′ . As a consequence, the unknown boundary displacements and tractions pertaining to the dynamic response of the body to given loading and initial conditions can be obtained as the displacements and tractions that make stationary the functional

F (p|Γ u , u|Γ p ) ≡ 1 2 L y, y -y, ḡ . (92)
Through discretization by boundary elements with interpolations as

p h (x; t) = Ψp(x; t)P on Γu × T u h (x; t) = Ψu(x; t)U on Γp × T ( 93
)
where T denotes the time interval of interest, like in Sec.

REGULARIZED SGBE FORMULATIONS

Regularization

The various fundamental solutions involved in SGBEM are singular when x = ξ. In linear elastostatics, their singularities have been specified by Eqs. ( 6) and [START_REF] Arnold | On the asymptotic convergence of collocation methods[END_REF]. For other classes of problems (potential theory, acoustics, elastodynamics, diffusion. . . ), the singular behavior of the kernels is similar. In particular it is well known that corresponding dynamic and static fundamental solutions have the same singularity (i.e. their difference is nonsingular). From Eqs. ( 6)- [START_REF] Arnold | On the asymptotic convergence of collocation methods[END_REF], it is apparent that the various Galerkin BIE formulations encountered in the previous sections cannot in general be implemented in a straightforward way. In particular, the integral of the so-called hypersingular kernel Gpp(ξ, x) over a single surface Γ can be defined only in the sense of Hadamard finite part or other similar limiting processes, which require C 1,α continuity of the source density function at the singular point. Its integral over the Cartesian product Γ × Γ is convergent for 3D, and still divergent for 2D. Moreover the strongly singular kernels Gup, Gpu are integrable over Γ only in the sense of the Cauchy principal value or other similar limiting processes, which require a C 0,α behavior of the density function at the singular point. Henceforth, C 0,α (Γ) is the set of functions f such that

∃α, 0 < α ≤ 1, ∃C > 0, ∀(x, y) ∈ Γ × Γ, |f (y) -f (x)| < C |y -x| α
and C 1,α (Γ) is the set of functions having C 0,α first-order derivatives.

In view of the above circumstances, it is convenient to rearrange somehow the SGBIE formulation in order to eliminate all difficulties related to integrability before performing the actual discretization. This step is usually known as regularization. This issue has been addressed in several ways.

The symmetric formulation of Balakrishna et al., Ref. 9, consists of using weighted residuals of regularized displacement and traction BIE; the integrability of kernel Gpp is then settled at the first (inner) integral level, at the expense of having to rely on C 1,α behavior for the unknown displacement. The somewhat involved resulting formulation otherwise allows for 2D and 3D implementation using any usual interpolation scheme.

In the 2D implementation of Sirtori et al., Ref. 154, the relevant double integrals are temporarily taken on two close but distinct curves Γ and Γ ⋆ , where Γ ⋆ surrounds Γ and is in one-to-one correspondence with Γ; then a piecewise linear discretization is introduced for Γ, Γ ⋆ and the unknown boundary variables and test functions. This particular choice of discretization allows to evaluate analytically all integrals while the discretized Γ h , Γ ⋆ h are kept distinct; then, the limit Γ ⋆ → Γ is taken analytically. In Ref. 154, a complex variable formalism was found convenient, while in Ref. [START_REF] Maier | On elastoplastic analysis by boundary elements[END_REF], where this procedure was extended to elastic-plastic 2D analysis, real variables have been preferred in view of future 3D implementations.

A recent generalization of the limiting approach of Ref. 154 and 98, applicable to either two-or three-dimensional situations, allows to formulate regularized SGBIE in full generality; in particular threedimensional situations fall within this framework. Essentially, the regularization technique consists of a limiting process Γ ⋆ → Γ carried out before any discretization. This is made possible by a key property of the Gpp kernel, which permits its reformulation as a repeated surface curl. Then,Stokes' formula is applied with respect to both x and ξ, resulting in the appearance of a new, weakly singular, kernel. This approach is presented below for elastic problems without body forces, following Ref. [START_REF] Bonnet | Regularized direct and indirect symmetric variational BIE formulations for three-dimensional elasticity[END_REF]. The SGBIE formulation comes e.g. from applying the Betti reciprocity identity and

lim x→Γ ⋆ Γ ⋆ Gup(x, ξ)d ⋆ (ξ) dΓ ξ = - 1 2 d ⋆ (x) + - Γ ⋆ Gup(x, ξ)d ⋆ (ξ) dΓ ξ (98)
denoting by the symbol -Γ a Cauchy principal value (CPV) integral. It follows that, in the limit Γ ⋆ → Γ: lim

Γ ⋆ →Γ Γ ⋆ Γ u T (x)Gpu(x, ξ)f ⋆ (ξ) dΓ ξ dΓx = 1 2 Γ u T f ⋆ dΓx + Γ - Γ u T (x)Gpu(x, ξ) dΓ ξ f ⋆ (ξ) dΓx ( 99
)
and

lim Γ ⋆ →Γ Γ Γ ⋆ p T (x)Gup(x, ξ)d ⋆ (ξ) dΓx dΓ ξ = - 1 2 Γ p T d ⋆ dΓx + Γ p T (x) - Γ ⋆ Gup(x, ξ)d ⋆ (ξ) dΓx dΓ ξ ( 100
)
At this point, the inner CPV integrals can be either evaluated by a direct method (Guiggiani and Gigante, Ref. 73), or regularized using the 'rigid-body' identity

- S Gpu(x, ξ) dΓ ξ = (κ(ξ) - 1 2 )1 (ξ ∈ S) (101)
where S is a closed surface, 1 denotes identity matrix and κ(ξ) = 1 or 0 depending on whether the unit normal n to S (implicitly present in the kernel Gpu(x, ξ)) is interior or exterior to the surface S. For instance: (i) if the domain boundary Γ is connected, one has κ = 1 for an exterior problem and κ = 0 for an interior problem; (ii) if the domain boundary Γ is made of two or more disjoint closed components, both cases κ = 1 or 0 arise. Besides, in 2D situations, the kernels Gup, Gpu can be recast as derivatives of other kernels with respect to the arc-length, as shown in Ref. [START_REF] Ghosh | A new boundary element method formulation for linear elasticity[END_REF]; hence regularization through integration by parts is available in this case, and has indeed been used in Refs. 60, 104, 154.

Hypersingular kernel

Turning now to the most singular kernel, Gpp, recall that the Kelvin displacement kernel Guu is such that

Guu(x, ξ) = Guu(ξ, x) = Guu(ξ -x) (102)
which in turn implies that

∇xGuu(x, ξ) = -∇ ξ Guu(x, ξ) = ∇Guu(ξ -x) (103)
where ∇f (ξx) denotes the gradient of f (z) taken at z = ξx. Indeed, from the Galerkin representation formula for solutions in isotropic elasticity (with shear modulus µ and Poisson ratio ν), one has

Guu = 2(1 -ν)∇ 2 F 1 -∇∇F (104)
where 1 denotes the identity matrix in I R 3 and the scalar function F (ξx) for 3D and 2D situations, respectively, reads:

F (ξ -x) = 1 16πµ(1 -ν) r F (ξ -x) = -1 16πµ(1 -ν)
r 2 ln r and are biharmonic: ∇ 4 F = 0 for any ξ = x. The ik-component of Gpp is then given by

[Gpp(x, ξ)] ik = -C ijab C kℓcd [Guu(x, ξ)] ac,bd nj (x)n ℓ (ξ) = Z ijkℓ (x, ξ)nj(x)n ℓ (ξ)
where C is the elastic tensor and f ,b (ξx) denotes the partial derivative ∂ ∂z b f (z) taken at z = ξx. Inserting Eq. ( 104) in the above formula, one obtains the expression of the Z ijkℓ (x, ξ) components in terms of F as follows:

1 µ 2 Z ijkℓ = -4ν[F ,ijkℓ -δijF ,ppkℓ -δ kℓ F,ppij ] -2(1 -ν)[2F ,ijkℓ -δ ik F ,ppjℓ -δ jℓ F ,ppik -δ iℓ F ,ppjk -δ jk F ,ppiℓ ] (105)
Then, δ being the Kronecker symbol and e the Ricci permutation tensor, using the algebraic identity eaepe cf p = δacδ efδ af δ be , one obtains eiepe jf q e kgr e ℓhs δpqδrsF ,ef gh

= F ,ijkℓ -δij F ,ppkℓ -δ kℓ F,ppij + δij δ kℓ ∇ 4 F
and similar relations under permutations of i, j, k, ℓ. Finally, since ∇ 4 F = 0 for any ξ = x, substitution of the above identities into Eq. ( 105) leads to the following expression of Gpp:

[Gpp(x, ξ)] ik = R ξ q R x s B iqks (x, ξ) (106)
in terms of the surface curl defined as Rcf (ξ) ≡ e abc na(ξ) ∂f ∂x b [START_REF] Maier | Symmetric Galerkin boundary element method for quasi-brittle fracture and frictional contact problems[END_REF] and in terms of the new kernel

B iqks = eiepe kgr µ 2 [4νδpqδrs + 2(1 -ν)(δprδqs + δpsδqr)] ∂ 2 F ∂ξe∂xg (108)
This kernel turns out to be only weakly singular (O(1/r) for 3D and O(ln r) for 2D problems). The surface curl, Eq. ( 107), is associated to the following form of the Stokes' formula:

S Rcf dS = ∂S f τc ds ( 109
)
where S is any regular surface and τ is the unit tangent to its edge ∂S, if present (i.e. if S is not closed). Thus, application of the above identity together with (106) leads to

Γ Γ ⋆ u T (x)Gpp(x, ξ)d ⋆ (ξ) dΓ ξ dΓx = Γ Γ ⋆ Rqui(x)B iqks (x, ξ)Rsd ⋆ k (ξ) dΓ ξ dΓx (110)
Note that no contour term appears since Γ and Γ ⋆ are closed surfaces and both u and d ⋆ are continuous. At this point, since the kernel B iqks is weakly singular, Eq. ( 110) is valid in the limit Γ ⋆ → Γ.

For two-dimensional problems, the above considerations about the hypersingular kernel Gpp still apply, with some simplifications. Let e3 be the direction of invariance, and define the unit tangent τ to the boundary curve Γ by τ = e3 ∧ n, or τ b = e ab3 na. Thus, Rcf , Eq. ( 107), reduces to δc3d/ds (s: arc length) and

[Gpp(x, ξ)] ik = d ds ξ d dsx B ik (x, ξ) B ik (x, ξ) = µ 2π(1 -ν) [ln rδ ik -r,ir ,k ] (111)
The above expression is used in Ref. 60.

To recapitulate, the final SGBIE formulation, ready for implementation, is obtained through the following steps: (i) substitute Eqs. [START_REF] Li | Symmetric weak-form integral equation method for three-dimensional fracture analysis[END_REF][START_REF] Liu | Plastic analysis by quasi-high order Galerkin BEM[END_REF], with f ⋆ = 0 on Γ ⋆ p and d ⋆ = 0 on Γ ⋆ u , into (94); (ii) take the limit expressions given by Eqs. ( 99), ( 100) and ( 110); (iii) split Γ into Γu, Γp and separate knowns ū, p and unknowns u, p.

It is worth noting that in general the unknown displacement u and the kinematical trial function d ⋆ , both defined on Γp, do not belong to the same function space, since they must realize continuous extensions on Γ of ū and the null function on Γp, respectively. Then, one way to guarantee a symmetric final formulation is to introduce a function û chosen so that û = ū on Γu, continuous on Γ and otherwise arbitrary, and define on Γp the new unknown v = uû; then both v and d ⋆ are continuous and vanish on Γu. For pure Neumann boundary conditions, the introduction of û, v is not necessary.

Assume that the boundary Γ, boundary variables u, p and trial functions d ⋆ , f ⋆ are discretized in the usual BEM fashion. Following the notations of Eqs. ( 18), [START_REF] Belytschko | Singular integration in variationally coupled FE-BE method[END_REF] of Sec. 2, û and v can be simply defined as the interpolations of the prescribed displacement nodal values and of the unknown ones, respectively, namely:

û = Ψu[0 T , Ū T ] T v = Ψu[X T u , 0 T ]
The final outcome of the regularization procedure is, of course, the linear system of equations [START_REF]Boundary Element Analysis of Plates and Shells[END_REF], i.e.

pp = Γp Γp RΨ T u (x)Bpp(x, ξ)RΨu(ξ) dΓx dΓ ξ ( 112
) Ĝup up = Γp Γu Ψ T u (x)Gpu(x, ξ)Ψp(ξ) dΓx dΓ ξ ( 113
) Ĝpu pu = Γu Γp Ψ T p (x)Gup(x, ξ)Ψu(ξ) dΓx dΓ ξ ( 114
) Ĝuu uu = Γu Γu Ψ T p (x)Guu(x, ξ)Ψp(ξ) dΓx dΓ ξ (115)
and the subvectors by

B f u = Γu κ(ξ)ū T (ξ)Ψp(ξ) dΓ ξ + Γ Γu [û T (x) -ûT (ξ)]Gpu(x, ξ)Ψp(ξ) dΓx dΓ ξ - Γp Γu pT (x)Guu(x, ξ)Ψu(ξ) dΓx dΓ ξ (116)
and

B d p = Γp (κ(ξ) -1)p T (ξ)Ψu(ξ) dΓ ξ - Γ Γp R ûT (x)Bpp(x, ξ)RΨu(ξ) dΓx dΓ ξ + Γp Γp pT (x)Gup(x, ξ)[Ψu(ξ) -Ψu(x)] dΓx dΓ ξ - Γp p(x) T Ψu(x) Γu Gup(x, ξ) dΓx dΓ ξ (117)
where κ(ξ) has the same meaning as in Eq. ( 101).

All single inner surface integrals in the above formulas are at most weakly singular and, hence, have a meaning in the usual Riemann sense. In particular, the integrals over Γp of Ĝup up and Ĝpu pu are at most weakly singular by virtue of the fact that the relevant interpolation functions Ψu are continuous and vanish along the separation curve between Γu and Γp.

It is worth stressing that Eqs. ( 112) to [START_REF] Nappi | An internal variable approach applied to boundary element elastic-plastic analysis[END_REF] hold for interior as well as exterior problems, and that they can accommodate nonconnected boundaries. They provide in all respects an adequate basis for a purely numerical implementation, i.e. relying solely on numerical quadrature techniques. This fact is very valuable for generalpurpose implementations to come, especially for three-dimensional SGBEM algorithms, because it enables one to use curved elements and high-degree interpolations, like in traditional collocation BEMs. Of course, the purely numerical approach to SGBEM implementation is by no means the only correct one, and various other approaches, mixing numerical and analytical quadrature, have been successfully developed, see e.g. Refs. 79 and 81.

Indirect SGBEM formulation

Following classical methods of potential theory, let Ω be bounded and introduce its unbounded complement Ω + = Ω∞ \ Ω. Consider an interior and an exterior problem simultaneously, with equal displacement data on Γu and opposite traction data on Γp: ū+ = ū and p+ + p = 0. Combining the SGBEM formulations for both problems, using Eqs. [START_REF] Maier | A Galerkin approach to boundary element elastoplastic analysis[END_REF] to [START_REF] Nappi | An internal variable approach applied to boundary element elastic-plastic analysis[END_REF], yields a new formulation, named indirect, in terms of new nodal unknowns Xϕ = X + u -Xu and Xπ = X + p + Xp. The corresponding linear system of equations is: Ĝuu

uu - Ĝup up - Ĝpu pu Ĝpp pp X π X ϕ = Lu u Lp p (118)

Crack problems

The regularization approach for cracked bodies follows the same steps as in Sec. 8.1, but with the Betti theorem and fictitious state defined by Eqs. ( 45), ( 46), [START_REF] Crouch | Boundary element methods in solid mechanics[END_REF]; an auxiliary surface S ⋆ must be introduced in order to consider the limiting situation S ⋆ → S. As a result, the new term B ψ ϕ in the right-hand side of ( 48) is given by

B ψ ϕ = S pT (ξ)Ψu(ξ) dΓx + Γ S R ûT (x)Bpp(x, ξ)RΨϕ(ξ) dΓx dΓ ξ - Γp S pT (x)Gup(x, ξ)Ψϕ(ξ) dΓx dΓ ξ (119)
whereas B d u , B p p are still given by Eqs. ( 116), [START_REF] Nappi | An internal variable approach applied to boundary element elastic-plastic analysis[END_REF].

Further comments on the regularization procedure

The kernel decomposition ( 106) is given by Nedelec in Ref. [START_REF] Nedelec | Integral equations with non integrable kernels[END_REF], see also Han,Ref. 78. For a scalar fundamental solution G of Laplace equation, such that ∆xG + δ(ξx) = 0, it takes the simpler form

G,ij(x, ξ)ni(x)nj(ξ) = -eiepni(x)e jf p nj (ξ)G ,ef (x, ξ)
For frequency-domain elastodynamics, the representation (104) of the fundamental displacement reads for 3D problems:

Guu = 2(1 -ν)[∇ 2 F + k 2 L F ]1 -∇∇F F = 1 4πµk 2 T (e ik L r -e ik T r ) 1 r (120)
with k 2 L = ρω 2 /(λ + 2µ), k 2 T = ρω 2 /µ, λ and µ being the Lamé constants. Then, an analysis similar to that conducted for the static case yields

[Gpp(x, ξ)] ik = R ξ q R x s B iqks (x, ξ) + 2(1 -ν)k 2 L C ijab C kℓcd δacF ,bd nj (x)n ℓ (ξ) + [2(1 -ν)(δ ik δ jℓ + δ jk δ iℓ ) + 4νδijδ kℓ ]∇ 4 F nj (x)n ℓ (ξ) (121)
where B iqks takes the form [START_REF] Maier | On Symmetrization in boundary element elastic and elastic-plastic analysis[END_REF] with F given by Eq. ( 120). One can show using series expansions that the second derivatives of F have a weak O(1/r) singularity; the factor of k 2 L is hence only weakly singular. Besides, one has

∇ 4 F = 1 4πµk 2 T (k 4 L e ik L r -k 4 T e ik T r ) 1 r
which is also weakly singular. A similar, more general, result has been achieved by Bécache et al., Ref. 17 for anisotropic elastodynamic fundamental solutions, both in frequency and time domain.

The above outlined regularization approach has been recently extended to Kirchhoff plate bending in Ref. [START_REF] Frangi | A Galerkin symmetric and direct BIE method for Kirchhoff elastic plates: formulation and implementation[END_REF]. In this work, extensive use is made of integration by parts; the most singular kernel present in the formulation, initially of order O(1/r 4 ), is recast in the fashion of (106) as a fourth-order derivative with respect to arc-length of another, O(ln r)-singular, kernel.

An interesting and promising alternative to regularization is represented by the fully numerical integration techniques developed in Refs. 1, 54 and 116

CONVERGENCE AND OTHER MATHEMATICAL FEATURES

The mathematical properties of the variational integral formulations, namely symmetry and sign-definiteness, prove quite useful to address general questions such as existence and uniqueness of solution to the continuous problem and convergence study of the corresponding discretized problem. Many fundamental boundary-value problems have been analyzed this way by Nedelec and co-workers, Refs. 11, 120-122. For instance, the variational formulation for the Dirichlet problem with boundary data ū ∈ H 1/2 (Γ) using a single-layer potential of unknown density q is

     bD(q, q ′ ) = ū, q ′ ∀q ′ ∈ H -1/2 (Γ) bD(q, q ′ ) = 1 4π Γ Γ q(ξ)q ′ (x) |x -ξ| dΓ ξ dΓx (122)
and the variational formulation for the Neumann problem with boundary data q ∈ H -1/2 (Γ) using a double-layer potential of unknown density q is

     bN (u, u ′ ) = q, u ′ ∀u ′ ∈ H 1/2 (Γ) bN (u, u ′ ) = 1 4π Γ Γ Riu(ξ)R ′ i u(x) |x -ξ| dΓ ξ dΓx (123)
where a double integration by parts regularization has been conducted and Riu is the surface curl defined by Eq. ( 107).

The mathematical definition of the function space H s (Γ), s > 0 is (see e.g. Ref. 166): Thanks to symmetry and sign-definiteness, the bilinear forms bD, bN are shown to be coercive, i.e. there exists positive constants βD, βN such that bD(q, q) ≥ βD |q| 2

H s (Γ) = f , |f | 2 H s (Γ) ≡ |α|≤[s] Γ |D α f | 2 dΓ + |α|=[s] Γ Γ |D α f (x) -D α f (ξ)| 2 |x -ξ| n-1+2(s-[s]) dΓx dΓ ξ < +∞ (124)
H -1/2 (Γ) bN (u, u) ≥ βN |u| 2 H 1/2 (Γ)/IR
Then, existence and uniqueness of the solution to problems [START_REF] Nedelec | Approximation des équations intégrales par éléments finis. Etude d'erreur[END_REF] or [START_REF] Pan | A symmetric boundary equation approach to transient poroelastic analysis[END_REF] follows directly from Lax-Milgram theorem. Moreover, using Lagrangian elements P k for the unknown density and Pm for the surface Γ, the following convergence results, in terms of the mesh size parameter h, are given in Ref. 122:

Γ u(ξ) -u h (P -1 (ξ) 2 dΓ ξ ≤ C1h m + C2[h k+1 + h m+1 ] Γ q(ξ) -q h (P -1 (ξ) 2 dΓ ξ ≤ C1h m + C2h k+1 (125)
(the orthogonal projection P onto Γ is used in order to 'transport' onto the exact boundary Γ the approximate solution q h or u h defined on the approximate boundary Γ h ). They suggest to choose interpolations so that m = k + 1 (e.g. flat triangles with piecewise constant unknowns) for the Dirichlet problem, and m = k (e.g. flat triangles with piecewise linear unknowns) for the Neumann problem. Convergence results of similar nature are known for SGBIE formulation of fundamental boundary value problems associated with e.g. Helmholtz or linear elasticity equations.

A more general approach to the mathematical study of Galerkin BIE formulations is found e.g. in Wendland,Refs. 166,[START_REF] Wendland | Boundary element methods for elliptic problems[END_REF] It is based on the fact that the various integral operators v → Gv:

Gv(ξ) ≡ Γ G(ξ, x)v(x) dΓx + (free term)
that arise in boundary integral formulations for elliptic problems, are particular instances of pseudo-differential operators, i.e. they admit a decomposition of the form

(Gv)(x(τ )) = Cv + IR n-1 IR n-1 g(τ , t, η).v(x(t))dt dη (126) with g(τ , t, η) = e 2iπ(τ -t).η a0(τ , η)Ψ(|t -τ |)
and where

t ∈ I R n-1 → x(t) is a local parametrization of Γ, Ψ(ρ)
is a C ∞ cut-off function with compact support and such that Ψ = 1 in a neighbourhood of ρ = 0. The principal symbol a0(t, η) is C ∞ with respect to η; moreover it is homogeneous of degree α:

a0(t, λη) = λa0(t, η)
for a certain value of α, called the order of G. In definition [START_REF] Panzeca | Boundary/field variational principles for the elastic-plastic rate problem[END_REF],

C : H s (Γ) → H s-α (Γ)
is a compact mapping for some s near α.

The integral operator G is then said to be strongly elliptic if and only if there exists a positive constant γ and a complex matrix-valued function Θ ∈ C ∞ (Γ) such that, for all x(τ ) ∈ Γ, all η ∈ I R n-1 with |η| = 1 and all χ ∈ C n , the principal symbol verifies

Re{χ T .[Θ(τ ).a0(t, λη)]. χ} ≥ γ |χ|
The principal symbol is strongly related to the degree of singularity of the integral operator. In particular, the integral operators with elastostatic kernels Guu, Gpp are found to be strongly elliptic (with Θ = 1) of order α = -1 and α = +1 respectively (see e.g. Eqs.

(2.1.18)-(2.1.24) of Ref. 167); this is also true of the corresponding dynamic integral operators.

A general convergence result is available for boundary element discretization of Galerkin formulations associated with strongly elliptic integral operators of order α, in connection with the use of the generic family of so-called S d,r h finite element approximations systems in (n -1) dimensions (Babuska and Aziz, Ref. 8), where the polynomial degree d of the interpolation functions and the regularity degree r < d -1 are nonnegative integers, and h denotes the maximum element size. This family of interpolations includes usual C 0 Lagrange interpolations Pm (for d = m + 1, r = 1) and piecewise constant interpolations on triangles (for d = 1, r = 0). The result states that, for α ≤ 2r + 1 (n = 2) or α ≤ 2r (n = 3) and for any real numbers t, s such that α ≤ t ≤ α/2 ≤ s ≤ d, the approximation v h to the true solution v of the Galerkin formulation for the integral equation Gv = b is such that

|v h -v| H t (Γ) ≤ ch s-t |v| H s (Γ) (127)
(see Eq. ( 124) for the definition of the H s norm). Estimates [START_REF] Panzeca | A BEM formulation of the dynamic elastic-plastic structural problem via variational principles[END_REF], neglecting the effect of geometrical modeling, are particular instances of the above result. Also, Eq. ( 127) is further generalized to systems of integral operators of different orders, so should also apply to SGBIE formulations of mixed boundary-value problems, although no emphasis is put on this issue in the references quoted herein. In contrast, convergence results concerning the more widely used collocation BEM formulations are scarce. They mostly address spline discretizations for 2D formulations, Refs. 7, 166, 167.

Few results have yet appeared concerning the mathematical study and actual implementation of adaptive SGBEM methods. A notable contribution, Ref. 156, addresses hp adaptive Galerkin BEM for two-and three-dimensional Dirichlet and Neumann problems for the Laplace equation on domains having polygonal or polyhedral boundaries, including screens (i.e. open surfaces across which the primary variable is allowed to jump) in unbounded domains. In this setting, the author establishes exponential convergence (with respect to the number of DOFs) of an hp adaptive strategy, and provides numerical experiments for both two-and three-dimensional examples.

Two-dimensional formulations are addressed in Ref. [START_REF] Wendland | A posteriori local error estimates of boundary element methods with some pseudo-differential equations on closed curves[END_REF], where estimates of the error I |uu h | 2 ds for a given interval I ⊂ ∂Ω are expressed in terms of a norm of the residual R = b-Gu over any interval J with I ⊂ J ⊂ ∂Ω, thus providing the basis for a posteriori error estimation using residuals. An a posteriori local error estimator based on the computation of the Galerkin residuals obtained by reversing the boundary condition types is introduced and numerically tested for 2D potential problems in Ref. [START_REF] Paulino | Galerkin residuals for error estimation and adaptivity in the symmetric Galerkin boundary element method[END_REF]. Two-dimensional hpand p-versions of the SGBEM are expounded and their faster rate of convergence is numerically experimented in Refs. 82, 83.

SENSITIVITY ANALYSIS FOR STRUCTURAL IDENTIFICATION AND OPTIMIZATION

Sensitivity analysis using boundary element methods is a growing area of investigation. In particular, BIE formulations for shape sensitivity of field variables or objective functions using analytical domain differentiation are frequently invoked in connection with shape optimization or inverse problems. In general, there is no privileged relationship between sensitivity analysis and Galerkin symmetric BIE formulations. Of course, SGBIE allow to formulate governing equations of adjoint problems or of field variables sensitivities, and thus provide one possible solution tool for the intermediary computations needed in most sensitivity analyses involving boundary value problems. However, one area of sensitivity analysis in which the SGBIE do play a privileged role is shape sensitivity problems in elasticity in the presence of cracks, as shown by the following considerations.

In shape sensitivity analysis, the adjoint variable approach is the most efficient one to objective function sensitivity evaluation, especially when the governing equations are linear. Usually (i.e. in the absence of geometrical singularities like cracks), it is easy to formulate such objective function sensitivities in terms of boundary integrals only, even when domain integral functionals are considered. However, the relevant derivations break down in the presence of cracks due to the appearance at some point of the calculation of nonintegrable crack front singularities (essentially, what happens is that Eulerian derivatives of the strain energy density are non-integrable at the crack front); see e.g. Ref. [START_REF] Bonnet | Sensitivity approach for internal crack or void perturbation using boundary integral equation and adjoint approach[END_REF]. Thus the combined application of BEM and classical adjoint variable techniques is problematic. In this case, the SGBIE framework allows to circumvent those difficulties.

For instance, let us consider again the cracked elastic solid defined in Sec. 4 with Neumann boundary conditions alone (i.e. given traction p on Γ, traction-free crack faces), together with an objective function J:

J = Γ f (u) dS (128)
For example, 2f = -p.u gives the potential energy and 2f = |u -umeas| 2 gives the L 2 -distance between computed (for an assumed crack S) and measured displacements, frequently used to solve crack identification problems. One has to assess the sensitivity of J to domain perturbations that affect the crack S but leave the external boundary Γ fixed. To this purpose, infinitesimal perturbations

using the notation: f,m ≡ ∂f /∂ym. The material derivative of a generic regular integrals over a moving surface S p is given by:

d dp S p f dS = S p ⋆ f +f Diθi dS ( 130
)
where Di = ∂i -ni∂n denotes the tangential part of a partial derivative (Diθi is the surface divergence divSθ of θ).

Domain formulation

Define the adjoint displacement field u a as the solution to the elastic problem with Neumann boundary data

pa = - ∂f ∂u on Γ pa = 0 on S (131)
Then, apply Betti reciprocity identity to the adjoint state u a and the material derivative ⋆ u (the latter is formally a displacement field that solves a certain elasticity problem with initial strain-like body forces defined in terms of the original elastic state (u, ε, σ)). Such application of Betti's identity leads to the result dJ dp = Ω σ T ε a div θσ T (∇u a .∇θ)σ aT (∇u.∇θ) dΩ (132) It turns out that the above expression can be recast as a boundaryonly formula using a series of integration by parts and the equilibrium relations satisfied by the initial and adjoint states. However, this manipulation breaks down in the presence of cracks because cartesian derivatives of the strain energy density are encountered in the process and are not integrable in the vicinity of the crack front. Then, the use of Eq. ( 132) in a BEM context becomes impractical.

For the special case of potential energy at equilibrium, u a = -1 2 u and Eq. ( 132) defines the so-called θ-method for computing the energy release rate in linear elastic fracture mechanics, Ref 53.

Boundary formulation.

Recall that the governing SGBIE formulation reads One sees here the essential role played by the symmetry of the governing integral operator: it allows to express the final result [START_REF] Polizzotto | Variational principles for boundary element formulations in structural mechanics[END_REF] without explicit reference to the lagrangian derivatives { ⋆ X T u , ⋆ X T ϕ } of the boundary unknowns. Equation [START_REF] Polizzotto | Variational principles for boundary element formulations in structural mechanics[END_REF] thus expresses the sensitivity of any objective function of the form [START_REF] Parreira | A numerical integration scheme for the Galerkin approach in Boundary Elements[END_REF] to domain perturbations in a boundary-only form that accommodates cracks. In particular, it can be applied to the special case of potential energy, leading to a θ-integral method for a boundary-only approach to the computation of energy release rate, see Ref. [START_REF] Bonnet | Stability of crack fronts under Griffith criterion: a computational approach using integral equations and domain derivatives of potential energy[END_REF].

COUPLING OF BEM WITH FEM

The BEM and the FEM possess specific computational features that make them complementary to each other. Whereas the FEM is well suited for problems with nonlinearities and inhomogeneities in bounded domains, the BEM turns out to be especially appropriate for linear homogeneous problems with unbounded domains or in the presence of singularities (e.g. crack problems). For this reason, often in engineering practice there is computational advantage in making use of both methods. Zienkiewicz et al. first in Ref. [START_REF] Zienkiewicz | The coupling of the finite element method and boundary solution procedures[END_REF] dealt with this topic in elasticity by dividing the problem domain into two subdomains, one for FE discretization, the other for BE discretization. Their main idea was to generate a boundary element node stiffness matrix to assemble with the other finite element matrices within a FEM package. Two difficulties arose in the above treatment, one concerned with the impossibility to obtain a symmetric stiffness matrix by means of the conventional direct BEM, the other concerned with the correct way to enforce the continuity conditions for displacements and tractions through the interface between the two subdomains. The 'forced symmetrization' of the indirect BEM stiffness matrix suggested in Ref. [START_REF] Zienkiewicz | The coupling of the finite element method and boundary solution procedures[END_REF] [START_REF] Reissner | On some variational theorems in elasticity[END_REF]. As to the interface continuity conditions, the most common strategy adopted consists in using the same modelling for the displacements and/or tractions on the two interface sides. The correct way these continuity conditions are to be complied with was discussed by Hsiao in a review paper, Ref. 84, devoted to the BEM/FEM coupling methods. BEM/FEM coupling via domain decomposition, including algorithms and parallel realizations, appear among the contributions in Refs. [START_REF]Boundary elements: implementation and analysis of advanced algorithms[END_REF][START_REF]Boundary element topics[END_REF].

A more general view point for the BEM/FEM coupling problem consists in adopting a single variational statement for the two subdomains such as to obtain the two sets of governing equations as the pertinent Euler-Lagrange equations with the interface continuity conditions as natural boundary conditions. Such a variational approach leads to consistent FE and BE discretizations with independent modellings on the interface. Examples of this approach were given by Be- The ability of the symmetric Galerkin BEM to provide the exact (symmetric) stiffness matrix associated to the relevant generalized boundary displacement can be shown as follows.

Let Ω denote again an elastic body, whose boundary Γ is closed. Let (u, p) denote any compatible (i.e. which is the trace on the boundary of an elastostatic state u, σ in Ω) pair of displacements and tractions on Γ. The elastic strain energy W for any such state is given by

W = 1 2 Ω σ T (x)ε(x) dΩx = 1 2 Γ p T (x)u(x) dΓx (138)
For simplicity, consider only discretized forms of (u, p) according to Eq. [START_REF] Belytschko | A variationally coupled finite element-boundary-element method[END_REF]. In particular, W is then given by

W = 1 2 P T Γ Ψ T p (x)Ψu(x) dΓx U (139)
for any compatible pair (U , P ). First, considering the traction p as induced by given displacement u (i.e. Γu = Γ, Γp = ∅), the discretized displacement and traction are related by the following matrix equation, which comes directly from Eq. ( 23) and the regularized expressions [START_REF] Michael | Galerkin formulation and singularity subtraction for spectral solutions of boundary integral equations[END_REF], [START_REF] Monegato | Numerical evaluation of hypersingular integrals[END_REF]

: ĜuuP = (Eup + Lup)U (140)
where

Ĝuu = Γ Γ Ψ T p (x)Guu(x, ξ)Ψp(ξ) dΓx dΓ ξ Eup = Γ κ(ξ)Ψ T p (ξ)Ψu(ξ) dΓ ξ (141) Lup = Γ Γ Ψ T p (x)Gup(x, ξ)[Ψu(ξ) -Ψu(x)] dΓx dΓ ξ
Then, considering the displacement u as induced by the given traction p (i.e. Γp = Γ, Γu = ∅), the following matrix equation also holds, as a direct consequence of Eqs. ( 23), [START_REF] Maier | A Galerkin approach to boundary element elastoplastic analysis[END_REF], [START_REF] Nappi | An internal variable approach applied to boundary element elastic-plastic analysis[END_REF]:

ĜppU = (Epu + Lpu)P (142)

having set Ĝpp = Γ Γ RΨ T u (x)Bpp(x, ξ)RΨu(ξ) dΓ ξ dΓx Epu = Γ (κ(x) -1)Ψ T u (x)Ψp(x) dΓx (143) Lpu = Γ [Ψu(x) -Ψu(ξ)]Gpu(x, ξ)Ψ T p (ξ) dΓx dΓ ξ
Finally, premultiplying Eq. ( 140) by P T and Eq. ( 142) by -U T and adding the resulting relations, we obtain:

P T ĜuuP -U T ĜppU = P T EupU -U T EpuP (144)
where use was made of the equality L T pu = Lup, flowing from Eqs. (141b) and (143b). At this point, one has simply to observe that

P T EupU -U T EpuP = P T (Eup -E T pu )U , Eup -E T pu = Γ [κ(x) -(κ(x) -1)]Ψ T p (x)Ψu(x) dΓx = Γ Ψ T p (x)Ψu(x) dΓx
Therefore Eq. ( 144), in view of Eq. ( 139), becomes

2W = P T ĜuuP -U T ĜppU (145)
in which vectors U and P must be compatible in the above sense. Now, solving Eq. ( 140) for P and substituting the result into the above equation ( 145) leads to:

W = 1 2 U T KU (146)
where:

K = (Eup + Lup) T Ĝ-1 uu (Eup + Lup) -Ĝpp (147)
On the other hand, solving Eq. ( 142) leads to the complementary strain energy interpretation of W , i.e.:

W = 1 2 P T KcP (148)
having set:

Kc = Ĝuu -(Epu + Lpu) T Ĝ-1 pp (Epu + Lpu)
In the foregoing analysis, the domain Ω may represent a subdomain of some region of space. The above results, Eqs. ( 146), [START_REF] Reissner | On some variational theorems in elasticity[END_REF] then provide the subdomain stiffness matrix, which can e.g. be added to the stiffness matrix for the complementary subdomain, obtained by a FE approach. Note that K and Kc are symmetric. Besides, either Ω or its complement with respect to Ω∞ is bounded, so that the above results allow to take care of an unbounded, linearly elastic, portion while a FE modelling is used for the bounded complement, e.g. because of nonlinear behavior. Note, however, that the actual computation of K or Kc needs the assembly of several integral operators over Γ and the inversion of either Ĝuu or Ĝpp (the latter being invertible only after provisions for removing rigid-body solutions have been made). Besides, equation ( 145) is only approximately true in a discretized setting, partly because [START_REF] Polizzotto | A consistent boundary/interior element method for evolutive elastic-plastic structural analysis[END_REF] and [START_REF] Polizzotto | A variational formulation of the BEM for elastic-plastic analysis[END_REF] cease to hold simultaneously in the discrete case (i.e. with the same U , P in both). However, the continuous counterpart of ( 145) is true. These considerations are analyzed in detail in Ref. [START_REF] Sirtori | Symmetric coupling of finite elements and boundary elements[END_REF].

Incidentally, the latter point gives an a posteriori justification of the use of discretized variables in the analysis. An explicit continuous expression for W in terms of either u or p on Γ is indeed not available except when special choices are made for the fundamental solutions. For example, if Guu is the Green's function of the domain Ω, i.e. if Guu(x, ξ) = 0 for any ξ ∈ Γ, then K = -Ĝpp.

Energy-based formulations of the BEM/FEM coupling problem in the framework of elastostatics were given in Ref. 144 through coupled variational principles that combine the boundary min-max principle of Refs. 133, 138 with classical principles of elasticity. These coupled principles take on the form of a saddle-point principle on making use of the minimum total potential energy or the the minimum complementary energy. In any case, the continuity conditions at the interface of the FE and BE subdomains are satisfied as natural conditions, which means that, in the numerical treatment of the problem, the trial functions need not satisfy the continuity conditions.

Coupling across interfaces between individually homogeneous SGBE-discretized subdomains of an inhomogeneous solid was formulated and implemented in Ref. 69 for two-dimensional static problems, and discussed on a variational basis in Ref. 100 for elastodynamics.

COMPUTER IMPLEMENTATIONS AND COMPUTATIONAL EXPERIENCES

It has been shown in Sec. 8 that SGBIE formulations are amenable to regularized form, in the sense that the various double integrations can be given the form of a weakly singular integral followed by a nonsingular integral. This provides a convenient basis for implementation. In fact, in the regularized SGBIE formulation, the kinematic (resp. static) boundary variables and trial functions can be modelled using a continuous (resp. piecewise continuous) interpolation. Thus, all the usual BEM discretization techniques are applicable. Let us then discuss the implementation of the SGBEM using standard BEM modelling. The boundary Γ is divided into elements Γe, each being mapped onto a reference (or 'master') element Me (either the square -1 ≤ η1, η2 ≤ 1 or the triangle η1, η2 ≥ 0, η1 + η2 ≤ 1 for three-dimensional problems). Nodes and shape functions are introduced in the usual fashion, so that the location of a point x on Γ is expressed in terms of vector shape functions Ψ k Γ (η) and nodal coordinates a k in the physical space of the problem, k running over the nodes of element Γe. The modelling of the geometry and the unknowns can thus be symbolically expressed, on a given element Γe, in terms of the intrinsic coordinates η as

x(η) = ΨΓ(η)A u(x) = Ψu(η)Xu (149) p(x) = Ψp(η)Xp
Note that the possibility of discontinuous unknown and trial tractions must be allowed for in the definition of Ψp, e.g. using nonconformal interpolation or multiple traction nodes where necessary. Then, upon substituting the above formulas, double integrals over element products Γa × Γ b must be computed and assembled, and three situations arise as follows (χ and η will denote the intrinsic coordinates on the master elements Ma and M b , respectively). (i) Disjoint elements: Γa ∩ Γ b = ∅. The two elements do not share a common side or vertex, so that the integral over either element is nonsingular. The numerical evaluation of the double element integral is straightforward using e.g. a Gaussian product rule. Like in more traditional BEMs, the number of Gauss points should be adjusted according to the relative proximity of the two elements (the number of Gauss points required increases as Γa

Algebraic manipulations allow to show that r can be expressed as

r 2 = b 2 r2 (a, b) (152)
with r(a, b) = 0 at the singularity χ = η where b = 0 as a consequence. The integral Iaa finally takes the form ( f denoting the function resulting from f after the change of variables):

Iaa ≡I r aa + I s aa = 1 0 db 1-b b-1 f (a, b) ln r2 (a, b) + f (a, -b) ln r2 (a, -b) da + 2 1 0 db 1-b b-1 [ f (a, b) + f (a, -b)] ln b da (153)
In Eq. ( 153), the inner integral in the first addend I s aa , as well as both integrals in the second one I r aa , can be evaluated using standard quadrature rules (e.g Gaussian), while the outer integral in I s aa should be computed using a logarithmic Gaussian rule, Ref. [START_REF] Stroud | Gaussian quadrature formulas[END_REF].

For another type of potentially singular integral, of the form

Jaa = +1 -1 +1 -1 f (η, χ)[Ψ(χ) -Ψ(η)] 1 r dη dχ
one observes, similarly to [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF], that:

Ψ(χ) -Ψ(η) = b Ψ(a, b)
with Ψ(a, b) = 0 at the singularity χ = η where, hence, b = 0. Hence Jaa is in fact nonsingular and can be computed using standard quadrature rules.

In three-dimensional situations, one has to consider singular integrals of the form ( dη and dχ denoting dη1 dη2 and dχ1 dχ2 respectively):

Iaa = Ma Ma f (η, χ) 1 r dη dχ (154)
which can be treated as a weakly singular integral with respect to χ followed by a regular integral with respect to η. In order to do so, for any fixed η, the parent element Ma is divided into up to four triangles, each defined by one element edge and the common vertex η, and the 1/r singularity is taken care of using an adequate change of variables. For instance, for quadrilateral elements, the triangle containing the edge χ1 = +1 is transformed into new intrinsic coordinates -1 ≤ a, b ≤ 1 according to:

χ1 = η1 + 1 2 (1 -η1)(a + 1) χ2 = η2 + 1 2 (b -1)(a + 1) (155)
Then, some algebra on the interpolation functions yields: 1 -η1 2 da db are nonsingular. Therefore, the corresponding χ-integral can be evaluated using Gaussian quadrature, in the (a, b) coordinates. Transformations similar to (155) must be performed separately for the other triangles. Note that it is possible to resort to a polar coordinate transformation on Ma, instead of (155); the resulting integral is again nonsingular due to similar considerations as above.

The subsequent, nonsingular, integration with respect to η is performed using Gaussian quadrature. This means, in particular, that all sample η points actually used e.g. to define [START_REF] Sirtori | Symmetric coupling of finite elements and boundary elements[END_REF], lie in the interior of Ma. (iii) Adjacent elements: Γa = Γ b , Γa ∩ Γ b = ∅. Unlike (i) and (ii), this case is specific to SGBEMs inasmuch no similar calculation arise in traditional collocation BEMs. This is, like (ii), a case of singular integration, but the singularity, weaker, occurs when x, ξ are located on the edge shared by the two elements. In two-dimensional situations, a typical element integral is

I ab = Γa Γ b g(x, ξ) ln r dΓx dΓ ξ = +1 -1 +1 -1 f (η, χ) ln r dη dχ ( 156
)
where the singularity r = 0 now occurs for η = 1, χ = -1. Introduce the transformation (151) if χη ≥ 0 (no singularity) and, otherwise, the new variables γ, δ such that:

η = γ(δ -1) + 1 , χ = γ(δ + 1) -1 , dη dχ = γ dδ dγ with -1 ≤ δ ≤ 1, 0 ≤ γ ≤ 1.
As a result, one obtains:

r = γ r(δ, γ) ; I ab =2 1 0 db 1-b b-1 f (a, b) ln r da + 2 1 0 γ ln γ dγ 1 -1 f (δ, γ) dδ (157) + 1 0 γ dγ 1 -1 ln r2 (δ, γ) f (δ, γ) dδ
where all integrations are performed using standard means, except for the γ-integration in the second integral, for which logarithmic Gaussian quadrature should be used instead.

A similar approach can be developed for adjacent integrations in three-dimensional problems, see in particular Ref. 5 where all cases of double singular integrations are addressed for the SGBIE formulation of the Neumann problem and in a rather general discretization setting.

Whatever the particular discretisation approach used, a system of linear equations of the form (23), i.e. It is important to note that the matrix is symmetric but not signdefinite (except for the special cases of pure Dirichlet and pure Neumann problems). Thus, Cholesky decomposition does not apply, except for the special cases. However, direct linear solvers for indefinite symmetric matrices are expounded in Ref. [START_REF] Golub | Matrix computations (second edition)[END_REF] and codes are provided e.g. in software libraries like LAPACK, Ref. 4. Such direct solvers, as well as the Choleski algorithm for sign-definite symmetric matrices, need about n 3 /6 arithmetic operations for n × n matrices; in contrast, direct solvers for unsymmetric n × n matrices need about n 3 /3 arithmetic operations.

Due to its particular structure, the coefficient matrix also lends itself to iterative solution techniques like biconjugate gradient or generalized minimal residual (GMRES). In such techniques, a natural preconditioner is the block-diagonal system of equations obtained by ignoring the non-diagonal blocks Ĝup up , Ĝup pu in Eq. ( 158), which makes advantageous use of the Choleski decompositions of the diagonal blocks Ĝuu uu and -Ĝpp pp . Iterative solvers are discussed in Refs. 71 (in general terms) and 75, 165 (in connection to SGBEM).

The techniques outlined above for the computation of double element integrals allow, in principle, to use any of the classical BEM discretization schemes, and accommodates curved elements and high-degree interpolations. In practice, general-purposes implementations of SGBEM methods based on purely numerical schemes are not commonly found in the literature. SGBEM implementations that follow the 'purely numerical' approach are presented in Refs. 60, 129 for two-dimensional elasticity and in Refs. 59 for bending of Kirchhoff elastic plates. Symbolic manipulations have been used in Refs. 9, 10 to implement two-dimensional SGBEM formulations with curved elements, for both thermal and elastic problems. An implementation for fracture analysis in plane orthotropic elasticity is reported in Ref. [START_REF] Gray | Symmetric Galerkin boundary integral fracture analysis for plane orthotropic elasticity[END_REF]. SGBEM implementations for 3D problems comparable to general-purpose 3D CBEMs are scarce. Ref. 95 presents a direct SGBEM implementation for three-dimensional linear fracture analysis of bounded elastic bodies, based on a formulation very close to Eq. ( 48) and using eight-noded boundary elements and numerical integration throughout; it is to our knowledge the best attempt in this direction.

In addition to general numerical integration techniques, many specific techniques have been developed and implemented. In particular, exact analytical element integrations for the case of straight line elements with piecewise linear or constant interpolation have been performed for two-dimensional problems, either using a complex variable formalism, Ref. 154, or directly in real variables, Ref 38. Also, exact expressions for some of the element integrals (mostly the singular ones, which are critical) have been employed e.g. in Ref. [START_REF] Ben Mariem | A new boundary finite element method for fluid-structure interaction problems[END_REF]. The first pand hp-versions of the SGBEM have been proposed, implemented and tested in Refs. [START_REF] Holzer | A p-version of the symmetric BEM[END_REF][START_REF] Holzer | The h-, pand hpversion of the BEM in elasticity: numerical results[END_REF] Perhaps the main obstacles to attempts at general implementations of SGBE methods, especially for three dimensions, are the complexity of the formulation (compared to traditional collocation BE methods) and the fact that a general regularization approach was not known until a few years ago. The general regularization issue is now well understood in principle, but many potential applications are yet to be made.

In addition to other desirable features provided by symmetry, as discussed in various passages of the present survey, the main computational gain provided by SGBE methods over traditional BEMs is related to the symmetry of the final system of equations, whereby a factor of two is gained in both solution computer time and computer storage. An objection frequently raised against SGBE methods is that the double integrals are expected to incur significantly higher demands in numerical quadrature and thereby slow down the assembly stage. However, this is not necessarily the case. Indeed, consider as an example a very fine mesh of N (say) quadrangular eight-noded curved elements, with the usual polynomial interpolation functions. From Euler formula, this mesh contains 3N + 2 nodes. Neglecting coincident, adjacent or otherwise singular integrations (whose relative proportion of the overall integration burden decreases as the mesh is refined) and assuming that a uniform number nG of Gauss points is used for each single element integral, a rough count shows that the Galerkin approach needs function evaluations at a total of N (N + 1)n 2 G /2 pairs of Gauss points for the (symmetric) matrix and N 2 n 2 G pairs of Gauss points for the right-hand side, whereas the collocation BIE needs about N (3N + 2)nG function evaluations for each of the two global matrices. This naive count gives a ratio of about nG/4 between the computational demands for Galerkin and traditional BEMs (note that by 'functions' we mean 3 × 3 matrices of functions). Given that four Gauss points are asymptotically (i.e. for pairs of very distant elements) sufficient in this context, the above count leads to a ratio of about unity. The integration task, although of course influenced by many other factors, can nevertheless be reasonably expected to have similar computational needs in both approaches. Numerical experiments presented in Kane's monograph, Ref. 87, show a factor of about two. Besides, the solution stage asymptotically dominates the assembly stage (O(N 3 ) vs. O(N 2 )), so that the factor two gained in matrix inversion becomes asymptotically decisive; this is also made apparent e.g. in Kane's numerical results and their discussion. In this respect, the availability of a proper scheme for automatic adjustment of Gauss points according to the ratio interelement distance to element typical length, such as that proposed for collocation BEMs in Ref. [START_REF] Rezayat | On time-harmonic elastic wave analysis by the boundary element method for moderate to high frequencies[END_REF], is crucial, the more so because in SGBEM integration times grow with the square of the typical number of Gauss points on an element.

Many of the applications of symmetric BE methods to date address rather specific situations. In particular, research effort has been devoted to problems in unbounded media. Scattering of acoustic or elastic waves by cracks is addressed in Refs. 3 (three-dimensional problems and frequency domain) and 16 (two-dimensional problems and time domain). Computations for wave radiation in exterior domains is addressed e.g. in Ref. [START_REF] Hamdi | Sound radiation from ducts: theory and experiment[END_REF]

SYMMETRIZATIONS BY ALTERNATIVE APPROACHES

The inability of the conventional BEMs to produce symmetric equation systems induced many researchers to carry out a 'forced symmetrization' of the unsymmetric equation sets, Refs. As noticed in Sec. 7, truly symmetric BEM formulations referred to as hybrid BEM (HBEM) formulations in the literature were obtained through energy-based approaches, all of which can basically be related to Reissner's boundary stationarity principle of Ref. 148. These HBEM formulations exhibit a number of features (e.g. simple integration over the boundary, symmetric positive-definite stiffness matrix for the boundary nodes), which make them an attractive and effective alternative to the SGBEM. HBEMs were presented in Refs. [START_REF] Belytschko | A variationally coupled FE-BE method for transient problem[END_REF][START_REF] De Figueiredo | A new hybrid displacement variational formulation of BEM for elastostatics[END_REF][START_REF] Dostal | Symmetric FE-BE coupling with iteration on the interface[END_REF][START_REF] Dumont | The hybrid boundary element method in elastostatics: Overview of the theory and examples[END_REF][START_REF] Leung | 2D elastostatic analysis by a symmetric BEM/FEM scheme[END_REF][START_REF] Schnack | A hybrid BEM model[END_REF][START_REF] Schnack | Modern coupling of FEM and BEM with mixed variational formulations[END_REF][START_REF] Zeng | Stable symmetric finite elementboundary integral coupling methods for fluid-structure interface problems[END_REF][START_REF] Zeng | Unified symmetric finite element and boundary integral variational coupling methods for linear fluid-structure interation[END_REF], with different motivations, starting points and purposes. In order to give a brief account of the essential features of HBEMs, the formulation of Ref. 52 is considered here as a representative one.

Reissner's boundary stationarity principle for elasticity, in an extended form used for hybrid finite element formulations, Ref. 161, considers the functional:

H = 1 2 Γ p T u dΓ - Γ (u -ũ) T p dΓ - Γp pT ũ dΓ - 1 2 Ω u T div σ dΩ - Ω bT u dΩ (159)
where the independent unknown fields are: σ and u in Ω, p and ũ on Γ, whereas p is the traction associated to σ on Γ, i.e. p = σ.n; The first addend on the r.h. side of eq. ( 161) is a quadratic form which, being representative of the strain energy stored in Ω as a consequence of the sources F , is positive for any F = 0, so that C is symmetric and positive definite, the product CF provides the generalized displacement vector dual of the source vector F . The Kuhn-Tucker conditions that characterize the stationarity of H of Eq. ( 161 where the matrix

K = M L -1 C(L -1 ) T M T
is the relevant symmetric positive-definite stiffness matrix.

The above formulation, extended to elastodynamics in Refs. 51, 62, has been applied to solve numerically sample problems with quite satisfactory results, Refs. [START_REF] Carini | Variational formulation of the boundary element method in transient heat conduction[END_REF][START_REF] Gaul | Hybrid symmetric boundary element method in elastodynamics[END_REF][START_REF] Panzeca | Boundary/field variational principles for the elastic-plastic rate problem[END_REF][START_REF] Panzeca | Dynamic plasticity analysis by boundary-interior element[END_REF][START_REF] Polizzotto | A step-wise variational approach to elastic-plastic analysis by boundary/interior elements[END_REF][START_REF] Polizzotto | Variational formulations for coupled BE/FE methods in elastostatics[END_REF]. Though an effective and computationally convenient alternative to SGBEM in elastic problems, HBEM reveals deficiencies in inelastic analysis and, more in general, as a tool for theoretical investigations. One of the reasons for this circumstance is that the HBEM makes use of the fictitious sources F which in general do not represent any real quantity of the related problem. In fact, the SGBEM turns out to be preferable for a wide class of inelastic problems, primarily as a consequence of the way the functional H, Eq. ( 159), is approximated within the HBEM.

It can be shown that, by interpreting H in a mechanically consistent way, the SGBEM can be generated from HBEM. To this purpose, let us set p = p on Γ, ũ = u on Γp, and ũ = ū on Γu; moreover, let equilibrium equations (160a) be satisfied in Ω. Thus, eq. ( 159) becomes where: -ḡ h , (h = u, p, σ), denote effects in Ω∞ due to the given external actions, Eqs. ( 16), ( 17) and (30); J h (h = u, p, σ), effects due to the unknown sources p ⋆ , -u ⋆ , i.e. Substituting Eq. ( 163) into Eq. ( 162), after some reordering, gives:

-

where C denotes a collection of integrals dependent on the external actions only. Noting that the expression within the square brackets vanishes by Betti's theorem (applied to Ω∞ subjected to the known and unknown sources, respectively), remembering Eqs. ( 81) and [START_REF] Becache | A variational boundary integral equation method for an elastodynamic antiplane crack[END_REF][START_REF] Becache | Regularization in 3D for anisotropic elastodynamic crack and obstacle problems[END_REF], one finally obtains:

-H = Π[p * , u * ] + C.
This result means that, to within the inessential constant C, the functional -H coincides with that pertaining to the boundary min-max principle of Sec. 7, which is known to constitute a specific energy basis for the SGBEM.

To conclude this section, other contributions to the issue of symmetric BEM formulations are briefly reviewed. In Ref. [START_REF] Kanarachos | On the symmetrization of the BEM formulation[END_REF], the authors simply start from the usual (domain) weak formulation of a linear (e.g. potential) problem. Both the physical unknown and its associated test function are modelled using shape functions that satisfy the relevant (homogeneous) partial differential equation (and thus are constructed using integral representation formulas), allowing to convert the energy bilinear form into a boundary integral. To avoid problems caused with singularities, the shape functions are defined using sources distributed on a curve external to the boundary. Also, setting up the 'stiffness' matrix then requires three-fold integrations. In Refs. 159 and 160, the minimum total and complementary energy principles of elasticity are employed and the field and boundary variables are suitably represented, through Kelvin kernels, in terms of discrete sources located outside the boundary. Depending on the kinds of unknown fields discretized and of sources employed, BEM formulations are obtained in the form of dual quadratic programming problems. These formulations are in some way similar to the HBEM formulations since they all use indirect representations of the unknown fields. In Ref. 34, a symmetric system of boundary integral equations is obtained, for linear elasticity, using representations of the kinematic and static field variables in terms of Helmholtz and Beltrami potentials, respectively (called 'conjugate variables' in that paper), and formulating the conditions under which the quadratic error in (elastic) constitutive equation, integrated over the domain occupied by the solid, is stationary.

FUTURE PROSPECTS AND CONCLUSIONS

Compared to the past and continuing developments of the traditional BEM in almost all areas of applied mathematics and engineering sciences, the SGBEM, as a relatively young branch of the BIE-based methodology, at the present time turns out to be little developed, especially in terms of numerical applications. The main reason for this is, in the authors' opinion, the conceptual, mathematical and computational difficulties implied by the hypersingular integrals featuring the method. Though several strategies have been, and are being, formulated and implemented in particular application contexts (see Secs 8 and 12), no general effective interpretations and consequent computational methods have been proposed yet for evaluating these hypersingular integrals with the desirable confidence. In the majority of application areas, the conventional BEM may still appear to be a computationally more attractive, easier to implement, user-friendlier alternative than the SGBEM.

However, as Francis Bacon wrote, "in more difficult things we should not expect that someone sows and harvests soon, but hope that they mature gradually" (verbatim: in rebus difficilioribus non est expectandum ut quis simul serat et metat, sed opus est ut per gradus maturescant). In the maturing process of the SGBEM, among the next stages at least the following advances are reasonably expected: a general method, mathematically justified in a way accessible to engineers, for the treatment of hypersingular integrands, as satisfactory as the Cauchy principal-value theory and possibly emerging as optimal in most cases from the diverse trends of nowaday research; criteria of algorithmic stability in time-stepping solutions of evolutive analyses, especially of those concerning nonlinear time-dependent problems, more particularly in the presence of constitutive material instabilities; coupling between SGBEM and FEM exploiting common frameworks of variational theorems; computer implementations apt to enrich the numerical experience (in particular in the areas where none is available so far) and to validate the reasonably wellfounded promises of SGBEM in comparative terms of accuracy and cost-effectiveness.

In Warner Koiter's words, "research in engineering sciences cannot be considered to be completed until its full numerical evaluation has been obtained, and the results have been presented in a form accessible to other engineers". In this direction, research on the SGBEM has undoubtedly still a long way to go, but, in the writers' opinion, it is likely to attract interest and contributions from more and more researchers in engineering mechanics and applied mathematics. Its motivation will not be merely the aesthetical search for symmetry or the dissatisfaction for its loss, as mentioned in the Introduction, but primarily by the ability of the SGBEM of preserving essential features of the original (continuum, differential) formulations of physical problems, its versatility and its special fitness to various problems of engineering interest, such as fracture and inelastic analysis and sensitivity analyses for structural optimization and identification.

  . Several contributions, e.g. those by Gu and Hew, Ref. 72, and Altiero and Gioda, Ref. 2, are within this line of thought. (C) Finally, if Ψ ⋆ i = Ψi because the interpolation functions are Dirac functions and only static sources are adopted (i.e. d ⋆ p = 0, d ⋆ u = 0), the starting point of the traditional BEM is recovered from Eq. (21). The following comments may supplement what precedes and provide links with what follows.

  been derived (and presented in Appendices) in the first papers on SGBEM on viscoelasticity, Ref.[START_REF] Carini | A variational formulation of viscoelastic analysis by boundary integral equations[END_REF] and poroelasticity, Ref. 123.

  various submatrices are given by Ĝpp

  occurs in that all double-integral terms on the right-hand side are cancelled out, because the kernels Gup, Gpu depend linearly on the exterior unit normal to the relevant domain, and Gpp depends quadratically on it.Indirect SGBE formulations are natural symmetrizations of indirect BEM obtained from e.g. single-layer representations for Dirichlet problems and double-layer representations for Neumann problems. Such, formulations have been applied to radiation and scattering problems. For instance, an early indirect SGBIE formulation for mixed boundary value problems in exterior acoustics was established in Ref.[START_REF] Hamdi | Une formulation variationnelle par équations intégrales pour la résolution de l'équation de Helmholtz avec des conditions aux limites mixtes[END_REF]; various symmetric indirect formulations for Dirichlet and Neumann problems are presented in Ref. 121.

(n = 2

 2 or 3: dimension of the geometrical space) where [s] is the integer part of s, α is a n-tuple of positive integers, |α| = |α1 + . . . |αn|| and D α f = (f ,α 1 ,...,f,α n ); the space H -s (Γ) is the dual space of H s (Γ) with respect to the L 2 scalar product (f , g) L 2 (Γ) = Γ f T g dΓ. The regularity of the functions of H s (Γ) increases with the value of s. It is interesting to note that for halfinteger values of s a kernel function similar to the fundamental solutions used in BEMs appears in the definition of the norm |f | H s (Γ) .

  of the form y p → y p+δp = y p + θ(y)δp + o(δp) with lim δp→0 |o(δp)| δp = 0 of a domain Ω p are introduced, where δp is the small variation of some shape parameter p upon which Ω is dependent, and θ is the domain transformation velocity. The material derivative ⋆ f (y) (also termed total or lagrangian derivative) of a generic field variable f (y p , p) in the domain transformation is defined as: ⋆ f (y, p) = lim δp→0 1 δp [f (y p+δp , p + δp)f (y p , p)] with y = y p = ∂f ∂p + f,mθm (y, p)

⋆

  the notations of Secs. 2,4 are used. Moreover, since the external boundary Γ is fixed, one has dJp dp = in general, this problem is not equivalent to one mentioned above with Neumann data[START_REF] Polizzotto | A unified treatment of shakedown theory and related bounding techniques[END_REF]. Besides, application of the lagrangian derivative operation to the SGBE formulation (133) leads to a governing linear system of equation for { the detailed expressions of the various d/dp come from application of formula (130) and the fact that neither Γ nor p are modified under the crack perturbation has been taken into account. Now, left-multiply Eq. (136) by { then subtract the resulting equations from one another. Taking advantage of the symmetric character of the matrix governing operator, one obtains the result

1 f

 1 , Γ b get closer to each other). (ii) Coincident elements: Γa = Γ b . In two-dimensional situations, a typical element integral is Iaa = Γa Γa g(x, ξ) ln r dΓx dΓ ξ (η, χ) ln r dη dχ (150) where r = |x -ξ| and the regular function f (η, χ) gathers interpolation functions, nonsingular parts of kernel functions and element jacobians. Then, following e.g. Refs. 129 or 60, one can introduce new intrinsic variables a, b: η = ab , χ = a + b dη dχ = 2 da db

2

 2 1)r(a, b) Ψu(χ) -Ψu(η) = 1 2 (a + 1) Ψ(a, b) where once again r(a, b) = 0 at the singularity χ = η, i.e. a = -1 there. Ψu(χ) -Ψu(η) dη = 1 r2 Ψ(a, b)

  , at least for boundary-value problems with mixed Dirichlet and Neumann data.

= 1 2 F

 2 σ and u are related to each other by Hooke's law, i.e. σ = Eε(u), where ε(u) = ∇ s u (∇ s : symmetric part of the gradient operator ∇); ũ is required to satisfy the essential boundary condition ũ = ū on Γu. It can be easily proved, Ref.[START_REF] De Figueiredo | A new hybrid displacement variational formulation of BEM for elastostatics[END_REF], that the Euler-Lagrange conditions that characterize the stationarity of H with respect to the above four fields read:div σ + b = 0 in Ω u = ũ, p = σ.n = p on Γ, p = p on Γp(160)With the HBEM, the functional H is discretized by representing u and p in terms of a discrete set of point sources and Kelvin's fundamental solutions, as well as by interpolating ũ and p from nodal values through suitable shape functions:u(x) = Gu(x)F , p(x) = Gp(x)F ũ(x) = Ψu(x)U , p(x) = Ψp(x)PMatrices Gu and Gp are constructed with the kernels Guu and Gup, respectively, and vector F collects point sources (i.e. concentrated forces) with location points over the boundary (or possibly outside of it). Assuming σ = GσF + σ b , where Gσ is constructed with the kernels Gσσ and σ b denotes some stress field in equilibrium with bin Ω, the functional H becomes:H T CF -F T LP + U T M P -U T lp -F T lΩ(161) having set C = Γ G T u Gp dΓ L = Γ G T u Ψp dΓ, M = Γ Ψ T u Ψp dΓ lp = Γp Ψ T p dΓ lΩ = 1 2 Ω G T u b dΩ

  ) read:CF = LP + lΩ, L T F = M T U M P = lpwhich are the algebraic equivalent of Eqs.[START_REF] Texeira De Freitas | A kinematic model for plastic limit analysis of solids by the boundary integral methods[END_REF], respectively. Slight manipulations of the above equations give Ku = l ≡ lp + M L -1 lΩ,

-H = 1 2 2 2 Ω

 222 Γu p T u dΓ -1 Γp u T p dΓ -Γu p T ū dΓ + Γp u T p dΓ + 1 bT u dΩ (162)where p = σ.n, σ = E∇ s u and u is a solution to the Navier equation of elasticity theory. Let the unknown fields of eq. (162) be represented by means of Kelvin's solutions for suitable sources applied in Ω∞, in which Ω is thought to be embedded. These sources consist of the following quantities: body forces b in Ω, simple layer (forces) p in Γp and double layer (relative displacements) -ū in Γu; simple layer (unknown forces) p ⋆ in Γu and double layer (unknown relative displacements) -u ⋆ in Γp. These produce the following effects:u = J u(p ⋆ , u ⋆ )ḡu , σ = J σ (p ⋆ , u ⋆ )ḡσ (in Ω) u = J u(p ⋆ , u ⋆ )ḡu , p = p ⋆ḡp (on Γu)(163)p = J p(p ⋆ , u ⋆ )ḡp u = u ⋆ḡu on Γp

J

  h (p ⋆ , u ⋆ ) = Γu G hu p ⋆ dΓ -Γp G hp u ⋆ dΓ, (h = u, p, σ).

  , Ref. 153; for some structures (beams and Kirchhoff plates) in 1985 by Hartmann et al., Ref. 79, and for elastoplastic solids in 1987 by Maier, Polizzotto, Ref. 112. Earlier symmetric formulations include one by Nedelec, Ref. 120 for scalar potential problems and another in a 1977 paper by Bui, Ref. 32, devoted to BIE for plane cracks under mode I loading.

  -Brebbia's 'dual reciprocity technique' (Ref. 118). While this technique has still to be tested in the context of concern herein, the SGBEM based on static kernels for inelastic dynamics has been enriched (Ref. 101) by theoretical results with computational interest, such as those established for SGBEM in quasi-static plasticity (Ref. 42) and mentioned in Sec. 3.3.

  a boundary stationarity principle for linear elasticity and a related Ritz type solution procedure. Zienkiewicz et al. employed in Ref. 176 the same stationarity principle to obtain a symmetric BEM formulation on the basis of a standard BEM procedure. Substantially equivalent formulations were proposed and numerically implemented by Schnack (Refs. 150, 151), Dumont (Refs. 56, 57), De Figuereido and Brebbia (Ref. 52), Davì (Ref. 50) and other authors (see Sec. 13).

  and are the equations governing the modelled boundary fields. The boundary min-max principle above has been extended as to include frictionless unilateral contact boundary conditions in Ref.[START_REF] Polizzotto | Variational boundary-integral-equation approach to unilateral contact problems in elesticity[END_REF] and deformation-theory plasticity in Refs. 140, 142, 143. Another, related, energy interpretation of the SGBEM as described in Sec. 2 is available. It is based on the fact, established in Sec. 11, that the strain energy W in Ω associated with any elastostatic state is given by

Ĝpp

pp P -Ĝpu pu U = Rp, Ĝup up P -Ĝuu uu U = Ru (85)

  -127, deserve further study. An alternative variational BEM formulation was given by Maier et al. in Ref. 99, 100 with the use of a boundary integral operator which is grounded on works on convolutive variational principles by Gurtin, Ref. 74, and Tonti, Ref. 162. Let us set:

  has been criticized by various researchers (see e.g. Ref. 113). Symmetric hybrid element coupling formulations (see Sec. 13) were given by Belytschko et al., Ref. 18, 19 and by others, Refs. 55, 93, 151, based on energy methods substantially equivalent to the Reissner boundary stationarity principle, Ref.

  lytschko et al. in Ref. 19, 20 and Zeng et al. in Refs. 173-175. The symmetric Galerkin BEM was employed to this purpose by Costabel in Refs. 45, Costabel and Stephan in Ref. 46, Hsiao in Ref. 84. Wendland in Ref. 166 provided error estimates of the discretization procedure and proved the related convergence features.

  for acoustics and Ref. 21 for electromagnetism. A SGBEM formulation for scattering by thin scatterers or scatterers with corners is investigated and numerically tested in Ref. 172. Other papers deal with fracture mechanics problems: in Ref. 72, a SGBE formulation for plane cracks in three-dimensional unbounded media is used to compute stress intensity factors; Yin and Ehrlacher, Ref. 171, present a complex-variable approach to SGBEM for fracture mechanics in plane elasticity; the domain derivative approach to the computation of energy release rate outlined in Sec. 10 is treated in Ref. 26; Xu et al. in Ref 169 present SGBEM-based numerical simulations of crack growth in fiber-reinforced composites. Most of the above-mentioned situations present relatively simple geometries; besides, they appear usually as either Dirichlet or Neumann problems (i.e. do not involve mixed-type boundary conditions).

  89, 90, 98, 176. Different procedures were proposed to different specific purposes, mainly for BEM applications to problems in which symmetry represents a feature of paramount importance, like BEM/FEM coupling, Refs.[START_REF] Kolino | Hybrid stress analysis of boundary and finite elements by superelement method[END_REF] 176 (Sec. 11), and plastic deformation bounding techniques, Ref. 105 (Sec. 6.1). These symmetrization procedures, though justified by their proposers with more or less appropriate arguments and numerical experiments, have been criticized for lack of firm rational basis, Refs. 14, 19, 94, 164, and it was even recognized, Refs.[START_REF] Belytschko | Singular integration in variationally coupled FE-BE method[END_REF][START_REF] Li | A new method for the coupling of finite element and boundary element discretized subdomains of elastic bodies[END_REF][START_REF] Tullberg | A critical study of different boundary element stiffness matrices[END_REF], that forced symmetrization may lead to worse numerical results than those obtainable by unsymmetric BEMs.

  H = 1 2 Γu p * T J u(p ⋆ , u ⋆ ) dΓ -1 2 Γp u * T J p(p ⋆ , u ⋆ ) dΓ

	-	p * T (ū + ḡu ) dΓ +	u * T (p + ḡp ) dΓ
		Γu				Γp
	+	1 2 Γu	p * T ḡu dΓ -	Γp	u * T ḡp dΓ
		-	ḡT p J u(p	⋆ , u ⋆ ) dΓ +	ḡT u J p(p	⋆ , u ⋆ ) dΓ
		Γu			Γp
		+	bT J u(p	⋆ , u
		Ω

⋆ ) dΩ + C

Appl Mech Rev vol 51, pp 669-704, 1998

e.g. in Ref.[START_REF]Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity[END_REF], Eqs. (2-3) take the following limiting forms when ξ = ξ -∈ Γ -:u ⋆ (ξ) = -1 2 d ⋆ (ξ) + Γ Guu(ξ, x)f ⋆ (x) dΓx + -Γ Gup(ξ, x)d ⋆ (x) dΓx(11)p ⋆ (ξ) = 1 2 f ⋆ (ξ) + -Γ Gpu(ξ, x)f ⋆ (x) dΓx + = Γ Gpp(ξ, x)d ⋆ (x) dΓx(12)where the symbols -Γ and = Γ denote Cauchy principal value and Hadamard finite part integrals respectively. Eqs.[START_REF] Bamberger | Formulation variationnelle espace-temps pour le calcul par potentiel retardé de la diffraction d'une onde acoustique (I)[END_REF][START_REF] Bamberger | Formulation variationnelle pour le calcul de la diffraction d'une onde acoustique par une surface rigide[END_REF] can then be substituted into Betti's equation[START_REF] Aimi | Numerical integration schemes for evaluation of (hyper)singular integrals in 2D BEM[END_REF]. written with b = b ⋆ = 0 and ϑ = ϑ ⋆ = 0, we obtain:Γ f ⋆T (x) -1 2 u(ξ) +

Γu, Γp, assume a prescribed self-equilibrated loading on the crack faces:p(x ± ) = ∓p x ∈ Γ d(conventionally, the surface Γ d is oriented so that the positive unit normal n is directed from S -to S + ). Consistently with the orientation convention on Γ d , the kinematic unknown on the crack surface is the COD w = u +u -, whereas the source distribution w ⋆ on Γ d is interpreted as a fictitious COD. Then, ignoring again domain external actions for brevity like in Sec. 2, a SGBIE formulation can be obtained by invoking Betti reciprocity identityΓ p T u ⋆u T p ⋆ dΓ ξ = Γ d p T w ⋆w T p ⋆ dΓ ξ(45)with a fictitious auxiliary state (u ⋆ , p ⋆ ) defined, over Ω∞ embedding Ω, in terms of densities d ⋆ , f ⋆ on Γ and w ⋆ on Γ d :u ⋆ (ξ) = Γ [Guu(ξ, x)f ⋆ (x) + Gup(ξ, x)d ⋆ (x)] dΓx + Γ d

Γ [p T u ⋆u T p ⋆ ] dΓ ξ = 0(94)with a fictitious auxiliary state (u ⋆ , p ⋆ ) defined in terms of densities d ⋆ , f ⋆ on the auxiliary surface Γ ⋆ exterior to Γ:u ⋆ (ξ) = Γ ⋆ [Guu(ξ, x)f ⋆ (x) + Gup(ξ, x)d ⋆ (x)] dΓx (95) p ⋆ (ξ) = Γ ⋆ [Gpu(ξ, x)f ⋆ (x) + Gpp(ξ, x)d ⋆ (x)] dΓx([START_REF] Liu | Plastic analysis by quasi-high order Galerkin BEM[END_REF]

Strongly singular kernels

Consider first the well-known limit-to-the-boundary relations:

BIOGRAPHICAL SKETCHES OF THE AUTHORS