Stéphanie Chaillat
email: stephanie.chaillat@ensta-paristech.fr

Samuel P Groth
email: samuelpgroth@gmail.com

Adrien Loseille
email: adrien.loseille@inria.fr

Metric-based anisotropic mesh adaptation for 3D acoustic boundary element methods

Keywords: Acoustic wave scattering, Boundary Element Method, fast BEMs, Mesh adaptation, Anisotropic meshes

This paper details the extension of a metric-based anisotropic mesh adaptation strategy to the boundary element method for problems of 3D acoustic wave propagation. Traditional mesh adaptation strategies for boundary element methods rely on Galerkin discretizations of the boundary integral equations, and the development of appropriate error indicators. They often require the solution of further integral equations. These methods utilise the error indicators to mark elements where the error is above a specified tolerance and then refine these elements. Such an approach cannot lead to anisotropic adaptation regardless of how these elements are refined, since the orientation and shape of current elements cannot be modified.

In contrast, the method proposed here is independent of the discretization technique (e.g., collocation, Galerkin). Furthermore, it completely remeshes at each refinement step, altering the shape, size, and orientation of each element according to an optimal metric based on a numerically recovered Hessian of the boundary solution. The resulting adaptation procedure is truly anisotropic and independent of the complexity of the geometry. We show via a variety of numerical examples that it recovers optimal convergence rates for domains with geometric singularities. In particular, a faster convergence rate is recovered for scattering problems with complex geometries.

Introduction

We consider the scattering of time-harmonic acoustic waves by three-dimensional obstacles embedded within an unbounded, homogeneous domain. Various numerical methods can be used to solve such a problem, however a natural and popular approach is to reformulate the problem as a boundary integral equation. The main advantage of such a reformulation is to restrict the computational domain to the boundary of the obstacle and to exactly fulfill the outgoing radiation condition. The numerical solution of boundary integral equations is known as the Boundary Element Method (BEM) (often called the Method of Moments in the electromagnetic community). Despite the reduction in dimension of the computational domain, the main drawback of the BEM is the fully-populated nature of the system matrix. The cost of BEM simulations is thus prohibitively high when large-scale problems are concerned. Hence there has been a great deal of work since the inception of the BEM to reduce its computational cost.

If no acceleration technique is used, the storage of such a system is O(N 2 ), where N is the number of degrees of freedom on the scatterer boundary (e.g., the number of nodes in the mesh of the discretized boundary). The cost of solving the dense system using a direct method such as Gaussian elimination requires O(N 3 ) flops, whereas solution via an iterative method such as GMRES requires O(N iter N 2 ), where N iter is the number of iterations. In the last decades, different approaches have been proposed to speed up the solution of dense systems. The best-known method is probably the fast multipole method (FMM) originally proposed by Greengard and Rokhlin [START_REF] Greengard | A fast algorithm for particle simulations[END_REF] which enables a fast evaluation of the matrix-vector product required by the iterative solver. Initially developed for N-body simulations, the FMM has since been extended to oscillatory kernels [START_REF] Darve | The fast multipole method: numerical implementation[END_REF][START_REF] Greengard | Accelerating fast multipole methods for the Helmholtz equation at low frequencies[END_REF]. Now it is widely used in many application fields and has shown its capabilities in the context of mechanical engineering problems solved with the BEM [START_REF] Chaillat | Recent advances on the fast multipole accelerated boundary element method for 3D time-harmonic elastodynamics[END_REF][START_REF] Takahashi | A wideband fast multipole accelerated boundary integral equation method for time-harmonic elastodynamics in two dimensions[END_REF]. An alternative approach designed for dense systems is based on the concept of hierarchical matrices (H-matrices) [START_REF] Bebendorf | Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems[END_REF]. The principle of H-matrices is to partition the initial dense linear system, and then reduce it to a data-sparse one by finding sub-blocks in the matrix that can be accurately approximated by low-rank matrices. The efficiency of hierarchical matrices relies on the possibility to approximate, under certain conditions, the underlying kernel function by low-rank matrices. The approach has been shown to be very efficient for asymptotically smooth kernels (e.g., the Laplace kernel) and efficient in a pre-asymptotic regime for oscillatory kernels such as Helmholtz or elastodynamic kernels [START_REF] Chaillat | Theory and implementation of H-matrix based iterative and direct solvers for Helmholtz and elastodynamic oscillatory kernels[END_REF].

Mesh adaptation is an additional technique to reduce the computational cost of a numerical method. The principle is to optimize (or at least improve) the positioning of a given number of degrees of freedom on the geometry of the obstacle, in order to yield simulations with superior accuracy compared to those obtained via the use of uniform meshes. Adaptation is particularly important for scattering obstacles that contain geometric singularities, i.e., corners and ridges, which lead to a rapid variation of the surface solution near these singularities. For such problems, meshes graded toward these singularities must be employed in order to yield accurate approximations. In addition, for wave scattering problems, we may exploit the directionality of the waves in order to further reduce the number of degrees of freedom. The best strategy to achieve these goals is via so-called "anisotropic" mesh adaptation for which an extensive literature exists for volume-based methods such as the finite element method and the discontinuous Galerkin method [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF]. However, there is relatively little research attention being paid to mesh adaptation in a boundary element context. One possible explanation is the large computational cost of standard BEMs. With the development of fast BEMs such as the Fast Multipole accelerated BEM (FM-BEM) [START_REF] Coifman | The fast multipole method for the wave equation: A pedestrian prescription[END_REF] or H-matrix accelerated BEM (H-BEM) [START_REF] Börm | Introduction to hierarchical matrices with applications[END_REF], the capabilities of the BEM are greatly improved such that efficient adaptive mesh strategies are needed not only to reduce further the computational cost, but also to certify the numerical results by assessing that the theoretical convergence order is observed during the computations.

In the BEM community, the majority of the research on mesh adaptation has been confined to isotropic techniques with a focus on the Laplace equation (see, e.g., the exhaustive review [START_REF] Feischl | Adaptive Boundary Element Methods: A Posteriori Error Estimators, Adaptivity, Convergence, and Implementation[END_REF]) and extensions to the Helmholtz equation being made only fairly recently [START_REF] Bakry | Fiabilité et optimisation des calculs obtenus par des formulations intégrales en propagation d'ondes[END_REF][START_REF] Bakry | A goal-oriented a posteriori error estimate for the oscillating single layer integral equation[END_REF][START_REF] Bakry | A new accurate residual-based a posteriori error indicator for the BEM in 2D-acoustics[END_REF]. These isotropic techniques are usually based on a posteriori error analysis from which error indicators are derived. An indicator is then used to steer the mesh refinement by systematically marking and refining only elements where the error is above a specified threshold -a process known as Dörfler marking [START_REF] Carstensen | Axioms of adaptivity[END_REF]. The derivation of appropriate local error estimators [START_REF] Faermann | Local a-posteriori error indicators for the Galerkin discretization of boundary integral equations[END_REF] is a significant challenge owing to the nonlocality of boundary integral operators. This difficulty is the main reason why adaptivity for the BEM is a much less well-explored research topic in comparison to adaptivity for the FEM where the relevant operators are local differential operators. In many works convergence rates for error estimates are proven rigorously, e.g., [START_REF] Ferraz-Leite | Simple a posteriori error estimators for the h-version of the boundary element method[END_REF][START_REF] Carstensen | An a posteriori error estimate for a first-kind integral equation[END_REF][START_REF] Gantumur | Adaptive boundary element methods with convergence rates[END_REF]. However, it is seen that Dörfler refinement techniques do not usually recover the optimal convergence rates for 3D problems with anisotropic features [START_REF] Aurada | Estimator reduction and convergence of adaptive bem[END_REF]. Anisotropic variants of this strategy have been considered in [START_REF] Aurada | Estimator reduction and convergence of adaptive bem[END_REF][START_REF] Ferraz-Leite | Simple a posteriori error estimators for the h-version of the boundary element method[END_REF] however with rectangular elements for cube or cube-like shapes (where all the ridges are right angles). For these shapes they obtain the optimal convergence rate, however for general shapes (or complex geometries), their approach would not perform as well. The additional drawback of previously published works is the problem-dependent or integral equation-dependent nature of the error estimates. Also, the error analysis of these methods requires a Galerkin discretization and hence a higher computational cost than, say, a collocation discretization.

The first novelty of the present work is the extension of metric-based anisotropic mesh adaptation (AMA) to the BEM. The metric-based AMA proposed in [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF][START_REF]Continuous mesh framework part II: validations and applications[END_REF] does not employ a Dörfler marking strategy but rather generates a sequence of non-nested meshes with a specified complexity (proportional to the number of vertices or elements). The different meshes are defined according to a metric field derived from the evaluation of the linear interpolation error of the (unknown) exact solution on the current mesh. From a theoretical point of view, a continuous metric is derived from the Hessian of the exact solution. From a practical point of view, an approximate metric is derived from the numerical solution only (obtained via the BEM on a mesh). This approximate Hessian is based on the extension of typical (volumetric) derivative recovery operators [START_REF] Vallet | Numerical comparison of some hessian recovery techniques[END_REF] to the case of numerical boundary solutions. In AMA, the size, shape, and orientation of elements are adjusted simultaneously. The advantages of this approach are that it is ideally suited to solutions with anisotropic features, it is independent of the underlying PDE and discretization technique (collocation, Galerkin, etc.), and it is inexpensive. The metric-based AMA approach, as outlined above, has never been applied to the BEM. The purpose of this paper is to detail and report on the first application of metric-based AMA within a boundary element setting. Furthermore, we address some issues encountered when using an iterative solver for the FM-BEM on the resulting refined anisotropic meshes. In particular, we present two simple techniques to reduce the number of GMRES iterations required to achieve convergence when anisotropic elements are contained in the mesh.

The second novelty of this work is the combination of two acceleration techniques, namely metric-based anisotropic mesh adaptivity (AMA) and Fast Multipole acceleration. If no fast BEM is used, the capabilities of anisotropic mesh techniques cannot be fully demonstrated for realistic large scale scattering scenarios. This original combination permits us to show the performance of AMA strategy for complex real-world scattering problems such as acoustic scattering from an aircraft (see Section 5).

The outline of the paper is as follows. In Section 2, we recall the boundary integral equation formulation for wave scattering problems, the discretization of the resulting boundary integral equations (BIEs), and the acceleration of the solution of the arising linear system via the FM-BEM. In Section 3, we outline the AMA approach and extend the relevant results from [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF][START_REF]Continuous mesh framework part II: validations and applications[END_REF] for the case of a 3D volume solution to a 3D surface solution. In Section 4, we demonstrate the capabilities of the FM-BEM solver and present results illustrating the utility of a sequence of meshes to reduce the number of iterations required by the iterative BEM solver. Finally, in Section 5, an array of numerical examples is presented demonstrating the effectiveness of the AMA approach for the boundary element method. We show that optimal convergence rates are achieved with AMA for scatterers with geometric singularities, as opposed to sub-optimal rates when employing uniform meshes. Further we see how the anisotropic meshing exploits the wave directionality for higher frequency problems, thereby reducing the number of mesh elements required. For large-scale real-world scattering scenarios, we observe enormous gains overs the employment of uniform meshes.

Standard and accelerated boundary element methods

Boundary integral equations for 3D Helmholtz equation

Consider a closed bounded domain Ω ⊂ R 3 with boundary Γ and let Ω + := R 3 \Ω denote the exterior scattering domain. Suppose that we have an incident plane wave u i with wavenumber k := 2π/λ (where λ is the wavelength) and direction d; this may be written as

u i (x) = e ikd•x , x ∈ R 3 .
u i is scattered by Ω leading to a scattered field u s which, when combined with u i , gives the total field u in Ω + . This total field is such that

∇ 2 u + k 2 u = 0 in Ω + . (1)
In this work, we consider only sound-soft (Dirichlet) boundary conditions (BCs)

u = 0 on Γ. (2)
Furthermore, we require that u s = uu i satisfies the Sommerfeld radiation condition

lim r→∞ ∂u s ∂r -iku s = 0, (3)
where r = |x|. We note that the adaptation method presented later is general in that it can be applied to problems with sound-hard (Neumann), mixed, or transmission boundary conditions, in addition to the sound-soft case considered here. The sound-soft BCs, when compared to the other kinds of BC, give rise to the most severe singularities in the boundary solution when corners and ridges are present. Hence it is for this case that mesh adaptivity is most imperative and therefore why it is our focus in this paper.

The first main difficulty arising in the numerical solution to this exterior boundary value problem is related to the unboundedness of the computational domain Ω + . Integral equation based methods are one of the possible tools to overcome this issue. The approach is based on the potential theory [START_REF] Colton | Integral equation methods in scattering theory[END_REF]. For any positive real number k, let Φ k (x, y) = e ik|x-y| 4π|x -y| be the fundamental solution of the 3D Helmholtz equation. The single-layer potential operator is defined by

Sφ(x) := Γ Φ k (x, y)φ(y)ds(y), x / ∈ Γ.
The trace of the single-layer potential is given by applying the exterior Dirichlet trace to S such that we have

(Sφ) |Γ = Sφ
where the boundary integral operator S is defined by

Sφ(x) := Γ Φ k (x, y)φ(y)ds(y), x ∈ Γ.
There exist various possible integral equations to obtain the Cauchy data. For the numerical examples considered in Section 5, we use the simplest choice for the sound-soft problems. We employ the Dirichlet formulation using only the single-layer operator

Sq = -u i , (4)
which is usually called the electric field integral equation (EFIE ) since, in the electromagnetic setting, the unknown represents the electric field. The solution of the integral equation ( 4) is the trace of q = -∂u ∂n on the boundary of the domain. The value in the domain is then obtained by using the following boundary integral representation

u s (x) = Sq(x), x ∈ Ω + . ( 5
)
When the problem at hand does not have a thickness, we require a slight modification of the BIE for these problems since here the unknown is the jump in ∂u/∂n across the surface, rather than ∂u/∂n itself. We denote this jump by ∂u ∂n . So the BIE to solve is [START_REF] Stephan | Boundary integral equations for screen problems in R 3[END_REF]:

-u i (x) = Γ Φ k (x, y) ∂u ∂n (y)ds(y), x ∈ Γ. ( 6
)
It is known [START_REF] Nédélec | Acoustic and Electromagnetic Equations: Integral Representations for Harmonic Problems[END_REF] that the Integral Equation ( 4) is not uniquely solvable for all wavenumbers k. Hence, in practice, it is often preferable to employ so-called "combined" formulations, such as the combined integral formulations [START_REF] Chandler-Wilde | Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering[END_REF], which are uniquely solvable for all k but are more computationally expensive. However, for the purpose of the present paper, the choice of integral equation is not particularly important since our focus is to demonstrate the application of AMA to BEMs. Thus we choose to employ the simplest BIE formulation in order to keep the computational cost at a minimum.

Boundary element discretization and Fast Multipole acceleration

The main ingredients of the BEM are a transposition of the concepts developed for the Finite Element Method [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF]. First, the numerical solution of the boundary integral equation ( 4) is based on a discretization of the surface Γ into N E isoparametric boundary elements of order one, i.e., three-node triangular elements. Each physical element E e on the approximate boundary is mapped onto a reference element ∆ e via an affine mapping

ξ ∈ ∆ e → y(ξ) ∈ E e , 1 ≤ e ≤ N e .
∆ e is the reference triangle in the (ξ 1 , ξ 2 )-plane. The N interpolation points y 1 , . . . , y N are chosen as the vertices of the mesh. Each unknown field q is approximated with globally continuous, piecewise-linear shape functions (v i (y)) 1≤i≤N : v i (y j ) = δ ij for 1 ≤ i, j ≤ N . A boundary element E e contains exactly 3 interpolation nodes (y e k ) 1≤k≤3 associated with 3 basis functions (v e k ) 1≤k≤3 . These basis functions are related to the canonical basis (v k ) 1≤k≤3 defined on the reference element ∆ e by v e k (y(ξ)) = vk (ξ). Each unknown field is approximated on the element E e by q(y)

≈ 3 k=1 q k v e k (y),
where q k denotes the approximation of the nodal value of q(y k ). To discretize the boundary integral equation ( 4), we employ the collocation approach. It consists in enforcing the equation at a finite number of collocation points x. To have a solvable discrete problem, one has to choose N collocation points. The N approximation nodes thus defined also serve as collocation points, i.e., (x i ) 1≤i≤N = (y j ) 1≤j≤N . In addition, a standard Gaussian quadrature formula is used to evaluate the integral in (4). This discretization process transforms (4) into a square complex-valued linear system of size N of the form Aq = b,

where the (N )-vector q collects the degrees of freedom (DOFs) while the (N )-vector b arises from the imposed incident wave field. Assembling the full dense matrix A [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF] requires the computation of all element integrals for each collocation point, thus requiring an O(N 2 ) computational time and memory. In addition, BEM matrix equations such as [START_REF] Bakry | A new accurate residual-based a posteriori error indicator for the BEM in 2D-acoustics[END_REF] are here solved iteratively using the GMRES algorithm. With reference to [START_REF] Bakry | A new accurate residual-based a posteriori error indicator for the BEM in 2D-acoustics[END_REF], each GMRES iteration requires one evaluation of Aq for a given q, a task requiring a computing time of order O(N 2 ). To lower this O(N 2 ) complexity, which is unacceptable for large BEM models, fast BEM solution techniques such as the Fast Multipole Method (FMM) must be employed. The goal of the FMM is to accelerate the evaluation of the matrix-vector product Aq for a given q, required at each iteration of an iterative solver applied to the BEM-discretized equations. Moreover, the governing BEM matrix is never explicitly formed, which leads to a storage requirement well below the O(N 2 ) memory required for holding it. Substantial savings in both computational time and memory requirements are thus achieved. In general terms, the FMM exploits a reformulation of the fundamental solution in terms of products of functions of x and of y, so that (unlike in the traditional BEM) integrations with respect to y can be reused when the collocation point x is changed. On decomposing the position vector r = yx into r = (yy 0 ) + r 0 -(xx 0 ), where x 0 and y 0 are two poles and r 0 = y 0x 0 (Figure 1a) and invoking a well-known plane wave decomposition, the Helmholtz fundamental solution is written as [START_REF] Darve | The fast multipole method: numerical implementation[END_REF]:

Φ k (x, y) = lim L→+∞ ŝ∈S e ikŝ.(y-y0) G L (ŝ; r 0 ; k)e -ikŝ.(x-x0) dŝ, ( 8
)
where S is the unit sphere of R 3 and the transfer function G L (ŝ; r 0 ; k) is defined in terms of Legendre polynomials P p and spherical Hankel functions of the first kind h

p by: (1) p (k|r 0 |)P p cos(ŝ, r 0 ) .

G L (ŝ; r 0 ; k) = ik 16π 2 0≤p≤L (2p + 1)i p h
It can be shown that Expression [START_REF] Bažant | Three-dimensional harmonic functions near termination or intersection of gradient singularity lines: a general numerical method[END_REF] is valid only for well-separated sets of collocation and integration points clustered around poles x 0 and y 0 .

In the single-level FMM, a 3D cubic grid of linear spacing d embedding the whole boundary Γ is then introduced to drive the computation (see Figure 1b). The FMM basically consists in using decomposition [START_REF] Bažant | Three-dimensional harmonic functions near termination or intersection of gradient singularity lines: a general numerical method[END_REF], with the poles x 0 and y 0 being chosen as the cell centres of the cells C x and C y , whenever x and y belong to non-adjacent cubic cells (i.e., cells that do not share a corner, Figure 2). The treatment of such Fast Multipole (FM) contributions exploits the plane wave decomposition (8) of the fundamental solution, truncated at a finite L and in a manner suggested by its multiplicative form. One notes that q → Aq is a discretized version of of the single-layer integral operator q → Sq. The FM-BEM essentially consists in an accelerated method for the evaluation of discretized integral operators on given densities. The efficient evaluation of Sq needed in (4) is decomposed into three steps:

1. Evaluation of the multipole moments for each cell C y R(ŝ; C y ) = Γ e ikŝ.(y-y0) q(y)ds(y).

Application of the truncated transfer functions

L(ŝ; C x ) = Cy / ∈A(Cx) G L (ŝ; r 0 )R(ŝ; C y ).
3. Numerical evaluation of the integration over the unit sphere with a quadrature rule (Sq) F M (x) ≈ q w q e -ikŝq.(x-x0) L(ŝ q ; C x ) . Conversely, when x and y belong to adjacent cells, traditional BEM evaluation methods are used instead.

To improve further the computational efficiency of the FM-BEM, standard (i.e., non-FMM) calculations must be confined to the smallest possible spatial regions while retaining the advantage of clustering the computation of influence terms into non-adjacent large groups whenever possible. This is achieved by recursively subdividing cubic cells into eight smaller cubic cells. New pairs of non-adjacent smaller cells, to which plane wave expansions are applicable, are thus obtained from the subdivision of pairs of adjacent cells. The cell-subdivision approach is systematized by means of an octree structure of cells. At each level , the linear cell size is denoted d . The level = 0, composed of only one cubic cell containing the whole surface Γ, is the tree root. The subdivision process is further repeated until the finest level = ¯ , implicitly defined by a preset subdivision-stopping criterion (d ¯ ≥ d min ), is reached. Level-¯ cells are usually termed leaf cells. This is the essence of the multi-level FMM, whose theoretical complexity is O(N log N ) per GMRES iteration both for computational time and memory requirements. More details on the implementation of the FM-BEM for waves can be found in [START_REF] Chaillat | A multi-level fast multipole BEM for 3-D elastodynamics in the frequency domain[END_REF].

Regularity of the boundary data

A standard rule of thumb for wave approximation problems is that between 6 and 10 degrees of freedom per wavelength are required to achieve "engineering accuracy" with numerical methods [START_REF] Marburg | Discretization requirements: how many elements per wavelength are necessary?[END_REF]. The most popular way to distribute these degrees of freedom is over a uniform mesh of the surface. However, for calculations involving scatterers with ridges, it is well known that this accuracy can be severely diminished if uniform meshes are employed [START_REF] Juhl | A note on the convergence of the direct collocation boundary element method[END_REF]. Furthermore, the convergence of BEMs for these problems will be suboptimal. This is due to the low regularity of the boundary data at ridges. In particular, at a ridge, we expect the acoustic pressure u to behave as [START_REF] Juhl | A note on the convergence of the direct collocation boundary element method[END_REF] u ∼ r π α as r → 0, [START_REF] Betcke | Numerical estimation of coercivity constants for boundary integral operators in acoustic scattering[END_REF] where r is the distance from the ridge, and α is the exterior angle. Hence the flow velocity will behave as

∂u ∂n ∼ r π α -1 as r → 0. ( 11
)
At corners where two or more ridges meet, the singularities in u and ∂u/∂n are even more severe and do not have simple expressions such as [START_REF] Betcke | Numerical estimation of coercivity constants for boundary integral operators in acoustic scattering[END_REF] and [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF] (see, e.g., [START_REF] Bažant | Three-dimensional harmonic functions near termination or intersection of gradient singularity lines: a general numerical method[END_REF]).

In order to accurately approximate the solutions of such scattering problems and recover the optimal convergence order, the mesh must be appropriately refined towards the ridges. We note that the optimal order of convergence in the L 2 -norm for P 1 discretizations is O(N -1 ) where N is the number of DOFs. Optimal a priori refinement strategies can be devised via a consideration of the polynomial interpolation of functions with the appropriate singularities [START_REF] Apel | Graded mesh refinement and error estimates for finite element solutions of elliptic boundary value problems in non-smooth domains[END_REF][START_REF] Graham | Hybrid Galerkin boundary elements on degenerate meshes[END_REF]. However, for complicated scattering geometries with vertices, designing an optimal mesh is challenging and often infeasible. Hence adaptive meshing strategies have proven popular. In the next section we describe the adaptive strategy employed here.

Metric-based anisotropic mesh adaptation

In this section, we derive new estimates for the metric-based anisotropic mesh adaptation procedure in the case of 3D surface solutions. Further details on the approach derivation for volume solutions and examples of application can be found in [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF][START_REF]Continuous mesh framework part II: validations and applications[END_REF], or in the comprehensive review [START_REF] Alauzet | A decade of progress on anisotropic mesh adaptation for computational fluid dynamics[END_REF].

The goal of mesh adaptation is to find the optimal mesh that permits the minimization of the approximation error due to the discretization scheme, i.e., the BEM in our context. To date, most mesh adaptation approaches for the BEM rely explicitly on the local error estimation of the boundary approximation to drive the adaptation procedure. The derivation of appropriate local error estimators is a significant challenge owing to the non-locality of boundary integral operators. This difficulty is the main reason why adaptivity for the BEM is a much less well-explored research topic in comparison to adaptivity for the FEM where the relevant operators are local differential operators.

The metric-based mesh adaptation method employed in this paper is also linked to the error in the surface approximation but in a much less explicit way. In addition, it is not linked to the underlying boundary integral operators. Rather, our method is connected to the linear interpolation error of the (unknown) exact solution on the mesh (interpolation with higher order functions has not yet been considered) within a specified norm. Then, given a number of degrees of freedom, meshes which minimize this interpolation error are designed in an anisotropic setting. Via arguments relating the BEM approximation error to this interpolation error, such meshes are seen to be close to optimal for BEM solutions. Many authors have considered such a technique, e.g., [START_REF] Chen | Optimal anisotropic meshes for minimizing interpolation errors in L p -norm[END_REF]. However this is the first time it is used in the boundary element context. In what follows, we focus on the description of metric-based mesh adaptation in the BEM setting, i.e., we assume that the solution is only provided on the boundary Γ.

Metric-based anisotropic mesh adaptation for a surface solution

A convenient framework to generate anisotropic meshes is that of Riemannian metric spaces. A Riemannian metric space is defined by a metric tensor M(x) (where M is a symmetric positive definite matrix), i.e., is a smoothly varying function of the physical variable x. In the context of anisotropic mesh generation, each mesh vertex x has an assigned value M(x) which dictates the size and orientation of adjacent elements. By generating a unit mesh in the corresponding Riemannian metric space, we obtain an anisotropic mesh refined in Euclidean space according to the metric M. This is the fundamental idea of metric-based mesh adaptation as introduced in [START_REF] George | Creation of internal points in Voronoi's type method. Control adaptation[END_REF].

A surface mesh, associated with a Riemannian metric space M = (M(x)) x∈Γ is a triangulation of the surface Γ. A triangle K, which is defined by its edges {e i } 3 i=1 , is unit with respect to M if the length of each edge is unit in this metric, i.e., if

||e i || M := t e i Me i = 1, for i = 1, 2, 3. (12)
At this point it is important to note that in the BEM context, the mesh is 3D but the elements are 2D. We need to introduce the operator R K that restricts a 3D quantity to the 2D local element frame. R K can be deduced from the edges of K, such that

R K (M -1/2 ) = t t 1 M -1/2 t 1 t t 1 M -1/2 t 2 t t 1 M -1/2 t 2 t t 2 M -1/2 t 2 ,
where

t 1 = e 1 e 1 and t 2 = [ t e 1 e 1 ] e 2 -[ t e 1 e 2 ] e 1 [ t e 1 e 1 ] e 2 -[ t e 1 e 2 ] e 1 .
If the three edges of K are of unit length with respect to the metric tensor M, then the area of K computed in the Riemannian (|K| M ) or Euclidean (|K|) space is, respectively,

p 3 2 1 |K| M = √ 3 4
, and

|K| = √ 3 4 det(R K (M -1/2 )).
The existence of a conforming unit mesh in which each triangle is perfectly unit with respect to a given Riemannian metric space is not guaranteed in general. Hence the objective of seeking a unit mesh must be relaxed somewhat. We seek a quasi-unit mesh with respect to M, which is a mesh composed of quasi-unit triangles.

A triangle K is quasi-unit with respect to M if 1 √ 2 ≤ ||e i || M ≤ √ 2, for i = 1, 2, 3, and |K| M = √ 3 4 . ( 13
)
These values come from practical considerations, see discussions in [START_REF] Loseille | Chapter 10 -unstructured mesh generation and adaptation[END_REF]. Using the above notion of a quasi-unit mesh, the generation of anisotropic meshes is simplified. The mesh generator is only required to generate a quasi-unit mesh in the prescribed Riemannian metric space, and this is shown to be always possible in [START_REF] Loseille | Chapter 10 -unstructured mesh generation and adaptation[END_REF]. The generated adapted mesh is anisotropic in Euclidean space.

The principle of AMA theory is then to minimize the error coming from the linear interpolation of the (unknown) exact solution. In the BEM context, we cannot directly apply the results proposed for 2D or 3D volume solutions as even though we work in 3D, we only have the solution on a surface. In particular, we cannot recover the Hessian of the numerical solution which is needed for the AMA procedure. One alternative would be to extend the theory for 3D volume solutions by considering the covariant derivatives. This approach would involve complex developments of differential geometry. In addition, we only focus on a linear approximation of the surface. These developments are mainly relevant for higher order approximations. A more natural and simple approach in the BEM context is to extend the surface solution in a tubular neighborhood, i.e., for a collection of points whose distance to the surface is smaller than ε. More precisely, to any solution u s defined on the surface Γ, we associate the lifting u defined on the tubular neighborhood Γ ε by u(x) = u s (P(x))

where P(x) is the projection onto the surface Γ. The function u is now a 3D volume solution and we can define its Hessian.

Controlling the linear interpolation error

We suppose that the solution u is sufficiently smooth to possess a second order Taylor approximation. For 3D volume solutions, the quadratic approximation u q of u is given by its truncated Taylor expansion centered at point a:

u(x) ≈ u q (x) := u(a) + t ∇u(a)(x -a) + 1 2 t (x -a)H(a)(x -a),
where H is a symmetric matrix representing the Hessian of the extended solution u.

We now have all the tools to derive an estimation of the interpolation error. We follow the same ideas than in [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF]. We note Π h the discrete linear interpolation operator, i.e., Π h u is the linear interpolant of u on the mesh H N (with N nodes). The starting point is to show that the interpolation error in the L 1 norm of the solution u on a triangular element K is

||u -Π h u|| L 1 (K) = |K| 24 | 3 i=1 t e i H e i | ≤ |K| 24
3 i=1 t e i |H| e i , (14)
where the e i are the edges of the triangle K provided in the local element frame and for every symmetric matrix H, |H| denotes the positive symmetric matrix deduced from H by taking the absolute values of its eigenvalues. The proof is based on the evaluation of the pointwise interpolation error within the element, i.e., (u -Π h u)(x) for x ∈ K. It only uses the mapping onto the reference element and the knowledge that the solution is a quadratic function. This error is then integrated over the element.

Proof: The reference element K ref is defined by its three vertices:

x 1 = t (0, 0), x 2 = t (1, 0), x 3 = t (0, 1).
The mapping onto the current element K is given by:

x = x 1 + B K x with B K = [x 2 -x 1 , x 3 -x 1 ] = [e 1 , e 2 ], x ∈ K, x ∈ K ref ,
where we have introduced the edges:

e 1 = x 2 -x 1 , e 2 = x 3 -x 1
, and e 3 = x 3x 2 . The quadratic function u reads in the frame of K ref :

u(x( x)) = 1 2 t x 1 H x 1 + 1 2 t x 1 H B K x + 1 2 t x t B K H x 1 + 1 2 t x t B K H B K x.
As we consider linear interpolation, the linear and constant terms of u(x( x)) are exactly interpolated such that only quadratic terms are kept. We introduce ũ(x) = 1 2

t x t B K H B K x, then it becomes (u -Π h u)(x) = (ũ -Π h ũ)(x).
We rewrite ũ in matrix form as

ũ(x( x)) = 1 2 t x ŷ t e
u(v i ) = u(v i ) for all i ∈ [1, 3], with x i the vertex i:        Π h ũ(x 1 ) = ũ(x(0, 0)) = 0 = c, Π h ũ(x 2 ) = ũ(x(1, 0)) = 1 2 ( t e 1 He 1 ) = a, Π h ũ(x 3 ) = ũ(x(0, 1)) = 1 2 ( t e 2 He 2 ) = b.
The solution of the previous linear system gives the final expression of Π h ũ:

Π h ũ(x( x)) = 1 2
( t e 1 He 1 ) x + ( t e 2 He 2 ) ŷ .

The exact pointwise interpolation error is then given by

(u -Π h u)(x( x)) = 1 2 [ ( t e 1 He 1 ) (x 2 -x) + ( t e 2 He 2 ) (ŷ 2 -ŷ) + 2 ( t e 1 He 2 ) xŷ ].
Classically, for every function F , its integration over K can be computed through its expression in K ref :

K F (x) dxdy = K ref F (x( x)) ||e 1 ∧ e 2 || dxdŷ = 2|K| K ref F (x( x)) dxdŷ.
Consequently, the interpolation error in the L 1 norm, evaluated by a direct integration of |(u -

Π h u)(x)|, is given by u -Π h u L 1 (K) = 2|K| 24 t e 1 H e 2 -t e 1 H e 1 -t e 2 H e 2 .
The cross term can be expressed only in terms of e i H e i for i = 1, .., 3 by remarking that e 1e 2 + e 3 = 0. It follows

2 t e 1 H e 2 = t e 1 H e 1 + t e 2 H e 2 -t e 3 H e 3 ,
and we deduce:

2 t e 1 H e 2 -( t e 1 H e 1 + t e 2 H e 2 ) = 3 i=1 t e i H e i .
When all the eigenvalues of H are greater than or equal to zero, we have

3 i=1 t e i H e i = 3 i=1 t e i H e i .
On the other hand, when H has positive and negative eigenvalues, we have only the following inequality

3 i=1 t e i H e i ≤ 3 i=1 t e i |H| e i .
This property is used to conclude the proof in the general case.

We consider now all elements K which are unit with respect to the metric M. The next step consists in introducing the metric M, i.e, in showing that

3 i=1 t e i H e i = 3 2 trace(R K (M -1 2 H M -1 2 )).
Proof: We consider first the simple case where R K (H) = I 2 and R K (M) = I 2 (where I 2 is the 2 × 2 identity matrix).

The regular triangle K 0 = (x 1 , x 2 , x 3 ) unit for I 2 is defined by the vertices:

x 1 = (0, 0) , x 2 = (1, 0) , x 3 = 1 2 , √ 3 2 .
We first show the following preliminary result: For every line (D) passing through one of the vertices of K 0 , the sum of the square lengths of the edges projected on (D) is invariant.

Without loss of generality, we assume that (D) passes through the vertex x 1 of K 0 . If (D) is defined by the vector n = (cos(v), sin(v)) , with v ∈ R, then the length of the three edges of K 0 projected on (D) are given by:

a = e 1 . n = cos(v), b = e 2 . n = 1 2 cos(v) + √ 3 2 sin(v), c = e 3 . n = -1 2 cos(v) + √ 3 2 sin(v). Such that a 2 + b 2 + c 2 = 3/2. When R K (M) is different from I 2 , we use R K (M - 1
2 ) that maps the unit ball of I 2 onto the unit ball of R K (M). (D) is now directed by one of the eigenvectors of R K (M), e.g., v j , and which is passing through x 1 . The lengths a, b and c are thus multiplied by h j , i.e., the size prescribed by R K (M) in the direction v j . Consequently, the square length of the edges projected on (D) are multiplied by h 2 j . It comes using the preliminary result:

3 i=1 ( t e i v j ) 2 = 3 2 h 2 j .
Considering this relation for the two principal directions of R K (M), we obtain:

3 i=1 e i 2 2 = 2 j=1 3 i=1 ( t e i v j ) 2 = 3 2 h 2 1 + h 2 2 = 3 2 trace(R K (M) -1 ). ( 15
)
We now consider the general case where a symmetric matrix R K (H) is involved in the estimation instead of R K (H) = I 2 . R K (H) being symmetric, it has two real eigenvalues (µ i ) i=1,2 along the principal directions (u i ) i=1,2 . If K is a regular triangle with edges (e i ) i=1,3 , according to the preliminary result, the sum of the projected square length of edges (e i ) i=1,3 on each principal direction u j is equal to 3/2. We deduce:

µ j 3 i=1 ( t e i u j ) 2 = 3 2 µ j , for j ∈ [1, 2] and 2 j=1 3 i=1 µ j ( t e i u j ) 2 = 3 2 (µ 1 + µ 2 ).
The previous equality reads:

3 i=1 t e i He i = 3 2 trace(R K (H)).
We deduce the general case by taking R K (M -1 2 HM -1 2 ) as a symmetric matrix. In that case, each edge e i of the regular triangle is mapped on ẽi = R K (M -1

2 )e i . The new triangle defined by edges (ẽ i ) i=1,3 is unit with respect to R K (M). This concludes the proof.

Using that |K| = √ 3 4 det(R K (M -1/2
)), the interpolation error ( 14) is finally expressed for a surface solution by: Importantly, [START_REF] Chaillat | Recent advances on the fast multipole accelerated boundary element method for 3D time-harmonic elastodynamics[END_REF] implies that the interpolation error does not depend on the element shape, i.e., on the mesh generator. The metric alone contains enough information to describe completely the linear interpolation error in the L 1 norm. In addition, this interpolation error can be expressed with continuous quantities only. Indeed M and H can be defined continuously anywhere. The restriction R K is replaced by its continuous counterpart R s , i.e., the tangent plane restriction. For a smooth surface, the operator R s is defined by using the local Frenet frame [t 1 , t 2 , n] (see Fig. 3). We thus define a continuous linear interpolation error, independent from the discrete mesh similarly to [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF] by removing the surface contribution of ( 16)

||u -Π h u|| L 1 (K) = √ 3 64 det(R K (M -1 2 )) |trace(R K (M -1 2 H M -1 2 ))|. (16)
2 ||u -Π h u|| L 1 (K) |K| = 1 8 |trace(R s (M -1 2 H M -1 2 ))| ≤ 1 8 trace(R s (M -1 2 |H| M -1 2 ))
Then, we define the upper bound of the continuous interpolation error by

e Rs(M) (x) = 1 8 trace(R s (M -1 2 |H| M -1 2 )). (17)
Finally, the continuous interpolation error calculated on the continuous mesh M is given by

E p (M) = Γ e Rs(M) (x) p dx, p ≥ 1. ( 18
)

Numerical validation on an analytical example

In this section, we validate the previous continuous analysis and the assumption to work in the tangent plane. To this end, we compare the continuous evaluation of the interpolation with the discrete evaluation in a simple analytic case. For the surface, we consider a quarter of a cylinder (S) defined by the following parametric equations:

  x y z   = σ(u, v) =   cos(u) sin(u) v   , with u ∈ [0, π/2] and v ∈ [0, 1]
. The Frenet frame is then given by n = t (x, y, 0), t 1 = t (0, 0, 1) and t 2 = t (-y, x, 0), so that the restriction R s is defined everywhere. We also consider a continuous description of the mesh given by the following metric:

M(x, y, x) = α   h -2 x h -2 y h -2 z   ,
with h x (x, y, z) = 0.2, h y (x, y, z) = 0.1(y + 1) + 0.05(x -1) and h z (x, y, z) = 0.2(z + 2). In the continuous mesh setting, we introduce the complexity N of a mesh. It is the continuous counterpart of the number of vertices and it governs the size of the meshes, i.e.,

N := Γ d(x)dx = Γ det(R s (M)(x)) dx,
where d = (h 1 h 2 ) -1 is the density and (h i ) i=1,2 are the local sizes along the principal directions of R s (M). In this analytical example, α is a parameter that controls the complexity of the continuous mesh.

Given the quadratic function u(x, y, z) = x 2y 2 + xz, we can compute the right hand side of (18) using a formal integration. We have : S e Rs(M) (a)dS = 0.01038536836π 8α .

We then compare this continuous estimation with a discrete one based on a sequence of meshes generated with α = [100, 500, 1000, 2000, 4000]. The interpolation error is computed with a numerical integration. Figure 4a represents the mesh for α = 100. In Figure 4b, we plot the continuous and discrete errors. We also give an upper and lower bound as the mesh generator is not able to generate a perfectly unit-mesh but only quasi-unit meshes. An edge is considered to be quasi-unit when its length computed in the metric is in

[1/ √ 2, √ 2
]. We observe that the discrete interpolation error fits the continuous one, and that the deviation is always contained within the lower and upper bounds (corresponding to the generation of a quasi-unit mesh instead of a perfectly unit one). This example confirms the relevance of the use of the tangent plane restriction. We now have all the tools to derive the mesh minimising the continuous interpolation error.

Deriving the optimal continuous mesh

The introduction of the continuous interpolation is very important to find the optimal mesh, i.e., the mesh minimizing the error introduced in the BEM solution. Indeed, the initial discrete optimization problem reads:

Find the surface mesh H opt N with N nodes, minimising the discrete interpolation error ||u -

Π h u|| L p (Γ h ) , (P d )
where Γ h is the discrete support representing Γ. This problem is generally intractable practically and would require some simplifications. By replacing the discrete interpolation error with its continuous equivalent, we do not have anymore to deal with discrete quantities. Using the mathematical tools defined on the continuous mesh space instead of the discrete one, the optimization problem is simplified. It is the originality of this work (similarly to [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF]), to solve this problem in a continuous setting. We consider in the following the continuous interpolation error controlled in L p norm (with p ≥ 1). The global optimization problem becomes to find the optimal continuous mesh M * minimizing the continuous interpolation error in the L p norm:

Given a complexity N , find

M * = argmin M Γ e Rs(M) (x) p dx 1/p . (P c )
The optimization problem (P c ) is well-posed and the interpolation error is expressed continuously without the need to define a discrete support. In the following section, (P c ) is solved in two steps.

We need first to introduce some concepts of the continuous mesh framework (detailed in [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF]). These meshes are represented by a Riemannian metric space

M = (R s (M)(x)) x∈Γ . Locally, the spectral decomposition of R s (M)(x) is R s (M)(x) = R(x)Λ(x) t R(x),
where R(x) is an orthonormal matrix that provides the local orientation given by the eigenvectors (v i (x)) i=1,2 . The (λ i (x)) i=1,2 are the corresponding eigenvalues. As a result, the diagonal matrix Λ is either

λ 1 (x) 0 0 λ 2 (x) or h -2 1 (x) 0 0 h -2 2 (x)
.

Introducing the density d = (h 1 h 2 ) -1 and the anisotropic quotients

r i = h 2 i (h 1 h 2 ) -1 , we can write R s (M)(x) = d(x)R(x) r -1 1 r -1 2 t R(x).
It follows from ( 17) that

e Rs(M) = 1 8 d -1 (a) 2 i=1 r i (a) t v i (a)|R s (H)|v i (a) ,
where v i denote again the eigenvectors of the metric R s (M).

Optimal density and anisotropic ratios for a fixed orientation. It is now possible to solve the continuous optimization problem. We introduce γ i = t v i (a)|R s (H)|v i (a), i.e., the directions of the continuous mesh. The first step in the solution of (P c ) is to fix γ 1 and γ 2 . Since the density is constant, we can write r 2 ∼ r -1 1 . The function e Rs(M) locally writes:

e Rs(M) = d -1 8 r 1 γ 1 + r -1 1 γ 2 ,
such that, we have to solve:

min (r 1 , d) Γ d -p r 1 γ 1 + r -1 1 γ 2 p
, under the linear constraint:

Γ d = N ,
where d is now an unknown of the optimization problem and the constraint is linear. In addition, the problem is uncoupled. It is thus natural to exhibit first the optimal anisotropic quotient r 1 and the optimal density in a second step.

To minimize E p (M) defined in [START_REF] Chaillat | Fast iterative boundary element methods for high-frequency scattering problems in 3D elastodynamics[END_REF], we use the calculus of variations. The classical Euler-Lagrange necessary condition states that for the variation δM = (δr 1 , δd), we have:

∀δr 1 , ∀δd with Γ δd = 0, δE p (M; δr 1 ) + δE p (M; δd) = 0. ( 19
) If we note β = r 1 γ 1 + r -1 1 γ 2 , ( 19
) for δd = 0 leads to δE p (M; δr 1 ) = Γ pd -p β p-1 γ 1 -r -2 1 γ 2 δr 1 = 0. (20)
It follows r 1 = γ 2 /γ 1 , r 2 = γ 1 /γ 2 and e Rs(M) becomes

e Rs(M) = d -1 8 r 1 γ 1 + r -1 1 γ 1 = 1 4 d -1 (γ 1 γ 2 ) 1 2 . ( 21
)
Then [START_REF] Chaillat | Theory and implementation of H-matrix based iterative and direct solvers for Helmholtz and elastodynamic oscillatory kernels[END_REF] for δr 1 = 0 leads to

δE p (M; δd) = Γ p 1 4 p d -p-1 (γ 1 γ 2 ) p 2 δd = 0 with the contraint Γ δd = 0. (22)
A condition to ensure ( 22) is given by d -p-1 (γ 1 γ 2 ) p 2 = C where C is a real constant. Using the constraint of the constant complexity N , the final expression of d is

d = N Γ (γ 1 γ 2 ) p 2(p+1) -1 (γ 1 γ 2 ) p 2(p+1) . (23)
Finally, the optimal continuous mesh M * obtained by solving (P c ) is described by

λ * i = (h * i ) -2 = N Γ (γ 1 γ 2 ) p 2(p+1) -1 (γ 1 γ 2 ) -1 2(p+1) γ i , and r * i = (γ 1 γ 2 ) 1 2 γ i ∀i = 1, 2 . (24)
Uniqueness. It is then possible to prove that the optimal continuous mesh M * given by ( 24) is the unique solution of (P c ) verifying E p (M * ) p < E p (M) p for all M having the same fixed (γ i ) i=1,2 .

Proof: Evaluating E p to the power p at the optimal point M * , it becomes

E p (M * ) p = 1 4 p N -p Γ (γ 1 γ 2 ) p 2(p+1) p+1 . (25)
A generic M of complexity N is given by its two anisotropic quotients (r i ) i=1,2 and its density d. To take into account the constraint on the density, the density is rewritten as d = N ( Γ f ) -1 f , where f is a strictly positive function. The error committed with a generic M is then

E p (M) p = N -p 8 p Γ f p Γ f -p (r 1 γ 1 + r 2 γ 2 ) p .
To prove E p (M * ) p < E p (M) p , we use the generalized arithmetic-geometric inequality which comes from the strict concavity of the logarithm function (since r 1 γ 1 = r 2 γ 2 , otherwise it would be the optimal continuous mesh):

ln 1 2 (r 1 γ 1 + r 2 γ 2 ) > 1 2 (ln (r 1 γ 1 ) + ln (r 2 γ 2 )) = ln r 1 2 1 γ 1 2 1 r 1 2 2 γ 1 2 2 = ln γ 1 2 1 γ 1 2 2 =⇒ (r 1 γ 1 + r 2 γ 2 ) p > 2 p γ 1 2 1 γ 1 2 2 p ,
where we have used that r 1 r 2 = 1. Finally, if we denote g = (γ 1 γ 2 ) p 2(p+1) , we have

         E p (M * ) p p+1 = 1 4 p p+1 N -p p+1 Γ g , E p (M) p p+1 > 1 4 p p+1 N -p p+1 Γ f p p+1 Γ f -p g p+1 1 p+1
.

Using the Hölder inequality, we have that

Γ f p p+1 Γ f -p g p+1 1 p+1 = Γ f p p+1 Γ g f p p+1 p+1 1 p+1 ≥ Γ f p p+1 g f p p+1 = Γ g , (26)
as 26) implies E p (M * ) < E p (M) for all M having the same fixed γ 1 , γ 2 . The optimal solution M * is thus unique.

p + 1 p ≥ 1, p + 1 ≥ 1 and p p + 1 + 1 p + 1 = 1 . Relation (
Optimal orientations. The last step consists in deriving the optimal directions of the continuous mesh M * . We know from ( 21) and ( 23) that for a given set of directions (v i ) i=1,2 , the optimal interpolation error reads:

E p (M * ) = N -1 4 Γ (det(G)) p 2(p+1) p+1 p
, where G is the diagonal matrix composed of γ 1 and γ 2 . It is then possible to minimize the optimal interpolation error by seeking the directions v i minimizing det(G) p 2(p+1) , or equivalently minimizing det(G). Geometrically, det(G) is the area of the rectangle defined by (v i ) i=1,2 and of length ( t v i |R s (H)|v i ) i=1,2 . We denote by (λ i ) i=1,2 the eigenvalues of |R s (H)| and (u i ) i=1,2 its principal directions. The length of a side computed with respect to |R s (H)| is:

t v i |R s (H)| v i = λ 1 t v i u 1 2 + λ 2 t v i u 2 2 .
We denote by θ the angle between u 1 and v 1 (see Fig. 5). The problem is now to minimize

( t v 1 |R s (H)| v 1 )( t v 2 |R s (H)| v 2 ) = (λ 1 cos 2 θ + λ 2 sin 2 θ)(λ 1 sin 2 θ + λ 2 cos 2 θ) = (λ 1 -λ 2 ) 2 sin 2 θ cos 2 θ.
It results that det(G) is minimal when the vectors (v i ) i=1,2 are aligned with the principal directions of |R s (H)|. Using [START_REF] Darbas | Combining analytic preconditioner and fast multipole method for the 3-D Helmholtz equation[END_REF], we conclude that the optimal metric both in sizes and directions is given by the following theorem.

✓ u 1 u 2 v 2 v 1
Figure 5: Determination of the optimal orientations of the continuous mesh.

Theorem 1. Let u be a twice continuously differentiable function and H its Hessian, the optimal continuous mesh M * minimizing Problem (P c ) reads locally:

R s (M * ) = D * det(|R s (H)|) -1 2(p+1) |R s (H)|, with D * = N Γ det(|R s (H)|) p 2(p+1) -1 . ( 27
)
M * is unique, locally aligned with the eigenvectors basis of R s (H) and has the same anisotropic quotients as R s (H).

In addition, M * provides an optimal explicit bound of the interpolation error in the L p norm (with p ≥ 1):

E p (M * ) = N -1 4 Γ det (|R s (H)|) p 2(p+1) p+1 p . ( 28
)
We now have all the tools to generate in practice a discrete mesh approximating the continuous optimal solution. We can use any metric-based adaptive mesh generators as soon as the generated mesh is quasi-unit. In addition, we provide bounds on the discrete interpolation error when the continuous mesh is projected into the space of discrete meshes by means of an adaptive mesh generator. If the mesh generator achieves a quasi-unit mesh with respect to M * , the following bounds for the discrete interpolation error are obtained:

1 2 E p (M * ) ≤ u -Π h u L p (Γ h ) ≤ 2 E p (M * ) ,
or equivalently,

1 8 N -1 Γ det (|Rs(H)|) p 2(p+1) p+1 p ≤ u -Π h u L p (Γ h ) ≤ 1 2 N -1 Γ det (|Rs(H)|) p 2(p+1) p+1 p . ( 29
)
In other words, M * defined by [START_REF] Feischl | Adaptive Boundary Element Methods: A Posteriori Error Estimators, Adaptivity, Convergence, and Implementation[END_REF] allows us to generate an optimal adapted mesh to control the interpolation error in the L p norm. One main advantage of this approach is to be completely generic, independent from the PDE and the numerical scheme.

Application to solutions given by a BEM approximation

The above results are valid for the linear interpolant Π h u of the exact solution u. However, neither Π h u or u are known. Only the approximation of u obtained by the BEM is at hand. In order to relate the interpolation error to the approximation error, we consider a reconstruction operator R h which is applied to the numerical approximation u h . This operator can be a recovery process, a hierarchical basis, or an operator connected to an a posteriori estimate. We make the assumption that the reconstruction R h u h is better than u h for a given norm || • ||, i.e., that

||u -R h u h || ≤ α||u -u h ||, where 0 ≤ α < 1.
Then, using the triangle inequality, we have

||u -u h || ≤ 1 1 -α ||u h -R h u h ||.
If Π h R h φ h = φ h for all φ h ∈ P 1 , i.e., the reconstruction operator is exact at the nodes, then

||u -u h || ≤ 1 1 -α ||R h u h -Π h R h u h ||.
a

Fictive points Surface points

Figure 6: Fictive points created to ensure the tubular neighborhood of a. The solution at blue node is set to the surface solution.

Then we have from ( 29) that

||u -u h || ≤ 1 2(1 -α) N -1 Γ det(|R s (H R h u h )|) p 2(p+1) p+1 p
.

Consequently we apply Estimate [START_REF] Ferraz-Leite | Simple a posteriori error estimators for the h-version of the boundary element method[END_REF] to R h u h . There are numerous techniques in the literature to approximate or "recover" the Hessian (see, e.g., [START_REF] Picasso | A numerical study of some hessian recovery techniques on isotropic and anisotropic meshes[END_REF][START_REF] Vallet | Numerical comparison of some hessian recovery techniques[END_REF]). From a practical point of view, we only need to recover a pointwise 3D Hessian at the nodes of the current surface mesh. The method we employ here is based on a local quadratic recovery and uses the boundary solution only. The lifting solution is deduced from u s , following the assumption of Section 3.1. The starting point is the quadratic approximation of the lifting function u at a point a = (x a , y a , z a ) located on the surface of the tubular neighborhood (such that we have P(a) = a):

u(x) ≈ u s (a) + t ∇u(a)dx + 1 2 t dx H(a) dx,
where dx = xa. The gradient of u being a 3D vector, it reads

t ∇ u(a) = (a b c).
Similarly, the Hessian is a 3 × 3 symmetric matrix:

H(a) =   d e f e g h f h i   .
At each surface mesh node, we thus have 9 unknowns to recover. The simple approach we use is a least square. We start with the patch composed of the elements adjacent to the node a, see Figure 6. The set of points is then augmented by creating points along the normal direction to the surface. This extruded geometry is depicted in Figure 6 and uses standard boundary layer mesh generation procedures [START_REF] Loseille | Robust Boundary Layer Mesh Generation[END_REF] (especially at complex corners or along ridges). The normal spacing is set to the minimal edge length of the patch of the current point. If the number of unknowns of this augmented system is not sufficient, an additional layer of neighbours is taken into account.

Refinement strategy

All the numerical examples shown in the following section use the same refinement strategy. Given a user-prescribed sequence of complexity N = [N 1 , . . . , N k ], we seek for the sequence of corresponding optimal meshes H = [H 1 , . . . , H k ]. This process is non-linear, i.e., both the mesh and the solution have to be converged. The following iterative algorithm is used to generate H i :

Start from the mesh H 0 i = H i-1 (or from the initial uniform mesh H 0 at iteration 0). for all k=0 : n iter do

Compute the approximation u h on the mesh H k i with the iterative FM-BEM solver (software COFFEE). Compute the recovered Hessian from u h and deduce Metric [START_REF] Feischl | Adaptive Boundary Element Methods: A Posteriori Error Estimators, Adaptivity, Convergence, and Implementation[END_REF] at the nodes of the mesh, with complexity N i . From [START_REF] Feischl | Adaptive Boundary Element Methods: A Posteriori Error Estimators, Adaptivity, Convergence, and Implementation[END_REF], a new quasi-unit mesh H k+1 i of complexity N i is generated with FEFLO.A [START_REF] Loseille | Chapter 10 -unstructured mesh generation and adaptation[END_REF]. end for For each case, the complexity sequence is always defined as N i+1 = 2 N i with n iter = 5.

Validation and enhancements of the iterative FM-BEM solver

In this section we first validate the efficiency of the FM-BEM solver at hand, COFFEE. We also describe several necessary enhancements to improve the efficiency of the FM-BEM solver when anisotropic elements are present in the mesh. This concerns the change of the initial guess for the GMRES solver along with the design of a simple but efficient preconditioner.

Validation of the fast BEM solver

The fast BEM solver for 3D Helmholtz problems, used in this article has been implemented in our in-house solver COFFEE [START_REF] Chaillat | A multi-level fast multipole BEM for 3-D elastodynamics in the frequency domain[END_REF]. Since so far the code has only been validated for 3D elastodynamic problems, this section is devoted to the verification of the accuracy and numerical efficiency of the solver for the scalar acoustic case of interest here. To this end, we consider the simple case of the diffraction of an incident plane wave by a unit sphere for which an analytical solution is known [START_REF] Bowman | Electromagnetic and acoustic scattering by simple shapes[END_REF].

In this section, we consider a uniform triangulated mesh with a density of points per wavelength equal to 10. The tolerance of the GMRES iterative solver is set to 10 -4 . In all the numerical examples presented, we use the CERFACS GMRES package [START_REF] Frayssé | Algorithm 842: A set of gmres routines for real and complex arithmetics on high performance computers[END_REF]. We report in this section only results for the EFIE with Dirichlet boundary conditions and similar results are obtained for the MFIE with Neumann boundary conditions. In Table 1, we report the number of degrees of freedom N together with the number of iterations required to achieve convergence in GMRES for both the standard BEM and the FM-BEM solvers. The two solvers are seen to converge with the same rate. Importantly, we were not able to solve problems larger than about 50 000 DOFs with the standard BEM due to time and memory constraints. We also report the relative L 2 -norm of the discrepancy between the computed solution and the analytical one, computed in the domain (on a circle of radius 2 enclosing the obstacle). Both solvers give the same accuracy. error (FM-BEM) 3.66 10 -3 2.01 10 -3 1.12 10 -3 8.09 10 -4 7.65 10 -4 Table 1: Validation of the FM-BEM solver for 3D Helmholtz equations: comparison of the accuracy and convergence of the iterative BEM and FM-BEM solvers.

In addition, to check the numerical efficiency of the solver, we report in Figure 7 the theoretical and numerical complexities in terms of computational time for one matrix-vector product. The black circles represent the computational time for the standard BEM and the red squares represent the time for the FM-BEM. As expect, the standard BEM is seen to be of the order of O(N 2 ) operations per iteration and the FM-BEM of the order of O(N log N ) operations per iteration.

Enhancements of the FM-BEM for anisotropic meshes and adaptivity

One important point in the adaptive process is the ability of the iterative solver to converge to the solution within the desired threshold. For small-scale examples, the FM-BEM is not necessary (only the standard BEM) and the number of GMRES iterations required to solve the matrix system is not a pressing issue. However, as soon as largescale problems are considered, the use of the FM-BEM is mandatory to compute the solution on a given mesh. The number of GMRES iterations is known to increase with the size of the problem [START_REF] Chaillat | Fast iterative boundary element methods for high-frequency scattering problems in 3D elastodynamics[END_REF]. The presence of anisotropic elements in the mesh also increases the number of GMRES iterations. As a result, the number of iterations may be a limiting factor when employing the FM-BEM on anisotropic meshes. Here we discuss two simple methods to reduce the number of GMRES iterations in the context of a mesh adaptation strategy. The first one is a diagonal scaling preconditioner and the second one is based on the use of an initial guess for the GMRES obtained from the previous level of adaptation. To illustrate the efficiency of these preconditioners, we consider the spherical resonator test case: a screen curved into a sphere with a circular aperture to create a resonator. The sphere has unit radius and the circular opening of radius π/10 has its centre on the positive x-axis, see Figure 13a. Figure 7: Comparison of the theoretical and numerical complexities in terms of computational time for one matrix-vector product. The black circles represent the computational time for the standard BEM and the red squares represent the time for the FM-BEM. In addition, we report the two theoretical complexities.

Preconditioning with a diagonal scaling. When performing BEM calculations on adapted meshes containing elements with drastically different sizes, the diagonal and off-diagonal entries will also vary drastically. A simple way to avoid this issue is to introduce a preconditioner which scales all the diagonal entries of the system matrix to 1. We first write the matrix system to be solved as Ax = b. We propose the following right-preconditioning strategy with the preconditioner M defined as

M ij = A ij , if i = j, 0, if i = j. (30)
Therefore the preconditioned system to be solved becomes:

AM -1 y = b, x = M -1 y. (31)
We report in Table 2 the number of GMRES iterations required to achieve convergence with a tolerance of 10 -4 with and without this preconditioner. We see that the gain in terms of number of iterations and computational times increases with the sizes of the meshes, to reach a reduction by a factor 4. Table 2: Spherical resonator: Number of GMRES iterations required to achieve convergence with a tolerance of 10 -4 both without and with the diagonal scaling preconditioner [START_REF] Gantumur | Adaptive boundary element methods with convergence rates[END_REF].

#DOFs

Using an appropriate initial guess for GMRES. A further advantage of an adaptive mesh refinement is the possibility to use a computed approximation at a previous level of refinement as an initial guess for the GMRES solver. This choice is expected to be more efficient than a guess chosen a priori, such as the zero vector or the incident wave evaluated on the surface. The simple modification consists in interpolating the previous solution onto the newly generated mesh (see the refinement strategy, Section 3.6). In our case, this is done as an additional post-processing step in FEFLO.A. In Table 3, we illustrate the improvement obtained with such an initial guess. More precisely, we consider the same resonator problem as above and solve the system at each step of refinement twice, once with the zero vector as the initial guess and once with the linear interpolation of the previous step numerical solution onto the current mesh as the initial guess. Again the tolerance for the GMRES solver is 10 -4 . The gain is less impressive than the one obtained with the diagonal scaling preconditioner, but an additional gain of 10% is obtained. To sum up, the preconditioning by scaling the diagonal entries of the matrix yields a large reduction in the number of GMRES iterations. Furthermore, this reduction is greater for meshes with larger anisotropic ratios. Next we saw that the initialization of GMRES using the previous solution also serves to save a number of iterations, a saving which increases with the size and degree of anisotropy of the mesh.

Validation of the adaptive mesh strategy for the FM-BEM

We illustrate in this section the capabilities of the adaptive mesh strategies for different 3D configurations: scattering of plane waves by open surfaces (a planar screen and a sphere with aperture), by a closed simple shape (a cube with cavity [START_REF] Darbas | Combining analytic preconditioner and fast multipole method for the 3-D Helmholtz equation[END_REF]), and by a large complex structure (an aircraft).

Planar screen

We consider the scattering of plane waves (with unit amplitude and direction d = (-1, 0, 0)) by a sound-soft unit screen centered at the origin, i.e., the screen occupies the region -1 2 , 1 2 × -1 2 , 1 2 . Such configurations lead to surface solutions with the most severe singular behavior, and are a good illustration of the capabilities of an adaptive mesh strategy.

It is well known (see, e.g., [START_REF] Juhl | A note on the convergence of the direct collocation boundary element method[END_REF]) that at the edge of a screen, we have ∂u ∂n ∼ (kr) -0.5 as kr → 0, [START_REF] Graham | Optimal panel-clustering in the presence of anisotropic mesh refinement[END_REF] where r is the distance from the edge. At the corner of the screen, this singularity is more severe, and was shown in [START_REF] Bažant | Three-dimensional harmonic functions near termination or intersection of gradient singularity lines: a general numerical method[END_REF] to take the form ∂u ∂n ∼ (kr) -0.704 as kr → 0.

The optimal convergence rate here (for our P 1 discretization) is O(N -1 ). Owing to the singular behavior, the approximated solution on a sequence of uniformly refined meshes should yield a sub-optimal convergence rate. We employ the adaptive strategy proposed in Section 3.6 for three wavenumbers (k = 5, 20, 80) and compare the results to those obtained with a sequence of uniformly refined meshes. In figures 8b-d, we present the relative L 2 errors in the scattered field evaluated on a circle in the (x, y)-plane of radius 2 (Figure 8a).

We observe that the use of the anisotropic adaptive strategy improves the computational efficiency for this scattering problem by almost an order of magnitude. For example, for k = 5, to achieve an accuracy of 1.37 10 -2 with a uniform refinement 3 622 DOFs are required while the anisotropic refinement produced an approximation with an accuracy of 1.78 10 -2 with only N = 271 DOFs. For k = 20, to achieve an accuracy of 1.12 10 -2 with a uniform refinement 13 929 DOFs are required while the anisotropic refinement produced an approximation with an accuracy of 1.06 10 only 801 DOFs. Finally, for k = 80, to achieve an accuracy of 1.12 10 -2 with a uniform refinement 55 224 DOFs are required while the anisotropic refinement produced an approximation with an accuracy of 7.9 10 -3 with only 4 672 DOFs. In addition, for all three wavenumbers the uniform refinement yields a convergence rate of approximately O(N -0.5 ) in the relative L 2 error. This sub-optimality is related to the prominent edge singularities of the form [START_REF] Graham | Optimal panel-clustering in the presence of anisotropic mesh refinement[END_REF]. While, with the adaptive mesh strategy, the optimal convergence rate of O(N -1 ) is recovered. In fact, it appears that we even improve slightly on this rate, especially at the higher wavenumbers. This improvement can be attributed to the exploitation of the directionality of the wave propagation.

O 1 √ N O 1 N N (#DOFs) L 2 error in domain Uniform Adapted (c) (d)
In Figure 9, we depict the final surface mesh for k = 5 and k = 20. For k = 5 (resp. k = 20), the mean anisotropic ratio is around 30 (resp. 60). If an isotropic error estimate was used, the obtained mesh would contained 55 608 DOFs, i.e., 15 times more DOFs than the anisotropic mesh for k = 20 at 5th iteration. This also emphasizes the advantage of using anisotropic estimates over isotropic estimates.

Finally, the level of anisotropy on the sides of the screen increases with the complexity, see Figure 10, showing the strong directionality of the error.

Cube with cavity

We next consider the scattering by a cube with cavity [START_REF] Darbas | Combining analytic preconditioner and fast multipole method for the 3-D Helmholtz equation[END_REF] which is the 3D analogue of the 2D trapping domain [START_REF] Betcke | Numerical estimation of coercivity constants for boundary integral operators in acoustic scattering[END_REF]. We use the adaptive strategy for the sound-soft problem and the incident direction is d i = (-2, -1, -1) throughout. The singularities in the boundary solution for this problem should be less severe than for the screen. At an edge, the unknown should behave as ∂u ∂n ∼ r -1/3 .

The general rule is that we expect the acoustic pressure u, to behave near an edge as u ∼ r π/α as r → 0, where r is the distance from the edge and α is the exterior angle subtended at the edge. Hence the flow velocity will behave as ∂u ∂n ∼ r π/α-1 as r → 0.

Again, the errors for the adapted and uniform meshes are computed on a circle of radius 2 and are shown in Figure 11b-d for the wavenumbers k = 5, 10, 20. Similar to the screen problem, the uniform refinement yields a suboptimal convergence rate. However, this time the rate is closer to O(N -2/3 ), i.e., superior to that for the screen due to a less severe singularity.

In Figure 12 we display the mesh at the fourth step of the adaptive process and the real part of the corresponding BEM approximation of the boundary solution. We note, on the top of the cube, the refinement towards the ridges as well as the exploitation of the directionality. To estimate the gain due to the level of anisotropy, we have generated the equivalent isotropic mesh having the minimal sizes of the anisotropic mesh. The isotropic mesh is composed of 19 409 DOFs, i.e., 4 times the number of DOFs in the equivalent anisotropic mesh.

Spherical resonator

We consider again the case of a spherical resonator (Figure 13a). Such objects are of interest in the context of anisotropic mesh adaptation because they require the retention of some level of isotropy in the refinement in order to maintain a good approximation of the surface by planar triangles. The incident wave is a plane wave travelling in the negative x-direction with a wavenumber of k = 10. In Figure 13b, we represent the modulus of the solution computed in a 2D planar slice of the domain (the screen in the (x, y) plane with z = 0) using the boundary integral representation [START_REF] Bakry | Fiabilité et optimisation des calculs obtenus par des formulations intégrales en propagation d'ondes[END_REF]. Since this example is an isotropic example, we perform an experiment demonstrating the capabilities of our mesh adaptation technique with three configurations: with full anisotropy (similarly to what is presented in previous examples), and two additional ones with different levels of isotropy enforced. The first level is called "semi-isotropy" and the second one is a fully isotropic case. For the semi-isotropy approach, the maximum dimension of the adapted triangles is limited accordingly to the complexity. Precisely, we enforce

max(h) = C √ N ,
where h is the maximum dimension of each triangle. This restricts the elongation of the elements but still permits some elongation controlled by the parameter C. In this work, C is set to 3.5 after some numerical experimentation.

In Figure 14, we give some meshes and corresponding solutions obtained during the adaptation process for the three strategies.

In Figure 15, we report the relative L 2 errors for uniform and (anisotropic, isotropic and semi-isotropic) adapted meshes with k = 5. Similarly to the case of the planar screen, a uniform refinement yields a convergence rate of O(N -1 2 ). The anisotropic and semi-isotropic adaptive strategies initially achieve a rate of O(N -3 2 ) which is super optimal, however this rate begins to slow as higher accuracies are reached. The reason is that the shape to which the refinement is leading is not the correct one due to the elongated elements, see Figure 14. By enforcing some isotropy the convergence rate of O(N -3

2 ) is maintained up to a larger mesh. However, it still slows once the relative error falls below approximately 10 -3 . On the other hand, for the fully isotropic refinement, which restricts the height to width ratio of the triangles, at first the convergence rate is similar to the one obtained by a uniform refinement then for accuracies smaller than 10 -2 , the convergence rate is O(N -3 2 ). Furthermore, it appears that the isotropic refinement maintains this rate asymptotically since it is also refining towards the spherical shape accurately as well as resolving the edge singularity. Nevertheless, it appears that if one is interested in achieving accuracies of the order of 10 -2 or 10 -3 , the anisotropic adaptive strategy is superior.

This spherical example illustrates the limitations of the current fully anisotropic mesh strategy. The potential anisotropy, provided by the error estimate, has to comply with the surface principal curvatures. In this example, the full anisotropic error estimate implies a deformed sphere that the solver is not able to recover. The semi-anisotropic approach reduces the anisotropy to limit this artifact. For more complex geometries, it is then necessary to intersect the metric provided by the error estimate with a geometric metric representing the surface curvatures, see [START_REF] Loseille | Chapter 10 -unstructured mesh generation and adaptation[END_REF]. We use such a metric for the following example. It is clear from this example that defining a high-order adaptive mesh strategy, for instance for P 2 interpolation, is necessary to recover anisotropy while preserving the geometry.

F15 aircraft

In this test case, we illustrate the robustness of the technique with respect to the complexity of the input geometry. We consider the geometry of an unarmed F15 aircraft, depicted in Figure 16a. The computation is started with the uniform mesh of Figure 16b composed of almost 80 000 vertices and 163 000 triangles. The incident plane wave has direction d i = (1/2, 1, 0) and wavenumber k = 6. In this numerical experiment, 15 adapted meshes are generated. The imposed complexity is kept constant for 3 iterations, and then it is multiplied by a factor 2. The prescribed complexities are 10 000, 20 000, 40 000, 80 000, 160 000 . The real part of the boundary solution is represented in Figure 17. We see the high directionality of the diffracted waves with different patterns for the front and back part of the aircraft. In Table 4, we report the sizes of the generated meshes. We observe that the size is almost doubled at each sequence of fixed complexity. Note that in this case, a surface metric is intersected with Metric [START_REF] Feischl | Adaptive Boundary Element Methods: A Posteriori Error Estimators, Adaptivity, Convergence, and Implementation[END_REF] with p = 2 to guarantee a good surface resolution. The expression of the surface metric is described in [START_REF] Loseille | Chapter 10 -unstructured mesh generation and adaptation[END_REF]. The meshing CPU time is the combined time taken to evaluate the error estimate, generate the adapted mesh, and interpolate the solution from the previous mesh to the current mesh. This interpolated solution is used as a restart for the acoustic solver. The number of solver iterations ranges from 300 to more than 700. The maximal number of iterations is obtained at the 9th iteration, where new wave reflections are captured in the mesh. For the subsequent iterations, the number of solver iterations is almost constant, showing that than no more new scales in the solution appear with smaller meshes. In term of sizes, we report in Table 5 the minimal size of the anisotropic mesh. To evaluate the gain of using an anisotropic approach, an isotropic and uniform mesh have been generated with the same level of discretization as its anisotropic counterpart. The factors provided in Table 5 are the multiplication factors needed to obtain the sizes of the isotropic and uniform meshes, respectively. We observe that if an isotropic approach was used, the size of the meshes would be approximately 3 times greater. And if uniform meshes were employed, the size of the meshes to achieve the same accuracy would be between 300 and 6 million times greater! For instance, the uniform mesh having the same accuracy than the 15th mesh is composed of more than 500 million vertices. This illustrates the optimality of the anisotropic approach in term of ratio accuracy/CPU time. Figure 18 shows closer views of the mesh (at the last iteration of the adaptive process) near the singularities of the geometry of the aircraft. Most of the anisotropy is reached in these areas when the diffracted waves have low regularity. Finally in Figure 19, we display the total field in the planes y = 0 and z = 0.

Conclusion

This paper details the first extension of a metric-based anisotropic mesh adaptation strategy to the boundary element method for 3D acoustic wave propagation problems. Following the works of [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF][START_REF]Continuous mesh framework part II: validations and applications[END_REF] for 3D volume solutions, we have derived the optimal continuous mesh minimizing the P 1 interpolation error for 3D surface solutions. The method proposed is independent of the discretization technique (e.g., collocation, Galerkin) and independent of the Table 5: F15 aircraft: Expected gain in terms of DOFs of using an anisotropic approach compared to isotropic or uniform refinement. The size of the corresponding isotropic or uniform mesh has to be multiplied by the factor given in the 3rd and 4th column respectively.

complexity of the geometry. It completely remeshes at each refinement step, altering the shape, size, and orientation of elements according to the optimal metric based on a numerically recovered Hessian of the boundary solution. The resulting adaptation is truly anisotropic and we have shown via a variety of numerical examples that it recovers optimal convergence rates for domains with geometric singularities: a planar screen, a sphere with aperture, a cube with cavity, and an F15 aircraft. We note that this optimal rate of convergence cannot be achieved when using uniform refinement since the rapid variation of the boundary solution near geometric singularities requires special meshing treatment. We observed that the sizes of the generated anisotropic meshes are relatively small, with a complexity increase of between 3 and 4 times, and up to millions of times being required to achieve the same level of accuracy with isotropic and uniform refinement, respectively. For the solver, high anisotropy tends to increase the number of iterations to solve the linear system. However, the mesh adaptation framework naturally provides a good initial guess for the iterative solver by simply interpolating the solution of the previous mesh onto the current mesh. This was shown to decrease the number of iterations at very little expense. We further showed that applying a simple diagonal scaling preconditioner decreases the number of iterations significantly for anisotropic meshes.

The first encouraging anisotropic BEM results presented in this work pave the way for further developments. The first important future direction is to improve the treatment of curved surfaces. It has been observed numerically that the level of anisotropy must be controlled to ensure a good approximation of the scatterer surface. For complex geometries, a more thorough method needs to be employed where two metrics are used, one based on the recovered Hessian and the other based on the surface curvature. Another approach is to propose an adaptation strategy for P k interpolation with k ≥ 2 on curved meshes. Another important improvement concerns the FMM. The version employed in this work is not optimal with respect to computational times for non uniform meshes. It would be of interest to employ an FMM which is tailored to the adapted meshes such as the one proposed in [START_REF] Graham | Optimal panel-clustering in the presence of anisotropic mesh refinement[END_REF].

Figure 1 :

 1 Figure 1: Fast Multipole Method: (a)Decomposition of the position vector and (b) 3D cubic grid embedding the boundary.

Cell C x Figure 2 :

 x2 Figure 2: FM-BEM: Definition of the adjacent cells.

Figure 3 :

 3 Figure 3: To extend the concept of metric-based error estimates, we introduce a surface restriction Rs on the local tangent plane.

Figure 4 :

 4 Figure 4: (a) Adapted mesh on the cylinder for α = 100 and (b) comparison between a continuous and discrete evaluation of the boundary interpolation error.

  N

Figure 8 :

 8 Figure 8: Scattering by a screen: Relative L 2 errors in the scattered field computed on the circle (a) for uniform and adapted meshes for three wavenumbers: (b)k = 5, (c) k = 20 and (d) k=80.

Figure 9 :

 9 Figure 9: Screen case: final adapted meshes for k = 5 (left) and k = 20 (right).

Figure 10 :

 10 Figure 10: Screen case: evolution of the anisotropic ratio from 2nd to 5th iteration for k = 20 at side x = 0.5.

Figure 11 :Figure 12 :

 1112 Figure 11: Scattering by a sound-soft cube with cavity: Relative L 2 errors for uniform and adapted meshes for (b) k = 5, (c) k = 10 and k = 20.

Figure 13 :

 13 Figure 13: Spherical resonator: (a)Uniform mesh of the spherical resonator and screen where the volume solution is computed; (b) corresponding modulus of total field in the (x, y)-plane.

Figure 14 :

 14 Figure 14: Spherical resonator: Adapted meshes with real part of the surface solution for (a) semi-isotropic, (b) fully anisotropic (c) fully isotropic adaptation strategies, and (d) mesh at 6th level of fully isotropic adaptation.

Figure 15 :

 15 Figure 15: Scattering by a spherical resonator: Relative L 2 errors for uniform and (anisotropic, isotropic and semi-isotropic) adapted meshes with k = 5.

Figure 16 :

 16 Figure 16: F15 aircraft: (a) geometry and (b) initial uniform mesh. Iteration Minimal Size Gain over isotropic Gain over uniform 3 9 • 10 -3 3.9 378 6 4 • 10 -3 4.9 2 535 9 1.3 • 10 -3 3.0 47 572 12 5.6 • 10 -4 3.2 377 246 15 1.8 • 10 -4 3.5 6 195 211

Figure 17 :

 17 Figure 17: F15 aircraft geometry : Solution on top and back for the last iteration of the adaptive process.

Figure 18 :

 18 Figure 18: F15 aircraft geometry : Closer view of the mesh at the last iteration of the adaptive process.

1

  He 1 t e 1 He 2 t e 2 He 1 t e 2 He 2

x ŷ , and then ũ in K ref reads: ũ(x( x)) = 1 2 ( t e 1 He 1 ) x2 + ( t e 2 He 2 ) ŷ2 + 2( t e 1 He 2 ) xŷ . ũ is now linearly interpolated on K ref . Its linear interpolant Π h ũ( x) writes a x + b y + c, where the coefficients (a, b, c) ∈ R 3 satisfy the following linear system ensuring Π h

Table 3 :

 3 Spherical resonator: number of iterations (with the diagonal scaling preconditioner) with and without the use of an interpolated initial guess.

Table 4 :

 4 F15 aircraft statistics of the adapted meshes, solver iteration and CPU times.

	Iteration # vertices # triangles # Solver iterations Meshing CPU Time Solver CPU time
	3	10 562	21 120	392	42 s	19 mn 47 s
	6	15 687	31 370	521	57 s	30 mn 36 s
	9	25 428	50 852	707	1 mn 34 s	52 mn 25 s
	12	43 217	86 430	593	2 mn 44 s	65 mn 38 s
	15	80 611	161 218	364	5 mn 13 s	70 mn 53 s

Acknowledgements

This work has been funded by the Agence Nationale de la Recherche (RAFFINE project, ANR-12-MONU-0021).