Applied Mathematical Implementation and computational aspects of a 3D elastic wave modelling by PUFEM

This paper presents an enriched ﬁnite element model for three dimensional elastic wave problems, in the frequency domain, capable of containing many wavelengths per nodal spacing. This is achieved by applying the plane wave basis decomposition to the three- dimensional (3D) elastic wave equation and expressing the displacement ﬁeld as a sum of both pressure (P) and shear (S) plane waves. The implementation of this model in 3D presents a number of issues in comparison to its 2D counterpart, especially regarding how S-waves are used in the basis at each node and how to choose the balance between P and S-waves in the approximation space. Various proposed techniques that could be used for the selection of wave directions in 3D are also summarised and used. The developed ele- ments allow us to relax the traditional requirement which consists to consider many nodal points per wavelength, used with low order polynomial based ﬁnite elements, and there- fore solve elastic wave problems without reﬁning the mesh of the computational domain at each frequency. The effectiveness of the proposed technique is determined by compar- ing solutions for selected problems with available analytical models or to high resolution numerical results using conventional ﬁnite elements, by considering the effect of the mesh size and the number of enriching 3D plane waves. Both balanced and unbalanced choices of plane wave directions in space on structured mesh grids are investigated for assessing the accuracy and conditioning of this 3D PUFEM model for elastic waves.


Introduction
Growing research activities have been taking place in various wave numerical modelling fields such as acoustics, surface water waves, radar waves, seismology, geophysics and biomedical ultrasound. The related problems in the frequency domain are modelled using mainly the Helmholtz equation [1][2][3] , Maxwell equation [4] and Navier equation [5][6][7][8] depending on the wave propagation medium and the type of application.
A number of different numerical methods have been used to solve such problems and the most commonly used is the finite element method (FEM) due to its flexibility in handling complex geometries and its ability to model different types of media. However, the use of low order polynomial based finite elements for solving wave problems requires fine mesh grids with many nodal points per wavelength. At medium and high frequencies, this leads to large memory requirements and high CPU time. Moreover, because of the so-called pollution error, even finer mesh grids would be required to maintain an acceptable accuracy.
The last two decades saw the development of numerical techniques aiming at modelling wave problems with improved accuracy and reduced computational cost, in comparison to the conventional FEM. These techniques involve the incorporation of a priori knowledge of the wave field in the numerical model. Indeed, the wave field is expanded into sets of analytical functions, for example, in the form of plane waves propagating in different directions. For the case of the Helmholtz equation, various methods have been successful in efficiently solving wave problems with coarse mesh grids in both 2D and 3D with reduced computational effort and better accuracy, in comparison to standard polynomial based FEM. Among these methods are the least-squares method [9] , the partition of unity finite element method (PUFEM) [10][11][12][13][14][15][16] , the ultra weak variational formulation (UWVF) [17,18] , Plane wave discontinuous Galerkin (PWDG) [19] , the generalized finite element method [20] and the discontinuous enrichment method (DEM) [21] . In the case of the boundary element method, a partition of unity boundary element method (PUBEM) has been developed for two and three dimensional Helmholtz problems [15,22] and an isogeometric wave-enriched method within the context of boundary elements (XIBEM) was presented for two and three dimensional Helmholtz problems [23,24] . Other solutions were also considered such as the oscillated finite element polynomials [25] and in a more recent work the stable discontinuous Galerkin method [ 26 ]. The variational theory of complex rays (VTCR) [27] and the phase reduction finite element method (PR-FEM) [28] are other dedicated strategies for solving short wave problems with coarse mesh grids.
In the current work, the PUFEM concept is extended to three-dimensional time-harmonic elastic wave problems. To the authors' knowledge, only the works [39,40] have dealt with three-dimensional plane wave enriched elements for elastic wave problems within DEM and UWVF frameworks, respectively. In both DEM [39] and UWVF [40] the approximating plane waves are defined at the element level and hence inter-element continuity must be enforced in the formulation. However, in the presented PUFEM elastic model, the plane waves defined at the nodes of the elements are multiplied by the Lagrange polynomial shape functions and hence continuity is automatically ensured.
In this model, the components of the displacement field are first expressed as sums of both pressure (P) and shear (S) wave contributions, like in the two-dimensional case [5] . The main difference with the three-dimensional case consists in the direction of S-wave induced displacement, which is orthogonal to the direction of propagation. While in two dimensions the orthogonal direction is well defined via a single vector, in three-dimensions the orthogonal direction is contained within the orthogonal plane to the propagation direction and therefore two vectors are required for its resolution. Furthermore, in two dimensions, the directions of the approximating plane waves are easily defined whether in a uniform or nonuniform distribution. This is not the case in three dimensions where the distribution of plane waves in space, either uniform or nonuniform, is not straightforward.
Indeed, in the implementation of the three-dimensional version of PUFEM the uniform distribution of directions in space is a challenging issue, unlike in the two-dimensional case. In fact, the uniform distribution of a set of points on a sphere has attracted the attention of researchers from various fields including meteorology and quantum theory. In this paper, different approaches for uniform distribution of directions in space are summarised and then some of them are used to obtain sets of linearly independent directions. Moreover, to keep the conditioning within acceptable limits, unbalanced numbers of directions for the pressure and shear waves are used.
The outline of the paper is as follows. In Section 2 , the model problem and the used notations are presented. Section 3 is dedicated to the variational formulation and Section 4 presents the PUFEM approximated solution. In Section 5 , various approaches to uniformly distribute sets of plane waves in space are presented. The PUFEM discrete model is given in Sections 6 and 7 presents a validation of the model against problems with analytical solutions such as progressive plane waves in an elastic infinite medium and elastic wave scattering by a spherical rigid body. For the latter problem, a comparison between PUFEM and the commonly used low order FEM is also attempted. Finally some concluding remarks are drawn in Section 6 .

Problem statement
Let be the spatial domain in R 3 occupied by an elastic medium. Let us denote by ( e 1 , e 2 , e 3 ) the cartesian vector system, and x = x 1 e 1 + x 2 e 2 + x 3 e 3 a generic point in R 3 . The propagation of elastic waves in a 3D homogeneous isotropic medium is governed by the wave equation [41] ρ∂ tt U − ∇ · σ(U ) = ρF , (1) where ρ is the density of the medium and F is a body force. The notation ∂ tt stands for the second partial derivative with respect to the time variable t . The stress tensor σ, evaluated at a displacement U = U 1 e 1 + U 2 e 2 + U 3 e 3 , is given by the classical Hooke's law σ(U ) = λ∇ · U I + μ ∇U + ( ∇U ) , (2) where I is the identity matrix in R 3 ×3 (3 × 3 matrix), λ and μ are the Lamé parameters of the elastic material, assumed constant, and ∇U = ( ∇ U 1 , ∇ U 2 , ∇ U 3 ) . We will denote here by ' ' the transpose of a given vector or tensor. The dot product ' ·' of ∇ with a tensor field ( However, for a vector field in R 3 , ' · ' is the usual dot product. The case of harmonic motion has practical interest. Suppose that F = exp ( i ωt) f , where ω is the circular frequency and i = √ −1 is the imaginary unit number. If U = exp ( i ωt) u is a solution of (1) then the strong form of the problem in terms of u can be rewritten as follows with the following boundary condition where n and t denote, respectively, the outward unit normal and tangent vectors to the boundary ( = ∂ ) and g is the source term. Details on how this boundary condition is deduced are given in reference [31] . Note that (5) is an absorbing boundary condition of the first order for elastic wave equations when g = 0 on . The aim of this work consists to validate the presented numerical model and therefore the source term g is introduced analytically in expression (5) to enable the solution of the problem. It is well known from the elastic wave theory [41] that a displacement field solution of (4) is a superposition of two types of purely oscillatory waves: P-wave (primary, or compressional) and S-wave (secondary, or shear). This can be demonstrated by writing The P and S waves travel with speeds The corresponding compressional wave number k P and shear wave number k S are given by These last two physical quantities are useful for the numerical modelling approach that we will discuss and investigate in this paper.

Variational formulation
Since PUFEM is a finite element type method, we will need to establish a variational formulation for the time-harmonic problem (4) and (5) . We first introduce the following Sobolev spaces Unless otherwise specified, standard notations of Sobolev spaces are adopted. Then the weak variational formulation of PUFEM can be described as: where v is the complex conjugate of a test function v in V and a, b are bilinear and linear forms defined respectively as: The product of two tensor fields A = (A 1 , A 2 , A 3 ) and B = (B 1 , B 2 , B 3 ) in R 3 ×3 is given by The curl operator ∇× is defined for a vector field v = (v 1 , v 2 , v 3 ) by Under the previous notations and definitions, the tensor σ can be written as follows where the matrix [ ∇ × u ] × is defined by with u = (u 1 , u 2 , u 3 ) . Let us notice that the product σ(u ) : ∇ v in (11) can be written, using (14) , under a form involving the operators ∇, ∇ · and ∇ × : Then thanks to the boundary condition (5) , we rewrite a and b of expression (11) in the following formulations: From the mathematical point of view, according to the boundary condition (5) on , the problem in (10) and (17) is well-posed and an existence and uniqueness result can be established by virtue of Fredholm's alternative theorem and continuation arguments [42,43] , under weak assumptions, i.e., let us consider that is a bounded Lipshitz domain,

PUFEM approximated solution
We now need to approximate the problem stated in expression (10) by the PUFEM, based on plane wave basis and polynomial shape functions. For this purpose, let us consider a finite element mesh containing n nodes, denoted z, z = 1 , n . We denote by { N z } the partition of unity by polynomial finite element shape functions, and respectively by m P and m S the numbers of approximating P and S plane waves. Following the two-dimensional PUFEM approximation [5] , the displacement field u is approximated by superimposing pressure and shear displacements as follows However, in three dimensions the shear wave induced displacement is resolved into two components as follows Therefore, using the plane wave enrichment approach, the approximated displacement field u h is then expressed in the following form The orthogonal unit S-vector d S , ⊥ , to the unit P-vector d , is resolved into two orthogonal vectors denoted by d l, 1 S, ⊥ and d l, 2 S, ⊥ and are defined by such that and then with · 2 denoting the L 2 -norm. Note that, in 3D, there is an infinite number of vectors that are orthogonal to a given vector. As a result, finding the perpendicular vector d × P to d l P can be chosen in different ways, one of them being where d l P = ( sin θ l cos φ l , sin θ l sin φ l , cos θ l ) , with θ l and φ l being the two angles of the spherical coordinate system with 0 ≤ θ l ≤ π and 0 ≤ φ l ≤ 2 π . In the case of d l P with the same components such as d l P = (a, a, a ) then the perpendicular vector is given by d × In expression (20) , h describes the computational mesh size. For simplicity of the notations, the dependency of the numerical solution u h on the numbers of approximating P and S plane waves, m P and m S , is not indicated.
The above approximation (20) can be derived from the Helmholtz decomposition theorem, as discussed in [30] . The unknowns are no longer the nodal values of the displacement u h , but are now the amplitudes A P z,l , A S, 1 z,l and A S, 2 z,l attached to a node z and corresponding to P and S plane waves travelling in the directions d l P and d l S , respectively.

Choice of plane wave directions in space
In two dimensions, the choice of the plane waves and their uniform spacing or clustering around directions of preference is a trivial problem. The move to the three dimensional case presents however some difficulties. In fact, in three dimensions, it is neither intuitive nor an easy task to define m uniform divisions of the 4 π solid angle, m being any integer number, a part from a few known solutions. In previous work [15] , the authors used a uniform boundary meshing of a cube. A reasonably well spaced set of points is defined by the vectors joining the centre of the cube to each node on the cube's boundary. However, this approach is also limited to a few special cases of m for which a boundary-meshed cube is available.
The uniform distribution of a set of points on a sphere attracted the attention of researchers from various fields. Early work in meteorological modelling, for example by Kurihara [44] , resulted in algorithms that produce distributions based on the mapping of a regular grid on an equilateral triangle to each octant of the sphere. However, this is limited to values of m = 4 p 2 + 4 p − 2 , where p + 1 is the number of grid points from the north pole to the equator. The icosahedral grid of Williamson [45] uses the vertices of the regular icosahedron with the possibility of subdividing each triangular face in a similar manner to the Kurihara grid [46] or in a symmetrical manner such as performed in reference [47] . All these possibilities produce distributions with limitations on the values of m . For example, in reference [47] , it is given as m = 10 p 2 + 2 where p = 2 q with q ≥ 0 being the number of equal intervals into which each side of the original icosahedral triangles is divided.
With the motivation to define sets of plane waves forming bases for the DEM, Grosu and Harari [48] corrected the technique of Beverly [49] finding equally spaced points on latitudinal lines of a sphere which are themselves equally spaced along the longitudinal lines. This technique provides sets of approximately uniform distributions of points on a sphere with more possibilities, in comparison to the previously cited approaches (see Table 1 ).
In the case of the UWVF [40] , a library was used in which the selection of m directions in space, with 4 ≤ m ≤ 130, was based on minimising the distances between points on a surface of a sphere. The recursive zonal equal-area partition proposed by Leopardi [50] offers the possibility to generate any number m of points on the sphere through the partitioning into regions of equal area and small equivalent diameter. In this approach, it is required to find the colatitudes of polar caps (zones between adjacent latitudinal lines), then determine the collar angle which provides the number of regions in each collar.
More recently, Peake et al. [51] presented a method of uniformly spacing a set of points on a sphere based on static equilibrium of charged particles and hence named it the Coulomb force method. This work is motivated by the need to use uniform wave directions in space within PUBEM. It is shown that the performance of the enriched elements, within PUBEM, is not sensitive to small variations in the uniformity of the plane wave direction distribution. The interested readers in methods dealing with uniform distribution of points on a sphere are referred to the references included in [51] .
Using sets of plane waves in space with approximately uniform distributions is of great importance to the proposed PUFEM model. To this end, two approaches for selecting plane waves in space are used in this work. The first approach consists to use the Cube-boundary Meshing (CM) which was previously used in PUFEM for 3D Helmholtz problems [15] .
where φ 0 is a constant angle between two consecutive latitudinal lines. Fig. 1 shows distribution examples of 58 points obtained by the cube-boundary meshing (left), the corrected Beverly method (middle) and the Leopardi method (right). It is obvious that the cube-boundary meshing distribution provides points relatively more spaced around the equator in comparison to the regions around the south and north poles, while the two other methods provide relatively uniform distributions. Moreover, the corrected Beverly method and the Leopardi method lead to very similar distributions for m = 6 , 12, 22, 34 ... etc. If the south cap, or the north cap, of the Leopardi method is rotated by a certain angle, both methods lead to the same distribution.

PUFEM discretization
In this work, the displacement field is approximated via propagating plane waves in space in which both balanced ( m P = m S ) and unbalanced ( m P = m S ) choices of plane waves are considered in expression (20) of the approximated displacement field. Let us set A finite element approximation of the variational problem (17) is derived by using N z p l , N z s 1 l and N z s 2 l : The main feature of the discrete scheme (27), (28) and (29) , thanks to the stress tensor form (14) and the numerical approximation of the displacement field (20) , is that the resulting matrices can be easily computed as in the case of the scalar Helmholtz equation [12] , in terms of the oscillatory shape functions The other matrices involved in (27), (28) and (29) can be written, in a similar way as for the mass matrix, in terms of the product of the new oscillatory shape functions. The numerical evaluation of the above integrals is performed with high-order Gauss-Legendre scheme.
A nodal point z , with z = 1 , n, has m P + 2 m S degrees of freedom. So the total number of degrees of freedom is n (m P + 2 m S ) , and the unknown amplitudes A P z,l , A S, 1 z,l and A S, 2 z,l for z = 1 , n can be stored in n blocks of vectors in R m P +2 m S of the ) . The global unknown amplitude vector can be stored as follows and the global matrix resulting from (27), (28) and (29) has the following form   Galerkin weighting is used in this work and hence the global matrix of the resulting system is symmetric and block banded. It is stored using a steering vector to locate the elements and the solution of the final system is obtained via a direct solver [52] based on LDL T decomposition where L T is the transpose of the lower triangular matrix L and D is a diagonal matrix.

Numerical analysis
In this section, the developed numerical model for selected elastic wave problems is investigated with respect to the computational mesh size h and the plane wave enrichment, in terms of the numbers m P and m S of the approximating P and S plane waves. In the convergence analysis considered below, special care is taken so that the incident wave direction d inc P is not too close to the plane wave directions used in the approximating model. The computational domain is defined by where u represents the analytical solution of the considered problem.
To quantify the number of degrees of freedom (DOF) per P and S wavelengths, respectively, the parameters τ P and τ S are introduced in the following forms [5] τ P = π k P ( n m P ) 1 / 3 and τ S = π where n , as mentioned before, is the number of nodal points in the mesh grid.
Let us denote by α the average rate of convergence with respect to h such that the L 2 error ε 2 = O (h α ) , where α is to be evaluated numerically via α = ( log ε (1 / 4) 2 − log ε (1 / 2) 2 ) / ( log h 1 / 4 − log h 1 / 2 ) , when considering mesh grids h 1/4 and h 1/2 , for example. Table 2 shows the number of shear wavelengths λ S per element size h of the considered mesh grids ( Fig. 2 ), for circular frequencies ω = 1, 10, 15 and 20. The comparison of the mesh size with the shear wavelength reveals that for low circular frequency ω, each finite element contains a fraction of λ S whereas at higher ω, each finite element becomes multi-wavelength sized. For example, for the mesh grid h 1 2 , h = 0 . 16 λ S at ω = 1 . However, at ω = 20 , h = 3 . 18 λ S . This latest remark will have a direct consequence on the number n gauss of integration points required for the numerical evaluation of the element matrices. Therefore, this number is adjusted in the three spatial directions, within each element, by adopting Table 3 Progressive plane wave problem: L 2 error in % for ω = 10 .
( m p , m s ) ( the empirical expression giving n gauss = 10 × h/λ S + 2 , to ensure enough integration points are used. It is possible to develop semi-analytical integration schemes such as done in [6,53] to avoid the burden of the high order scheme, but this is not the objective of the current work. All computations are carried out in Fortran with double-precision complex numbers on Workstation Intel Core i7 3.5 GHz with 32 GB RAM.

Progressive plane wave test problem
Let us first consider a simple mathematical model consisting of a progressive plane wave in a three dimensional infinite elastic medium. This problem has an analytical solution and the expression of the displacement field is given by where d inc P is chosen such that (d inc ) = ( sin β 1 cos β 2 , sin β 1 sin β 2 cos β 1 ) , β 1 and the two components d 1 S, ⊥ and d 2 S, ⊥ of vector d S orthogonal to the direction d inc P are given by where (d inc P ) × is defined in expression (24) . The source term g of expression (5) is calculated from the traction σ( u ) n on , evaluated with the analytical solution u given in expression (37) . Note that this is only possible when exact boundary conditions are prescribed on .  (38) .
From the numerical results it is obvious that h -refinement improves the accuracy of the scheme for all cases of frequencies and plane wave enrichments. Increasing the numbers of approximating plane waves, for a fixed h , does also improve the accuracy of the plane wave enrichment model. In fact, this is expected as both the mesh refinement and the increase in the number of approximating plane waves lead to higher numbers τ P and τ S of DOF per P and S wavelengths, respectively.
Nevertheless, two remarks are worth mentioning. The first one concerns the cases of coarse mesh grids, in comparison to the wavelength, with low numbers of approximating plane waves. In such cases, the L 2 error is high and the reason lies in the fact that the discretization level represented by the numbers of DOF per P-and S-wavelength, τ P and τ S , is very low. For example, for the mesh grid h 1 2 at ω = 15 , the L 2 error is very high, ε 2 = 44 . 4% , for (m P , m S ) = (10 , 10) for which (τ P , τ S ) = (3 . 1 , 1 . 8) . For a higher enrichment, (m P , m S ) = (26 , 26) the L 2 error is a lot lower in comparison to the previous case, ε 2 = 1 . 13% , due to the increase of the discretization level to (τ P , τ S ) = (4 . 3 , 2 . 5) . When the latter is further increased to (τ P , τ S ) = (5 . 6 , 3 . 2) , for (m P , m S ) = (58 , 58) , the L 2 error is even lower, ε 2 = 0 . 097% , and then for (m P , m S ) = (98 , 98) the L 2 error reduces to ε 2 = 0 . 019% .
The second remark is related to the average rate of convergence α. From the three Tables 2 -5 , for fixed numbers ( m P , m S ) of approximating plane waves, the results show an exponential decrease of the L 2 error with respect to h -refinement. In general, the average rate α increases with the increase of the numbers of approximating plane waves, such as shown for the case of ω = 20 . However, continuing to increase the number of approximating plane waves, for a fixed circular frequency, does not necessarily lead to an increase in the rate of convergence. This is the case for ( m P , m S ) = (98, 98), at ω = 10 and 15 where the rate has decreased from 7.8 to 5.4 and from 5.9 to 5.1, respectively. This is most likely due to the ill conditioning issue, which is an inherent feature of the plane wave basis finite element approaches. This aspect is investigated in the next test example dealing with elastic wave scattering by a rigid sphere.

Wave scattering test problem
The problem of harmonic waves scattered by a spherical rigid body contained in a three dimensional homogenous and isotropic infinite elastic space is considered. An upward P-wave with amplitude 0 , for which the the potential is given by where 0 = i k p , impinges on a spherical body of radius a . The incident displacement field is obtained from The total displacement field u is a superposition of the incident u in and scattered u sc fields, u = u in + u sc . In fact, a P-wave impinging on a rigid body would lead to both scattered P and S waves, in an infinite elastic medium. The total displacement can be written in terms of its components with respect to the spherical unit vectors e r , e θ and e φ with where P ν are Legendre polynomials for ν = 0 , ∞ . The constants A ν and B ν are chosen such that u | r= a = 0 representing the total reflection condition on the boundary of the spherical rigid body. This is obtained via for all ν in N . The associated stress field σ( u ) is then obtained from the displacement field (42) in terms of its spherical and The symbols ε p(0) i , ε p(1) i and ε s (1) i , for i = 1 , 9 , are coefficients obtained form the various derivatives. The radial spherical functions are defined as, see [54] : where i = 0 or i = 1 such that ξ p(0) ν (k P r) = J ν (k P r) is the spherical Bessel function of first kind and order n and P 0 n ( cos θ ) is Legendre's polynomial of order n and degree m = 0 . A = (−i ) ν (2 n + 1) and ξ p (1) The stress field is then used in the evaluation of the source term g of expression (5) to be used in the boundary condition of the elastic wave scattering problem. The series used in expressions (43) and (44) are truncated when a sufficient number N t of terms is included. This number is taken such that N t ≥ k S R where R is the radius of the most distant point of the considered computational domain [55] . For illustration purpose, Fig. 3 shows contour plots of the components of | Re ( u h )| in the computational domain for the case of ω = 40 when a P-wave travelling upward along the x 3 -axis impinges on a spherical rigid object of unit radius. The centre of the rigid scatterer coincides with the centre of the Cartesian system and hence the considered computational domain is in the vicinity of the scattering object. The scattered field contains both body waves, P and S, which interfere around the hard scatterer. The first and second components are opposite each other and the third component has higher values due to the incident wave propagating along x 3 -axis. To obtain the contour plots of Fig. 3 , the computed nodal values from the problem solution are used in conjunction with the model of expression (20) .
The same numerical analysis is followed for this elastic wave scattering problem considering the same parameters already used for the case of the progressive plane wave test problem but with more options regarding the used sets of approximating plane waves. Indeed, in this case, plane wave distributions obtained form the Cube-boundary Meshing (CM) and the Corrected Beverly (CB) method are used. The two approaches, CM and CB, offer more flexibility in the choice of uniformly distributed plane waves in space and even when they allow the same number of plane waves the CM and CB distributions are different. Tables 6 -8   approximating plane waves are increased, for ω = 10 , 15 and 20, respectively. They also give the average rate of convergence computed numerically from the obtained L 2 errors.
From Tables 6 -8 , it can be seen that, in general, the obtained results lead to similar conclusions already drawn from the case of the progressive plane wave test problem. Basically, the L 2 error improves as the numbers ( m P , m S ) of the approximating P and S plane waves increase for all cases of mesh grids and frequencies presented here. However, we notice in certain cases after reaching a low L 2 error a further increase in the number of enriching plane waves does not necessarily improve the accuracy such as reflected in the case of h 1 4 and ω = 10 when increasing ( m P , m S ) form (58, 58) to (106, 106). In all these cases, the L 2 error remained unchanged at ε 2 = 0 . 0 0 0 02% . For frequencies ω = 15 and ω = 20 , this is again seen for the mesh grid h 1 4 when ( m P , m S ) increased from (98, 98) to (106, 106) and for both enrichments the L 2 error stagnated at ε 2 = 0 . 0 0 0 04% and 0.0 0 0 06%, respectively. Moreover, in the case of frequency ω = 20 and the mesh grids h 1 and h 1 2 , for high numbers of enriching plane waves the quality of the results deteriorates instead of improving, or at least stagnating. This is seen when ( m P , m S ) increased from (98, 98) to (106, 106), for which the L 2 error increased from 4.00% to 7.19% in the case of mesh grid h 1 and increased slightly from 0.0022% to 0.0024% for the mesh grid h 1 2 . This is known to be caused by the ill conditioning aspect, which is investigated here.
The results overall also confirm the exponential convergence with respect to h -refinement. It is noticed that the average rate of convergence increases with the increase of the numbers of approximating plane waves up to a certain level but continuing to increase ( m P , m S ), for a fixed circular frequency, does not necessarily lead to an increase of the rate of convergence. In fact, it could even decrease such as for the cases of ( m P , m S ) > (58, 58), for ω = 10 , in which α decreased from     6.1 to 4.7, 4.0 and 3.5. It can also be seen that the L 2 error stagnated at the value of 0.0 0 0 06% for the mesh size h 1 4 . As will be shown next, this is obviously due to the ill conditioning of the plane wave enrichment model which is known to occur for relatively small size elements, in comparison to the wavelength, with high numbers of approximating plane waves.
The behaviour of the condition number denoted κ = W W −1 , with respect to the 1-norm, is numerically investigated for the considered circular frequencies ω = 10, 15 and 20, mesh grids h 1 , h 1 2 and h 1 4 , and various plane wave enrichments.
The package MF71 of the Harwell Subroutine Library [56] , which provides an estimate of the 1-norm for sparse or not explicitly available matrices, is used with our LDL T linear solver to evaluate the condition number. Fig. 4 clearly shows that, for fixed circular frequency ω and plane wave enrichment m = m P = m S , the conditioning deteriorates with h -refinement. In general, log κ increases linearly with m as shown for the mesh grids h 1 and h 1 2 but when the condition number is relatively high, which is the case of mesh grid h 1 4 , the machine precision does not allow the computation of the true condition number anymore ( Fig. 4 ). However, it is obvious that it decreases with the increase if the circular frequency for which the elements are relatively of large size in comparisons to the wavelengths. For example for the mesh grid h 1 4 , κ ≈ 10 22 for ω = 10 , while it is about 10 18 and 10 16 for ω= 15 and 20, respectively. This shows that the proposed numerical model is more suitable for the high frequency range where coarse mesh grids, in comparison to the wavelength, are used with relatively high numbers of approximating plane waves. For such cases, it is obvious from the results of Fig. 4 that the condition number remains within acceptable limits. Previous work [6] has shown that a careful choice of the numbers of approximating P and S plane waves may improve the accuracy of the PUFEM and its conditioning. This is due to the fact that k S is always greater than k P and hence the multi-wavelength mesh size h contains more S wavelengths than P wavelengths. Therefore, it makes sense to choose more S plane waves than P plane waves in the approximating sets, m P < m S . It was suggested to choose m S / m P of the same order of k S / k P [6] . See also references [31] and [57] .
In this case, the scattering problem with ω = 20 is solved with various combinations of m P and m S plane waves. The The results of Fig. 6 show the behaviour of the L 2 error and the condition number for various combinations of m P and m S plane waves, with the total number of approximating plane waves being m P + 2 m S . It is obvious again, as previously mentioned, that increasing the number of approximating plane waves in general improves the quality of the results but leads to an increase of the condition number. Moreover, it is also revealed that balanced choices of P and S plane waves do not lead to the best quality of results. This is clearly shown by the cases of m S = 102 and 128. For the first case, m S = 102 , the lowest L 2 error is provided for the total number of plane waves equal to 288 leading to m P = 84 . For the second case, m S = 128 , the lowest L 2 error is obtained for m P + 2 m S = 322 giving m P = 66 . In both cases, m S / m P > 1 and more precisely m S /m P = 1 . 21 and 1.93, respectively. It is worth noting that the suggestion to choose m S / m P of the order of k S / k P ≈ 1.73 was drawn from past work for indication [6] . But given that the element sizes are not constant throughout the mesh grid it is anticipated to see variations in the values of the ratio m S / m P . As for conditioning, it was expected to increase with the increase of the total number of approximating plane waves, as the mesh size and the frequency are both fixed.

FEM and PUFEM solutions for elastic short wave scattering problems
The plane wave basis finite element model presented in this work is compared to polynomial based finite elements for the solution of the scattering problem dealt with above for the range of frequencies ω = 1, 5, 10 and 15. This is achieved by comparing the accuracy of the results based on the L 2 error, the total number of degrees of freedom (totdof) required in the solution as well as the total number of storage locations (totsys) of the final system matrix to be solved. The resolution for the tri-quadratic case. The refinement approach ensures that each S-wavelength is modelled by a sufficient number of elements reflected in the discretization level τ S .
Before considering the results, it is worth mentioning that this is not a fair comparison since PUFEM is a high order approach while for FEM only low order elements are used. However, from the practice point of view, when solving engineering problems it is very common to use low order elements, mainly linear and quadratic, and if these low order elements are enriched via the incorporation of a priori knowledge of the wave field in the numerical model, such as done here with the use of enriching P and S plane waves, improved quality results with low discretization levels would be achieved, as will be shown next.
The results of Table 9 show that the PUFEM approach leads to better quality results as well as to significant reductions of the total number of degrees of freedom (totdof), the total number of storage locations (totsys) and the number of degrees of freedom per wavelength ( τ S ). For example, for ω = 1, mesh grid h 1 combined with a plane wave enrichment with m = 10 and a discretization level τ s = 20, PUFEM leads to a very low L 2 error of 0.003%. For the same frequency, FEM1 produced an error of 0.33% with τ s = 40 . 77 while requiring more than double in terms of totdof and totsys. At other circular frequencies, the performance of the PUFEM model is even more obvious as it leads to relatively very low L 2 errors with discretization levels as low as 2.9, for example in the case of ω = 15 and m = 98 . For the latter case, the L 2 error ε 2 = 0 . 189% while with FEM1 the most refined mesh leads to ε 2 = 1 . 66% with a discretization level of 14 DOF/ λ S and very high numbers related to totdof and totsys. Using PUFEM, offers the possibility to further enrich the model with m = 154 and hence decrease the L 2 error to 0.0199% or even use the mesh grid h 1 2 with m = 58 to achieve ε 2 = 0 . 0251% . In both cases, the discretization level with respect to the S-wavelength remains around 3.3. Furthermore, with PUFEM, there is also room to increase the circular frequency of the problem and achieve good quality results with practical parameters for totdof and totsys. However, in the case of FEM, further refining the mesh would lead to more than a quarter of a million entries for totdof and around 2 billions for totsys, which becomes impractical to run on most personal computers. In such cases, computers with high specifications would be required. Table 9 also shows that using tri-quadratic finite elements (FEM2) improves the quality of the results in comparison to tri-linear elements (FEM1). The tendency shows that it also reduces the required total number of degrees of freedom (totdof) and the total number of storage locations (totsys). Continuing to increase the order of the polynomial based FEM may lead to significantly better quality results and lower computational requirements in terms of totdof, totsys and τ S . As already mentioned, comparing PUFEM to FEM1 and FEM2 is not fair as significantly better results could be obtained using elements of higher order than quadratic. This should be investigated but, in practice, often linear and quadratic elements are used. Moreover, the incorporation of the wave numbers k P and k S into the enrichment model of expression (20) , via P and S plane waves, incorporates valuable knowledge in the presented finite element model about the wave problem at hand. This latter aspect may play a key role in making PUFEM competitive in solving elastic wave problems at high frequencies.
Last, it is important to add that for PUFEM a high integration scheme is used to evaluate the highly oscillatory integrands of the element matrices. This uses thousands of integration points at high frequencies. So, at this stage, the computational time is mostly taken by the assembling process while the solution part represents only a fraction of the total time thanks to the drastic reduction of the required total number of degrees of freedom. However, developing semi-analytical integration schemes in three dimensions similar to those presented in the two-dimensional case [6,53] would significantly reduce the computational effort and would make the PUFEM model even more attractive.

Conclusions
In this work, plane wave enriched finite elements are developed within the framework of PUFEM for the solution of three dimensional elastic wave problems. These elements are capable of containing many wavelengths per nodal spacing and consequently allow the relaxation of the traditional requirement of several nodal points per wavelength, used in low order polynomial based FEM. This is achieved by expanding the displacement field into discrete series of displacements with respect to many directions corresponding to P and S plane waves, each propagating at a specified angle in the three-dimensional space. The displacement due to a P plane wave lies in its direction of propagation whereas the displacement due to the S plane wave is normal to its propagation direction and hence it is resolved into two components contained in the plane normal to the S plane wave direction.
It is shown that the proposed approach provides better quality results with significantly reduced requirements in terms of the total number of degrees of freedom and total number of storage locations due to the fact that the elements may contain many wavelengths per nodal spacing. However, in order to make the current model competitive, the high order numerical integration issue must be addressed by developing fast integration schemes similar to those produced for the two-dimensional elastic wave problems [6] .
It is shown that using an unbalanced choice of the plane wave enrichment leads to better quality results and lower condition number in comparison to the balanced choice option. Moreover, flexibility in the choice of the plane wave directions in space is an important aspect for the good performance of the method. This could be overcome by using some of the many strategies already developed in other research fields, some of them referenced herein, and providing uniform distributions for any chosen number of directions.
Last, ill conditioning being an inherent feature of the plane wave enrichment, even if it is possible to choose parameters in terms of the mesh size and plane wave enrichment for a given frequency such that the condition number is kept within acceptable limits, developing suitable preconditioners as well as using iterative solvers for large problems would be of practical interest.