1. Introduction
We consider the efficient numerical solution of electromagnetic wave propagation modeled by timedependent Maxwell’s equations in second order form
Here is the electric field, and
are the symmetric and positive definite permittivity and permeability tensors, and
is 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 timedomain 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 timestepping schemes, they lead to very efficient and accurate numerical approximations. Nontrivial modifications are, however, required to guarantee stability of the schemes in the case of nonrectangular grids and discontinuous or anisotropic coefficients [19], 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, masslumping 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 masslumping can be interpreted as inexact numerical integration. In this spirit, masslumping for finite element methods for Maxwell’s equations on quadrilateral and hexahedral grids has been investigated in [5]. In fact, a close relation exists between finite difference schemes [20, 21] and low order finite element approximations with masslumping; we refer to [4, 13] for details.
In order to obtain the full geometric flexibility of finite element approximations, we here consider masslumping 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 [9] and first order convergence has been established. Related methods have been proposed in [3] in the context of the finite integration technique, but no convergence analysis is given there. A masslumping strategy based on an extension of the lowest order elements has been proposed by Elmkies and Joly [10] 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 masslumping 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 masslumping, 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 masslumping 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 masslumping, whose number increases with higher order of approximation. We therefore expect that discontinuous Galerkin methods
[12], which also have (block) diagonal mass matrices, become more efficient for higher polynomial degree. A thorough comparison of finite elements with masslumping and discontinuous Galerkin schemes is given in [11] 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
(1) 
on some bounded Lipschitz domain . For ease of presentation, we complement (1) by homogeneous boundary and initial conditions
(2)  
(3) 
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)  
(5) 
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 semigroup theory [14, 18], and solutions of (1) are characterized by the variational principle
(6) 
for all and ; see [16] for details. Here and below, denotes the scalar product of .
2.1. Space discretization
Let denote a shaperegular tetrahedral mesh of the domain . We denote by the spaces of piecewise smooth functions over the mesh , and write
(7) 
for the space of piecewise polynomial functions whose restrictions to any element belong to the Nédélec space ; see [1]. We look for approximations for in a finite element space satisfying

.
There exist locally defined projection operators such that
(8)  
(9) 
Due to shaperegularity 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.
For the numerical approximation of problem (1)–(3), we then consider inexact Galerkin finite element methods of the following form.
Problem 2.1.
Find such that and
(10) 
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 wellposedness of the semidiscrete 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
(11) 
By choosing any basis for the finite dimensional space , we can transform the discrete variational equation (10) into a linear system
(12) 
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 PicardLindelöf theorem. This also implies the wellposedness 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.
Theorem 2.2.
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 timediscretization 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
(14) 
Here is the approximation for the semidiscrete solution at time resulting from time discretization, and is the time step size. Furthermore,
(15) 
are the usual central difference quotients of first and second order. Moreover, let
(16) 
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 [13] for details. The following error estimates can then be proven via energy arguments.
Theorem 2.4.
A detailed proof of this result will again be given in the appendix.
2.3. Standard Nédélec elements
To illustrate the applicability of the convergence results above, let us briefly discuss the approximation of (1)–(3) using second order Nédélec elements and inexact numerical integration. We set
(17) 
As inexact scalar product for Problem 2.1, we choose
(18) 
with evaluated by appropriate numerical quadrature.
Lemma 2.5.
Proof.
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. ∎
Remark 2.6.
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 [16], 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 masslumping
As observed in [4, 10], masslumping via numerical quadrature relies on the following key ingredients to be satisfied on every element :

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.
We refer to Section 4 for details. An appropriate quadrature rule is given by [10]
(19) 
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
Lemma 3.1.
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 [10] 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 [10], which we call the element in the sequel. Let , denote the barycentric coordinates of the tetrahedron , and consider the four functions
(20) 
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
(21) 
By elementary arguments, we can verify the following properties.
Lemma 3.2.
Proof.
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.
Lemma 3.3.
and .
This means that the four additional basis functions , are independent but one specific linear combination leads to a curlfree 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.
Lemma 3.4.
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
Proof.
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. ∎
Remark 3.5.
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.
In the sequel, we utilize an additional divergence preserving projection operator. Let be the canonical projection operator for the space , see [1, Sec 2.5] or [2], and recall that
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.
Lemma 3.6.
Proof.
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 divergencepreserving 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 –. ∎
Remark 3.7.
The element with masslumping thus yields second order convergence provided that the exact solution is divergence free. From (1), we see that
Hence the assumption holds, if and for all . In this case, Lemma 3.6 guarantees second order convergence, which explains the good numerical results obtained in [4, 10].
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
(22) 
and define the modified element by
By elementary computations, one can again verify the following properties.
Lemma 3.8.
and , i.e., the four additional functions have linear independent curls. Moreover, is integrated exactly by the quadrature rule (19), i.e., for all .
As approximation space for Problems 2.1 and 2.3, we now consider
Based on the previous results, we can prove the following assertion.
Lemma 3.9.
Proof.
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 blockdiagonal 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 nonzero 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.
Comments
There are no comments yet.