A macro‐element strategy based upon spectral finite elements and mortar elements for transient wave propagation modeling. Application to ultrasonic testing of laminate composite materials

Proposing efficient numerical modeling tools for high‐frequency wave propagation in realistic configurations, such as the one appearing in ultrasonic testing experiments, is a major challenge, especially in the perspective of inversion loops or parametric studies. We propose a numerical methodology addressing this challenge and based upon the combination of the spectral finite element method and the mortar element method. From a prior decomposition of the scene of interest into “macro‐elements,” we show how one can improve the performances of the standard finite element procedures in terms of memory footprint and computational load. Additionally, using this decomposition, we are able to efficiently reconstruct important modeling features on‐the‐fly, such as orientations of anisotropic materials or splitting directions of perfectly matched layers formulations, altogether in a robust and efficient manner. We believe that this strategy is particularly suitable for parametric studies and sensitivity analysis. We illustrate our strategy by simulating the propagation of an ultrasonic wave into an immersed and curved anisotropic laminate 3D specimen flawed with an internal circular delamination of varying size, thus showing the efficiency and the robustness of our approach.

of direct simulations that are generally embedded within an inversion algorithm. Hence, proposing useful modeling tools should be performed following a primary goals: (G1) The forward solver should be able to efficiently handle parametric variations of the configuration of interest.
In the scope of this communication, the presented work is driven by a second major goal, which results from our wish to make the modeling tools accessible by a nonspecialist end user, eg, through a commercial software such as the CIVA platform. 12 (G2) No specific hardware architecture should be required from the end user in order to support the simulation's runs.
In the context of high-frequency wave propagation, apart from dedicated analytic solvers as proposed, eg, in the work of Shen and Giurgiutiu, 13 one can consider without being fully comprehensive at least two principal families of generic modeling strategies. First, the class of asymptotic methods 14,15 relying on an asymptotic expansion of the wave field w.r.t. the frequency is widely spread in the NDT 16,17 and seismic waves 18,19 modeling communities. This type of methods provides an efficient and meaningful way to represent the propagation of waves over large distances. However, due to the truncation of the asymptotic development, they fail to represent various, and potentially important, phenomena such as diffraction by small flaws or reflections at critical angles. Second, a full-wave modeling of the propagation phenomena can be obtained using numerical methods, such as finite element methods, 20,21 finite difference methods, 22 boundary element methods, 23 or more recent developments 24,25 based upon these standard strategies. Nevertheless, as they are based upon discretization steps related to the wavelength scale, they often require dedicated solvers based upon specific hardware architectures 26,27,28,29,30 in order to achieve sufficient performances.
In our work, we propose a UT modeling numerical tool based upon a "macro-element" strategy that satisfies both our primary goals (G1) and (G2). This strategy relies on the assumption that the geometry associated to the UT configuration can be decomposed into subdomains. Each subdomain is defined as the transformation of a reference macro-element, the unit cube in 3D, and is assigned to a unique formulation (eg, acoustics, elastodynamics, or absorbing layers, with potentially inhomogeneous material properties). In order to provide an efficient fully discrete propagator, each macro-element is subdiscretized, depending on the estimated wavelengths of interest, using the spectral finite element method. 20,21,31 The communication between the various formulations assigned to the macro-elements is handled using the mortar element method. 32,33,34 In this context, a first novelty presented in this communication consists in making the most of this macro-element strategy so that (1) we significantly increase the performances of standard finite element procedures in both memory and CPU load; (2) we are able to conduct an efficient and lightweight on-the-fly reconstruction of relevant modeling features, such as material anisotropic orientations. Furthermore, since a core challenge of this strategy is to handle the communication between macro-elements, we propose a specific coupling formalism based upon the mortar element method leading to a fully explicit global scheme.
The outline of this paper is as follows. In Section 2, we recall the main components of the discrete propagators built upon spectral finite elements for each formulation of interest. This enables us to detail in Section 3 the macro-element strategy and how it addresses goals (G1) and (G2). In Section 4, we show how the mortar element method can be used in practice to efficiently connect neighboring macro-elements with potentially different formulations. In Section 5, we illustrate our strategy in the context of the UT of curved carbon fiber reinforced polymer (CFRP) composite structures in a 3D setting. We finally give some conclusions and perspectives in Section 6.

DISCRETE PROPAGATORS BASED UPON SPECTRAL FINITE ELEMENTS
We consider the propagators for the fluid, solid, and corresponding perfectly matched absorbing layers (PMLs). They form the minimal set of formulations that one needs to handle in order to address most of the UT configurations. In this section, we recall the main components of the finite element method employed to propose discrete propagators for every formulations, in order to clearly list the benefits brought by the macro-element strategy in the next section. For further details on finite elements, we invite readers to refer to other works 20,21,35,36 for a comprehensive presentation in the case of acoustic and elastodynamics and to the works of Bécache et al, 37 Joly, 38 and Demaldent and Imperiale 39 for the PML formulations.

Fluid, solid, and PML formulations
In the following, we consider Ω ⊂ R d , with d = 2 or d = 3, a bounded computational domain bearing one of the possible formulations. We represent by Ω its boundary and by n the corresponding exterior normal vector field. We denote the coordinates in Ω by x = (x 1 … x d ) ⊺ ∈ R d , and by , the gradient of a scalar and vector field, respectively. Accordingly, the divergence operators of vector and tensor fields are denoted by so that the scalar Laplacian operator simply reads Δ x u = ∇ x · (∇ x u). In the case of a fluid domain, the material is characterized by its mass density , assumed constant for simplicity, and its sound velocity c. The wave equation satisfied by the acoustic pressure p reads along with given initial conditions. The weak formulation associated to (1) consists in finding p ∈ V = H 1 (Ω) for any time t > 0, such that, for any test function p * ∈ V, where m(·, ·) and k(·, ·) are two bilinear forms defined by In the case of a solid domain, the displacement field satisfies the following field equation: with being the mass density, and being the stress tensor field. We complete (4) with relevant initial conditions. The linearized Green-Lagrange tensor and the stress tensor are linked through a linear constitutive law We assume that the fourth-order tensor  can be decomposed into a constant tensor  * defined in a local orthonormal basis so that  satisfies the standard symmetry and positivity conditions. The corresponding weak formulation is similar to (2) where, for any test function * ∈ V = [H 1 (Ω)] d , the bilinear forms are In addition to these natural formulations, we are interested in formulations addressing the challenge of numerically representing infinite propagation areas. Building such formulations is an active field of research that goes beyond the scope of this paper. Hence, without going into details, we use as is a class of absorbing layers referred to as PMLs and, more specifically, the formulation proposed in the work of Demaldent and Imperiale. 39 In the acoustic case, it is a first-order in time formulation based upon split pressure variables {p i } d i=1 and a velocity variable v, satisfying The formulation (7) is completed with relevant initial conditions and boundary conditions. In (7), are the absorbing coefficients. In practice, these directions can be the canonical basis of R d or potentially another set of varying orthonormal directions. In the former case, we refer to the work of Demaldent and Imperiale 39 for the discussion of the "perfectly matched" condition of this formulation. We also refer to the aforementioned work 39 for the expression of the corresponding weak formulation. For solid domains, we consider the similar split first-order formulation. In the scope of this paper, we cast aside the issue of PML stability, and we assume that the material properties  * in (5) allow for a stable use of this formulation. Readers may refer to the works of Bécache et al 37 and Joly 38 and references therein for more details on this matter.

Generalities on spectral finite elements
In this section, we propose to go through the key notions of the finite element discretization. To simplify the presentation, we consider mainly the scalar case and systems of the form of (2). This will enable us to introduce the principal notations in order to detail the benefits of the macro-element strategy later on.

Fully discrete scheme
We consider a Galerkin approximation space V h ⊂ V generated by the basis functions We define the matrix form of the discrete (in space) counter part of (2) as In (8), − → P ∈ R N h is the vector regrouping the coefficients of the discrete solution on the basis of V h . The matrices M and K are the so-called mass and stiffness matrices defined as The fully discrete scheme is obtained by employing an explicit second-order leapfrog time scheme which is stable upon the following Courant-Friedrichs-Lewy (CFL) condition (see, eg, the work of Joly 40 or Appendix A) where r(·) represents the spectral radius of a matrix. Note that this discretization procedure trivially extends to the elastodynamic case with relevant changes in the bilinear forms defined in (6). For PML formulations, we also use a time scheme consistent of order 2, such as in the work of Demaldent and Imperiale. 39

Spectral finite elements and mass lumping
Lagrange finite element methods 35 propose a specific construction of the approximation space V h based upon a given mesh  h of the computational domain. In the following, we consider meshes of hexahedral (or quadrilateral) elements satisfying the standard conformity constraints. The space V h is defined as the space of functions that are globally continuous and that have a local, ie, per element, polynomial representative in a reference element. More precisely, assuming that each element K ∈  h is the result of a diffeomorphism F K applied onto a reference cube (or square)K, we have, for any v h ∈ V h , In (11),  k is the space of polynomials expressed as the tensor product of one-dimensional polynomial spaces, and we denote byn h its dimension We define on the reference element a set ofn h nodesΞ = {̂i}n h i=1 , and we denote by {̂i}n h i=1 the associated local Lagrange polynomial basis. From the transformation of these local nodes by every In (12), we have implicitly introduced the "local-to-global" index mapping Note that, in (11), we allow functions in V h to be locally represented by polynomials with orders potentially different in every directions. This particularity will be referred to as "anisotropic" orders of approximation. We will see in Section 5 that anisotropic orders enables a local adaptation of the discretization, in order to efficiently take into account thin layers of materials.
In the context of hexahedral (or quadrilateral) elements, a convenient way to define a proper set of Lagrange nodeŝ Ξ is by tensor product of d one-dimensional point distributions. The spectral finite element method is based upon Gauss-Lobatto points for defining these one-dimensional point distributions. Due to their optimal convergence, as the order of approximation increases, 31,41 they have received a significant amount of interest in numerous fields of application. In the context of transient wave propagation modeling, spectral elements are particularly interesting since they allow for consistent "mass lumping" at any order of approximation. 20,21,36,42 Broadly, this approach aims at approximating the mass matrix in (9) by a diagonal one, in order to obtain a fully explicit scheme. This approximation of the mass matrix is performed using a specific quadrature formula for computing the integrals appearing in the local mass matrices. The main conditions to achieve consistent mass lumping are the following 21,36,43 : i. the quadrature points and the nodes must coincide to obtain diagonal local mass matrices; ii. for stability reasons, the quadrature weights need to be strictly positive; iii. the quadrature formula must be exact at least for polynomials of order 2k − 2, in order to be as consistent as the case of an exact integration.
For hexahedral (or quadrilateral) elements, these three specific conditions are only satisfied by the Gauss-Lobatto quadrature formula. Spectral finite elements and mass lumping extend to elastodynamics 20,36,44 and PML formulations. 37,39

Local stiffness operations
In practice, the diagonal mass matrix is assembled and stored in an initializing step of the numerical scheme. Hence, the main computational load in (9) comes from the application, at each time step, of the stiffness matrix to a finite element vector. Traditionally, the application of the global bilinear form is decomposed into the application of element-wise forms. Upon the assumption of a constant mass density of the fluid, we have After a change of variable and using the local polynomial representatives of v h and w h , as in (11), the local stiffness operator readsk where J K = | det(∇xF K )|. For conciseness, we introduce the notation where co(·) is the cofactor matrix. As for the mass matrix, we use a Gauss-Lobatto quadrature formula to approximate the integral in (14). We denote by { q }n q q=1 and {x q }n q q=1 ⊂K the quadrature weights and points. The local stiffness operators (14) are approximated byk the expression of the corresponding gradient at the quadrature points reads Introducing the matrix ∈  d×n h ,n h (R) such that the local finite element vectors representingv K and its gradient are given by Regrouping in a matrix  K ∈  d×n h ,d×n h (R) the evaluation of (15) at every local nodes times the quadrature weights, the local stiffness matrix (16) readsk This decomposition extends to the elastic stiffness operator (6) and, for the sake of conciseness, is left to the reader. Similar arguments can also be used in the context of PML formulations in order to derive factorized forms of the bilinear operators appearing in the weak form of (7).

Primary division of the computational domain
In the following, we consider a more complex configuration where, within the bounded domain Ω, multiple formulations (fluid, solid, or PML) can be considered. We assume that the computational domain is decomposed into N subdomains, and that each subdomain M i , referred to as a macro-element, is built from a polynomial transformation i of a reference cube (or square)M, This division of Ω is driven by the following constraints: • the macromesh is a conform subdivision of Ω; • to each interface between two different formulations corresponds an interface, or a set of interfaces, between two MEs or more; • the number of macro-elements should be minimal.
In Figure 1, we illustrate a macro-element decomposition of a curved plate having a through-hole at its center. In this example, the macromesh is made of twelve MEs, deformed from the unit cube in order to fit the geometry of the hole.
In some cases, we may consider an additional set of polynomial transformations i of the reference macro-element. These transformations can be used to represent the computational domain prior to some global deformation Ψ. Let us denote by Ω * the domain in this prior state such that In Figure 1, we illustrate the case where Ω * represents the plane geometry prior to a simple cylindrical deformation. In concrete applications, we might be interested in global deformations Ψ that can be a combination of cylindrical deformations, a spline deformation, or potentially more complex alterations of Ω * emanating from the mechanical response of the specimen to a potential loading. The effect of this deformation onto a macro-element M * i can be represented by the pair of polynomial transformations i and , and we will see in Section 3.4 how we can use this representation to incorporate important modeling features into our numerical strategy, such as the evolution of the directions of anisotropy after deformation of the material.

Subdivision of the macro-elements
We seek, for performance purposes, a macro-element decomposition with the lowest number of subdomains. This prior division of the computational domain is unfit for adequate numerical procedures at the wavelength scale, and we resort to an internal subdivision of each macro-element M i . It corresponds to a grid of hexahedral (or quadrilateral) cells deformed by i . Let us denote by h,i and  h,i the grid prior to and after its deformation so that, formally, Each grid cellK ∈ h,i is obtained from the reference element by a simple affine transformation FK. For instance, for d = 3, it reads where {O p (K)} d p=1 and {h p (K)} d p=1 are the origin and the lengths of the cell, respectively. With these notations, for each element K ∈  h,i , the transformation F K introduced in Section 2.2.2 is expressed as To each element in  h,i , we associate a spectral finite element with anisotropic orders. We denote by N h,i the number of nodes and by the coordinates of the nodes. Note that (18) and (20) naturally extend to the polynomial transformation i of the macro-element so that we can define a mesh  * h,i of M * i , composed of the nodes Ξ * i . In practice, the refinement of the grid depends, for precision purposes, on an a priori estimation of the wavelengths of interest and on the maximal deformation of the grid. Note that, in (19), the lengths may vary from one cell to another, as long as h,i remains a conform grid. This will be of major importance when dealing with stratified materials in Section 5.

Assumption of conform interfaces
Throughout this paper, we make the assumption that the union of every macro-element subdivisions generates a conform mesh of the computational domain after discarding redundant coordinates at interfaces. In particular, let M i and M j be two macro-elements having a face in common, say Γ, then there exists a one-to-one mapping between the coordinates of the nodes Ξ i | Γ and Ξ j | Γ . Let us denote by M h the number of nodes on the interface Γ, we can introduce the "interface-to-volume" index mappings Alternatively, these mappings may be represented by the rectangular matrices This aspect, referred to in the following as the conform interface assumption, will have important consequences when using the mortar element method in Section 4.

Benefits of the macro-element strategy
In order to cope with our objective (G2), we favor low-memory strategies. In this regard, we consider "unassembled" operations for representing the application of the stiffness matrix to an input finite element vector. In essence, it corresponds to applying every local stiffness matrices, described in Section 2.2.3, to local input vectors. In the simplified configuration where the subdomains bear an acoustic formulation, we propose the pseudocode in Algorithm 1 for performing these operations. In this pseudocode, are the output vectors storing the results of the application of the stiffness matrices. We also consider in Algorithm 1 a standard coloring of each element in  h,i so that no elements in a color group have nodes in common.
These unassembled operations are rather natural in the context of high-frequency wave propagation modeling. However, our approach enables some significant improvements, mainly due to the fact that the underlying subdiscretization of each macro-element is related to a reference grid inM. In particular, while in a general case of unstructured meshes, the "local-to-global" mapping G defined in (13) is usually stored, it can be recomputed at no expense within each macro-element, thus sparing important memory load. Additionally, the matrix  K in (15), built from the transformation of the current element, can also be recomputed on-the-fly. To do so, we simply apply the transformation (20) to the local nodes coordinatesΞ and compute the gradient of the transformation using ⊺ , defined in (17). Thus, most of the memory footprint of a macro-element comes from the three finite element vectors of the time scheme and the diagonal mass matrix. Note that in this context, the major part of the workload in Algorithm 1 is centered on computing gradients of local finite element vectors. In practice, we use a specific data structure in order to improve this operation (readers may refer to the work of Carrascal-Manzanares et al 45 for more details). On top of improving the memory load, one can also increase the CPU performance by applying the local stiffness operations in parallel. Indeed, due to the definition of the elements coloring, the third loop in Algorithm 1 is naturally fit for parallelism. Since the topology of the mesh of a single macro-element is essentially a grid, this coloring is trivial and optimal in the sense that (1) each color group has the same number of elements and (2) the number of different color is minimal.
Another interesting aspect of the macro-element strategy resides in its capacity to reconstruct efficiently and on-the-fly specific local vector fields. An important example is when a constant unitary vector, say v * , is defined on a macro-element M * i prior to its deformation by Ψ. In this context, we can retrieve its local variation within the final macro-element through A natural application is the case of fiber materials, where specific material properties are defined in a privileged direction, the fiber, while being isotropic in the orthogonal plane. One may typically define the constant fiber orientation in Ω * and compute the result of its transformation from (23). From the local fiber orientation, we can reconstruct any two orthogonal directions and define a local isotropic transverse constitutive law using (5). In Section 5, we present numerical results using this computational strategy to model curved laminate composites. A similar configuration is when a constant unitary vector, sayṽ, is defined in the reference macro-elementM. We can retrieve, within each macro-element, its local deformation from This can typically be used to compute splitting directions of PML systems of the form of (7). Operating in this fashion, we encompass in the same framework curved and straight PML domains with an automatic and lightweight definition of the local directions. Low memory footprint, parallel operations on standard CPU, and automatic reconstruction of relevant modeling components form the main benefits of our approach that enable us to address our primary goals (G1) and (G2). These advantages can be significant especially for 3D configurations, but they require a "minimal" high-order hexahedral (or quadrilateral) decomposition of Ω, which is known to be a major challenge in the field of computational geometry. However, since our goal is to achieve parametric studies, we do not intend to provide efficient numerical procedures for fully generic configurations. In most of the practical cases, we restrict ourselves to parametric geometries that can be obtained from case-dependent meshing procedures. Readers may refer to other works 46,47,48,49 for some examples in the context of UT modeling.

COUPLING MACRO-ELEMENT FORMULATIONS USING MORTAR ELEMENTS
One of the challenge of our approach is to enable the communication between neighboring subdomains. Upon the conform interface assumption stated in Section 3.3, the difficulty essentially lies in coupling different formulations within the same framework. To address this issue, we rely on the mortar element method. 50,51,52,53 This domain decomposition method has been successfully used as a mean of incorporating different and independent space discretizations within a global numerical scheme. A significant amount of research (see, for instance, other works 32,33,34,54,55 ) has been dedicated to the analysis and application of mortar elements for elastodynamics and wave propagation problems. In our work, we use this method for our specific objectives, and we show how one can minimize its computational cost by employing a lumping technique, as the one described in Section 2.2.2.

Schur complement
To start with, we recall the Schur complement method, which will be used to devise a global algorithm for coupling macro-element formulations. We consider two vectors of unknown: − → X ∈ R N , referred to as the volume unknown, and − → L ∈ R M , referred to as the interface unknown. We assume that they satisfy the following linear system: where − → F and − → G are given right-hand sides. The matrices D and D , of dimension N × N and M × M respectively, are diagonal. The matrices B and C of dimension N × M may be referred to as transmission (from interface to volume) matrices. The Schur complement method for solving this type of linear system reads: Step 1. Preprocessing. We compute an auxiliary variable Step 2. Computation. Defining the Schur complement matrix and assuming that this matrix is invertible, we solve the following linear system Step 3. Postprocessing. The volume unknown is obtained by Note that, in these three steps, the main issues are the invertibility of the Schur complement matrix S in (27), and, if it is invertible, the efficient computation of the interface unknown solution of (28).

Two-domain problems
To start with, we consider the case of two subdomains, ie, Ω = M 1 ∪ M 2 , connected at an interface Γ satisfying the conform interface assumption, detailed in Section 3.3. We denote by the normal vector field of the interface oriented from M 1 to M 2 .

The illustrative example of fluid-fluid coupling
We assume that both M 1 and M 2 support a fluid formulation (1), satisfied by the corresponding pressure fields p 1 and p 2 with mass densities 1 and 2 and sound velocities c 1 and c 2 . In this case, the mortar element method introduces a Lagrange multiplier ∈ W = H − 1 2 (Γ), satisfying 56 and weakly imposes the continuity condition on p 1 and p 2 at the interface. We define, for i = 1, 2, the space V i = H 1 (M i ), and the product space X = V 1 × V 2 × W. The weak formulation of the coupling scheme aims at finding the solution (p 1 , p 2 , ) ∈ X for a time t > 0 such that, for any (p * 1 , p * 2 , * ) ∈ X, we have In (30), we have assumed free-surface boundary conditions on the remaining faces of the macro-elements. In addition, m i (·, ·) and k i (·, ·) are the bilinear forms defined as in (3), with corresponding material properties i and c i . Note that, in this formulation, we have introduced in the weak continuity constraint a penalization term proportional to a (small) positive constant coefficient ≥ 0 and the second time-derivative of the Lagrange multiplier. This specific choice of the penalization term is justified in Appendix A, where we consider the fully discrete energy norm of the coupling system. We define V h,i ⊂ V i the finite element approximation spaces, detailed in Section 2.2.2, and we denote by the associated basis functions. The approximation space of W is defined by taking the trace space of V h,1 onto the interface, namely, As a consequence, there is a trivial mapping between the basis function of W h and the basis functions of V h,1 . Upon the conform interface assumption, this trivial mapping extends to basis functions of V h,2 , and using the "interface-to-volume" mapping (21), we have In the context of the mortar element method, our choice of discrete space for the Lagrange multipliers is nonstandard and usually prohibited, since it may lead, in general configurations, to an unsatisfied discrete inf-sup condition. 57,58 It is for this particular reason that the penalization term was introduced in (30), in order to retrieve the stability of the discrete coupling scheme. Applying a second-order time-scheme, such as the one presented in (9), and a centered scheme for the coupling terms lead to the following fully discrete coupling scheme In (32) From the choice of the discrete space for the Lagrange multipliers, we have where  1 Γ (·) is the "interface-to-volume" index mapping defined in (21). Hence, applying the same consistent mass lumping technique as for the volume mass matrix, presented in Section 2.2.2, we obtain a diagonal interface mass matrix. The transmission matrices C i ∈  N h,i ,M h (R) are expressed as and, using the matrix representation of the "interface-to-volume" index mappings (22), can be expressed as Remark 1. In Appendix A, we detail the arguments proving that, in order for the discrete system (32) to be stable, the time step must respect the CFL condition In other words, the introduction of the Lagrange multipliers does not modify the stability condition on the time step deriving from both subdomains, which is a major asset of this domain decomposition method.
The discrete scheme (32) can be expressed as where the right-hand sides are System (37) is a specific case of system (25) with unknowns and matrices Hence, the Schur complement matrix reads Since the interface mass matrix and the volume mass matrices are diagonal, using the decomposition (35) of the transmission matrices, and remarking that is the identity matrix, one can verify that the Schur complement matrix is diagonal This lumped Schur complement matrix is invertible for ≥ 0, and its memory footprint is very limited. Hence, the coupling procedures represent little computational overheads.

The case of fluid-solid coupling
We now consider the case where M 1 supports a fluid formulation (1) and M 2 supports a solid formulation (4). We introduce two Lagrange multipliers v and s used to solve the normal velocity constraint and the normal stress constraint, where v 1 is the velocity associated to the fluid acoustic pressure in M 1 . The first relation can be written in terms of the pressure unknown since In (41), m 2 (·, ·) and k 2 (·, ·) are the bilinear forms defined in (6). Compared to (30), no penalization term is needed since the normal velocity and normal stress weak continuity relations directly involve coercive terms w.r.t. the Lagrange multipliers. We define the discrete space W h as in (31). Using an order 2 centered time scheme for the complementary terms appearing in (41), we obtain the following fully discrete scheme: In (42), the matrices M Γ and C 1 are defined in (33) and (34), respectively. The matrix C 2, ∈  N h,2 ,M h (R) is defined as are the (vectorial) basis functions generating V h,2 . Using similar energy arguments as in Appendix A, one can prove that the CFL condition on the time step, for the discrete system (42) to be stable, is the lowest condition of the two subdomains computed independently. The numerical scheme (42) is equivalent to where the right-hand sides related to the volume unknowns are expressed as while the right-hand sides associated to the discrete Lagrange multipliers are System (43) is a specific case of system (25) with unknowns .
The diagonal and extradiagonal block matrices are ) . .
By renumbering the unknown vector − → L in the following way and applying the same operation in the Schur complement matrix, we observe that S is block diagonal with each block defined as Note that, to obtain (44), we used the fact that the interface mass matrix is diagonal and the expression (35) of the transmission matrices. The determinant of each block is strictly positive, thus the computation step (28) is well-posed. In practice, we store the inverse of each 2 × 2 block in an initializing step, and we apply each inverse in parallel in order to solve (28) at each time step, thus limiting the cost overhead of the coupling procedures.

Coupling with PML formulations
Coupling with PML formulations can be handled in the same way. To illustrate this case, we suppose that M 1 holds an acoustic formulation and M 2 an acoustic PML formulation such as the one given in (7). We introduce the Lagrange multiplier The first relation can be expressed in terms of pressure variables as in Section 4.2.2. Assuming that where is the set of splitting directions in (7), then we can express the second relation in the equivalent form Note that the assumption (45) is satisfied in most practical cases. Using the relevant Green formula in both systems, we can introduce the Lagrange multiplier in both weak forms and retrieve the expression provided in the work of Demaldent and Imperiale. 39 After time and space discretization, we obtain the following Schur complement matrix are the (constant) absorption coefficients in each splitting direction. Readers can refer to the aforementioned work 39 for the details of the numerical scheme, not recalled here for the sake of conciseness. Note that, as for the previous case (40), we obtain, in the conform case, a lumped Schur complement matrix.

N-domain problems
So far, we have considered two-domain problems exclusively and we have shown how, by using the natural trace of the volume discrete space, we can obtain lumped Schur complement matrices leading to very efficient coupling between subdomains. Traditionally, this specific choice of discrete space for the Lagrange multipliers is prohibited since it may entail ill-posed Schur complement matrices, even in the conform case. To circumvent this difficulty, numerous possibilities of adequate discrete spaces for the multipliers have been proposed (see, eg, other works 33,34,51,52 ). In our work, we rely on an alternate method based upon a penalization strategy, suggested in the work of Brezzi and Fortin. 58 This enables us to recover the invertibility of the Schur complement matrices, while keeping the computational performances of lumped matrices.

Canonical example of an ill-posed lumped Schur complement matrix
To illustrate the loss of invertibility of the Schur complement matrix, we consider a canonical example of four subdomains put together in a square-shape fashion, as depicted in Figure 2. In this example, the computational domain is composed of four macro-elements We denote by = {a, b, c, d} a generic interface index and by i, j ∈ {1, 2, 3, 4} two subdomain indexes. To each interface Γ , we associate a normal vector field , oriented as depicted in Figure 2. We assume that each subdomain bears an acoustic formulation with identical mass densities and sound velocities. The Lagrange multiplier associated to each interface is defined as in (29). Additionally, we assume that the macro-elements are of the same size and that their corresponding meshes have the same number of elements and the same (isotropic) order of approximation. Similarly to (37), we can write in this context a system of the form of (25). The volume and interfaces unknowns are The matrix D is the concatenation of four identical diagonal mass matrices M, and D is the concatenation of four identical interface penalization matrices M Γ . We denote by W h, the discrete space for the Lagrange multipliers associated to the interface Γ , such that The transmission matrices are where ∀i = 1, … , 4 and ∀ ∈ {a, b, c, d} Note that, in (46), the sign of the block matrices can be directly deduced from the chosen orientation of the normal at each interface. The Schur complement matrix appearing in the computation step (28) is symmetric and reads where we denote by (·) the nonzero extradiagonal blocks deduced from symmetry and by S ,i the diagonal blocks such that For a fixed i = 1, … , 4 and ( , ) ∈ {a, b, c, d} 2 such that ≠ , blocks of the form C ⊺ i M −1 i C i represent the interactions between the Lagrange multipliers associated to the interfaces Γ and Γ through the macro-element M i . Hence, from the decomposition (35), one can see that these extradiagonal blocks have only one nonzero value at the cross-point. Assuming an identical numbering of the interfaces going from the cross-point to the free extremity, we denote by the index of the node concerned with the cross-point and the rest of interface indexes in Γ , respectively. Following this notation, we propose the renumbering of the interface unknown so that the Schur complement matrix takes the form of where S × is the 4 × 4 cross-point matrix expressed as The expression of each extradiagonal value of this cross-point matrix is In the specific case of our canonical example, the mass and interface mass matrices are identical. We denote by m × and m × Γ their corresponding values at the cross-point index. After factorization, the cross-point matrix takes the simpler form of and has a nonempty kernel generated by the constant vectors Thus, we can see the necessity of the penalization strategy since the Schur complement matrix S is ill-posed for = 0. It should be noted that the simple form (50) of the kernel of the cross-point matrix is linked to the specific parameters of our canonical configuration. In more generic cases, we expect this kernel to vary depending on the discretization patterns and the material properties of the macro-elements.

Schur complement matrix for N-domain problems
We consider a specific data structure, referred to as the "skeleton" of the macromesh, regrouping the information of the incident faces (or edges) at every cross points, and the incident macro-elements at every interfaces. Using this "skeleton," we can assemble the set of cross-point matrices, denoted by {S × r } N × r=1 , where N × is the number of cross-points in the macromesh. If d = 2, the cross-points are restricted to points with multiple incident edges, while, if d = 3, the cross-points are spread across edges with multiple incident faces. The form of the cross-point matrices is similar to (49) with dimensions equal to the number of incident faces (or edges). We also assemble the Schur complement matrices set on the interior nodes. Let  be the set of interface indexes, such that for any ∈ , the interface Γ is shared by two macro-elements M i( ) and M j( ) . These interior Schur complement matrices are denoted by {S • i( ) ( ), } ∈ , with expressions deduced from (47). Finally, generalizing the previous renumbering (48), we can express the Schur complement matrix for N-domain problems in the general form of .
In practice, the interior Schur matrices are diagonal and every cross-point matrices are independent matrices of reasonable size so that we can compute and store their inverse in an initializing step. Applying their inverse at every time steps is performed in parallel, leading to a limited cost overhead.
Remark 2. The penalization procedure is a compromise between stability and consistency. In order to limit spurious consistency errors, we set = (Δt 3 ) which is below the consistency errors of the numerical time schemes used for the time-discretization of the various formulations. Additionally, since only the cross-point matrices are ill-posed, we use an "inhomogeneous" penalization parameter set to zero at the interior nodes of the interfaces.

NUMERICAL SIMULATION OF IMMERSED CURVED CFRP IN 3D
As a numerical illustration, we propose to model the UT of a 4 mm thick immersed laminate composite material. The specimen is composed of 16 isotropic transverse plies of 235 m thickness. The anisotropy direction of the plies, representing the fiber orientation, changes from 0 = (1 0 0) ⊺ to 90 = (0 1 0) ⊺ successively from one ply to another. Between each ply and on the top face of the structure, we consider thin intermediate epoxy layers of 15 m represented as isotropic materials. Overall, the stratification is made of 32 layers. Using Voigt notation to represent the constitutive law  * in (5), the material properties associated to each ply and epoxy layers are presented in Table 1. The surrounding fluid is water with a mass density of w = 1.0 g · cm −3 and a sound velocity of c w = 1.483 mm· s −1 .
We assume that the structure is subject to a cylindrical deformation with respect to the y-axis with a curvature center positioned at 10 mm below the specimen. We consider an incident pressure field p inc as a plane wave with a spatial Gaussian window where x 0 = (0 0 3) ⊺ is the incident wave origin, d = (0 0 −1) ⊺ is the direction of the propagation, the standard deviation of the Gaussian window is set to G = 1.0, and we define the excitation signal with F = 3 MHz, t 0 = 0.75 s, and f = 0.2. In the following subsections, we use our numerical model to compute the total field, ie, the sum of the incident field and its interactions with the stratified specimen. Due to the localized spatial

FIGURE 3
Numerical domain for the healthy immersed curved Carbon Fiber Reinforced Polymer. PML, perfectly matched absorbing layer support of the incident field, we truncate the numerical domain in the tangent plane of the specimen so that its lengths are reduced to 6 mm in both directions. To avoid spurious reflections, we surround the solid area with 1-mm solid PMLs. Above and below, we append two 1-mm fluid subdomains surrounded with additional 1-mm thick fluid PMLs. To sum up, prior to its cylindrical deformation, the overall numerical domain has a size of 8 × 8 × 8 mm 3 , and the disposition of the various formulations is depicted in Figure 3. The incident field is evaluated at the upper interface between fluid and fluid PML subdomains.
Remark 3. In Appendix B, we show how the transmission conditions are adequately taken into account in the mortar element method. As a result, above the entry surface in the fluid PMLs, we solve the field equations for the diffracted field and below for the total field. In practice, this enables us to restrict the surrounding fluid subdomains to a limited size while avoiding spurious interferences between the incident and total fields.

Simulation in the case of a healthy curved CFRP
To start with, we consider the case of a healthy specimen without internal flaw, as presented in Figure 3, for which we have a total of 45 macro-elements. Note that, due to the cylindrical transformation, the constitutive law varies locally and is expressed as in (5). In this case, the local fiber orientation of each ply is computed using (24), withṽ being either 0 or 90 . In this configuration, the smallest wavelength is the one related to the fluid material, which is about 0.5 mm. Hence, in the tangent plane, we perform a subdiscretization of each macro-element so that we can insure two elements of order 4 per wavelength. In the thickness of the solid material, since the size of each layer is smaller than half of the wavelength, we use a specific discretization pattern and use one element of order 3 for the each ply and one element of order 2 for each epoxy layer. Note here the benefits of the anisotropic orders in the finite element discretization, depicted in Section 2.2.
In Table 2, we gather the various characteristics of the simulation. We differentiate the number of nodes in the final (subdiscretized) finite element space from the number of degrees of freedom (DoFs). The latter is obtained by taking into account the dimension of the unknown at each node, which is 1 for fluid, 3 for solid, 6 for fluid PML, and 18 for solid PML subdomains. It is important to emphasize that, even though we are facing a reasonably large number of DoFs within a locally varying anisotropic context, the memory footprint of the overall numerical scheme is quite low (less than 1 GB). Additionally, for reaching the maximal time corresponding to the propagation of the ultrasonic beam through the structure, the computational time is about 10 minutes. The simulation was carried out on a biprocessor computer Intel ® Xeon ® Processor E5-2687W v2, 2 × 8 cores at 3.40 GHz. For a more detailed examination of the performances of the proposed numerical solver, including a comparison with a commercial finite element software, readers may refer to the dedicated communication. 59  We present in Figure 4 the snapshots of the solution at various time steps. In these snapshots, we observe the propagation of the primary ultrasonic beam influenced by the anisotropy of the structure and the generation of structural noise coming from the multiple reflections at the layer interfaces. Additionally, we consider an observation point P = (0, 0, 2.95) ⊺ located above the structure, and an observation point Q = (0, 0, −2.95) ⊺ located below the structure (see Figure 3). We plot the pressure field obtained at the observation point P in Figure 5A. We see three main contributions: (P.I) is the incident field defined in (51); (P.II) is the reflection of the incident field at the structure's upper boundary; (P.III) is the backwall echo. The signal between (P.II) and (P.III) is the phenomena referred to as the structural noise. In Figure 5B, we plot the signal at point Q and we observe a main contribution (Q.I), which is the transmitted wave that traveled through the complete structure.

Illustration of a variation of the configuration with a circular delamination flaw
Starting from the configuration described in the previous section, we introduce a circular delamination flaw located at the center of the structure. We add the necessary macro-elements so that the flaw geometry corresponds to an interface between two macro-elements. Between these two MEs, we discard the mortar element coupling scheme in order to incorporate the effect of the flaw in the numerical scheme, which is equivalent to duplicate and detach interfacial nodes. We show this configuration in Figure 6. In this more intricate case, the local anisotropy orientation of each ply is computed using (23), with v * being either 0 or 90 .
In Table 3 Depending on the radius, we adapt the geometries of the macro-elements; hence, the number of nodes varies from one simulation to another. It should be noticed that, due to the inclusion of the flaw geometry and the conformity constraints that need to be satisfied, the final mesh incorporates elements of different sizes. As a consequence, the time step has decreased compared to the healthy case, and the overall computation time has increased in order to reach an identical maximal time of 5.5 s. We present in Figure 7 the snapshots of the solution at the same time steps than in Figure 4. For this configuration, we have fixed the radius of the circular delamination flaw to r = 0.75 mm. We observe the diffracted wave due to the presence of the flaw. This contribution is clearly seen in Figures 8A and 8B where we plot the pressure field at two observations points above and below the flawed structure. The contributions (P.I), (P.II), and (P.IV) are identical to the healthy case, context, this strategy can be used to represent the propagation of the incident field, without having to incorporate the complete propagation area in the numerical domain. When modeling realistic UT configurations, it is also required of the forward solver to represent the electrical signal V R captured by the receiving piezoelectric transducer. In this context, the pioneer works 1,2,60,61 express the received electrical signal as a function of the so-called reciprocity quantity , namely, where T is the maximal time, · ⋆ · is the convolution in time domain, and  is a convolution kernel representing the sensitivity of the transducer acquisition chain. The reciprocity quantity  reads where p tot, E is the total pressure field computed from the incident field generated by the emitting transducer, and p inc, R is a "virtually emitted" field from the receiving transducer. Note that Γ is a boundary surrounding the structure, or the flaw, and in the context of UT of composite materials, we generally choose the entry surface. In practice, the numerical solver is used to compute the quantity p tot, E , as presented in the numerical examples of this section. Hence, by appending the computation of the reciprocity signal to the finite element computations, our approach is compatible with this type of hybrid methods. For more details and illustrations of numerical results obtained using the development version of the CIVA software, 12 readers may refer to prior communications. 46,47,49,62

CONCLUSION AND PERSPECTIVES
In this work, we have proposed a specific numerical solver for wave propagation modeling designed for UT configurations and particularly suited for stratified composite materials. Our approach is focused on addressing, in an efficient and robust manner, parametric variations of the configuration, a fundamental prerequisite for inversion loops or sensibility studies. To reach this objective, we have detailed a macro-element strategy which is based upon a decomposition of the configuration of interest. Each subdomain, or macro-element, in this decomposition is associated to a specific wave propagation formulation, from which we derive a discrete propagator using the spectral finite element method. Making the most of the a priori information embedded in the scene decomposition, we are able to improve the performances of the standard finite element operations in terms of memory footprint and computational load. Additionally, we have shown how we can extract from the macro-elements' deformations important modeling components such as the anisotropy orientations or the splitting directions of PMLs, altogether in a lightweight and efficient fashion. The transmission conditions between each subdomain are solved using the mortar element method, which is sufficiently flexible to take into account, in the same formalism, the various set of formulations arising in numerical UT modeling. In particular, we have depicted how, using the traces of the volume finite element spaces and a penalization strategy, one can apply a lumping integration technique on the Schur complement matrices in order to significantly reduce the cost of the coupling method, in the case of conform interfaces. To illustrate the efficiency of this approach, we have proposed a 3D configuration of a curved composite materials, flawed with an internal circular delamination of varying radius. A first and important perspective of this work would be to incorporate, within the macro-element strategy, the case of nonconform interfaces. Indeed, in numerous UT configurations, the wavelength ratio may vary strongly between different inclusions of materials. A typical example being the case of immersed specimen, as the one proposed in the numerical experiments of this communication. In these configurations, using the ability of the mortar element method to deal with independent subdomain discretizations, we could significantly improve the efficiency of the global numerical scheme. However, the main issue, in this context, is to solve the interface linear system, which cannot be lumped, in an efficient manner. A potential lead to address this difficulty could be to rely on the conformity of the macromesh and to only consider "controlled" nonconformities, where we have an a priori assumption on the refinement ratio between two adjacent discretizations. In addition to spatial nonconformities, we could also consider the class of methods that allows to include, within the same numerical scheme, different time discretizations. We can mention, for instance, the local time-stepping method 63 or the locally implicit method. 64 Another possible way to improve the range of configurations efficiently handled by our approach could be to use the mortar element method to "implicitly" incorporate thin layers of materials, such as the epoxy layers in the presented numerical illustrations. We could, for instance, consider the standard spring-mass model 65 or a more complex asymptotic propagator associated to the thin layer, such as the one proposed in the work of Bonnet et al. 66

ACKNOWLEDGEMENT
Authors would like to thank Dr Sébastien Imperiale of MΞDISIM project-team at Inria Paris-Saclay for the useful insights on the penalized form of the mortar element method.

A.1 Discrete energy norm of second-order schemes
To start with, we recall a standard stability result (see the work of Joly 40 for more details). Considering the second-order fully discrete scheme (9), one can verify that the functional satisfies a specific conservation property. Indeed, multiplying (9) by 1 To turn this conservation property into a stability result, we need to prove that the functional is positive for any solution vector. While it is self-evident for the first term in (A1), the second one needs some further manipulations. Remarking that Hence, upon the CFL condition (10), the functional is positive and represents a discrete energy norm.

A.2 Stability of the fluid-fluid coupling system
We consider the case of system (32). For i = 1, 2, let {E Δt Using the discrete form of the continuity relation in (32), we obtain