We consider the efficient numerical solution of electromagnetic wave propagation modeled by time-dependent Maxwell’s equations in second order form
Here is the electric field, and
are the symmetric and positive definite permittivity and permeability tensors, andis the current density. Conduction currents can be included by setting , where are the impressed source currents and is the electric conductivity.
Today’s industry standard for solving Maxwell’s equations in time-domain are the finite difference time domain method and the finite integration technique [20, 21] which provide second order approximations on rectangular grids and for smooth and isotropic coefficients . Due to the underlying explicit time-stepping schemes, they lead to very efficient and accurate numerical approximations. Non-trivial modifications are, however, required to guarantee stability of the schemes in the case of non-rectangular grids and discontinuous or anisotropic coefficients , which in general also leads to reduced convergence orders.
A flexible alternative is provided by Galerkin approximations based on conforming finite elements, for which a rigorous stability and convergence analysis is possible under rather general assumptions. The standard finite element approximation with Nédélec elements of order
, leads to error estimates of the form
which are optimal in view of the approximation properties of these spaces; see [15, 16, 17] for details. A major drawback of standard finite element approximations for wave propagation problems, however, is that due to the required conformity of the basis functions, the linear systems
arising from discretization in space have a mass matrix which is sparse, but does not have a sparse inverse. This prohibits an efficient time integration by explicit time stepping schemes.
In order to overcome this source of inefficiency, mass-lumping strategies can be applied, which aim at replacing the mass matrix by a (block) diagonal approximation , in such a way that the overall accuracy of the approximation is not reduced. A systematic analysis of such schemes is possible, if mass-lumping can be interpreted as inexact numerical integration. In this spirit, mass-lumping for finite element methods for Maxwell’s equations on quadrilateral and hexahedral grids has been investigated in . In fact, a close relation exists between finite difference schemes [20, 21] and low order finite element approximations with mass-lumping; we refer to [4, 13] for details.
In order to obtain the full geometric flexibility of finite element approximations, we here consider mass-lumping for Maxwell’s equations on tetrahedral meshes, for which only few results are available. Lowest order Nédélec elements of type one and two has been proposed in  and first order convergence has been established. Related methods have been proposed in  in the context of the finite integration technique, but no convergence analysis is given there. A mass-lumping strategy based on an extension of the lowest order elements has been proposed by Elmkies and Joly  and first order convergence has been illustrated by a numerical dispersion analysis. Second order convergence has been observed by the authors for an appropriate extension of the second order Nédélec element , which we call element in the following; we refer to [4, 10] for details.
As a first result of this paper, we will prove that
the element with mass-lumping yields in fact second order convergence, if , but in general, only first order convergence can be obtained.
Since is satisfied when , our analysis also explains the good convergence behavior observed in the numerical tests in [4, 10]. Our proof of second order convergence when is based on a detailed analysis of quadrature errors, which also provides insight into the cause for the convergence order reduction in the general case. This allows us to
propose a modification of the element which, together with appropriate mass-lumping, leads to second order convergence in the general case.
In fact, only one of the basis functions of the element has to be slightly changed. In summary, we obtain a second order inexact Galerkin approximation of Maxwell’s equations with block diagonal mass matrix with the same accuracy and flexibility of standard finite element approximations.
The focus of this manuscript lies on the second order approximations for Maxwell’s equations, but the basic arguments can in principle also be used for the construction and analysis of mass-lumping schemes for other equations and approximations of higher order. Some ideas in these directions will be discussed at the end of the manuscript. Let us note that some additional degrees of freedom are required for mass-lumping, whose number increases with higher order of approximation. We therefore expect that discontinuous Galerkin methods, which also have (block) diagonal mass matrices, become more efficient for higher polynomial degree. A thorough comparison of finite elements with mass-lumping and discontinuous Galerkin schemes is given in  in the context of elastodynamics.
The remainder of the manuscript is organized as follows: In Section 2, we briefly summarize some results about the discretization of electromagnetic wave propagation problems by inexact Galerkin methods in space and explicit time integration scheme. A convergence analysis is given under some simple abstract conditions, which can easily be verified for particular approximations. As an example, in Section 2.3, we apply the results to the standard element. In Section 3.1, we then analyze the effect of inexact numerical integration for the element. Sections 3.2 and 3.3 then contain our main results: We first show that the element leads to second order convergence, if , and then propose and analyze the new element, which leads to second order convergence in the general case. Some numerical tests are presented in Section 5 for illustration of our theoretical results. In Section 6, we briefly review the main ingredients that are required to obtain higher order approximations or discretizations for other types of equations. Detailed proofs for some technical lemmas and a list of basis functions for the element are provided in the appendix.
2. Inexact Galerkin approximations
We consider Maxwell’s equations
on some bounded Lipschitz domain . For ease of presentation, we complement (1) by homogeneous boundary and initial conditions
More general boundary and initial conditions and also lower order terms in (1) can be considered with minor modifications. To avoid technicalities, we assume that
|, are symmetric positive definite matrices, and that||(4)|
Piecewise smooth coefficient functions and more general right hand sides can again be considered with minor modifications. Under assumptions (4)–(5), the existence of a unique solution to (1)–(3) can be proven by semi-group theory [14, 18], and solutions of (1) are characterized by the variational principle
for all and ; see  for details. Here and below, denotes the scalar product of .
2.1. Space discretization
Let denote a shape-regular tetrahedral mesh of the domain . We denote by the spaces of piecewise smooth functions over the mesh , and write
for the space of piecewise polynomial functions whose restrictions to any element belong to the Nédélec space ; see . We look for approximations for in a finite element space satisfying
There exist locally defined projection operators such that
Due to shape-regularity of the mesh the constant can be chosen independent of the element . We further denote by the -projection operator to piecewise polynomials of order and note that
The definition of the
-projection naturally extends to vector valued functions.
Find such that and
for all and . Here denotes a suitable approximation for the scalar product on , which is part of the definition of the method.
In order to ensure the well-posedness of the semi-discrete problem, we require that
the bilinear form defines a scalar product on and
for some positive constants .
As a consequence of this assumption and the positivity of , we can estimate
By choosing any basis for the finite dimensional space , we can transform the discrete variational equation (10) into a linear system
describing the evolution of the coordinate vector representing the finite element function . Due to condition (A1), the mass matrix in (12) is symmetric positive definite, and existence of a unique solution for given initial values follows from the Picard-Lindelöf theorem. This also implies the well-posedness of Problem 2.1. As a second ingredient, we assume that
the inexact scalar product is sufficiently accurate, i.e. there exists a constant such that
where represents the quadrature error and is a suitable projection operator satisfying
Based on these general assumptions, we can prove the following convergence result.
The result follows from standard energy arguments and assumptions (A0)–(A2). For convenience of the reader, a complete derivation is given in the appendix.
2.2. Time discretization
Let us also consider the time-discretization by a typical explicit scheme, which will be used in our numerical tests.
Problem 2.3 (Fully discrete scheme).
Set , and then determine for by the variational equations
Here is the approximation for the semi-discrete solution at time resulting from time discretization, and is the time step size. Furthermore,
are the usual central difference quotients of first and second order. Moreover, let
In order to ensure the stability of the fully discrete scheme, we require that
the time step is chosen such that
which can be interpreted as an abstract CFL condition; see  for details. The following error estimates can then be proven via energy arguments.
A detailed proof of this result will again be given in the appendix.
2.3. Standard Nédélec elements
As inexact scalar product for Problem 2.1, we choose
with evaluated by appropriate numerical quadrature.
It suffices to verify assumption (A2) with and . Let denote the -projection onto . Since is constant, , and the local quadrature rule is exact for polynomials of degree , we have
This already concludes the proof of the lemma. ∎
If the quadrature formula is exact for polynomials of degree , then the quadrature error is zero and the method of Problem 2.1 coincides with the standard finite element approximation of second order , which is also included in our analysis. The previous lemma shows that some amount of inexact numerical integration is allowed without degrading the second order convergence.
3. Second order finite elements with mass-lumping
three degrees of freedom are required for every quadrature point;
sufficiently many quadrature have to be located at the boundary in order to allow for appropriate continuity of the associated basis functions.
where are the vertices and are the face midpoints of the tetrahedron ; see Figure 1 for a graphical illustration. By elementary computations, one can verify
The quadrature rule (19) is exact for polynomials of degree .
From Lemma 2.5, we infer that the Nédélec elements of second order with (19) as quadrature rule yield second order convergence. To satisfy condition , Elmkies and Joly proposed an extension of the space by four additional basis functions. In the following two sections, we analyze the approach of  and show that the method is second order convergent only in particular cases. We then propose a modification that yields second order convergence in general.
3.1. First order convergence of the element
We start by defining the extension of the Nédélec finite element space proposed in , which we call the element in the sequel. Let , denote the barycentric coordinates of the tetrahedron , and consider the four functions
with in circular permutation. Each of these functions can be associated to the face opposite to the vertex ; see Figure 1 for an illustration. Note that and has zero tangential trace on and, therefore, its extension by zero lies in . Following [4, 10], we define
and we introduce the corresponding global approximation space
By elementary arguments, we can verify the following properties.
Assumptions (A0) and (A1) hold by construction. To verify (A2) with and , we choose and use that and Lemma 3.1. This yields . ∎
Let us note that the additional basis functions (20) are cubic polynomials, which forced us to use in the estimate of the quadrature error. As we will indicate by numerical tests, the assertions of Lemma 3.2 are sharp, i.e., in general, the method only provides first order convergence. In the following section, we will prove, however, that second order convergence can be obtained in special situations.
3.2. Second order convergence for the element
Let us start by summarizing some additional properties of the finite element space.
This means that the four additional basis functions , are independent but one specific linear combination leads to a curl-free function. Since
the sum of the basis functions is the problematic linear combination, since it is the gradient of the bubble function . Based on this observation and Lemma 3.3, we can split any function into
with , and . Moreover, this splitting is unique and direct, which allows to prove the following assertions.
Let be split into as above. Then
with a constant depending only on the shape of the element . For the second component, we further have
The first two estimates follow from linear independence of the basis functions and a mapping argument. The last is based on the fact that the curls of the three basis functions spanning these components are linearly independent. As a consequence, defines a norm on the space , and one can see that , hence also defines a norm on this three dimensional subspace. The fact that and can be chosen depending only on the shape of the element follows from the usual mapping argument and the uniform shape regularity of the mesh. ∎
The above splitting generalizes directly to the whole space, i.e., any function defined in (21) can be split uniquely into
with components given by
Since we assumed uniform shape regularity of the mesh, we can take the same constant in Lemma 3.4 and thus obtain global estimates
which follow by summation of the local estimates of the previous lemma.
Moreover, the following approximation error estimates hold
By we define the corresponding global projection operator for piecewise smooth functions . With the help of this projection operator and the splitting of the test space , we can now establish the following improved estimates for the quadrature error.
We choose and let be the splitting of defined above. Then
Due to the exactness of the quadrature rule stated in Lemma 3.1, we have . For the second term, we get
By Lemma 3.1, the last two terms vanish identically, and we obtain
where we used the approximation properties of and the third estimate of Remark 3.5 in the last step. This is the required estimate for the second component. From the assumption and the divergence-preserving property of the projection operator , we infer that
A close inspection of the quadrature rule (19) shows that
and as defined above, i.e. the quadrature rule also integrates one additional forth order polynomial exactly. As a consequence, , and the proof is concluded by summing up the estimates for the terms –. ∎
The element with mass-lumping thus yields second order convergence provided that the exact solution is divergence free. From (1), we see that
3.3. A modification of the element
From the error analysis presented in the previous section, one can see that the problematic component in the estimate for the quadrature error is that related with the test function which cannot be controlled via its curl nor is integrated with sufficient accuracy in the general case. We therefore replace the basis functions in the element by
and define the modified element by
By elementary computations, one can again verify the following properties.
and , i.e., the four additional functions have linear independent curls. Moreover, is integrated exactly by the quadrature rule (19), i.e., for all .
Based on the previous results, we can prove the following assertion.
Conditions (A0) and (A1) follow by construction. For the proof of property (A2), we note that can now be split into with components
We can then continue as in the proof of Lemma 3.6 without taking the third component in the decomposition of the quadrature error into consideration. ∎
In the next section, we define a local basis for the space , which together with the quadrature rule (19) leads to a block-diagonal mass matrix . We thus obtain an efficient method of second order accuracy.
4. Definition of the basis functions
We start by defining basis functions for the standard Nédélec space . We only consider a single element and omit the corresponding subscript in the following. Let , denote the barycentric coordinates of the element associated to the vertices , which are ordered with respect to their global index in the mesh. Furthermore, let , denote the midpoint of the face opposite to the vertex , see Figure 1. For each edge spanned by vertices and with we define two basis functions of the form leading to
Any basis function associated to an edge has non-zero tangential trace only on the respective edge and vanishes in all but one vertex. For a face with vertices , , , , , we define two basis functions of the form , viz.
These functions vanish identically at all vertices and on one of the faces. The functions form a basis for the Nédélec space . Note that three basis functions can be associated to each vertex while only two functions are associated to every face midpoint; see Figure 1 for an illustration.