Partially discontinuous nodal finite elements for $H(\operatorname{curl})$ and $H(\operatorname{div})$

We assume that $\Omega$ is a polygonal or polyhedral domain. Given a mesh, we will use $\mathcal{V}$ to denote the set of vertices, $\mathcal{E}$ for edges, $\mathcal{F}$ for faces, and $\mathcal{T}$ for the 3D simplices. We use $V$, $E$, $F$, and $T$ to denote the number of vertices, edges, faces, and tetrahedra, respectively. From Euler’s formula, one has $V-E+F=1$ in 2D and $V-E+F-T=1$ in 3D for contractible domains.

For a triangle $f$, we denote its Clough-Tocher split as $f_{\operatorname{CT}}$, which is the collection of three triangles obtained by connecting the vertices to an interior point of $f$. For a tetrahedron $K$ in $\mathcal{T}$, the Worsey-Farin split of $K$, denoted by $K_{\operatorname{WF}}$, is obtained by adding a vertex to the interior of $K$, connecting this vertex to the vertices of $K$; then adding an interior vertex to each face of $K$, and connecting this point to the vertices of the face and the interior point (see Figure 1). Therefore, each face is divided into three triangles like the Clough-Tocher split. The set of the Worsy-Farin splits of all the tetrahedra in $\mathcal{T}$ is referred to as the Worsey-Farin mesh and denoted by $\mathcal{T}_{\operatorname{WF}}$. We refer to [21, 23, 32] for more details.

We use $\bm{\nu}_{f}$ and $\bm{\tau}_{f}$ to denote the unit normal and tangential vectors of a simplex $f$, respectively. In 2D, the tangential and normal directions of an edge are uniquely defined up to an orientation. For edges in 3D there are one unit tangential and two linearly independent unit normal directions, and for faces in 3D there are one unit normal and two linearly independent unit tangential directions. We will write $\bm{\tau}_{e}$, $\bm{\nu}_{e,i}$ and $\bm{\nu}_{f}$, $\bm{\tau}_{f,i}$, $i=1,2$ for these cases, respectively.

In our discussions, $C^{r}(\mathcal{V})$ includes piecewise smooth functions for which derivatives up to order $r$ coincide at vertices. Similarly we can define $C^{r}(\mathcal{E})$ and $C^{r}(\mathcal{F})$, $C^{r}(\mathcal{T})$ for the space of functions with $C^{r}$ continuity on the edges, faces, and 3D cells, repectively.

In $d$D ($d=2,3$), we define the gradient operator $\operatorname{grad}=\nabla$ and the $\operatorname{div}$ operator $\operatorname{div}=\nabla\cdot$. In 3D, for $\bm{u}=(u_{1},u_{2},u_{3})^{\mathrm{T}}$, we define

$$\operatorname{curl}\bm{u}=(\partial_{x_{2}}u_{3}-\partial_{x_{3}}u_{2},\partial_{x_{3}}u_{1}-\partial_{x_{1}}u_{3},\partial_{x_{1}}u_{2}-\partial_{x_{2}}u_{1})^{\mathrm{T}}.$$

In 2D, we define

$$\operatorname{rot}\bm{u}=\partial_{x_{1}}u_{2}-\partial_{x_{2}}u_{1}\text{ for }\bm{u}=(u_{1},u_{2})^{\mathrm{T}}\text{ and }\operatorname{curl}u=(\partial_{x_{2}}u,-\partial_{x_{1}}u)\text{ for a scalar }u.$$

We review some basic facts about the de Rham complex; further details can be found, for instance, in [5, 13]. The de Rham complex consists of differential forms and exterior derivatives. In 3D, the de Rham complex on $\Omega$ reads

(2.1)

with $\operatorname{curl}\operatorname{grad}=\operatorname{div}\operatorname{curl}=0$. A complex is called exact at a space if the kernel of the subsequent operator is identical to the range of the previous one. For example, we say that the complex (2.1) is exact at $H(\operatorname{curl};\Omega)$ if $\mathcal{N}(\operatorname{curl},H(\operatorname{curl}))=\mathcal{R}(\operatorname{grad},H^{1})$. Here for a linear operator $d:X\to Y$, $\mathcal{N}(d,X)\subset X$ denotes the kernel of $d$ and $\mathcal{R}(d,Y)\subset Y$ denotes the range. The complex (2.1) is called exact if it is exact at all the spaces.

If there are finite element spaces $\Sigma_{h}\subset H^{1}(\Omega)$, $V_{h}\subset H(\operatorname{curl};\Omega)$, $X_{h}\subset H(\operatorname{div};\Omega)$, and $W_{h}\subset L^{2}(\Omega)$ such that

(2.2)

is a complex, we say that (2.2) is a finite element subcomplex of (2.1). In this case, a necessary condition for the exactness is the following dimension condition

(2.3)

$\displaystyle\dim(\mathbb{R})-\dim\Sigma_{h}+\dim V_{h}-\dim X_{h}+\dim W_{h}=0.$

As an example, the Lagrange element space, Nédélec element space of the first kind [33], Raviart-Thomas (RT) element space [33], and discontinuous element space fit into a discrete de Rham complex. The Lagrange element space, Nédélec element space of the second kind [34], Brezzi-Douglas-Marini (BDM) element space [34], and discontinuous element space also fit into a discrete de Rham complex. For a detailed discussion of these finite element spaces, see, e.g., [9].

We use $\mathcal{P}_{p}$ to denote the polynomial space of degree $p$ (a nonnegative integer). Denote

$$\operatorname{NED}_{p}^{1}=(\mathcal{P}_{p-1})^{3}+S_{p}\text{ with }S_{p}=\{q\in(\mathcal{P}_{p})^{3}:q\cdot\bm{x}=0\}$$

as the shape function spaces of the first-kind Nédélec element with polynomial degree $p$, and

$$\operatorname{RT}_{p}=(\mathcal{P}_{p-1})^{3}+\bm{x}\mathcal{P}_{p-1},$$

as the shape function space of the RT element with the polynomial degree $p$. The shape function space of the second-kind Nédélec element and the BDM element with the polynomial degree $p$ are

$$\operatorname{NED}_{p}^{2}=[\mathcal{P}_{p}]^{3}\text{ and }\operatorname{BDM}_{p}=[\mathcal{P}_{p}]^{3},$$

respectively. We denote the Nédélec finite element space of the first and second kind with the polynomial degree $p$ on mesh $\mathcal{T}$ as $\operatorname{NED}_{p}^{1}(\mathcal{T})$ and $\operatorname{NED}_{p}^{2}(\mathcal{T})$, respectively.

In 2D, the de Rham complex reads

or

In this section, we briefly review the constructions in 2D in [16] and some questions that were left open there.

In 2D, the Stenberg $H(\operatorname{div})$ element can be obtained by breaking the tangential component of a vector-valued Lagrange element on each interior edge of the mesh, i.e., the tangential trace of the functions in the Stenberg $H(\operatorname{div})$ element space from two sides of an edge can be different. Therefore, the continuity of the Stenberg element is in between the corresponding vector-valued $C^{0}$ Lagrange element and the discontinuous element (see Figure 2). The Stenberg element fits in a discrete de Rham complex in 2D, which begins with the Hermite element taking vertex derivatives as degrees of freedom [16].

The 2D Stenberg $H(\operatorname{div})$ element can be extended to 3D (see [16, Section 3.5] and $\mathcal{W}_{p+1}^{2}$ in (4.1)). However, the complex containing it was not discussed in [16].

The construction in [16] contains two complexes in 3D: one begins with the Hermite element, and another begins with a $C^{0}$ element with second order supersmoothness at the vertices (formally, $C^{2}(\mathcal{V})$, although the function is not $C^{2}$). The latter can be viewed as a 3D extension of the Falk-Neilan complex for the Stokes problem [19] in terms of supersmoothness. The $H(\operatorname{curl})$ element in the former complex does not admit nodal basis, while the $H(\operatorname{curl})$ element in the latter complex can be represented in terms of scalar Hermite bases. Nevertheless, it was left open in [16] to construct a sequence containing an $H(\operatorname{curl})$ element with Lagrange bases. We solve these problems in this paper.

A natural candidate for a partially discontinuous nodal $H(\operatorname{curl})$ element is obtained by breaking the normal components of the vector-valued Lagrange finite element on each face (at the same time, the $C^{0}$ continuity at vertices and on edges is retained):

$$Z_{h}^{p}:=\operatorname{NED}_{p}^{2}(\mathcal{T})\cap C^{0}(\mathcal{V})\cap C^{0}(\mathcal{E}).$$

This element is in between the $C^{0}$ vector-valued Lagrange element and the Nédélec element in terms of continuity. Nevertheless, this element unlikely fits in a finite element de Rham complex due to the following observation. In most existing cases, the “restriction” (precise definition discussed in [17]) of a higher-dimensional finite element (complex) to a lower-dimensional cell is a lower-dimensional finite element (complex). This holds for the standard finite element differential forms [3], as well as all the examples in [16]. Therefore, to construct a complex in 3D, we want the restrictions of the involved spaces and the degrees of freedom on each face to form a 2D complex as well. Nevertheless, this is unlikely for the case of $Z_{h}^{p}$. Because the restriction of $Z_{h}^{p}$ on each face is a 2D Lagrange element, which unlikely fits in a de Rham complex, as reflected in the deficiency of $Z_{h}^{p}$ (at least for low order cases) as a discretization for the Stokes problem, c.f., [9, 24, 39].

Although the above observation already suggests that $Z_{h}^{p}$ may not be a good candidate for an $H(\operatorname{curl})$ nodal element, one may still ask if $Z_{h}^{p}$, as a space in between the vector-valued Lagrange element (which produces spurious solutions) and the Nédélec element (which excludes spurious solutions), suffers from spurious solutions for eigenvalue problems. To the best of our knowledge, the emergence of spurious solutions for eigenvalue problems is not predicted by theory in most cases. Therefore we answer this question with numerical tests.

We apply the $H(\operatorname{rot})$ version of the Stenberg element in 2D defined by

$$Y_{h}^{p}:=\operatorname{NED}_{p}^{2}(\mathcal{T})\cap C^{0}(\mathcal{V}),$$

and $Z_{h}^{p}$ in 3D to solve the following Maxwell eigenvalue problem on the domain $\Omega=(0,\pi)^{d}$, $d=2,3$:

(3.1)		$\displaystyle\operatorname{curl}\operatorname{curl}\bm{u}=\lambda\bm{u}\ \text{in }\Omega,$
		$\displaystyle\operatorname{div}\bm{u}=0\text{ in }\Omega,$
		$\displaystyle\bm{u}\times\bm{n}=0\text{ on }\partial\Omega.$

In 2D, the eigenvalues of this problems are $m^{2}+n^{2}$ [3], and the corresponding eigenvectors are

$$\bm{u}(x_{1},x_{2})=\left(\begin{matrix}n\sin(mx_{1})\cos(nx_{2})\\ -m\cos(mx_{1})\sin(nx_{2})\end{matrix}\right),$$

where $m,n$ are nonnegative integers, not both 0. In 3D, the eigenvalues of this problems are $m^{2}+n^{2}+l^{2}$, and the corresponding eigenvectors are

$$\bm{u}(x_{1},x_{2},x_{3})=\left(\begin{matrix}(l-n)\cos(mx_{1})\sin(lx_{3})\sin(nx_{2})\\ (m-l)\cos(nx_{2})\sin(mx_{1})\sin(lx_{3})\\ (n-m)\cos(lx_{3})\sin(mx_{1})\sin(nx_{2})\end{matrix}\right),$$

where $m,n,l$ are nonnegative integers, not all 0.

We first apply the Stenberg $H(\operatorname{rot})$ finite element space $Y_{h}^{p}$ to solve the problem (3.1) in 2D. The numerical results are shown in Table 1.

From the table, we see that spurious solutions do not appear. Therefore the vector-valued Lagrange elements are remedied by breaking the normal degrees of freedom on edges. Then we use $Z_{h}^{p}$ to solve the problem (3.1) for $d=3$. Table 2 shows the numerical eigenvalues around 3. In this case, spurious solutions occur. Therefore relaxing the normal continuity on faces (equivalently, enriching face-based bubbles) on the vector-valued Lagrange element is not enough to give a proper discretization of the Maxwell equations.

Now both heuristic argument and numerical experiments inspire us to seek a different construction other than $Z_{h}^{p}$. We cannot further relax the vertex and edge continuity, which is required to ensure the existence of nodal bases. Motivated by the heuristic argument on the face modes, now we plan to construct an element such that its restriction on each face is a proper velocity finite element for the Stokes problem. Recall that the $\mathcal{P}_{2}$-$\mathcal{P}_{1}$ Scott-Vogelius pair on the Clough-Tocher split of triangular meshes (see Figure 3) leads to a complex and is thus stable [7, 36]. This motivates us to use the Worsey-Farin split of a tetrahedron [32, 45], which, restricted on each face of the macro tetrahedron, is a Clough-Tocher split.

We consider a complex for $p\geq 3$ as follows:

(3.2)

where each space is constructed on the Worsey-Farin split. The Worsey-Farin $C^{1}$ element [45] requires a mesh alignment, the subdivision point on each face of the macro tetrahedron lies on the line connecting the subdivision points of the two macro tetrahedra sharing the face. Nevertheless, such a condition is not required in the construction below.