Anton Arnold

Sjoerd Geevers

Ilaria Perugia

Dmitry Ponomarev

AN ADAPTIVE FINITE ELEMENT METHOD FOR HIGH-FREQUENCY SCATTERING PROBLEMS WITH SMOOTHLY VARYING COEFFICIENTS

Keywords: Helmholtz equation, scattering problem, time-domain wave problem, limiting amplitude principle, adaptive FEM, front-tracking mesh Mathematics Subject Classification 35J05, 35L05, 65M60, 65M50

HAL is

Introduction

We consider time-harmonic wave scattering problems in inhomogeneous media with smoothly varying material properties. Such problems become notoriously hard to solve when the angular frequency ω is large. To obtain a given accuracy with a standard finite difference or finite element method requires at least O(ω d ) degrees of freedom [START_REF] Babuška | Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers[END_REF], where d denotes the number of space dimensions. On top of that, standard iterative solvers and multigrid methods break down or converge slowly for high frequencies [START_REF] Ernst | Why it is difficult to solve Helmholtz problems with classical iterative methods[END_REF].

One way to reduce the computational complexity for large frequencies is by combining finite element methods with asymptotic methods [START_REF] Giladi | A hybrid numerical asymptotic method for scattering problems[END_REF][START_REF] Nguyen | A phase-based hybridizable discontinuous Galerkin method for the numerical solution of the Helmholtz equation[END_REF]. An asymptotic method, such as the geometrical optics method, is used to determine the wave propagation directions and a plane-wave finite element method is then used to solve the scattering problem. A drawback of this approach is that standard geometrical optics does not account for diffracted fields and incorporating diffraction phenomena typically requires an ad hoc approach. Instead of using an asymptotic method, one can also extract the dominant wave propagation directions from the solution of a lower-frequency scattering problem [START_REF] Fang | Learning dominant wave directions for plane wave methods for high-frequency Helmholtz equations[END_REF].

Another way to reduce the computational complexity is by using a classical finite difference or finite element method with a standard iterative solver, but in combination with a sweeping preconditioner [START_REF] Engquist | Sweeping preconditioner for the Helmholtz equation: hierarchical matrix representation[END_REF][START_REF] Engquist | Sweeping preconditioner for the Helmholtz equation: moving perfectly matched layers[END_REF][START_REF] Stolk | An improved sweeping domain decomposition preconditioner for the Helmholtz equation[END_REF]; for a more recent overview of several sweeping preconditioning methods, see [START_REF] Gander | A class of iterative solvers for the Helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz methods[END_REF]. The number of iterations then becomes nearly independent of the frequency, which means that the computational complexity scales almost as O(ω d ).

Instead of solving the scattering problem directly in the frequency domain, one can also solve the scattering problem in the time domain. The time-harmonic solution can be obtained from a solution to a time-dependent wave equation by exploiting the limiting amplitude principle [START_REF] Ladyzhenskaya | On the principle of limit amplitude[END_REF][START_REF] Vainberg | Principles of radiation, limit absorption and limit amplitude in the general theory of partial differential equations[END_REF][START_REF] Morawetz | The limiting amplitude principle[END_REF] or by applying a Fourier transformation in time. Classical time domain methods are the finite difference and finite element time domain methods [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[END_REF][START_REF] Taflove | Computational electrodynamics: the finite-difference time-domain method[END_REF]. Timedomain methods that are specifically devised for solving frequency-domain problems include the controllability method [START_REF] Bristeau | Controllability methods for the computation of timeperiodic solutions; application to scattering[END_REF][START_REF] Glowinski | A mixed formulation and exact controllability approach for the computation of the periodic solutions of the scalar wave equation. I. Controllability problem formulation and related iterative solution[END_REF], with its spectral version [START_REF] Heikkola | Controllability method for the Helmholtz equation with higher-order discretizations[END_REF] and its extensions [START_REF] Grote | On controllability methods for the Helmholtz equation[END_REF][START_REF] Grote | Parallel controllability methods for the Helmholtz equation[END_REF], the WaveHoltz method [START_REF] Appelo | Waveholtz: Iterative solution of the Helmholtz equation via the wave equation[END_REF], and the time-domain preconditioner of [START_REF] Stolk | A time-domain preconditioner for the Helmholtz equation[END_REF].

While both frequency-domain and time-domain methods are commonly used in practice and are expected to remain relevant in the future, this paper will focus on the time-domain approach. Some of the advantages of time-domain methods are that they are inherently parallel and straightforward to implement, without the need of storing Krylov subspaces and matrix factorisations and without the need of implementing linear solvers and moving absorbing boundary layers. However, for classical finite difference and finite element time domain methods, the number of degrees of freedom is at least O(ω d ) and the number of time steps is at least O(ω) due to the CFL condition, resulting in a computational complexity of at least O(ω d+1 ). In this paper, we present an adaptive finite element time-domain method that reduces the average number of degrees of freedom per time step to almost O(ω d-1 ), resulting in a computational cost that scales almost as O(ω d ).

The main idea of adaptive finite element time-domain methods is to automatically adapt the mesh over time in such a way that fine elements are used near the wave front and coarser elements are used away from the wave front. The mesh adaptation algorithms are typically driven by a posteriori error estimators/indicators. Adaptive finite element methods for the wave equation were studied for a conforming finite element discretisation in space combined with a discontinuous Galerkin discretisation in time [START_REF] Johnson | Discontinuous Galerkin finite element methods for second order hyperbolic problems[END_REF][START_REF] Li | Implementation and adaptivity of a space-time finite element method for structural dynamics[END_REF][START_REF] Thompson | Adaptive space-time finite element methods for the wave equation on unbounded domains[END_REF], the Crank-Nicholson scheme in time [START_REF] Bangerth | Finite element method for time dependent scattering: nonreflecting boundary condition, adaptivity, and energy decay[END_REF][START_REF] Gorynina | Time and space adaptivity of the wave equation discretized in time by a second-order scheme[END_REF][START_REF] Gorynina | An easily computable error estimator in space and time for the wave equation[END_REF], the implicit Euler scheme in time [START_REF] Bernardi | Time and space adaptivity for the second-order wave equation[END_REF][START_REF] Georgoulis | A posteriori L ∞ (L 2 )-error bounds for finite element approximations to the wave equation[END_REF], and the leap-frog scheme in time [START_REF] Georgoulis | A posteriori error estimates for leap-frog and cosine methods for second order evolution problems[END_REF]. Adaptive finite element methods based on a discontinuous Galerkin discretisation in space were presented and analysed in [START_REF] Adjerid | A posteriori finite element error estimation for second-order hyperbolic problems[END_REF][START_REF] Adjerid | A posteriori discontinuous finite element error estimation for two-dimensional hyperbolic problems[END_REF][START_REF] Adjerid | Asymptotically exact a posteriori error estimates for a one-dimensional linear hyperbolic problem[END_REF]. An anisotropic adaptive mesh refinement algorithm was studied in [START_REF] Picasso | Numerical study of an anisotropic error estimator in the L 2 (H 1 ) norm for the finite element discretization of the wave equation[END_REF]. Adaptive finite element schemes tailored for a given goal functional were studied in [START_REF] Bangerth | Finite element approximation of the acoustic wave equation: Error control and mesh adaptation[END_REF][START_REF] Bangerth | Adaptive Galerkin finite element methods for the wave equation[END_REF]. Adaptive finite element schemes have also been studied for other time-dependent problems such as the Stefan problem [START_REF] Nochetto | An adaptive finite element method for two-phase Stefan problems in two space dimensions. I. Stability and error estimates[END_REF] and nonlinear wave equations [START_REF] Alauzet | Transient fixed point-based unstructured mesh adaptation[END_REF]. Here, we present a new adaptive finite element time-domain method that is tailored for efficiently solving time-harmonic scattering problems and has the following distinguishing features:

• The source term is obtained from an incoming plane wavelet with compact support in the direction of propagation of width O(ω -1 ). • The adapted meshes are obtained from a set of nested meshes that are defined a priori.

• The mesh is not updated at each time step, but only after every m time steps. In the numerical examples, we have m ∼ 10 -100. • The time-harmonic field is obtained using an adapted algorithm for computing the Fourier transform in time that exploits the reduced number of degrees of freedom of the adapted meshes. In the numerical section, we present an implementation of the method using a fully explicit conforming finite element time-stepping scheme combined with the perfectly matched layer of [START_REF] Grote | Efficient PML for the wave equation[END_REF].

The paper is organised as follows: In Section 2, we explain how we solve the time-harmonic scattering problem using a time-domain approach. The adaptive finite element method for solving the time-dependent problem and the adapted method for computing the Fourier transform in time are then given in Section 3. Details of the numerical implementation and several numerical examples are given in Section 4. Finally, our findings are summarised in Section 5.

Solving the time-harmonic scattering problem in the time domain

We are interested in solving the time-harmonic wave scattering problem in d dimensions governed by the Helmholtz equation

-ω 2 (U S + U I ) -β -1 ∇ • (α∇(U S + U I )) = 0 in R d , (1a) [far field radiation condition on U S ], (1b)
where U S = U S (x) is the scattered wave field that needs to be resolved, U I = U I (x) is a given incoming plane wave, ω > 0 is the angular frequency, ∇ and ∇• denote the gradient and divergence operator, respectively, and α = α(x) ≥ α min > 0 and β = β(x) ≥ β min > 0 are two material parameters that are assumed to vary smoothly in space. Our main interest is the case where ω is large. We assume that there exists a bounded domain Ω in ⊂ R d such that α and β are constant, say α = α 0 and β = β 0 , in the exterior domain Ω ex := R d \ Ω in . We also assume that the incoming plane wave is of the form U I (x) = e iω(r•x)/c 0 , with i the imaginary unit (i 2 = -1), r a unit direction vector, and c 0 := α 0 /β 0 the wave propagation speed in the exterior domain.

We can rewrite (1) as

-ω 2 U S -β -1 ∇ • (α∇U S ) = F in R d , (2a)
[far field radiation condition on U S ], (2b)

with F := ω 2 U I + β -1 ∇ • (α∇U I ). Note that F = 0 in Ω ex .
Remark 2.1. In acoustic scattering, U S is the scattered pressure field, α is the reciprocal of the mass density ρ = ρ(x) of the medium, and β is the reciprocal of ρc 2 , with c = c(x) being the wave propagation speed of the medium.

A common way to obtain the scattered wave field U S is by solving a wave scattering problem in the time domain for a time-harmonic source term of the form F (x)e -iωt . If the limiting amplitude principle is valid [START_REF] Odeh | Principles of limiting absorption and limiting amplitude in scattering theory. II. The wave equation in an inhomogeneous medium[END_REF][START_REF] Èidus | The principle of limit amplitude[END_REF], the time-dependent scattered wave field converges to U S (x)e -iωt as t tends to infinity; see Appendix A.

Alternatively, we can compute the scattered wave field u S (x, t) for a suitable source term f (x, t) with compact support in space and time. Let F t denote the Fourier transform with respect to time, namely

F t [ϕ](ω ) = R e -iω t ϕ(t) dt. If F t [f ](•, -ω) = F , then it follows from the limiting amplitude principle that U S = F t [u S ](•, -ω); see Lemma A.2 in Appendix A.
Our proposed numerical method is based on this latter approach. In particular, we solve the scattered wave field u S (x, t) corresponding to a single incoming plane wavelet u I (x, t), defined such that U I = F t [u I ](•, -ω), and then compute U S = F t [u S ](•, -ω). By plane wavelet we mean a plane wave with compact support in the propagation direction. We describe this approach in detail in the three steps below.

Step 1. Defining the incoming plane wavelet. We consider an incoming plane wavelet of the form u I (x, t) := ωψ(ω(t -(r • x)/c 0 )), where ψ = ψ(ξ) is some real-valued, smooth function with supp(ψ) = [-ξ 0 , ξ 0 ], where ξ 0 > 0 is some constant independent of ω, and such that its Fourier transform

F[ψ] satisfies F[ψ](-1) = ξ 0 -ξ 0 e iξ ψ(ξ) dξ = 1.
The incoming wave field is thus a traveling plane wavelet of amplitude O(ω) and with a support of width 2c 0 ξ 0 ω -1 = O(ω -1 ). The Fourier transform of u I is given by

F t [u I ](x, ω) = F t [ωψ(ω(t -(r • x)/c 0 ))](x, ω) = ωF t [ψ(ω(t -(r • x)/c 0 ))](x, ω) = ωe -iω(r•x)/c 0 F t [ψ(ωt)](x, ω) = e -iω(r•x)/c 0 F t [ψ(t)](x, ω/ω) = e -iω(r•x)/c 0 F[ψ](ω/ω) and since F[ψ](-1) = 1, we therefore have F t [u I ](•, -ω) = U I . An illustration of u I is given in Figure 1.
Step 2. Solving a wave scattering problem in the time-domain. Having defined the incoming plane wavelet u I , we next solve the scattered wave field u S (x, t) given by the wave equation

∂ 2 t (u S + u I ) -β -1 ∇ • (α∇(u S + u I )) = 0 in R d × (t 0 , ∞), (3a)
[zero initial conditions on u S at t = t 0 ], (3b) where ∂ 2 t denotes the second-order time derivative, and t 0 := inf x∈Ω in (r • x)/c 0 -ξ 0 ω -1 is the time when the incoming plane wavelet first enters Ω in . We can rewrite this equation as

∂ 2 t u S -β -1 ∇ • (α∇u S ) = f in R d × (t 0 , ∞), (4a)
[zero initial conditions on u S at t = t 0 ], (4b) where f := -∂ 2 t u I + β -1 ∇ • (α∇u I ). Note that f has only support in Ω in × (t 0 , t f ), where t f := sup x∈Ω in (r • x)/c 0 + ξ 0 ω -1 is the time when the incoming plane wavelet leaves Ω in . Furthermore, since

F t [u I ](•, -ω) = U I , it follows that F t [f ](•, -ω) = F .
Step 3. Reconstructing the Helmholtz solution using a Fourier transform.

Since F t [f ](•, -ω) = F , it follows that, after extending u S by zero in R d × (-∞, t 0 ), we have F t [u S ](•, -ω) = U S .
Having computed the scattered wave field u S given by (4), we can thus reconstruct the solution to the Helmholtz equation given in (2) by computing the limit

U S = F t [u S ](•, -ω) = lim t→∞ t t 0 e iωτ u S (•, τ ) dτ. ( 5
)
The advantage of this approach is that the scattered wave field u S corresponds to an incoming plane wavelet of amplitude O(ω) and with a support of width O(ω -1 ). Motivated by geometric optics, we expect that, when ω is large and when there are no trapping modes, the solution u S is a travelling wave that has a steep gradient near the wave front, and a small gradient that is more or less independent of ω everywhere else; see also Figure 6 in Section 4. In the finite element approximation of u S , we can exploit this property and significantly reduce the computational cost by using an adaptive, time-dependent spatial mesh, where a fine mesh is used near the wave front and a coarser mesh, with mesh width independent of ω, is used elsewhere. Remark 2.2. We can readily extend the approach for a wave scattering problem that includes a sound-soft scatterer Ω sc . The wave scattering problem is then given by equation (1), but The approach remains identical, except that in equations (3) and (4) in Step 2, the spatial domain is R d \ Ω sc instead of R d and the additional boundary condition is of the form

u S = -u I on (t 0 , ∞) × ∂Ω sc .
Remark 2.3. With a slight modification, the approach can also be applied to a wave scattering problem of the form in (2) with F = F (x) an arbitrary external source term with bounded support. In Step 1, we then only need to define ψ, but not u I . In Step 2, we then solve equation (4) for f (x, t) = ωψ(ωt)F (x), with t 0 := -ξ 0 ω -1 the earliest time when f (•, t) is non-zero. Note that f only has support in supp(F ) × (-ξ 0 ω -1 , ξ 0 ω -1 ).

Step 3 remains unaltered. If F is a source term with very local support in space, such as a point source, we again expect that we can solve the time-domain problem efficiently using an adaptive mesh.

A description of the adaptive finite element method is given in the following section.

An adaptive finite element method

We aim to solve the wave equation given in (4) by using a finite element method with an adapted spatial mesh that is constantly updated over time. We simultaneously update the right-hand side of (5) by applying a discretised Fourier transform in time that exploits the reduced number of degrees of freedom of the adapted spatial meshes.

The main idea of the adaptive finite element method is as follows: we consider a finite computational domain Ω ⊃ Ω in with an absorbing boundary and split the time domain into small intervals (T j-1 , T j ) of length O(ω -1 ). Let u h denote the finite element approximation to u S . At the beginning of each time interval, we construct an adapted mesh T j of the domain Ω based on the current discrete approximation u h (•, T j-1 ). We then project u h (•, T j-1 ) into the finite element space of T j and solve the discretised wave equation on T j for the time interval (T j-1 , T j ).

For constructing T j , we aim for the coarsest possible mesh on which u S can still be approximated accurately during the time interval (T j-1 , T j ). To construct such a mesh, we assume that the current discrete approximation u h (•, T j-1 ) is an accurate approximation of u S (•, T j-1 ) and take into account that u S travels during the time interval (T j-1 , T j ) and is generated by the incoming wave u I . Let c max := sup x∈Ω α(x)/β(x) denote the maximum wave propagation speed. The construction of T j consists of the following steps:

• Coarsen. Coarsen the mesh T j-1 and project u h (•, T j-1 ) onto the discrete space associated with the coarser mesh. • Estimate. Compute the projection error, i.e. compute the difference between u h (•, T j-1 ) and its projection onto the coarser mesh. • Mark 1. Mark all elements of the coarser mesh where the projection error is above a certain threshold. Also mark all elements of the coarser mesh that overlap with the support of the incoming wave u I (•, T j-1 ). • Mark 2. Mark all elements of the coarser mesh that are within a distance c max (T j -T j-1 ) of elements that were marked in the first round. • Refine. Refine the coarser mesh at all marked elements to obtain T j .

To accurately approximate the Fourier transform of u h in time, we normally need to sample u h at all the degrees of freedom of a uniformly fine mesh at each time step. By computing the Fourier transform using an adapted algorithm, we can significantly reduce the average number of sampling points per time step.

A detailed description of the complete method is given in the following subsections.

3.1.

A finite element method with a time-dependent mesh. Consider the wave scattering problem given in (4). To approximate the scattered wave field u S using a finite element method, we consider a bounded polygonal domain Ω ⊃ Ω in , and impose an absorbing boundary condition or add an absorbing boundary layer at ∂Ω. Since, however, the main steps of our adaptive finite element method do not depend on the type of boundary condition, we consider in this section a zero Dirichlet boundary condition u S | ∂Ω = 0 in order to simplify the presentation. We thus consider a wave equation of the form

∂ 2 t u -β -1 ∇ • (α∇u) = f in Ω × (t 0 , ∞), (6a) u(•, t 0 ) = ∂ t u(•, t 0 ) = 0 in Ω, (6b) u = 0 on ∂Ω × (t 0 , ∞). ( 6c
)
Let T up = O(ω -1 ) be the time after which the mesh is updated and define T j := t 0 + jT up . Also, let T j denote the simplicial/square/cubic mesh of Ω used during the time interval (T j-1 , T j ]. For any mesh T , let U T denote the corresponding finite element space, given by

U T := {u ∈ H 1 0 (Ω) | u • φ E ∈ Û for all E ∈ T },
with φ E : Ê → E the affine element mapping, Ê the reference element, and Û the polynomial reference space. Also, let L T : U T → U T denote the discretisation of the spatial operator u → -β -1 ∇ • (α∇u) for a given mesh T . For simplicity, we consider here the operator L T defined such that

(βL T u, w) = (α∇u, ∇w) ∀w ∈ U T ,
where (•, •) denotes the L 2 (Ω) or L 2 (Ω) d inner product. A slightly different discretisation that allows for an explicit expression of L T is given in Section 4.1.2. Finally, let Π j denote a projection operator that projects into the space U T j . The semi-discrete finite element formulation can be stated as follows: for j = 1, 2, . . . , find

u j : [T j-1 , T j ] → U T j such that ∂ 2 t u j + L T j u j = f T j in Ω × (T j-1 , T j ), (7a) u j (•, T j-1 ) = Π j u j-1 (•, T j-1 ), (7b) ∂ t u j (•, T j-1 ) = Π j ∂ t u j-1 (•, T j-1 ), (7c) with u 0 (•, T 0 ) ≡ 0, ∂ t u 0 (•, T 0 ) ≡ 0, and f T j (•, t) ∈ U T j a discretisation of f (•, t).
For the time discretisation, we consider here the central difference scheme, although the adaptive method can also be applied to other time integration schemes. Let ∆t = T up /m denote the time step size, with m > 0 some positive integer. Also, let t n := t 0 + n∆t and u n j := u j (•, t n ). We approximate ∂ 2 t u j (•, t n ) by the central difference scheme

D 2 t u n j := u n+1 j -2u n j + u n-1 j ∆t 2 . (8)
Now, let n j := mj. The fully discrete finite element formulation can be stated as follows: for j = 1, 2, . . . and for n :

n j-1 -1 ≤ n ≤ n j , find u n j ∈ U T j such that D 2 t u n j + L T j u n j = f T j (•, t n ) for n : n j-1 ≤ n ≤ n j -1, (9a) u n j = Π j u n j-1
for n = n j-1 -1 and n = n j-1 , (9b) with u 0 0 ≡ 0 and u -1 0 ≡ 0. We can rewrite equation (9a) as

u n+1 j = -u n-1 j + 2u n j + ∆t 2 (-L T j u n j + f T j (•, t n )). ( 10
)
An extended time stepping scheme, that takes into account an absorbing boundary layer and the discretisation of f , is given in Section 4.1.2.

In practice, we cannot solve the wave equation for t → ∞, but have to stop at some finite time t stop = T jstop . To determine j stop , we use a stopping criterion of the form

sup x∈Ω |u n j j (x)| ≤ 0 , (11)
where 0 > 0 is an a priori chosen threshold value. In other words, we stop the computations when the scattered wave field is close to zero, which means it has almost completely left the computational domain.

An overview of how to implement the adaptive finite element method is given in Algorithm 1. Here, we use the following functions:

• T j = updateMesh(T j-1 , u n j-1 , t n ): computes the new mesh T j given the current mesh T j-1 , the current discrete wave field u n j-1 , and the current time t n . A detailed description of updateMesh and how to choose the initial mesh T 0 is given in Section 3.2 below.

• u n j := project(u n j-1 , T j-1 , T j ): computes the projection u n j = Π j u n j-1 .

• u n+1 j = doTimeStep(u n j , u n-1 j
, T j , t n ): computes the wave field at the next time step u n+1 j using the formula in [START_REF] Bangerth | Finite element approximation of the acoustic wave equation: Error control and mesh adaptation[END_REF].

• Stop(u n j , T j , t n ): returns true if t n > t f and if the stopping criterion given in ( 11) is satisfied. Returns false otherwise.

3.2.

Adapting the mesh. To construct adapted meshes T j , we define a priori a set of nested meshes {T 1 , T 2 , . . . , T K }, K ≥ 2, where T 1 is the coarsest mesh with a mesh width h 1 independent of ω, and T K is the finest mesh with a mesh width h K of order ω -1 or less. We assume that, in case of no mesh adaptation, the mesh T K is sufficiently fine for solving the wave equation with the desired accuracy.

We construct adapted meshes T j from elements in K k=1 T k . The initial mesh T 0 is chosen as the finest mesh T K . The algorithm T j = updateMesh(T j-1 , u n j-1 , t n ) for updating the mesh is given in Algorithm 2, which consists of the following functions:

• P = getParentElements(T ): returns the set of all elements in K-1 k=1 T k that are coarser than those of the given mesh T , i.e. it returns

P = {E ∈ K-1 k=1 T k \ T | E ⊃ E for some E ∈ T }.
• T = getChildElements(P): returns the mesh T ⊂ K k=1 T k , given its parent elements P = getParentElements(T ). In particular, getChildElements returns T = E∈P getSubElements(E) \ P.

Algorithm 1 solving the wave equation using a time-dependent mesh procedure solveWaveEquation T h ← T 0 set initial mesh u h ← 0, u new h ← 0, and u old h ← 0 initialise wave field for j = 1, 2, . . . do n = n j-1 at this point,

u h = u n j-1 , u old h = u n-1 j-1 , T h = T j-1 if stop(u h , T h , t n ) then return u h end if T new h ← updateMesh(T h , u h , t n ) T new h ← T j u h ← project(u h , T h , T new h ) u h ← u n j u old h ← project(u old h , T h , T new h ) u old h ← u n-1 j T h ← T new h T h ← T j for = 0, 1, 2, . . . , m -1 do n ← n j-1 + at this point, u h = u n j , u old h = u n-1 j u new h ← doTimeStep(u h , u old h , T h , t n ) u new h ← u n+1 j u old h ← u h u h ← u new h
end for end for end procedure

• T E = getSubelements(E): returns, for a given element E ∈ T k with k ≤ K -1,
the set of the elements in T k+1 that are a subset of E. • P * j-1 = markElements(P j-1 , u n j-1 , t n ): returns the set of all elements in P j-1 that need to be refined. We mark all elements E ∈ P j-1 that also have subelements in P j-1 . For the elements E ∈ P j-1 that have no further subelements in P j-1 , we only mark those for which the function needsRefinement(E, u n j-1 , t n ) returns true. Pseudocode of the function markElements is given in Algorithm 3.

• needsRefinement(E, u n j-1 , t n ): returns true if and only if E ∩ support(u I (•, t n )) = ∅, (12)
or

η E := sup x∈E |u n j-1 (x) -Π E u n j-1 (x)| > η 0 , (13)
where Π E denotes a projection operator that projects into the discrete space of E and η 0 > 0 is some threshold value defined a priori.

• P j = markNearbyElements(P * j-1 ): for k = 1, 2, . . . , K -1, returns all elements in T k that are within a distance c max T up of an element in

P * j-1 ∩ T k . The distance between two elements E 1 and E 2 is defined as dist(E 1 , E 2 ) := inf x∈E 1 ,y∈E 2 |x -y|.
An illustration of the mesh adaptation algorithm is given in Figure 8 below.

3.3.

Computing the Fourier Transformation. We can approximate the Fourier transform in (5) by a discrete Fourier transformation:

U S = F t [u S ](•, -ω) = ∞ t 0 e iωt u S (•, t) dt ≈ ∞ n=1 ∆te iωt n u S (•, t n ).
Furthermore, we can approximate u S (•, t n ) by the finite element approximation u n h , where u n h := u n j for n : n j-1 + 1 ≤ n ≤ n j and where u n j is the solution to the fully discrete problem Algorithm 2 update the mesh function updateMesh(T j-1 , u n j-1 , t n ) P j-1 ← getParentElements(T j-1 ) P * j-1 ← markElements(P j-1 , u n j-1 , t n ) P j ← markNearbyElements(P * j-1 ) T j ← getChildElements(P j ) return T j end function Algorithm 3 mark elements for refinement [START_REF] Bangerth | Finite element method for time dependent scattering: nonreflecting boundary condition, adaptivity, and energy decay[END_REF]. We assume that u S (•, t n ) ≈ 0 in Ω for n > n stop := n jstop . We then obtain the approximation

function markElements(P j-1 , u n j-1 , t n ) P * j-1 ← ∅ initialise the set of marked elements for E ∈ P j-1 do if getSubelements(E) ∩ P j-1 = ∅ then P * j-1 ← P * j-1 ∪ E mark E else if needsRefinement(E, u n j-1 , t n ) then P * j-1 ← P * j-1 ∪ E mark E end if end for return P * j-1 end function formulated in
U S ≈ U nstop h := nstop n=1 ∆te iωt n u n h .
We can compute U nstop h by setting U 0 h ≡ 0 in Ω and by using the recursive relation

U n h = U n-1 h + ∆te iωt n u n h for n = 1, 2, . . . , n stop . (14)
To accurately approximate U S on the entire computational domain Ω, we need to compute U nstop h on a globally fine mesh T K , which means that, if we use the formula in [START_REF] Èidus | On the principle of limiting absorption[END_REF], we would need to evaluate u n h at each time step at all the degrees of freedom of T K . Since u n h is only known at the degrees of freedom of some adapted mesh T j , this means we would need to interpolate u n h at the degrees of freedom of T K at each time step. We can significantly reduce the average number of interpolation points per time step in the computation of U nstop h by using an adapted space-time mesh.

Let Q be a space-time mesh for the space-time domain Ω × (T 0 , T jstop ) with space-time elements Q of the form

Q = E Q × (T j Q,0 , T j Q,1
). We choose the space-time elements such that

E Q ∈ T j for j : j Q,0 < j ≤ j Q,1 . Let χ Q (x, t) be the characteristic function given by χ Q (x, t) := χ E Q , t ∈ (T j Q,0 , T j Q,1 ], 0, otherwise, χ E (x) := 1, x ∈ E 0, otherwise.
A discrete version of χ E is given in Section 4.1.2 below. We have the partition of unity property

Q∈Q χ Q (x, t n ) = 1 for a.e. x ∈ Ω, ∀n : 1 ≤ n ≤ n stop . ( 15
)
We can therefore write

U nstop h = nstop n=1 ∆te iωt n u n h (15) = nstop n=1   Q∈Q χ Q (•, t n )   ∆te iωt n u n h = Q∈Q nstop n=1 χ Q (•, t n )∆te iωt n u n h =: Q∈Q ∆U Q . ( 16
)
Note that ∆U Q has support only in E Q . With a slight abuse of notation, we let ∆U Q also denote its restriction to E Q . Which definition is used will be clear from the context.

Let

n Q,0 := n j Q,0 and n Q,1 := n j Q,1 . We can compute ∆U Q by first setting ∆U Q = ∆U n Q,1 Q , where ∆U n Q,1 Q is computed by setting ∆U n Q,0 Q = 0 in E Q and by using the recursive relation ∆U n Q = ∆U n-1 Q + ∆te iωt n u n h | E Q for n : n Q,0 + 1 ≤ n ≤ n Q,1 .
Note that, in order to compute ∆U Q , we only need the values of u n h at the degrees of freedom corresponding to the spatial element E Q . To compute U nstop h

=

Q∈Q ∆U Q , we need to interpolate ∆U Q at the degrees of freedom of the finest mesh T K at E Q . However, we only need to do this once for each space time element.

To minimise the computational cost, we choose the space-time elements as large as possible. This means that we choose the time intervals (T Q,0 , T Q,1 ) as large as possible, namely such that E Q ∈ T j for all j : j Q,0 < j ≤ j Q,1 and such that E Q / ∈ T j for j = j Q,0 and j = j Q,1 + 1. An illustration of a space-time mesh is given in Figure 4 below.

In practice, we do not need to construct the space-time mesh explicitly. Let Q(j, E) denote the unique space-time element Q ∈ Q such that E = E Q and j Q,0 < j ≤ j Q,1 . We define ∆U j,E := ∆U Q(j,E) and ∆U n j,E := ∆U n Q(j,E) . We also define

U n j h := Q∈Q: j Q,1 ≤j ∆U Q .
We have that U 0 h ≡ 0 and we can compute U nstop h by computing, for j = 1, 2, . . . , j stop , the following:

U n j h = U n j-1 h + Q∈Q: j Q,1 =j ∆U Q , or, equivalently, U n j h = U n j-1 h + E∈T j \T j+1 ∆U j,E , (17)
with T jstop+1 := ∅. We can compute ∆U j,E = ∆U n j j,E by first setting ∆U 0 0,E ≡ 0 for all E ∈ T 0 and by computing, for j = 1, 2, . . . , j stop , the following:

∆U n j,E = ∆U n-1 j,E + ∆te iωt n u n h | E ∀E ∈ T j & n : n j-1 + 1 ≤ n ≤ n j , (18a) ∆U n j-1 j,E ≡ 0 ∀E ∈ T j \ T j-1 , (18b) ∆U n j-1 j,E = ∆U n j-1 j-1,E ∀E ∈ T j ∩ T j-1 . (18c)
Therefore, letting ∆U n j := {∆U n j,E } E∈T j , the algorithm for computing

U nstop h
is given in Algorithm 4, which consists of the following functions:

• U h = updateFT(U h , ∆U n j-1 , T j-1 , T j ): computes U h | E ← U h | E + ∆U n j-1,E for all E ∈ T j-1 \ T j . • ∆U n j = initialiseNewIncrements(∆U n j-1 , T j-1 , T j ): sets ∆U n j,E ← 0 for all E ∈ T j \ T j-1 and ∆U n j,E ← ∆U n j-1,E for all E ∈ T j ∩ T j-1 . • ∆U n j = updateIncrements(∆U n-1 j , u n h , T j , t n ): computes U n j,E = U n-1 j,E +∆te iωt n u n h | E for all E ∈ T j .
Algorithm 4 compute the discrete Fourier transform with respect to time

function computeFT({u n h } nstop n=1 , {T j } jstop j=1 ) U h ← 0 initialise Fourier transform ∆U h ← 0 initialise increments for j = 1, 2, . . . , j stop do n ← n j-1 at this point, U h = U n j-2 h , ∆U h = ∆U n j-1 U h ← updateFT(U h , ∆U h , T j-1 , T j ) U h ← U n h ∆U h ← initialiseNewIncrements(∆U h , T j-1 , T j ) ∆U h ← ∆U n j for = 1, 2, . . . , m do n ← n j-1 + ∆U h ← updateIncrements(∆U h , u n h , T j , t n ) ∆U h ← ∆U n j end for end for U h ← updateFT(U h , ∆U h , T jstop , ∅) U h ← U nstop h
return U h end function 3.4. Overview of the complete algorithm. To approximate the solution U S to the Helmholtz equation given in (2), we approximate the solution u S to the time-dependent wave equation in (4) using Algorithm 1 and then approximate the Fourier transform F t [u S ](•, -ω) using Algorithm 4. We can solve the wave equation and compute the Fourier transform simultaneously, resulting in Algorithm 5. This last algorithm gives a complete overview of the proposed method for solving the Helmholtz equation.

Numerical Examples

We present numerical examples for wave scattering problems in 1 and 2 spatial dimensions. Details of the finite element discretisation are provided for the 2-dimensional case. The discretisation for the 1-dimensional case can be readily deduced from the 2-dimensional case. All the numerical experiments presented in this section have been carried out in MATLAB R2017a.

4.1.

Absorbing boundary layer and finite element discretisation. We start by specifying how we impose an absorbing boundary layer and then provide details of the finite element discretisation. 4.1.1. Absorbing boundary layer. We consider the wave equation in ( 6) for a rectangular domain and add an additional absorbing boundary layer. Let Ω 0 = (-L 1 , L 1 ) × (-L 2 , L 2 ) ⊃ Ω in be the region of interest, and let Ω ABL be an additional absorbing boundary layer of width W surrounding Ω 0 . We define the computational domain by

Ω = Ω 0 ∪ Ω ABL = (-L 1 -W, L 1 + W ) × (-L 2 -W, L 2 + W ).
We apply the perfectly matched layer that was introduced in [START_REF] Grote | Efficient PML for the wave equation[END_REF] and further analysed in [START_REF] Kaltenbacher | A modified and stable version of a perfectly matched layer technique for the 3-d second order wave equation in time domain with an application to aeroacoustics[END_REF][START_REF] Baffet | Energy decay and stability of a perfectly matched layer for the wave equation[END_REF]. The resulting wave equation is given by

∂ 2 t u + (ζ 1 + ζ 2 )∂ t u + ζ 1 ζ 2 u -β -1 ∇ • (α(∇u + s)) = f in Ω × (t 0 , ∞), (19a)
∂ t s + Z 1 s + Z 2 ∇u = 0 in Ω × (t 0 , ∞), (19b) u(•, t 0 ) = ∂ t u(•, t 0 ) = 0 in Ω, (19c) s(•, t 0 ) = 0 in Ω, (19d) u = 0 on ∂Ω × (t 0 , ∞), (19e)
Algorithm 5 solving the Helmholtz equation using a time-domain approach procedure solveHelmholtzEquation T h ← T K set initial mesh u h ← 0, u new h ← 0, and u old h ← 0 initialise wave field U h ← 0 initialise Fourier transform ∆U h ← 0 initialise increments for j = 1, 2, . . . do n = n j-1 at this point,

T h = T j-1 , u h = u n j-1 , u old h = u n-1 j-1 also, at this point, U h = U n j-2 h , ∆U h = ∆U n j-1 if stop(u h , T h ) then U h ← updateFT(U h , ∆U h , T h , ∅) U h ← U n h return U h end if T new h ← updateMesh(T h , u h , t n ) T new h ← T j U h ← updateFT(U h , ∆U h , T h , T new h ) U h ← U n h ∆U h ← initialiseNewIncrements(∆U h , T h , T new h ) ∆U h ← ∆U n j u h ← project(u h , T h , T new h ) u h ← u n j u old h ← project(u old h , T h , T new h ) u new h ← u n-1 j T h ← T new h T h ← T j for = 0, 1, 2, . . . , m -1 do n ← n j-1 + at this point, u h = u n j , u old h = u n-1 j , ∆U h = ∆U n j u new h ← doTimeStep(u h , u old h , T h , t n ) u new h ← u n+1 j u old h ← u h u old h ← u n j u h ← u new h u h ← u n+1 j ∆U h ← updateIncrements(∆U h , u h , T h , t n+1 ) ∆U h ← ∆U n+1 j
end for end for end procedure with

Z 1 := ζ 1 0 0 ζ 2 and Z 2 := ζ 1 -ζ 2 0 0 ζ 2 -ζ 1 .
Here, s = s(x) = (s 1 (x), s 2 (x)) is an auxiliary vector field that has support only in Ω ABL , and

ζ 1 = ζ 1 (x 1 ) ≥ 0 and ζ 2 = ζ 2 (x 2 )
≥ 0 are two additional parameters that are nonzero only in Ω ABL . In particular, ζ 1 (x 1 ) is nonzero only for

x 1 : L 1 ≤ |x 1 | ≤ L 1 + W and ζ 2 (x 2 ) is nonzero only for x 2 : L 2 ≤ |x 2 | ≤ L 2 + W .
4.1.2. Finite element discretisation. To construct the adapted meshes, we use a set of nested rectangular meshes {T k } K k=1 that are constructed as follows. First, we define a sequence of mesh widths

h 1 > h 2 > • • • > h K . Mesh T k is then constructed using elements of size h k,1 × h k,2 , where h k,i = h k in the region |x i | < L i and h k,i = h K in the region |x i | > L i , for i = 1, 2.
In other words, T k is a Cartesian mesh of width h k in the main domain Ω 0 , but has thin elements of the finest resolution h K in the absorbing boundary layer Ω ABL . An illustration of the nested meshes is given in Figure 2.

To construct the conforming finite element space U T for a given mesh T , we use tensorproduct polynomials of degree at most p, namely the polynomial reference space is Û = span{x k 1 x 2 | k, ≤ p}. As degrees of freedom, we use the values at the nodes X T , where X T consists of the (p + 1) × (p + 1) tensor-product Gauss-Lobatto points of each element E ∈ T . In case T has a hanging node, X T contains the Gauss-Lobatto points corresponding to the

T k , k = 1, 2, 3, with {h k } = { 2 5 , 1 5 , 1 10 }, Ω 0 = (-1, 1) 2 , W = 1 5 .
coarsest edge adjacent to the hanging node, but not the points corresponding to the finer edges adjacent to the hanging node.

For the discretisation of s, we use the discontinuous finite element space S 2 T , where S T is given by

S T := {s ∈ L 2 (Ω) | s| Ω 0 ≡ 0 and s • φ E ∈ Û for all E ∈ T : E ⊂ Ω ABL }.
Let s E := s| E . The degrees of freedom of s ∈ S 2

T are given by s E (x) for all E ∈ T and all x ∈ X E , where X E := X {E} denotes the set of nodes on E.

We define the discretisation L T (u, s)

: (U T , S 2 T ) → U T of the spatial operator (u, s) → -β -1 ∇ • (α(∇u + s)) such that (βL T (u, s), w) T ,M L = (α(∇u + s), ∇w) T ∀w ∈ U T .
Here, (•, •) T denotes the approximation of the L 2 inner product (•, •) using the tensor-product (p + 1)-point Gauss-Lobatto quadrature rule for each element in T . Furthermore, (•, •) T ,M L denotes the approximation of (•, •) using a mass-lumping technique, i.e.

(u, w)

T ,M L = x∈X T u(x)w(x)σ x,T ,
where σ x,T := Ω w x,T (y) dy and w x,T ∈ U T denotes the nodal basis function corresponding to node x. We can give an explicit expression for L T :

L T (u, s)(x) = α(∇u + s), ∇w x,T T β(x)σ x,T ∀x ∈ X T .
For the projection operators, let T j-1 + T j denote the mesh constructed from elements of T j-1 and T j by always selecting the finest elements. We define the projection operators Π j : U T j-1 → U T j and Π S j : S T j-1 → S T j in such a way that

(Π j u, w) T j ,M L = (u, w) T j-1 +T j ∀w ∈ U T j , (Π S j s, w) T j = (s, w) T j-1 +T j ∀w ∈ S T j .
Furthermore, let E ∈ T k with k ≤ K -1, and let T E = getSubelements(E). We define the projection operator Π E : U E := U {E} → U T E , used for the refinement criterion in [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF], in such a way that

(Π E u, w) {E} = (u, w) T E ∀w ∈ U E .
We can give explicit expressions for these projection operators:

Π j u(x) = u, w x,T j T j-1 +T j σ x,T j ∀x ∈ X T j , (Π S j s) E (x) = s, w x,E T j-1 +T j σ x,E ∀E ∈ T j and x ∈ X E , Π E u(x) = (u, w x,E ) T E σ x,E ∀x ∈ X E .
Here σ x,E := E w x,E (y) dy and w x,E ∈ S E denotes the discontinuous nodal basis function corresponding to element E and node x.

For the time discretisation, let s n := s(•, t n ). We define u n+1/2 := 1 2 (u n + u n+1 ) and s n+1/2 := 1 2 (s n + s n+1 ). We also define the discrete second-order time derivative D 2 t u n j as in ( 8) and we define the discrete time derivatives D 2t u n and D t s n+1/2 as

D 2t u n := u n+1 -u n-1 2∆t , D t s n+1/2 := s n+1 -s n ∆t .
The fully discrete finite element formulation can then be stated as follows: for j = 1, 2, . . . , find u n j ∈ U T j for n : n j-1 -1 ≤ n ≤ n j and s n j ∈ S T j for n : n j-1 ≤ n ≤ n j such that

D 2 t u n j + (ζ 1 + ζ 2 )D 2t u n j + ζ 1 ζ 2 u n j + L T j (u n j , s n j ) = f T j (•, t n ) at all x ∈ X T j , (20a) D t s n+1/2 j,E + Z 1 s n+1/2 j,E + Z 2 ∇u n+1/2 j | E = 0 at all x ∈ X E , ∀E ∈ T j , (20b)
for n : n j-1 + 1 ≤ n ≤ n j , and

u n j = Π j u n j-1
for n = n j-1 -1 and n = n j-1 ,

s n j = Π S j s n j-1
for n = n j-1 , with u 0 0 ≡ 0, u -1 0 ≡ 0, and s 0 ≡ 0. We can rewrite [START_REF] Gander | A class of iterative solvers for the Helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz methods[END_REF] as

u n+1 j = z1 u n-1 j + z2 u n j + ∆t 2 (-L T j (u n j , s n j ) + f T j (•, t n )) z3 at all x ∈ X T j (21a) s n+1 j,E = Z-1 3 Z1 s n j,E + Z2 ∇u n+1/2 j | E at all x ∈ X E , ∀E ∈ T j , (21b)
where

z1 := -1 + 1 2 ∆t(ζ 1 + ζ 2 ), Z1 := I - 1 2 ∆tZ 1 , z2 := 2 -∆t 2 ζ 1 ζ 2 , Z2 := -∆tZ 2 , z3 := 1 + 1 2 ∆t(ζ 1 + ζ 2 ), Z3 := I + 1 2 ∆tZ 1 ,
with I ∈ R 2×2 the identity matrix.

For the discretisation of the source term f , recall that

f = -∂ 2 t u I + β -1 ∇ • (α∇u I ). One can check that -∂ 2 t u I + β -1 0 ∇ • (α 0 ∇u I ) ≡ 0 and therefore we have -β 0 β -1 ∂ 2 t u I + β -1 ∇ • (α 0 ∇u I ) ≡ 0. We can therefore write f = -(β -β 0 )β -1 ∂ 2 t u I + β -1 ∇ • ((α -α 0 )∇u I ).
We can discretise the time-and spatial derivatives in a similar way as before. The discrete source term f T (•, t n ) ∈ U T can then be given by

f T (x, t n ) := - β(x) -β 0 β(x) D 2 t u n I (x) - (α -α 0 )∇u n I , ∇w x,T T β(x)σ x,T ∀x ∈ X T ,
where u n I := u I (•, t n ). Note that, due to this discretisation, the discrete source term is still zero in the exterior domain Ω ex and for all t n > t f .

To discretise the characteristic function χ E for a square element

E = (x 1,E , x 1,E + h E ) × (x 2,E , x 2,E + h E ), we define χ E (x) =      1, x 1 ∈ [x 1,E , x 1,E + h E ) ∪ ({x 1,E + h E } ∩ {L 1 + W }) and x 2 ∈ [x 2,E , x 2,E + h E ) ∪ ({x 2,E + h E } ∩ {L 2 + W }), 0, otherwise. (22)
Then the partition of unity property ( 15) is valid for every x ∈ Ω (not just a.e. x ∈ Ω) and, in particular, for all nodes x ∈ X T K . 4.1.3. Algorithm for solving the Helmholtz equation. To solve the Helmholtz equation given in (2), we use Algorithm 5 with a few small modifications given below. These modifications take the absorbing boundary layer into account and ensure that all the steps are fully computable.

• At the start, we also initialise the auxiliary variables s h ← 0.

• After we compute u h ← project(u h , T h , T new h ), we also compute the auxiliary variable s h ← projectS(s h , T h , T new h ), where the function

s n j = projectS(s n j-1 , T j-1 , T j ) computes s n j = Π S j s n j-1 . • Instead of computing u new h ← doTimeStep(u h , u old h , T h , t n ), we now compute the pair (u new h , s h ) ← doTimeStep(u h , u old h , s h , T h , t n ), where we use the modified function (u n+1 j , s n+1 j ) = doTimeStep(u n j , u n-1 j , s n j , T j , t n ) that computes u n+1 j
and s n+1 j with the formulae in [START_REF] Georgoulis | A posteriori L ∞ (L 2 )-error bounds for finite element approximations to the wave equation[END_REF].

• For the stopping criterion given in [START_REF] Bernardi | Time and space adaptivity for the second-order wave equation[END_REF], we do not take the supremum over all x ∈ Ω, but instead we compute the supremum over all nodes x ∈ X T j . Similarly, for the refinement criterion given in (13), we do not take the supremum over all x ∈ E, but only over all nodes x ∈ X T E , where T E = getSubelements(E).

• For the function

U h = updateFT(U h , ∆U n j-1 , T j-1 , T j ), we now compute U h (x) ← U h (x) + χ E (x)∆U n j-1,E (x
) for all x ∈ X T K ∩ E and all E ∈ T j-1 \ T j , with χ E defined as in [START_REF] Georgoulis | A posteriori error estimates for leap-frog and cosine methods for second order evolution problems[END_REF]. 4.2. Numerical example in 1D. As a first numerical example, we consider the 1D domain Ω = (-1, 1), with spatial parameters

α(x) = 1 + 3(1 -2x) 2 (1 + 2x) 2 , x ∈ (-1 2 , 1 2 ), 1, otherwise, β(x) = 1,
and an incoming plane wave U I (x) = e iωx/c 0 , with c 0 = 1. An illustration of α is given in Figure 3.

For the numerical approximation, we consider an incoming wavelet of the form u I (x, t) = ωψ(ω(t -x/c 0 )), with

ψ(ξ) =    (ξ -π) 4 (ξ + π) 4 3840π(21 -2π 2 ) , ξ ∈ (-π, π), 0, otherwise.
An illustration of ψ is given in Figure 3.

At the boundary, we apply a perfectly matched layer of width W = c 0 π/ω, i.e. of width half a wave length. The equations for the absorbing boundary layer in 1D are the same as those in 2D given in [START_REF] Fang | Learning dominant wave directions for plane wave methods for high-frequency Helmholtz equations[END_REF], but with ζ 2 ≡ 0 and with a scalar field s(x) instead of s(x). We choose the damping parameter ζ 1 (x) = ζ(x) as in [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF], where

ζ(x) =    | log(R)| 3 2W W -(1-|x|) W 2 , |x| > 1 0, otherwise, (23)
with R = 10 -10 the expected artificial reflection.

For the time step size, we choose ∆t = T up /m, with m the smallest positive integer such that

∆t = T up m ≤ c CF L 2 λ * max ,
where the CFL number is chosen as c CF L = 0.9, and where λ * max is an upper bound of the largest eigenvalue λ max of the discrete spatial differential operator β -1 ∇ • α∇, given by

λ * max := 1 h 2 K α max β min sup u∈ Û \{0} (∇u, ∇u) { Ê} (u, u) { Ê} ,
where α max := sup x∈Ω α(x), and β min := inf x∈Ω β(x). A smaller value of T up results in fewer neighbouring elements being marked for refinement, and therefore in fewer degrees of freedom on average. However, it also means that the mesh needs to be updated more frequently. In other words, choosing T up very small might slow down the method due to computational overhead, while choosing T up very large might render the method less efficient due to many additional neighbouring elements being refined. In the numerical examples, we choose T up as half a time-period, i.e. T up = πω -1 .

For the spatial discretisation, we use quadratic elements (so degree p = 2). We compute for the time interval (t 0 , t stop ) = (t 0 , T jstop ), where t 0 = -0.5-πω -1 and where j stop is determined using the stopping criterion in [START_REF] Bernardi | Time and space adaptivity for the second-order wave equation[END_REF] with 0 = ω/100. For the mesh refinement criterion in [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF], we use η 0 = ω/100. An overview of the L 2 (Ω 0 ) error err 2 := U S -U h (•, T jstop ) Ω 0 and average number of degrees of freedom n DOF := 1 jstop jstop j=1 |X (T j )| is given in Table 1 for the adapted finite element method. The results are compared with the classical finite element method using a uniform mesh of width h = h K and using the same polynomial degree, time step size, and stopping time as for the adapted finite element method. From this table, we can see that the error of the adapted finite element method and classical finite element method behave very similarly, whereas the average number of degrees of freedom grows at a significantly slower rate for the adapted finite element method as the frequency ω increases. In particular, Tables 1 and4 1. Estimated L 2 (Ω 0 ) error and average number of degrees of freedom for the quadratic adapted (AFEM) and classical (FEM) finite element approximation to the 1D Helmholtz problem for different angular frequencies ω. To estimate the error, we take the numerical approximation on a uniform mesh of width h K /2 as reference solution.

To compute the errors in Table 1, we use the discrete solution obtained on a uniform mesh of width h K/2 as reference solution. However, this error does not take into account the error produced by the absorbing boundary layer or the truncation of the wave field at time t stop . To measure these errors, we compute the following numerical approximations and reference solutions.

• U h : the adapted finite element approximation considered in Table 1.

• U 1 : the reference solution used in Table 1.

• U 2 : similar to U 1 , but using the exact absorbing boundary condition ∂ t u + c 0 ∂ x u = 0 on ∂Ω 0 instead of an absorbing boundary layer. • U 3 : similar to U 2 , but using a time interval (t 0 , t 0 + 100) instead of (t 0 , t stop ) (t stopt 0 < 2.5 for all cases in Table 1). The difference U h -U 1 indicates the error due to the spatial and time discretisation, U 1 -U 2 indicates the error due to the absorbing boundary layer, and U 2 -U 3 indicates the error due to the truncation in time. An overview of these errors is given in Table 2. From this table, we can see that the error is dominated by the discretisation error, whereas the errors due to the absorbing boundary layer and truncation in time are negligible. We will use this as a motivation to also estimate the error by U h -U 1 Ω 0 in the 2D case. An illustration of u h for the case ω = 40π and an illustration of the corresponding spacetime mesh as described in Section 3.3 are given in Figure 4. 4.3. Numerical example in 2D: incoming plane wave. For the first 2D example, we consider a domain Ω 0 = (-1, 1) 2 , with spatial parameters and an incoming plane wave U I (x 1 , x 2 ) = e iωx 1 /c 0 , with c 0 = 1. An illustration of α is given in Figure 5.

ω U h -U 1 Ω 0 U 1 -U 2 Ω 0 U 2 -U 3 Ω 0 10π 1.
α(x 1 , x 2 ) = 1 + 3(1 -2 x 2 1 + x 2 2 ) 2 (1 + 2 x 2 1 + x 2 2 ) 2 , x 2 1 + x 2 2 ≤ 1 2 , 1, otherwise, β(x 1 , x 2 ) = 1,

Figure 5

For the numerical approximation, we consider an incoming wavelet of the form u I (x 1 , x 2 , t) = ωψ(ω(t -x 1 /c 0 )), with ψ as in the 1D example. At the boundary, we apply a perfectly matched layer of width W = c 0 π/ω, i.e. of width half a wave length. We choose the damping parameters ζ 1 (x 1 ) = ζ(x 1 ) and ζ 2 (x 2 ) = ζ(x 2 ), with ζ(x) as in the 1D example.

We compute for the time interval (t 0 , t stop ) = (t 0 , T jstop ), where t 0 = -0.5 -πω -1 . For the spatial discretisation, we use biquadratic elements. The time step size ∆t and the parameters η 0 and 0 are chosen in the same way as in the 1D case. We compute the L 2 (Ω 0 ) error and compare the results with a classical finite element method with a uniform mesh in the same way as we did for the 1D case. The results are presented in Table 3. From this table, we can see that, just as in the 1D case, the errors of the adaptive finite element method and classical finite element method behave similarly, whereas the average number of degrees of freedom grows at a significantly slower rate for the adapted finite element method as ω increases. In particular, Table 4 illustrates that for the adaptive finite element method, the average number of degrees of freedom grows almost linearly with ω instead of as ω 2 .

Snapshots of the time-dependent wave field for different frequencies are shown in Figure 6. This figure shows that the wavefront gets sharper as ω increases, while away from the wavefront, e.g. on the left of x 1 = 0, the wave field is similar for different frequencies. An 4. Estimated growth rate of the average number of degrees of freedom n DOF with respect to the frequency ω for the adaptive finite element method. The results correspond to the incoming plane wave problem in 1D and 2D. illustration of the total time-harmonic field and the error for the case ω = 40π is given in Figure 7. In other words, level(P)(x) returns the level of the finest element in P that contains x. If no element in P contains x, then level(P)(x) returns 0. An illustration of the mesh adaptation algorithm is given in Figure 8. Bottom: parent elements level(P j-1 ), marked elements level(P * j-1 ), and marked neighbour elements level(P j ). Here, P j-1 , P * j-1 , and P j are defined as in Algorithm 2.

F of the form

F (x 1 , x 2 ) := 200 λ 2 f 0 ρ(x 1 , x 2 ) 1 2 λ f 0 (ξ) := (ξ 2 -1) 4 , |ξ| ≤ 1, 0, otherwise,
where λ = 2π/k 0 denotes the wave length in the exterior domain, k 0 = ω/c 0 denotes the wave number in the exterior domain, ρ(x 1 , x 2 ) = (x 1 -x 1,0 ) 2 + (x 2 -x 2,0 ) 2 denotes the distance to (x 1,0 , x 2,0 ), and (x 1,0 , x 2,0 ) := (0.5, 0.5) is the position of the point source. Note that the support of F is centered at (x 1,0 , x 2,0 ) and has a diameter of one wave length λ. The domain, absorbing boundary layer, and spatial parameters α and β are chosen as in the previous example.

For the numerical approximation, we use the same spatial discretisation and parameters η 0 and 0 as in the previous example. The initial time is t 0 = -πω -1 . Similar to the previous examples, we compare the numerical results of the adaptive method with the classical finite element method for different frequencies. The results are presented in Table 5. Again, the accuracy of the adaptive and classical method are comparable, whereas the average number of degrees of freedom is significantly smaller for the adaptive method and grows almost linearly with ω instead of quadratically. Table 5. Estimated L 2 (Ω 0 ) error and average number of degrees of freedom for the biquadratic adapted (AFEM) and classical (FEM) finite element approximation to the 2D point source Helmholtz problem for different angular frequencies ω. To estimate the error, we take the numerical approximation on a uniform mesh of width h K /2 as the exact solution.

An illustration of the adaptive method is given in Figure 9. The left image of Figure 9 shows a sharp circular wave front generated by the point source and an adapted mesh that is only refined near this wave front. The right image shows the approximated time-harmonic wave field. For the numerical test, we only consider the case ω = 30π. We use two nested meshes with mesh size {h k } = { 1 10 , 1 150 } and with biquadratic elements, a mesh update time T up = 1/30, a mesh refinement criterion ξ 0 = 1 100 ω and a stopping criterion 0 = 5 100 ω. The initial time is t 0 = 0.4 -πω -1 , the number of time steps between each mesh update is m = 20, and the stopping criterion is triggered at t stop = T jstop , with j stop = 332. We compare the results of the adaptive finite element method and classical finite element method as in the previous examples. The results are listed in Table 6. An illustration of the adapted mesh and of the computed time-harmonic wave field is also given in Figure 10. Due to trapping, it takes much longer before the wave field vanishes, which makes time-domain approaches less efficient. Furthermore, as illustrated in Figure 10, a large region around the trapping area requires a fine mesh, which makes the adaptive method less efficient. The average number of degrees of freedom is still smaller than for the classical finite element method, since the adaptive method correctly determines in which part of the domain the wave field is active.

ω = 30π
err 2 n DOF AFEM 1.77e-01 9.95e+04 FEM 1.57e-01 3.31e+05 Table 6. Estimated L 2 (Ω 0 ) error and average number of degrees of freedom for the biquadratic adapted (AFEM) and classical (FEM) finite element approximation to the 2D trapping mode Helmholtz problem for angular frequencies ω = 30π. To estimate the error, we take the numerical approximation on a uniform mesh of width h K /2 = 1/300 as the exact solution.

Conclusion

We considered the time-harmonic acoustic scattering problem with smoothly varying coefficients for an incoming plane wave of angular frequency ω. The proposed method consists of solving the wave equation in the time domain for a single incoming plane wavelet using an adaptive mesh. The time-harmonic solution is then recovered by computing the Fourier transform in time using an adaptive algorithm that exploits the reduced number of degrees of freedom corresponding to the adapted meshes. We compared our adaptive finite element method to a standard classical finite element time domain method and show that the accuracy is comparable, whereas the average number of degrees of freedom for our adaptive method grows at a significantly smaller rate as the frequency ω increases. In particular, numerical examples indicate that the average number of degrees of freedom for the adaptive finite element method scales almost like O(ω d-1 ), with d the number of dimensions in space, instead of O(ω d ). Numerical examples also demonstrate that our method can be extended to include external source terms and sound-soft scatterers. The method, however, provides only a limited advantage in the presence of trapping modes.

Proof. Fix ω > 0 and let U be the solution to [START_REF] Glowinski | A mixed formulation and exact controllability approach for the computation of the periodic solutions of the scalar wave equation. I. Controllability problem formulation and related iterative solution[END_REF] with source term F = Fω . We need to show that Ũω = U . To do so, define G(t) := H(t)e -iωt , where H(t) denotes the Heaviside step function (H(t) = 1 for t ≥ 0 and H(t) = 0 for t < 0). Also, let * t denote the convolution operator with respect to time. If we apply G * t to (25), we obtain

∂ 2 t (G * t u) -β -1 ∇ • (α∇(G * t u)) = (G * t f ) in R d × (t 0 , ∞), (27a)
[zero initial conditions on (G * t u) at t = t 0 ]. (27b) Since f has finite support in time, we have that f (•, t) = 0 for all t > T for some T > t 0 . Therefore, (G * t f ) = e -iωt Fω for t > T . It then follows from the limiting amplitude principle that (G * t u)(•, t) converges to e -iωt U as t → ∞ in U = L 2 (Ω). In other words, lim t→∞ e iωt (G * t u)(•, t) = U . By definition of Ũω , we also have that lim t→∞ e iωt (G * t u) = Ũω and hence, Ũω = U .

Figure 1 .

 1 Figure 1. Illustration of the incoming plane wavelet u I (x, t) = ωψ(ω(t -(r • x)/c 0 )) at time t = -0.5 (left) and t = 0.5 (middle), and its Fourier transform F t [u I ](•, -ω) (right). We chose ω = 20π, r = (1, 0), c 0 = 1, and ψ defined as in Section 4.2.

Figure 2 .

 2 Figure 2. Illustration of Tk , k = 1, 2, 3, with {h k } = { 2 5 , 1 5 , 1 10 }, Ω 0 = (-1, 1) 2 , W = 1 5 .

Figure 3 .

 3 Figure 3. Illustration of parameter α(x) (left) and wavelet ψ(ξ) (right).

Figure 4 .

 4 Figure 4. Illustration of u h and the space-time mesh described in Section 3.3 for the 1D test case with ω = 40π. In the actual algorithm, we never explicitly compute the space-time mesh.

Figure 6 .

 6 Figure 6. Snapshot of u h at time t = t 0 + 0.70 for the 2D plane wave test case with ω = 20π (left), ω = 40π (middle), and ω = 80π (right).

4. 4 .

 4 Numerical example in 2D: point source. As a next 2D numerical example, we consider the Helmholtz equation of the form in (2) with a slightly smeared out point source

Figure 7 .

 7 Figure 7. Total time-harmonic field U tot = U h + U I (left) and error |U S -U h | (right) for the 2D plane wave test case with ω = 40π.

Figure 8 .

 8 Figure 8. Top: illustration of wave field u j-1 and the updated wave field u j at time T j = t 0 + 0.65, for the 2D plane wave test case with ω = 40π using three levels of nested meshes with {h k } = { 1 5 , 1 20 , 1 200 }. Bottom: parent elements level(P j-1 ), marked elements level(P * j-1 ), and marked neighbour elements level(P j ). Here, P j-1 , P * j-1 , and P j are defined as in Algorithm 2.

Figure 9 .

 9 Figure 9. Snapshot of the wave field u h at time T j = t 0 + 1 (left) and timeharmonic field U h (right) for the 2D point source test case with ω = 40π.

Figure 10 .

 10 Figure 10. Snapshot of the wave field u h at time T j = t 0 + 2 (left) and total time-harmonic field U tot = U h + U I (right) for the 2D trapping mode test case with ω = 30π.

  illustrate that the average number of degrees of freedom is almost independent of ω for the adaptive finite element method in 1D.

				AFEM			FEM
	ω	{h k }	T up m j stop	err 2	n DOF	err 2	n DOF
	10π	{ 1 5 , 1 50 }	1 10	28 20 1.40e-02 1.14e+02 1.00e-02 2.21e+02
	20π { 1 5 , 1 10 , 1 100 } 1 20	28 36 1.67e-02 1.43e+02 1.71e-02 4.21e+02
	40π { 1 5 , 1 20 , 1 200 } 1 40	28 63 3.92e-02 1.78e+02 3.93e-02 8.21e+02
	80π { 1 5 , 1 40 , 1 400 } 1 80	28 123 6.49e-02 2.12e+02 5.92e-02 1.62e+03
	Table

Table 2 .

 2 Estimated L 2 (Ω 0 ) error due to the spatial and time discretisation U h -U 1 Ω 0 , due to the absorbing boundary layer U 1 -U 2 Ω 0 , and due to the truncation in time U 2 -U 3 Ω 0 for the 1D test cases.

		40e-02	2.42e-04	4.50e-03
	20π	1.67e-02	2.00e-04	3.45e-03
	40π	3.92e-02	1.38e-04	3.94e-03
	80π	6.49e-02	1.21e-04	2.12e-03

Table 3 .

 3 Estimated L 2 (Ω 0 ) error and average number of degrees of freedom for the biquadratic adapted (AFEM) and classical (FEM) finite element approximation to the 2D plane wave Helmholtz problem for different angular frequencies ω. To estimate the error, we take the numerical approximation on a uniform mesh of width h K /2 as the exact solution.

				AFEM			FEM
	ω	{h k }	T up m j stop	err 2	n DOF	err 2	n DOF
	10π	{ 1 5 , 1 50 }	1 10	39 22 7.47e-03 1.93e+04 4.97e-03 4.88e+04
	20π { 1 5 , 1 10 , 1 100 } 1 20	39 39 8.71e-03 4.15e+04 8.43e-03 1.77e+05
	40π { 1 5 , 1 20 , 1 200 } 1 40	39 73 1.77e-02 9.33e+04 1.74e-02 6.74e+05
	80π { 1 5 , 1 40 , 1 400 } 1 80	39 139 3.48e-02 2.09e+05 3.44e-02 2.63e+06
				1D		2D
		ω	n DOF ratio rate	n DOF ratio rate
		10π 1.14e+02	1.93e+04
		20π 1.43e+02 1.25 0.33 4.15e+04 2.15 1.10
		40π 1.78e+02 1.24 0.31 9.33e+04 2.25 1.17
		80π 2.12e+02 1.19 0.25 2.09e+05 2.24 1.16

Table

*A. Arnold, S. Geevers, and I. Perugia have been funded by the Austrian Science Fund (FWF) through the project F 65 "Taming Complexity in Partial Differential Systems". I. Perugia has also been funded by the FWF through the project P 29197-N32. A. Arnold and D. Ponomarev were supported by the bi-national FWF-project I3538-N32.

Appendix A. Limiting amplitude principle Let U = U (x) be the solution to the Helmholtz problem

[far field radiation condition on U ], (24b) and let u = u(x, t) be the solution to the wave problem in the time domain

[zero initial conditions on u at t = t 0 ], (25b) with spatial parameters α = α(x) and β = β(x), source terms f = f (x, t) and F = F (x), and frequency ω > 0.

We assume that α(x) ≥ α min and β(x) ≥ β min for x ∈ R d , and that α(x) ≡ α 0 and β(x) ≡ β 0 for x ∈ R d \Ω in , where Ω in is a bounded domain and α min , β min , α 0 , and β 0 are positive constants. We also assume that f (•, t) and F are supported within Ω in .

Let U be a Hilbert space on a bounded domain Ω with Ω ⊃ Ω in . The limiting amplitude principle states that, if f is of the form f (x, t) = F (x)e -iωt , then u(•, t) converges in U to U e -iωt as t tends to infinity. In particular, we can define the limiting amplitude principle as follows:

Definition A.1 (limiting amplitude principle). Let u be the solution to the wave equation given in [START_REF] Gorynina | An easily computable error estimator in space and time for the wave equation[END_REF], with f (x, t) := F (x)e -iωt for all t > T for some T > t 0 , and let U be the solution to the Helmholtz equation given in [START_REF] Glowinski | A mixed formulation and exact controllability approach for the computation of the periodic solutions of the scalar wave equation. I. Controllability problem formulation and related iterative solution[END_REF]. The limiting amplitude principle states that [START_REF] Gorynina | Time and space adaptivity of the wave equation discretized in time by a second-order scheme[END_REF] lim

The following is known about the validity of the limiting amplitude principle, with supp(F ) ⊂ Ω ⊂ R d .

• For d = 3, ( 26) was derived in [START_REF] Tamura | Resolvent estimates at low frequencies and limiting amplitude principle for acoustic propagators[END_REF] for [START_REF] Gorynina | Time and space adaptivity of the wave equation discretized in time by a second-order scheme[END_REF] of the limiting amplitude principle is not valid [14, Sect. 3]; a modified form is currently under investigation and will be presented in a forthcoming paper. Whenever the limiting amplitude principle is valid for U = L 2 (Ω), we have the following result.

Lemma A.2. Let u be the solution to [START_REF] Gorynina | An easily computable error estimator in space and time for the wave equation[END_REF] with a source term f that has compact support in space and time. Extend u and f by zero to R d × (-∞, t 0 ) and define, for any frequency ω > 0, Ũω := F t [u](•, -ω) and Fω := F t [f ](•, -ω), where F t denotes the Fourier transform with respect to time. If the limiting amplitude principle is valid for U = L 2 (Ω), then, for any bounded domain Ω ⊂ R d with supp(f (•, t)) ⊂ Ω, we have that Ũω is the solution to the Helmholtz equation given in [START_REF] Glowinski | A mixed formulation and exact controllability approach for the computation of the periodic solutions of the scalar wave equation. I. Controllability problem formulation and related iterative solution[END_REF] with source term F = Fω .


Valérie Rolle
email: valerie.rolle@univ-nantes.fr

Sébastien Galliot
email: sebastien.galliot@univ-amu.fr

Benjamin Daniels

Joana Ferreira Pires

Keywords: Body modification, Tattooing, Rite of passage, Subjectivation, Deviance

Any individual or collective action with the purpose of intentionally changing all or part of the body can be termed "body modifications." This complicates the definition of a precise category, as all societies practice some degree of body modification albeit with extreme diversity in terms of both its visible form and its place within its social setting. Consequently, the inclusion of techniques such as body painting or plastic surgery within this generic definition is open to debate. Technical criteria, such as the depth and irreversibility of the act, do not significantly refine this definition. The decline of some irreversible, visible or mutilating practices (genital mutilation in Pohnpei, neck extension among the Kayan people and lip plates among the Surma) and recent expansion of other reversible and/or superficial practices (bodybuilding, fasting and tattooing) suggests that body modifications should be reconsidered from a more socio-anthropological than technical standpoint. Whereas observation and study of physical appearance predates the emergence of an anthropology of the body, it is the progressive consideration of the body as an interface and an operator for social relations (particularly since Mauss, 1973[START_REF] Mauss | Techniques of the body[END_REF]) that led to an understanding of the cultural basis of body modification. Nevertheless, analytical frameworks have still been adjusted throughout history.

From "primitive" to "governmentality"

From the 16 th century, European explorations in Africa, the Americas, Asia, and Oceania expanded and the bodies of indigenous peoples were observed through an ethnocentric lens imbued with Christian moral values. "Barbarian," "savage," "pagan": body modifications were interpreted from a utilitarian standpoint [START_REF] Kant | Critique of Judgment[END_REF](Kant [1790]], Darwin 1871) but very often simply disqualified with no further scrutiny. In the 19 th century, their study was developed due to the naturalistic paradigm specific to Western science and the taxonomic approach that intended to describe diverse lifestyles and human "races" via a study of morphological characteristics. Given this framework comprising physical anthropology and criminology, physical characteristics deviating from aesthetic standards inherited from Greek statuary were analyzed as divergent, inferior. The nature of body modifications, whether deliberate, accidental, or caused by environmental factors, was thus of little importance. In the early 20 th century, criticisms of evolutionary theories led anthropologists to study social institutions and practices as interdependent parts of a cultural system. Body modifications were henceforth considered an integral part of a society as a whole. Detailed ethnographic studies defined them as life crisis rituals or the corporeal expression of the esthetic values of a given society [START_REF] Ashley-Montagu | The Origin of Subincision in Australia[END_REF][START_REF] Bohannan | Beauty and scarification among the Tiv[END_REF][START_REF] Jeffreys | The winged solar disk or Ibo Itεi facial scarification[END_REF][START_REF] Van Gennep | Les rites de passage [The rites of passage[END_REF]. As demonstrated in the review articles of Benson (in [START_REF] Caplan | Written on the Body[END_REF] and [START_REF] Schildkrout | Inscribing the Body[END_REF], many ethnographic works describe numerous potential interventions on the skin and body as surfaces and materials molded by socialization processes. Conversely, more directly theoretical approaches to the body and its transformation remain rare, if not marginal. In classical anthropological works, the surface of the body is similar to a canvas, on which it is less the subject choosing to share its inner self with the society than the society acting towards a subject that must be transformed by its peers in order to be extracted from its primal state and take an acknowledged role within the group [START_REF] Clastres | De la torture dans les sociétés primitives[END_REF][START_REF] Godelier | La production du corps. Approches anthropologiques et historiques[END_REF][START_REF] Bourdieu | Les rites comme actes d'institution » [Rites as institutional acts[END_REF]. All body modifications thus convey a belonging to a given culture, which can be deciphered if the codes of that culture are known. Non-Western body modifications, studied in more recent monographs focus primarily on tattooing revivals or contemporary skin markings. They are described as beautification techniques, ritualized stages of the life cycle but also as empowerment practices, or markers of belonging to a given clan, ethnic group or society [START_REF] Mallon | Tatau: A history of Samoan Tattooing[END_REF][START_REF] Krutak | Tattoo Traditions of Native North America[END_REF][START_REF] Kuwahara | Tattoo. An anthropology[END_REF][START_REF] Salvador-Amores | Tapping ink, tattooing identities[END_REF], Te Awekotuku & Nikora 2011). In these diverse contexts, individuals are endowed with greater autonomy insofar as they convey more conscious political statements as well as postcolonial and patrimonial discourses, thereby contributing to reduce the so-called "great divide" between modern and pre-modern societies. Talal [START_REF] Asad | Remarks on the Anthropology of the Body[END_REF] stated that classic anthropology studied only body symbols (particularly within ritual contexts) at the expense of the perceptions and instruments that act directly on participants by way of their bodies. In his classic study, Turner (1980) demonstrated the importance of modifications to the skin and body surface as part of a progressive socialization. Between the lines, he initiated an exploratory study of the transformative acts performed on the body and their effects on the subject. The major epistemological obstacle to an anthropological study of body modifications is thus the individualistic nature of the body and the inaccessibility of the mind that inhabits it, as though the body and its perception were intimate and incommunicable. Alfred [START_REF] Gell | Wrapping in Images[END_REF] built on this critically using the psychoanalytical works of Didier Anzieu (1989) on the "skin-ego" (the development of subjectivity and basic relational skills through the skin as a medium) and examined processes that pre-date body modifications as a symbol, exploring how the thematic possibilities of the underlying technical schema of tattooing (wounding-scabbing-permanent ornament) were exploited in different ways by different social regimes. This approach inspired research on technical and instrumental aspects of body modification [START_REF] Galliot | Ritual efficacy in the making[END_REF]. Beyond these redevelopments, the works of Michel [START_REF] Foucault | Les techniques de soi » [Technologies of the Self[END_REF] remain decisive for the anthropological approach to body modifications. The notion of "governmentality," as a combination of power and subjectivity, repositions the social dimension of the body, the subject, and the individual into a biopolitical and historical reality. For Foucault, the subjectification process lies, on the one hand, in a set of social and technological devices (armies, schools, prisons, hospitals) that completely confiscate the availability of the body and its acts through some kind of capillarity via multiple micropowers, and, on the other hand, in "technologies of the self" with which subjects respond to this biopower and establish themselves as social beings. This association of "technologies of power" with "technologies of the self" can explain body modifications as diverse as bodybuilding, tattooing and plastic surgery. This approach of the body as a strategic space of intervention for political subjectivation is found in anthropological works such as [START_REF] Bloch | From blessing to violence. History and ideology in the circumcision ritual of the Merina of Madagascar[END_REF] or [START_REF] Warnier | The Pot-King: The Body and Technologies of Power[END_REF] and has guided more recent studies, such as those labeled "somatechnics" (Murray & Sullivan 2009). The criticism expressed by Nikki [START_REF] Sullivan | The Somatechnics of Bodily Inscription: Tattooing[END_REF] towards studies interpreting tattooing as an "external expression of the inner self," more generally be applied to approaches that consider body modifications as an attempt of personal narration. In that manner, the subject tries to define itself visà-vis others by exploiting the iconic power of a graphic expression with varying degrees of coding and idiosyncrasy sometime confining to floating signifier. This semiotic paradigm tends to be predominant in approaches to body modifications as negative, pathological, and compulsive markings (most psychology and criminology studies) as well as in those that perceive them as positive, deliberate, and liberating acts (most contemporary studies of body modifications).

From "narrative" to "industry"

In 1989, Clinton Sanders published a groundbreaking ethnographic work on contemporary Western tattooing, analyzing the practice -already partially rehabilitated -from an interactionist angle: how do tattooists and the tattooed pre-empt the social risks of acquiring a potentially discrediting marking by promoting a positive image of their occupation and themselves? Subsequent works developed the reflections initiated by Sanders, documenting the process of social rehabilitation of the tattoo, its appropriation by a recomposed population still partly consisting of sub-cultures labeled as "deviant," and the professionalization bordering or "artification" of the field. In the 1980s, the "Tattoo Renaissance" (see Fleming in [START_REF] Caplan | Written on the Body[END_REF] in the West involved improvements and changes in terms of technique (needle assemblies, colors), health and safety (cleaning and sterilization), iconography (increasingly diverse styles, customized designs) and representations (a history grounded on immemorial tribal roots, a redefinition of tattoo as a form of self-expression and artistic practice). The expansion of this new material and representational economy within and without the tattooing world (by way of international conventions, specialist magazines, cyberspace, as well as general news, literature, and media) had non-negligible effects on the social make-up of the populations involved. While tattooing never left the male-dominated, blue-collar, deviant milieux in which it was popularized in the West over the 20 th century, it additionally conquered the mainstream with an increasingly female clientele over the 1990s. Margo DeMello (2000) noted a polarization and hierarchization of the "tattoo community," whereby some tattoos were relegated to the working class (particularly bikers) while others were prized as a middle-class conception of the art. This normalization of the practice ensured its respectability: those with tattoos state that they very carefully considered their project and chose the best tattooist for a custom, unique, and artistic design that expressed their personality [START_REF] Irwin | Legitimating the First Tattoo: Moral Passage through Informal Interaction[END_REF]. Michael [START_REF] Atkinson | Tattooed: the Sociogenesis of a Body Art[END_REF], inspired by the figurational sociology of Norbert Elias, analyzed the enthusiasm for tattooing as views on the body began to change. Acceptance of tattoos in contemporary Western societies is based in a conception of identity dominated by self-expression. However, the reconciliation of body projects remains just as guided by the norms prevalent in the social spaces (friends, romance, family, work) through which the tattooed move. With one exception [START_REF] Sardo | De chair, d'encre et de quotidien. Une ethnographie du corps tatoué » [Flesh, ink, and everyday life. An ethnography of tattooed bodies[END_REF], this interpretation of the practice contrasts with French post-modern anthropological approaches in which the body is studied as space for individuals to express their "sovereignty," where the tattoo becomes an "intimate rite of passage" resulting in "rebirth" (Le Breton 2003). Anglophone works in the field conceptualize bodies more as a space of power, consistent with works in Cultural Studies, Gender studies and deviance. The commodification of tattoos as a rebellious symbol of identity would include them in the "supermarket of style" (as defined by Polhemus), particularly when the level of commitment to the practice is low [START_REF] Sweetman | Anchoring the (Postmodern) Self? Body Modification, Fashion and Identity[END_REF]. Nonetheless, the inevitable pain associated with acquiring a tattoo does distinguish it from consumption products displayed as "pure signs". Tattoos are thus an "ironic fad" [START_REF] Kosut | An Ironic Fad: The Commodification and Consumption of Tattoos[END_REF]), as they are permanent and incorporated while also conveying meaning and statements. Body modifications are sported as forms of self-affirmation bordering on empowerment, and may therefore reflect both an adherence to and rejection of dominant norms. Where the literature generally agrees that tattoos are consistent with female [START_REF] Atkinson | Pretty in Ink: Conformity, Resistance, and Negotiation in Women's Tattooing[END_REF] and male [START_REF] Gill | Body Projects and the Regulation of Normative Masculinity[END_REF] norms, feminist literature is keen to study their subversive potential (as defined by Butler). In a review of two centuries of tattooing among women, Margot [START_REF] Mifflin | Bodies of Subversion. A Secret History of Women and Tattoo[END_REF] demonstrated the extent to which the practice reflects an appropriation of the female body by men, although it became a means to reappropriate one's own body during feminist movements of the 1970s. The case of heavily tattooed women" clearly testifies to this contradiction: whereas tattoos allow an individual to reclaim her body from a subjective standpoint, the societal result is a re-establishment of gender norms. The "monster beauty" [START_REF] Braunberger | Revolting Bodies: The Monster Beauty of Tattooed Women[END_REF] of "heavily tattooed women" may be considered a source of distinction within the tattoo world. However, as described by Beverley [START_REF] Thompson | Covered in Ink: Tattoos, Women, and the Politics of the Body[END_REF], this alternative beauty is often subject to sanctions mitigated by compensatory (good behavior) or avoidant (tattoo covering) strategies in family, professional and public spaces. This construction and management of a transgressive body has also been observed in intersectional frameworks [START_REF] Santos | The Chicana Canvas: Doing Class, Gender, Race and Sexuality Through Tattooing[END_REF], see Hutter, Castellani & Thompson in Kloss 2020). The skin thus appears as a surface marked with "texts of race, gender, sexuality, class, and age, before it is marked with ink" (MacCormack 2006, 59). Commitment to "non-mainstream body modifications" often reflects adherence to a non-conformist lifestyle [START_REF] Wojcik | Punk and Neo-Tribal Body Art. Jackson: University of Mississippi Press. Figure captions Tables Related Articles (See Also) Art, Anthropology of; Cross-Cultural Aesthetics; Authenticity, Cultural; Display[END_REF]. Following on from Howard Becker's notion of "career deviance," Angus Vail (1999) examined the processes of affinity, affiliation, and meaning that allow the tattooed to become "tattoo collectors" and manage their ensuing deviant identity in spaces outside the tattoo world. Although "radical body modifications" (such as genital piercings, branding, cutting, and burning) are less common, such practices arise from a similar process of affiliation to sub-cultures that sometimes claim to be disaffiliated with mainstream society and its norms [START_REF] Myers | Nonmainstream Body Modification: Genital Piercing, Branding, Burning, and Cutting[END_REF]). However, protests by these sub-cultures (such as those self-proclaimed or labelled "radical queers", "modern primitives", and "cyberpunks") have little impact in subverting the social order, as their practices remain pathologized as a form of "self-mutilation". Moreover, their conception of identity and body modification is aligned with the Cartesian ("mind over matter") and neocolonial (authenticity of so-called "primitive" practices) ideologies of Western societies [START_REF] Pitts | The Flesh. The Cultural Politics of Body Modification[END_REF]. Unlike body modifiers, tattooists may claim the recognition of their craftsmanship as an art form thanks to their graphical skills. For Katherine [START_REF] Irwin | Saints and Sinners: Elite Tattoo Collectors and Tattooists as Positive and Negative Deviants[END_REF], this quest endeavors to acquire a privileged position within the tattoo world, which prefers to remain at the fringes of society by affirming a form of "positive deviance". This "internal artification" [START_REF] Kosut | An Ironic Fad: The Commodification and Consumption of Tattoos[END_REF]) went hand-in-hand with a stratification of the profession based on a generational gap and difference in educational background (with a loud minority being first trained in an art school). Anyhow, the profession is still a provision of services, which requires professional tattooists to run their studio as a business, negotiate the production of visual elements from a primarily uninitiated client base, and found and maintain an occupation governed by informal exchanges yet subject to a partial, forced professionalization process due to its regulation [START_REF] Rolle | L'art de tatouer. La pratique d'un métier créatif[END_REF]. David [START_REF] Lane | The Other End of the Needle: Continuity and Change among Tattoo Workers[END_REF] focused in particular on cultural transformations of the profession due to the redefined limits of and hierarchies within this group which is still predominantly organized around the model of "traditional tattoo shops."

Article ID e.g. abc0001, article title abc0002, article title abc0003, article title


Dong Keun Oh
email: spinhalf@kfe.re.kr

On Static H-formulation Adapted to the Critical Current Model of HTS Conductors

Keywords: Critical Current, Numerical Simulations, Hformulation, Superconducting Cables

A new attempt in H-formulation is presented, as an alternative to the static formulation of magnetic vector potential, to be applied to the self-consistent critical current model of HTS (High Temperature Superconducting) devices. Avoiding the intrinsic difficulty of the illposedness in static condition, the H-formulation is carefully implemented by critical combination of the known ideas which includes the BEM (Boundary Element Method) of the Laplace equation. On the merit of the vectorial field model, the efficacy of the alternative H-formulation is discussed not only in terms of optimal size of the air mesh, but also about modeling the background field, as improved to be easier than the previous work.

INTRODUCTION

The self-consistent critical current model has been recognized for its potential use to inspect the design of HTS (High Temperature Superconducting) devices by evaluating their current capability in terms of critical current, as typical in recent works, of the 2G REBCO (2 nd Generation Rare Earth Barium Copper Oxide) coated tapes [START_REF] Zermeño | A self-consistent model for estimating the critical current of superconducting devices[END_REF][START_REF] Zermeño | Open-source codes for computing the critical current of superconducting devices[END_REF]. Proposed to be an agile method of 2-d simplification, the idea has been founded upon the condition of static current, as equivalent to the uniform electric field (E-field) over the cross section or each slice of the cross sectional planes [START_REF] Zermeño | A self-consistent model for estimating the critical current of superconducting devices[END_REF].

Yielding rather efficient evaluations, a sensible trick of Pvalue indicator, which represents the domain-wise E-field uniformity, has been introduced not to solve the nonlinear E-J relationship directly, but to lead to the steady state solution restricted by and the domain-wise uniformity P [START_REF] Zermeño | A self-consistent model for estimating the critical current of superconducting devices[END_REF]. Under the idea, a numerically sound magnetic formulation should be laid, as the steady state transport is understood from the asymptotic behavior of magnetic induction. Thus, the equation of magnetic vector potential was favorably chosen for the underlaid magnetic formulation, because its static limit guarantees a well-posed problem by itself.

In this study, an alternative approach based on the field strength is attempted, instead of the magnetic vector potential, to substitute to the previous solution scheme based on the P-value. Taking advantage of direct access to the field quantity of interest, the alternative H-formulation is motivated to investigate the expected capabilities.

II. MODEL DESCRIPTION

In order to carry out the attempt, a particular variant of H-formulation has to be implemented to avoid the intrinsic difficulty of ill-posed problem, which appears in the static limit of H-formulation. In this sense, a legacy idea of penalization is invoked to prepare the H-formulation in weak form [START_REF] Bandelier | A new hformulation for nonlinear magnetostatics in R/sup 3[END_REF][START_REF] Bossavit | A mixed FEM-BIEM method to solve 3-D eddy-current problems[END_REF]. Substituting the penalized magnetic vector potential , the weak form is closed on the domain of interest , as written in terms of the field strength , the current density , and the arbitrary test field [START_REF] Bandelier | A new hformulation for nonlinear magnetostatics in R/sup 3[END_REF]. Besides, the boundary integral along is described by the reduced vector potentials, i.e., and , as defined for the background field to be , and the H-field on the boundary to be .

(

) 1
To complete the formulation, the second equation of simple identity emerges, which leads the free boundary condition constrained by the stray field in the outside of , virtually solving the Laplace equation of the reduced potentials through its Green's function kernels [START_REF] Bossavit | A mixed FEM-BIEM method to solve 3-D eddy-current problems[END_REF][START_REF] Costabel | Principles of boundary element methods[END_REF].

III. NUMERICAL MODEL

Following the formulation in Eq. 1, the numerical model of magnetic field is composed on the FEM (Finite Element Method) model of curl-free (Nédélec's) vectorial shape function [START_REF] Nédélec | Mixed finite elements in IR3[END_REF], and the BEM (Boundary Element Method) of simple linear shape function [START_REF] Costabel | Principles of boundary element methods[END_REF]. By virtue of the recent development of high-level modeling tool, the mathematical formulations are synthesized by the FreeFEM++'s lexicons in a straightforward manner [START_REF] Hecht | New development in FreeFem++[END_REF], and composed to the matrix formula (Eq. 2) by discretization, as codes in simple idioms, in terms of the unknowns, i.e., of the H-field and of the reduced vector potential on the boundary, [START_REF] Zermeño | Open-source codes for computing the critical current of superconducting devices[END_REF] where the current density and the external magnetic field are given to the load vector on the right hand side of Eq. 2. The matrices are obtained from the numerical integrations over the shape functions. Thus, and are matrices induced from the integrals on the domain , and and are from the boundary integrals along , coupling the vectorial (domain) elements and the boundary elements. Applying the Green's function formulae of the Laplace equation, is prepared as a dense matrix of so called BEM discretization [START_REF] Costabel | Principles of boundary element methods[END_REF].

Through the final equation (Eq. 2) of block matrix, one is able to evaluate the nonlinear solution of from the Pvalue based solution scheme, to be absolutely parallel to the previous work. Thus, the FreeFEM++ code is written as an iterative routine to make the domain wise P-value indicators converged [START_REF] Zermeño | Open-source codes for computing the critical current of superconducting devices[END_REF].

IV. BENCHMARK RESULT

As a benchmark test, the Roebel cable model of critical current is brought, which was introduced as the primary case

J c ( ⃗ B ) ⃗ A ⃗ H ⃗ A = 1 ν ( ⃗ J -⃗ ∇ × ⃗ H ) Ω ⃗ H ⃗ J ⃗ H′ ∂Ω a a 0 ⃗ H 0 ⃗ ∇ × ( a 0 ẑ ) ⃗ ∇ × (a ẑ ) + ⃗ H 0 ∫ Ω ⃗ H ⋅ ⃗ H′ d x 2 + 1 μ 0 ν ∫ Ω ( ⃗ ∇ × ⃗ H -⃗ J ) ⋅ ( ⃗ ∇ × ⃗ H′ ) d x 2 + ∮ ∂Ω ( a + a 0) ⃗ H′ ⋅ τ d l = 0 ∮ ∂Ω ∂a ∂n a′ d l + ∮ ∂Ω ( ⃗ H ⋅ τ -⃗ H 0 ⋅ τ )a′ d l = 0 Ω h a + a 0 [ [F] [C] T [C] [B] ] [ h a + a 0 ] = [ [V] ⋅ j (h) 2[C] ⋅ h 0 ] j h 0 [F] [V] Ω [C] [C] T ∂Ω [B] j (h)
in the reference [START_REF] Zermeño | A self-consistent model for estimating the critical current of superconducting devices[END_REF][START_REF] Zermeño | Open-source codes for computing the critical current of superconducting devices[END_REF]. For a better comparison, the same mesh structure of 10 HTS strands is applied to the alternative model of H-formulation [START_REF] Zermeño | Open-source codes for computing the critical current of superconducting devices[END_REF], while the outer diameter of air mesh is 13 times smaller, thanks to the BEM technique [START_REF] Bossavit | A mixed FEM-BIEM method to solve 3-D eddy-current problems[END_REF][START_REF] Costabel | Principles of boundary element methods[END_REF], than the original A-formulation. In spite of such downsizing in diameter, the total number of DoF is rather increased, but just becomes comparable to the reference of A-formulation at the application of the vectorial element (the Nédélec's element) [START_REF] Nédélec | Mixed finite elements in IR3[END_REF]; in detail, the A-formulation model has 6516 vertices and 26435 DoFs, and the H-formulation has 4111 vertices and 40950.

The benchmark test is carried out with respect to two different modes of criterion, one is in average over the 10 HTS tapes (AVG), and the other is the same P but at maximum among those ones (MAX) [START_REF] Zermeño | Open-source codes for computing the critical current of superconducting devices[END_REF]. As presented in the figures, the field profile of the H-field model just reproduces the reference model of A-formulation. In Table I, the result is summarized, in good agreements with the reference, arguing the potential of H-formulation as a practical alternative to the A-formulation. Computing the vectorial field quantity directly from the discretized equation, the developed method has its own merits, for instance, the background field is imposed explicitly into the model, whereas the original formulation cannot help representing the background effect indirectly as far field vector potential on the infinite boundary. Moreover, the far field approximation itself becomes unnecessary for the boundary constraint, because the BEM part virtually evaluates the stray field in the outside, so that it is possible to reduce the redundancy of air mesh significantly. On those positive features, the technical considerations in plan are discussed as following;

• A COMSOL complement of the same H-formulation is highly recommendable, taking into account its popular use in the HTS modeling community.

• Extending the developed method is also considerable, for instance, to the 3-d model, as well as to the hysteresis loss evaluation.

• Another means to impose the external field is feasible on the main domain of vectorial elements instead of the coupling integral along the boundary; it deserves to evaluate the cons and pros.

VI. CONCLUSION

The idea of H-formulation is successfully adapted to the outline of the previous method of critical current evaluation. In a straightforward manner, the developed mathematical formulations is translated to the numerical codes in the recentest FreeFEM++ platform. Through the benchmark result, the potential merits are discussed as a practical alternative to the established method for HTS modeling.

P = 1

 1 critical current model of the Roebel cable without background field; the field maps are obtained from the H-formulation model (upper), and from the reference model in A-formulation (lower).

2 .

 2 The critical current model of the Roebel cable with the background field of 0.1 T to y-axis; the field maps are obtained from the Hformulation model (upper), and from the reference model in Aformulation (lower).

ACKNOWLEDGMENT

This work was supported by the Korean Ministry of Science, ICT.