Direct serendipity and mixed finite elements on convex quadrilaterals

The classical serendipity and mixed finite element spaces suffer from poor approximation on nondegenerate, convex quadrilaterals. In this paper, we develop families of direct serendipity and direct mixed finite element spaces, which achieve optimal approximation properties and have minimal local dimension. The set of local shape functions for either the serendipity or mixed elements contains the full set of scalar or vector polynomials of degree r, respectively, defined directly on each element (i.e., not mapped from a reference element). Because there are not enough degrees of freedom for global H1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^1$$\end{document} or H(div)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H(\text {div})$$\end{document} conformity, exactly two supplemental shape functions must be added to each element when r≥2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r\ge 2$$\end{document}, and only one when r=1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=1$$\end{document}. The specific choice of supplemental functions gives rise to different families of direct elements. These new spaces are related through a de Rham complex. For index r≥1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r\ge 1$$\end{document}, the new families of serendipity spaces DSr+1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr {DS}}_{r+1}$$\end{document} are the precursors under the curl operator of our direct mixed finite element spaces, which can be constructed to have reduced or full H(div)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H(\text {div})$$\end{document} approximation properties. One choice of direct serendipity supplements gives the precursor of the recently introduced Arbogast–Correa spaces (SIAM J Numer Anal 54:3332–3356, 2016. https://doi.org/10.1137/15M1013705). Other fully direct serendipity supplements can be defined without the use of mappings from reference elements, and these give rise in turn to fully direct mixed spaces. Our development is constructive, so we are able to give global bases for our spaces. Numerical results are presented to illustrate their properties.

constructive, so we are able to give global bases for our spaces. Numerical results are presented to illustrate their properties.

Introduction
Serendipity finite elements S r (Ê) can be defined on a rectangleÊ [19,23,46]. Over a rectangular mesh, they merge together to form H 1 conforming spaces of scalar functions. These finite elements have the distinction of having a minimal number of degrees of freedom (DoFs) for the given order of approximation r +1 in L 2 . Similarly, the Brezzi-Douglas-Marini mixed finite elements BDM r (Ê) [20] are defined so that they merge together on a rectangular mesh into H (div) = v ∈ (L 2 ) 2 : ∇ · v ∈ L 2 conforming spaces of vector functions. These finite elements also have a minimal number of DoFs for the same order of approximation.
Both of these elements appear in the periodic table of the finite elements as given by Arnold and Logg [13] (where they are denoted S r Λ 0 and S r Λ 1 , respectively). They should be studied together, since they are related by a de Rham complex [5,11,12] which implies that BDM r (Ê) = curl S r +1 (Ê) ⊕ xP r −1 (Ê), where P r −1 (Ê) are polynomials of degree r − 1.
In this paper, we define new families of (we call them direct) serendipity and mixed finite elements on a general nondegenerate, convex quadrilateral E. These new finite elements generalize the complex (1), and they maintain H 1 or H (div) conformity, provide optimal order approximation properties, and possess a minimal number of DoFs.

Existing finite elements
The serendipity finite elements on rectangles S r (Ê), especially the 8-node biquadratic (r = 2) and the 12-node bicubic (r = 3) elements, have been well studied for many years. They appear in almost any introductory reference on finite elements, e.g., [19,23,46], and they are provided by software packages both in academia [27] and industry [35]. Compared with the full tensor product Lagrange finite elements P r ,r (Ê), serendipity finite elements use fewer degrees of freedom, and they are usually more efficient in terms of the number of local computations performed. It was not until recently, however, that a general definition of the serendipity finite element spaces of arbitrary order on rectangles in any space dimension was given by Arnold and Awanou [6,7] (see also [32]).
The serendipity finite element spaces work very well on computational meshes of rectangular elements, but it is well known that their performance is degraded on quadrilaterals when the space is mapped from a rectangle, when r ≥ 2. This is not the case for tensor product Lagrange finite elements [8,36,38]. To be more precise, mapped serendipity elements of index r do not approximate to optimal order r + 1 on E, but the image of the full space of tensor product polynomials P r ,r (Ê) maintains accuracy on E. We note that Rand, Gillette, and Bajaj [40] recently introduced a new family of Serendipity finite elements based on generalized barycentric coordinates of index r = 2 that is accurate to order three on any convex, planar polygon. A generalization to any order of approximation was given by Floater and Lai [31], but on quadrilaterals, they require dim P r + r shape functions, which is more than the minimal required when r > 2.
There are many families of mixed finite elements on rectangles, beginning with those of Raviart and Thomas [41] and generalized by Nédélec [39]. These and the BDM r finite elements are extended to quadrilaterals using the Piola transform [41,48]. For most spaces, this creates a consistency error and consequent loss of approximation of the divergence [1,9,16,21,48].
The construction of mixed finite elements on quadrilaterals that maintain optimal order accuracy is considered in many papers. Most address only low order cases (see, e.g., [15,17,22,26,37,43,44]). The exceptions we are aware of are the families of finite elements of Arnold, Boffi, and Falk (ABF r (E)) [9], Siqueira, Devloo, and Gomes [45], and Arbogast and Correa (AC r (E) and AC red r (E), written in this paper as AC r r (E) and AC r −1 r (E), respectively) [1]. The ABF elements are defined for rectangles and extended to quadrilaterals in the usual way (i.e., by mapping via the Piola transformation). They rectify the problem of poor divergence approximation by including more degrees of freedom in the space, so that approximation properties are maintained after Piola mapping. The spaces of [45] also involve the Piola map, but in a unique way. They also add shape functions to their space to obtain accuracy. The AC elements use a different strategy. These elements are defined by using vector polynomials directly on the element (i.e., without being mapped) and supplemented by two vector shape functions defined on a reference square and mapped via Piola. The AC spaces have minimal local dimension with respect to the requirement of global H (div) conformity and optimal order of approximation.

New finite elements
In this paper, we introduce new families of direct serendipity and mixed finite elements that have optimal approximation properties of all orders r = (0, )1, 2, . . . and maintain minimal local dimension. They are direct in the sense that the shape functions contain a full set of polynomials of degree r defined directly on the element, as in the AC spaces. Because there are not enough degrees of freedom to achieve H 1 or H (div) conformity over meshes of quadrilaterals, supplemental functions need to be added to each element. These supplemental functions can be defined in many ways as we will see, and each choice gives rise to a different family of direct finite elements. These families are novel, and each finite element of a given index r in the family is a novel finite element, except possibly for the low order elements (r ≤ 2) of some specific families, since as we noted above, many low order finite elements are known to exist.
The families of direct serendipity elements have the same number of degrees of freedom as the corresponding classical serendipity element, and they take the form DS r (E) = P r (E) ⊕ S DS r (E), r ≥ 1. (2) Each family is defined by the choice of the two supplemental functions (or one, if r = 1) spanning S DS r (E). We give a very general, but explicit, construction for these supplements. They can be defined directly on E, or they can be defined onÊ and mapped to E. In fact, we will construct a nodal basis over E, and these nodal basis functions can be merged together to give an H 1 conforming global basis. When r = 1 we obtain new elements akin to barycentric coordinates (i.e., they are linear on edges and sum to one, but they are not necessarily positive everywhere).
There are two classes of families of direct mixed elements, which correspond to reduced and full H (div)-approximation. For index r , a vector function is approximated to order r +1 accuracy, but the divergence of the vector is approximated to order r −1 or r for reduced and full H (div)-approximation, respectively. Each class of direct mixed elements has the same optimal number of degrees of freedom as the AC elements of that class. They take a form similar to (2), which is for the reduced and full H (div)-approximation spaces, respectively, whereP r are homogeneous polynomials of degree r . Again, each family is defined by the choice of the two supplemental functions spanning S V r (E). When r = 0, we also construct V 0 0 (E) = P 2 0 (E) ⊕ xP 0 (E) ⊕ S V 0 (E) with a single supplement. The serendipity and mixed families are related by de Rham theory: We define one family of direct serendipity elements that is the precursor of the reduced and full AC spaces. We also define many fully direct serendipity elements that use no mappings to define S DS r (E), which in turn generate new reduced and full direct H (div) approximation mixed spaces that use no mappings whatsoever.

Outline of the paper
We set some basic notation in the next section. We construct new families of direct serendipity elements in Sect. 3 for any index r ≥ 2, and in Sect. 4 for index r = 1. Our development is constructive, and results in a local nodal basis. These serendipity spaces do not involve mappings from a reference element. In Sect. 5, we define direct serendipity spaces based on supplements that are mapped from a reference element (when r ≥ 1). We discuss the stability and approximation properties of the new direct serendipity elements in Sect. 6. In Sect. 7, we turn our attention to the construction of direct mixed finite elements through a de Rham complex. We recover the spaces AC r −1 r (E) and AC r r (E) before defining new direct mixed finite element spaces, which do not require mappings from a reference element. Implementation is discussed, either using the hybrid method, or using an H (div)-conforming global basis, which is constructed. We discuss the stability and approximation properties of the new mixed elements in Sect. 8. Some numerical results illustrating the performance of our new direct serendipity and mixed finite elements appear in Sect. 9. Finally, a summary of our results and conclusions are given in the last section of the paper. Moreover, a second de Rham complex involving the gradient and curl operators provides new H (curl) = v ∈ (L 2 ) 2 : curl v ∈ L 2 elements as well.

Some notation
Let P r (ω) denote the space of polynomials of degree up to r on ω ⊂ R d , where d = 0 (a point), 1, or 2. Recall that LetP r (ω) denote the space of homogeneous polynomials of degree r on ω. Then Let the element E ⊂ R 2 be a closed, nondegenerate, convex quadrilateral. By nondegenerate, we mean that E does not degenerate to a triangle, line segment, or point. We choose to identify the edges and vertices of E adjacently in the counterclockwise direction, as depicted in Fig. 1. Let the edges of E be denoted e i , i = 1, 2, 3, 4, and the vertices be x v,1 = e 1 ∩ e 2 , x v,2 = e 2 ∩ e 3 , x v,3 = e 3 ∩ e 4 , and x v,4 = e 1 ∩ e 4 . Let ν i denote the unit outer normal to edge e i , and let τ i denote the unit tangent vector of e i oriented in the counterclockwise direction, for i = 1, 2, 3, 4. Let the reference elementÊ be [−1, 1] 2 . Define the bilinear and bijective map F E :Ê → E that maps the vertices ofÊ to those of E, oriented with (−1, −1) being mapped to x v,1 . We define the linear polynomial λ i (x) giving the distance of x ∈ R 2 to edge e i opposite the normal direction. It is where x * i ∈ e i is any point on the edge. If x is in the interior of E, these functions are strictly positive, and each vanishes on the edge which defines it.
We denote by F 0 E the pullback map associated with F −1 E ; that is, F 0 E is the map taking a functionφ defined onÊ to a function φ defined on E by the rule where x = F E (x). We denote by F 1 E the Piola map taking a vector functionψ ψ ψ defined onÊ to a vector function ψ ψ ψ defined on E by the rule where D F E (x) is the Jacobian matrix of F E and J E is its absolute determinant. Recall Ciarlet's definition [23] of a finite element.

Definition 1 (Ciarlet 1978) Let
1. E ⊂ R d be a bounded closed set with nonempty interior and a Lipschitz continuous boundary, 2. P be a finite-dimensional space of functions on E, and 3. N = {N 1 , N 2 , . . . , N dim P } be a basis for P .
Then (E, P, N ) is called a finite element.
Our task is to define the shape functions P and the degrees of freedom (DoFs) N . The DoFs give a basis for P provided that we have unisolvence of the shape functions (i.e., for φ ∈ P, N j (φ) = 0 for all j implies that φ = 0). To achieve optimal approximation properties, we will require that P ⊃ P r (E) for each index r . That is, the polynomials will be directly included within the function space, and hence we call our new finite elements direct serendipity and direct mixed finite elements.
Let Ω ⊂ R 2 be a connected, polygonal domain with a Lipschitz boundary (i.e., Ω has no slits), and let T h be a conforming finite element partition or mesh of Ω into nondegenerate, convex quadrilaterals of maximal diameter h > 0. The DoFs must be defined so that the shape functions on adjoining elements merge together. For serendipity spaces, we want the global space to reside in H 1 (Ω), so the elements must merge continuously across each edge e. For mixed spaces, the vector variable must lie in H (div; Ω), which means that the normal components (fluxes) of the vectors on an edge e in adjacent elements must be continuous. Table 1 Geometric decomposition and number of degrees of freedom (DoFs) associated to each geometric object of a quadrilateral for a serendipity element of index r ≥ 2.

Dimension
Object name  Object count  DoFs per object  Total DoFs   0  Vertex  4  1  4 1

Direct serendipity elements when r ≥ 2
We present our definition of direct serendipity elements when r ≥ 2 in this section. The case r = 1 requires special treatment, and will be given below in Sect. 4. Our dual objectives are that P r (E) ⊂ DS r (E) and that shape functions on adjoining elements merge continuously, i.e., so the space over Ω satisfies DS r (Ω) ⊂ H 1 (Ω). These objectives require us to consider the lower dimensional geometric objects within E (as in [6]). The minimal number of DoFs associated to each lower dimensional object must correspond to the dimension of the polynomials of degree r that restrict to that object. These numbers are given in Table 1. A quadrilateral has 4 vertices, 4 edges, and one cell of dimension 0, 1, and 2, respectively. Each vertex requires dim P r (R 0 ) = 1 DoF, the interior of each edge requires dim P r −2 (R) = r − 1 DoFs (not counting the two vertices), and the interior of each DoFs (not counting the edge and vertex DoFs). There are cell DoFs only if r ≥ 4, but the formula works for r ≥ 2. The total number of DoFs is then and so to define DS r (E), we will supplement P r (E) ⊂ DS r (E) with the span of two linearly independent functions, φ s,1 (x) and φ s,2 (x). We have many choices for the supplemental functions, the span of which is denoted S DS r (E) = span{φ s,1 , φ s,2 }. Each choice gives rise to a distinct family of direct serendipity elements of index r ≥ 2; that is, the shape functions (P in Definition 1) are We define the DoFs (N in Definition 1) as a set of nodal functionals N i defined at a nodal point x node j , i.e., As depicted in Fig. 2, for vertex DoFs, the nodal points are exactly the vertices x v,1 , 3 , and x v,4 of E. We have choices for the location of the rest of the nodal points.  For edge DoFs, we simply choose nodal points so that they, plus the two vertices, are equally distributed on each edge. There are r − 1 nodal points on the interior of each edge, which can be denoted x e,i j , j = 1, . . . , r − 1, for nodal points that lie on edge e i , i = 1, 2, 3, 4, ordered in the counterclockwise direction. The interior cell DoFs can be set, for example, on points of a triangle T strictly inside E, where the set of nodal points is the same as the nodes of the Lagrange element of order r − 4 on the triangle T . We denote the interior nodal points as x E,i , i = 1, . . . , 1 2 (r − 2)(r − 3). The total number of nodal points is indeed D r .
We will define a basis for the shape functions P dual to N , which will give unisolvence and a properly defined finite element. Such shape functions are called nodal basis functions. For a nodal point x node j , they have the property that N j (ϕ i ) = ϕ i (x node j ) = δ i j , the Kronecker delta. But first, we define the supplemental functions.

Supplemental functions
As stated above, we define distinct families of direct serendipity finite elements depending on a choice of two supplemental functions. These will be defined by a choice of four functions, two of which are linear polynomials, denoted λ 24 (x) and λ 13 (x). The other two functions should be continuous on E (and so bounded), and they are denoted R 13 (x) and R 24 (x). The supplemental functions are then φ s,1 = λ 2 λ 4 λ r −2 24 R 13 and φ s, When r = 2, λ 24 and λ 13 are not needed. The linear function λ 24 is defined by its zero set line L 24 . As shown in Fig. 3, let L 1 and L 3 be the infinite lines containing the edges e 1 and e 3 , respectively. We require that L 24 is chosen to intersect both L 1 and L 3 . Moreover, when e 1 and e 3 are not parallel, L 1 and L 3 intersect in a point x 13 = L 1 ∩ L 3 , and we also require that x 13 / ∈ L 24 (i.e., so that λ 24 (x 13 ) = 0). In a similar way, L 13 , the zero set line of λ 13 , is chosen to intersect the lines L 2 and L 4 extending e 2 and e 4 , respectively, and when they are not parallel, L 13 must avoid the intersection point x 24 = L 2 ∩ L 4 . Let ν 24 and ν 13 denote a unit normal to L 24 and L 13 , and let x * 24 and x * 13 denote any point on these lines, respectively. Then we define This definition is very general, but it is sufficient to provide a well defined finite element; however, accuracy considerations require a more restrictive definition, such as that given in Lemma 1 (where L 24 should intersect e 2 and e 4 , and L 13 should intersect e 1 and e 3 ). We remark that a simple choice is to take although the normalization is not strictly necessary. Later in (97) we will need the tangent vectors τ 24 and τ 13 , defined counterclockwise from the normals ν 24 and ν 13 , respectively.
The continuous functions R 13 and R 24 are defined to satisfy the properties They are ±1 on opposite edges, but arbitrary on the other two edges. A smoothness requirement will be added later in Lemma 1. A simple choice is to take the rational functions R simple 13 and R simple 24 (note that the denominators do not vanish on E). One could also use the bilinear map Theorem 1 Let E be a nondegenerate, convex quadrilateral. Suppose that λ 13 and λ 24 are linear functions with zero lines that intersect, for λ 13 , the lines containing e 1 and e 3 , and for λ 24 , the lines containing e 2 and e 4 , but avoiding the intersection points if they exist. Suppose also that the bounded functions R 13 and R 24 are continuous and satisfy (17). Then for r ≥ 2, with nodal DoFs defined by (12), is a well defined finite element. Moreover, DS r (E) has the minimal possible dimension (10) needed for H 1 conformity, and so it is a direct serendipity finite element.
It remains only to prove that the DoFs are unisolvent (i.e., that N is a basis for the dual space). We will do this by constructing explicitly a nodal basis, which could be used in practical implementations if one wished to do so.

Remark 1
If E is a rectangle, the classic serendipity spaces arise from our construction using a specific set of choices. For simplicity, assume E =Ê = [−1, 1] 2 is the square oriented with e 1 ⊂ {(−1,ŷ) :ŷ ∈ R} being on the left (every other rectangle is affine equivalent toÊ). To recover the classic serendipity space S r (E), take in DS r (E) λ 13 =ŷ + c 13 and λ 24 =x + c 24 for some constants c 13 and c 24 , and take R 13 =x and R 24 =ŷ, which is the simple choice (18). All other choices give new serendipity finite elements on the square, to the best of our knowledge. However, since these are defined in a nonsymmetric way, they are probably of little interest.

Proof of Theorem 1: construction of nodal basis functions
We define so that R i is 1 on edge e i , 0 on the opposite edge, and arbitrary on the other two edges. Note that

Interior cell nodal basis functions
For the entire cell E, we have interior shape functions only when r ≥ 4 (recall Table 1). These shape functions are and they vanish on all four edges (i.e., at all edge and vertex nodes). Let {φ E,i } ⊂ P r −4 be a nodal basis for the cell nodes Our interior cell nodal basis functions are then

Edge nodal basis functions
We construct ϕ e,11 (x), which is 1 at x e, 11 and vanishes at all other nodal points. The construction of the other edge nodal basis functions is similar.
This function vanishes on all edges but e 1 . Let We want q to vanish at the nodes x e,1i for i = 2, . . . , r − 1; that is, we want p to satisfy the conditions (note that λ 3 (x e,1i ) = 0). These r − 2 conditions uniquely define p along the edge e 1 , i.e., they definep ∈ P r −3 (e 1 ) as a function of t = x · τ 1 . The coefficients can be found using Newton's divided difference interpolation formulas for the points t i = x e,1i · τ 1 (i = 2, . . . , r − 1) and values given on the right hand side of (25). For example, the low order cases arẽ 13 ) We define p(x) by extendingp(t) to E constantly along perpendicular lines, i.e., Our construction will succeed provided q(x e,11 ) = 0. So suppose to the contrary that q(x e,11 ) = 0. We restrict x to L 1 (the line extending e 1 , as in Fig. 3) and let t = x · τ 1 . Conversely, given t, there is a unique x ∈ L 1 such that We have assumed thatq(t 1 ) = 0, soq(t i ) = q(x e,1i ) = 0 for all i = 1, . . . , r −1. That is,q(t) is a polynomial of degree r − 2 vanishing at r − 1 points, and so it vanishes identically. We have two cases to consider. First, suppose that the lines through e 1 and e 3 intersect at x 13 (see Fig. 3). Now is a contradiction, since λ 3 (x 13 ) = 0 and λ 24 (x 13 ) = 0 by our choice of this linear function. Second, suppose that the lines through e 1 and e 3 are parallel. Then λ 3 | e 1 = α > 0 is a strictly positive constant, and so This is clearly a contradiction, since the zero line of λ 24 is transverse to e 3 (again by our choice) leading us to conclude thatλ r −2 24 must have strict degree r − 2. We have concluded that q(x e,11 ) = 0, and so also φ e,11 (x e,11 ) = 0. We complete the construction by defining 11 ) .

Vertex nodal basis functions
For the vertices, r ≥ 2, so we can define the shape functions wherein we interpret indices modulo 4. These four functions vanish at all of the edge nodes, and φ v,i (x v, j ) = 0 if i = j and is positive otherwise. The nodal basis functions are then This completes the construction of the D r = dim P r (E) + 2 (recall (10)) nodal basis functions for DS r (E). This also completes the proof of Theorem 1.

Implementation as an H 1 -Conforming Space
On the mesh T h of Ω, the global direct serendipity finite element space of index r is Because our elements are polynomials of degree r on the edges, they merge together continuously, provided that their edge and vertex DoFs match on element boundaries. We constructed a local nodal basis for DS r (E), r ≥ 2; that is, one for which every basis function vanishes at all but one nodal point, and it equals to one at this point. These local basis functions (after extension by zero outside the element), merge together continuously to give a global nodal basis for In this way we construct global nodal basis functions, each equal to one at a nodal point and zero at all the other nodal points (so far for r ≥ 2, but also for r = 1, after the construction of the next section is complete). element S 1 (E) is the tensor product space of bilinear functions P 1,1 (Ê) onÊ mapped to E by F 0 E , and, in fact, has the form of a direct serendipity space with only one supplemental function.
for some supplemental function R that reduces to a linear function on each edge of E and satisfies As just noted, the mapped function gives the usual parametric serendipity space, and it is linear on each edge of E and satisfies (34).
Proof We first show that (33) with any such R will give a well defined serendipity finite element DS 1 (E) by showing that it has a nodal basis. The nodal basis is constructed by first defining the linear functions with zero lines corresponding to the diagonals of the element E. Let ν d,1 be either unit normal to the first diagonal joining x v,1 with x v,3 , and let ν d,2 be either unit normal for the second diagonal joining x v,2 with x v,4 . Then let The nodal basis functions are Note that there is no division by zero, so these functions are well defined, and that they are in DS 1 (E) = P 1 ⊕ span{R}, as required.
For the direct implication of the theorem, every well defined direct serendipity space DS 1 (E) has a nodal basis of dimension four (recall Table 1) and contains P 1 (E).
is in the space and cannot be in P 1 (E), we conclude that DS 1 (E) has the form (33).
One way to define R is to use generalized barycentric coordinates (GBCs) [29,30,33,40]. There are many types of GBCs, including Wachspress, mean value, Sibson, and harmonic coordinates. The functions ϕ i (x) : E → R, i = 1, 2, 3, 4, are GBCs if they satisfy the two properties: These four functions are linearly independent, they are linear on each edge e of ∂ E, their span includes P 1 (E), and they form a nodal basis with respect to the vertices, i.e., ϕ i (x v, j ) = δ i j for all i, j [30,33]. So by their definition, their span is a direct serendipity finite element, and moreover they constitute the nodal basis.
However, we do not require the non-negativity property. Functions satisfying linear completeness for L ∈P 1 (E) are called homogeneous coordinates, and they were completely characterized in [30] in terms of areas of triangles. After normalization, these give direct serendipity finite elements. In the technical sense of their characterization (which requires the choice of four functions), no one can find new spaces. However, our new characterization of the spaces (33) can be used to give an alternate construction that is based on simple linear functions.
Our idea is to construct R inside DS 2 (E) so that R satisfies (34). There are many ways to define DS 2 (E), so we get many R's, a different one for each choice of the space DS 2 (E). Let ϕ (2) e,i1 (x) be the edge nodal basis function for the node x (2) e,i1 in DS 2 (E), i = 1, 2, 3, 4. It is quadratic on each edge. Let which is quadratic as well, and then define r (x (2) e,i1 ) ϕ (2) e,i1 (x). Then Restricted to the edges, R is nominally quadratic, but it reduces to a linear polynomial on each edge, since R is 1 at one vertex, −1 at the other, and vanishes at the midpoint, i.e., at x (2) e,i1 for all i.

Serendipity supplements based on mapping from a reference element
There are other ways to define serendipity finite elements. In this section, we define the supplemental space on the reference elementÊ and map it to E using the bilinear map.
This construction gives us direct serendipity elements with mapped supplements. For r ≥ 2, we must show unisolvence with this supplemental space. As before, we show this property by constructing a nodal basis.
Since the supplements are not used in the construction of interior cell nodal basis functions, the definition (23) continues to be valid. Moreover, once the edge nodal basis functions are defined, the vertex nodal functions are defined by (30). Thus, we need only construct the edge nodal basis functions. As in Sect. 3.2.2, we discuss only the nodal basis function ϕ e,11 (x) at nodal point x e,11 , since the other edge nodal basis functions are constructed similarly. (19) and which is a nonlinear function. However, because F E is a bilinear map, on the edges e 1 and e 3 , λ * 13 is linear and F 0 E (1 −x 2 2 ) is quadratic (but these may be different linear and quadratic functions on the two edges).
The function λ * 24 has the zero set being the line joining the center of e 1 to the center of e 3 . Let x * 24 be any point on this line and ν 24 denote a unit normal to the line. Define the linear function which mimics λ * 24 in the sense that both are linear on e 1 and e 3 , and they have the same zero set. Then there are constants a i = 0 of the same sign such that In fact, is quadratic and vanishes at the ends of e 1 and e 3 , so for constants b i > 0, defined by considering the center points of e 1 and e 3 .
This function vanishes on e 2 and e 4 . It also vanishes on e 3 , since there R 13 = 1, and so Restricted to e 1 , we have R 13 = −1 and is a well defined finite element.

Approximation properties of DS r
In this section, we develop the stability and convergence theory for our new direct serendipity finite elements. For the most part, we work over the entire domain Ω, with the assumption that diam(Ω) = 1 for simplicity. To obtain global approximation properties, we need to assume that the mesh is uniformly shape regular [34,pp. 104-105], which means the following.

Definition 2
For any E ∈ T h , denote by T i , i = 1, 2, 3, 4, the sub-triangle of E with vertices being three of the four vertices of E. Define the parameters A shape regular mesh has a bound on the number of quadrilaterals that can share a single vertex.
In Sects. 3 and 4 , we constructed local and global nodal bases for DS r (E) and DS r = DS r (Ω) for r ≥ 1; that is, one for which every basis function vanishes at all but one nodal point, and it equals to one at this point. In this section, we denote the set of global nodal basis functions as {ϕ 1 , . . . , ϕ N r }, corresponding to global nodal points {a 1 , . . . , a N r }, respectively, where N r = dim DS r .

An interpolation operator mapping onto DS r
We first construct an interpolation operator mapping onto DS r following Scott and Zhang [42]. For each node a i in the interior of some element E ∈ T h , we set K i to be (the closed set) E, and we call such a node an interior node. For each node a i in the interior of edge e of T h (i.e., not at the vertices), we set K i = e (a closed set). These nodes are called edge nodes. For each node a i being a vertex of T h , we choose K i to be any fixed edge e containing a i , with the restriction that if a i ∈ ∂Ω, then e ⊂ ∂Ω. Note that e is chosen from among multiple edges. These nodes are called vertex nodes.
We define a special L 2 -dual nodal basis {ψ 1 , . . . , ψ N r } as follows. For each node a i , we denote the total number of nodes in K i as n i , and then denote these nodes in where we use a slight abuse of notation in that dx should be dσ for edge and vertex nodes. Finally, for the node a i , we take ψ i = ψ i,1 . (As described, this construction is highly redundant. Since it is used only for theoretical purposes, we do not explore its efficient implementation.) For any node a i giving rise to K i and ψ i , This is easily seen as follows. If a i is an interior node, then this expression is exactly (52) (since the latter expression holds for all nodes on E). If a i is an edge or vertex node, then when a j is also an edge or vertex node, (53) is (52), and when a j is an interior node, ϕ j vanishes on the edges of T h (and thus on K i ).
We can now define an interpolation operator I r h : W l p (Ω) → DS r by where 1 ≤ p ≤ ∞ and l > 1/ p (but l ≥ 1 if p = 1). By the trace theorem, the nodal values , even when K i is an edge. Note that I r h depends on our choice of K i , but we suppress this fact in the notation for simplicity. Because (53) holds, this operator is a projection on DS r .

Boundedness of the interpolation operator I r h
To obtain approximation properties, the interpolation I r h needs to be a bounded operator. Scott and Zhang's proof of this fact in their situation [42] does not hold directly for our construction, since we need to use non-affine mappings from the reference element to the actual elements, and we use non polynomial shape functions. We give a proof of boundedness based on a continuous dependence argument.
On an element E ∈ T h , it is clear that the linear functions λ i defined in (7) depend only on x and the vertices of E, and that this dependence is continuous. For simplicity of discussion, we use the simple choices given in (16) and (18) and consider only the direct serendipity finite elements based on this choice. We will remove this restriction at the end of the section. However, for these simple choices, it is clear that these four functions used to define the serendipity finite elements are continuously differentiable functions of x and the vertices of E.
We need to fix the domain to the reference elementÊ = [−1, 1] 2 . For any E ∈ T h , let H = H E be the maximal edge length of E. (By shape regularity, H E is comparable to h E , the diameter of E.) Since our finite element construction is invariant under translation and rotation, for any E ∈ T h we can assume that x v,1 = (0, 0) and x v,2 = (H , 0), as depicted in Fig. 4. Furthermore, we can scale E by 1/H to define a local reference elementẼ = E/H with vertices (0, 0), (1, 0), (v 1 , v 2 ), and  (v 3 , v 4 ). We can view E as the image of the bilinear map F E = H FẼ defined in Sect. 2, which is a continuous function ofx ∈Ê and the vertices of E (since E is nondegenerate). Boundedness of this mapping is well-known for shape regular meshes. To be more precise, let F K i = F E | K i when K i = E or K i is an edge of some E. By the uniform shape regularity of the mesh, for some constant C independent of K i [23,Theorem 4.3.3], These are the properties of the mapping needed in the argument of Scott and Zhang [42]. Further, since the bilinear mappings F E are defined on nondegenerate quadrilaterals, F E and F −1 E are smooth, and similar bounds hold for higher order derivatives. In particular, higher order derivatives of FẼ and F −1 E are uniformly bounded. The edge nodal points have been placed uniformly on each edge. We need to fix a place for the interior nodal points, so that their positions vary continuously with the location of the vertices. Recall that the interior nodal points correspond to Lagrange nodal points for P r −4 . We can fix these on a triangle with vertices at, say, The set of admissible ξ is bounded, since no side ofẼ has length greater than 1. Moreover, the shape regularity constraints are given by (albeit complicated) continuous functions of (v 1 , v 2 , v 3 , v 4 ) involving maximal inscribed circles being required to lie in the closed interval [σ * , ∞). Therefore the set of admissible ξ is a closed, and hence compact, set. We conclude that each ϕẼ j and its derivatives are bounded uniformly with respect to the shape ofẼ. That is, there exists a constant C = C(σ * , m, q) independent of ξ such that where σ * is the shape regularity parameter. Since locally ϕẼ j (x/H ) = ϕ E j (x), we also have a bound for the global nodal basis functions, namely, We also need to show that the dual basis functions are bounded such that for any node a i , which is a result analogous to [42,Lemma 3.1]. Here, the dual nodal basis functions defined for nodes inẼ as defined in Sect. 6.1. These are also continuous functions of ξ (and possibly the corresponding values of ξ for its neighboring elements, due to the treatment of vertex nodes). By a similar continuity and compactness argument, we conclude that ψẼ j is bounded uniformly with respect to the shape ofẼ and its neighbors. When a i is an interior node,K i =Ẽ, so and so (59) holds. When a i is an edge or vertex node,K i and K i are edges which are affinely related to the reference interval [−1, 1]. Moreover, the functions in question, when restricted to the edge, are polynomials. Thus Scott and Zhang's argument holds directly, and so (59) holds in general. This and (57) lead to the conclusion that the interpolation operator is bounded.
We can extend the proof to more general λ 24 and λ 13 . These linear functions are defined by their zero lines, which are defined by two points each. However, we have a restriction that for λ 24 , say, the zero set line L 24 must intersect both L 1 and L 3 , but not at their intersection when they are not parallel. This choice of four points (two for each of λ 24 and λ 13 ) could be added as parameters to the variable ξ . However, then the restriction implies that ξ does not vary over a compact set. So we must restrict the choice of these four points. We make a simple (and practical) requirement that the zero set line L 24 intersects e 1 and e 3 , and L 13 intersects e 2 and e 4 . We actually choose four points (which are added into ξ ), one on each fixed and closed edge ofÊ, and map them through F E to define λ 24 and λ 13 . In this way, λ 24 and λ 13 are defined as continuous functions of ξ , and ξ varies over a compact set, and the argument above continues to hold.
We can generalize the possible R 13 and R 24 that can be used as well. Of course they must satisfy the requirement (17), but they must also be uniformly differentiable functions of ξ . Such is the case for the mapped choice R

Lemma 1
Let v ∈ W l p (Ω), with 1 ≤ p ≤ ∞ and > 1/ p (or ≥ 1 if p = 1). Let T h be uniformly shape regular (Definition 2) with shape regularity parameter σ * . For every E ∈ T h , suppose that the basis functions of DS r (E) are constructed using λ 24 and λ 13 such that the zero set L 24 intersects e 1 and e 3 , and L 13 intersects e 2 and e 4 . Moreover, assume that R 13 and R 24 are uniformly differentiable functions of the vertices of E up to order m. Then for r ≥ 1, E ∈ T h , 1 ≤ q ≤ ∞, and any nonnegative integer m, where E * = F∈T h , F∩E =∅ F and | · | W k p is the seminorm of k-th order derivatives. We remark that the mapped supplements (43) vary continuously with the element shape, and so satisfy the lemma above.

Approximation properties of the DS r spaces
We use the Bramble-Hilbert lemma [18] in the form developed by Dupont and Scott in [28] (see also [19]). A domain ω is star-shaped with respect to a ball B r of radius r if for all x ∈ ω, the closed convex hull of {x}∪ B r is a subset of ω. Let r max = sup{r : ω is star-shaped with respect to B r } and h = diam(ω), and define the chunkiness parameter of ω by h/r max . Then the Bramble-Hilbert Lemma, i.e., that has a constant C that depends continuously on the chunkiness parameter. On the local reference element ω =Ẽ (actually its interior), the chunkiness parameter varies continuously in a compact set due to the shape regularity assumption, so the constant C has an upper bound independent of the vertices ofẼ. Combining Lemma 1 and the Bramble-Hilbert lemma (61), we derive our theorem for local and global error estimation using E * defined just after (60).

Construction of direct mixed finite elements using a de Rham complex
The de Rham complex of interest here is where the curl (or rot) of a scalar function φ( . From left to right, the image of one linear map is the kernel of the next. On rectangular elements, it is known [6,7] that the serendipity space S r +1 is the precursor of the Brezzi-Douglas-Marini space BDM r [20] for r ≥ 1; that is, on the reference squareÊ, (1) holds.

Reduced and full AC spaces
We have the following extension of (1) to quadrilateral elements E. The direct serendipity spaces DS mapped r (49) using the mapped supplements (43) is the precursor of the reduced H (div)-approximating Arbogast-Correa space AC r −1 r [1], r ≥ 1, defined on meshes of convex quadrilaterals: Moreover, the full H (div)-approximating space AC r r , for r ≥ 0, satisfies This observation is clear once one realizes three sets of facts. First, the direct serendipity elements based on (43) have the structure Second, the AC elements have the structure where F 1 E is the Piola mapping (9) fromÊ to E. Finally, we have the fairly well-known Helmholtz-like decomposition (see, e.g., [1]) the relation between the curl operator and the bilinear and Piola maps and the fact that the div operator takes xP s one-to-one and onto P s for any s ≥ 0. Now we see from (73) that, when r ≥ 1, and so is in the kernel of the operator div. Finally, (72) says that These spaces satisfy the exact sequence properties of the de Rham complex (65)-(66). For r = 0, it is easy to check that S 1 (E) = DS 1 (E) (see (32)) precedes the element AC 0 0 (E) in the de Rham sequence (66).

Direct mixed finite elements when r ≥ 1
According to [1], a reduced or full H (div)-approximating mixed finite element space defined directly on a quadrilateral E of minimal local dimension takes the form (P in Definition 1) where, in that paper, the choice of S V r (E) is given by taking (71). However, it is noted that other supplemental functions could be used [1,near (3.15)]. Their normal components must lie in P r (e i ) on each edge e i and, if they are mapped by the Piola transform, they must contain a nontrivial component of the DoFs of curlx r +1 1x 2 and curlx 1x r +1 2 . For r ≥ 1, the supplemental space S V r (E) must be of dimension two and linearly independent of P 2 r (E), so that dim V r −1 r (E) = (r + 2)(r + 1) + 2 = r 2 + 3r + 4, dim V r r (E) = (r + 2)(r + 1) + (r + 1) + 2 = r 2 + 4r + 5. (79) As given in [1], the DoFs (N in Definition 1) for ψ ψ ψ ∈ V s r (E), s = r − 1, r , are given (after fixing a basis for the test functions) by where dσ is the one dimensional surface measure and the H 1 (E) and H (div; E) bubble functions, for r ≥ 3, are Note that the number of DoFs agrees with the dimension of the space (79). Unlike the construction given in [1], which only considered mixed finite elements, we use here the de Rham theory to construct a mixed finite element space V s r based on a well defined direct serendipity space. For r ≥ 1, we have de Rham complexes for both reduced and full direct H (div)-approximating mixed elements: for any variant of our new direct serendipity spaces. We define spaces of vector functions according to these de Rham complexes, using the fact that the div operator takes xP k one-to-one and onto P k . These spaces are and they have the form (77)-(78), provided we define Theorem 4 Let E be a nondegenerate, convex quadrilateral and DS r +1 (E) = P r +1 (E) ⊕ S DS r +1 (E) be a well defined direct serendipity finite element for r ≥ 1.

with DoFs defined by (80)-(82) are well defined finite elements. Moreover, for s = r −1, r, V s r (E) has the minimal possible dimension (79) needed for H (div) conformity and the property that div V s r (E) = P s (E).
Proof It remains only to show that in fact these spaces are unisolvent for the DoFs (80)-(82); that is, that these spaces are well defined mixed finite elements. So suppose that ψ ψ ψ ∈ V s r (E), s = r − 1, r , and has vanishing DoFs. By construction, where φ ∈ DS r +1 (E) and ψ ψ ψ d ∈ xP s . The DoFs (81), with (80), imply that for q ∈ P s (E), Since ∇ · ψ ψ ψ d ∈ P s (E), we conclude that ∇ · ψ ψ ψ d = 0, and further that ψ ψ ψ d = 0. We next observe the well known fact that tangential derivatives of functions along the edges of E map by the curl operator to normal components; that is, if we define the (counterclockwise) unit tangential vector then for φ ∈ DS r +1 (E), Therefore, the DoFs (80) imply that for any edge e i and p ∈ P r (e i ), and we conclude that ∇φ · τ = 0 on ∂ E, and so φ = c is constant on Thus ψ ψ ψ = 0, and the proof of unisolvence is complete.
We can use any of our direct serendipity spaces to define direct mixed finite elements V r −1 r (E) or V r r (E). As we saw earlier, if we use DS mapped r +1 (E), we recover the known finite elements AC r −1 r (E) and AC r r (E). However, the direct serendipity spaces of Sect. 3 give new (direct) mixed finite elements.

Direct mixed finite elements when r = 0
When r = 0, there are no reduced H (div) approximation finite element spaces. The full H (div) approximation direct elements have dimension 4, and they take the form where R was defined in Sect. 4 above as F 0 E (x 1x2 ) or in (42) using the special construction in DS 2 (E). The former is the space AC 0 0 [1], and the later appears to be a new mixed finite element satisfying (93)

Implementation using the hybrid mixed method
The mixed space of vector functions V s r over Ω is defined by merging continuously the normal fluxes across each edge e of the mesh T h . That is, for r ≥ 0, s = r − 1, r , s ≥ 0, This space is normally paired with a space approximating scalar functions. When, say, solving a second order elliptic partial differential equation in mixed form, the vector functions V s r are paired with the scalar functions These scalar spaces are the divergences of the corresponding vector function spaces. However, in practical implementation, the hybrid form of the mixed method is often used [10]. In that case, the elements V s r (E) are simply concatenated, and no globally merged basis is required. The Lagrange multiplier space, used to enforce the normal flux continuity through an additional equation, is simply To represent the vector functions in V s r as we presented them, it would appear that we need to apply the curl operator extensively. However, this is not the case, since the vector polynomials in these spaces are clear, so we need only apply the curl operator to the supplemental functions. Even then, taking a curl can be avoided in some cases. As we saw earlier, curl S DS ,mapped r +1 (E) gives the supplements for the known finite elements AC s r (E), s = r − 1, r , which are computed in (71) using the Piola transform rather than the curl operator.
For example, suppose we use the direct serendipity elements from Sect. 3 (so r ≥ 1), with supplements defined by (14). We note that where j = 1, 2, 3, 4 or j = (24), (13 The supplemental vector functions are then

Implementation as an H(div)-conforming mixed space
The H (div) spaces V r −1 r and V r r defined in (94) can be given a global basis with the normal components of the shape functions merged continuously on each edge of the mesh T h . We present a method that uses the nodal basis functions of the serendipity space of index r + 1 preceding V r −1 r or V r r in the de Rham complex. Since the H (div) bubble functions B V r (83) have no normal flux, they can be handled easily. For each E ∈ T h , when r ≥ 3 one can define the global basis functions using the interior cell nodal basis functions (23) with a superscript as a reminder that the index of the direct serendipity space is r + 1. The shape functions arising from the edge nodal basis functions, like (28), also present no particular difficulty. For each edge e of the mesh shared by elements E k and E (with e being locally edge i 1 and i 2 , respectively), when r ≥ 1 one can define These global basis functions are in H (div) because the serendipity elements are continuous, so the tangential derivatives (i.e., the flux-see (91)) agree across e. Moreover, the average normal flux vanishes on each edge. The most delicate global basis functions to construct are those for which the average normal flux is a constant on each edge e of the mesh. For each r ≥ 0, these are given primarily, but not completely, by taking curls of the serendipity vertex nodal basis functions (30). However, each of these functions has normal flux on two edges, and there are only three linearly independent functions per element (since the kernel of curl consists of the constant functions). For each element E having e as an edge, we need to consider the vector functions P 2 0 (E)⊕xP 0 (E) with these curls. The proper construction requires some work on each element E of the mesh. We begin the construction by = δ ik and its restriction to each edge of E is a linear function. To be precise, .
We also define ψ ψ ψ * * v, Finally, for each edge e of the mesh which is edge i of element E, we define Finally, when s ≥ 1 we define the global basis functions associated to the nonconstant divergences. These functions are local to each element E ∈ T h . Working on E, we begin with the functions xP * s (E), where P * s (E) = s k=1P k (E) ⊂ P s (E). Take p i (x) in a basis for P * s (E), so i = 1, . . . , 1 2 (s + 2)(s + 1) − 1. We must remove the normal flux on ∂ E from x p i (x). We do this using (99) and (100) by defining and setting the coefficients α j,k on each edge e j so that where c j = x · ν j | e j is a constant. The coefficients can be found in a number of ways, including a straightforward application of linear algebra requiring the solution of four small (r + 1) × (r + 1) linear systems. An alternative can be given once one realizes that on edge e j , for k ≥ 1, and the coefficients can be read off by substituting in the Lagrange points t = /(r +1) for = 1, . . . , r + 1. The global basis is now fully defined.

Stability and approximation properties for V s r
In this section, we develop the stability and convergence theory for our new direct mixed finite elements. Again, we assume that diam(Ω) = 1 for simplicity. As was done by Raviart and Thomas [41] for their mixed spaces, we can define a projection operator π : H (div; Ω) ∩ (L 2+ (Ω)) 2 → V s r , s = r − 1, r , where > 0. The operator π is pieced together from locally defined operators π E . Following [1], we define π E v in terms of the DoFs (80)-(82). The operator π satisfies the commuting diagram property [25], which is to say that where P W s is the L 2 -orthogonal projection operator onto W s = ∇ · V s r . To show that certain important properties of π E do not depend on the vertices of E except for scale, we work over the scaled elementẼ, which was introduced in Sect. 6.2.
We need to fix a basis for the DoFs (80)-(82), so letp j (t) =t j andp j,k (x) =x j 1x k 2 where j, k are integers. For anyṽ ∈ H 1 (Ẽ), we define the linear functional NẼ β (ṽ) with index β as follows: Denote the set of all possible indices of β as B. By a continuity and compactness argument, similar to that given in Sect. 6.2, there exists a constant C > 0, independent of the vertices ofẼ, such that ∀β ∈ B, where · j,ω is the norm of W j 2 (ω) = H j (ω). By unisolvence, there exists a basis {φ φ φẼ β , β ∈ B} for V s r (Ẽ), such that Then πẼ can be defined as Note that φ φ φẼ β varies continuously with respect to vertices ofẼ, and so there exists a constant C > 0, such that Combining (105) and (107), there exists a constant C > 0, independent of the vertices ofẼ, such that By the boundedness of πẼ in H 1 , the Bramble-Hilbert Lemma (61), and usual scaling arguments, there exists a polynomialp ∈ P 2 where the constants do not depend on the vertices ofẼ. The L 2 -orthogonal projection P W s gives optimal approximation, and so also ∇ · (u − π u) = ∇ · u − P W s ∇ · u is optimally approximated. We have shown the following.

Lemma 2 With the assumptions of Lemma 1, there is a constant C
where s = r − 1 ≥ 0 and s = r ≥ 1 for reduced and full H (div)-approximation, respectively. Moreover, the discrete inf-sup condition holds for some γ > 0 independent of T h and h > 0.

Numerical results
We test our finite elements on the Dirichlet problem where the second order tensor a(x) is uniformly positive definite and bounded, and f ∈ L 2 (Ω). The problem can be written in the weak form: where (·, ·) is the L 2 (Ω) inner product. Setting we also have the mixed weak form: Find u ∈ H (div; Ω) and p ∈ L 2 (Ω) such that These weak forms give rise to finite element approximations. In view of Theorem 3 and Lemma 2, it is well known that the following theorem holds [19,21].

Theorem 5
With the assumptions of Lemma 1, there exists a constant C > 0, independent of T h and h > 0, such that where p h ∈ DS r (Ω) ∩ H 1 0 (Ω) approximates (115). Moreover, with s = r − 1, r, where (u h , p h ) ∈ V s r × W s approximates (117)-(118). We consider the test problem (113)-(114) defined on the unit square Ω = [0, 1] 2 with the coefficient a being the 2 × 2 identity matrix, i.e., we solve the Poisson equation. We use the method of manufactured solutions, taking the exact solution to be u(x, y) = sin(π x) sin(π y) and the source term is f (x, y) = 2π 2 sin(π x) sin(π y).
Solutions are computed on three different sequences of meshes. The first sequence, T 1 h , is a uniform mesh of n 2 square elements (two sets of parallel edges per element). The second sequence, T 2 h , is a mesh of n 2 trapezoids of base h and one pair of parallel edges of size 0.75h and 1.25h, as proposed in [8]. The third sequence, T 3 h , is chosen so as to have no pair of edges parallel. The first 4 × 4 meshes for each sequence are shown in Fig. 5. Finer meshes are constructed by repeating the same pattern over the domain.
Our computer program uses the deal.II library [14]. The integrals must be approximated using quadrature rules, since we have nonpolynomial basis functions. In general, one can use a rule based on squares mapped to the quadrilateral, or one can cut each quadrilateral into two triangles and use a rule suitable for triangles. To accurately approximate the bilinear form, the order of the quadrature rule should be at least 2r . Construction of the finite element basis requires some computation, as described in Sects. 3 and 7. If one uses parallel computing, the time cost for these routines can be scaled nearly perfectly in parallel, since they basically involve only local computations. Moreover, in a time dependent problem, the basis needs to be computed only once. We find that reducing the global number of degrees of freedom in a serendipity space versus a tensor product space, even at the expense of a slightly more expensive basis construction, is worthwhile [2,47].

Direct serendipity spaces
We present convergence studies for the fully direct serendipity spaces DS r using the elements defined in (11) and (14), i.e., the ones that use no mappings. We compare with the two spaces of elements given by mapping the local serendipity spaces S r (Ê) and the tensor product spaces P r ,r (Ê) to the mesh elements (the latter are simply called the P r ,r spaces).
We take the simple choice (16) for λ 24 and λ 13 . We do not quite use the simple choice (18) for R 13 and R 24 . Instead, we let

4), and then set
These functions are constant on each opposite pair of edges, but not ±1 (this makes no difference to the development presented above).
For an n × n mesh, the total number of degrees of freedom for P r ,r is (nr + 1) 2 = O(r 2 n 2 ), and for S r and DS r it is dim(S r ) = dim(DS r ) = (number of vertices) + (number of edges)(r − 1) Therefore, the total number of degrees of freedom for a serendipity space is asymptotically about half the size of that for a tensor product space of the same order.

Convergence order versus maximal element diameter h
We report the L 2 and H 1 -seminorm errors and convergence orders for r = 2, 3, 4, 5 as h decreases (i.e., as n increases) on mesh sequence T 1 h in Tables 2 and 3. The direct serendipity space DS r and the regular serendipity space S r coincide on the square mesh T 1 h . All three spaces show an (r + 1)-th order convergence in the L 2 -norm and an r -th order convergence in the H 1 -seminorm, as we should expect from theory (see Theorem 5). Our results are fully consistent with the recently reported results in [24] for the standard serendipity elements on rectangles. Tables 4 and 5 show the errors and orders of convergence for the trapezoidal mesh sequence T 2 h . The tensor product spaces P r ,r and DS r achieve the expected optimal convergence rates: (r + 1)-th order in the L 2 norm and r -th order in the H 1 -seminorm. The regular serendipity spaces S r have worse than optimal convergence rates in both norms (as was also observed in [8]). The errors and convergence rates for P r ,r , DS r , and S r on mesh sequence T 3 h are similar to those on T 2 h , so we omit showing them. Finally, we present the errors and convergence rates for the serendipity spaces DS mapped r using elements defined in (49), which have supplements mapped from the reference element. The results for r = 2, 3, 4, 5 on mesh T 2 h are shown in Table 6. One can see optimal convergence rates and errors that compare favorably with those of the fully direct spaces in Tables 4 and 5, although the latter are perhaps slightly better.
In all cases, for a fixed value of n, the errors for P r ,r are smaller than that for DS r . One might have expected to see such a result, since P r ,r has many more degrees of freedom at its disposal to better approximate the solution. However, the serendipity spaces are advantageous in terms of memory usage and time to solution (see [24]), so if one is given a single fixed mesh for solving the problem, it is more efficient to use the serendipity spaces. It is not uncommon in engineering applications to be given a fixed mesh, such as in petroleum engineering applications where the mesh is determined by the subsurface geology and geostatistical issues.

Convergence order versus the number of DoFs
If the mesh is not viewed as being fixed but can be varied as one wishes, it would be more appropriate to compare the errors with respect to the number of unknowns (i.e., the total number of DoFs). This is done in Fig. 6 for square meshes, where we plot the log (base 10) of the error versus the log of the square root of the number of DoFs,   which is a surrogate for h. One can see that the serendipity space is more efficient for r = 2, the two spaces are about equally efficient for r = 3, and the tensor product space is more efficient for r ≥ 4. This behavior was also observed on rectangular meshes in [24] for the standard serendipity spaces. Similar results are seen in Fig. 7 for the trapezoidal meshes T 2 h . This phenomenon can be explained as follows. Suppose that the solution u is analytic, so that locally on an element E, we may view it as an infinite dimensional polynomial. The DS r (E) finite elements are able to match only the monomials in u up to order r (see [28]). However, the P r ,r (E) finite elements are able to match many more monomials, up to the highest degree monomial x r y r . When r = 2, P r ,r (E) can only match the three extra monomials x 2 y, x y 2 , and x 2 y 2 , and no improvement in accuracy over DS r (E) was seen in our examples. However, as r increases, P r ,r (E) can match on the order of 1 2 r 2 more monomials than DS r (E), and some improvement in accuracy was seen.  To test this explanation, we changed the manufactured true solution to be a single monomial u(x, y) = x 5 y k , with k = 0, 1, 2, 3, 4, 5, and tested the case r = 4. Simple polynomials should better approximate the case k = 0, while tensor product approximations should handle the "diagonal" case k = 5 better. This is indeed the case, as shown in Fig. 8. The more advantageous choice of finite element changes from DS r to P r ,r as k increases, with fairly equivalent results seen for k = 2.
The next series of tests considers the manufactured true solution u(x, y) = log(x + ky + 1), which is analytic with all orders of monomials in its Taylor expansion. As k increases, the relative strength of the "diagonal" monomials decrease, and so we should expect that DS r performs better. This is the case, as shown in Fig. 9, where we test r = 2, 3, 4, 5, k = 1, 2, 4, and use a mesh generated by a random perturbation of the vertices of a rectangular mesh. We see that P r ,r is more efficient when k = 1,  DS r is more efficient when k = 4, and they are fairly equivalent when k = 2. Similar results are seen if one uses the trapezoidal mesh T 2 h , and other analytic functions such as square root, exponential, and trigonometric functions.
The previous test motivates us to consider a typical groundwater or petroleum engineering problem related to a linear flood. In such a problem, a pressure difference is imposed from left to right. In a homogeneous medium, the pressure field would vary linearly, but in a heterogeneous porous medium, the pressure is perturbed from a clean linear variation. To simulate such a flood, we imposed the manufactured true solution u(x, y) = 1 + 0.1 sin π cos(π x) cos(3π y/2) (1 − x). The convergence of the error versus the number of DoFs is shown in Fig. 10. In this test, DS r outperforms P r ,r for all r = 2, 3, 4, 5.
We conclude that if one is given a fixed number of DoFs (rather than a fixed mesh), it is difficult to predict whether it is more efficient in terms of accuracy to use P r ,r or DS r . However, the P r ,r elements do not lead to mixed finite elements with the appropriate order of convergence [1,8]. So even though there are cases where one might prefer to use P r ,r , it is nevertheless critical to understand the DS r finite elements, since they are needed to define the direct mixed finite elements in the de Rham sequence (84). We give numerical examples of these spaces next.

Fully direct mixed finite elements on quadrilaterals
We verify the convergence rates for the new fully direct mixed finite elements V s r × W s derived in Sect. 7.2. These are implemented without the use of any mapping from the reference element. We take the simple choices (16) and (18), although the choice (123) provides similar results. We apply the hybrid form of the mixed finite element method (Sect. 7.4). The errors and the orders of convergence when r = 1, 2 on mesh T 2 h are presented in Table 7. Again, results are similar on T 3 h meshes. As the theory predicts, the scalar p, the vector u, and the divergence ∇ · u obtain order of approximation s + 1, r + 1, and s + 1, respectively, for the reduced (s = r − 1) and full (s = r ) H (div)-approximation spaces.
Our numerical test agrees with that taken in [1], where results for the full and reduced AC spaces and the mapped BDM spaces appear. Results for our fully direct mixed spaces agree very closely with the results for the AC spaces, and these far exceed the performance of the mapped BDM spaces.

Summary and conclusions
We defined finite elements on a nondegenerate, convex quadrilateral E of minimal local dimension (i.e., with a minimal number of DoFs) using a different construction technique than what has been used by other authors. Our approach allows us to define entire families of finite elements, one for each index of approximation r . Moreover, the choices made in the construction give rise to different families of finite elements. To the best of our knowledge, they are all novel, except in three instances. The classic serendipity spaces on rectangles are recovered by our techniques as noted in Remark 1 (and then also the mixed elements that these give rise to, i.e., to the BDM family of mixed elements on rectangles). The two families of AC mixed finite elements are known [1], and these can be recovered by our construction as discussed in Sect. 7.1. It is also possible that for certain families, and then only for low order finite elements (with r ≤ 2), that our construction leads to known finite elements. It is difficult to determine what choices we would need to make in the construction of our families of finite elements that would reproduce these known elements, and we did not attempt to do that here.
We defined a wide variety of direct serendipity elements on a nondegenerate, convex quadrilateral E and proved their convergence properties (Theorems 1, 2, and 3 ). They have the form The supplemental space S DS r (E) has dimension 2 when r ≥ 2 and 1 when r = 1. It can be defined by a choice of four functions. Referring to Fig. 1, the linear functions λ 24 and λ 13 are arbitrary except that the zero line of λ 24 should intersect the edges e 2 and e 4 , and λ 13 should intersect the edges e 1 and e 3 . The (most likely nonlinear) functions R 13 and R 24 should be continuously differentiable and reduce to −1 (or any negative constants) on e 1 and e 2 , respectively, and 1 (or any positive constants) on e 3 and e 4 , respectively. For example, one can take the simple choices λ 24 = λ 2 − λ 4 , λ 13 The case r = 1 is special, but when r ≥ 2, the fully direct supplemental space is then A supplemental space can also be defined using the bilinear map betweenÊ = [−1, 1] 2 and E as S DS ,mapped We defined a wide variety of direct mixed elements on E and proved their convergence properties (Theorem 4 and Lemma 2). The de Rham theory is useful in this regard, and the elements take the form for reduced and full H (div)-approximation spaces, respectively, where If (127) is used, one obtains the AC spaces [1] (these are are not fully direct since the supplemental functions are mapped from the reference element). Otherwise, the choices we make in the construction lead to different families of finite elements, and they are the first families of fully direct mixed spaces with a minimal number of DoFs to be defined on quadrilaterals.
The direct serendipity and mixed elements can be merged to create H 1 (Ω) and H (div; Ω) conforming spaces, respectively, on a mesh T h of nondegenerate, convex quadrilaterals. If the meshes are shape regular as h → 0 and the four functions in (125) are chosen to be continuously differentiable with respect to the vertices of the element, then the spaces have both optimal approximation properties and minimal local dimension.
Numerical results were presented to illustrate their performance. All of the spaces attained the theoretical optimal convergence rates. On a fixed mesh, the direct serendipity spaces are more efficient in terms of accuracy than finite elements based on mapping tensor product polynomials to quadrilaterals. When comparing the total number of DoFs versus accuracy, it is not clear whether the finer mesh direct serendipity or coarser mesh tensor product spaces will be more efficient. In any case, the direct serendipity spaces are needed to define the direct mixed elements.
As a final result, we make a simple observation. The de Rham complex (64) in R 2 is sometimes described as where the curl of a vector function ψ ψ ψ = ψ 1 , ψ 2 is the scalar curl ψ ψ ψ( This form of the complex is just a rotation of (64). It gives us reduced and full direct H (curl)-approximating elements where These spaces can be arranged to form an exact sequence that discretizes the de Rham complex, namely for any variant of our new direct serendipity spaces. We can merge these elements globally by matching tangential components across the edges. We close by mentioning two open problems relating to this work. First, we plan in future work to extend the ideas presented in this paper to define families of direct serendipity and mixed finite elements on convex polygons [4]. We expect that a similar construction can be given, although the number of supplemental functions will depend on the number of sides N of the polygon E N and the index r . Our construction of serendipity elements for N = 4 was to define them for r ≥ 2 and then define the elements for r = 1 using the r = 2 elements. We expect that a similar development can be used for polygons.
A second open problem is to determine how to construct general families of direct finite elements on cuboidal hexahedra (and perhaps even on convex polytopes), perhaps by generalizing the work [3] on mixed elements. The cubical case shows that the supplemental functions of the serendipity elements are used to separate DoFs on both faces and edges. This mixing of concerns makes the construction on hexahedra a bit more involved. Moreover, for H 1 conformity, the traces of the elements on the faces cannot reduce to just polynomials. If the scalar serendipity spaces can be defined, one cannot discover the mixed finite elements immediately. The relevant de Rham sequence is and its discrete analogue is the exact sequence where the space of vector-valued functions X must be determined.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.