Xavier Antoine

Bertrand Thierry

The

Introduction to the user guide

What is µ-diff?

The toolbox µ-diff has been developed for solving two dimensional acoustic multiple scattering problems by disks based on boundary integral equations. The toolbox is a set of Matlab functions that compute accurately and efficiently the standard integral operators with in addition pre-and post-processing facilities. When the boundary integral formulation has been written mathematically, the coding part can be easily done in µ-diff. Solving the linear system is realized within the Matlab framework. No particular computational skill is needed, except basic notions in Matlab. The boundary integral operators are spectrally discretized in Fourier basis, thanks to the circular shape of the obstacles. This leads to an accurate solution and a smaller size linear systems to solve, compared to usual Boundary Element Methods (BEMs). For high-frequency and/or a large number of scatterers, an adapted algorithm is moreover used that takes advantage of the the particular structure of the matrix (block Toeplitz).

What is not µ-diff?

Even if µ-diff is designed to solve multiple scattering problem, it is not a black-box that only solves a particular problem, i.e. µ-diff does not a priori fix the integral equation formulation to be used for solving a particular problem. This theoretical aspect is rather let to the user. However, for users who are not fluent with boundary integral equations, this user guide provides a rapid survey on boundary integral equations (see chapter 1). Furthermore, some examples of efficient solutions based on integral equations in classical cases (Dirichlet/Neumann boundary condition, penetrable scatterers) are detailled. We expect that these elements will help a nonspecialist of boundary integral equation to use µ-diff for solving multiple scattering problems.

How to use µ-diff?

It is difficult to describe the way to use µ-diff without speaking about integral equations. To simplify the presentation, we provide a basic example. Let us assume that a (circular) sound-soft scatterer Ω -is placed within a homogeneous medium which is illuminated by an incident plane wave. The mathematical problem is to compute the scattered field solution to the Helmholtz equation in the propagation domain

       (∆ + k 2 )u = 0 in R 2 \ Ω - u = -u inc
on ∂Ω - + Sommerfeld's radiation condition at infinity A possible integral representation of the scattered field u is the following

u(x) = ∂Ω - G(x, y)ρ(y)dy, ∀x ∈ R \ Ω -,
where G is the Green function and ρ is the density, solution of a boundary integral equation, which can be e.g. the EFIE (see chapter 1)

Lρ = -u inc , with Lρ(x) = ∂Ω - G(x, y)ρ(y)dy, ∀x ∈ ∂Ω -.
The problem is hence reduced to solving this boundary integral equation: knowing ρ leads to the possible computation of u at any point x of the propagation domain. Note that, in addition, the far-field can be obtained easily thanks to its integral representation. The theoretical part of the problem ends here and the computational part is handled now by µ-diff through the four following steps 1. Pre-processing (a) Creating the obstacles (b) Building the right-hand side (-u inc )

2. Assembling the matrix of the integral operator(s)

3. Solving the resulting linear system 4. Post-processing: computing the far-field, near-field (at a point of the domain or on a grid), or any other physical quantity of interest . . .

Most of the steps call some µ-diff functions except for the linear system solution which uses the already existing and optimized Matlab's functions for the direct (backslash operator \ ) and iterative solvers (GMRES, . . . ). This is summarized as a diagram in Figure 1.

To whom µ-diff is designed to?

The toolbox µ-diff is designed for any scientist who needs an easy-to-use and efficient way to either solve accurately the acoustic multiple scattering problem or to compute boundary integral operators for a collection of disks. Indeed, the µ-diff toolbox can be used

• as a solver of the multiple scattering problem,

• as a framework for more theoretical studies on boundary integral equations.

What does this user guide contain?

The first chapter recalls some well-known properties of the boundary integral equations. Even if this part is mostly theoretical, some practical aspects of boundary integral equations are also provided for the impenetrable case and some examples of robust integral equations are given for the penetrable case. Moreover, some examples of usage are provided with the µ-diff toolbox.

Chapter 2 explains the theory on which the µ-diff toolbox is based. More precisely, the boundary integral operators are projected in the Fourier bases. Their associated matrices are derived, in addition to other physical quantities of interest such as the far-or the near-fields. The µ-diff toolbox is detailed in chapter 3: each function is explained and examples are provided. Simple example are provided in chapter 4. Let us note that the list of available functions is given in the appendix A, in the alphabetic order, and arranged by folder location in appendix B.

Other tests can be launched to check that all evaluations of the integral operators works correctly

BenchmarkNeumann; BenchmarkPenetrable; BenchmarkCalderon;

5. If everything went right: congratulation! Your µ-diff installation is working! The µ-diff toolbox is based on integral equations and uses the four standard boundary integral operators. The derivation and the choice of the integral formulation is let to the user, even if some of them are given below. In this chapter, we provide all the necessary mathematical background to solve time-harmonic wave scattering problems by (penetrable or impenetrable) circular cylinders based on integral equations. The only mathematical effort requires is that the user must of course choose the integral formulation that he wants to solve: when the integral formulation is written, then it can be solved by using µ-diff. This chapter begins by presenting the potential theory and the four classical boundary integral operators with their main properties. The case of the scattering by disks is then studied and the boundary integral operators are projected in the Fourier bases, leading to infinite matrices but with some analytic expressions of their coefficients. We then discuss the finite dimensional approximation. The chapter concludes with the expressions of both the near-and far-fields, and the projections of the right-hand side onto the Fourier bases (incident wave).

Standard integral equation formulations in acoustic scattering

We present a way to derive standard direct integral equations in the case of non penetrable obstacles. Even if µ-diff can be used to solve the penetrable case, studying the impenetrable case is a suitable way to introduce the boundary integral operators and their properties. This section is strongly inspired by [START_REF] Antoine | Integral Equations and Iterative Schemes for Acoustic Scattering Problems[END_REF][START_REF] Bendali | Boundary integral equations methods in acoustics in computational acoustic scattering[END_REF][START_REF] Darbas | Préconditionneurs Analytiques de type Calderòn pour les Formulations Intégrales des Problèmes de Diffraction d[END_REF][START_REF] Thierry | Analyse et Simulations Numériques du Retournement Temporel et de la Diffraction Multiple[END_REF].

Scattering problems

Let us consider a homogeneous, isotropic and non dissipative medium filling the whole space R 2 and containing an open and bounded set Ω -, possibly not connected but such that each component is itself simply connected. Let Γ be its boundary and n the unit normal vector outwardly directed to Ω -. The domain of propagation is denoted by Ω + = R 2 \ Ω -. When illuminated by an incident time-harmonic wave u inc , the obstacle Ω -generates a scattered wave field u that is solution to the boundary-value problem

       ∆u + k 2 u = 0 in Ω + u = -u inc on Γ + u is outgoing to Ω - (1.1)
(The time-dependent form of the wave is assumed to be e -iωt , the wavenumber k := 2π/ω is supposed to be real and positive.) The first equation of system (1.1) is the well-known Helmholtz equation. For the sake of conciseness, we consider a Dirichlet boundary condition on Γ. The Neumann case will be studied later. The condition u is outgoing to Ω -means that u satisfies the Sommerfeld's radiation condition at infinity

lim x →+∞ x 1/2 ∇u • x x -iku = 0, where x = x 2 1 + x 2 2
is the euclidian norm of R 2 . Since u inc is a solution of the Helmholtz equation in R 2 , then the total field u T = u + u inc is solution of the following problem

       ∆u T + k 2 u T = 0 in Ω + u T = 0 on Γ +(u T -u inc ) is outgoing to Ω - (1.2)
It is known that, in a suitable mathematical framework, the Dirichlet scattering boundary-value problems is well-posed [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF]. More precisely, we have Theorem 1.1. The problems (1.1) and (1.2) are uniquely solvable.

Volume and boundary integral operators

Let the volume single-layer integral operator L be defined by (see e.g. [START_REF] Mclean | Strongly elliptic systems and boundary integral equations[END_REF]Theorem 6.12])

L : H -1/2 (Γ) -→ H 1 loc (R 2 ) ρ -→ L ρ, ∀x ∈ R 2 , L ρ(x) = Γ G(x, y)ρ(y) dΓ(y),
and the volume double-layer integral operator M by

M : H 1/2 (Γ) -→ H 1 loc (R 2 \ Γ) λ -→ M λ, ∀x ∈ R 2 \ Γ, M λ(x) = - Γ ∂ ny G(x, y)λ(y) dΓ(y),
where the spaces

H -1/2 (Γ), H 1 (R 2 \ Γ), H 1 loc (R 2 \ Γ) are the usual Sobolev spaces and the two-dimensional Green function G(• , •) is given by ∀x, y ∈ R 2 , x = y, G(x, y) = i 4 H (1) 0 (k x -y ). (1.
3)

The function H

(1) 0 is the first-kind Hankel function of order zero. The way the integral operators on Γ must be understood is through the duality product between the Sobolev space H 1/2 (Γ) and its dual space H -1/2 (Γ). However, when the data (u inc and Γ) are smooth enough, then the scattered field u is also regular and the dual product can be identified to the (non-hermitian) inner product on L 2 (Γ)

f, g H -1/2 ,H 1/2 = Γ f (x)g(x) dΓ(x).
This identification is considered throughout this paper. Let us now define the trace γ ± 0 and the normal derivative trace γ ± 1 operators (see [13, Appendix A]), where the plus/minus signs specify whether the trace is taken from the inside of Ω + /Ω -. First, the trace operators γ ± 0 :

H 1 (Ω ± ) → H 1/2 (Γ) are defined so that, if v ∈ C ∞ (Ω ± ), then γ ± 0 v(x) = lim z∈Ω ± →x v(z),
for almost every x ∈ Γ. By introducing the space

H 1 (Ω ± ; ∆) := {v ∈ H 1 (Ω ± ); ∆v ∈ L 2 (Ω ± )} and the linear operators γ ± * : H 1/2 (Γ) → H 1 (Ω ± ) such that γ ± 0 γ ± * ϕ = ϕ, for all ϕ ∈ H 1/2 (Γ), the normal traces γ ± 1 : H 1 (Ω ± ; ∆) → H -1/2 (Γ) can be defined [13, Equation (A.28)] ∀v ∈ H 1 (Ω ± ; ∆), ∀ϕ ∈ H 1/2 (Γ), γ ± 1 v, ϕ H -1/2 (Γ),H 1/2 (Γ) := ∓ Ω ± ∆v(x)w(x) dx + Ω ± ∇v(x) • ∇w(x) dx , (1.4)
where w := γ ± * ϕ (and thus satisfies γ ± 0 w = ϕ). Since the quantities involved in a scattering problem do not belong to H 1 (Ω + ) but rather to H 1 loc (Ω + ), the exterior trace and normal derivative trace operators are naturally extended as γ + 0 :

H 1 loc (Ω + ) → H 1/2 (Γ) and γ + 1 : H 1 loc (Ω + ; ∆) → H -1/2 (Γ) by γ + 0 (v) = γ + 0 (vv ) and γ + 1 (v) = γ + 1 (vv )
, where v is an arbitrarily compactly supported and infinitely differentiable function on Ω + which is equal to 1 in a neighborhood of Γ, and where

H 1 loc (Ω + ; ∆) := {v ∈ H 1 loc (Ω + ); ∆v ∈ L 2 loc (Ω + )}.
Remark that, when the function v is sufficiently smooth, then its normal derivative trace γ ± 1 v, given by (1.4), belongs to L 2 (Γ) and can be written as γ ± 1 v(x) = lim z∈Ω ± →x ∇v(z) • n(x), for almost every x on Γ. Note also that, the single-and double-layer potentials, introduced previously, belong not only to

H 1 loc (Ω + ) H 1 (Ω -) but also to H 1 loc (Ω + ; ∆) H 1 (Ω -; ∆) (see e.g. [13, §2.2]
). Some well-known properties of the single-and double-layer potentials are summarized in the two following propositions. Their proof can be found for example in [20, Theorems 7.5 and 9.6] for proposition 1.2 and in [START_REF] Mclean | Strongly elliptic systems and boundary integral equations[END_REF]Theorem 6.12] for proposition 1.3. Proposition 1.2. For any densities ρ ∈ H -1/2 (Γ) and λ ∈ H 1/2 (Γ), the single-layer potential L ρ and double-layer potential M λ are outgoing solutions to the Helmholtz equation in R 2 \ Γ. Moreover, the scattered field u, solution to (1.1), can be written as

∀x ∈ Ω + , u(x) = -L (∂ n u| Γ )(x) -M (u| Γ )(x).
We also have Proposition 1.3. The trace and the normal derivative trace of the operators L and M are given by the following relations

γ ± 0 L ρ = Lρ, γ ± 0 M λ = ∓ 1 2 I + M λ, γ ± 1 L ρ = ∓ 1 2 I + N ρ, γ ± 1 M λ = Dλ, (1.5)
where I is the identity operator and, for x ∈ Γ, ρ ∈ H -1/2 (Γ) and λ ∈ H 1/2 (Γ), the four boundary integral operators are defined by

L : H -1/2 (Γ) -→ H 1/2 (Γ), Lρ(x) = Γ G(x, y)ρ(y)dΓ(y), M : H 1/2 (Γ) -→ H 1/2 (Γ), M λ(x) = - Γ ∂ ny G(x, y)λ(y)dΓ(y), N : H -1/2 (Γ) -→ H -1/2 (Γ), N ρ(x) = Γ ∂ nx G(x, y)ρ(y)dΓ(y) = -M * ρ(x), D : H 1/2 (Γ) -→ H -1/2 (Γ), Dλ(x) = -∂ nx Γ ∂ ny G(x,
-∆v = k 2 v in Ω -, v = 0 (resp. ∂ n v = 0) on Γ, (1.7)
admits non-trivial solutions. Then, the operator L (resp. D) realizes an isomorphism from

H -1/2 (Γ) into H 1/2 (Γ) (resp. from H 1/2 (Γ) into H -1/2 (Γ)) if and only if k ∈ F D (Ω -) (resp. k ∈ F N (Ω -)).
These irregular frequencies k of F D (Ω -) (resp. of F N (Ω -)) are exactly the square-roots of the eigenvalues of the Laplacian operator (-∆) for the homogeneous interior Dirichlet (resp. Neumann) problem. In the multiple scattering case, that is when Ω -= M p=1 Ω - p is multiply connected, the following equalities hold true

F D (Ω -) = M p=1 F D (Ω - p ) and F N (Ω -) = M p=1 F N (Ω - p ). (1.8)
Throughout the paper, F DN (Ω -) denotes the set of all irregular frequencies

F DN (Ω -) = F D (Ω -) F N (Ω -).
(1.9)

Direct integral equations

Generalities

This section details the way of deriving direct integral equations, as described in [START_REF] Antoine | Integral Equations and Iterative Schemes for Acoustic Scattering Problems[END_REF][START_REF] Bendali | Boundary integral equations methods in acoustics in computational acoustic scattering[END_REF][START_REF] Darbas | Préconditionneurs Analytiques de type Calderòn pour les Formulations Intégrales des Problèmes de Diffraction d[END_REF][START_REF] Thierry | Analyse et Simulations Numériques du Retournement Temporel et de la Diffraction Multiple[END_REF]. This approach is nonstandard but has advantages that appear later in the user guide in section 1.2.2. The principle is to write the total field u T as a linear combination of a single-and double-layer potentials

u T (x) = L ρ(x) + M λ(x) + u inc (x), ∀x ∈ Ω + , ( 1.10)
where (λ, ρ) is now the unknown of the problem. Thanks to proposition 1.2, such an expression ensures that both u T is solution of the Helmholtz equation in Ω + and (u T -u inc ) is outgoing. Following [START_REF] Bendali | Boundary integral equations methods in acoustics in computational acoustic scattering[END_REF], an integral equation is said to be direct when the densities (λ, ρ) have a physical meaning. Indeed, for these integral equations, they are exactly the Cauchy data

(-u T | Γ , -∂ n u T | Γ )
. However, this is not a choice but a consequence of the construction of the integral equation. In electromagnetic scattering, direct and indirect integral equations are more often referred to as respectively field and source integral equations (see e.g. Harington and Mautz [START_REF] Harrington | H-field, E-field and combined field solution for conducting bodies of revolution[END_REF][START_REF] Mautz | A combined-source solution for radiation and scattering from a perfectly conducting body[END_REF] or Borel [START_REF] Borel | Résolution des équations intégrales pour la diffraction d'ondes acoustiques et électromagnétiques. Stabilisation d'algorithmes itératifs et aspects de l'analyse numérique[END_REF]). From now on, the problem, with unknown (λ, ρ), has only one equation given by the Dirichlet boundary condition on Γ. To obtain a second equation, a fictitious interior wave u - T , defined in Ω -, is introduced through

u - T (x) = L ρ(x) + M λ(x) + u inc (x), ∀x ∈ Ω -. (1.11)
Remark that, on the one hand u - T is a solution of the Helmholtz equation in Ω -and, on the other hand, due to the trace relations (1.5), the couple of unknowns (λ, ρ) satisfies the well-known jump-relations

λ = u - T | Γ -u T | Γ , ρ = ∂ n u - T | Γ -∂ n u T | Γ .
(1.12)

As the wave u - T is fictitious, it does not act on the solution u T of the scattering problem. As a consequence, the boundary condition on Γ imposed to u - T has no influence on u T . Let this constraint be represented by an operator A such that u - T is the solution of the following interior problem ∆u

- T + k 2 u - T = 0 in Ω -, Au - T = 0 on Γ. (1.13)
To build a direct integral equation, the operator A is chosen such that the field u - T vanishes in Ω -. Supposing that such an operator exists, then the following equalities hold on the boundary Γ u -

T | Γ = 0, ∂ n u - T | Γ = 0.
Consequently and thanks to the Dirichlet boundary condition u T | Γ = 0, the jump relations (1.12) read as

λ = 0, ρ = -∂ n u T | Γ ,
Therefore, both the fictitious interior field u - T and the total field u T can be composed by a single-layer potential only

u T (x) = L ρ(x) + u inc (x), ∀x ∈ Ω + , u - T (x) = L ρ(x) + u inc (x), ∀x ∈ Ω -.
The unknown ρ is finally obtained through the solution of the (direct) integral equation Au - T = 0, which can be written as AL ρ = -Au inc .

(1.14)

Both the expression and the nature of the integral equation (1.14) depend on the boundary condition imposed to u - T , represented here by the operator A. The next steps introduce three usual direct integral equations. The proofs are not provided and can be found for example in [START_REF] Bendali | Boundary integral equations methods in acoustics in computational acoustic scattering[END_REF] or [START_REF] Thierry | Analyse et Simulations Numériques du Retournement Temporel et de la Diffraction Multiple[END_REF].

EFIE (Electric Field Integral Equation)

For this first integral equation, the operator A is the interior trace operator γ - 0 on Γ. Thanks to the continuity on Γ of the single-layer integral operator L (see equation ( When k ∈ F D (Ω -), the integral operator L is no longer bijective but is still one-to-one. It can be shown that the kernel of the operator L is a subset of the kernel of the operator L . Consequently, for every solution ρ of the EFIE, the associated single-layer potential L ρ + u inc is still the solution of the scattering problem (1.2).

MFIE (Magnetic Field Integral Equation)

Another possibility is to choose A = γ - 1 , the interior normal derivative trace. Using the traces formulas (1.5), the integral equation (1.14) becomes For any irregular frequency k of F N (Ω -), the operator 1 2 I + N is no longer one-to-one. In that case and unlike the EFIE, the single-layer potential L ρ + u inc based on a solution ρ of the MFIE is not guaranteed to be the solution of the scattering problem (1.2).

1 2 I + N ρ = -∂ n u inc | Γ . ( 1

CFIE (Combined Field Integral Equation)

To avoid the irregular frequencies problem, Burton and Miller [START_REF] Burton | The application of integral equation methods to the numerical solution of some exterior boundary-value problems[END_REF] considered a linear combination of the EFIE and the MFIE by imposing a Fourier-Robin boundary condition to u - T on the boundary Γ

A = (1 -α)γ - 1 + αηγ - 0 , (1.17)
with 0 < α < 1 and (η) = 0, (1.18) where (η) is the imaginary part of the complex number η. Hence, the boundary integral equation (1.14) reads as

(1 -α) 1 2 I + N + αηL ρ = -(1 -α)∂ n u inc | Γ + αηu inc | Γ . (1.19)
This Combined Field Integral Equation (CFIE, denomination of Harrington and Mautz [START_REF] Harrington | H-field, E-field and combined field solution for conducting bodies of revolution[END_REF] in electromagnetism) or Burton-Miller integral equation [START_REF] Burton | The application of integral equation methods to the numerical solution of some exterior boundary-value problems[END_REF] is well-posed for any frequency k.

Proposition 1.7. For any k > 0 and for any complex-valued numbers α and η satisfying condition (1.18), the function L ρ + u inc is the solution of the scattering problem if and only if ρ is the solution of the CFIE (1.2).

Brakhage-Werner indirect integral equation

The indirect Brakhage-Werner Integral Equation (BWIE) is derived now. Let us first start by introducing some notations. The total field u T is here sought as a linear combination of a singleand a double-layer potentials applied to the density ψ ∈ H 1/2 (Γ)

u T = u inc + L BW ψ,
where the operator L BW with parameter η BW is given by

∀x ∈ R d \ Γ, L BW ψ(x) = (-η BW L -M )ψ(x) = Γ ∂ ny G(x, y) -η BW G(x, y) ψ(y) dΓ(y), (1.20)
with (η) = 0. The integral equation is obtained by applying the exterior trace γ + 0 to u T . Indeed, the Dirichlet boundary condition γ + 0 u T = 0 and the trace relations (1.5) directly give the Brakhage-Werner integral equation with ψ as unknown

L BW ψ = -u inc | Γ , ( 1.21)
with

L BW = -ηL -M + 1 2 I .
This second-kind integral equation does not suffer from irregular frequency [START_REF] Brakhage | Über das Dirichletsche Aussenraumproblem für die Helmholtzsche Schwingungsgleichung[END_REF].

Proposition 1.8. For all k > 0, the quantity L BW ψ + u inc is the solution to (1.2) if and only if ψ is the solution of the Brakhage-Werner integral equation (1.21).

A numerical study concerning the optimal choice of the parameter η BW (see relation (1.20)) is proposed in [START_REF] Kress | On the condition number of boundary integral operators for the exterior Dirichlet problem for the Helmholtz equation[END_REF] in the case of a single spherical or circular obstacle of radius R. For a Dirichlet boundary condition, the choice η BW = i/2 max(1/R, k) leads to a reasonable condition number of the matrix of the linear system associated to the Brakhage-Werner integral equation, for a sufficiently high frequency. Recent works have been developed on how to choose this parameter for much more general domains, see for example [12, §6] and [13, §5.1] for the case of large k and [8, §2.6 and §2.7] for the case of small frequency k. Note also that, according to [13, Remark 2.24], these results apply to both L BW and the CFIE operator, since when α = 1/2, these operators are adjoints (up to a factor of 1/2) in the real L 2 inner product. Other generalizations of these equations, when η BW is an operator, are available for example in [START_REF] Alouges | A stable well-conditioned integral equation for electromagnetism scattering[END_REF][START_REF] Antoine | Alternative integral equations for the iterative solution of acoustic scattering problems[END_REF][START_REF] Antoine | Generalized combined field integral equations for the iterative solution of the three-dimensional Helmholtz equation[END_REF]. They provide a clear background concerning the generalization and improvement of both the CFIE and BWIE.

Neumann boundary condition

Consider now the scattering of a wave by a sound-hard obstacle (Neumann boundary condition)

       ∆u + k 2 u = 0 in Ω + , ∂ n u = -∂ n u inc on Γ, + u is outgoing.
The scattered field u is sought in the form of a linear combination between a single-and a double-layer potentials

u(x) = L ρ(x) + M λ(x), x ∈ Ω + .
To build the direct integral equations, let us introduce the interior fictitious wave u - T = L ρ + M λ + u inc defined in Ω -. On Γ, a boundary condition is hence enforced to the field u - T such that it vanishes in Ω -. In that case the jump relations (1.12) lead to

λ = -u T | Γ , ρ = -∂ n u T | Γ = 0.
The Neumann boundary condition then makes the density ρ vanishing. The interior u - T and exterior u T fields are obtained through a double-layer potential

u - T (x) = M λ(x) + u inc (x) ∀x ∈ Ω -, u T = M λ(x) + u inc (x) ∀x ∈ Ω + .
As previously, the boundary condition applied to u - T is the integral equation to solve. Here are listed the four integral equations obtained in addition with their properties

Dλ = -∂ n u inc | Γ .
This first-kind integral equation is well-posed and equivalent to the scattering problem as long as k is not an irregular frequency for the interior homogeneous Neumann problem.

If k ∈ F N (Ω -) however, then after reconstruction, the solution obtained by the EFIE is correct.

• MFIE: applying a Dirichlet boundary condition to the field u -

T u - T | Γ = 0,
leads to the MFIE for the Neumann boundary condition

1 2 I + M λ = -u inc | Γ . (1.23)
This second-kind Fredholm integral equation (the operator M is compact) is well-posed.

It is equivalent to the scattering problem if k is not an interior resonance for the Dirichlet boundary-value problem. In that case, the solution obtained from the MFIE is not correct and must not be used for a practical computation.

• CFIE: applying the mixed boundary condition (1.17) to the wave u -

T (1 -α)∂ n u - T | Γ + αηu - T | Γ = 0, with 0 < α < 1 et (η) = 0,
leads to the CFIE for a Neumann boundary condition

(1 -α)D + αη 1 2 I + M λ = -(1 -α)∂ n u inc | Γ + αηu inc | Γ . (1.24)
This first-kind integral equation is uniquely solvable for any wavenumber k and is equivalent to the scattering problem.

• BWIE: the wave u T is represented as

u = (-L -ηM )ψ,
with (η) = 0. Applying the boundary condition on Γ gives the BWIE

-ηD + 1 2 I -N ψ = -∂ n u inc | Γ . (1.25)
This first-kind integral equation is also uniquely solvable for every k > 0 and equivalent to the scattering problem.

Summary

The

∈ F D (Ω -) (resp. F N (Ω -)) yes MFIE 2 nd kind k ∈ F N (Ω -) (resp. F D (Ω -)) no CFIE 2 nd (1 st ) kind k > 0 yes BWIE 2 nd (1 st ) kind k > 0 yes

Multiple scattering case 1.2.1 A more explicit writing of the integral equation formulations

The domain Ω -is now supposed to be a collection of

M disjoint bounded open sets Ω - p of R 2 , p = 1, . . . , M , such that each domain R 2 \ Ω - p is connected.
In the user guide, single-scattering designates scattering in a medium containing only one scatterer whereas multiple scattering is used for a medium containing more than one obstacle. We focus on the multiple scattering case, i.e. M ≥ 2. Since Ω -is composed of M disjoint obstacles Ω - p , p = 1, . . . , M , the single-layer volume integral operator L can be written as the sum of M operators L q , q = 1, . . . , M , defined by

L q : H -1/2 (Γ q ) -→ H 1 loc (R 2 ) ρ q -→ L q ρ q , ∀x ∈ R 2 , L q ρ q (x) = Γq G(x, y)ρ q (y) dy. (1.26)
Therefore, the single-layer potential can be decomposed as follows

∀ρ ∈ H -1/2 (Γ), L ρ = M q=1
L q ρ q , with ρ q = ρ| Γq .

(1.27)

By introducing the operators L p,q , for p, q = 1, . . . , M , defined on H -1/2 (Γ q ) by

∀ρ q ∈ H -1/2 (Γ q ), L p,q ρ q = (L q ρ q )| Γp , that is ∀ρ q ∈ H -1/2 (Γ q ), ∀x ∈ Γ p , L p,q ρ q (x) = Γq G(x, y)ρ q (y) dy, (1.28)
then the EFIE can be written in the following matrix form

      L 1,1 L 1,2 . . . L 1,M L 2,1 L 2,2 . . . L 2,M . . . . . . . . . . . . L M,1 L M,2 . . . L M,M             ρ 1 ρ 2 . . . ρ M       = -       u inc | Γ 1 u inc | Γ 2 . . . u inc | Γ M      
.

(1.29)

The same procedure can be applied to the other volume operator M and the other boundary integral operators M, N and D to obtain

L p,q : H -1/2 (Γ q ) -→ H 1/2 (Γ p ), Lρ q (x) = Γq G(x, y)ρ q (y)dΓ q (y), M p,q : H 1/2 (Γ q ) -→ H 1/2 (Γ p ), M λ q (x) = - Γq ∂ ny G(x, y)λ q (y)dΓ q (y), N p,q : H -1/2 (Γ q ) -→ H -1/2 (Γ p ), N ρ q (x) = Γq ∂ nx G(x, y)ρ q (y)dΓ q (y), D p,q : H 1/2 (Γ q ) -→ H -1/2 (Γ p ), Dλ q (x) = -∂ nx

Γq

∂ ny G(x, y)λ q (y)dΓ q (y).

Single-scattering preconditioning

Let us introduce the single-scattering operator of the EFIE L as the diagonal part of the operator L defined by

L =       L 1,1 0 . . . 0 0 L 2,2 . . . 0 . . . . . . . . . . . . 0 0 . . . L M,M       . (1.30)
Let k be a wavenumber k that is not an interior resonance (k ∈ F D (Ω -)(Ω)), implying that L and L are both invertible. The single-scattering preconditioner then simply consists in L

L -1 =       (L 1,1 ) -1 0 . . . 0 0 (L 2,2 ) -1 . . . 0 . . . . . . . . . . . . 0 0 . . . (L M,M ) -1      
.

The associated preconditioned EFIE becomes

L -1 Lρ = -L -1 u inc | Γ , (1.31)
where the operator L -1 L has the following matrix form

L -1 L =       I (L 1,1 ) -1 L 1,2 . . . (L 1,1 ) -1 L 1,M (L 2,2 ) -1 L 2,1 I . . . (L 2,2 ) -1 L 2,M . . . . . . . . . . . . (L M,M ) -1 L M,1 (L M,M ) -1 L M,2 . . . I       , (1.32)
with I the identity operator. Note that this preconditioning accelerates the convergence rate of an iterative solver, like e.g. the GMRES. The single scattering preconditioner can be applied to the other integral equations and the following result holds [START_REF] Thierry | A remark on the single scattering preconditioner applied to boundary integral equations[END_REF] Proposition 1.9. When preconditioned by their single-scattering preconditioners, the EFIE, MFIE and CFIE become identical and similar to the preconditioned BWIE (equal up to an invertible operator). In other words, after being preconditioned, the four integral equations lead to the same convergence rate when an iterative solver is applied.

This proposition shows in particular that there is no need in computing the single-scattering preconditioned versions of all the integral equations. Because the resulting operator exhibits a good convergence rate for multiple scattering, it is hard coded in µ-diff for both the Dirichlet and the Neumann cases. They are based on the single-layer potential for the Dirichlet case and the double-layer potential for the Neumann case (EFIE, MFIE or CFIE in both cases).

Mixing Dirichlet and Neumann boundary conditions

Let us assume that the collection of M obstacles is such that there is M D sound-soft obstacles and M N sound-hard scatterers, with M D + M N = M . Without loss of generality, let the sound-soft obstacles Ω - p be numbered for p = 1, . . . , M D , and the sound-hard ones for p = M D + 1, . . . , M D + M N . The problem then reads as

         (∆ + k 2 )u = 0 in Ω + , u = -u inc on Γ p , p = 1, . . . , M D , ∂ n u = -∂ n u inc on Γ p , p = M D + 1, . . . , M D + M N , + u is outgoing to Ω -.
The field is then represented as a linear combination of single-and double-layer potentials

u = M q=1 L q ρ q + M q=1 M q λ q ,
where L q is given by (1.33) and M q by

M q : H 1/2 (Γ q ) -→ H 1 loc (R 2 ) λ q -→ M q λ q , ∀x ∈ R 2 , M q λ q (x) = Γq ∂ ny G(x, y)λ q (y) dy. ( 1

.33)

A possible way to solve this problem is to apply the CFIE on both collections. To derive the CFIE, a fictitious field is introduced and the mixed condition (1.17) is applied on it. Following the theory previously described, this would lead to

λ q = 0 for q = 1, . . . , M D ρ q = 0 for q = M D + 1, . . . , M D + M N ,
and the scattered field u is represented as

u = M D q=1 L q ρ q + M D +M N q=M D +1
M q λ q .

To simplify, let us assume now that there are only two obstacles, the first one being sound-soft and the second one being sound-hard. The extension to M various scatterers is straightforward. The boundary condition (1.17) applied to the interior fictitious field leads to the integral operator

A 1 L 1 A 1 M 2 A 2 L 1 A 2 M 2 ,
where A p represents the boundary condition (1.17) on Γ p only

A p = (1 -α)γ - 1,p + αηγ - 0,p .
Applying the trace formulas leads to the system

    (1 -α) I 2 + N 1,1 + αηL 1,1 (1 -α)D 2,1 + αηM 2,1 (1 -α)N 1,2 + αηL 1,2 (1 -α)D 2,2 + αη I 2 + M 2,2     ρ 1 λ 2 = b 1 b 2 , with b p = (1 -α)u inc | Γp + αη∂ n u inc | Γp .
For M obstacles, the same procedure can be applied to obtain

I 2 + A ϕ = b, (1.34) where
• the matrix A of size M × M is such that

A p,q = (1 -α)N p,q + αηL p,q , if q ≤ M D , (1 -α)D p,q + αηM p,q , if q > M D , ( 1.35)
• the unknown vector ϕ = (ϕ q ) 1≤q≤M collects the contributions of the elementary densities

ϕ q = ρ q , if 1 ≤ q ≤ M D , λ q , if M D + 1 ≤ q ≤ N D , • the right-hand side b = (b p ) 1≤p≤M is given by b p = (1 -α)u inc | Γp + αη∂ n u inc | Γp .

The penetrable case 1.4.1 The boundary-value problem

Let assume that the obstacles are now penetrable with a contrast function n such that

n(x) = n p if x ∈ Ω - p , 0 otherwise,
where n p are constant for each obstacle. The boundary-value problem to solve is given by: find the total field u T such that

             (∆ + k 2 )u T = 0 in Ω + , (∆ + k 2 n 2 )u T = 0 in Ω -, u T = 0 on Γ, ∂ n u T = 0 on Γ, (u T -u inc ) is outgoing to Ω -.

An example of integral equation for the penetrable case

To rewrite this problem as a boundary integral equation, let the field u T be rewritten as

u T = u + + u inc in Ω + , u - in Ω -,
where u + and u -are solutions of the coupled problems

             (∆ + k 2 )u + = 0 in Ω + , (∆ + k 2 n 2 )u -= 0 in Ω -, (u + -u -)| Γ = -u inc | Γ on Γ, (∂ n u + -∂ n u -)| Γ = -∂ n u inc | Γ on Γ, u + is outgoing to Ω -.
Following the EFIE procedure, let the two fields u + and u -be expressed as a single-layer potential only

     u + (x) = L + ρ + (x) = Γ G + (x, y)ρ + (y) dy in Ω + , u -(x) = L -ρ -(x) = Γ G -(x, y)ρ -(y) dy in Ω -,
where

G ± = i 4 H (1) 0 (k ± x -y )
corresponds to the Green function with respectively the wavenumber k + := k and k -:= kn (note that k -can differ from one obstacle to another). Applying the transmission boundary conditions for the trace and normal derivative trace on Γ leads to the following couple system of boundary integral equations

  L + -L - - I 2 + N + - I 2 -N -   ρ + ρ - = -u inc | Γ -∂ n u inc | Γ . (1.36)
The upper index ± appearing on the boundary integral operators specifies which Green function is used in their expression. For instance, we have

L ± ρ = Γ G ± (x, y)ρ(y) dy.
It can be proved that this integral equation is well-posed if k is not an irregular frequency of the interior homogeneous Dirichlet problem. The example BenchmarkPenetrable.m located in Examples/Benchmark solves the penetrable scattering problem by using this integral equation. Actually, the program has been written in the dielectric case (not the acoustic one), and thus, the transmission condition on the normal derivative is there slightly modified to allow a possible jump in the normal derivative

(∂ n u + - µ µ 0 ∂ n u -)| Γ = -∂ n u inc | Γ ,
where µ 0 is the magnetic permeability of the medium and µ = µ p in Ω - p is the magnetic permeability of the domain Ω - p . The associated integral equation is then obtained by modifying the quantity

- I 2 -N - by µ - I 2 -N -
in the matrix of (1.36).

Calderón projectors

Let us consider the Helmholtz representation of the interior field u - T = u -and the exterior one

u + T = u inc -M + (u + | Γ ) -L + (∂ n u + | Γ ), u - p = u - T | Ωp = M - p (u -| Γ ) + L - p (∂ n u -| Γ ), ∀p = 1, . . . , M,
where again the ± upper index specifies the interior or exterior Green function. Note that the interior total field u - T is equal to the interior field u -. However, in the propagation domain, the total field u + T is obtained by summing the scattered wave u + and the incident wave u inc . The idea is to obtain two sets of M equations, one taken from the interior and the other one from the exterior, and next to combine them. Let us first consider the M interior problems and take the trace and normal derivative trace of the interior fields (see Eq. (1.5) for the trace formulae)

∀p = 1, . . . , M,        u - p | Γp = I 2 + M - p (u - p | Γp ) + L - p (∂ n u - p | Γp ), ∂ n u - p | Γp = D - p (u - p | Γp ) + I 2 + N - p (∂ n u - p | Γp ).
Then, by introducing the following quantities

∀p = 1, . . . , M, U - p = u - p | Γp ∂ n u - p | Γp ,
the previous set of equations can be rewritten as

∀p = 1, . . . , M, A - p,p - I 2 U - p = 0,
where

A - p,p = M - p,p L - p,p D - p,p N - p,p
.

The operator L - p,p are the operators L p,p with the interior Green functions (same for the three others)

L - p,p ρ p = Γp G -(k x -y )ρ p dΓ p | Γp .
Note that the Calderón projector P - p , defined by

P - p = I 2 + A p ,
is here hidden. This set of equations can be rewritten in the following matrix form

          - I 2 + A - 1,1 0 0 . . . 0 0 - I 2 + A - 2,2 0 . . . 0 . . . . . . . . . . . . . . . 0 0 0 . . . - I 2 + A - M,M                 U - 1 U - 2 . . . U - M       =       0 0 . . . 0      
.

(1.37) Now, a second set of equations can be obtained, thanks to the exterior part. Following the same idea, the trace and normal derivative trace of u + T on Γ p can be computed. However, the M contributions appear to be here

∀p = 1, . . . , M,          u + T | Γp = u inc | Γp + I 2 u + | Γp - M q=1 M + p,q (u + | Γq ) + L + p,q (∂ n u + | Γq ) , ∂ n u + T | Γp = ∂ n u inc | Γp + I 2 ∂ n u + | Γp -D + p,q (u + | Γq ) + N + p,q (∂ n u + | Γq ) ,
which can be rewritten as

∀p = 1, . . . , M, I 2 U + p + M q=1 A + p,q U + q = 0,
with the following quantities

U + p = u + T | Γp ∂ n u + T | Γp and A + p,q = M + p,q L + p,q D + p,q N + p,q
.

The operator L - p,q is given by (same for the three others)

L - p,q ρ q = Γq G -(k x -y )ρ q dΓ q | Γp .
Finally, this second set of equations can also be rewritten in the following matrix form

          I 2 + A + 1,1 A + 1,2 A + 1,3 . . . A + 1,M A + 2,1 I 2 + A + 2,2 A + 2,3 . . . A + 2,M . . . . . . . . . . . . . . . A + M,1 A + M,2 A + M,3 . . . I 2 + A + M,M                 U + 1 -U inc 1 U + 2 -U inc 2 . . . U + M -U inc M       =       0 0 . . . 0       . ( 1

Spectral formulation used in µ-diff

We consider that the scatterers are some circular cylinders. In this case, we can explicitly compute the boundary integral equations in some Fourier bases, leading therefore to an efficient computational spectral method when used in conjunction with numerical linear algebra methods (direct or iterative solvers).

Notations and Fourier bases

Let us consider an orthonormal system (O,

--→ Ox 1 , --→ Ox 2
). We assume that the scattering obstacle Ω -is the union of M disks Ω - p , for p = 1, . . . , M , of radius a p and center O p . We define Γ p as the boundary of Ω - p and by Γ = ∪ p=1...M Γ p the boundary of Ω -. The unit normal vector n to Ω -is outgoing. An illustration of the notations is reported in Figure 2.1. For any p = 1, . . . , M , we introduce b p as the vector between the center O p and the origin O

x b pq O p O q b q b p Ω - p Ω - q O x 1 x 2 r p (x) r q (x) α pq r(x) α q θ(x) α p θ q (x) θ p (x)
Figure 2.1: Illustration of the notations for two disks Ω - p and Ω - q and a point x ∈ Ω + . Furthermore, any point x is described by its global polar coordinates

r(x) = Ox, r(x) = r(x) , θ(x) = Angle( --→ Ox 1 , r(x)),
or by its polar coordinates in the orthonormal system associated with the obstacle Ω - p , with p = 1, . . . , M ,

r p (x) = O p x, r p (x) = r p (x) , θ p (x) = Angle( ---→ O p x 1 , r p (x)).
Let us now build a basis of L 2 (Γ) to approximate the integral operators. To this end, we first construct a basis of L 2 (Γ p ) associated with Ω - p , for p = 1, . . . , M . If the circle Γ p has a radius one and is centered at the origin, then a suitable basis of L 2 (Γ p ) is the spectral Fourier basis of functions (e imθ ) m∈Z . We adapt this basis to the general case where a p = 1 by introducing, on one hand, the functions (ϕ m ) m∈Z defined on R 2 by: ∀m ∈ Z, ∀x ∈ R 2 , ϕ m (x) = e imθ(x) , and, on the other hand, the functions (ϕ p m ) 1 p M, m∈Z given by

∀p = 1, . . . , M, ∀m ∈ Z, ∀x ∈ Γ p , ϕ p m (x) = ϕ m (r p (x)) 2πa p = e imθp(x) 2πa p .
For p = 1, . . . , M , the family (ϕ p m ) m∈Z forms an orthonormal basis of L 2 (Γ p ) for the hermitian inner product (•,

•) L 2 (Γp) ∀f, g ∈ L 2 (Γ p ), (f, g) L 2 (Γp) = Γp f (x)g(x)dΓ p (x).
To build a basis of L 2 (Γ), we introduce the functions Φ p m of L 2 (Γ) as the union of these M families

∀p, q = 1, . . . , M, ∀m ∈ Z, Φ p m | Γq = 0 if q = p, ϕ p m if q = p. The family B = {Φ p m , m ∈ Z, p = 1, .
. . , M }, also called Fourier or spectral basis, is a Hilbert basis of L 2 (Γ) for the usual scalar product (•, •) L 2 (Γ) .

Integral operators -integral equations for a cluster of circular cylinders

In view of a numerical procedure, µ-diff can use for example the weak formulation of the EFIE (see Eq. 1.29, page 22) in L 2 (Γ) based on the Fourier basis B    Find ρ ∈ H -1/2 (Γ) such that for any p = 1, . . . , M, and m ∈ Z,

(Lρ, Φ p m ) L 2 (Γ) = -u inc | Γ , Φ p m L 2 (Γ) .
Since u inc is assumed to be smooth enough (typically C ∞ ) and that Γ is C ∞ , then the scattered wavefield is also C ∞ (Ω + ) and the density ρ is (at least) in H 1/2 (Γ). Therefore, ρ can be developed in B as

ρ = M q=1 n∈Z ρ q n Φ q n
and the weak form of the EFIE is

      
Find the Fourier coefficients ρ q n ∈ C, for q = 1, . . . , M , and n ∈ Z, such that,

∀p = 1, . . . , M, ∀m ∈ Z, M q=1 n∈Z ρ q n (LΦ q n , Φ p m ) L 2 (Γ) = -u inc | Γ , Φ p m L 2 (Γ)
.

This formulation can be written under the following matrix form L ρ = U, where the infinite matrix representation L = ( L p,q ) 1 p,q M and the infinite vectors ρ = ( ρ p ) 1 p M and U = ( U p ) 1 p M are defined by blocks as

L =       L 1,1 L 1,2 . . . L 1,M L 2,1 L 2,2 . . . L 2,M . . . . . . . . . . . . L M,1 L M,2 . . . L M,M       , ρ =       ρ 1 ρ 2 . . . ρ M       , U =       U 1 U 2 . . . U M       , ( 2.1)
with, for any p, q = 1, . . . , M , and m, n

∈ Z: L p,q m,n = (LΦ q n , Φ p m ) L 2 (Γ) , ρ p m = ρ p m and U p m = -u inc | Γ , Φ p m L 2 (Γ)
. For the other integral formulations (section 1.1.3) or even for any other boundary condition, the expressions of the three boundary integral operators M , N and D are needed. Therefore, to compute an integral equation, we introduce the infinite matrices M = ( M p,q ) 1 p,q M , N = ( N p,q ) 1 p,q M and D = ( D p,q ) 1 p,q M , with the same block structure as L (see equation (2.1)). For p, q = 1, . . . , M , the coefficients of the infinite matrices M p,q , N p,q and D p,q are defined for any indices m and n in Z by

M p,q m,n = (M Φ q n , Φ p m ) L 2 (Γ) , N p,q m,n = (N Φ q n , Φ p m ) L 2 (Γ) , and D p,q m,n = (DΦ q n , Φ p m ) L 2 (Γ) .
For a numerical implementation, we can explicitly compute [START_REF] Antoine | Wide frequency band numerical approaches for multiple scattering problems by disks[END_REF][START_REF] Thierry | Analyse et Simulations Numériques du Retournement Temporel et de la Diffraction Multiple[END_REF] the matrix blocks L p,q , M p,q , N p,q and D p,q involved in L, M, N and D, for p, q = 1, . . . , M . To this end, we introduce the infinite diagonal matrices J p , (d J) p , H p and (d H) p , with general terms, for m ∈ Z, J p mm = J m (ka p ), (d J) p mm = J m (ka p ), H p mm = H (1) m (ka p ), (d H) p mm = H (1) m (ka p ).

In addition, let I p be the infinite identity matrix, and, for q = p, the infinite separation matrix S p,q between the obstacles Ω - p and Ω - q , defined by S p,q = ( S p,q m,n ) m∈Z,n∈Z and S p,q m,n = S mn (b pq ) = H

(1) m-n (kb pq )e i(m-n)α bq .

Under these notations, we rewrite the blocks L p,q , M p,q , N p,q and D p,q of the infinite matrices L, M, N and D under the matrix form, for any p, q = 1, . . . , M ,

• L p,q =        iπa p 2 J p H p , if p = q, iπ √ a p a q 2 J p ( S p,q ) T J q , if p = q,
(2.2)

• M p,q =        - 1 2 I p - iπka p 2 J p (d H) p = 1 2 I p - iπka p 2 (d J) p H p , if p = q, - ikπ √ a p a q 2 J p ( S p,q ) T (d J) q , if p = q, (2.3) • N p,q =        1 2 I p + iπka p 2 J p (d H) p = - 1 2 I p + iπka p 2 (d J) p H p , if p = q, ikπ √ a p a q 2 (d J) p ( S p,q ) T J q , if p = q, ( 2.4)
• D p,q =          iπk 2 a p 2 (d J) p (d H) p , if p = q, - ik 2 π √ a p a q 2 (d J) p ( S p,q ) T (d J) q , if p = q, ( 2.5)
where ( S p,q ) T is the transpose matrix of the separation matrix S p,q . The integral equations involve the trace or normal derivative trace of the incident wavefield on Γ. We have already introduced the infinite vector U of the coefficients of u inc | Γ in the Fourier basis. We then define similarly the infinite vector d U = (d U p ) 1 p M of the coefficients of the normal derivative trace

∂ n u inc | Γ , such that ∀p = 1, . . . , M, ∀m ∈ Z, ( d
U) p m = ∂ n u inc | Γ , Φ p m L 2 (Γ)
.

Finally, the density changes according to the integral equation and most particularly with respect to the boundary condition. To keep the same notations as previously, we introduce the densities λ and ψ (used in the BWIE) that are expanded in the Fourier basis as

Single-scattering preconditioned integral equations

The EFIE preconditioned by its single scattering component (see Section 1.2.2), given by

L -1 Lρ = L -1 u inc | Γ ,
can also be computed analytically in the Fourier bases. Indeed, let L -1 L be the matrix associated to the operator L -1 L, then

∀p, q = 1, . . . , M, ( L -1 L) p,q =        I p if p = q, a q a p ( H p ) -1 ( S p,q ) T J q otherwise. (2.6)
For the Neumann case, the preconditioned integral equation (EFIE) reads as

D -1 Dλ = D -1 ∂ n u inc | Γ .
The matrix representation D -1 D of D -1 D is then given by ∀p, q = 1, . . . , M, ( D

-1 D) p,q =        I p if p = q, - a q a p ((d H) p ) -1 ( S p,q ) T (d J) q otherwise. (2.7)
As highlighted by Proposition 1.9, there is no need to compute the preconditioned versions of the other integral equations as they lead to the same operator (up to an invertible operators, for BWIE). To solve sound-hard or sound-soft scattering, using the above integral equations appears to be a suitable choice.

Projection of the incident waves in the Fourier basis

To fully solve one of the integral equations (EFIE, MFIE, CFIE or BWIE), we need to compute the Fourier coefficients of the trace and normal derivative traces of the incident wave. We give the results for both an incident plane wave and a pointwise source term (Green's function).

For an incident plane wave, the following proposition holds [START_REF] Antoine | On the numerical approximation of high-frequency acoustic multiple scattering problems by circular cylinders[END_REF].

Proposition 2.1. Let us assume that u inc is an incident plane wave of direction β, with β = (cos(β), sin(β)) and β ∈ [0, 2π], i.e.

∀x ∈ R 2 , u inc (x) = e ikβ•x .
Then we have the following equalities

U p m = u inc | Γ , Φ p m L 2 (Γ) = d p m J m (ka p ), ( d
U) p m = ∂ n u inc | Γ , Φ p m L 2 (Γ) = kd p m J m (ka p ), with d p m = 2πa p e ikβ•bp e im(π/2-β) .
Let us consider now an incident wave emitted by a pointwise source located at s ∈ Ω + , i.e. the wave u inc is the Green's function centered at s. The Fourier coefficients of the trace and normal derivative trace of u inc on Γ are then given by the following proposition [START_REF] Thierry | Analyse et Simulations Numériques du Retournement Temporel et de la Diffraction Multiple[END_REF].

Proposition 2.2. Let s ∈ Ω + . We assume that the incident wave u inc is the Green's function centered at

s ∀x ∈ R 2 \ {s}, u inc (x) = G(x, s) = i 4 H (1) 0 (k x -s ).
The Fourier coefficients in B of the trace and normal derivative trace of the incident wave on Γ are respectively given by

U p m = u inc | Γ , Φ p m L 2 (Γ) = iπa p 2 J m (ka p )H (1) m (kr p (s)) Φ p m (s)
and

(d U) p m = ∂ n u inc | Γ , Φ p m L 2 (Γ) = k iπa p 2 J m (ka p )H (1)
m (kr p (s)) Φ p m (s).

Near-field evaluation

Outside the obstacles By using the Graf's addition theorem [START_REF] Martin | Multiple Scattering. Interaction of Time-Harmonic Waves with N Obstacles[END_REF][START_REF] Thierry | Analyse et Simulations Numériques du Retournement Temporel et de la Diffraction Multiple[END_REF], we can compute the expression of the single-and double-layer potentials at a point x located in the propagation domain Ω + .

Proposition 2.3. Let ρ ∈ L 2 (Γ) and µ ∈ H 1/2 (Γ) be two densities admitting the following decompositions in the Fourier basis

B ρ = M p=1 m∈Z ρ p m Φ p m and λ = M p=1 m∈Z λ p m Φ p m .
Then, for any point x in the domain of propagation Ω + , the single-layer potential reads

L ρ(x) = M p=1 m∈Z ρ p m L Φ p m (x) = M p=1 m∈Z ρ p m iπa p 2 J m (ka p )H (1) m (kr p (x)) Φ p m (x), (2.8)
and the double-layer potential can be expressed as

M λ(x) = M p=1 m∈Z λ p m M Φ p m (x) = - M p=1 m∈Z λ p m iπka p 2 J m (ka p )H (1) m (kr p (x)) Φ p m (x).
(2.9) Proposition 2.3 implies that, for any x in Ω + ,

u(x) = L ρ(x) + M λ(x) = M p=1 m∈Z iπa p 2 ρ p m J m (ka p ) + λ p m J m (ka p ) H (1) m (kr p (x)) Φ p m (x).

Inside the obstacles

Similarly, the potentials can be computed inside the obstacles, which is useful for penetrable obstacles for instance. In this case, only the contribution of the current obstacle is taken into account Then, for any point x inside the obstacle Ω p , the single-layer potential reads

u -(x) = L p ρ p + M p λ p , ∀x ∈ Ω p .
L ρ(x) = m∈Z ρ p m L Φ p m (x) = M p=1 m∈Z ρ p m iπa p 2 H (1) m (ka p )J m (kr p (x)) Φ p m (x), (2.10)
and the double-layer potential can be expressed as

M λ(x) = M p=1 m∈Z λ p m M Φ p m (x) = - M p=1 m∈Z λ p m iπka p 2 H (1) m (ka p )J m (kr p (x)) Φ p m (x).
(2.11)

Far-field and Radar Cross Section (RCS)

For computing the far-field pattern, let us recall that the scattered field u admits the following Helmholtz's integral representation: u = L ρ + M λ, where ρ and λ are two unknown densities.

In the polar coordinates system (r, θ) and by using an asymptotic expansion of u when r → +∞, the following relation holds [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF] ∀θ ∈ [0, 2π],

u(r, θ) = e ikr r 1/2 [a L (θ) + a M (θ)] + O 1 r 3/2 ,
where a L and a M are the radiated far-fields for the single-and double-layer potentials, respectively, defined for any angle θ of [0, 2π] by

         a L (θ) = 1 √ 8kπ e iπ/4 Γ e -ikθ•y ρ(y)dΓ(y), a M (θ) = 1 √ 8kπ e iπ/4 Γ - ik y θ • ye -ikθ•y λ(y)dΓ(y),
with θ := (cos(θ), sin(θ)). In addition, the Radar Cross Section (RCS) is defined by ∀θ ∈ [0, 2π], RCS(θ) = 10 log 10 2π |a L (θ) + a M (θ)| 2 (dB).

To optimize the far-fields computation, these relations can be written thanks to the inner product between two infinite vectors. Indeed, let us introduce a L = (( a L ) p ) 1 p M and a M = (( a M ) p ) 1 p M , where ( a L ) p and ( a M ) p are given by: ∀p = 1, . . . , M ,

         ( a L ) p = ( a L ) p m m∈Z , ( a L ) p m = ie -iπ/4 √ a p 2 √ k e -ibpk cos(θ-αp) J m (ka p )e im(θ-π/2) , ( a M ) p = ( a L ) p m m∈Z , ( a L ) p m = ie -iπ/4 ka p 2 e -ibpk cos(θ-αp) J m (ka p )e im(θ-π/2) .
(2.12) Then, we obtain the following: a L (θ) = ( a L ) T ρ and a M (θ) = ( a M ) T λ.

Finite-dimensional approximations and numerical solutions proposed in µ-diff

We now have all the ingredients to numerically solve the four integral equations EFIE, MFIE, CFIE and BWIE, for sound-soft obstacles. In fact, any integral equation for any boundary condition can be solved according to the previous developments. In practice, the resulting infinite Fourier system needs to be truncated to get a finite dimensional problem: we must pass from a sum over m ∈ Z to a finite number of Fourier modes that depends on ka p , p = 1, ..., M . Let us consider e.g. the EFIE, the extension to the other boundary integral operators being direct. The EFIE is given by equation (2.1): L ρ = -U. To truncate each Fourier series associated with (Φ p m ) m∈Z for the obstacle Ω - p , we only keep 2N p + 1 modes in such a way that the indices m of the truncated series satisfy: ∀p = 1, . . . , M , -N p m N p . The truncation parameter N p must be fixed large enough, with N p ka p , for p = 1, ..., M . An example [START_REF] Antoine | On the numerical approximation of high-frequency acoustic multiple scattering problems by circular cylinders[END_REF][START_REF] Antoine | Wide frequency band numerical approaches for multiple scattering problems by disks[END_REF] is: N p = ka p + C p , where C p weakly grows with ka p . A numerical study of the parameter N p is proposed in [START_REF] Antoine | On the numerical approximation of high-frequency acoustic multiple scattering problems by circular cylinders[END_REF][START_REF] Antoine | Wide frequency band numerical approaches for multiple scattering problems by disks[END_REF] where the following formula leads to a stable and accurate computation

N p = ka p + 1 2 √ 2 ln(2 √ 2πka p ε -1 ) 2 3 (ka p ) 1/3 + 1 , (2.13)
where ε is a small parameter (related to the relative tolerance required in the iterative Krylov subspace solver used for solving the truncated linear system (2.14), see [START_REF] Antoine | On the numerical approximation of high-frequency acoustic multiple scattering problems by circular cylinders[END_REF][START_REF] Antoine | Wide frequency band numerical approaches for multiple scattering problems by disks[END_REF]). The resulting linear system writes

Lρ = -U, (2.14)
where we introduced the block matrix L = (L p,q ) 1 p,q M , the vectors ρ = (ρ p ) 1 p M and U = (U p ) 1 p M as

L =       L 1,1 L 1,2 . . . L 1,M L 2,1 L 2,2 . . . L 2,M . . . . . . . . . . . . L M,1 L M,2 . . . L M,M       , ρ =       ρ 1 ρ 2 . . . ρ M       , U =       U 1 U 2 . . . U M       . (2.15)
For p, q = 1, . . . , M , the complex-valued matrix L p,q is of size (2N p + 1) × (2N q + 1) and its coefficients L p,q m,n are: L p,q m,n = L p,q m,n , for m = -N p , . . . , N p , n = -N q , . . . , N q . The complexvalued components of the vector ρ p = (ρ p m ) -Np m Np of size 2N p + 1 are the approximate Fourier coefficients ρ p m of ρ. For the sake of clarity, we keep on writing: ρ p m = ρ p m = ρ p m , for all m = -N p , . . . , N p . The complex-valued vector U p = (U p m ) -Np m Np is composed of the 2N p +1 Fourier coefficients of the trace of the incident wave on Γ, i.e.

U p m = U p m = u inc | Γ , Φ p m L 2 (Γ) , ∀m = -N p , . . . , N p . If N tot = M p=1 (2N p + 1)
denotes the total number of modes, the size of the complex-valued matrix L is then N tot × N tot . More generally, all the boundary integral operators can be truncated according to this process. Concerning the notations, it is sufficient to formally omit the tilde symbol ∼ over the quantities involved in sections (2.1.2) to (2.1.6). Since the four finite-dimensional matrices L, M, N and D that respectively correspond to the four boundary integral operators L, M , N and D can be computed, the linear systems that approximate the EFIE, MFIE, CFIE and BWIE can be stated. For example, the CFIE leads to (with 0 α 1 and (η) = 0)

αηL + (1 -α) I 2 + N ρ = -αηU -(1 -α)dU. (2.16)
Let us remark that the matrix obtained after discretization is always a linear combination of the four integral operators L, M, N, D and the identity matrix I. As a consequence, for a given integral equation, the resulting matrix is of size N tot × N tot and has the same block structure as e.g. L (see equation (2.15)). The finite-dimensional linear system (2.14) (or (2.16)) is accurately solved in µ-diff by using the Matlab direct solver or a preconditioned Krylov subspace linear solver that uses fast matrix-vector products based on Fast Fourier Transforms (FFTs), the choice of the linear algebra strategy (direct vs. iterative) depending on the configuration with respect to ka p and M . The use of FFTs is made possible since the off-diagonal blocks of the integral operators can be written as the products of diagonal and Toeplitz matrices [START_REF] Antoine | On the numerical approximation of high-frequency acoustic multiple scattering problems by circular cylinders[END_REF][START_REF] Antoine | Wide frequency band numerical approaches for multiple scattering problems by disks[END_REF] (see e.g. the matrices S p,q m,n in section 2.1.2). In addition, low memory is only necessary when ka p is large enough since the storage of the Toeplitz matrices can be optimized. This resulting storage technique is called sparse representation in µ-diff, in contrast with the dense (full) storage of the complex-valued matrices. Let us assume that a p ≈ a, for 1 p M . In terms of storage, the dense version of a matrix requires to store about 4M 2 [ka] 2 coefficients (assuming that N p are fixed by formula (2.13), and [r] denotes the integer part of a real number r) while the sparse storage needs about 4M 2 [ka] complex-valued coefficients. In terms of computational time for solving the linear system, the direct (multi-threaded) gaussian solver included in Matlab leads to a cost that scales with O(M 3 (ka) 3 ). For the preconditioned iterative Krylov subspace methods (i.e. restarted GMRES)), the global cost is O(M 2 ka log 2 (ka)), the converge rate depending on the physical situation and robustness of the preconditioner. From these remarks, we deduce that an iterative method can be an efficient alternative to a direct solver for large wavenumbers ka, but also for large M . We refer to [START_REF] Antoine | On the numerical approximation of high-frequency acoustic multiple scattering problems by circular cylinders[END_REF][START_REF] Antoine | Wide frequency band numerical approaches for multiple scattering problems by disks[END_REF] for a thorough computational study of the various numerical strategies. A few examples in µ-diff are provided with the toolbox. Finally, the postprocessing formulas (near-and far-fields quantities) clearly inherits of the truncation procedure (see sections 2.1.5 and 2.1.6).

Generalities

An effort has been made to keep the same notations between the mathematical framework and the µ-diff toolbox. Some explanations (goal of the function and input/output arguments) about a function name_of_the_function can be obtained by typing the classical Matlab command window help name_of_the_function;

Because µ-diff includes all the integral operators that are needed in scattering (traces and normal derivative traces of the single-and double-layer potentials), a large class of scattering problems can be solved as long as the shape of the obstacles is circular. Concerning the geometrical configurations, any deterministic or random distribution of disks is possible. Finally, µ-diff includes post-processing facilities like e.g.: surface and far-fields computations, total and scattered exterior (near-field) visualization. . . We now introduce the µ-diff toolbox based on Matlab by explaining the main predefined functions and their relations with the previous mathematical derivations. Section 3.2 gives a presentation and introduces the notations. Section 3.3 shows how to define the scattering configuration (geometry and physical parameters). Section 3.4 presents the way the integral equations must be defined and solved. Finally, section 3.5 describes the data post-processing. The µ-diff toolbox is organized following the five subdirectories (see also Figure 3.1 for a diagram of the directory arborescence)

• mudiff/PreProcessing/: pre-processing data functions (incident wave and geometry) (section 3.3).

• mudiff/IntOperators/: functions for the four basic integral operators (dense and sparse structure) used in the definition of the integral equations to solve (section 3.4).

• mudiff/PostProcessing/: post-processing functions of the solution (trace and normal derivative traces, computation of the scattered/total wavefield at some points of the spatial domain or on a grid, far-field and RCS) (section 3.5).

• mudiff/Common/: this directory includes functions that are used in µ-diff but which does not need to be known from the standard user point of view.

• mudiff/Examples/: various scripts are presented for the user in standard configurations.

In addition, the µ-diff user guide, license and credits can be found under the directory mudiff/Doc/.

Common argument and notations

A large number of arguments of the µ-diff functions are similar. For the sake of conciseness and if nothing is specified, the following arguments refer to the ones given in the Tables below, where the indices p and q vary from 1 to M (or N_scat in the µ-diff scripts). In addition, any function or value specific to µ-diff is written with the following font.

Geometry

Name Size Description

N_scat

Pre-processing

The pre-processing in µ-diff consists in defining the right-hand side (or the incident wave) and the geometry (the obstacles). The associated functions are located respectively in the folders mudiff/PreProcessing/IncidentWave and in mudiff/PreProcessing/Geometry.

Geometry: creating the obstacles

The obstacles are stored in memory as a row vector a containing the radii of the disks and a [2 × N_scat] matrix O with the centers of the disks (O(1,p=x-and O(2,p=y-coordinates of the center of the p th disk). These two arrays can be created either manually or by using some built-in functions.

Manually

The disks can be created manually by simply creating the two variables O and a containing respectively the coordinates of the disks and their radii. For example, for three obstacles placed at points (-1, 2), [START_REF] Antoine | Integral Equations and Iterative Schemes for Acoustic Scattering Problems[END_REF][START_REF] Antoine | Integral Equations and Iterative Schemes for Acoustic Scattering Problems[END_REF] where the arguments are

• bx: distance separating two centers in the x 1 -direction (be careful with the radii of the disks to avoid overlapping!).

• by: distance separating two row of obstacles in the x 2 -direction.

• Nx: number of disks in a row.

• Ny: number of rows.

For both functions, the vector of radii must be built separately and manually. For a set of unitary disks, the following command is used

a = ones(size(O,2));
The available options are (same for RectangularLattice and TriangularLattice)

• RectangularLattice(..., 'Origin', Ostart) places the first disk at Ostart position, Ostart being a [2x1] vector. (Default: [0; 0]).

• RectangularLattice(..., 'Centered', Ocenter) center the collection on the point Ocenter, Ocenter being a [2x1] vector. (Default: none but replace Ostart value).

• RectangularLattice(..., 'Direction', DIR) places the row in the increasing (DIR = +1) or decreasing x 2 -direction (DIR = -1). The default value is: DIR = +1.

Random placement

The toolbox µ-diff provides a function In that case, CreateRandomDisks builds N_scat disks with unit radius in the box. The function takes care to not overlap the disks. Note that it is possible that the function does not succeed to place the obstacles (e.g. if the user specifies too many obstacles in a box which is not large enough). This is the reason why a security test has been set: only 500 possible placements are allowed per disk.

The function comes along with a large set of optional arguments where each additional argument is optional (but the order must be kept (amin must be set, then amax, etc...!) and given by Variable Type Default Description amin scalar 1 minimal (random) radius of the obstacles allowed amax scalar 1 maximal (random) radius of the obstacles allowed dmin scalar realmin minimal distance allowed between two obstacles (not between the centers!). Setting ≤ 0 value will set dmin to realmin (i.e. ignore it) dmax scalar realmax maximal distance allowed between two obstacles (not between the centers!). The maximal distance is efficiently reached! Setting ≤ 0 value will set dmax to realmax (i.e. ignore it)

O_avoid [2 × N ] []
center of N hole(s) where the obstacles must not overlap.

Useful for example for the points source location

a_avoid [1 × N ] [] radii of the N holes dmin_avoid [1 × N ] []
minimal distance between an obstacle and a hole

The holes , represented by the * _avoid arguments, are circular regions where the obstacles must not overlap, for example where a point source is emitting a wave. Example 3: now consider that a point source is located at (2, 2) and that the obstacles must be separated from the source by at least a distance equal to 0.3. Then, the * _avoid arguments can be used and the resulting function call is

[O, a] = CreateRandomDisks(-10, 10, -10, 10, 7, 0.1, 0.5, 0.1, -1, [2;2], 0.3);

The disk centered at (2, 2) with radius 0.3 is then avoided. A second option is to set a_void to zero and set the minimal distance dmin_avoid to 0.3

[O, a] = CreateRandomDisks(-10, 10, -10, 10, 7, 0.1, 0.5, 0.1, -1, [2;2], 0, 0.3); Remark 3.1. To check if a disk is correctly placed, CreateRandomDisks calls the CheckPlacement function which can also be useful for a user who is placing the obstacles manually.

Removing disks

The function RemoveDisk aims to remove some disks of the geometrical configuration, either disk by disk, by row, by column or by radius. This can be useful for example to delete a row of disks. Here is the syntax The resulting vector is such that M_modes(p)=N p , where N p satisfies (2.16), with obviously a minimal value equal to 0 (which consists in only one mode). If k is a vector (which corresponds to one wavenumber per obstacle), then FourierTruncation uses formula (2.16) with k(p) as the wavenumber. The following options are moreover available

• M_modes = FourierTruncation(..., 'Min', MIN);
To force a minimal value: M_modes(p) is then either the min value between MIN and formula (2.16).

• M_modes = FourierTruncation(..., 'Tol', TOL);

The tolerance, set by default to 10 -10 , is then set to TOL.

Incident waves

Generalities

Two different incident waves are available in the µ-diff toolbox: the plane wave and the point source wave. The user can build his own incident wave. They are all located in the directory PreProcessing/IncidentWave/.

As explained in section 2.1.4, a right-hand side b is decomposed by blocks, each of these blocks representing one obstacle: b = (b p ) p=1,...,M . A different condition can be applied on two different obstacles (e.g. Dirichlet on Ω 1 and Neumann on Ω 2 ) or a different integral equation can be considered on each obstacle (e.g. EFIE on Ω 1 and MFIE on Ω 2 ). To this end, µ-diff builds each block separately thanks to the function BlockIncidentWave which computes the vector b p . According to the input data, the function builds one of the available right-hand sides described in Table 3.1.

On the other hand, the common function IncidentWave computes the whole vector b. For all obstacles p, the function BlockIncidentWave is called and the whole vector is assembled.

In addition, for all the incident waves, an interface function is available. These easy-to-use interfaces build the whole vector on only one pattern (trace of plane wave, normal derivative of a point source wave, . . . ). They are located in the interface/ directory and their names allow an easy interpretation (see also table 3.1, column µ-diff name): PlaneWave, PointSource, DnPlaneWave,. . . . Let us recall that the help informations of the interface functions contains the mathematical description of the incident wave. A plane wave is defined by e ik(cos(β)x 1 +sin(β)x 2 2 DnPlaneWave beta_inc normal derivative of a plane wave of angle of incidence DnPlaneWavePrecond beta_inc same as DnPlaneWave but multiplied by the inverse of the double-layer block diagonal operator Table 3.1: Right-hand sides already coded in µ-diff.

The two main functions, BlockIncidentWave and IncidentWave, are now detailed.

BlockIncidentWave

This function computes the block vector of the opposite of the coefficients of an incident wave, either the trace of the normal derivative trace, on one of the obstacles, in the Fourier bases. Its syntax is Bp = BlockIncidentWave(Op, ap, Np, k, TypeOfWave, Param);

where TypeOfWave is a scalar value specifying the incident wave (see table 3.1) and Param is the parameter of the wave: angle of direction, position of a point source, . . . . The returned value Bp is a column vector of length 2Np + 1.

IncidentWave

The function call is the following

B = IncidentWave(O, a, M_modes, k, TypeOfWave, Param)
The resulting vector B is of size M p=1 (2M_modes(p) + 1). The value Param is the same as for BlockIncidentWave whereas TypeOfWave is a vector of size M , where TypeOfWave(p) is the fixed choice for the block b p . In other word, the block b p is built by calling the function with the arguments BlockIncidentWave(Op, ap, Np, k, TypeOfWave, Param);. To simplify, if TypeOfWave is a scalar value, then it is considered as a vector with the same scalar value.

Example 1: Building a vector associated to the trace of an incident plane wave of direction beta_inc is done thanks to the following command (the 1 argument refers to as PlaneWave) In other words, this builds the vector ( [START_REF] Antoine | Alternative integral equations for the iterative solution of acoustic scattering problems[END_REF][START_REF] Antoine | Generalized combined field integral equations for the iterative solution of the three-dimensional Helmholtz equation[END_REF] can be translated as PointSource, DnPointSource)

b = -u inc | Γ 1 -∂ n u inc | Γ 2 ,
where

u inc = i 4 H (1) 0 (k x -x s ), with x s = [1, 2].
Remark 3.2. Remember that the resulting vector corresponds to the opposite of the trace or normal derivative trace!

Integral operators

Generalities

The functions defining the integral operators are available in the directory IntOperators/ which has the Dense/ and Sparse/ subdirectories for the dense (matrix) and sparse (@function) representations of the four basic integral operators used in scattering, i.e. L, M, N and D, given in their infinite dimensional operator versions by respectively (2.2), (2.3), (2.4) and (2.5). Preconditioned versions of the operators by their single-scattering operators [START_REF] Thierry | A remark on the single scattering preconditioner applied to boundary integral equations[END_REF] are also defined. Following Proposition 1.9, only L -1 L (EFIE Dirichlet) and D -1 D (EFIE Neumann) are provided since they are the only ones needed for the Dirichlet and Neumann problems. Two different types of storage can be used within the µ-diff toolbox: dense and sparse. The assembly of the matrix is almost the same in both cases. The way µ-diff is developed is close to the mathematics for the assembly process, in the sense that the matrix A of an integral operator is built block-by-block (A p,q ). For both storages, two main functions exist: a block function (BlockIntegralOperator and SpBlockIntegralOperator) and a global matrix function (IntegralOperator and SpIntegralOperator) which assemble all the blocks into a matrix. This separation allows the user to either build a simple matrix for one operator or to construct a more complex matrix where each block represents a different operator or a linear combination of them. Let us also note that all the operators have interface functions (located in Dense/Interface or Sparse/Interface folders). For example, SingleLayer is the interface function of IntegralOperator to build the single-layer operator, and BlockSingleLayer the one for BlockIntegralOperator. The same applies for the sparse storage with a prefix Sp: SpSingleLayer, SpBlockSingleLayer. The dense storage should be preferred for solving small scale problems (most particularly low frequency problems), when the memory storage and the CPU cost is not a problem, where the sparse storage must be used when the limits are reached. For the sparse storage, which is detailed later, the block matrices are stored through vectors and the matrix-vector product is then fast when an iterative Krylov subspace solver is used. This is essentially a suitable strategy when a sufficiently high number of modes N p is required per obstacle, i.e. for large enough wave numbers k. The drawback is that some instabilities may arise in the numerical process if the truncations of the Fourier series are not done correctly. The formula (2.16) provides a stable result following [START_REF] Antoine | On the numerical approximation of high-frequency acoustic multiple scattering problems by circular cylinders[END_REF].

Available integral operators and numbering

As for the incident wave, for each operator, there exists a function that builds the whole operator for the multiple scattering problem, both for the dense and sparse storages. The available operators are listed in Table 3.2, with their unique identifier (integer), their associated name in µ-diff (useful to get their interface functions) and their definitions.

Dense storage

Building a block A p,q

An elementary block matrix in a global matrix can be created by using the following function The Weight argument is optional and set to 1 by default. The quantity TypeOfOperator specifies the integral operator to compute, thanks to the numbering of Table 3

.5 × I 1,1 + M 1,1 L 1,2 D 2,1 0.5 × I 2,2 + N 2,2
can be computed by using a three-dimensional array and the null operator

Sparse storage

Storing the matrices in a sparse way leads to a significant reduction in memory storage and the ability to get access to fast matrix-vector product evaluations. The linear system must then be solved by an iterative solver since the global matrix is no longer built. Combining linearly the matrices is still possible but this is realized during the computation of the matrix-vector products. Indeed, summing two matrices for two different integral operators is not guaranteed to keep the particular matrix structure. If the problem involves (at least) two different integral operators, they must hence be computed separately.

Compressed storage

To explain how the matrix is stored by using the µ-diff sparse storage, let us consider that A is a matrix corresponding to one of the four integral operators L, M, N or D. The matrix A has the following special structure (for p, q = 1, . . . , M and p = q)

• A p,p is diagonal,

• A p,q is full and can be written as A p,q = A p,q L T p,q A p,q R , where A p,q L and A p,q R are diagonal and called respectively the left and right parts, and T p,q = (T p,q m,n )), with

T p,q m,n = iπe i(n-m)αpq H (1)
n-m (kb pq ), is a Toeplitz matrix: T p,q m,n = T p,q m+1,n+1 . The idea is that the diagonal matrices can be stored as vectors containing the diagonal elements. A matrix-vector product between a diagonal matrix and a vector is then simply an elementby-element multiplication. The Toeplitz matrix can also be stored in a compressed form as a vector and the matrix-vector product is handled by a cross-correlation based on the Fast Fourier Transform (FFT). The diagonals parts are stored in the three-dimensional arrays A L and A R ( L for Left part and R for Right part). The left part A L also includes the diagonal part of the matrix

A L (:, p, q) = diag(A p,p ) if p = q, diag(A p,q L ) otherwise, A R (:, p, q) = 0 if p = q, diag(A p,q R ) otherwise.
The matrix T p,q of size (2N p + 1) × (2N q + 1) is then also compressed as a Toeplitz matrix following Clearly, the root vector containing the first row and the first column is enough to rebuild the matrix. This root vector, called A p,q M , is of size (2N p + 1) + (2N q + 1) -1 = 2N p + 2N q + 1 and such that

T p,q =       
A p,q M =                     t (2Nq+1)+(2Np+1)-1 t (2Nq+1)+(2Np+1)-2
. . .

t (2Nq+1)+2 t (2Nq+1)+1 t 1 t 2 t 3 . . . t (2Nq+1),1                     .
Finally, all these root vectors are stored in the three-dimensional array A M A M (:, p, q) = 0 if p = q, A p,q M otherwise.

The global matrix A is thus stored through three different parts: A L , A M and A R . From the computer point of view, as a three-dimensional array has a fixed size in each direction, the second and third dimensions of the arrays are both of size N_scat and the length in the first dimension is given by 3D-array Length in the dimension. . . 1 2 3

A L max p (2N p + 1) N_scat N_scat A M max p max q =p [(2N p + 1) + (2N q + 1) -1] N_scat N_scat A R max p (2N p + 1)
N_scat N_scat

With this constraint, the vector A L (:, p, q) can be larger than 2N p + 1 and thus A p,q L must be extracted: A p,q L = A L (1 : 2N p + 1, p, p). The same occurs for A M and A R . Finally, these three-dimensional arrays are merged into a cell A, representing the matrix A, such that, by using the Matlab notations

A{1} = A L , A{2} = A M , A{3} = A R .

Fast matrix-vector products

The matrix-vector product between A and a vector X is divided into different elementary operations. Let us consider Y = AX and more particularly the p th component

Y p = q A p,q X q .
By using the previous notations, we have

Y p = A p,p L X p + q =p
A p,q L (A p,q M (A p,q R X q )).

Since A p,p L , A p,q L and A p,q R are diagonal matrices (that are stored as vectors), the corresponding matrix-vector products are easy to compute. The only difficulty concerns A p,q M Z q (where Z q = A p,q R X q ). This can however be achieved efficiently. Indeed, the discrete cross-correlation product (xcorr) between A p,q M and Z q gives W q = xcorr(Z q , A p,q M ), where the bar denotes the complex conjugate of a complex number. The result W q is then extracted from W q by W q = W q (2N q + 1 : 2N q + 1 + 2N p ).

The matrix-vector product between A and X is then done in a fast way (cross-correlation is efficiently evaluated through the FFT).

Assembling the matrix

The assembly process is very similar to the dense one, the difference being the prefix Sp. The block function is then called by SpBlockIntegralOperator(Op, ap, Np, Oq, aq, Nq, Nmax, k, TypeOfOperator, Weight);

The function returns three vectors corresponding to A p,q L , A p,q M and A p,q R . The major difference here is that no linear combination of operators can be done during the assembly process for the sparse representation. The quantity TypeOfOperator cannot hence be a vector but must be a scalar. The assembly process of the global matrix can be done through A = SpIntegralOperator(O, a, M_modes, k, TypeOfOperator, Weight);

The result is a Matlab cell of three components, where each component is a three-dimensional array. The quantity TypeOfOperator can be a scalar or a matrix (not a three-dimensional array nor a vector!). If TypeOfOperator is a matrix, then the block A p,q is assumed to be of type TypeOfOperator(p,q). The linear combination of operators can still be done, but must be specified in the matrix-vector product (see below). Note that a function exists to add the (weighted-)identity to a sparse operator (see below).

Radar Cross Section (RCS)

The Radar Cross Section (RCS) can be computed either from the far-field with FarField_to_RCS or directly from the density with RCS. In µ-diff, the RCS σ is obtained by (F being the far-field) where the arguments are exactly the same as for the far-field. In fact, RCS first calls FarField and then FarField_to_RCS.

σ = 10 log 10 (2π|F | 2 ).
Remark 3.4. The RCS function is also interfaced with the two following ready-to-use functions RCSSingleLayer and RCSDoubleLayer. Both of them are located in the directory called PostProcessing/FarField/Interface.

Near-field: inside and outside the obstacles

The near-field of a potential is given by equations (2.8) and (2.9) outside the obstacles and by (2.10) and (2.11) inside. In the same way, µ-diff separates the outside and the inside computations. Let us first explain how to compute the field outside the obstacles. The user only needs to create a mesh by using e.g. meshgrid and launch ExternalPotential through the syntax U = ExternalPotential(X, Y, O, a, M_modes, k, Density, Weight, OPTIONS);

The resulting matrix U has zero value inside the obstacle. By default, the computation is not realized on the boundary of the obstacles (this can be changed thanks to the options). The arguments are

• Density and Weight are exactly the same as for FarField,

• the quantities X and Y are the grid points which are created by the Matlab meshgrid function,

• the OPTIONS must be chosen from the options of the potential, see below (after the interior potential).

Let us consider now the computation inside the obstacle Ω p . For the interior field, the only contribution is assumed to be of the form

L p ρ p + M p λ p .
In other words, the contributions of the other obstacles are not taken into account. In the same way as for the external potential, the computation is realized in µ-diff by the function call

The Dirichlet boundary-value problem

Let us consider the scattering problem by a collection of sound-soft obstacles

    
(∆ + k 2 )u = 0, in Ω + , u = -u inc , on Γ, u outgoing,

The Neumann boundary-value problem

Let us now consider the sound-hard scattering problem The post-processing is based on the double-layer potential (compared to Dirichlet, the modification is realized in the RCS function, the last argument [START_REF] Alouges | A stable well-conditioned integral equation for electromagnetism scattering[END_REF]0] is then replaced by

Figure 1 :

 1 Figure 1: Schematic structure of a µ-diff script.

  History of convergence of the GMRES

Figure 2 : 1

 21 Figure 2: 4 of the different figures that pop up after launching BenchmarkDirichlet;. The first figure shows the two obstacles and the second one the total field. The two others figures, (c) and (d) represent respectively the radar cross section for the different boundary integral equations and the history of convergence of the GMRES solver.

1. 1 . 4 21 1. 2 22 1. 3 23 1. 4

 14212223234 Brakhage-Werner indirect integral equation . . . . . . . . . . . . . . . . 19 1.1.5 Neumann boundary condition . . . . . . . . . . . . . . . . . . . . . . . . 20 1.1.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Multiple scattering case . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.2.1 A more explicit writing of the integral equation formulations . . . . . . 21 1.2.2 Single-scattering preconditioning . . . . . . . . . . . . . . . . . . . . . . Mixing Dirichlet and Neumann boundary conditions . . . . . . . . . The penetrable case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 1.4.1 The boundary-value problem . . . . . . . . . . . . . . . . . . . . . . . . 24 1.4.2 An example of integral equation for the penetrable case . . . . . . . . . 25 1.4.3 Calderón projectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

Theorem 1 . 4 .

 14 y)λ(y)dΓ(y).(1.6) All along the user guide, the boundary integral operators are written with a roman letter (e.g. L) whereas the volume integral operators are denoted by a calligraphic letter (e.g. L ). According to [21, Theorems 3.4.1 and 3.4.2], the boundary integral operators L and D are invertible, providing k is not an irregular frequency. More precisely, we have Let F D (Ω -) (resp. F N (Ω -)) be the countable set of positive wavenumbers k accumulating at infinity such that the interior homogeneous Dirichlet (resp. Neumann) problem

•

  EFIE: applying a homogeneous Neumann boundary condition to the wave u - T ∂ n u - T | Γ = 0, (1.22) leads to the EFIE for the Neumann boundary condition

  Finally, we set: λ = ( λ p ) 1 p M and Ψ = ( Ψ p ) 1 p M , where each block λ p = ( λ p m ) m∈Z and Ψ p = ( Ψ p m ) m∈Z is defined by: ∀m ∈ Z, λ p m = λ p m and Ψ p m = ψ p m .

Proposition 2 . 4 .

 24 Let ρ ∈ L 2 (Γ) and µ ∈ H 1/2 (Γ) be two densities admitting the following decompositions in the Fourier basis B

3 . 4 3 . 6

 3436 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.2 Common argument and notations . . . . . . . . . . . . . . . . . . . . 40 3.3 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.3.1 Geometry: creating the obstacles . . . . . . . . . . . . . . . . . . . . . . 42 3.3.2 Truncation of the Fourier series . . . . . . . . . . . . . . . . . . . . . . . 46 3.3.3 Incident waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 Integral operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.4.1 Generalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.4.2 Available integral operators and numbering . . . . . . . . . . . . . . . . 50 3.4.3 Dense storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.4.4 Sparse storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.5 Post-Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Examples available in µ-diff . . . . . . . . . . . . . . . . . . . . . . . . 62

[ 1 ×[ 2 ×[ 1 ×[ 1 × 1 ][ 1 × 1 ]Figure 3 . 1 :[ 1 × 1 ] 4 H ( 1 )0∂ 6 L- 1 L 7 D- 1 D

 12111113111416171 Figure 3.1: Arborescence of µ-diff toolbox. 41

  and (-15, 10) with respective radii 0.1, 0.5 and 10, one gets O = [-1, 5, 2 ; -15, 5, 10]; a = [0.1, 0.5, 10];Periodic placementTwo built-in functions are available with the toolbox to periodically place some disks, with a rectangular or a triangular lattice, as shown on Figure3.2. The two functions must be called as follows, first for the rectangular lattice O = RectangularLattice(bx, by, Nx, Ny, OPTIONS); and second for the triangular lattice O = TriangularLattice(bx, by, Nx, Ny, OPTIONS);

Figure 3 . 2 :

 32 Figure 3.2: Rectangular (left) and triangular (right) lattice with bx = 3, by = 4, Nx = 5, Ny=6.

[

  O, a] = CreateRandomDisks(xmin, xmax, ymin, ymax, N_scat, amin, amax, dmim, dmax, O_avoid, a_avoid, dmin_avoid, dmax_avoid);

Example 1 :Example 2 :

 12 creating N_scat random disks with random radii is realized by the command [O, a] = CreateRandomDisks(xmin, xmax, ymin, ymax, N_scat, amin, amax); building 7 obstacles in the box [-10, 10] × [-10, 10] with radii between 0.1 and 0.5. The disks must be separated by a minimal distance equal to 0.1 and without maximal value. The command is then [O, a] = CreateRandomDisks(-10, 10, -10, 10, 7, 0.1, 0.5, 0.1, -1);

Figure 3 . 3 :

 33 Figure 3.3: Periodic placement with a row and a column deleted.

Value µ-diff function name Param Description 1

 1 PlaneWave beta_inc trace of a plane wave of angle of direction beta_inc.

4 H

 4 wave emitted by a point source placed at XS. Such a wave is defined in µ-diff by i of the wave emitted by a point source placed at XS 5 PlaneWavePrecond beta_inc same as PlaneWave but multiplied by the inverse of the single-layer block diagonal operator (see section 2.1.3 on the single-scattering preconditioner) 6

BExample 2 :

 2 = IncidentWave(O, a, M_modes, k, 1, beta_inc); or by using the interface function B = PlaneWave(O, a, M_modes, k, beta_inc); For two obstacles and a point source centered at (1, 2), if b 1 is the trace of the wave and b 2 is the normal derivative trace, then B = IncidentWave(O, a, M_modes, k, [3;4],[1;2]);

7

 7 PrecondNeumanndouble-layer preconditioned by its diagonal

  BlockIntegralOperatorApq = BlockIntegralOperator(Op, ap, Np, Oq, aq, Nq, k, TypeOfOperator, Weight);

1 :Example 2 :Example 3 :

 123 .2. If TypeOfOperator Example The single-layer potential L is L = IntegralOperator(O, a, M_modes, k, 2); The MFIE operator 0.5 × I + N for a Dirichlet boundary condition (1.16) is A_MFIE = IntegralOperator(O, a, M_modes, k, [1,4], [0.5,1]); The following two-obstacles operator is 0

Example 4 :

 4 TypeOfOp = zeros(2,2,2); Weight = zeros(2,2,2); TypeOfOp(:, 1,1) = [1, 3]; Weight(:, 1,1) = [0.5, 1]; %block A_{1,1} TypeOfOp(:, 1,2) = [2, 0]; Weight(:, 1,2) = [1, 0]; %block A_{1,2} TypeOfOp(:, 2,1) = [5, 0]; Weight(:, 2,1) = [1, 0]; %block A_{2,1} TypeOfOp(:, 2,2) = [1, 4]; Weight(:, 2,2) = [0.5, 1]; %block A_{2,2} A = IntegralOperator(O, a, M_modes, k, TypeOfOp, Weight); The Brakhage-Werner Integral Equation (BWIE) for a Dirichlet problem (1.21) is solved by writing k = 1; beta_inc = pi; eta = i/k; Uinc = PlaneWave(O, a, M_modes, k, beta_inc); A = IntegralOperator(O, a, M_modes, k, [1, 2, 3], [0.5, -eta_BW, -1]); psi = A \ Uinc;

Example 1 :Example 2 :Example 1 :

 121 Creation of µ-diff sparse matrix of the single-layer operator L L = SpIntegralOperator(O, a, M_modes, k, 2);L is now a cell with three components! Note the interface function can also be used L = SpSingleLayer(O, a, M_modes, k); For two obstacles, building the following matrixM 1,1 0.5 × L 1,2 1.5 × D 2,1 2.5 × N 2,2can be done with the commands TypeOfOp = zeros(2,2); Weight = zeros(2,2); TypeOfOp(1,1) = [3]; Weight(1,1) = [1]; %block A_{1,1} TypeOfOp(1,2) = [2]; Weight(1,2) = [0.5]; %block A_{1,2} TypeOfOp(2,1) = [5]; Weight(2,1) = [1.5]; %block A_{2,1} TypeOfOp(2,2) = [4]; Weight(2,2) = [2.5]; %block A_{2,2} A = SpIntegralOperator(O, a, M_modes, k, TypeOfOp, Weight); Assembling: adding the identity It is possible to add the identity (multiplied by a constant) to a sparse operator A = SpAddIdentity(A, alpha, M_modes); which simply returns A = A + αI. If L is the sparse representation of the single-layer operator, then L = SpIntegralOperator(O, a, M_modes, k, 2); alpha_L = SpAddIdentity(L, 0.5, M_modes); computes 0.5 × L. and Weight(p,2)=γ p . If Weight is a row of length 2, then Weight(1)=η p and Weight(2)=γ p for all p (the coefficient are the same for all obstacles). With these notations, the resulting far-field is N_scat p=1 Weight(p, 1)L Density(:,1) + Weight(p,2)M Density(:,2).

Example 1 :Example 2 :Example 3 :Remark 3 . 3 .

 12333 the Dirichlet EFIE (1.15): the scattered field u reads as u = L ρ and the far-field is then computed by F = FarField(O, a, M_modes, k, theta, rho, [1,0]); the Dirichlet BWIE (1.21): the scattered field u is given by u = (-ηL -M )ψ and the corresponding far-field by eta = i/k; F = FarField(O, a, M_modes, k, theta, psi, [-eta,-1]); if the densities ρ and λ are known and if the user needs to compute for example M p=1 pL p ρ p -pM p λ p , then this can be done by the function call F = FarField(O, a, M_modes, k, theta, [rho, lambda], [[1:p].', -[1:p].']); Let us remark that the FarField function is also interfaced with the ready-to-use function: FarFieldSingleLayer and FarFieldDoubleLayer. These two functions are located in the PostProcessing/FarField/Interface directory.

  The computation is realized by the function callR = FarField_to_RCS(F);where F has been computed by the FarField function. Otherwise, the RCS can be computed directly from the densities by R = RCS(O, a, M_modes, k, theta, Density, Weight);

0

  History of convergence of the GMRES, k = 1, without restart Iteration number GMRES residual (log) History of convergence of the GMRES (d) Absolute value of the scattered field (e) Absolute value of the total field (f) Real part of the total field

Figure 4 . 1 :

 41 Figure 4.1: Different results for a Dirichlet problem and an incident plane wave solved using µ-diff. The first figure shows the three obstacles and the second one the radar cross section obtained for the different integral equation (all superimposed). The next figure (c) represents the GMRES history of convergence of the different integral equations. The figure (d) shows the absolute value of the scattered field, and the last two, (e) and (f), represent respectively the absolute value and the real part of the total field.

0

  History of convergence of the GMRES, k = 1, without restart Iteration number GMRES residual (log) History of convergence of the GMRES (d) Absolute value of the scattered field (e) Absolute value of the total field (f) Real part of the total field

Figure 4 . 2 :

 42 Figure 4.2: Different results for a Dirichlet problem for a wave emitted by a point source solved using µ-diff. The first two pictures show respectively the obstacles and the radar cross section, while the third one shows the GMRES history of convergence. The absolute value of the scattered field is presented on subfigure (d), and the total field on figures (e) (absolute value) and (f) (real part).

  k 2 )u = 0, in Ω + , ∂ n u = -∂ n u inc , on Γ, u outgoing.An efficient solution to this problem is given for example by a preconditioned integral equation for sound-hard obstacles. Here, we only present this solution but the extension to other kinds of integral equations is direct. The µ-diff script is close to the one developed for the Dirichlet problem, only the two following functions must be modified: PrecondDirichlet is replaced by PrecondNeumann and the right-hand side PlaneWavePrecond is now given by DnPlaneWavePrecond. For the Neumann problem, the preconditioned boundary integral equation is based on the double-layer representationu = M λ, D -1 Dλ = -D -1 ∂ n u inc . % Three unit disks O = [-5, 0, 5; -2, 0, 2]; a = [1, 1, 1]; %Set the parameters... k = 1; %wavenumber beta_inc = 0; %incident angle %Fourier series truncation parameter M_modes = FourierTruncation(a, k, 'Min', 1); %Right-hand side DnUincPrecond = DnPlaneWavePrecond(O, a, M_modes, k, beta_inc); %Matrix of the system (the two following lines are the same) APrecond = PrecondNeumann(O, a, M_modes, k); %Solving (here, direct) lambda = APrecond \ DnUincPrecond;

* 2 *

 2 pi/360; %Radar Cross Section associated with the double-layer potential myRCS = RCS(O, a, M_modes, k, theta_RCS_rad, lambda, [0,1]); plot(theta_RCS, myRCS, 'k'); With the previous parameter, we obtain the result shown on figure 4.3 for an incident plane wave of direction 0 and on figure 4.4 for a point source centered on [-10, 0]. Only the case of the preconditioned integral equation has been shown.

0

  History of convergence of the GMRES, k = 1, without restart Iteration number GMRES residual (log) Precond (c) History of convergence of the GMRES (d) Absolute value of the scattered field (e) Absolute value of the total field (f) Real part of the total field

Figure 4 . 4 :•

 44 Figure 4.4: Neumann boundary value problem with a point source emitter solved using the single scattering preconditioned integral equation. In the respective order is shown: the obstacles, the radar cross section, the history of convergence of the GMRES and then, the absolute value of the scattered field, the absolute and real part of the total field.

  , which is well-conditioned but can provide inaccurate solution. These two integral equations are solved in the example function BenchmarkCalderon. The matrix A ±

	Chapter 2
	Multiple scattering by disks:
	approximation method in µ-diff
	.38)
	To obtain only one set of equations, systems (1.37) and (1.38) must be combined. A first-kind
	integral equation is obtained by doing (1.37) + (1.38) , which is known to have a poor condition
	number but is robust, and a second-kind integral equation can be obtained by computing
	(1.37) -(1.38) p,q
	CalderonProjector.

is computed by the BlockCalderonProjector function and the whole matrix (A ± p,q ) p,q by Contents 2.1 Spectral formulation used in µ-diff . . . . . . . . . . . . . . . . . . . . 29 2.1.1 Notations and Fourier bases . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.1.2 Integral operators -integral equations for a cluster of circular cylinders 31 2.1.3 Single-scattering preconditioned integral equations . . . . . . . . . . . . 33 2.1.4 Projection of the incident waves in the Fourier basis . . . . . . . . . . . 33 2.1.5 Near-field evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.1.6 Far-field and Radar Cross Section (RCS) . . . . . . . . . . . . . . . . . 35 2.2 Finite-dimensional approximations and numerical solutions proposed in µ-diff . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

Table 3 .

 3 2: Available integral operators in µ-diff, their (unique) identifier, function name (interface) and the mathematical definition. The zero operator function does not have an interface function and the sparse version is obtained by adding the prefix Sp (SpSingleLayer, SpDnDoubleLayer,. . . ).The block interface functions are also prefixed by Block (SpBlockSingleLayer, BlockSingleLayer,. . . ).

2 The Neumann boundary-value problem . . . . . . . . . . . . . . . . . 4.3 Mixing Dirichlet and Neumann boundary conditions . . . . . . . . . 4.4 Penetrable case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  where U has zero value outside the obstacles and by default on the boundaries too. The arguments are exactly the same as ExternalPotential. Finally, the common options for ExternalPotential and InternalPotential are described below. For the sake of clarity, the examples are shown with only ExternalPotential even if they also apply for InternalPotential Some examples are available in the Examples/ directory and especially in Examples/Benchmark, where each file is independent and solves a different problem (Dirichlet, Neumann or penetrable). They should also be launched to check that µ-diff is correctly installed on your computer and provides the expected results. In the next chapter, we present some standard examples of µ-diff scripts. Writing and solving the BIE using µ-diff . . . . . . . . . . . . . . . . . . 4.4.4 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The aim of this chapter is to provide some examples of multiple scattering problems solved by the µ-diff toolbox. The impenetrable case with a Dirichlet, a Neumann or a mixed of both boundary conditions set on the boundaries of the obstacles are fully treated. For penetrable obstacles, an example of the implementation of equation (1.36) is provided.

	Chapter 4
	Simple examples of multiple
	scattering problems solved with
	µ-diff
	U = InternalPotential(X, Y, O, a, M_modes, k, Density, Weight, OPTIONS);
	• ExternalPotential(..., 'Verbosity', VERBOSITY)
	VERBOSITY is a scalar, which when set to 0, stop displaying message (default = 1)

• ExternalPotential(..., 'OnBoundary',

ONBOUNDARY)

if ONBOUNDARY is set to 1 then the boundary is also covered (beware of the jump relation. . . ) (default = 0)

3.6 Examples available in µ-diff

Contents 4.1 The Dirichlet boundary-value problem . . . . . . . . . . . . . . . . . 4.1.1 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.2 The case of the EFIE . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.3 The case of the MFIE . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.4 The case of the CFIE . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.5 The case of the single-scattering preconditioned integral equation . . . . 4.1.6 The case of the Brakhage-Werner integral equation . . . . . . . . . . . . 4.1.7 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.8 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.9 Point source wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.4.1 Integral equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2 A more complex geometry . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.3

b p = OO p , b p = b p , α p = Angle( --→ Ox 1 , b p ),and, for q = 1, . . . , M , with q = p, b pq as the vector between the centers O q and O pb pq = O q O p , b pq = b pq , α pq = Angle( --→ Ox 1 , b pq ).

toolbox µ-diff requires the installation of the Matlab software (http://www.mathworks. com/), version 2011

[O,a] = RemoveDisk(O_old, a_old, ...);

where O_old and a_old are the centers and radii of the current geometry. Without optional argument, the function has no effect. The available arguments are • [O,a] = RemoveDisk(..., 'X', [X1, X2, .., XN]);

Remove all the disks with a center which has an abscissa equal to X1, X2, . . . , or XN.

• [O,a] = RemoveDisk(..., 'Y', [Y1, Y2, ..., YN]); Remove all the disks with a center which has an ordinate equal to Y1, Y2, . . . , or YN.

• Remove all the disks with a radius equal to a1, a2, . . . , or aN.

• [O,a] = RemoveDisk(..., 'Verbosity', VERBOSITY); set VERBOSITY to 0 to avoid display message, to 1 to only show results, and to > 1 to see everything (default).

• [O,a] = RemoveDisk(..., 'Tol', TOL); Tolerance used for the conditional statement (default 10 -10 ).

Example 1: Remove all the obstacles that are either on the row of with abscissa equal to 1 or on the column with an ordinate equal to 2.5

[O,a] = RemoveDisk(O_old, a_old, 'X', 1, 'Y', 2.5);

Example 2: Remove the obstacles centered at [START_REF] Antoine | On the numerical approximation of high-frequency acoustic multiple scattering problems by circular cylinders[END_REF][START_REF] Antoine | Integral Equations and Iterative Schemes for Acoustic Scattering Problems[END_REF] and [START_REF] Antoine | Alternative integral equations for the iterative solution of acoustic scattering problems[END_REF][START_REF] Antoine | Generalized combined field integral equations for the iterative solution of the three-dimensional Helmholtz equation[END_REF] [O,a] = RemoveDisk(O_old, a_old, 'XY', [2, 3; 5, 4]);

Example 3: Create a periodic placement of 11 × 11 unit disks, separated by a distance equal to 1, then remove the middle line and the middle column, as shown on Figure 3.3. The central disk is moreover centered at (0, 0).

Truncation of the Fourier series

To help the user, the formula (2.16) has been coded in µ-diff. The values of N p are stored in a Matlab vector of size [1 × N_scat] called M_modes in µ-diff. Computing M_modes only involves the radii of the disks and the wavenumber k, which is assumed to be created by the user is a scalar (e.g. = 2), then the resulting matrix Apq is the elementary matrix of the associated operator (e.g. L p,q ). If TypeOfOperator is a row (e.g. [START_REF] Alouges | A stable well-conditioned integral equation for electromagnetism scattering[END_REF][START_REF] Antoine | Alternative integral equations for the iterative solution of acoustic scattering problems[END_REF]), the sum of the two operators is computed (e.g. I p,q +M p,q ). Finally, the Weight quantity (of the same size as TypeOfOperator) is the constant that is used to multiply the block and hence

where Operator( ) is one of the integral operators.

Example 1: Build the single-layer block L p,q Apq = BlockIntegralOperator(Op, ap, Np, Oq, aq, Nq, k, 1);

Example 2: Build the whole matrix 0.5 × I p,q + N p,q , appearing in the MFIE (1.16) Apq = BlockIntegralOperator(Op, ap, Np, Oq, aq, Nq, k, [START_REF] Alouges | A stable well-conditioned integral equation for electromagnetism scattering[END_REF][START_REF] Antoine | Generalized combined field integral equations for the iterative solution of the three-dimensional Helmholtz equation[END_REF], [0.5, 1]);

Example 3: Compute the sum of the four blocks: 0.5×L p,q +1.5×M p,q +2.5×N p,q +3.5×D p,q Apq = BlockIntegralOperator(Op, ap, Np, Oq, aq, Nq, k, [START_REF] Antoine | On the numerical approximation of high-frequency acoustic multiple scattering problems by circular cylinders[END_REF][START_REF] Antoine | Alternative integral equations for the iterative solution of acoustic scattering problems[END_REF][START_REF] Antoine | Generalized combined field integral equations for the iterative solution of the three-dimensional Helmholtz equation[END_REF][START_REF] Antoine | Integral Equations and Iterative Schemes for Acoustic Scattering Problems[END_REF], [0.5, 1.5, 2.5, 3.5]);

Assembling the matrix A Now that the construction of an elementary block matrix is well-understood, building a global matrix is easy thanks to the common function A = IntegralOperator(O, a, M_modes, k, TypeOfOperator, Weight);

As for BlockIntegralOperator, the quantity Weight is optional and set to 1 by default. Roughly speaking, IntegralOperator creates a loop on all the obstacles p and q, launches the following command Apq = BlockIntegralOperator(O(:,p), a(p), M_modes(p), O(:,q), a(q), M_modes(q), k, Tpq, Wpq);

and places Apq in the expected elementary block matrix of the global matrix. The quantities Tpq and Wpq are given by TypeOfOperator and Weight such that (Weight is of the same size as TypeOfOperator and so Wpq follows the same rules as Tpq)

• If TypeOfOperator is a scalar then Tpq = TypeOfOperator.

• If TypeOfOperator is a row or a column vector then Tpq is an array given by Tpq =

TypeOfOperator.

• If TypeOfOperator is a matrix then Tpq = TypeOfOperator(p,q).

• If TypeOfOperator is a three-dimensional array, then Tpq is an array given by Tpq(:) =

TypeOfOperator(p,q,:).

Sparse matrix-vector product

A matrix-vector product Y = AX is done thanks to Y = SpMatVec(X, M_modes, ListOfOperators, Weight);

where Weight is optional. ListOfOperators is a Matlab cell (not an array! Make use of {•} instead of [•]) containing all the sparse matrices that the user wants to involve in the computational process. The linear combination of the operators is done at that time for each matrix-vector product.

Solving a linear system

Now that the matrices and the right-hand side have been built, the sparse linear system can be solved iteratively. If one uses for example the GMRES solver [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF], the syntax is then. . .

Post-Processing

Now that the system has been solved, the next step is to display some results. The µ-diff toolbox proposes some post-processing features such as computing the far-field (fast) or the near-field of the wave on a grid (Matlab meshgrid) or only at some points. In addition, other possibilities are offered such as drawing the disks on a figure or the incident field. Both the near-and farfield computations can be done for a linear combination of a single-and double-layer potentials L ρ + M λ, only one of them (e.g. L ρ) or with the same density (L + M )ψ, with a single command line. All post-processing functions are located in the PostProcessing/ directory.

Far-field

The far-field of the combination of the single-and double-layer potentials where the parameters are

• theta is the vector of receiving angles (in radians)

• Density is either ρ, λ or both. If Density is a column vector then Density(:,2) = Density(:,1) and only one density is used both for the single-and the double-layer potentials contribution. If Density has 2 columns, then Density(:,1) is considered to be the single-layer density and Density(:,2) the double-layer density.

• Weight is the weight vector to apply to the M volume integral operators (η p and γ p ). Example 2: consider two obstacles and compute

inside the obstacle, where the densities have already been computed

Incident wave

The incident waves can also be computed on a grid by using where, like for the IncidentWave function from the pre-processing part (see §3. The last line could also be (switch 'PlaneWave' with its identifier 1)

Uinc = IncidentWaveOnGrid(X, Y, k, 1, beta_inc);

Geometry: drawing the obstacles and creating mask matrices

The µ-diff toolbox also proposes some functions for the post-processing of the obstacles. Drawing the circular cylinders can be done by calling the following function If X and Y are obtained by meshgrid, let us define the mask matrix S by

This class of matrices is handy when plotting the potentials and is actually used by the two functions ExternalPotential and InternalPotential. The function MaskMatrixObstacles which provides the matrix M has the following syntax

Let us explain how to extract the boundary points. In the same way as for the MaskMatrixObstacles function, BoundaryOfObstacles leads to a matrix G where

Similarly, G is computed by

with Ω -= M p=1 Ω - p . We propose to solve this problem through various integral equations: the EFIE (1.15), the MFIE (1.16), the CFIE (1. [START_REF] Mautz | A combined-source solution for radiation and scattering from a perfectly conducting body[END_REF]) and the single-scattering preconditioned integral equations (1.32). We show how to use both the full and sparse storages of the matrices. Before starting, we recommend to use the single-scattering preconditioned integral equation as presented in §4.1.5. Indeed, the resulting system is well-posed and is well-conditioned leading to an efficient solution by a Krylov subspace iterative solver.

Pre-processing

Let us first consider a collection of three sound-soft unit circular cylinders. The wavenumber is k = 2π and the direction of incidence of the wave is β = 0 degree. The resulting µ-diff pre-processing code for setting these parameters is then %% Pre-processing % Three unit disks O = [-5, 0, 5; -2, 0, 2]; a = [1, 1, 1]; %Set the parameters... k = 1; %wavenumber beta_inc = 0; %incident angle (plane wave case) %Fourier series truncation parameter M_modes = FourierTruncation(a, k, 'Min', 1);

For each integral equation, we now present the assembly process, the computation of the solution and finally the post-processing of the computed wave fields. The common pre-processing part is the one described above. All the functionalities presented here are also available in the file BenchmarkDirichlet.m which is located in the Examples/Benchmark folder.

The case of the EFIE

This integral formulation reads as

where the first line is the integral equation representation of the exterior wavefield u and the second one is the surface integral equation to solve.

Dense storage

In µ-diff, the surface single-layer operator L and the incident plane wave field u inc | Γ are predefined quantities. If the full storage of the integral equation is used, then the direct solution of the resulting linear system can be obtained by the standard backslash Matlab operator \ %Right-hand side (plane wave) Uinc = PlaneWave(O, a, M_modes, k, beta_inc); %% Assembling %Matrix of the system (the two following lines are the same) L = SingleLayer(O, a, M_modes, k); %% Solving (here, direct) rho = L \ Uinc;

Sparse storage

For the sparse storage version, only the assembly process of the single-layer matrix and the system solution need to be modified as follows %Matrix of the system (the two following lines are the same) SpL = SpSingleLayer(O, a, M_modes, k); %% Solving (here, direct) rho = gmres(@(X)SpMatVec(X, M_modes, SpL), Uinc);

The case of the MFIE

The resolution of the scattering problem by the MFIE (1.16) leads to the integral equation representations

Dense storage

The MFIE operator

can be computed thanks to the frontal function IntegralOperator with two arguments: the type of the operators (for the identity operator and the double-layer potential operator N , see Table 3.2) and their associated weights (0.5 and 1). The post-processing part is exactly the same as for the EFIE since the surface equation is based on the volume single-layer integral representation.

Sparse storage

The sparse storage version is almost the same as for the dense storage except for assembling the matrix and solving the linear system. Indeed, the matrices I and N cannot be computed by the same function since the sparse function representations SpIntegralOperator and IntegralOperator cannot be summed together. It is however possible to add the identity to a sparse operator thanks to SpAddIdentity SpN = SpDnSingleLayer(O, a, M_modes, k); %Add I/2 to N: SpA_MFIE = SpAddIdentity(SpN, 0.5, M_modes) rho = gmres(@(X)SpMatVec(X, M_modes, SpA_MFIE), DnUinc);

The case of the CFIE

Let us now consider the well-posed and well-conditioned CFIE (see also Eq. (1.19))

Here, we fix the parameters to α = 0.5 and η = i/k.

Dense storage

The operator

Sparse storage

The sparse storage version changes compared to the dense one: the operators I/2 + N and L are computed separately and merged during the matrix-vector products. This is done in the SpMatVec function SpL = SpSingleLayer(O, a, M_modes, k); SpN = SpDnSingleLayer(O, a, M_modes, k); SpA_MFIE = SpAddIdentity(SpN, 0.5, M_modes) %% Solving and combining operators: rho = gmres(@(X)SpMatVec(X,M_modes,{SpL, SpA_MFIE}, [alpha * eta, 1-alpha]), B_CFIE);

The case of the single-scattering preconditioned integral equation

We strongly recommend to use the single-scattering preconditioned version of the EFIE, which is rigorously the same as the MFIE and CFIE and, up to an invertible operator, to any other boundary integral equation (see Proposition 1.9). The EFIE version is available in µ-diff and is represented as

Dense storage

In µ-diff, the quantity -L -1 u inc | Γ is provided by PlaneWavePrecond whereas L -1 L is obtained with PrecondDirichlet. The syntax for the dense version is then the following

Sparse storage

The sparse storage is here almost the same as for the dense version thanks to SpPrecondDirichlet SpPrecond = SpPrecondDirichlet(O, a, M_modes, k); %% Solving and combining operators: rho = gmres(@(X)SpMatVec(X,M_modes, SpPrecond), UincPrecond);

The case of the Brakhage-Werner integral equation

The Brakage-Werner integral equation for the Dirichlet problem (1.21) reads as

Dense storage

As in the previous chapter, we make here use of the common function IntegralOperator:

Sparse storage

For the sparse storage, the three operators I, L and N must be computed separately and merged during each matrix-vector products. Note that, in fact, the operators I and N can be merged together, as shown below using SpAddIdentity: SpL = SpSingleLayer(O, a, M_modes, k); %L SpM = -SpDoubleLayer(O, a, M_modes, k); %-M SpIminusM = SpAddIdentity(SpM, 0.5, M_modes); %0.5I -M %% Solving and combining operators: psi_BW = gmres(@(X)SpMatVec(X,M_modes, {SpL, SpIminusM}, [-eta_BW, 1]), Uinc);

Post-processing

The EFIE, MFIE, the CFIE and their preconditioned version share the same integral representation of the scattered field u as a single-layer potential only: u = L ρ. Their post-processing operations are the same. For the Brakhage-Werner integral equation, the post-processing is different since u = (-ηL -M ) ψ.

Radar Cross Section (RCS) and far field

Once the surface wavefield has been computed, the far field can be calculated by the following µ-diff commands, here for a discretization of [0, 2π[ with a step of 1 degree: The RCS can also be computed simply by: RCSSingleLayer = FarFieldToRCS(FSingleLayer); RCSBWIE = FarFieldToRCS(FBWIE);

or, if the far fields are not computed:

Near fields

The scattered field can also be computed on a grid. Even with a vectorization of the code, this functionality remains slow and cpu consuming. For the EFIE/MFIE and CFIE, only the single-layer potential must be computed whereas for the Brakhage-Werner integral equation, the double-layer must also be computed. In all the case, the computations are done in the exterior of the obstacles and ExternalPotential will do the job properly: The total field u T = u + uinc can also be calculated, thanks to IncidentWaveOnGrid which computes the incident wave on a grid, even inside the obstacles (on these point, 0 value must be set, using MaskMatrixObstacles).

Results

The figure 4.1 presents the results obtained with the considered configuration, that is 3 unit disks placed on (-5, -2), (0, 0) and (5, 2) with k = 1 and an incident plane wave of direction 0. On the figure is shown the obstacles, the history of convergence of the GMRES for the 5 integral equations for their dense and sparse storage, the radar cross section and also the near-fields.

Point source wave

Our toolbox µ-diff also provides a right hand side derived from a wave emitted by a point source. To solve that kind of problem, the right hand side is built thanks to PointSource and DnPointSource instead of respectively PlaneWave and DnPlaneWave (there is currently no equivalent of PlaneWavePrecond). For example, for the EFIE and a point source centered on (-10, 0), the right hand side must be computed as: The other change appears in the post-processing when computing the total field on a grid. Indeed, the incident wave is different and IncidentWaveOnGrid must get the right argument:

The rest of the code remains unchanged and the results obtained with a point source are shown on figure 4.2.

Mixing Dirichlet and Neumann boundary conditions

Let us consider the following situation where we mix Dirichlet and Neumann boundary conditions. The scatterer is composed of M D sound-soft and M N = M -M D sound-hard obstacles, leading to the scattering problem

For this problem, the preconditioned integral is not directly available. We can apply a Combined Field Integral Equation for the mixed problem (see equation (1.34)) and written as

where the matrix A is given by (1.35). We have

The matrix is particularly easy to build with µ-diff thanks to the frontal function IntegralOperator. To this end, two three-dimensional arrays, Assembling and Weight, are built such that

The indices in Assembling corresponds to the indices of the boundary integral operators (2 = L p,q , 3 = M p,q , 4 = N p,q , 5 = D p,q ). The assembling process is realized by IntegralOperator. The post-processing is then done by specifying to µ-diff how the density must be used: for the first M D obstacles, a single-layer potential is used, while for the others, a double-layer potential is required. The RCS function can simply do that. It just needs an array TypeOfOp of size N_scat × 2 such that

In a µ-diff script, this means Apply the single-layer potential (multiplied by 1) for the first M D part of the density and a double-layer potential (multiplied by 1) for the others .

Penetrable case 4.4.1 Integral equation

Let us now consider the case of penetrable obstacles, with a frequency k - p (possibly different) in each scatterer Ω - p , for p = 1, . . . , M . The problem then reads as:

where

Solving this problem can be done thanks to the integral equation presented in section 1.36 where the interior total field u -(equal also to the total field) and the exterior scattered field u + are represented as a single-layer potential:

The associated integral equation (1.36) is recalled to be

A more complex geometry

As already explained in §3.3.1, µ-diff also proposes some tools to place the obstacles and in a particular order, leading to original geometry, as the one shown previously on figure 3.

Writing and solving the BIE using µ-diff

Let us solve the problem using a dense storage and first write the matrix of the system, which consists in assembling the different operator into a larger matrix: Finally, we solve the system, here using a direct solver. The two densities ρ + and ρ -are here also extracted to simplify: rho = A\B; %% Extracting the "exterior" and "interior" densities rho_plus = rho(1:sum_modes); rho_minus = rho(sum_modes+1:end); clear rho;

Post-processing

Radar cross section

We compute now the radar cross section, which is based on the (exterior) single-layer potential: The radar cross section is shown on figure 4.5(b).

Near field

Let us now compute the near-fields by first building the grid: The quantity U is the scattered field outside the obstacles and the total field inside in the obstacles. The option OnBoundary is set to compute the potential also on the boundary, in case there is some points of the grid on it. To get the total field, we need to add the incident wave to U, outside the obstacles:

Appendix A

List of the µ-diff functions (alphabetical order)

IntOperators/Dense