Construction of Supplemental Functions for Direct Serendipity and Mixed Finite Elements on Polygons

Abstract

New families of direct serendipity and direct mixed finite elements on general planar, strictly convex polygons were recently defined by the authors. The finite elements of index r are $H^1$ and $H\left(\mathrm{div}\right)$ conforming, respectively, and approximate optimally to order $r+1$ while using the minimal number of degrees of freedom. The shape function space consists of the full set of polynomials defined directly on the element and augmented with a space of supplemental functions. The supplemental functions were constructed as rational functions, which can be difficult to integrate accurately using numerical quadrature rules when the index is high. This can result in a loss of accuracy in certain cases. In this work, we propose alternative ways to construct the supplemental functions on the element as continuous piecewise polynomials. One approach results in supplemental functions that are in $H^p$ for any $p\ge 1$. We prove the optimal approximation property for these new finite elements. We also perform numerical tests on them, comparing results for the original supplemental functions and the various alternatives. The new piecewise polynomial supplements can be integrated accurately, and therefore show better robustness with respect to the underlying meshes used.

1. Introduction

There has been strong interest in using polygonal and polyhedral meshes when solving certain types of problems via the finite element method. For just a few examples, we note problems in solid mechanics [1,2], elasticity [3,4], fracture mechanics [5,6,7], thin plates [8], shells [9], porous media [10], topology optimization [11,12,13], and finding eigenvalues [14]. In fact, polygonal meshes are an important motivation for the development and use of methods beyond the classic finite element method, which include, for example, the discontinuous Galerkin methods (including weak Galerkin [15] and ultra-weak methods [16,17,18]), mimetic methods [19,20,21], and virtual element methods [22,23,24,25].

Classic conforming finite element methods have also been developed for use on polygonal meshes, and especially for quadrilateral meshes. Approaches taken include the use of maps from reference finite elements [26,27,28], restriction to low order elements [29,30,31,32], the use of macro-elements [33], basis function enrichment [34,35,36], and construction using barycentric coordinates [9,37,38,39]. Ideally, we would have families of conforming finite elements defined for any order of accuracy. These would possess a minimal number of degrees of freedom (DoFs) subject to conformity and accuracy constraints. Finite elements based on the use of non-affine maps from reference finite elements display degraded accuracy. Accuracy is restricted if only low order elements are defined. Macro-elements, basis function enrichment, and the use of barycentric coordinates in higher order cases results in finite elements with an excess number of DoFs.

Families of conforming finite elements defined on polygons that maintain both accuracy and a minimal number of DoFs have appeared recently [40,41,42,43] (as well as some finite elements in three dimensions [44,45,46]). The approach taken is to begin with the space of polynomials ${\mathbb{P}}_r\left(E\right)$ of degree up to r defined directly on the physical element E to achieve accuracy of order $r+1$. To achieve conformity, one then adds in a space of supplemental functions. A basis for the supplemental functions must have certain properties on $\partial E$, but they must be defined over all of E by filling in the interior. The “supplemental function space” is sometimes called the “filling space”.

In this paper, we discuss the construction of the supplemental functions, in the context of the finite elements developed by the current authors in [43], which are called direct finite elements. Let the element $E=E_N\subset {\mathbb{R}}^2$ be a closed, nondegenerate, convex polygon with $N\ge 3$ edges. The direct serendipity finite elements of index $r\ge 1$ are $H^1$-conforming and take the form

$${\mathcal{DS}}_r\left(E_N\right)={\mathbb{P}}_r\left(E_N\right)\oplus {\mathbb{S}}_r^{\mathcal{DS}}\left(E_N\right),$$

(1)

where ${\mathbb{S}}_r^{\mathcal{DS}}\left(E_N\right)$ is the space of supplemental functions. The direct mixed finite elements are $H\left(\mathrm{div}\right)$-conforming and take two forms,

$$\begin{array}{c}{\mathbf{V}}_r^r\left(E_N\right)={\mathbb{P}}_r^2\left(E_N\right)\oplus \mathbf{x}{\tilde{\mathbb{P}}}_r\left(E_N\right)\oplus {\mathbb{S}}_r^{\mathbf{V}}\left(E_N\right),\\ {\mathbf{V}}_r^{r-1}\left(E_N\right)={\mathbb{P}}_r^2\left(E_N\right)\oplus {\mathbb{S}}_r^{\mathbf{V}}\left(E_N\right),\end{array}$$

(2)

for full ($r\ge 0$) and reduced ($r\ge 1$) $H\left(\mathrm{div}\right)$-approximation, respectively, where ${\tilde{\mathbb{P}}}_r\left(E_N\right)$ are the homogeneous polynomials of (exact) degree r. These two finite elements are related to each other by the finite element exterior calculus [47] through the de Rham complex

$$\mathbb{R}\hookrightarrow H^1\stackrel{\mathrm{curl}}{\to }H\left(\mathrm{div}\right)\stackrel{\mathrm{div}}{\to }{L}^2⟶0,$$

(3)

resulting in, for $s=r-1,r(s\ge 0)$,

$$\mathbb{R}\hookrightarrow {\mathcal{DS}}_{r+1}\left(E_N\right)\stackrel{\mathrm{curl}}{\to }{\mathbf{V}}_r^s\left(E_N\right)\stackrel{\mathrm{div}}{\to }{\mathbb{P}}_s\left(E_N\right)⟶0.$$

(4)

The consequence is that

$$\begin{array}{c}{\mathbf{V}}_r^r\left(E_N\right)=\mathrm{curl}{\mathcal{D}\mathcal{S}}_{r+1}\left(E_N\right)\oplus \mathbf{x}{\mathbb{P}}_r,\\ {\mathbf{V}}_r^{r-1}\left(E_N\right)=\mathrm{curl}{\mathcal{DS}}_{r+1}\left(E_N\right)\oplus \mathbf{x}{\mathbb{P}}_{r-1},\end{array}$$

(5)

and, therefore,

$${\mathbb{S}}_r^{\mathbf{V}}\left(E_N\right)=\mathrm{curl}{\mathbb{S}}_{r+1}^{\mathcal{DS}}\left(E_N\right).$$

(6)

The original construction of supplemental functions made use of rational functions (see (21)), which are difficult to numerically integrate accurately. As a consequence, when solving a partial differential equation using direct finite elements, the quadrature error may be significant, leading to poor overall approximation of the solution. This was observed in [43], although in that paper, the degradation in the approximation was attributed to poor mesh quality. While mesh quality remains an important ingredient in finite element analysis, quadrature approximation is also a critical component, especially for high order methods.

In this work, we introduce two constructions of the supplemental functions ${\mathbb{S}}_r^{\mathcal{DS}}\left(E_N\right)$ which involve using continuous piecewise polynomials. Such constructions are motivated by the work of Kuznetsov and Repin [33], and suggested by the work of Cockburn and Fu [41]. These new supplemental functions can then be accurately integrated by quadrature rules. (A similar, but more complex, construction in three dimensions for cuboidal hexahedra is discussed in [46]).

In the next two sections we present some basic notation and review the general definition of the original direct serendipity and direct mixed finite elements, which have supplemental functions that are $C^{\infty }$-smooth. Our new families of direct finite elements based on piecewise continuous supplemental functions are given in Section 4. We give two constructions of the supplemental functions, so that one set lies in $H^1$ and the other in $H^p$ for any integer $p\ge 1$. The approximation properties of these new direct finite elements are given in Section 5. The results are optimal, up to the bounding constant. The proof follows that given in [43], and we concentrate on the modifications that are required to handle the new supplements. In Section 6, we present numerical tests that compare the errors and convergence rates of the new and original direct finite elements. We conclude the paper in Section 7.

2. Notation

We choose to identify the edges and vertices of $E_N$ adjacently in the counterclockwise direction, as depicted in Figure 1 (throughout the paper, we interpret indices modulo N). Let the edges of $E_N$ be denoted $e_i$, $i=1,2,\dots ,N$, and the vertices be ${\mathbf{x}}_{v,i}=e_i\cap e_{i+1}$. Let ${\nu }_i$ denote the unit outer normal to edge $e_i$, and let ${\tau }_i$ denote the unit tangent vector of $e_i$ oriented in the counterclockwise direction, for $i=1,2,\dots ,N$.

For any two distinct points ${\mathbf{y}}_1$ and ${\mathbf{y}}_2$, let $\mathcal{L}[{\mathbf{y}}_1,{\mathbf{y}}_2]$ be the line passing through ${\mathbf{y}}_1$ and ${\mathbf{y}}_2$, and take $\nu [{\mathbf{y}}_1,{\mathbf{y}}_2]$ to be the unit vector normal to this line interpreted as going from ${\mathbf{y}}_1$ to ${\mathbf{y}}_2$ and then spinning 90 degrees in the clockwise direction (i.e., pointing to the right). Then we define a linear polynomial giving the signed distance of $\mathbf{x}$ to $\mathcal{L}[{\mathbf{y}}_1,{\mathbf{y}}_2]$ as

$$\lambda [{\mathbf{y}}_1,{\mathbf{y}}_2]\left(\mathbf{x}\right)=-(\mathbf{x}-{\mathbf{y}}_2)·\nu [{\mathbf{y}}_1,{\mathbf{y}}_2].$$

(7)

To simplify the notation for linear functions that will be used throughout the paper, let ${\mathcal{L}}_i=\mathcal{L}[{\mathbf{x}}_{v,i-1},{\mathbf{x}}_{v,i}]$ be the line containing edge $e_i$ and let ${\lambda }_i\left(\mathbf{x}\right)$ give the distance of $\mathbf{x}\in {\mathbb{R}}^2$ to edge $e_i$ opposite the normal direction, i.e.,

$$\begin{array}{cc}\hfill {\lambda }_i\left(\mathbf{x}\right)& =\lambda [{\mathbf{x}}_{v,i-1},{\mathbf{x}}_{v,i}]\left(\mathbf{x}\right)=-(\mathbf{x}-{\mathbf{x}}_{v,i})·{\nu }_i,i=1,2,\dots ,N.\hfill \end{array}$$

(8)

These functions are strictly positive in the interior of $E_N$, and ${\lambda }_i$ vanishes on the edge $e_i$.

3. Direct Serendipity and Mixed Finite Elements

The general development of direct serendipity and mixed finite elements is given in [43]. The definition of the supplemental space ${\mathbb{S}}_r^{\mathcal{DS}}\left(E_N\right)$ in (1) is key to the construction. For completeness, we review the definitions of these direct finite elements here.

When $N=3$ (triangles), the direct serendipity supplemental space ${\mathbb{S}}_r^{\mathcal{DS}}\left(E_3\right)$ is empty. When $N\ge 4$ and $1\le r<N-2$, the direct serendipity spaces ${\mathcal{DS}}_r\left(E_N\right)$ are defined as subspaces of ${\mathcal{DS}}_{N-2}\left(E_N\right)$ by the rule

$$\begin{array}{c}\hfill {\mathcal{DS}}_r\left(E_N\right)=\left\{\phi \in {\mathcal{DS}}_{N-2}\left(E_N\right){:\phi |}_e\in {\mathbb{P}}_r\left(e\right)\mathrm{for}\mathrm{all}\mathrm{edges}e\mathrm{of}{E}_N\right\}.\end{array}$$

(9)

Therefore, we only need to understand ${\mathbb{S}}_r^{\mathcal{DS}}\left(E_N\right)$ for $r\ge N-2$ and $N\ge 4$.

To define the supplemental basis functions, two series of choices must be made for each $i,j$ such that $1\le i<j\le N$ and $2\le j-i\le N-2$ (i.e., $e_i$ and $e_j$ are nonadjacent edges). First, as shown in Figure 2, one must choose two distinct points ${\mathbf{x}}_{i,j}^1\in {\mathcal{L}}_i$ and ${\mathbf{x}}_{i,j}^2\in {\mathcal{L}}_j$ that avoid the intersection point ${\mathbf{x}}_{i,j}={\mathcal{L}}_i\cap {\mathcal{L}}_j$, if it exists. Then let

$${\lambda }_{i,j}\left(\mathbf{x}\right)=\lambda [{\mathbf{x}}_{i,j}^1,{\mathbf{x}}_{i,j}^2]\left(\mathbf{x}\right)=-(\mathbf{x}-{\mathbf{x}}_{i,j}^2)·{\nu }_{i,j},{\nu }_{i,j}=\nu [{\mathbf{x}}_{i,j}^1,{\mathbf{x}}_{i,j}^2],$$

(10)

be the linear function associated to the line ${\mathcal{L}}_{i,j}=\mathcal{L}[{\mathbf{x}}_{i,j}^1,{\mathbf{x}}_{i,j}^2]$. Second, $R_{i,j}$ must be chosen to satisfy

$$\begin{array}{cccccc}\hfill R_{i,j}{\left(\mathbf{x}\right)|}_{e_i}& =-1,\hfill & \hfill & \hfill & \hfill R_{i,j}{\left(\mathbf{x}\right)|}_{e_j}& =1.\hfill \end{array}$$

(11)

The supplemental space for $r\ge N-2\ge 2$ is of the form

$$\begin{array}{cc}\hfill {\mathbb{S}}_r^{\mathcal{DS}}\left(E_N\right)& =\mathrm{span}\left\{{\varphi }_{s,i,j}:1\le i<j\le N,2\le j-i\le N-2\right\},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\varphi }_{s,i,j}& =\left(\prod _{k\ne i,j}{\lambda }_k\right){\lambda }_{i,j}^{r-N+2}{R}_{i,j}.\hfill \end{array}$$

(13)

3.1. Direct Serendipity Finite Elements

Every shape function of the direct serendipity finite element ${\mathcal{DS}}_r\left(E_N\right)$ is a sum of a polynomial and a linear combination of the supplemental functions, as in (1). To implement them, one must define the DoFs. For example, for $\psi \in {\mathcal{DS}}_r\left(E_N\right)$, one can take

$$\begin{array}{cccc}\hfill & \psi \left({\mathbf{x}}_{v,i}\right),\hfill & & \forall i=1,2,\dots ,N,\hfill \end{array}$$

(14)

$$\begin{array}{cccc}\hfill & {\int }_{e_i}\psi pd\sigma ,\hfill & & \forall p\in {\mathbb{P}}_{r-2}\left(e_i\right),i=1,2,\dots ,N,\hfill \end{array}$$

(15)

$$\begin{array}{cccc}\hfill & {\int }_{E_N}\psi qdx,\hfill & & \forall q\in {\mathbb{P}}_{r-N}\left(E_N\right),\hfill \end{array}$$

(16)

where $d\sigma$ is the one dimensional surface measure. Alternatively, one can use nodal DoFs (i.e., evaluation at a node point) in place of (15) and/or (16). For the former, on each edge $e_i$, its corresponding edge nodes are $r-1$ points such that they, along with the two vertices, are equally distributed on $e_i$. For the latter, the interior cell nodes can be set to be the Lagrange nodes of order $r-N$ of a triangle that lies strictly inside $E_N$.

The basis of ${\mathcal{DS}}_r\left(E_N\right)$ corresponding to the DoFs can be constructed. Given a computational mesh of convex polygons ${\mathcal{T}}_h$ over a domain $\Omega$, the basis can be simply pieced together to form a global $H^1$-conforming basis of the space ${\mathcal{DS}}_r(\Omega )\subset H^1(\Omega )$.

3.2. Direct Mixed Finite Elements

As discussed in the introduction, full, ${\mathbf{V}}_r^r\left(E_N\right)$, and reduced, ${\mathbf{V}}_r^{r-1}\left(E_N\right)$, $H\left(\mathrm{div}\right)$-approximating mixed finite element spaces follow from a de Rham complex (3), where the direct serendipity finite elements serve as the precursor (4). The supplemental space is related to ${\mathbb{S}}_{r+1}^{\mathcal{DS}}\left(E_N\right)$ by the simple Formula (6).

The DoFs for these spaces (with $s=r\ge 0$ or $s=r-1\ge 0$) can be taken to be

$$\begin{array}{cccc}\hfill & {\int }_{e_i}\mathbf{v}·{\nu }_{i}pd\sigma ,\hfill & & \forall p\in {\mathbb{P}}_r\left(e_i\right),i=1,2,\dots ,N,\hfill \end{array}$$

(17)

$$\begin{array}{cccc}\hfill & {\int }_{E_N}\mathbf{v}·\nabla qdx,\hfill & & \forall q\in {\mathbb{P}}_s\left(E_N\right),q\mathrm{not}\mathrm{constant},\hfill \end{array}$$

(18)

$$\begin{array}{cccc}\hfill & {\int }_{E_N}\mathbf{v}·\mathit{\psi }dx,\hfill & & \forall \mathit{\psi }\in {\mathbb{B}}_r^{\mathbf{V}}\left(E_N\right),\mathrm{if}r\ge N-1,\hfill \end{array}$$

(19)

where the $H^1\left(E_N\right)$ and $H(\mathrm{div};E_N)$ bubble functions, for $r\ge N-1$, are

$${\mathbb{B}}_{r+1}\left(E_N\right)={\lambda }_1{\lambda }_2\dots {\lambda }_N{\mathbb{P}}_{r-N+1}\left(E_N\right)\mathrm{and}{\mathbb{B}}_r^{\mathbf{V}}\left(E_N\right)=\mathrm{curl}{\mathbb{B}}_{r+1}\left(E_N\right).$$

(20)

Given the mesh ${\mathcal{T}}_h$ over $\Omega$, one constructs the basis and the $H\left(\mathrm{div}\right)$-conforming global space ${\mathbf{V}}_r^s(\Omega )\subset H\left(\mathrm{div}\right)$ (see [43] for details). As an alternative, when solving partial differential equations, one can use the hybrid form of the method [48], which does not require the construction of global basis functions.

4. Piecewise Continuous Supplements

In [43], $R_{i,j}$ satisfying (11) on $E_N$, for $1\le i<j\le N$, $2\le j-i\le N-2$, was taken to be the simple rational function

$$R_{i,j}^{\mathrm{rational}}\left(\mathbf{x}\right)=\frac{{\lambda }_i\left(\mathbf{x}\right)-{\lambda }_j\left(\mathbf{x}\right)}{{\lambda }_i\left(\mathbf{x}\right)+{\lambda }_j\left(\mathbf{x}\right)}.$$

(21)

These rational functions are smooth over the element. We now give new direct serendipity and mixed finite elements by providing an alternate construction of $R_{i,j}$ as a piecewise continuous polynomial defined over a sub-partition of $E_N$. We present two strategies, the first of which is convenient for the construction of continuous supplemental functions in $H^1\left(E_N\right)$, and the second for constructing smoother supplemental functions in $H^p\left(E_N\right)$ for integer $p\ge 1$.

4.1. Supplemental Functions in $H^1\left(E_N\right)$

Our first strategy for constructing $R_{i,j}$ requires a sub-triangulation of the element $E_N$, and we present two natural choices. The first sub-triangulation is depicted in Figure 3 and denoted as ${\mathcal{T}}^n\left(E_N\right)$. One picks a vertex ${\mathbf{x}}_{v,n}$ and divides $E_N$ into $N-2$ sub-triangles. The sub-triangles are $T_m^n$ with vertices ${\mathbf{x}}_{v,n}$, ${\mathbf{x}}_{v,m}$, and ${\mathbf{x}}_{v,m+1}$, where $m=n+1,\dots ,n+N-2$. For the second sub-triangulation, depicted in Figure 4 and denoted as ${\mathcal{T}}^{{\mathbf{x}}_c}\left(E_N\right)$, one picks a point ${\mathbf{x}}_c$ in the interior of $E_N$ and divides it into N sub-triangles. Now the sub-triangles are $T_m^{{\mathbf{x}}_c}$ with vertices ${\mathbf{x}}_c$, ${\mathbf{x}}_{v,m}$, and ${\mathbf{x}}_{v,m+1}$, where $m=1,2,\dots ,N$. We use the centroid of the element for ${\mathbf{x}}_c$.

Let the piecewise polynomial function space of degree s corresponding to each sub-triangulation be

$$\begin{array}{cccc}& {\mathcal{P}}_s\left({\mathcal{T}}^n\left(E_N\right)\right)\hfill & & =\{f\in C^0\left(E\right)|f{|}_{T_m^n}\in {\mathbb{P}}_s\left(T_m^n\right),m=n+1,\dots ,n+N-2\},\hfill \end{array}$$

(22)

$$\begin{array}{cccc}& {\mathcal{P}}_s\left({\mathcal{T}}^{{\mathbf{x}}_c}\left(E_N\right)\right)\hfill & & =\{f\in C^0\left(E\right)|f{|}_{T_m^{{\mathbf{x}}_c}}\in {\mathbb{P}}_s\left(T_m^{{\mathbf{x}}_c}\right),m=1,2,\dots ,N\}.\hfill \end{array}$$

(23)

We construct $R_{i,j}$ in ${\mathcal{P}}_1\left({\mathcal{T}}^n\left(E_N\right)\right)$ or ${\mathcal{P}}_1\left({\mathcal{T}}^{{\mathbf{x}}_c}\left(E_N\right)\right)$, depending on which of the two sub-triangulations is used, such that

$$\begin{array}{c}\hfill R_{i,j}{|}_{e_i}=-1,R_{i,j}{{|}_{e_j}=1,R_{i,j}|}_{v_k}=0,\forall k\ne i-1,i,j-1,j.\end{array}$$

(24)

by using interpolation at the vertices of the sub-triangles. If the sub-triangulation is chosen to be ${\mathcal{T}}^n\left(E_N\right)$, the restrictions (24) uniquely specify all the vertex values. However, if the triangulation is ${\mathcal{T}}^{{\mathbf{x}}_c}\left(E_N\right)$, the center value is not determined, so we assign $R_{i,j}\left({\mathbf{x}}_c\right)=0$.

Our construction has $R_{i,j}$ being $-1$ on $e_i$ and 1 on $e_j$ as required by (11). Moreover, $R_{i,j}\in H^1\left(E_N\right)$. After constructing the supplemental functions in (13) with this $R_{i,j}$, each ${\varphi }_{s,i,j}$ is in ${\mathcal{P}}_{r+1}\left({\mathcal{T}}^n\left(E_N\right)\right)$ or ${\mathcal{P}}_{r+1}\left({\mathcal{T}}^{{\mathbf{x}}_c}\left(E_N\right)\right)$, and therefore also in $H^1\left(E_N\right)$.

4.2. Supplemental Functions in $H^p\left(E_N\right)$

We now present the second of our two strategies for constructing $R_{i,j}$ for two nonadjacent edges $e_i$ and $e_j$. Recall that ${\lambda }_k\left(\mathbf{x}\right)$ is the linear polynomial giving the (signed) distance to the line ${\mathcal{L}}_k$ extending edge $e_k$. When $e_i$ and $e_j$ are parallel, we simply define $R_{i,j}$ as the linear polynomial

$$R_{i,j}=\frac{{\lambda }_i-{\lambda }_j}{{\lambda }_i\left({\mathbf{x}}_{v,j}\right)}.$$

(25)

When $e_i$ and $e_j$ are not parallel, we first define a sub-partition of $E_N$ by adding a single extra line ${\mathsf{\ell }}^{i,j}$ through a point ${\mathbf{x}}^{i,j}$ as depicted in Figure 5. The point ${\mathbf{x}}^{i,j}$ is chosen so that it is closer to ${\mathcal{L}}_j$ than the endpoints of $e_i$, i.e.,

$$\begin{array}{c}\hfill {\lambda }_j\left({\mathbf{x}}^{i,j}\right)\le min\{{\lambda }_j\left({\mathbf{x}}_{v,i-1}\right),{\lambda }_j\left({\mathbf{x}}_{v,i}\right)\}.\end{array}$$

(26)

The line ${\mathsf{\ell }}^{i,j}$ passes through ${\mathbf{x}}^{i,j}$ and is parallel to $e_j$. This line divides $E_N$ into $E_N^{i,j,1}$ near $e_i$ and $E_N^{i,j,0}$ near $e_j$, i.e.,

$$\begin{array}{cc}\hfill E_N^{i,j,1}& =E_N\cap \left\{\mathbf{x}\right|{\lambda }_j\left(\mathbf{x}\right)\ge {\lambda }_j\left({\mathbf{x}}^{i,j}\right)\},\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill E_N^{i,j,0}& =E_N\cap \left\{\mathbf{x}\right|{\lambda }_j\left(\mathbf{x}\right)<{\lambda }_j\left({\mathbf{x}}^{i,j}\right)\}.\hfill \end{array}$$

(28)

Let ${\nu }^{i,j}=-{\nu }_j$ be the unit normal vector of ${\mathsf{\ell }}^{i,j}$ pointing into $E_N^{i,j,1}$, and let ${\tau }^{i,j}={\tau }_j$ be a unit tangent vector.

We next construct the function ${\rho }^{i,j}$, which is 1 on edge $e_i$ and 0 on edge $e_j$. It is defined piecewise on the sub-partition of $E_N$ as

$$\begin{array}{c}\hfill {\rho }^{i,j}\left(\mathbf{x}\right)=\left\{\begin{array}{cc}1,\hfill & \mathbf{x}\in E_N^{i,j,1},\hfill \\ {\displaystyle 1-{\left(1-\frac{{\lambda }_j\left(\mathbf{x}\right)}{{\lambda }_j\left({\mathbf{x}}^{i,j}\right)}\right)}^p,}\hfill & \mathbf{x}\in E_N^{i,j,0},\hfill \end{array}\right.\end{array}$$

(29)

where $p\ge 1$ is an integer. The function is continuous, since ${\lambda }_j\left(\mathbf{x}\right)={\lambda }_j\left({\mathbf{x}}^{i,j}\right)$ on ${\mathsf{\ell }}^{i,j}$ implies that ${\rho }^{i,j}{|}_{{\mathsf{\ell }}^{i,j}}=1$ in either case of the definition. Moreover, in the tangential direction,

$$\begin{array}{c}\hfill \frac{\partial {\rho }^{i,j}}{\partial {\tau }^{i,j}}{|}_{{\mathsf{\ell }}^{i,j}}=0,\end{array}$$

(30)

and, in the normal direction,

$$\begin{array}{c}\hfill \frac{\partial {\rho }^{i,j}}{\partial {\nu }^{i,j}}\left(\mathbf{x}\right)=\left\{\begin{array}{cc}0,\hfill & \mathbf{x}\in E_N^{i,j,1},\hfill \\ {\displaystyle \frac{p}{{\lambda }_j\left({\mathbf{x}}^{i,j}\right)}{\left(1-\frac{{\lambda }_j\left(\mathbf{x}\right)}{{\lambda }_j\left({\mathbf{x}}^{i,j}\right)}\right)}^{p-1},}\hfill & \mathbf{x}\in E_N^{i,j,0},\hfill \end{array}\right.\end{array}$$

(31)

which is continuous for $p>1$, so ${\rho }^{i,j}\in C^1\left(E_N\right)$. By iterating the argument, we have that ${\rho }^{i,j}\in C^{p-1}\left(E_N\right)$ and so also in $H^p\left(E_N\right)$ for $p>1$. If $p=1$, ${\rho }^{i,j}$ is continuous, so it is in $H^1\left(E_N\right)$.

Finally, after constructing both ${\rho }^{i,j}$ and ${\rho }^{j,i}$, we define

$$R_{i,j}={\rho }^{j,i}-{\rho }^{i,j},$$

(32)

which is $-1$ on $e_i$, 1 on $e_j$. Moreover, $R_{i,j}\in H^p\left(E_N\right)$. The supplemental functions in (13) constructed with this $R_{i,j}$ lie in $H^p\left(E_N\right)$.

We end this section with two specific examples, using the sub-partitions shown in Figure 6, which divide $E_N$ by N lines. The first example has a sub-partition based on the midpoints ${\mathbf{x}}_{e,i}^{\mathrm{M}}$ of the edges $e_i$, $i=1,2,\dots ,N$, and gives rise to the spaces denoted ${\mathcal{DS}}_r^{\mathrm{M}}\left(E_N\right)$ and ${\mathbf{V}}_r^{\mathrm{M},s}\left(E_N\right)$. We compute the minimal distance of the midpoints to the edges, i.e.,

$$h_{\mathrm{M}}=\underset{1\le i\le N,k=i\pm 1}{min}{\lambda }_i\left({\mathbf{x}}_{e,k}^{\mathrm{M}}\right).$$

(33)

Then for any two non parallel and nonadjacent edges $e_i$ and $e_j$, simply take the partition line ${\mathsf{\ell }}^{i,j}$ to be the line parallel to $e_j$ that is the fixed distance $h_{\mathrm{M}}>0$ away and intersects $E_N$.

The second specific example uses a sub-partition based on trisecting each edge, resulting in the points, for edge $e_i$, being denoted counterclockwise as ${\mathbf{x}}_{e,i,k}$ for $k=1,2$, $i=1,2,\dots ,N$. In this case, we simply take ${\mathbf{x}}^{i,j}$ to be the closest of these points to ${\mathcal{L}}_j$, omitting ${\mathbf{x}}_{e,j,1}$ and ${\mathbf{x}}_{e,j,2}$. We denote the resulting spaces ${\mathcal{DS}}_r^{\mathrm{T}}\left(E_N\right)$ and ${\mathbf{V}}_r^{\mathrm{T},s}\left(E_N\right)$.

5. Approximation Properties

We discuss now the global approximation properties for our direct finite element spaces. The results of [43] do not directly apply here because there it was assumed that the functions $R_{i,j}$ are smooth on the element. Consider a collection of meshes ${\mathcal{T}}_h$ of convex polygons partitioning a domain $\Omega$, where $h>0$ is the maximal element diameter.

We need to make the usual assumption that our collection of meshes is uniformly shape regular [49] (pp. 104–105). For any $E_N\in {\mathcal{T}}_h$, let $h_{E_N}$ be its diameter. Denote by $T_i$, $i=1,2,\dots ,N(N-1)(N-2)/6$, the sub-triangle of $E_N$ with vertices being three of the N vertices of $E_N$, and define

$${\rho }_{E_N}=2\underset{1\le i\le N(N-1)(N-2)/6}{min}\left\{\mathrm{diameter}\mathrm{of}\mathrm{the}\mathrm{largest}\mathrm{circle}\mathrm{inscribed}\mathrm{in}{T}_i\right\}.$$

(34)

The shape regularity parameter of the single mesh ${\mathcal{T}}_h$ is

$${\sigma }_{{\mathcal{T}}_h}=\underset{E_N\in {\mathcal{T}}_h}{min}\frac{{\rho }_{E_N}}{h_{E_N}}.$$

(35)

Assumption 1. 

The collection of meshes ${\left\{{\mathcal{T}}_h\right\}}_{h>0}$ is  uniformly shape regular.  That is, the shape regularity parameters are bounded below by a positive constant: there exists ${\sigma }_{\ast }>0$, independent of ${\mathcal{T}}_h$ and $h>0$, such that the ratio

$$\frac{{\rho }_{E_N}}{h_{E_N}}\ge {\sigma }_{\ast }>0\mathit{for}\mathit{all}{E}_N\in {\mathcal{T}}_h,h>0.$$

(36)

We also require some mild restrictions on the construction of ${\mathbb{S}}_r^{\mathcal{DS}}\left(E_N\right)$.

Assumption 2. 

For every $E_N\in {\mathcal{T}}_h$, assume that the functions of ${\mathbb{S}}_r^{\mathcal{DS}}\left(E_N\right)$ are constructed using ${\lambda }_{i,j}$ such that the zero set ${\mathcal{L}}_{i,j}$ intersects $e_i$ and $e_j$. Moreover, suppose that $R_{i,j}\in H^p\left(E_N\right)$ for some $p\ge 1$ and that the sub-partitions introduced in Section 4 for their construction depend continuously on the vertices of $E_N$.

The continuous dependence requirement of the sub-partitions is met if we systematically choose the points ${\mathbf{x}}_c$ in Section 4.1 (say as the centroid) and ${\mathbf{x}}^{i,j}$ satisfying (26) in Section 4.2 (say be taking ${\mathbf{x}}^{i,j}$ as the closer endpoint of $e_i$ to ${\mathcal{L}}_j$, or so that ${\lambda }_j\left({\mathbf{x}}^{i,j}\right)=\frac{1}{2}min\{{\lambda }_j\left({\mathbf{x}}_{v,i-1}\right),{\lambda }_j\left({\mathbf{x}}_{v,i}\right)\}$).

We state first the approximation result for ${\mathcal{DS}}_r(\Omega )$.

Theorem 1. 

Let $r\ge 1$, $1\le p\le \infty$, and $\mathsf{\ell }>1/p$ (or $\mathsf{\ell }\ge 1$ if $p=1$). If Assumptions 1 and 2 hold (so the basis functions are in $H^p$ on each element), then there exists a constant $C>0$, such that for all functions $v\in W^{\mathsf{\ell },p}(\Omega )$,

$$\begin{array}{c}\hfill \underset{v_h\in {\mathcal{DS}}_r(\Omega )}{inf}\parallel v-v_h{\parallel }_{W^{m,p}(\Omega )}\le C{h}^{\mathsf{\ell }-m}{\parallel v\parallel }_{W^{\mathsf{\ell },p}(\Omega )},0\le \mathsf{\ell }\le r+1,m=0,1.\end{array}$$

(37)

Proof. 

The methodology of the proof follows [42,43]. The key difference is that we must relax the smoothness requirement made on the supplemental functions. We highlight the differences, and leave the reader to consult [42,43] for some of the details.

Given a mesh ${\mathcal{T}}_h$, we construct an interpolation operator ${\mathcal{I}}_h^r:W_p^l(\Omega )\to {\mathcal{DS}}_r$ as a generalization of that defined in [50]. To do so, we use a nodal set of DoFs for the finite elements, and identify global nodal points $a_i$, $i=1,2,\dots ,dim{\mathcal{DS}}_r$. These nodal points must be chosen systematically with respect to the vertices of the mesh, so they depend continuously on them. The global nodal basis function for $a_i$ is denoted ${\phi }_i$.

A geometry object $K_i$ is associated to each $a_i$. If $a_i$ lies in the interior of some element, we choose the element to be $K_i$. Otherwise, we choose an edge containing $a_i$ to be $K_i$, where we additionally ask that $K_i\subset \partial \Omega$ if $a_i\in \partial \Omega$. We use these to define the dual basis ${\psi }_i$ with respect to $L^2\left(K_i\right)$, $i=1,2,\dots ,dim{\mathcal{DS}}_r$. The corresponding interpolation operator ${\mathcal{I}}_h^r:W_p^l(\Omega )\to {\mathcal{DS}}_r$ is then

$${\mathcal{I}}_h^{r}v\left(\mathbf{x}\right)=\sum _{i=1}^{dim{\mathcal{DS}}_r}\left({\int }_{K_i}{\psi }_i\left(\mathbf{y}\right)v\left(\mathbf{y}\right)dy\right){\phi }_i\left(\mathbf{x}\right).$$

(38)

There are two essential steps towards showing the approximation property. First, the nodal basis functions are bounded,

$$\underset{1\le i\le dim{\mathcal{DS}}_r(\Omega )}{max}\underset{E\in {\mathcal{T}}_h}{max}\left|\right|{\phi }_i{\left|\right|}_{W_q^m\left(E\right)}\le C,$$

(39)

and, second, the dual basis functions are bounded up to a scaling factor,

$$\left|\right|{\psi }_i{\left|\right|}_{L^{\infty }\left(K_i\right)}\le C{h}_{K_i}^{-dim{K}_i}.$$

(40)

We show the necessary boundedness by mapping the elements and using a continuity and compactness argument.

As depicted in Figure 7, to each element $E_N$, we associate a regular polygon (equilateral and equiangular) ${\widehat{E}}_N$. We can then define a map ${\mathbf{F}}_{E_N}:{\widehat{E}}_N\to E_N$ as a composition of a map that changes the geometry but not the size to ${\tilde{E}}_N$, and then a scaling map (see [43] for precise details).

Define the nodal basis functions ${\phi }_i^{{\tilde{E}}_N}$ on ${\tilde{E}}_N$. It is enough to show the boundedness of their $W_q^m$ norms. Although they are no longer smooth functions, compared to [42,43], they are continuous on ${\tilde{E}}_N$, and smooth on all the subregions generated by the sub-partition. Moreover, by assumption the sub-partition is required to depend continuously on the vertices of ${\tilde{E}}_N$. Therefore, ${\phi }_i^{{\tilde{E}}_N}$ will still depend continuously on $\widehat{\mathbf{x}}={\mathbf{F}}_{{\tilde{E}}_N}^{-1}\left(\tilde{\mathbf{x}}\right)$ and the vertices of ${\tilde{E}}_N$, which vary in a compact set. We conclude that the nodal basis functions are bounded in $W_q^m$ norm. The boundedness of ${\psi }_i$ in the $L^{\infty }$ norm can be shown in a similar way. □

For the mixed finite elements, we have the following result, wherein we see projection operators $\pi :H(\mathrm{div};\Omega )\cap {\left(L^{2+ϵ}(\Omega )\right)}^2\to {\mathbf{V}}_r^s$, $s=r-1,r$, where $ϵ>0$, and ${\mathcal{P}}_{W_s}$, the $L^2$-orthogonal projection operator onto $W_s=\nabla ·{\mathbf{V}}_r^s$.

Theorem 2. 

Let $r=s\ge 0$ or $r\ge 1$, $s=r-1$. If Assumptions 1 and 2 hold, then there is a constant $C>0$, such that

$$\begin{array}{cccc}\hfill {\parallel \mathbf{v}-\pi \mathbf{v}\parallel }_{L^2(\Omega )}& \le {C\parallel \mathbf{v}\parallel }_{H^k(\Omega )}{h}^k,\hfill & & k=1,\dots ,r+1,\hfill \end{array}$$

(41)

$$\begin{array}{cccc}\hfill \parallel p-{\mathcal{P}}_{W_s}{p\parallel }_{L^2(\Omega )}& \le {C\parallel p\parallel }_{H^k(\Omega )}{h}^k,\hfill & & k=0,1,\dots ,s+1,\hfill \end{array}$$

(42)

$$\begin{array}{cccc}\hfill {\parallel \nabla ·(\mathbf{v}-\pi \mathbf{v})\parallel }_{L^2(\Omega )}& \le {C\parallel \nabla ·\mathbf{v}\parallel }_{H^k(\Omega )}{h}^k,\hfill & & k=0,1,\dots ,s+1,\hfill \end{array}$$

(43)

where $s=r-1\ge 0$ and $s=r\ge 1$ for reduced and full $H\left(\mathrm{div}\right)$-approximation, respectively. Moreover, the discrete inf-sup condition

$$\underset{{\mathbf{v}}_h\in {\mathbf{V}}_r^s}{sup}\frac{(w_h,\nabla ·{\mathbf{v}}_h)}{\parallel {\mathbf{v}}_h{\parallel }_{H\left(\mathrm{div}\right)}}\ge \gamma {\parallel w_h\parallel }_{L^2(\Omega )},\forall w_h\in W_s,$$

(44)

holds for some $\gamma >0$ independent of $h>0$.

For the proof, we define the projection operator $\pi$ by piecing together local operators ${\pi }_E$ that are defined in terms of the DoFs (17)–(19). The approximation properties given in [42,43] hold with a similar proof, using now that the subregions generated by the sub-partition depend continuously on the vertices of the element.

6. Numerical Results

We present numerical experiments for our new finite elements as applied to Poisson’s equation

$$\begin{array}{cccc}\hfill -\nabla ·(\nabla p)& =f\hfill & & \mathrm{in}\Omega ,\hfill \end{array}$$

(45)

$$\begin{array}{cccc}\hfill p& =0\hfill & & \mathrm{on}\partial \Omega ,\hfill \end{array}$$

(46)

where $f\in L^2(\Omega )$. The corresponding weak form finds $p\in H_0^1(\Omega )$ such that

$$\begin{array}{c}\hfill (\nabla p,\nabla q)=(f,q),\forall q\in H_0^1(\Omega ),\end{array}$$

(47)

where $(·,·)$ is the $L^2(\Omega )$ inner product. Setting

$$\mathbf{u}=-\nabla p,$$

(48)

we have the mixed weak form, which finds $\mathbf{u}\in H(\mathrm{div};\Omega )$ and $p\in L^2(\Omega )$ such that

$$\begin{array}{cccccc}\hfill & (\mathbf{u},\mathbf{v})-(p,\nabla ·\mathbf{v})\hfill & & =0,\hfill & & \forall \mathbf{v}\in H(\mathrm{div};\Omega ),\hfill \end{array}$$

(49)

$$\begin{array}{cccccc}\hfill & (\nabla ·\mathbf{u},w)\hfill & & =(f,w),\hfill & & \forall w\in L^2(\Omega ).\hfill \end{array}$$

(50)

These weak forms naturally give rise to finite element approximations. According to Theorems 1 and 2, the following convergence analysis holds by a standard argument [27,51].

Theorem 3. 

If Assumptions 1 and 2 hold, then there exists a constant $C>0$, independent of ${\mathcal{T}}_h$ and $h>0$, such that for $r\ge 1$,

$$\begin{array}{cc}\hfill \parallel p-p_h{\parallel }_{m,\Omega }& \le C{h}^{\mathsf{\ell }+1-m}{\left|p\right|}_{\mathsf{\ell }+1,\Omega },\mathsf{\ell }=0,1,\dots ,r,m=0,1,\hfill \end{array}$$

(51)

where $p_h\in {\mathcal{DS}}_r(\Omega )\cap H_0^1(\Omega )$ approximates (47). Moreover, with $r=s\ge 0$ or $r\ge 1$, $s=r-1$,

$$\begin{array}{cccc}\hfill \parallel \mathbf{u}-{\mathbf{u}}_h{\parallel }_{0,\Omega }& \le {C\parallel \mathbf{u}\parallel }_{k,\Omega }{h}^k,\hfill & & k=1,\dots ,r+1,\hfill \end{array}$$

(52)

$$\begin{array}{cccc}\hfill \parallel p-p_h{\parallel }_{0,\Omega }& \le {C\parallel \mathbf{u}\parallel }_{k,\Omega }{h}^k,\hfill & & k=1,\dots ,s+1,\hfill \end{array}$$

(53)

$$\begin{array}{cccc}\hfill \parallel \nabla ·(\mathbf{u}-{\mathbf{u}}_h){\parallel }_{0,\Omega }& \le {C\parallel \nabla ·\mathbf{u}\parallel }_{k,\Omega }{h}^k,\hfill & & k=0,1,\dots ,s+1,\hfill \end{array}$$

(54)

where $({\mathbf{u}}_h,p_h)\in {\mathbf{V}}_r^s \times W_s$ approximates (49) and (50).

We perform our tests on a unit square domain $\Omega ={[0,1]}^2$, and take the source term $f\left(\mathbf{x}\right)=2{\pi }^{2}sin\left(\pi x_1\right)sin\left(\pi x_2\right)$, so the exact solution is $u(x_1,x_2)=sin\left(\pi x_1\right)sin\left(\pi x_2\right)$. We consider five types of supplemental spaces. The original direct serendipity and mixed finite element spaces will be denoted ${\mathcal{DS}}_r^{\mathrm{R}}$ and ${\mathbf{V}}_r^{\mathrm{R},s}$, respectively. These use supplements based on the rational functions (21).

For the $H^1$ supplemental functions introduced in Section 4.1, there are two varieties. Denote the space using supplemental functions that are constructed based on the vertex sub-triangulation as ${\mathcal{DS}}_r^{\mathrm{V}}$ and its corresponding mixed spaces as ${\mathbf{V}}_r^{\mathrm{V},s}$, and those based on the center point sub-triangulation as ${\mathcal{DS}}_r^{\mathrm{C}}$ and ${\mathbf{V}}_r^{\mathrm{C},s}$, respectively. The spaces based on the $H^p$ supplements were described in Section 4.2 and denoted ${\mathcal{DS}}_r^{\mathrm{M}},{\mathbf{V}}_r^{\mathrm{M},s}$ and ${\mathcal{DS}}_r^{\mathrm{T}},{\mathbf{V}}_r^{\mathrm{T},s}$.

6.1. The Meshes Used

Approximate solutions are computed on a sequence of Voronoi meshes ${\mathcal{T}}_h^2$ generated by the package PolyMesher [52]. Each mesh has $n^2$ elements, which are generated with $n^2$ random initial seeds and up to ${10}^4$ smoothing iterations to improve the shape regularity.

For comparison to the results appearing in [43], we use the same mesh sequence ${\mathcal{T}}_h^2$ with $n=6$, 10, 14, 18, and 22. We show the meshes for $n=6$ and $n=18$ in Figure 8. The shape regularity parameters are given in

Table 1. Shape regularity parameters for each mesh ${\mathcal{T}}_h^2$.

	$\mathit{n}\mathbf{=}\mathbf{6}$	$\mathit{n}\mathbf{=}\mathbf{10}$	$\mathit{n}\mathbf{=}\mathbf{14}$	$\mathit{n}\mathbf{=}\mathbf{18}$	$\mathit{n}\mathbf{=}\mathbf{22}$
${\sigma }_{{\mathcal{T}}_h}$	0.180	0.115	0.161	0.127	0.150

. Note that the $n=10$ and $n=18$ meshes are the least regular.

In [43] it was observed numerically that the $n=18$ mesh performed well for the original direct finite elements (using rational supplemental functions) when $r=2,3,4$, but had a degraded convergence rate when $r=5$. The problem was resolved by removing short edges from the $n=18$ mesh. However, as we will see in this section, the problem is actually due to inaccurate numerical quadrature of the rational supplemental functions, which only showed up in those tests for the more refined mesh (i.e., not $n=10$) and higher values of r.

6.2. Results for Direct Serendipity Spaces

We present and compare the results of the numerical tests performed for ${\mathcal{DS}}_r^{\mathrm{R}}$, ${\mathcal{DS}}_r^{\mathrm{V}}$, ${\mathcal{DS}}_r^{\mathrm{C}}$, ${\mathcal{DS}}_r^{\mathrm{M}}$, and ${\mathcal{DS}}_r^{\mathrm{T}}$, where $r=2,3,4,5$. We take $p=1$ in Section 4.2 for the construction of ${\mathcal{DS}}_r^{\mathrm{M}}$ and ${\mathcal{DS}}_r^{\mathrm{T}}$, because it gives better results than a larger p in most cases. According to Theorem 3, we expect all those spaces to have the convergence rates $r+1$ for $L^2$ errors and r for $H^1$-seminorm errors. As

Table 2. $L^2$ errors and convergence rates for ${\mathcal{DS}}_r^{\mathrm{R}}$ and ${\mathcal{DS}}_r^{\mathrm{C}}$.

	${\mathcal{DS}}_2^{\mathrm{R}}$		${\mathcal{DS}}_2^{\mathrm{C}}$		${\mathcal{DS}}_3^{\mathrm{R}}$		${\mathcal{DS}}_3^{\mathrm{C}}$
n	error	rate	error	rate	error	rate	error	rate
10	2.160$\times {10}^{-4}$	3.45	2.144$\times {10}^{-4}$	3.50	8.859$\times {10}^{-6}$	4.34	1.031$\times {10}^{-5}$	4.35
14	7.329$\times {10}^{-5}$	3.16	7.165$\times {10}^{-5}$	3.21	2.175$\times {10}^{-6}$	4.11	2.518$\times {10}^{-6}$	4.13
18	3.452$\times {10}^{-5}$	2.95	3.409$\times {10}^{-5}$	2.92	7.927$\times {10}^{-7}$	3.96	8.964$\times {10}^{-7}$	4.05
22	1.863$\times {10}^{-5}$	3.47	1.841$\times {10}^{-5}$	3.46	3.555$\times {10}^{-7}$	4.51	4.045$\times {10}^{-7}$	4.48
	${\mathcal{DS}}_4^{\mathrm{R}}$		${\mathcal{DS}}_4^{\mathrm{C}}$		${\mathcal{DS}}_5^{\mathrm{R}}$		${\mathcal{DS}}_5^{\mathrm{C}}$
n	error	rate	error	rate	error	rate	error	rate
10	3.467$\times {10}^{-7}$	5.69	3.972$\times {10}^{-7}$	6.11	1.133$\times {10}^{-8}$	6.97	1.730$\times {10}^{-8}$	6.61
14	5.644$\times {10}^{-8}$	5.31	6.622$\times {10}^{-8}$	5.24	1.202$\times {10}^{-9}$	6.57	1.964$\times {10}^{-9}$	6.37
18	1.530$\times {10}^{-8}$	5.12	1.823$\times {10}^{-8}$	5.06	4.376$\times {10}^{-10}$	3.97	4.134$\times {10}^{-10}$	6.12
22	5.314$\times {10}^{-9}$	5.95	6.239$\times {10}^{-9}$	6.03	8.905$\times {10}^{-11}$	8.95	1.243$\times {10}^{-10}$	6.76

and

Table 3. $H^1$-seminorm errors and convergence rates for ${\mathcal{DS}}_r^{\mathrm{R}}$ and ${\mathcal{DS}}_r^{\mathrm{C}}$.

	${\mathcal{DS}}_2^{\mathrm{R}}$		${\mathcal{DS}}_2^{\mathrm{C}}$		${\mathcal{DS}}_3^{\mathrm{R}}$		${\mathcal{DS}}_3^{\mathrm{C}}$
n	error	rate	error	rate	error	rate	error	rate
10	3.561$\times {10}^{-3}$	2.32	3.552$\times {10}^{-3}$	2.36	1.933$\times {10}^{-4}$	3.13	2.390$\times {10}^{-4}$	3.11
14	1.683$\times {10}^{-3}$	2.19	1.660$\times {10}^{-3}$	2.23	6.724$\times {10}^{-5}$	3.09	8.343$\times {10}^{-5}$	3.08
18	1.018$\times {10}^{-3}$	1.97	1.013$\times {10}^{-3}$	1.94	3.144$\times {10}^{-5}$	2.98	3.783$\times {10}^{-5}$	3.10
22	6.712$\times {10}^{-4}$	2.34	6.696$\times {10}^{-4}$	2.33	1.730$\times {10}^{-5}$	3.36	2.114$\times {10}^{-5}$	3.27
	${\mathcal{DS}}_4^{\mathrm{R}}$		${\mathcal{DS}}_4^{\mathrm{C}}$		${\mathcal{DS}}_5^{\mathrm{R}}$		${\mathcal{DS}}_5^{\mathrm{C}}$
n	error	rate	error	rate	error	rate	error	rate
10	8.530$\times {10}^{-6}$	4.55	1.027$\times {10}^{-5}$	4.91	3.103$\times {10}^{-7}$	5.73	4.394$\times {10}^{-7}$	5.57
14	1.973$\times {10}^{-6}$	4.29	2.439$\times {10}^{-6}$	4.21	4.625$\times {10}^{-8}$	5.57	7.098$\times {10}^{-8}$	5.34
18	6.952$\times {10}^{-7}$	4.09	8.785$\times {10}^{-7}$	4.01	2.646$\times {10}^{-8}$	2.19	1.981$\times {10}^{-8}$	5.01
22	2.969$\times {10}^{-7}$	4.78	3.689$\times {10}^{-7}$	4.88	5.973$\times {10}^{-9}$	8.37	7.233$\times {10}^{-9}$	5.66

suggests, the convergence rates at $n=10,14,22$ for ${\mathcal{DS}}_r^{\mathrm{R}}$ are all slightly better than optimal in this test. However, we can observe a slower convergence rate at $n=18$ for ${\mathcal{DS}}_r^{\mathrm{R}}$. (We interject that convergence rates are computed as improvement from the previous mesh in our sequence.)

In

Table 4. Errors and convergence rates at $n=18$ computed from the previous step $n=14$, for the direct serendipity spaces ${\mathcal{DS}}_r^{\mathrm{R}}$, ${\mathcal{DS}}_r^{\mathrm{V}}$, ${\mathcal{DS}}_r^{\mathrm{C}}$, ${\mathcal{DS}}_r^{\mathrm{M}}$, and ${\mathcal{DS}}_r^{\mathrm{T}}$.

	$\mathit{r}\mathbf{=}\mathbf{2}$		$\mathit{r}\mathbf{=}\mathbf{3}$		$\mathit{r}\mathbf{=}\mathbf{4}$		$\mathit{r}\mathbf{=}\mathbf{5}$
	Error	Rate	Error	Rate	Error	Rate	Error	Rate
$L^2$ errors and convergence rates
${\mathcal{DS}}_r^{\mathrm{R}}$	3.452$\times {10}^{-5}$	2.95	7.927$\times {10}^{-7}$	3.96	1.530$\times {10}^{-8}$	5.12	4.376$\times {10}^{-10}$	3.97
${\mathcal{DS}}_r^{\mathrm{V}}$	3.554$\times {10}^{-5}$	2.89	1.073$\times {10}^{-6}$	3.87	2.108$\times {10}^{-8}$	4.83	4.637$\times {10}^{-10}$	5.92
${\mathcal{DS}}_r^{\mathrm{C}}$	3.409$\times {10}^{-5}$	2.92	8.964$\times {10}^{-7}$	4.05	1.823$\times {10}^{-8}$	5.06	4.13$\times {10}^{-10}$	6.12
${\mathcal{DS}}_r^{\mathrm{M}}$	6.820$\times {10}^{-5}$	2.88	1.697$\times {10}^{-6}$	3.85	3.095$\times {10}^{-8}$	4.80	6.11$\times {10}^{-10}$	6.02
${\mathcal{DS}}_r^{\mathrm{T}}$	7.072$\times {10}^{-5}$	2.92	1.866$\times {10}^{-6}$	4.04	3.367$\times {10}^{-8}$	4.91	5.83$\times {10}^{-10}$	6.07
$H^1$-seminorm errors and convergence rates
${\mathcal{DS}}_r^{\mathrm{R}}$	1.018$\times {10}^{-3}$	1.97	3.144$\times {10}^{-5}$	2.98	6.952$\times {10}^{-7}$	4.09	2.646$\times {10}^{-8}$	2.19
${\mathcal{DS}}_r^{\mathrm{V}}$	1.059$\times {10}^{-3}$	1.89	4.199$\times {10}^{-5}$	2.93	9.959$\times {10}^{-7}$	3.85	2.184$\times {10}^{-8}$	4.80
${\mathcal{DS}}_r^{\mathrm{C}}$	1.013$\times {10}^{-3}$	1.94	3.783$\times {10}^{-5}$	3.10	8.785$\times {10}^{-7}$	4.01	1.981$\times {10}^{-8}$	5.01
${\mathcal{DS}}_r^{\mathrm{M}}$	1.976$\times {10}^{-3}$	1.88	6.334$\times {10}^{-5}$	3.01	1.590$\times {10}^{-6}$	3.95	3.076$\times {10}^{-8}$	4.83
${\mathcal{DS}}_r^{\mathrm{T}}$	2.059$\times {10}^{-3}$	1.94	6.895$\times {10}^{-5}$	3.16	1.710$\times {10}^{-6}$	4.04	3.008$\times {10}^{-8}$	4.87

, we compare the results for the $n=18$ mesh of ${\mathcal{DS}}_r^{\mathrm{R}}$, ${\mathcal{DS}}_r^{\mathrm{V}}$, ${\mathcal{DS}}_r^{\mathrm{C}}$, ${\mathcal{DS}}_r^{\mathrm{M}}$, and ${\mathcal{DS}}_r^{\mathrm{T}}$. On the one hand, the results suggest that the new spaces are all approximately optimal for $r=5$, which is an obvious improvement compared to ${\mathcal{DS}}_5^{\mathrm{R}}$. On the other hand, the errors for $r=2,3,4$ of the new spaces are slightly worse than those of ${\mathcal{DS}}_r^{\mathrm{R}}$, and among all the new spaces, ${\mathcal{DS}}_r^{\mathrm{C}}$ has the best performance in error. We conclude that ${\mathcal{DS}}_r^{\mathrm{C}}$ shows the best overall performance among all the spaces considered.

We suggest that the reason for such an observation is that the dominant errors for $r=5$ are from the numerical quadrature applied to the integration of rational functions, especially on the elements that are less shape regular. However, for $r=2,3,4$, the new supplements, as piecewise polynomials, cannot approximate the shape of a smooth function as well as the original rational supplements, especially those from ${\mathcal{DS}}_r^{\mathrm{M}}$ and ${\mathcal{DS}}_r^{\mathrm{T}}$, of which $R_{i,j}$ for $e_i\nparallel e_j$ are flat in the middle and oscillate near the boundary. In contrast, the supplements from ${\mathcal{DS}}_r^{\mathrm{V}}$ and ${\mathcal{DS}}_r^{\mathrm{C}}$ are more reasonably shaped, and those from ${\mathcal{DS}}_r^{\mathrm{C}}$ are better since its partition has sub-triangles that are more shape regular (as was shown in Figure 3 and Figure 4). This argument is also supported by the observation that the results are usually worse if we take larger p for ${\mathcal{DS}}_r^{\mathrm{M}}$ and ${\mathcal{DS}}_r^{\mathrm{T}}$, where the shape of the supplements are even worse.

6.3. Results for Direct Mixed Spaces

We perform numerical tests for ${\mathbf{V}}_r^{\mathrm{R},s}$, ${\mathbf{V}}_r^{\mathrm{V},s}$, ${\mathbf{V}}_r^{\mathrm{C},s}$, ${\mathbf{V}}_r^{\mathrm{M},s}$, ${\mathbf{V}}_r^{\mathrm{T},s}$, for the full $H\left(\mathrm{div}\right)$-approximation spaces where $r=s=0,1,2,3,$ and the reduced $H\left(\mathrm{div}\right)$-approximation spaces where $r=1,2,3$, and $s=r-1$. Since those mixed spaces are constructed from corresponding direct serendipity spaces ${\mathcal{DS}}_{r+1}$, it is natural that we find the comparison of the results similar to the small r cases discussed in Section 6.2. For all the spaces, we can observe the convergence rates approximately optimal in general but the errors are slightly worse for $n=18$, especially when $r=s=3$, as shown in

Table 5. Errors and convergence rates at $n=18$ computed from the previous step $n=14$, for the reduced $H\left(\mathrm{div}\right)$-approximation spaces ${\mathbf{V}}_r^{\mathrm{R},r-1}$, ${\mathbf{V}}_r^{\mathrm{V},r-1}$, ${\mathbf{V}}_r^{\mathrm{C},r-1}$, ${\mathbf{V}}_r^{\mathrm{M},r-1}$, and ${\mathbf{V}}_r^{\mathrm{T},r-1}$.

	$\mathbf{\parallel }\mathit{p}\mathbf{-}{\mathit{p}}_{\mathit{h}}\mathbf{\parallel }$		$\mathbf{\parallel }\mathit{u}\mathbf{-}{\mathit{u}}_{\mathit{h}}\mathbf{\parallel }$		$\mathbf{\parallel }\mathit{\nabla }\mathbf{·}\mathbf{(}\mathit{u}\mathbf{-}{\mathit{u}}_{\mathit{h}}\mathbf{)}\mathbf{\parallel }$
	Error	Rate	Error	Rate	Error	Rate
$r=1$, reduced $H\left(\mathrm{div}\right)$-approximation
${\mathbf{V}}_1^{\mathrm{R},0}$	7.039$\times {10}^{-2}$	1.01	5.428$\times {10}^{-3}$	1.98	6.988$\times {10}^{-2}$	0.99
${\mathbf{V}}_1^{\mathrm{V},0}$	7.039$\times {10}^{-2}$	1.01	5.429$\times {10}^{-3}$	1.98	6.988$\times {10}^{-2}$	0.99
${\mathbf{V}}_1^{\mathrm{C},0}$	7.039$\times {10}^{-2}$	1.01	5.443$\times {10}^{-3}$	1.98	6.988$\times {10}^{-2}$	0.99
${\mathbf{V}}_1^{\mathrm{M},0}$	7.039$\times {10}^{-2}$	1.01	5.366$\times {10}^{-3}$	1.98	6.988$\times {10}^{-2}$	0.99
${\mathbf{V}}_1^{\mathrm{T},0}$	7.039$\times {10}^{-2}$	1.01	5.362$\times {10}^{-3}$	1.98	6.988$\times {10}^{-2}$	0.99
$r=2$, reduced $H\left(\mathrm{div}\right)$-approximation
${\mathbf{V}}_2^{\mathrm{R},1}$	2.614$\times {10}^{-3}$	1.96	8.492$\times {10}^{-5}$	2.92	2.614$\times {10}^{-3}$	1.96
${\mathbf{V}}_2^{\mathrm{V},1}$	2.614$\times {10}^{-3}$	1.96	8.876$\times {10}^{-5}$	2.89	2.614$\times {10}^{-3}$	1.96
${\mathbf{V}}_2^{\mathrm{C},1}$	2.614$\times {10}^{-3}$	1.96	8.850$\times {10}^{-5}$	2.87	2.614$\times {10}^{-3}$	1.96
${\mathbf{V}}_2^{\mathrm{M},1}$	2.614$\times {10}^{-3}$	1.96	8.895$\times {10}^{-5}$	2.85	2.614$\times {10}^{-3}$	1.96
${\mathbf{V}}_2^{\mathrm{T},1}$	2.614$\times {10}^{-3}$	1.96	8.973$\times {10}^{-5}$	2.85	2.614$\times {10}^{-3}$	1.96
$r=3$, reduced $H\left(\mathrm{div}\right)$-approximation
${\mathbf{V}}_3^{\mathrm{R},2}$	6.515$\times {10}^{-5}$	2.96	1.887$\times {10}^{-6}$	3.90	6.515$\times {10}^{-5}$	2.96
${\mathbf{V}}_3^{\mathrm{V},2}$	6.515$\times {10}^{-5}$	2.96	1.931$\times {10}^{-6}$	3.89	6.515$\times {10}^{-5}$	2.96
${\mathbf{V}}_3^{\mathrm{C},2}$	6.515$\times {10}^{-5}$	2.96	1.911$\times {10}^{-6}$	3.89	6.515$\times {10}^{-5}$	2.96
${\mathbf{V}}_3^{\mathrm{M},2}$	6.515$\times {10}^{-5}$	2.96	2.007$\times {10}^{-6}$	3.81	6.515$\times {10}^{-5}$	2.96
${\mathbf{V}}_3^{\mathrm{T},2}$	6.515$\times {10}^{-5}$	2.96	2.105$\times {10}^{-6}$	3.83	6.515$\times {10}^{-5}$	2.96

and

Table 6. Errors and convergence rates at $n=18$ computed from the previous step $n=14$, for the full $H\left(\mathrm{div}\right)$-approximation spaces ${\mathbf{V}}_r^{\mathrm{R},r}$, ${\mathbf{V}}_r^{\mathrm{V},r}$, ${\mathbf{V}}_r^{\mathrm{C},r}$, ${\mathbf{V}}_r^{\mathrm{M},r}$, and ${\mathbf{V}}_r^{\mathrm{T},r}$.

	$\mathbf{\parallel }\mathit{p}\mathbf{-}{\mathit{p}}_{\mathit{h}}\mathbf{\parallel }$		$\mathbf{\parallel }\mathit{u}\mathbf{-}{\mathit{u}}_{\mathit{h}}\mathbf{\parallel }$		$\mathbf{\parallel }\mathit{\nabla }\mathbf{·}\mathbf{(}\mathit{u}\mathbf{-}{\mathit{u}}_{\mathit{h}}\mathbf{)}\mathbf{\parallel }$
$\mathit{n}$	Error	Rate	Error	Rate	Error	Rate
$r=0$, full $H\left(\mathrm{div}\right)$-approximation
${\mathbf{V}}_0^{\mathrm{R},0}$	7.030$\times {10}^{-2}$	1.01	2.701$\times {10}^{-2}$	1.10	6.988$\times {10}^{-2}$	0.99
${\mathbf{V}}_0^{\mathrm{V},0}$	7.027$\times {10}^{-2}$	1.01	3.095$\times {10}^{-2}$	1.03	6.988$\times {10}^{-2}$	0.99
${\mathbf{V}}_0^{\mathrm{C},0}$	7.028$\times {10}^{-2}$	1.01	2.951$\times {10}^{-2}$	1.03	6.988$\times {10}^{-2}$	0.99
${\mathbf{V}}_0^{\mathrm{M},0}$	7.027$\times {10}^{-2}$	1.01	3.065$\times {10}^{-2}$	0.93	6.988$\times {10}^{-2}$	0.99
${\mathbf{V}}_0^{\mathrm{T},0}$	7.026$\times {10}^{-2}$	1.01	3.163$\times {10}^{-2}$	0.92	6.988$\times {10}^{-2}$	0.99
$r=1$, full $H\left(\mathrm{div}\right)$-approximation
${\mathbf{V}}_1^{\mathrm{R},1}$	2.614$\times {10}^{-3}$	1.96	4.895$\times {10}^{-4}$	2.19	2.614$\times {10}^{-3}$	1.96
${\mathbf{V}}_1^{\mathrm{V},1}$	2.614$\times {10}^{-3}$	1.96	5.542$\times {10}^{-4}$	2.13	2.614$\times {10}^{-3}$	1.96
${\mathbf{V}}_1^{\mathrm{C},1}$	2.614$\times {10}^{-3}$	1.96	5.226$\times {10}^{-4}$	2.17	2.614$\times {10}^{-3}$	1.96
${\mathbf{V}}_1^{\mathrm{M},1}$	2.614$\times {10}^{-3}$	1.96	7.505$\times {10}^{-4}$	2.08	2.614$\times {10}^{-3}$	1.96
${\mathbf{V}}_1^{\mathrm{T},1}$	2.614$\times {10}^{-3}$	1.96	7.917$\times {10}^{-4}$	2.15	2.614$\times {10}^{-3}$	1.96
$r=2$, full $H\left(\mathrm{div}\right)$-approximation
${\mathbf{V}}_2^{\mathrm{R},2}$	6.515$\times {10}^{-5}$	2.96	8.818$\times {10}^{-6}$	3.10	6.515$\times {10}^{-5}$	2.96
${\mathbf{V}}_2^{\mathrm{V},2}$	6.515$\times {10}^{-5}$	2.96	1.887$\times {10}^{-5}$	2.92	6.515$\times {10}^{-5}$	2.96
${\mathbf{V}}_2^{\mathrm{C},2}$	6.515$\times {10}^{-5}$	2.96	1.526$\times {10}^{-5}$	3.03	6.515$\times {10}^{-5}$	2.96
${\mathbf{V}}_2^{\mathrm{M},2}$	6.515$\times {10}^{-5}$	2.96	2.801$\times {10}^{-5}$	2.49	6.515$\times {10}^{-5}$	2.96
${\mathbf{V}}_2^{\mathrm{T},2}$	6.515$\times {10}^{-5}$	2.96	3.010$\times {10}^{-5}$	2.67	6.515$\times {10}^{-5}$	2.96
$r=3$, full $H\left(\mathrm{div}\right)$-approximation
${\mathbf{V}}_3^{\mathrm{R},3}$	1.182$\times {10}^{-6}$	3.99	2.144$\times {10}^{-7}$	3.65	1.182$\times {10}^{-6}$	3.99
${\mathbf{V}}_3^{\mathrm{V},3}$	1.182$\times {10}^{-6}$	3.99	3.324$\times {10}^{-7}$	3.43	1.182$\times {10}^{-6}$	3.99
${\mathbf{V}}_3^{\mathrm{C},3}$	1.182$\times {10}^{-6}$	3.99	2.933$\times {10}^{-7}$	3.50	1.182$\times {10}^{-6}$	3.99
${\mathbf{V}}_3^{\mathrm{M},3}$	1.183$\times {10}^{-6}$	3.99	1.254$\times {10}^{-6}$	3.10	1.182$\times {10}^{-6}$	3.99
${\mathbf{V}}_3^{\mathrm{T},3}$	1.183$\times {10}^{-6}$	3.99	1.547$\times {10}^{-6}$	3.61	1.182$\times {10}^{-6}$	3.99

. All spaces perform similarly well, although ${\mathbf{V}}_r^{\mathrm{R},s}$ usually performs best in these tests. Among the new spaces, ${\mathbf{V}}_r^{\mathrm{C},s}$ performs a bit better, and it gives results close to those of ${\mathbf{V}}_r^{\mathrm{R},s}$. For reference, we provide the numerical results for ${\mathbf{V}}_r^{\mathrm{C},s}$ in

Table 7. Errors and convergence rates in $L^2$ for ${\mathbf{V}}_r^{\mathrm{C},r-1}$.

	$\mathbf{\parallel }\mathit{p}\mathbf{-}{\mathit{p}}_{\mathit{h}}\mathbf{\parallel }$		$\mathbf{\parallel }\mathit{u}\mathbf{-}{\mathit{u}}_{\mathit{h}}\mathbf{\parallel }$		$\mathbf{\parallel }\mathit{\nabla }\mathbf{·}\mathbf{(}\mathit{u}\mathbf{-}{\mathit{u}}_{\mathit{h}}\mathbf{)}\mathbf{\parallel }$
$\mathit{n}$	Error	Rate	Error	Rate	Error	Rate
$r=1$, reduced $H\left(\mathrm{div}\right)$-approximation
10	1.290$\times {10}^{-1}$	1.24	1.775$\times {10}^{-2}$	2.29	1.260$\times {10}^{-1}$	1.15
14	9.109$\times {10}^{-2}$	1.02	9.024$\times {10}^{-3}$	1.98	9.001$\times {10}^{-2}$	0.98
18	7.039$\times {10}^{-2}$	1.01	5.443$\times {10}^{-3}$	1.98	6.988$\times {10}^{-2}$	0.99
22	5.736$\times {10}^{-2}$	1.15	3.630$\times {10}^{-3}$	2.28	5.708$\times {10}^{-2}$	1.14
$r=2$, reduced $H\left(\mathrm{div}\right)$-approximation
10	8.635$\times {10}^{-3}$	2.23	5.210$\times {10}^{-4}$	3.28	8.634$\times {10}^{-3}$	2.23
14	4.308$\times {10}^{-3}$	2.04	1.841$\times {10}^{-4}$	3.04	4.308$\times {10}^{-3}$	2.03
18	2.614$\times {10}^{-3}$	1.96	8.850$\times {10}^{-5}$	2.87	2.614$\times {10}^{-3}$	1.96
22	1.715$\times {10}^{-3}$	2.37	4.772$\times {10}^{-5}$	3.47	1.715$\times {10}^{-3}$	2.37
$r=3$, reduced $H\left(\mathrm{div}\right)$-approximation
10	3.881$\times {10}^{-4}$	3.38	2.021$\times {10}^{-5}$	4.39	3.881$\times {10}^{-4}$	3.38
14	1.384$\times {10}^{-4}$	3.02	5.151$\times {10}^{-6}$	4.00	1.384$\times {10}^{-4}$	3.02
18	6.515$\times {10}^{-5}$	2.96	1.911$\times {10}^{-6}$	3.89	6.515$\times {10}^{-5}$	2.96
22	3.507$\times {10}^{-5}$	3.48	8.432$\times {10}^{-7}$	4.60	3.507$\times {10}^{-5}$	3.48

and

Table 8. Errors and convergence rates in $L^2$ for ${\mathbf{V}}_r^{\mathrm{C},r}$.

	$\mathbf{\parallel }\mathit{p}\mathbf{-}{\mathit{p}}_{\mathit{h}}\mathbf{\parallel }$		$\mathbf{\parallel }\mathit{u}\mathbf{-}{\mathit{u}}_{\mathit{h}}\mathbf{\parallel }$		$\mathbf{\parallel }\mathit{\nabla }\mathbf{·}\mathbf{(}\mathit{u}\mathbf{-}{\mathit{u}}_{\mathit{h}}\mathbf{)}\mathbf{\parallel }$
$\mathit{n}$	Error	Rate	Error	Rate	Error	Rate
$r=0$, full $H\left(\mathrm{div}\right)$-approximation
10	1.281$\times {10}^{-1}$	1.20	6.389$\times {10}^{-2}$	1.54	1.260$\times {10}^{-1}$	1.15
14	9.086$\times {10}^{-2}$	1.01	3.832$\times {10}^{-2}$	1.50	9.001$\times {10}^{-2}$	0.98
18	7.028$\times {10}^{-2}$	1.01	2.951$\times {10}^{-2}$	1.03	6.988$\times {10}^{-2}$	0.99
22	5.731$\times {10}^{-2}$	1.15	2.145$\times {10}^{-2}$	1.79	5.708$\times {10}^{-2}$	1.14
$r=1$, full $H\left(\mathrm{div}\right)$-approximation
10	8.635$\times {10}^{-3}$	2.23	2.003$\times {10}^{-3}$	2.65	8.634$\times {10}^{-3}$	2.23
14	4.308$\times {10}^{-3}$	2.04	9.076$\times {10}^{-4}$	2.32	4.308$\times {10}^{-3}$	2.03
18	2.614$\times {10}^{-3}$	1.96	5.226$\times {10}^{-4}$	2.17	2.614$\times {10}^{-3}$	1.96
22	1.715$\times {10}^{-3}$	2.37	3.320$\times {10}^{-4}$	2.55	1.715$\times {10}^{-3}$	2.37
$r=2$, full $H\left(\mathrm{div}\right)$-approximation
10	3.881$\times {10}^{-4}$	3.38	1.007$\times {10}^{-4}$	3.47	3.881$\times {10}^{-4}$	3.38
14	1.384$\times {10}^{-4}$	3.02	3.303$\times {10}^{-5}$	3.26	1.384$\times {10}^{-4}$	3.02
18	6.515$\times {10}^{-5}$	2.96	1.526$\times {10}^{-5}$	3.03	6.515$\times {10}^{-5}$	2.96
22	3.507$\times {10}^{-5}$	3.48	7.889$\times {10}^{-6}$	3.71	3.507$\times {10}^{-5}$	3.48
$r=3$, full $H\left(\mathrm{div}\right)$-approximation
10	1.294$\times {10}^{-5}$	4.60	3.537$\times {10}^{-6}$	4.84	1.294$\times {10}^{-5}$	4.60
14	3.270$\times {10}^{-6}$	4.03	7.157$\times {10}^{-7}$	4.68	3.270$\times {10}^{-6}$	4.03
18	1.182$\times {10}^{-6}$	3.99	2.933$\times {10}^{-7}$	3.50	1.182$\times {10}^{-6}$	3.99
22	5.219$\times {10}^{-7}$	4.60	1.301$\times {10}^{-7}$	4.57	5.219$\times {10}^{-7}$	4.60

.

7. Conclusions

We reviewed the construction of direct serendipity and mixed finite elements on non-degenerate, planar convex polygons. The direct serendipity finite element spaces are of the form

$$\begin{array}{c}\hfill {\mathcal{DS}}_r\left(E_N\right)={\mathbb{P}}_r\left(E_N\right)\oplus {\mathbb{S}}_r^{\mathcal{DS}}\left(E_N\right).\end{array}$$

(55)

The full and reduced $H\left(\mathrm{div}\right)$-approximation mixed finite element spaces are obtained from a de Rham complex, where the direct serendipity finite elements serve as a precursor. The mixed spaces are of the form

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{V}}_r^r\left(E_N\right)={\mathbb{P}}_r^2\left(E_N\right)\oplus \mathbf{x}{\tilde{\mathbb{P}}}_r\left(E_N\right)\oplus {\mathbb{S}}_r^{\mathbf{V}}\left(E_N\right),\\ {\mathbf{V}}_r^{r-1}\left(E_N\right)={\mathbb{P}}_r^2\left(E_N\right)\oplus {\mathbb{S}}_r^{\mathbf{V}}\left(E_N\right),\end{array}\end{array}$$

(56)

where

$${\mathbb{S}}_r^{\mathbf{V}}\left(E_N\right)=\mathrm{curl}{\mathbb{S}}_{r+1}^{\mathcal{DS}}\left(E_N\right).$$

(57)

We presented two approaches to construct the supplemental functions in ${\mathbb{S}}_r^{\mathcal{DS}}\left(E_N\right)$ as piecewise polynomials. The first approach divides a polygonal element $E_N$ into sub-triangles, and constructs the supplements as continuous piecewise polynomials that lie in $H^1\left(E_N\right)$. The second approach has a more complicated subdivision of $E_N$ that needs to be treated carefully. However, it provides a framework for constructing supplements that lie in $H^p\left(E_N\right)$ for any $p\ge 1$.

The approximation properties of the new finite elements were proved under the regularity assumption of the mesh sequences and some mild restrictions on the construction.

We performed numerical tests on a randomly generated mesh sequence and compared results for five different ways of constructing the supplemental functions, including the original construction using smooth but rational functions. The comparison suggested that it is better to use the piecewise polynomial supplements rather than the rational supplements for higher order r. Although the rational supplements are smooth and so tend to approximate better, noticeable errors could be seen due to inaccurate numerical integration, especially on meshes with short edges. Among the new spaces, it was found that the spaces with supplements based on the center point sub-triangulation (23) performed best.