An automatic perfectly matched layer for acoustic finite element simulations in convex domains of general shape

This article addresses the efficient finite element solution of exterior acoustic problems with truncated computational domains surrounded by perfectly matched layers (PMLs). The PML is a popular nonreflecting technique that combines accuracy, computational efficiency, and geometric flexibility. Unfortunately, the effective implementation of the PML for convex domains of general shape is tricky because of the geometric parameters that are required to define the PML medium. In this work, a comprehensive implementation strategy is proposed. This approach, which we call the automatically matched layer (AML) implementation, is versatile and fully automatic for the end‐user. With the AML approach, the mesh of the layer is extruded, the required geometric parameters are automatically obtained during the extrusion step, and the practical implementation relies on a simple modification of the Jacobian matrix in the elementwise integrals. The AML implementation is validated and compared with other implementation strategies using numerical benchmarks in two and three dimensions, considering computational domains with regular and nonregular boundaries. A three‐dimensional application with a generally shaped domain generated using a convex hull is proposed to illustrate the interest of the AML approach for realistic industrial cases.

propagation of waves. Many approaches have been proposed and studied over the time. The basic impedance boundary condition is the simplest and cheapest approach, but it provides a rather poor accuracy. High-order absorbing boundary conditions (e.g., , infinite elements (e.g., , and perfectly matched layers (e.g.,  constitute families of more accurate techniques. In particular, the perfectly matched layers (PMLs) are easy to use and accurate for a limited computational cost, which makes them attractive for large-scale industrial simulations. They are frequently used in the aeroacoustic community. 17,18 With the PML technique, the computational domain is surrounded with an artificial layer where the outgoing waves are damped. The governing equations are modified in such a way that any outgoing wave is perfectly transmitted from the domain to the layer, whatever the angle of incidence, which is the key property of this approach. The PML has been introduced by Bérenger 12 for time-dependent simulations of electromagnetic waves with rectangular truncated domains. It has been rapidly applied to other physical waves and to time-harmonic problems. In the frequency domain, the PML equations are simply obtained by performing a complex coordinate stretch in the governing equations, 19,20 which has been reinterpreted as a modification of the metric tensor 21,22 and as a modification of the material properties. 23,24 Mathematical analyses of the Helmholtz problem with a PML have been proposed in References 22,[25][26][27][28][29]. The selection of the absorbing function of the PML, which is critical for practical applications, has been studied in References 14,[30][31][32] In particular, it has been shown in References 14,30,32 that specific unbounded absorbing functions provide high-fidelity solutions without requiring the tuning of parameters, in contrast with commonly used polynomial functions, for which some free parameters must be adjusted.
Most of the PML formulations have been derived for truncated domains with simple shapes, such as cuboids, spheres, and cylinders. With the standard approach, the complex coordinate stretch is performed in a specific coordinate system, that is, the Cartesian, spherical and cylindrical coordinate systems, respectively (see, e.g., References 27,33,34). Nevertheless, in many applications, an unnecessarily large free-space region must be discretized between the PML and the object under investigation should the truncated domain have one of these shapes. In order to improve the geometric flexibility of the PML, Teixeira and Chew 33,34 proposed a conformal PML that is suited to convex truncated domains with regular boundary. They used a specific local curvilinear coordinate system based on a local parametrization of the boundary, with normal and tangent coordinates. The conformal PML was studied for the Helmholtz equation in Reference 22. The approach has been applied to time-domain simulations with a discontinuous Galerkin scheme 35 and to time-harmonic simulations with Bernstein-Bézier finite elements. 36 Unfortunately, the conformal PML requires the exact knowledge of some geometric information that may not be available in practical cases, which limits the applicability of the approach. In three dimensions, the principal curvatures, the principal directions, and the normal of the surface between the domain and the layer must be known, as well as the distance between any point of the layer and that surface. These data are easily derived for simple geometries, such as spherical and ellipsoidal surfaces, but it is not the case for more complicated surfaces. The data could be estimated a priori by computational procedures, but it would be at the price of additional processing steps. The conformal PML is a priori not suited to polyhedric domains with nonregular boundary, which is another drawback of the approach. Indeed, in some applicative cases, only polyhedric meshed domains are available.
To address truncated domains of general shape, few alternative PML realizations have been investigated. Zschiedrich et al. 37 proposed a two-dimensional PML implementation for polygonal domains where the complex coordinate stretch is performed in a local prismatoidal coordinate system. In References 13,38, Ozgun and Kuzuoglu proposed a locally conformal PML finite-element implementation suited to convex domains of any shape. With their approach, the distance to the interface is estimated by solving a minimization problem, and the coordinate stretch is performed directly on the coordinates of the mesh nodes. That approach is very flexible and easy to implement, since it does not require further geometric data or any modification of the equations. However, the practical realization must be carefully handled, because it implies involved numerical approximations that can generate local instabilities. 39 To alleviate that issue, an alternative estimate procedure for the distance function has been proposed in References 39,40. In this work, we propose a novel finite element PML approach-which we call the automatically matched layer (AML) implementation-suited to convex domains of general shape with regular and nonregular boundaries. For convex domains with regular boundaries, this approach can be considered as a specific finite element implementation of the conformal PML that requires less geometric data than the original formulation. With the AML approach, the layer is generated at solver level by extruding the surface mesh of the domain boundary. The extrusion provides all the geometric data that are required, avoiding the need for supplementary preprocessing steps or costly optimization procedures. Then, the coordinate stretch is embedded into the finite element scheme thanks to a rather simple modification of the Jacobian transformation on the reference element in the elementwise integrals. The accuracy of the obtained layer is improved in comparison with the approach based on complexified mesh nodes. In particular, unbounded absorbing functions, which are proved to be very efficient, 14,30,32 can be used with the AML approach, in contrast to the approach of Ozgun and Kuzuoglu. Therefore, the whole strategy avoids some issues of the existing approaches while retaining the ease of implementation in existing codes and the accuracy by the use of unbounded absorbing functions.
This article is structured as follows. In Section 2, the Helmholtz problem and the classical conformal PML are introduced, together with a description of the coordinate system associated with the interface between the domain and the layer. The AML approach, with the mesh extrusion and the PML implementation, is explained in Section 3. In particular, the links and the differences with the existing approaches are discussed. In Sections 4 and 5, numerical results are proposed to validate the approach and to illustrate the use on a realistic 3D application. Finally, a conclusion is proposed in Section 6. Geometric data for reference benchmarks are provided in the Appendix.
Notation. Scalars are denoted with italic lower-case letters (e.g., k and r). Vectors are denoted with italic bold lower-case letters (e.g., x and n). Column matrices with the components of vectors in given coordinate systems are denoted with bold lowercase letters (e.g., x = [x 1 x 2 x 3 ] ⊤ and n = [n 1 n 2 n 3 ] ⊤ ). Matrices are denotes with bold uppercase letters (e.g., J).

CONFORMAL PML FOR THE HELMHOLTZ EQUATION
We consider a general time-harmonic wave propagation problem defined on the unbounded domain R d , with the dimension d = 2, 3. The field u(x) is the solution of the Helmholtz equation with the Sommerfeld boundary condition at infinity, where k is the (constant) wavenumber an f (x) is a source term with compact support on the region Ω dom ⊂ R d . For finite element simulations, the domain Ω dom is meshed and extended with a PML Ω pml to represent the propagation of waves toward the infinity. The interface between the domain and the layer is noted Γ int . The exterior border of the layer is noted Γ ext . We take the convention that the time-dependence of the fields is e −ı t , where is the angular frequency and t is the time.
The PML equations are obtained from the original equations by stretching a spatial coordinate in the complex plan, which introduces a directional damping of waves in the layer. For rectangular and cuboid domains, the stretching is performed in the Cartesian coordinate system. For circular and spherical domains, the radial coordinates are stretched in the corresponding coordinate systems. 27 This approach has also been used with elliptical and ellipsoidal domains. 41 It is however rather complicated to apply it for domains with more general shapes, because it requires the use of a global coordinate system associated with the interface Γ int . With the conformal PML, 22,33,34,42 the coordinate stretch is performed in a curvilinear coordinate system based on surfaces parallel to the interface Γ int . The layer then has a constant thickness , and the exterior border and the interface are parallel. This approach is more flexible and exhibits some features which simplify the practical finite element implementation.
The coordinate stretch used to derive the conformal PML is described in Section 2.1, together with the curvilinear coordinate system. The PML equation for the Helmholtz problem is derived in Section 2.2. The variational formulation is given in Section 2.3. Only the three-dimensional case is described, but it can be straightforwardly adapted to the two-dimensional case.

Coordinate stretch and curvilinear coordinates associated with int
The PML equation will be derived by performing a complex coordinate stretch along the normal to the interface Γ int between the domain and the layer. The domain is assumed to be convex. In this section, the interface is assumed to be sufficiently regular.
For each point x of the layer Ω pml , we consider the closest point p belonging to the interface Γ int . The distance between x and the interface Γ int is given by r ∶= dist(x, Γ int ) = ||x − p|| ∈]0, [. The point x can then be represented in a unique way as where n(p) is the unit normal of Γ int at p pointing toward the exterior of Ω dom . The complex coordinate stretch then consists in replacing the real coordinate r in Equation (3) with the complex functionr(r) defined as where (r) ≥ 0 is the so-called absorption function. Because of the coordinate stretch, any outgoing wave is damped in the normal direction n(p), that is, withx(r, p) ∶= p +r (r)n(p) and cos = n ⋅ k∕k.
To derive the PML equation, the coordinate stretch is performed in a specific curvilinear coordinate system ( 1 , 2 , 3 ) associated with the interface Γ int . For each point x( 1 , 2 , 3 ) of the layer, the coordinate 1 = r is the distance to the closest point belonging to Γ int , and the coordinates 2 and 3 are provided by a local parametrization of Γ int . The local parametrization is chosen in such a way that where t 2 ( 2 , 3 ) and t 3 ( 2 , 3 ) are the principal directions and 2 ( 2 , 3 ) and 3 ( 2 , 3 ) are the principal curvatures of the surface Γ int at p( 2 , 3 ) (see, e.g., References 2,35 for further details). We can then rewrite Equation (3) as and Equation (4) as̃1 The set of coordinates ( 1 , 2 , 3 ) forms an orthogonal curvilinear coordinate system, and the set of vectors (e 1 , e 2 , e 3 ) := (n, t 2 , t 3 ) constitutes an orthonormal frame. The two-dimensional version of the coordinate system is illustrated in Figure 1. The Jacobian of the transformation from the curvilinear coordinates ( 1 , 2 , 3 ) to the Cartesian coordinates (x 1 , x 2 , x 3 ) is obtained by differentiating Equation (7) with respect to i , which gives with the scale factors defined as h i ( 1 , 2 , 3 ) ∶= || x∕ i ||. For the system considered here, the scale factors are (see, e.g., Reference 33) Writing Equation (9) in Cartesian coordinates gives the Jacobian matrix where J x∕ is a 3 × 3 matrix and e i , n, and t i are column matrices with the Cartesian components of e i , n, and t i , respectively.

Derivation of the PML equation
The PML equation is straightforwardly obtained in the curvilinear coordinate system by writing explicitly the Helmholtz equation in curvilinear coordinates and by replacing 1 with̃1 ( 1 )in the obtained equation. Nevertheless, for the practical implementation in standard finite element codes, it is convenient to have the PML equation in Cartesian coordinates. In Reference 43, Matuszyk and Demkowicz described an elegant method to derive the PML equation and the variational formulation of the problem in Cartesian coordinates, when the complex stretching is performed in another coordinate system. Following their approach, the PML equation reads with the operator ∇ x , the matrix pml , and the scalar pml defined as where the matrix J pml is the Jacobian of the transformation from the Cartesian coordinates (x 1 , x 2 , x 3 ) to the complexified Cartesian coordinates (x 1 ,x 2 ,x 3 ) (i.e., the Cartesian components ofx) defined as It is well known that the PML equation (13) can be interpreted as an anisotropic absorber with specific complex material properties (see, e.g., some references in the electromagnetic community 23,24,33,40 ). The PML can also be interpreted as a space with a complex metric tensor. This point of view is discussed by Teixeira and Chew in Reference 34.
In the remainder of this section, we derive an explicit form for the Jacobian matrix (15). Using the chain rule, this matrix is factorized as This factorization highlights successive changes of variables: the transformation from Cartesian coordinates to curvilinear coordinates in the complex space, the use of the complex stretch to come back in the real space, and the transformation from curvilinear coordinates to Cartesian coordinates in the real space.
• The first matrix corresponds to J x∕ , defined in Equation (12), where the coordinate 1 is replaced with̃1 ( 1 ). In J x∕ , only the scale factors h 2 and h 3 depend on 1 . Using Equations (8) and (11), we can write where we have introduced the so-called stretching functions The first matrix can then be written as • The second matrix is the Jacobian of the transformation from the real curvilinear coordinates to the stretched curvilinear coordinates. Becausẽ1 where we have introduced the supplementary stretching function • The last matrix is the inverse of the Jacobian J x∕ , which reads This matrix verifies It contains all the information relative to the geometry and the complex stretch.

Variational formulation
The truncated problem consists in solving the Helmholtz equation (1) in the domain Ω dom , the PML equation (13) in the layer Ω pml , continuity conditions at the interface Γ int and a boundary condition at the exterior border Γ ext . In this work, we use a homogeneous Neumann condition on Γ ext . The region containing both the domain and the layer is denoted Ω. It is defined as an open set that verifies Ω = Ω dom ∪ Ω pml .
Since the PML equation is a Helmholtz-type equation, the variational formulation of the problem is obtained rather immediately with the standard approach. It reads Using the definitions of pml and pml , the formulation can be rewritten as This second formulation could also be obtained by performing the complex coordinate stretch and the change of variables directly on the classical bilinear form of the Helmholtz equation. 43 It highlights the important role of the Jacobian matrix J pml .

AN AUTOMATIC FINITE ELEMENT IMPLEMENTATION OF THE PML
The variational formulations (Equations 27 and 28) can be straightforwardly discretized with H 1 -conforming finite elements by using the standard approach. Nevertheless, the direct implementation of the obtained finite element schemes for generally shaped computational domains is more difficult because of the geometric data that are required. For a direct implementation, the material properties ( pml and pml ) or, equivalently, the Jacobian matrix J pml must be known at every point of the layer. They depend on the distance r, the normal direction n, the tangential directions t 2 and t 3 , and the curvatures 2 and 3 . Unfortunately, these geometric data are available only for simple geometries, such as circular, elliptical or spherical domains. For more general domains, they could be provided by CAD tools or by numerical procedures (see, e.g., Reference 8) at the price of additional preprocessing steps. Another issue is related to the approximation of the geometry by the mesh. Even with curvilinear elements, it could be impossible to align the finite element mesh with the interface Γ int and the exterior border Γ ext . A typical example is a circular domain meshed with straight-sided elements, which leads to a polygonal meshed interface. The definition of the geometric parameters then is ambiguous: should the distance r start at the meshed interface or at the circular interface? Finally, when only a mesh of the domain is available, without a reference geometry, we possibly have to address polyhedral boundaries, which falls a priori outside the scope of application of the conformal PML presented in the previous section.
With the proposed AML implementation, only the meshed interface is required. As in the previous section, the computational domain is assumed to be convex, but the meshed interface could be irregular. Our approach is based on a mesh extrusion generating the layer, and a local curvilinear coordinate system ( 1 , 2 , 3 ) associated with the meshed interface Γ int,h . By connecting directly the local coordinates ( 1 , 2 , 3 ) to the reference coordinates (u 1 , u 2 , u 3 ) used to compute the elementwise integrals, the issue of geometric data requirements is fully alleviated.
The mesh extrusion, the reference coordinate system and the standard finite element scheme are described in Section 3.1. The coordinate stretch and the curvilinear coordinates associated with the meshed interface Γ int,h are discussed in Section 3.2. The effective computation of the elementwise integrals for elements in the layer is discussed in Section 3.3.

Mesh extrusion and geometric data
We consider a general mesh for the convex domain Ω dom,h . The nature of the elements inside the domain Ω dom,h does not matter to present our approach, but we assume that the restriction of these elements on its exterior boundary, denoted Γ int,h , gives a conformal surface mesh made of linear straight elements (N geo = 1) or quadratic curvilinear elements (N geo = 2). Our approach could be extended to curvilinear elements of higher order or NURBS-enhanced methods, but this goes beyond the scope of this article.
To generate the mesh of a layer, Ω pml,h , the vertices and the second-order nodes belonging to the surface mesh Γ int,h are extruded along a direction n h , which should correspond to the exterior normal. Nevertheless, the definition of the

F I G U R E 2 Mesh extrusion for three cases in two dimensions with linear elements (cases (a) and (b)) and quadratic curvilinear
elements (case (c)). The mesh vertices and the second-order nodes are represented. The mesh extrusion is performed twice here (N pml = 2). The interpolated distance field r h ∈ [0, pml ] and extrusion directions n h are also shown exterior normal of Γ int,h is ambiguous at the points where the surface is not regular. This is the case, for instance, at the vertices and the edges of polyhedral surfaces. The extrusion direction n h is then chosen following the empirical rules: • for a given vertex of the surface, n h is the average of the normal of the elements touching the vertex; • for a second-order node defined in the middle of an edge or a face, n h is the average of the extrusion directions defined on the vertices touching the edge or the face.
The two-dimensional version of this strategy is illustrated in Figure 2. It is performed N pml times with a constant extrusion distance h pml . The obtained layer is structured and made of quadrangular elements (in two dimensions) or prismatic and hexahedral elements (in three dimensions) and has a thickness of pml = N pml h pml .
For each point x h ⊂ Ω pml,h , we shall need values for the extrusion direction n h , for the closest point belonging to the interface p h ⊂ Γ int,h and for the distance to the interface r h . During the mesh extrusion, these geometric data can be recorded at the nodes. Indeed, for each extruded node, n h and p h are the same as for the origin node belonging to the interface. For a node generated at the nth extrusion, the distance r h is equal to nh pml . At every node, we have To get values for the geometric data inside the elements we use a polynomial extrapolation. Let us consider a hexahedral element D e of the layer. The reference hexahedron is defined as where (u 1 , u 2 , u 3 ) are the reference coordinates. In the reference coordinate system, the polynomial representation of the position vector reads are the N geo th-order Lagrange functions and x e n 1 ,n 2 ,n 3 is the mesh node of the physical hexahedron D e corresponding to the set of indices (n 1 , n 2 , n 3 ). Similarly, we can obtain the extrusion direction n e , the closest point p e and the numerical distance r e at every point of the element D e by extrapolating the data recorded at the nodes. The extrusion direction and the closest points do not vary in the stretching direction, while the numerical distance varies only in that direction. In addition, the numerical distance varies linearly. If the reference coordinate u 1 is aligned with the stretching direction, we can then write with if the element D e has been generated at the nth extrusion, with n = 1 … N pml . Figure 2 presents the numerical distance field obtained using this strategy in two dimensions.

Coordinate stretch and curvilinear coordinates associated with int,h
Similarly to the approach of Section 2.1, a complex coordinate stretch is performed along the extrusion direction n h . The main differences here are that the interface Γ int,h can be irregular and n h is an interpolated field.
To perform the complex coordinate stretch in the extrusion direction, we introduce a local curvilinear coordinate system associated with the interface Γ int,h and based on the reference coordinate system. The first coordinate is the numerical distance. The second and third curvilinear coordinates are the second and third reference coordinates, which constitute a local parametrization of the interface. We then have e 1 (u 1 ) = r e (u 1 ), e 2 (u 2 ) = u 2 and e 3 (u 3 ) = u 3 . The position vector can be rewritten with a dependence on the curvilinear coordinates, Replacing the coordinate e 1 with the stretched coordinatẽe 1 ∕ık leads to the complex vector Thanks to these formula, we can obtain Jacobian matrices to switch between the Cartesian coordinate system and the curvilinear coordinate system in the real and complex spaces.

Finite element scheme and automatic PML implementation
The variational formulation (28) is solved by using an H 1 -conforming finite element method. The numerical solution u h is constructed on the mesh with hierarchical Lobatto shape functions (see, e.g., References 44,45).
Following the classical approach, the discrete variational formulation reads where a Jacobian matrix J pml,h = ( has to be defined. We could derive explicitly J pml,h by performing successive changes of variables between the Cartesian coordinates and the curvilinear coordinates in the complex and real spaces, following the strategy used in Section 2.1. Nevertheless, the computation is simplified by combining J pml,h with the Jacobian matrix of the mapping with the reference element, which is generally used to compute elementwise integrals. Let us consider an entry (A, B) of the matrix of the global system obtained by the discretization of the problem. The elementwise integral corresponding to the contribution of an element D e ⊂ Ω pml,h in the variational formulation (35) reads where Ψ A and Ψ B are global basis functions (defined on Ω h ) associated with the global unknowns u A and u B . We assume that A and B are such that D e belongs to the support of Ψ A and Ψ B . Using the reference mapping between the physical element D e and the reference element D ref , the integral can be rewritten as ] ⊤ is the gradient in the reference coordi- where the combined Jacobian matrix is defined as J ∶= J pml,h J ref .
The practical computation of J is rather straightforward. First, we refactorize the matrix J as with changes of variables between the complex Cartesian coordinates, the real curvilinear coordinates associated with Γ int,h and the reference coordinates. Using Equation (34), we easily derive the first Jacobian matrix, Using the definition of the curvilinear coordinates, we have with u 1 1 = h pml ∕2. Finally, we get Thanks to Equation (42), implementing the PML with the AML approach simply consists in adding an imaginary part to the Jacobian matrix used in the elementwise integrals. This part only requires a spatial representation of the distance r e (u 1 ) and the extrusion direction n e (u 2 , u 3 ), which the nodal values are automatically obtained during the mesh extrusion and then extrapolated on the layer, as explained in Section 3.2.
Although our approach is based on a complex stretch in a coordinate system associated with the mesh, it can be related to the conformal PML. Indeed, using Equations (21) and (22), we have, for the conformal PML, We obtain the corresponding Jacobian for the AML implementation (Equation 40) by taking ( 1 , 2 , 3 ) = ( 1 , 2 , 3 ) with the approximations n ≈ n e , 2 t 2 ≈ 2 n e and 3 t 3 ≈ 3 n e over each element D e of the layer. In the cases where the interface is smooth, the AML corresponds to a conformal PML with approximate geometric parameters. For nonsmooth interfaces, the AML can be seen as an approximate conformal PML with effective parameters at the corners. This approach is similar to the one proposed in Reference 8 for high-order absorbing boundary conditions, where effective geometric parameters were obtained with an empirical procedure to deal with corners. Despite the relation between the conformal PML and the AML, the mathematical analyses performed for the conformal PML do not apply straightforwardly here because of the geometric approximations. A comprehensive analysis should take into account errors due to the use of meshed geometries and approximate geometric parameters which depend on these meshes. This analysis is out of the scope of this article and is left for future work.
The approach proposed by Ozgun and Kuzuoglu in References 13,38, which consists in replacing the coordinates of the mesh nodes with the stretched coordinates, is a priori simpler than the approach described here. With that approach, the Jacobian matrix used in the elementwise integrals would be wherex e corresponds to a polynomial representation of the stretched coordinates. Because of the definition of the stretched position (Equation 34), the real part of this matrix is identical to the one of the matrix obtained with our approach (Equation 42), but the imaginary part is different. With the approach of Ozgun and Kuzuoglu, the stretched coordinates are interpolated with polynomial shape functions, and then differentiated numerically in Equation (44). Nevertheless, both steps introduce approximation errors, which can generate local instabilities. 39 The integrated absorbing function f ( 1 )that is in the imaginary part of Equation (34) is not represented accurately with the polynomial shape functions. For instance, the widely used quadratic absorbing function ( 1 ) = 2 1 requires third-degree polynomial shape functions for an exact representation, and a supplementary degree for an accurate numerical differentiation. The very efficient unbounded absorbing functions introduced by Bermúdez et al. 14 cannot be represented accurately with polynomials. By contrast, with our approach, the absorbing function ( 1 )and the integrated function f ( 1 )are evaluated exactly in the Jacobian matrix (Equation 42) without any polynomial interpolation and numerical derivation. The unbounded absorbing functions can be used efficiently.

NUMERICAL VALIDATION AND COMPARISON
In this section, the efficiency of the PML implementations is studied for truncated computational domains with regular and nonregular borders in two and three dimensions. The AML approach is compared with the direct implementation of the PML (when the exact geometric data are available) and to the implementation based on complex mesh nodes proposed by Ozgun and Kuzuoglu. 13,38 The reference benchmarks and the PML parameters used for the numerical simulations are detailed in Section 4.1. The validation and comparison results for truncated domains with regular and nonregular borders are presented in Sections 4.2 and 4.3, respectively.

Reference benchmarks and PML parameters
The numerical simulations represent the scattering of the plane wave u inc (x) = e ζkx with the propagation directionk = e x by a sound-hard object centered at the origin. In two dimensions, the scattering object is the disk of radius a sca = 1 and the resulting scattered field is where(r, )are the polar coordinates, J m is the mth-order Bessel function, H (1) m is the mth-order first-kind Hankel function, and m is the Neumann function which is equal to 1 for m = 0 and 2 otherwise. In three dimensions, the scattering object is the sphere of radius a sca = 1 and the resulting scattered field is where r = ||x|| is the radial coordinate, j m is the mth-order spherical Bessel function, H (1) m is the mth-order first-kind spherical Hankel function, and P m is the mth-order Legendre polynomial.
Truncated computational domains with different shapes are considered. For each case, the domain is meshed with the mesh generator Gmsh. 46 A layer surrounding the domain is generated by extrusion at run time by the solver, which is a Matlab code. The mesh of the domain is made of triangular or tetrahedral elements, and the extruded layer is composed of quadrangular or prismatic elements in two and three dimensions, respectively. Quadratic curvilinear elements are used for the accurate geometric representation of the problem (N geo = 2). The finite element solution is computed with hierarchical polynomial shape functions of maximal degree p = 1, 2 or 3 (see, e.g., References 44,45). The Neumann boundary condition n u = − n u inc is prescribed on the boundary of the scattering object, and a homogeneous Neumann condition is used on the exterior boundary of the extruded layer.
The accuracy of the PML critically depends on the absorbing function ( 1 ), the thickness of the layer and the spatial discretization (see, e.g., References 14,31,32). These parameters must then be chosen with great care. For a given discretization, increasing the thickness of the layer generally improves the accuracy, but it also increases the number of degrees of freedom and the size of the algebraic system. Therefore, in order to limit the computational cost, the layer should be as thin as possible with an optimized absorbing function.
Polynomial absorbing functions are widely used. These functions ensure a progressive damping of the outgoing waves in the layer. In particular, the cubic function is a very frequent choice, even though a supplementary free parameter is introduced. An approach to select the parameter consists in choosing a priori the reflection coefficient for outgoing waves with normal incidence. The reflection coefficient of the planar PML reads where inc is the angle of incidence. For the cubic function, the constant amplitude should then be = (2∕ ) ( ln R −1 0 ) , where and R 0 are given. Nevertheless, must not be too large to avoid any significant dispersion error due to the spatial discretization. 31,32 For practical applications, the parameter selection is a critical issue.
In an alternative approach, Bermúdez et al. 14 studied unbounded absorbing functions which cancel the reflection coefficient (48). In particular, the hyperbolic function provides good accuracy without requiring the tuning of supplementary free parameters for finite element simulations. 14,32 Nevertheless, because this function is singular on the exterior border of the layer, the numerical approximation of the elementwise integrals requires some care. For the AML implementation and the direct implementation of the conformal PML, we have used the hyperbolic function and Gauss-Legendre quadratures without any difficulty. Nevertheless, the implementation with complex mesh nodes has not provided satisfactory results with this absorbing function. That implementation is not suited to this kind of absorbing function (see discussion at the end of Section 3.3). In all the cases, the thickness of the layer is equal to N pml h pml , where N pml is the number of extrusions and h pml is the extrusion distance. The distance h pml is taken equal to h, the characteristic size of the mesh cells in the computational domain. The influence of N pml on the accuracy is studied in the following sections.

Numerical results for computational domains with regular boundary
The AML implementation is compared with the other approaches for two computational domains with regular boundary: the disk of radius a = 1.1 and the ellipse of semi-axes a x = 1.6 and a y = 1.1. The scattering disk of unit radius is placed in the middle of these domains. The simulations are performed with the wavenumber k = 25, the polynomial degrees p = 1, 2 and 3, the characteristic size of the mesh cells h = p ∕d , the wavelength = 2 ∕k and the resolution rate d = 20. The implementations are tested with different layer thicknesses (or, equivalently, different numbers of mesh extrusions N pml ) and both hyperbolic and cubic absorbing functions.
For a quantitative comparison of the implementations, we consider the relative L 2 -error between the numerical solution u num and the reference solution (45) on the computational domain, We also consider the relative L 2 -error between the reference solution and its L 2 -projection onto the finite element space, By Céa's lemma, this error corresponds to the best numerical solution that can be reached on each mesh, whatever the boundary treatment.
The AML implementation and the direct implementation of the conformal PML are firstly tested with the hyperbolic absorbing function hyp (49). The exact geometric parameters used for the direct implementation of the conformal PML are provided in the Appendix. The complex mesh nodes approach is not considered because it does not provide satisfactory results with hyp . In Figure 3, the relative L 2 -error is plotted as a function of the layer thickness for both implementations, with different polynomial degrees and both circular and elliptical domains. In all the cases, the AML implementation and the direct implementation give the same errors. The error is obviously larger with thin layers and small polynomial degrees. For these cases, the numerical dispersion is significant. Increasing the layer thickness decreases the error until a plateau that is close to the best interpolation error for p = 2 and 3. With p = 1, the plateau is substantially higher than the best interpolation error. The level of this plateau is higher in the elliptical case.
The PML implementations are now tested with the cubic absorbing function cub (47). When using this function, the parameter must be chosen. Here, it is taken equal to (2∕ ) , where the tuning parameter R 0 has to be adjusted. We first study the parameter selection with the direct implementation of the PML in the circular case. Figure  are similar with the hyperbolic and cubic functions for layers with several cells. The cubic function is outperformed by the hyperbolic function for layers with one or two mesh cells. In the remaining of the article, the cubic function is always used with R 0 = 10 −6 . Figure 5 shows the error curves obtained with the AML implementation and the complex mesh nodes approach equipped with the cubic function cub . The results of the direct implementation, not shown, are nearly identical to those with the AML. For thin layers, the errors are always smaller with the AML than with the complex mesh nodes approach. As discussed at the end of Section 3.3, approximation errors are introduced with the second approach because of the interpolation of the complex coordinates on the finite element mesh. This interpolation, which is the main difference between both implementations, explains the higher errors observed with the second one. In all the cases, increasing the (A) Circular domain-Absorbing function σ cub layer thickness decreases the errors until plateaus, which are at the same levels for all the implementations and both absorbing functions. These plateaus vary only with the polynomial degree p. In a nutshell, the direct implementation of the conformal PML and the AML implementation give the same results, which validates our approach. The implementation based on complex mesh nodes provides similar results, but only for large layers with a polynomial absorbing function that must be tuned.

Numerical results for computational domains with corners
The PML implementations are studied and compared for computational domains having corners with right and nonright angles. To analyze the influence of the angle on the error, polygonal and polyhedral domains are considered with different numbers of edges and faces, respectively. The circular and spherical domains, which are limit cases, are considered as well. In two dimensions, the scattering disk of unit radius is placed in the middle of polyhedral domains of midradius 1.65. In three dimensions, the scattering sphere of unit radius and polyhedral domains of midradius 2 are considered. In both cases, the polynomial degree is p = 2 and the resolution rate is d = 20. The wavenumber is k = 25 in two dimensions and k = 10 in three dimensions.
Only the AML implementation and the complex mesh nodes approach, which can address any kind of convex truncated domains, are systematically compared. The standard PML is suited only to domains with smooth borders and right angles. In the latter case, the coordinate stretch is performed in two or three directions, which must be orthogonal. Strategies to address nonright angles have been proposed for absorbing layers that are defined differently. [47][48][49] Nevertheless, these strategies do not apply to standard PML formulations, and they involve specific modifications of the computational scheme at the corners. By contrast, the implementations considered here correspond to the standard PML on the regular parts of the domain boundary, with a corner treatment implicitly embedded in the coordinate stretch based on the mesh extrusion. Figure 6 shows snapshots of the numerical solutions and the error distributions, together with the meshes, for numerical simulations performed with polygonal and circular domains. The AML implementation was used with N pml = 4 mesh extrusions and the absorbing function hyp . For the circular case, the error distribution is clearly dominated by the dispersion error of the finite element scheme. For the squared and pentagonal domains, the levels of error remain low, but structures in the error distributions indicate that spurious errors are generated at the boundary of the domain, and likely at the corners. Finally, the largest levels of error are observed with the triangular domain. For that case, the extruded mesh cells close to the corners have the worst aspect ratio. In addition, the exterior border is rather close to the scatterer. This is the most challenging case for the PML implementations.
For a quantitative comparison of the implementations, the L 2 -error is computed for different configurations with the polygonal and circular domains. We consider the AML implementation with hyp , and the complex mesh nodes approach with cub . As a reference, the standard Cartesian PML is tested with both absorbing functions for the configurations with the square domain. The relative error is plotted as a function of the layer thickness for the different configurations in Figure 7. Let us note that the relative projection error has nearly the same value for all the computational domains. Different observations can be made: • When the hyperbolic function is used (Figure 7(A)), the errors obtained with the Cartesian PML are very close to the projection error, and do not vary with the layer thickness. With the AML implementation, the errors are slightly larger for the circular and polyhedral domains with obtuse angles. The errors significantly larger, by at least one order of magnitude, in the triangular case. This corroborates the preliminary qualitative study.
• When the cubic function is used (Figure 7(B)), the errors obtained with the Cartesian PML remain close to the projection error, but only for thick layers (i.e., with N pml ≥ 5). Either the parameter R 0 is not optimal, or the cubic function cannot provide results as good as the hyperbolic function for thin layers. With the complex mesh nodes approach, the errors are similar or slightly higher than the errors obtained with the Cartesian PML.
• In all the cases, increasing the layer thickness decreases the errors until they plateau. The plateaus are at the same levels for the AML implementation and the complex mesh nodes approach, despite the different absorbing functions. The plateau is much higher in the triangular case.
To summarize, both the AML implementation and the complex mesh nodes approach are equivalent for large layers. The results are satisfactory, except for the triangular case, where the border is close to the scattering object and the angles F I G U R E 6 Polygonal computational domains: mesh (left), real part of the numerical solution (center) and corresponding error distribution (right) obtained with the AML implementation. The layer is generated with N pml = 4 mesh extrusions and the hyperbolic absorbing function is used are acute. For thin layers, which allow for a reduction of the computational cost, combining the AML with the hyperbolic function is the most efficient approach.
The implementations are finally tested in three dimensions, with tetrahedral, cubic, octahedral, icosahedral, and spherical domains. In Figure 8, the error curves are shown for the AML implementation with hyp and the complex mesh nodes approach with cub . The error curves corresponding to the Cartesian PML with the cubic domain are shown as well. These results confirm the observations made for the two-dimensional results. Both the AML implementation and the complex mesh nodes approach are efficient for sufficiently large layers. For thin layers, using the AML implementation with the hyperbolic absorbing function is the best strategy. For the icosahedral and the spherical cases, only one mesh extrusion is necessary.

APPLICATION WITH AN AUTOMATICALLY GENERATED CONVEX DOMAIN
To illustrate the interest of the AML approach for realistic applications, we consider the scattering of a plane wave by a three-dimensional scattering object: the shark submarine. The original geometry 1 and the mesh of the submarine are shown in Figure 9. The submarine is placed in the middle of a computational domain where the pressure field is computed. In order to minimize the number of degrees of freedom, the exterior surface of this computational domain is generated automatically by using a complex hull algorithm. The exterior surface mesh is then extruded to generate the absorbing layer. The final model is set up and ran using the simulation package Simcenter 3D (version 2019.2) 50 developed by Siemens Industry Software, in which the AML implementation has been developed.
We describe step-by-step the process to get the mesh of the computational domain and the layer. First, a surface mesh of six-noded triangular elements is generated on the original CAD geometry of the submarine, which is of unit total length. The characteristic length of the elements is h ≈ 15 × 10 −3 . Then, an approximate convex hull of this surface mesh is computed by applying the quickhull algorithm. 51 To create the exterior surface of the computational domain, the convex hull is scaled in the exterior normal direction, and it is homogeneously triangulated with the element length h. The minimal distance between the submarine and the scaled surface is equal to 2h. In a third step, a quadratic tetrahedral mesh of the volume between the submarine and the scaled surface is generated. The final mesh of the domain, shown in the bottom right part of Figure 9, is made of 302 803 10-noded tetrahedral elements and 439 884 nodes. As expected, close to the nose and at the fine tips, the exterior surface of the domain is only 2-elements away from the submarine surface. Finally, the mesh of the PML region, composed of 10-noded prismatic elements, is extruded by the solver at run time. The required geometric data are generated automatically (see discussion in Section 3.1). The extrusion length is h pml = h, and the layer thickness is pml = N pml h pml . Only the AML implementation and the hyperbolic absorbing function are considered.
The simulations were performed for two incident plane waves hitting the submarine frontally (k inc = [0 0 − 1] ⊤ ) and with an oblique incidence ( 2), respectively, with the wavenumber k = 40, the polynomial degree p = 2 and the resolution rate d = 20. The real part of the scattered fields obtained in the computational domain is represented in Figure 10 for both incident fields in the case N pml = 5. The orientation of the global axes is visible on the figures. The computations were performed on a DELL Precision 7520 laptop, with 3.10 GHz clock speed and 32 Gb of RAM. The most computationally intensive operations were the system assembly and the factorization, which required respectively 4.5 and 62 seconds on four threads. The factorization required 11.6 Gb of memory.
We first study the influence of the PML on the accuracy by comparing simulations performed with three layer thicknesses: N pml = 1, 5 and 10. Figure 11 shows the directivity of the scattered field on a circle of radius 2 in the xz-plane, for both incidence directions and each layer thickness. For each incidence, the results obtained with the different thicknesses are very close. For layers with N pml = 5 and 10, the curves are superimposed. Therefore, a thin absorbing layer, with a thickness of only few mesh cells, is sufficient to get accurate results, although the interface between the domain and the layer is very close to the scatterer.
The computational cost, in terms of runtime and memory storage, is directly related to the total number of degrees of freedom, and then to the volume of the computational domain. Using the automatically generated convex hull to define the domain, instead of a more standard geometry (e.g., cuboid, cylinder, or spheroid), reduces the volume of the computational domain, yielding significant savings on the computational resources. To illustrate this, we have generated cuboidal and cylindrical computational domains aligned with the main axis of the submarine. For each geometry, the minimal distance between the exterior boundary and the submarine is equal to 2h, and the exterior surface mesh is extruded to generate a layer with N pml = 5. The characteristics of the resulting finite element models are compiled in Table 1.
The model with the convex envelope is significantly less computationally intensive than the cuboidal and cylindrical models. The volume and the number of degrees of freedom in the computational domain are reduced by approximately a factor of three and two, respectively. This difference may be explained by the presence of several mesh refinement regions close to the scatterer, which capture the small geometric details of the submarine. The PML unknowns represent a significant proportion of the total number of degrees of freedom in all the cases: 50%, 44%, and 40% for the convex, cuboidal, and cylindrical cases, respectively. The use of a convex envelope allows to reduce the size of the resulting global system of around 40% as compared with a canonical truncation. This results into a significant memory footprint reduction of 60% (resp. 50%) in comparison to the cylindrical (resp. cuboidal) case. Let us note that, the layer where generated by Note: Model characteristics for truncated domains with different shapes. The distance between the interface and the scatterer is equal to 2h and an extruded layer with N pml = 5 is generated. The volume of the domain corresponds to the sum of the volumes of the mesh cells.
extrusion, even at the corners, using the procedure described in Section 3.1. With the standard Cartesian and cylindrical PML approaches, where the complex stretch is performed in orthogonal directions, the number of elements and the number of degrees of freedom would be larger at the corners, resulting in a larger computational cost.
In conclusion, it is important to limit the volume of the computational domain and the thickness of the layer to control the computational cost, especially in three dimensions. This advocates using the AML implementation of the PML, which brings geometric flexibility for the definition of the computational domain, and which allows to use efficient unbounded absorbing functions, such as the hyperbolic one.

CONCLUSION
In this article, a comprehensive strategy to implement the PML is proposed for the accurate finite element solution of scattering problems with generally shaped convex computational domains. The approach, which we call the AML implementation, is fully automatic: the end-user only needs to generate the mesh of the main domain, and to choose the layer thickness. The layer does not have to be modeled with a CAD tool, and no geometric information on the layer is required a priori. The approach is very versatile. It can be used with high-order finite element schemes, with unbounded absorbing functions, and for convex domains with regular and nonregular borders. The AML implementation relies on a mesh extrusion generating the layer. Geometric data recorded during the extrusion step are interpolated on the finite element space, providing the required geometric parameters. No further geometric data are required. A rather simple modification of the Jacobian matrix in the elementwise finite element integrals is performed to progressively attenuate the solution inside the layer.
For computational domains with regular borders, the AML can be seen as a specific implementation of the conformal PML. With the conformal PML, a complex coordinate stretch is performed in a local coordinate system associated with the border of the domain. The direct implementation is rather easy thanks to complex material parameters. Unfortunately, the definition of these material parameters required the explicit knowledge of geometric parameters at all points of the layer (i.e., the distance to the border of the domain, the stretching direction, the principal curvatures, and the principal directions), which can be difficult to obtain in practical cases (e.g., with complicated geometries or when only a mesh of the domain is available). With the AML approach, rather than using empirical approaches to estimate the geometric parameters, the complex coordinate stretch is directly performed in an interpolated extrusion direction. This simple strategy allows to alleviate the issue of the accurate estimation of the geometric parameters. Numerical results have shown that the AML is equivalent to the direct implementation of the conformal PML in terms of accuracy for academic settings where the geometric parameters are known.
The proposed approach can deal with computational domains having corners, where the geometric parameters are normally not well defined. The interpolated geometric parameters are then used at the corners as an empirical strategy. The numerical results have shown that the method is very accurate for corners with obtuse angles, and that the accuracy decreases with more acute angles. For settings with right angles, the AML is not equivalent to the Cartesian PML, but the error remains at an acceptable level for the considered cases.
In term of easiness of implementation, the AML has been compared with the direct PML implementation and the implementation based on complex mesh nodes proposed by Ozgun and Kuzuoglu. 13,38 Contrarily to the direct PML implementation, the AML approach does not require the explicit definition of the principal curvatures and the principal directions, though they could be estimated with the interpolated geometric parameters. The minimal geometric knowledge is also a feature of the complex mesh nodes implementation, but that approach cannot be applied with unbounded absorbing functions, which is a severe limitation.
The AML implementation of the PML, combined with the hyperbolic absorbing function, is an accurate, automatic, and computationally efficient approach to solve large-scale scattering problems. The geometric flexibility allows for the use of a convex hull to minimize the volume of the computational domain, thereby reducing the size of the problem and the computational cost.

APPENDIX. GEOMETRIC DATA AND JACOBIAN MATRICES FOR 2D REFERENCE CASES
In this appendix, geometric data and Jacobian matrices used for the direct implementation of the conformal PML are provided. In the general two-dimensional case, the Jacobian matrix of the PML transformation, its inverse and its determinant read, respectively, where e 1 and e 2 constitute an orthonormal frame, and s 1 and s 2 are the stretching functions.
and the stretching functions are with the polar coordinates (r, ) and the stretched coordinate 1 = r − R.
For an elliptical domain of semiaxes a x and a y , the Cartesian components of the orthonormal directions are where p = [p x , p y ] is the closest point belonging to the boundary of the domain. One has p 2 x ∕a 2 x + p 2 y ∕a 2 y = 1. The stretching functions are with the radius of curvature a x a y (A7) and the stretched coordinate = dist(x, p).