A second order finite element method with mass lumping for Maxwell's equations on tetrahedra

02/13/2020 ∙ by Herbert Egger, et al. ∙ Technische Universität Darmstadt 0

We consider the numerical approximation of Maxwell's equations in time domain by a second order H(curl) conforming finite element approximation. In order to enable the efficient application of explicit time stepping schemes, we utilize a mass-lumping strategy resulting from numerical integration in conjunction with the finite element spaces introduced by Elmkies and Joly. We prove that this method is second order accurate if the true solution is divergence free, but the order of accuracy reduces to one in the general case. We then propose a modification of the finite element space, which yields second order accuracy in the general case.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

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, 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 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 [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, 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 [5]. 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 [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 mass-lumping 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 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

[12], 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 [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


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 [16] 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 [1]. 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.

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


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.

Theorem 2.2.

Let (A0)–(A2) hold and let denote a sufficiently smooth solution of (1)–(3). Then the inexact Galerkin approximation of Problem 2.1 satisfies


with rate and constant

and constant only depending on the domain, the shape-regularity of the mesh, and the bounds for the coefficients.

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 [13] for details. The following error estimates can then be proven via energy arguments.

Theorem 2.4.

Let (A0)–(A3) hold and let denote a sufficiently smooth solution of (1)–(3). Then the discrete approximations , of Problem 2.3 satisfy

for all with rate and constant from Theorem 2.2 and

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


As inexact scalar product for Problem 2.1, we choose


with evaluated by appropriate numerical quadrature.

Lemma 2.5.

Assume that is exact for polynomials of degree and chosen such that (A1) is valid. Then assumption (A2) holds with and . As a consequence, the estimates of Theorems 2.2 and 2.4 hold with rate , i.e., the method is second order accurate.


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 mass-lumping

As observed in [4, 10], mass-lumping 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]


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


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.

Lemma 3.2.

Let be defined by the quadrature rule (19). Then assumptions (A0)–(A2) hold with , , and . As a consequence, the estimates of Theorems 2.2 and 2.4 hold with , i.e., the method is first order accurate.


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 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.

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


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.

Let . Then assumptions (A0)-(A2) hold with , , and some depending only on the shape regularity of the mesh. As a consequence, the estimates of Theorems 2.2 and 2.4 hold with rate , i.e., the method is second order accurate if the exact solution is divergence free.


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 . ∎

Remark 3.7.

The element with mass-lumping 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


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.

Assumptions (A0)–(A2) hold with , , and only depending on the shape-regularity of the mesh. Hence the assertions of Theorems 2.2 and 2.4 hold with rate , i.e., the method is second order accurate.


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.