Parametric finite elements, exact sequences and perfectly matched layers

The paper establishes a relation between exact sequences, parametric finite elements, and perfectly matched layer (PML) techniques. We illuminate the analogy between the Piola-like maps used to define parametric H1-, H(curl)-, H(div)-, and L2-conforming elements, and the corresponding PML complex coordinates stretching for the same energy spaces. We deliver a method for obtaining PML-stretched bilinear forms (constituting the new weak formulation for the original problem with PML absorbing boundary layers) directly from their classical counterparts.

responding PML via complex stretching. If nothing else, this helps in coding the PML for multiphysics problems involving a simultaneous discretization using H 1
The perfectly matched layer (PML) was introduced by Bérenger [1] as a technique for reflectionless absorbing of electromagnetic (EM) waves at the boundary of the domain of interest. Absorption of the waves inside of the PML was achieved by matching material conductivities, and by splitting the EM field components into subcomponents. In [2], Chew and Weedon proposed an alternative formulation of the PML through a complex coordinate stretching of the spatial variables of the original differential equations posed in the frequency domain. The approach was later reinterpreted in terms of an analytic continuation [3,4], and motivated the extension of the PML to curvilinear coordinates, conformal mesh terminations, and more general media. Simultaneously, the new technique was successfully applied to acoustics [5], elastic wave propagation [6], and poroelasticity [7].
Sacks et al. [8] developed a non-split version of the PML for Maxwell's equations, preserving the original form of the equations (the so called Maxwellian PML) at the expense of introducing anisotropic material properties. In consequence, the resulting field in the PML can be interpreted as a physical field in the aforementioned anisotropic medium. This is the manifestation of the metric invariance of Maxwell's equations.
For more general hyperbolic equations (e.g. elastodynamics), due to the presence or absence of specific symmetries in the physical fields in the PML introduced medium, the interpretation of the PML in terms of an anisotropic material seems to be impossible. A rigorous proof of stability for such systems remains open [9].
In the series of papers [10][11][12], Teixeira and Chew showed that an appropriate change of field variables (analogous to Piola transforms for H(curl)-and H(div)-conforming fields) combined with complex coordinate stretching leads to Maxwellian PML for general curvilinear coordinates and anisotropic linear media. A similar approach was used for elastic wave propagation in [13,14], where the H(div)-conforming Piola transform was applied to the stress tensor.
Finally, in [15], Teixeira and Chew proposed a generalized geometric viewpoint to the PML obtained via complex coordinate stretching in the Fourier domain. It appears that the analytic continuation is equivalent to a complexification of the space metric tensor. Since the spatial differential operators depend on the metric, the PML amounts to a modification of these operators due to the new complex metric. This enables the construction of more elaborate versions of PML, provided one can determine the transformation matrix relating original and complexified metric tensors, which in turn is used to define stretched differential operators or field variables.
In this paper, we present a unified framework for a systematic derivation of complex stretchings for all fields and exterior differentiation operators present in the classical grad-curl-div exact sequence. We expose the algebraic structure of the construction for arbitrary curvilinear systems of coordinates, and establish practical rules for transformation of bilinear forms corresponding to classical weak formulations. As a result of the construction, the PML-stretched operators and spaces preserve the exact sequence property.
The content of the paper is as follows. We begin in Sect. 2 with a review of the notion of the exact sequence and its connection with a proper definition of FE spaces for parametric elements through the use of Piola-like transforms. Section 3 describes the PML technique, and reinterprets +it in terms of the new Piola-like complex transforms. In Sect. 4, we illustrate the general technique with two specific examples of calculations of Jacobian matrices describing the PML transformations. Finally, in Sect. 5, we show how to derive the bilinear forms corresponding to various PMLs by a direct "stretching" of the classical bilinear forms for several standard variational formulations involving the exact sequence energy spaces. We conclude with a summary in Sect. 6.

Exact sequence and parametric elements
The following energy spaces, used henceforth, are implied by physics: Along with operators of grad, curl, and div, for a simply connected domain Ω, the spaces form the well-known exact sequence shown below. The sequence can be reproduced at the level of a single element or a Finite Element (FE) mesh using various versions of Lagrange, Nédélec, Raviart-Thomas, and discontinuous elements, possibly with a variable order (the hp elements), see [16,17]. Figure 1 presents the lowest order elements on a tetrahedral element. The corresponding degrees-of-freedom for H(curl) elements are defined using tangential components along edges, and for H(div) spaces, using normal components over faces; hence, the frequently-used names of edge and face elements.
The de Rham complex (1) combines the continuous and discrete versions of the sequence with specially constructed Projection-Based Interpolation (PBI) operators that make the diagram commute [18].
Notation. We use superscripts for contravariant and subscripts for covariant components of vectors and tensors. The Latin alphabet is used for coordinates and indices in Cartesian coordinate systems, and the Greek alphabet for their curvilinear counterparts.

Cartesian coordinates
We reproduce the reasoning from [17, p. 34]. Given a bijective smooth map T y :Ω y → x ∈ Ω (alternatively x i = x i (y j ), i, j = 1, 2, 3), transforming a master element Ω ⊂ R 3 onto a physical element Ω ⊂ R 3 , and master element H 1 -conforming shape functionsû(y), we define the H 1 -conforming shape functions on the physical element as compositions of the map T y and functions defined on the master element, In order to preserve the exact sequence property, we have to define the H(curl)-, H(div)-, and L 2 -conforming shape functions (and thus the corresponding spaces) consistently with the transformation rules for the differential operators. The transformation rule for gradients follows from the chain differentiation formula, and, therefore, we define the H(curl)-conforming shape functions on the parametric element as follows or, more precisely, In turn, for the curl operator we have where i jk is the permutation Levi-Civita symbol. But, from the definition of the determinant of a 3 × 3 matrix, we have where J y is the Jacobian of the map Consequently, This leads to the definition of the H(div)-conforming shape functions, Finally, using the differentiated identity (2), we have which establishes the transformation rule for the L 2 -conforming shape functions, We summarize the results in the table below, This can be rewritten in the form of geometric relations Defining finite element spaces W hp , Y hp , V hp , and Q hp , and using the transformation rules listed above, we preserve the exact sequence property.
Transformations of tangent and normal vectors. Assume we are given a parameterized curve y i (s). We define a tangent vector, depending on a parameter s, The vectort is transformed into a tangent vector t via a mapping T y as follows thus we have a geometric transformation for tangent vectors A vectorn normal to the surface, defined by a pair of noncollinear tangent vectorst 1 andt 2 , can be expressed bŷ The transformed normal vector can be calculated as follows where we use the following identity, valid for any 3×3 matrix M 2.2 Curvilinear to Cartesian coordinates Assume now that we have a curvilinear system of coordinates Denoting by e i the unit vectors for Cartesian system x i , we start with the formulas for the basis and cobasis vectors, This implies relations between curvilinear and Cartesian components of vectors, We record now the formulas for differential operators Here ε αβγ = αβγ /J ξ denotes contravariant components of the third-order alternating tensor, and J ξ = det(∂ x i /∂ξ α ).
Notice that the formulas favor using covariant components for gradients and H(curl)-conforming fields, but contravariant components for the curls and H(div)-conforming fields. For Cartesian coordinates, the new formulas reduce to the original ones.

Curvilinear to curvilinear coordinates
Assume now that we have a curvilinear system of coordinates y j = y j (η β ) with covariant and contravariant bases b α and b β , and the map x i (y j ) is given through curvilinear coordinates ξ α and η β , namely Figure 2 presents more clearly the above transformations. Functions u, E, V, q are given in terms of x i or ξ α , and u,Ê,V,q are functions of y j or η β , respectively. Using the formulas for grad, curl, and div in curvilinear coordinates, we can rewrite now the Piola transforms in terms of the curvilinear coordinates ξ α and η β . Transformation for the components of the gradient implies the transformation rule for vectors in the space H(curl) Thus, for the curl, using previously derived relationships, we have where we defined the new alternating 3rd order tensorε δβγ = δβγ /J η , and denoted J η = det(∂ y j /∂η β ). The notation for the Jacobians follows the notation for the corresponding transformations. This leads to the formulas for curls and vectors in the space H(div) Finally, for the div operator, we have We summarize the formulas for the Piola transforms in the curvilinear coordinates in the table below where J y = J ξ det(∂ξ α /∂η β ) J −1 η . Again, note that in formulas (8), gradients and H(curl)-conforming fields are represented in terms of their covariant components but curls and H(div)-conforming fields are represented in terms of their contravariant components.
For Cartesian coordinates, transformations (8) reduce to the original transformations (3). Both sets of formulas are coordinate representations of geometric relations (4) with the appropriately defined Jacobian of the transformation. In the most general case of curvilinear coordinates, the Jacobian J is represented in terms of mixed components, namely with its determinant given by,

Perfectly matched layer
In this Section, we review in detail the concept of complex stretching behind the perfectly matched layer technique applied to boundary-value problems involving grad, curl, and div operators. We reinterpret the method in terms of new complex Piola-like transforms, showing the analogy between the PML transforms and the transforms for conforming parametric elements.

PML in Cartesian coordinates
We shall start again with Cartesian coordinates. Consider a boundary-value problem formulated in an unbounded domain. Construction of the PML, based on the analytic extension argument, is done as follows.
Step 1: Assume or prove analyticity of the solution in one or more arguments. Typical examples include exterior wave propagation problems. Dependent upon the material properties of the medium extending to infinity, we may assume analyticity in one or more variables. For instance, for a 2D problem involving horizontal layers, we expect analyticity in x (direction of the layers) but not in y (direction perpendicular to the layers). Typically, the solution is not even C 1 across the layers.
Step 2: Extend involved operators to the complex domain.
Once we have assumed (proven) that the solution is analytic in a specific variable, we can complexify the variable. The real domain for the variable being complexified is replaced with a complex domain. The original function u = u(x i ) is replaced with its unique analytical extensionũ =ũ(x i ). Without any loss of generality, we will assume that we have complexified all variables x i . The corresponding derivatives with respect to the complexified variable are now understood in the complex sense. Differential operators like grad, curl, or div remain unchanged except that the complex derivatives replace the real derivatives.
Step 3: Define complex stretching, i.e. a map from a new real domain into the complex domain.
In the special case of stretching each coordinate separately, we havẽ The goal is now to replace the original solution In the domain of interest, say in a box |x j | < a, the mapx(x) = x, so the compositionũ•x coincides with u, the solution remains unchanged. Outside of the domain of interest, the stretching (11) is constructed in such a way that, for outgoing waves, the compositionũ •x diminishes exponentially with growing x j , but for incoming waves, it grows exponentially in x j . The situation resembles the behavior of the solution in the Laplace domain, although there we complexify the frequency, not a particular coordinate. Consequently, a radiation condition that is supposed to eliminate waves coming from infinity, after the stretching can be replaced with a simple Dirichlet condition at infinity. In practice, we truncate the domain at a finite length, e.g. at |x j | = b > a, making sure that the exponentially decreasing compositionũ•x reaches a machine zero at the truncating boundary. Finally, if the solution is analytic only in some variables, we do not stretch the variables in which we do not have an analytic extension.
Step 4: We identify differential equations satisfied by the stretched solution. This amounts simply using the chain rule for the compositionũ •x. For instance, Notice that ∂ũ/∂x i is the complex derivative. It is essential to work with the formulation based on a system of first-order equations (see the Remark 2 below).
Step 5: Stretch the metric (Jacobian). Once we have determined the "stretched" first-order equations, we proceed with the derivation of the appropriate weak formulation. We follow the choices made in the derivation of the original formulation, and relax or keep in the strong form the same equations, comp. [19,17]. When integrating over the stretched domain, we need to account for the Jacobian of the stretching, Without this step, we might lose the symmetry of the bilinear form.

Remark 1
Steps 4 and 5 can be replaced with a direct modification (stretching) of the original bilinear form, including the complexification of the Jacobian. This is perhaps the easiest way to derive the PML modified bilinear forms in practice. This ends the whole procedure for H 1 -conforming discretization. For problems involving H(curl)-, H(div), or L 2 -conforming discretizations, one extra step is necessary and illuminating this step is the main purpose of writing this contribution.
Step 6: Use complex Piola transforms to identify new unknowns in the PML domain. This step is necessary to preserve the H(curl), H(div), and L 2 energy spaces, see e.g. [17] for a detailed discussion in the context of Maxwell's equations.

Remark 2
It is natural to derive the PML-stretched equations using the system of first-order equations. This point was already made in [17] where several classical weak formulations have been derived without even using of secondorder PDEs. One practical advantage of the approach lies in the fact that one can clearly see which equations have been relaxed (integrated by parts) and which have been kept in the strong form. Another important factor is the simplicity. We never have to work with the second-order differential operator, a non-trivial burden in context of curvilinear coordinates.
The same comments apply to the PML-stretched equations. Additionally, one can clearly see the need for introducing the stretched Jacobian, necessary for preserving the symmetry of the bilinear form. The ultimate variational formulation for the stretched equations can be obtained directly by "stretching the original variational formulation," one of the points made in this paper.

PML in arbitrary coordinates
We shall consider now the most general possible case of complex coordinate stretching in all coordinates, where n = 1, 2, 3 is the dimension of the space. Obviously, this covers all possible cases of partial stretchings in less than n variables. Here x = (x 1 , . . . , x n ) stands for the ultimate Cartesian coordinates. The corresponding Jacobian matrix will be denoted bỹ withJ −1 denoting the inverse matrix of the complex-valued Jacobian matrix J. Suppose now that the original formula for the bilinear form or a first-order differential operator includes gradient ∇ x u where u is the unknown solution or a test function. Upon switching to complex extensionũ =ũ(x), the gradient is replaced with ∇xũ. Composing ∇xũ with complex stretching, we obtain back a function of real coordinates x, Denoting the modified function u :=ũ •T by the same letter as for the original function, the chain formula implies that or, using the indexless notation, which is equivalent to This is the formula we need to use when computing the bilinear form in the PML domain. The logic of the exact sequence implies that we have to use the same transformation for H(curl)-conforming fields, whereẼ is the complex extension of an original vector field E and, denoted by the same symbol, the field on the right-hand side is the new unknown in the stretched part of the domain.
Following the logic of exact sequence we arrive at formulas fully analogous to the Piola pull-backs used in the definition of parametric elements whereJ denotes the Jacobian, i.e. the determinant of the matrixJ.

Examples of Jacobians for coordinate stretchings
We present here examples of transformations and derivations of Jacobian matrices for two representative cases defined in curvilinear coordinates.

Example of a problem stretched in curvilinear but solved in Cartesian coordinates
Complex stretching is done frequently in curvilinear coordinates. The most common example perhaps involves spherical coordinates in 3D where the stretching takes place only for the radial coordinate r . We start with the comment that even for curvilinear stretchings it may be convenient to use Cartesian formulas (14). This is a consequence of the fact that FE computations take place, most of the time, in Cartesian coordinates. We will illustrate this point with a spherical stretching example. We use the standard spherical coordinates, and their analogous relations for complex variables ⎧ ⎨ ⎩ x = r sin ψ cos θ y = r sin ψ sin θ z = r cos ψ ⎧ ⎨ ⎩x =r sin ψ cos θ y =r sin ψ sin θ z =r cos ψ Notice that we have not put the tilde over ψ and θ as we postulate analyticity and stretch in the r variable only. After stretching, we get, ⎧ ⎨ ⎩x =r (r ) sin ψ cos θ y =r (r ) sin ψ sin θ z =r (r ) cos ψ It is easy now to compute the Jacobian matrix, where Letting the computer do the algebra for us is the best way to avoid algebraic mistakes.

Example of a problem stretched and solved in curvilinear coordinates
However, if a boundary problem is formulated from the very beginning in curvilinear coordinates (e.g. axisymmetric problems), we need the equivalent of formulas (14) in curvilinear coordinates. Geometric formulas (14) remain valid except that we need complex equivalents of formulas (9) and (10), Note that, along with scalarsξ,J ξ , one needs also to complexify vectors a α .
Cylindrical coordinates. We use the standard cylindrical coordinates, and their complex counterparts Notice that we have not put the tilde over θ as we postulate analyticity and stretch in the r and z variables only. Additionally, we define transformationT : ξ →ξ as follows ⎧ ⎨ ⎩r =r (r ) θ = θ z =z(z)

Transformation of weak formulations
We show now how to obtain the PML-stretched versions of weak formulations by a direct application of the proposed transformation rules on the classical bilinear forms corresponding to the problem of interest. The rules are as follows: 1. Start with the original (classical) weak formulation and identify the energy spaces for all the quantities present in the governing equations, 2. Calculate the stretched Jacobian for designed PML, using Eqs. (13) or (15), depending on the selected system of coordinates, 3. In the weak form, replace each field (e.g. u, E, V, q) with its stretched counterpart (ũ,Ẽ,Ṽ,q), 4. Replace each differential operator with its stretched counterpart, 5. Multiply the integrand by Jacobian (determinant of the Jacobian matrix)J in each volume integral, 6. Replace the normal vector n by its stretched counterpart n in each surface integral, 7. Apply transformations (14) to express the integrand in terms of the original field variables and classical differential operators.
Application of this method to the classical bilinear form results in the final PML-stretched weak formulation. This can be also achieved by performing all the classical steps of development of a weak formulation, starting from the PML-stretched differential equations. The former approach involves much less effort.

Examples of transformation of bilinear forms
Linear Acoustics. The governing equations for acoustics are obtained by linearization of the isentropic form of the compressible Euler equations around the hydrostatic equilibrium. The resulting system of equations in the frequency domain has the form iωp + ρ f c 2 f ∇ x · u = 0 (conservation of mass) iωρ f u + ∇ x p = 0 (conservation of momentum) (17) where p and u are fluid pressure and velocity (perturbations), respectively, c f is the speed of sound in the fluid, and ρ f is fluid density.
A formulation in which the continuity equation is treated in a weak sense, and the momentum equation point-wise leads to a classical variational formulation for the scalar Helmholtz equation posed in H 1 (Ω) space with the bilinear form whereÃ =JJ −1J−T is a symmetric matrix, and bothÃ and k 2 fJ can be treated as the new material parameters. Similarly, reverting treatment of the Eq. (17), we arrive at the alternative bilinear form for acoustics, posed in the PML stretching results iñ In both cases, we recovered the initial form of the governing equations, changing only material parameters.
Maxwell's Equations. The governing equations for EM wave propagation in the frequency domain have the form where E and H stand for electric and magnetic field, σ, μ, and ε are electric conductivity, permeability, and permittivity, respectively, and I i stands for prescribed impressed electric current.
The classical formulation for that problem gives the following bilinear form posed in H(curl, Ω) space Thus, once again, we recovered the initial form of the governing equations.
Linear Elastodynamics. The governing equations for linear elastodynamics in the frequency domain have the form where u i are components of the displacement, σ stands for the Cauchy stress tensor, ρ is the solid density, and C is the 4th-order elastic tensor. The classical bilinear form for that problem, posed in Understandingṽ i, j as the components of the second-order tensor ∇xṽ, PML stretching results iñ whereC i pkq =JJ −1 pjJ −1 ql C i jkl andJ ρ are the new material parameters. The form of the equation seems to be apparently preserved; however, we have lost the important properties of the elastic tensor. For example, if the Jacobian matrixJ is diagonal (when each coordinate is stretched separately, see Eq. (12)) we havẽ It is easy to see that the new tensorC preserves the major symmetry, i.e.C i jkl =C kli j , but both the minor symmetries are lost, namelyC i jkl =C jikl andC i jkl =C i jlk .

Examples of transformation of bilinear forms for coupled problems
When a coupled problem is solved, a need for the calculation of boundary integrals coupling two physical domains within the PML layer arises. In those cases, one has to transform the coupling equation into its stretched form, and proceed in the analogous way as for standard bilinear forms.
Coupled acoustics-elastic wave propagation. There are two coupling conditions that need to be considered (t is the traction vector defined on the boundary): n · u f = iωn · u s (normal velocity continuity), − pn = t (traction continuity). , where u f and u s dentote fluid velocity and solid displacement, respectively. The boundary integral for the pressure formulation of linear acoustics has the form − iωρ f qn · u f d .
The boundary integral for elastodynamics has the form t · v d .
Inserting the stretched form of the traction continuity condition (−pñ =t), we obtain However, when the coupling interface is perpendicular to the PML-stretching direction (Fig. 3), which is often the case, and we consider simple coordinate stretching (12) resulting in diagonal Jacobian matrix, the termJJ −T n simplifies Thus, the only change in the classical surface coupling terms is the presence of the stretched Jacobian determinant.

Summary
We have shown a close correspondence between the transformations used in defining conforming parametric elements originating in the exact sequence, and the complex transformations defining the reinterpreted PML stretchings for problems involving grad, curl, and div operators. We have formulated a procedure for developing PMLstretched weak formulations by a direct transformation of the original bilinear forms. This approach requires a precise identification of the involved energy spaces. The method has been illustrated with several examples. It is important to note that the results can also be obtained by: (a) replacing the differential operators in the governing differential equations with the stretched ones, (b) applying the proposed transformation rules both to the differential equations and the field variables, (c) multiplying with the test functions and integrating over the problem domain and, finally, (d) performing integration by parts. Obviously, such an approach is more complicated and prone to making mistakes.