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

: New families of direct serendipity and direct mixed ﬁnite elements on general planar, strictly convex polygons were recently deﬁned by the authors. The ﬁnite elements of index r are H 1 and H ( div ) 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 deﬁned directly on the element and augmented with a space of supplemental functions. The supplemental functions were constructed as rational functions, which can be difﬁcult 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 ≥ 1. We prove the optimal approximation property for these new ﬁnite 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.

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 P r (E) 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 ∂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 ⊂ R 2 be a closed, nondegenerate, convex polygon with N ≥ 3 edges.The direct serendipity finite elements of index r ≥ 1 are H 1 -conforming and take the form DS r (E N ) = P r (E N ) ⊕ S DS r (E N ), where S DS r (E N ) is the space of supplemental functions.The direct mixed finite elements are H(div)-conforming and take two forms, for full (r ≥ 0) and reduced (r ≥ 1) H(div)-approximation, respectively, where Pr (E N ) 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 resulting in, for s = r − 1, r (s ≥ 0), The consequence is that and, therefore, 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 S DS r (E N ) 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 ∞ -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 ≥ 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.

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, . . ., N, and the vertices be x v,i = e i ∩ e i+1 .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, . . ., N.
e 4 e 5 e 1

Figure 1.
A pentagon E 5 , with edges e i , outer unit normals ν i , tangents τ i , and vertices x v,i .
For any two distinct points y 1 and y 2 , let L[y 1 , y 2 ] be the line passing through y 1 and y 2 , and take ν[y 1 , y 2 ] to be the unit vector normal to this line interpreted as going from y 1 to 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 x to L[y 1 , y 2 ] as To simplify the notation for linear functions that will be used throughout the paper, let ] be the line containing edge e i and let λ i (x) give the distance of x ∈ R 2 to edge e i opposite the normal direction, i.e., These functions are strictly positive in the interior of E N , and λ i vanishes on the edge e i .

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 S DS r (E N ) 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 S DS r (E 3 ) is empty.When N ≥ 4 and 1 ≤ r < N − 2, the direct serendipity spaces DS r (E N ) are defined as subspaces of DS N−2 (E N ) by the rule DS r (E N ) = ϕ ∈ DS N−2 (E N ) : ϕ| e ∈ P r (e) for all edges e of E N .
Therefore, we only need to understand S DS r (E N ) for r ≥ N − 2 and N ≥ 4.
To define the supplemental basis functions, two series of choices must be made for each i, j such that 1 ≤ i < j ≤ N and 2 ≤ j − i ≤ N − 2 (i.e., e i and e j are nonadjacent edges).First, as shown in Figure 2, one must choose two distinct points x 1 i,j ∈ L i and x 2 i,j ∈ L j that avoid the intersection point x i,j = L i ∩ L j , if it exists.Then let be the linear function associated to the line L i,j = L[x 1 i,j , x 2 i,j ].Second, R i,j must be chosen to satisfy The supplemental space for r ≥ N − 2 ≥ 2 is of the form

Direct Serendipity Finite Elements
Every shape function of the direct serendipity finite element DS r (E N ) 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 ψ ∈ DS r (E N ), one can take where dσ 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 DS r (E N ) corresponding to the DoFs can be constructed.Given a computational mesh of convex polygons T h over a domain Ω, the basis can be simply pieced together to form a global H 1 -conforming basis of the space DS r (Ω) ⊂ H 1 (Ω).

Direct Mixed Finite Elements
As discussed in the introduction, full, V r r (E N ), and reduced, V r−1 r (E N ), H(div)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 S DS r+1 (E N ) by the simple Formula (6).
The DoFs for these spaces (with s = r ≥ 0 or s = r − 1 ≥ 0) can be taken to be where the H 1 (E N ) and H(div; E N ) bubble functions, for r ≥ N − 1, are Given the mesh T h over Ω, one constructs the basis and the H(div)-conforming global space V s r (Ω) ⊂ H(div) (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.

Piecewise Continuous Supplements
In [43], R i,j satisfying ( 11) on E N , for 1 ≤ i < j ≤ N, 2 ≤ j − i ≤ N − 2, was taken to be the simple rational function 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 (E N ), and the second for constructing smoother supplemental functions in H p (E N ) for integer p ≥ 1.

Supplemental Functions in H 1 (E N )
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 T n (E N ).One picks a vertex x v,n and divides E N into N − 2 sub-triangles.The sub-triangles are T n m with vertices x v,n , x v,m , and x v,m+1 , where m = n + 1, . . ., n + N − 2. For the second sub-triangulation, depicted in Figure 4 and denoted as T x c (E N ), one picks a point x c in the interior of E N and divides it into N sub-triangles.Now the sub-triangles are T x c m with vertices x c , x v,m , and x v,m+1 , where m = 1, 2, . . ., N. We use the centroid of the element for x c .
A sub-triangulation based on a common fixed vertex.Shown is T 5 (E 5 ) using the fixed vertex x v,5 .
Let the piecewise polynomial function space of degree s corresponding to each subtriangulation be We construct R i,j in P 1 (T n (E N )) or P 1 (T x c (E N )), depending on which of the two subtriangulations is used, such that by using interpolation at the vertices of the sub-triangles.If the sub-triangulation is chosen to be T n (E N ), the restrictions (24) uniquely specify all the vertex values.However, if the triangulation is T x c (E N ), the center value is not determined, so we assign R i,j (x c ) = 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 ∈ H 1 (E N ).After constructing the supplemental functions in (13) with this R i,j , each φ s,i,j is in P r+1 (T n (E N )) or P r+1 (T x c (E N )), and therefore also in H 1 (E N ).

Supplemental Functions in H p (E N )
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 λ k (x) is the linear polynomial giving the (signed) distance to the line L k extending edge e k .When e i and e j are parallel, we simply define R i,j as the linear polynomial When e i and e j are not parallel, we first define a sub-partition of E N by adding a single extra line i,j through a point x i,j as depicted in Figure 5.The point x i,j is chosen so that it is closer to L j than the endpoints of e i , i.e., The line i,j passes through x i,j and is parallel to e j .This line divides N near e i and E i,j,0 N near e j , i.e., Let ν i,j = −ν j be the unit normal vector of i,j pointing into E i,j,1 N , and let τ i,j = τ j be a unit tangent vector.
A sub-division of E 6 using the line i,j .
We next construct the function ρ 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 where p ≥ 1 is an integer.The function is continuous, since λ j (x) = λ j (x i,j ) on i,j implies that ρ i,j | i,j = 1 in either case of the definition.Moreover, in the tangential direction, and, in the normal direction, which is continuous for p > 1, so ρ i,j ∈ C 1 (E N ).By iterating the argument, we have that ρ i,j ∈ C p−1 (E N ) and so also in Finally, after constructing both ρ i,j and ρ j,i , we define which is −1 on e i , 1 on e j .Moreover, R i,j ∈ H p (E N ).The supplemental functions in (13) constructed with this R i,j lie in H p (E N ).
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 x M e,i of the edges e i , i = 1, 2, . . ., N, and gives rise to the spaces denoted DS M r (E N ) and V M,s r (E N ).We compute the minimal distance of the midpoints to the edges, i.e., Then for any two non parallel and nonadjacent edges e i and e j , simply take the partition line i,j to be the line parallel to e j that is the fixed distance h 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 x e,i,k for k = 1, 2, i = 1, 2, . . ., N. In this case, we simply take x i,j to be the closest of these points to L j , omitting x e,j,1 and x e,j,2 .We denote the resulting spaces DS T r (E N ) and V T,s r (E N ).
x v,5 x v,1 x v,2 x v,3 x v,4 x M e,4 h M x v,5 x v,1 x v,2 x v,3 x v,4 x e,5,2 x e,3,1 x e,4,1 x e,3,2 x e,4,2 Figure 6.The two sub-partitions of E 5 used for constructing specific H p supplemental functions.The left one is used for DS M r and V M,s r , where h M = λ 5 (x M e,4 ).The right one is used for DS T r and V T,s r , where the closest trisection points are marked in red.

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 T h of convex polygons partitioning a domain Ω, 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 ∈ T h , let h E N be its diameter.Denote by T i , i = 1, 2, . . ., 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 {diameter of the largest circle inscribed in The shape regularity parameter of the single mesh T h is Assumption 1.The collection of meshes {T h } h>0 is uniformly shape regular.That is, the shape regularity parameters are bounded below by a positive constant: there exists σ * > 0, independent of T h and h > 0, such that the ratio We also require some mild restrictions on the construction of S DS r (E N ).
Assumption 2. For every E N ∈ T h , assume that the functions of S DS r (E N ) are constructed using λ i,j such that the zero set L i,j intersects e i and e j .Moreover, suppose that R i,j ∈ H p (E N ) for some p ≥ 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 x c in Section 4.1 (say as the centroid) and x i,j satisfying (26) in Section 4.2 (say be taking x i,j as the closer endpoint of e i to L j , or so that λ j (x i,j ) = We state first the approximation result for DS r (Ω).
Theorem 1.Let r ≥ 1, 1 ≤ p ≤ ∞, and > 1/p (or ≥ 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 ∈ W ,p (Ω), 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 T h , we construct an interpolation operator I r h : W l p (Ω) → 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, . . ., dim 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 ϕ 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 ⊂ ∂Ω if a i ∈ ∂Ω.We use these to define the dual basis ψ i with respect to L 2 (K i ), i = 1, 2, . . ., dim DS r .The corresponding interpolation operator I r h : W l p (Ω) → DS r is then There are two essential steps towards showing the approximation property.First, the nodal basis functions are bounded, and, second, the dual basis functions are bounded up to a scaling factor, 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) ÊN .We can then define a map F E N : ÊN → E N as a composition of a map that changes the geometry but not the size to ẼN , and then a scaling map (see [43] for precise details).
(0, 0) Ê5 (1, 0) x x v,3 x v,4 x v,5 x Figure 7.An element E 5 ∈ T h is shown on the right-hand side in its translated and rotated local coordinates.It is the image of a regular reference polygon Ê5 on the left-hand side.The map is decomposed into one that changes the geometry but not the size F Ẽ5 : Ê5 → Ẽ5 , and a scaling map x → H x.
Define the nodal basis functions ϕ ẼN i on ẼN .It is enough to show the boundedness of their W m q norms.Although they are no longer smooth functions, compared to [42,43], they are continuous on Ẽ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 ẼN .Therefore, ϕ ẼN i will still depend continuously on x = F −1 ẼN ( x) and the vertices of ẼN , which vary in a compact set.We conclude that the nodal basis functions are bounded in W m q norm.The boundedness of ψ i in the L ∞ norm can be shown in a similar way.
For the mixed finite elements, we have the following result, wherein we see projection operators π : H(div; Ω) ∩ (L 2+ (Ω)) 2 → V s r , s = r − 1, r, where > 0, and P W s , the L 2 -orthogonal projection operator onto 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 h > 0.
For the proof, we define the projection operator π by piecing together local operators π 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.

Numerical Results
We present numerical experiments for our new finite elements as applied to Poisson's equation p = 0 on ∂Ω, (46) where f ∈ L 2 (Ω).The corresponding weak form finds p ∈ H 1 0 (Ω) such that where (•, •) is the L 2 (Ω) inner product.Setting we have the mixed weak form, which finds u ∈ H(div; Ω) 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].
We perform our tests on a unit square domain Ω = [0, 1] 2 , and take the source term f (x) = 2π 2 sin(πx 1 ) sin(πx 2 ), so the exact solution is u(x 1 , x 2 ) = sin(πx 1 ) sin(πx 2 ).We consider five types of supplemental spaces.The original direct serendipity and mixed finite element spaces will be denoted DS R r and V R,s r , 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 DS V r and its corresponding mixed spaces as V V,s r , and those based on the center point sub-triangulation as DS C r and V C,s r , respectively.The spaces based on the H p supplements were described in Section 4.2 and denoted DS M r , V M,s r and DS T r , V T,s r .

The Meshes Used
Approximate solutions are computed on a sequence of Voronoi meshes T 2 h 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 T 2 h 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.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.

Results for Direct Serendipity Spaces
We present and compare the results of the numerical tests performed for DS R r , DS V r , DS C r , DS M r , and DS T r , where r = 2, 3, 4, 5.We take p = 1 in Section 4.2 for the construction of DS M r and DS T r , 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 Tables 2 and 3 suggests, the convergence rates at n = 10, 14, 22 for DS R r are all slightly better than optimal in this test.However, we can observe a slower convergence rate at n = 18 for DS R r .(We interject that convergence rates are computed as improvement from the previous mesh in our sequence.)In Table 4, we compare the results for the n = 18 mesh of DS R r , DS V r , DS C r , DS M r , and DS T r .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 DS R 5 .On the other hand, the errors for r = 2, 3, 4 of the new spaces are slightly worse than those of DS R r , and among all the new spaces, DS C r has the best performance in error.We conclude that DS C r 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 DS M r and DS T r , of which R i,j for e i ∦ e j are flat in the middle and oscillate near the boundary.In contrast, the supplements from DS V r and DS C r are reasonably shaped, and those from DS C r are better since its partition has sub-triangles that are more shape regular (as was shown in Figures 3 and 4).This argument is also supported by the observation that the results are usually worse if we take larger p for DS M r and DS T r , where the shape of the supplements are even worse.We perform numerical tests for V R,s r , V V,s r , V C,s r , V M,s r , V T,s r , for the full H(div)approximation spaces where r = s = 0, 1, 2, 3, and the reduced H(div)-approximation spaces where r = 1, 2, 3, and s = r − 1.Since those mixed spaces are constructed from corresponding direct serendipity spaces 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 Tables 5 and 6.All spaces perform similarly well, although V R,s r usually performs best in these tests.Among the new spaces, V C,s r performs a bit better, and it gives results close to those of V R,s r .For reference, we provide the numerical results for V C,s r in Tables 7 and 8.
Table 5. Errors and convergence rates at n = 18 computed from the previous step n = 14, for the reduced , and V

Conclusions
We reviewed the construction of direct serendipity and mixed finite elements on nondegenerate, planar convex polygons.The direct serendipity finite element spaces are of the form DS r (E N ) = P r (E N ) ⊕ S DS r (E N ).
The full and reduced H(div)-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 where S V r (E N ) = curl S DS r+1 (E N ).
We presented two approaches to construct the supplemental functions in S DS r (E N ) as piecewise polynomials.The first approach divides a polygonal element E N into subtriangles, and constructs the supplements as continuous piecewise polynomials that lie in H 1 (E N ).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 (E N ) for any p ≥ 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.

Table 1 .
Shape regularity parameters for each mesh T 2 h .

Table 2 .
L 2 errors and convergence rates for DS R r and DS C r .

Table 3 .
H 1 -seminorm errors and convergence rates for DS R r and DS C r .

Table 6 .
Errors and convergence rates at n = 18 computed from the previous step n = 14, for the full H

Table 7 .
Errors and convergence rates in L 2 for V C,r−1 r r , and V T,r r .p− p h u − u h ∇ • (u − u h )

Table 8 .
Errors and convergence rates in L 2 for V C,r r .
p − p h u − u h ∇ • (u − u h )