A family of C1 quadrilateral finite elements

We present a novel family of C1 quadrilateral finite elements, which define global C1 spaces over a general quadrilateral mesh with vertices of arbitrary valency. The elements extend the construction by Brenner and Sung (J. Sci. Comput. 22(1-3), 83-118, 2005), which is based on polynomial elements of tensor-product degree p ≥ 6, to all degrees p ≥ 3. The proposed C1 quadrilateral is based upon the construction of multi-patch C1 isogeometric spaces developed in Kapl et al. (Comput. Aided Geometr. Des. 69, 55–75 2019). The quadrilateral elements possess similar degrees of freedom as the classical Argyris triangles, developed in Argyris et al. (Aeronaut. J. 72(692), 701–709 1968). Just as for the Argyris triangle, we additionally impose C2 continuity at the vertices. In contrast to Kapl et al. (Comput. Aided Geometr. Des. 69, 55–75 2019), in this paper, we concentrate on quadrilateral finite elements, which significantly simplifies the construction. We present macro-element constructions, extending the elements in Brenner and Sung (J. Sci. Comput. 22(1–3), 83–118 2005), for polynomial degrees p = 3 and p = 4 by employing a splitting into 3 × 3 or 2 × 2 polynomial pieces, respectively. We moreover provide approximation error bounds in L∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$L^{\infty }$\end{document}, L2, H1 and H2 for the piecewise-polynomial macro-element constructions of degree p ∈{3,4} and polynomial elements of degree p ≥ 5. Since the elements locally reproduce polynomials of total degree p, the approximation orders are optimal with respect to the mesh size. Note that the proposed construction combines the possibility for spline refinement (equivalent to a regular splitting of quadrilateral finite elements) as in Kapl et al. (Comput. Aided Geometr. Des. 69, 55–75 30) with the purely local description of the finite element space and basis as in Brenner and Sung (J. Sci. Comput. 22(1–3), 83–118 2005). In addition, we describe the construction of a simple, local basis and give for p ∈{3,4,5} explicit formulas for the Bézier or B-spline coefficients of the basis functions. Numerical experiments by solving the biharmonic equation demonstrate the potential of the proposed C1 quadrilateral finite element for the numerical analysis of fourth order problems, also indicating that (for p = 5) the proposed element performs comparable or in general even better than the Argyris triangle with respect to the number of degrees of freedom.


Introduction
Since the origin of the finite element method (FEM), there has been a significant interest for the construction of C 1 spaces. The main motivation comes from the discretization of fourth-order problems, where a direct Galerkin formulation requires C 1 global smoothness of the approximating functions. Indeed, the first constructions have been plate elements: the Argyris triangle [1] and the Bell triangle [4]. Both elements require polynomials of degree p ≥ 5, and are additionally C 2 at the vertices. While the normal derivative along an edge is of degree p − 1 for the Argyris element, its degree reduces to p − 2 for the Bell element. This leads for instance in case of polynomial degree p = 5 to the fact that the Argyris triangular space possesses six degrees of freedom for each vertex and one degree of freedom for each edge, while the Bell triangular space just has six degrees of freedom for each vertex and no additional degrees of freedom for the edges. For more details on the Argyris and Bell triangular element as well as on other C 1 triangular finite elements, we refer to the books [8,13]. C 1 finite element spaces of lower polynomial degree are in general based on splines, which are constructed over general triangulations, see [37].
The design of C 1 finite elements over quadrilateral meshes is in general more challenging compared to the case of triangular meshes, in particular with respect to the selection of the degrees of freedom. It is of course always possible to split a quadrilateral into several triangles and employ constructions over triangulations. One such possibility is the de Veubeke-Sander element as developed in [16,48], which is a piecewise polynomial element of degree p = 3 over a quadrilateral split into four triangles by splitting along the diagonals. However, C 1 constructions for "true" quadrilateral elements, i.e., elements having local spaces defined by bivariate polynomials of maximum degree p mapped onto a quadrilateral using a bilinear mapping, are significantly more complicated. Only quite recently, such constructions were developed in [5,9,38]. For certain special cases, constructions are well-known, such as the Bogner-Fox-Schmit element [7], which is a simple bivariate Hermitetype C 1 construction that works for low polynomial degrees such as p = 3, but is limited to tensor-product meshes. In contrast, the C 1 elements [5,9,38] are applicable to more general quadrilateral meshes, but require a polynomial degree p ≥ 6 in case of [9] and a polynomial degree p ≥ 5 (for some specific settings just p = 4) in case of [5,38]. The degrees of freedom for the finite element space [9] are selected similar to the Argyris triangular finite element space [1] by enforcing additionally C 2 -continuity at the vertices.
In contrast, the functions in [5,38] are just C 1 at the vertices and the degrees of freedom are defined by means of the concept of minimal determining sets (cf. [37]), which is a common strategy for the construction of C 1 splines over triangular meshes, see also [37]. A different but related problem is the construction of C 1 function spaces over general quadrilateral meshes for the design of surfaces, such as in [21,43,44,47]. The methods are based on the concept of geometric continuity [45], which is a well-known tool for generating smooth complex surfaces in computer aided geometric design.
An alternative to FEM is the use of isogeometric analysis (IgA), which was introduced in [23], and employs the same spline function space for describing the physical domain of interest and for representing the solution of the considered PDE, see, e.g., [15,23] for more details. In case of a single patch geometry, this allows the direct discretization of fourth order PDEs [52], such as Kirchhoff-Love shells, e.g., [35,36], the Navier-Stokes-Korteweg equation, e.g., [20], problems of strain gradient elasticity, e.g., [18,42], or the Cahn-Hilliard equation, e.g., [19], by just using C 1 splines. In case of multi-patch geometries with possibly extraordinary vertices, i.e., vertices with a patch valency different to four, the design of smooth spline spaces is challenging and is the topic of current research.
Depending on the used type of parametrizations for the single patches of the given unstructured quadrilateral mesh, different techniques for the design of a C 1 spline space over this mesh have been developed. Possible examples in the case of planar, unstructured quadrilateral meshes are to use C 1 multi-patch parametrizations with a singularity at an extraordinary vertex, e.g., [41,53], multi-patch parametrizations which are C 1 except in the vicinity of an extraordinary vertex, e.g., [32][33][34]40], or multi-patch parametrizations which have to be just C 0 at all interfaces, e.g., [6, 11, 12, 14, 26-28, 30, 31, 39]. For more details about existing C 1 constructions for unstructured quadrilateral meshes, we refer to the recent survey articles [29] and [24]. Beside this, in [10,49], different approaches for the construction of smooth spline functions of degree p are presented, which are C s (1 ≤ s ≤ p − 1) everywhere, except in the vicinity of an extraordinary vertex, where they are just C 0 .
In this work, we present a family of C 1 quadrilateral finite elements, that are the low-degree (for p ∈ {3, 4, 5}) and spline counterparts of the quadrilateral finite elements proposed in [9] (for p ≥ 6) by Brenner and Sung. The interest for the low-degree case is that the computational cost for the linear system formation (due to numerical quadrature) and solution (that depends on the matrix conditioning) is more favorable. The spline constructions of arbitrary degree p ≥ 3 allow for simple, more efficient refinement rules that yield nested spaces. Such constructions based on spline refinement may be utilized to obtain efficient multigrid solvers as in [51].
The proposed quadrilateral (macro-)elements possess similar degrees of freedom as the classical C 1 Argyris triangle [1]. An advantage of the quadrilateral construction over the triangular one is the simpler extension to the lower polynomial degrees p = 3 and p = 4 by just using tensor-product splines without the need of special splits for the mesh elements.
While in [30] the optimal approximation properties of the C 1 isogeometric spline space is just numerically shown, in this work, the optimal approximation order of the C 1 quadrilateral space is proven. The isogeometric C 1 multi-patch spaces from [30] depend, along each interface, on the two neighboring patch parametrizations. Hence, basis functions restricted to one patch may also depend on the parametrizations of the neighboring patches. In the most general setting a purely local construction is not possible. However, when restricting to bilinearly mapped patches/elements, the construction becomes local to each element. Moreover, we explicitly give, for some particular cases, the Bézier or spline coefficients of the basis functions by simple formulas depending only on the vertices of the quadrilateral. Several numerical tests of solving the biharmonic equation also show the potential of the C 1 quadrilateral spaces for the numerical analysis of fourth order PDEs.
The outline of this paper is as follows. Section 2 introduces the quadrilateral mesh which will be used throughout the paper. In Sections 3-4, the construction of the C 1 quadrilateral is described, focusing first on the bi-quintic polynomials, and then generalizing to splines, which allow the use of the lower polynomial degrees p = 3 and p = 4. Section 3 also discusses the connection of the C 1 quadrilateral with two wellknown triangular finite elements, namely with the Argyris triangle [1] and with the Bell triangle [4]. Section 5 describes the global spaces constructed from the previous quadrilateral elements. In Section 6, we analyze the approximation properties of the C 1 quadrilateral space. Section 7 discusses the spline refinement of the quadrilateral elements introduced in Sections 3-4. Then, Section 8 describes the design of local basis functions of the C 1 quadrilateral space for the case of polynomials and its extension for the case of splines, respectively, giving for the low-degree p ∈ {3, 4, 5} cases the explicit Bézier and spline coefficients of the basis functions. In Section 9, numerical benchmarks on the biharmonic equation with different quadrilateral meshes are presented. Finally, the isoparametric extension of the C 1 quadrilateral, and its relation to the isogeometric element of [30], is briefly discussed in Section 10. We conclude the paper in Section 11.

Quadrilateral mesh
We consider planar domains that allow meshing by quadrilaterals. Note that a generalization to domains with curved boundaries is possible with some additional care. We refer the reader to [5,28], where such discretizations were developed, see also Section 10.
Let Ω ⊂ R 2 be an open, planar and connected region, with Lipschitz boundary. We assume that Ω allows a quadrangulation as defined below. This is the case if the boundary is piecewise linear, including all inner boundaries if Ω is not simply connected. The coordinates in physical space are given as ( consisting of a set of elements Q, edges E and vertices V satisfying the following properties.
-Each vertex is a point in the plane, that is V ⊂ R 2 .
-Each edge ε ∈ E is an open segment, and there exist two vertices v, v ∈ V such that -Each element Q ∈ Q is a convex, non-empty, open quadrilateral; there exist four vertices v (1) , . . . , v (4) ∈ V, four edges ε (1) , . . . , ε (4) ∈ E, with edge ε (k) connecting v (k) with v (k+1) (modulo 4), given in counter-clockwise order; and the element admits a parametrization by a mapping F Q : Q → Q which is bilinear on Q = [0, 1] 2 , precisely: Fig. 1 The last condition means that there are no hanging vertices in the quadrilateral mesh, i.e., all neighboring quadrilaterals share an entire edge or a vertex in their closure. An edge ε ∈ E of the mesh can be assigned either to the interior or to the boundary of the domain, if ε ⊂ Ω or ε ⊂ ∂Ω, respectively. Each inner edge is shared by exactly two elements, i.e., there exist exactly two Q, Q ∈ Q, such that ε = ∂Q ∩ ∂Q , whereas each boundary edge belongs to one element only, i.e., there exists exactly one Q ∈ Q such that ε ⊂ ∂Q ∩ ∂Ω. Moreover, we introduce the following notation. Fig. 1 we have, e.g.,v (1)

Definition 1 (Pre-images of points and edges) For every point
Given a Q ∈ Q, we introduce the following notation: We denote by Fig. 1 Visualization of the mapping F Q for a quadrilateral Q, with vertices, edges, parameter domain and local coordinates the vector corresponding to the edge ε (k) , and define with indices modulo 4. We have that a (k) = J(v (k) ), where Furthermore, h ε (k) = t (k) denotes the length of the corresponding edge ε (k) , and ρ Q its minimum angle defined as follows: if Q is a degenerate quadrilateral (that is a triangle), then ρ Q = 0; otherwise, ρ Q is the minimum of the angles of the four triangles that are formed by the edges and the diagonals of Q. We assume however that each Q ∈ Q is non-degenerate, indeed the following holds.

2)
and where c 1 (ρ Q ), c 2 (ρ Q ) ∈ R are constants that depend only on ρ Q and are strictly positive for ρ Q > 0.
Since J is a linear polynomial, its extrema are attained at the vertices of Q, that is min(J) = min{a (k) , k = 1, . . . , 4}. Since the lower bound (2.4) follows from Eqs. 2.2 and 2.3.
The condition ρ Q > 0 also implies that the parametrization F Q : Q → Q is regular, that is its inverse F −1 Q : Q → Q has bounded derivatives too. This follows from Eq. 2.4. Similar to [9], we assume that the quadrilateral mesh M is shape regular, that is

C 1 Brenner-Sung quadrilateral element
In the following, we recall the definition of the Brenner-Sung quadrilateral element from [9, Section 6], extend it to the lower degree cases, and define the associated piecewise polynomial C 1 space over the domain of interest Ω, given a quadrilateral mesh M. In our presentation, we loosely follow the style of [8,13]. The Brenner-Sung quadrilateral of degree p ≥ 5 is constructed from bi-quintic polynomials with normal derivatives across interfaces that are polynomials of degree p − 1. For p = 5, the degrees of freedom are given as C 2 -data at the vertices, normal derivatives at the edge midpoints, as well as interior point evaluations. This is in accordance with the degrees of freedom of the Argyris triangle, see [1]. We denote with Q p the space of bivariate polynomials of maximum degree p and with P p the space of polynomials of total degree p, either uni-or bivariate, depending on context.
In the next subsections, we introduce the local spaces and degrees of freedom corresponding to a single quadrilateral Q. One set of degrees of freedom is the normal derivative at the edge midpoint, where we use the following notation. For every edge ε ∈ E between vertices v and v , let m ε = 1 2 v + 1 2 v be the edge midpoint. Moreover, let n ε be its unit normal vector and ∂ n ε be the normal derivative of a function defined on Ω across the edge ε. Here, we assume that the direction of the normal is fixed for every edge of the mesh M.

Local space and degrees of freedom, p = 5
Given a quadrilateral Q ∈ Q we define the local function space and the local degrees of freedom as follows.
Definition 2 (Brenner-Sung quadrilateral for p = 5) Given a quadrilateral Q with vertices v (1) , v (2) , v (3) and v (4) and edges ε (1) , ε (2) , ε (3) and ε (4) follow-ing [13], we define the Brenner-Sung quadrilateral of degree p = 5 as (Q, P 5 Q , Λ 5 Q ), with (3.1) and The set of face points is given as See Fig. 2 for a visualization of the local degrees of freedom of the Brenner-Sung quadrilateral. Representations of first and second derivatives of ϕ in terms of local coordinates on the reference element can be found in Appendix A.
The unisolvency of the degrees of freedom Λ 5 Q for the space P 5 Q follows from the basis construction in Section 8.
We have already observed the similarity of this construction with the Argyris triangle. Let us recall the definition of the Argyris triangle for p = 5 as given in [1,13]. Definition 3 (Argyris triangle for p = 5) Given a triangle T with vertices v (1) , v (2) and v (3) and edges ε (1) , ε (2) and ε (3) , we define the Argyris triangle as (T , P A T , Λ A T ), with P A T = P 5 and Λ A T = Λ 0,T ∪ Λ 1,T , with Hence, the degrees of freedom for the Brenner-Sung quadrilateral (p = 5) and Argyris triangle (p = 5) are the same, except for the additional point evaluations at face points in the quadrilateral case. In addition, the traces as well as normal derivatives along edges are the same in both elements, i.e., for p = 5 traces are quintic polynomials and normal derivatives are quartic polynomials. The degrees of freedom for the Argyris triangle are visualized in Fig. 3 (left).
In addition, the condition that the normal derivative along an edge is of degree 4, is similar to the condition on the Bell triangular element [4], a quintic element, where normal derivatives are assumed to be polynomials of degree 3, thus eliminating the normal derivative degrees of freedom and resulting in 18 degrees of freedom per triangle.
Definition 4 (Bell triangle for p = 5) Given a triangle T with vertices v (1) , v (2) and v (3) and edges ε 1 , ε 2 and ε 3 , we define the Bell triangle as (T , P B T , Λ B T ), with and The degrees of freedom for the Bell triangle are visualized in Fig. 3 (right). Both triangle elements possess variants of higher degree, see [1,4,13,37]. For triangular elements, constructions of smooth spaces for lower degrees are usually based on special splits, such as the Clough-Tocher or Powell-Sabin 6-or 12-splits. Unlike the triangular case, in the quadrilateral case, variants of lower degree are relatively straightforward: they are based on regular splits of the quadrilaterals as visualized in Fig. 4 and follow from the spline constructions developed in [30].

Local space and degrees of freedom, p ≥ 6
Given a quadrilateral Q ∈ Q we define the local function space and the local degrees of freedom for p ≥ 6 as follows.
Definition 5 (Brenner-Sung quadrilateral [9]) Given a quadrilateral Q with vertices v (1) , v (2) , v (3) and v (4) and edges ε (1) , ε (2) , ε (3) and ε (4) we define the Brenner-Sung Here, F ε (k) = F Q |ε(k), and the set of face points is given as As one can easily see, Definition 5 covers also the case of Definition 2. Obviously, we have the following. The degrees of freedom Λ p Q are unisolvent for the space P p Q . Indeed, one can show (see [9]) that the dimension of P p Q is given by dim(Q p ) = (p + 1) 2 minus the number of constraints from (∂ n ε (k) ϕ • F Q )|ε(k) ∈ P p−1 , which are one per edge, that is four. Then the dimension of P p Q is (p + 1) 2 − 4 and equals the cardinality of Λ p Q .

Quadrilateral spline macro-element
In the following we extend the construction of Section 3 on quadrilaterals to lower degrees p = 3 and p = 4 using a split into sub-elements, that is we assume that the parameter domain Q is split into sub-elementsq ∈ s κ ( Q), with In this section, we focus on the coarsest sub-element split, that is κ = 3 and κ = 2 for p = 3 and p = 4, respectively, while the general case will be adressed in Section 7.
Definition 6 (C 1 quadrilateral macro-element for p ∈ {3, 4}) Given a quadrilateral Q with vertices v (1) , v (2) , v (3) and v (4) and edges ε 1 , ε 2 , ε 3 and ε 4 we define the For p = 4, the set of face points is given as As for p = 5, the degrees of freedom Λ p Q completely determine the functions from the space P p Q and the dimension is given by dim(P This follows as a special case of Lemma 7.
In Fig. 4, we visualize the polynomial sub-elements from Eq. 4.2 and local degrees of freedom from Eq. 4.3 for p ∈ {3, 4}.

Global space and degrees of freedom
In this section, we describe the global space and set of degrees of freedom from the local spaces and degrees of freedom defined above, with focus on the low-degree cases p ∈ {3, 4, 5}.
Definition 7 (Global degrees of freedom) Let p ∈ {3, 4, 5}. Given a quadrilateral mesh M, we have the degrees of freedom Λ p , given as The global degrees of freedom in Definition 7 together with the finite element descriptions in Definitions 2 and 6 determine a global space S p (M) ⊂ C 1 (Ω).
and let M of Ω and let the space S p (M) be given by the degrees of freedom Λ p as in Definition 7, with where the local spaces P p Q are given as in Definition 2 or 6, respectively. Then, the global space satisfies S p (M) ⊂ C 1 (Ω) and we have Proof Note that the piecewise polynomial space P p Q from Definition 6 covers also the polynomial case p = 5 for κ = 1, where Consequently, since F Q |ε is a linear function, we have ϕ| ε∩q ∈ P p , ϕ| ε ∈ C p−1 , Hence, ϕ| ε is a piecewise polynomial of degree p. Thus, ϕ| ε has 6 degrees of freedom, which follows from Eq. 5.1. Value, first and second derivative (in direction of the edge) of ϕ at the two vertices of ε are determined by the C 2 -data. The function ϕ| ε is thus completely determined by the C 2 -data. This is independent of the element Q, Q under consideration. Hence, we have ϕ ∈ C 0 (Ω). Moreover, due to Eq. 5.2, the function ∂ n ε ϕ| ε is a piecewise polynomial of degree p − 1, having 5 degrees of freedom, independent of Q, Q . Of those 5 degrees of freedom, the C 2 -data at the vertices determine two each, whereas one is determined by ∂ n ε ϕ(m ε ). Hence, ϕ| ε and ∂ n ε ϕ| ε are completely determined by the global degrees of freedom and ϕ ∈ C 1 (Ω). What remains to be shown is that dim(S p (M)) = |Λ p |. Its proof follows directly from a simple counting argument.
We have presented Lemma 1 and its proof purely in terms of a finite element setting, considering the local spaces and global degrees of freedom. See [30,Section 4] for a more general statement on spline patches. Note that the space S p (M) is C 2 at all vertices by construction.

Remark 1
Since both the degrees of freedom Λ p Q as well as the definition of the local space P p Q depend on derivatives in normal direction, the proposed Brenner-Sung quadrilaterals (including the macro-element variants) are not affine invariant, as the Argyris triangle, which possesses an affine invariant space, but no affine invariant degrees of freedom.

Approximation properties
In this section, we prove local and global approximation estimates, where the error is measured only in the norms of interest · L ∞ , · L 2 , and · H , with ∈ {1, 2}, for simplicity. For the notation concerning Sobolev spaces, we follow [8].
Given a convex quadrilateral Q ∈ Q, the main ingredient to prove the local approximation estimate is the projector Π P p The existence of such a basis is a consequence of the unisolvence of the set of degrees of freedom Λ p Q . A key property for the approximation result is the basis stability stated below.
where · is any of the norms of interest.
Proof Each basis function β λ Q can be obtained by imposing the conditions to belong to the space and to be in duality to the degrees of freedom In parametric coordinates, it means that β λ Q • F Q defined on Q is a polynomial (for p ≥ 5, it belongs to Q p ) or piecewise polynomial (for p ∈ {3, 4}, its restriction to each subelementq ∈ s 6−p ( Q) belongs to Q p ) that fulfills (6.2)-(6.3). These conditions above involve the first and second derivatives of the inverse parametrization F −1 Q that are well defined and bounded on Q thanks to Eq. 2.4. Recalling the expression (2.5) of ∇F Q , the first and second derivatives of F −1 Q are rational polynomials in x 1 and x 2 and depend continuously on the parameters t (1) , . . . , t (4) . Therefore, β λ Q only depends on t (1) , . . . , t (4) , the dependence is continuous and the parameters belong to the compact set t (1) , . . . , t (4) thanks to Proposition 1. Continuity and compactness give the existence of a maximum of β λ Q which only depends on h Q , ρ Q and on p.
The bound from Lemma 2 yields the following local stability of the projector.

Lemma 3
Let Q ∈ Q be any convex quadrilateral. There exists a constant C > 0, dependent on h Q , ρ Q , and p such that for all ψ ∈ C 2 (Q), where · is one of the norms of interest and is the supremum of all derivatives up to second order on Q.
Proof Thanks to Lemma 2, we have that and then we use the obvious continuity |λ(ψ)| ≤ ψ C 2 (Q) .
The next two Lemmata, from [8], concern standard Sobolev inequalities and standard polynomial approximation over Q ∈ Q.

Lemma 4 [8, Lemma 4.3.4]
Let Q ∈ Q be any convex quadrilateral. There exists a constant C SI > 0, dependent on h Q and ρ Q such that for all ψ ∈ H 4 (Q), we have ψ ∈ C 2 (Q) and where Π P p ϕ is the averaged Taylor polynomial of degree p of ϕ over B.
The last property we need is that the quadrilateral element space under consideration contains the polynomials of total degree p.
Proof Let ψ ∈ P p , we need to show that ψ •F Q ∈ Q p and (∂ n ε (k) ψ •F Q )|ε(k) ∈ P p−1 for all ε (k) , according to Eqs. 3.1 and 4.2. Note that we do not need to consider the sub-elements separately, as ψ is a global polynomial. The composition of a polynomial of total degree p with a bilinear function always results in a polynomial of maximum degree p; hence, we have ψ • F Q ∈ Q p . Moreover, the directional derivative ∂ n ε (k) ψ is a polynomial of total degree p − 1, restricted to an edge yields a univariate polynomial of degree p − 1, which gives (∂ n ε (k) ψ • F Q )|ε(k) ∈ P p−1 since F Q |ε(k) is a linear parametrization.
We can now state and prove the local approximation estimate.
Theorem 1 Let Q ∈ Q be a convex quadrilateral. There exists a constant C > 0, dependent on ρ Q and on p, such that for 0 ≤ ≤ 2, 4 ≤ m ≤ p + 1 and for all ϕ ∈ H m (Q) we have Proof The proof follows the proof of [8,Theorem 4.4.4]. We can assume h Q = 1, since the general case and the role of h Q in the estimates follow by an homogeneity argument. We have Applying the bounds from Lemmata 3, 4 and 5, we obtain The L ∞ -estimate follows the same idea as the H -estimates, where a bound of the form is needed together with estimates similar to Lemma 4 and 5. Note that in case of the L ∞ estimate we only need m ≥ 3, see again [8,Theorem 4.4.4].
From this local error estimate, a global estimate follows straightforwardly. Let be the global projector defined as where each β λ ∈ S p (M) satisfies λ(β λ ) = 1 and λ (β λ ) = 0 for all λ = λ ∈ Λ p . By definition of the local and global spaces and degrees of freedom, the global projector and the local projector fulfill, for any Q ∈ Q, For a given local functional λ Q (·) = λ(·| Q ), we have β λ Q = β λ | Q . Hence, the support of β λ is given by all elements on which λ is defined, i.e., one element for all face point evaluations, two neighboring elements for all edge midpoint evaluations and, in case of vertex degrees of freedom, all elements around the vertex.

Quadrilateral macro-element: spline refinement
We can extend the definitions of Sections 3 and 4 to certain B-spline based macroelements of any degree p ≥ 3. For simplicity, and to highlight especially the low degree construction, we introduce the details only for p = 3. The higher degree cases can be constructed analogously. For B-spline macro-elements the degrees of freedom are given as C 2 -data at the vertices, normal derivative and point data at certain points along the edges, as well as suitably many interior functions that have vanishing values and gradients at all element boundaries. In such a setting, refinement can be performed either by splitting the macro-elements or by knot insertion within every macro-element. Note that, in the construction below, the continuity within the macro-element is of order 1. For p ≥ 3 this can be generalized to any order of continuity r, with 1 ≤ r ≤ p − 2.
We assume that every quadrilateral Q is split into κ × κ elements by mapping a regular split of the parameter domain Q = [0, 1] 2 using F Q . Let S p,r κ be the univariate B-spline space of degree p and regularity r over the interval [0, 1] split into κ polynomial segments of the same length, i.e., having the knot vector 0, . . . , 0, . . , κ + p − 1, be the corresponding Greville abscissae for the first and second knot vector, respectively. Note that the Greville abscissae γ i corresponding to a given knot vector (ξ 0 , . . . , ξ N ) are defined as knot averages γ i = (ξ i+1 + · · · + ξ i+p )/p for i = 0, . . . , N − p − 1.
As in Definition 2, we can define the local function space and the local degrees of freedom, where we need to assume κ ≥ 3 in order to be able to split the vertex degrees of freedom.
Definition 8 (C 1 quadrilateral spline macro-element for p = 3) Given a quadrilateral Q with vertices v (1) , v (2) , v (3) and v (4) we define the cubic C 1 quadrilateral spline macro-element as (Q, P Q,κ , Λ Q,κ ), with Here, r (k) , and the set of face points is given as We have the following lemma, which is related to the results developed in [30, Section 4].
The final step is then to show that, given an arbitrary ϕ ∈ P Q,κ , if every functional in Λ Q,κ vanishes on ϕ then ϕ = 0. If all derivatives up to second order of ϕ at the corners vanish, then the trace of ϕ is set to zero due to ϕ(r (k) j 0 ) = 0, for 1 ≤ j 0 ≤ κ − 3. This is an interpolation property for S 3,2 κ with zero end values, first and second derivatives and the interpolation is well-defined since r (k) j 0 are suitable Greville points. Analogously, the normal derivative of ϕ has to be zero since ∂ n ε (k) ϕ(q (k) j 1 ) = 0, for 1 ≤ j 1 ≤ κ − 2 (this is an interpolation property for S 2,1 κ with zero end values and first derivatives). Once ϕ and ∇ϕ are null on ∂Q, all functionals in Λ * 2,Q,κ evaluating to zero finally implies ϕ = 0.
It follows immediately from Lemma 7 that the dimension of the space can be determined completely by counting, having six degrees of freedom per vertex, 2κ − 5 degrees of freedom per edge and (2κ − 2) 2 degrees of freedom inside the element. Thus, we need 2κ − 5 ≥ 0, yielding the constraint κ ≥ 3.
Similar to Lemma 1 defining the global space S 3 (M), we can define a global space S 3 κ (M) using the spline macro-elements as given in Definition 8. Note that by definition, S 3 3 (M) = S 3 (M). It is now possible to refine the S 3 κ (M) in two different ways, either by splitting each quadrilateral element and thus refining the mesh M, or by performing spline refinement, i.e., increasing κ. We have the following inclusions of spline spaces, which follow directly from the definitions of the local spaces.
The inclusion S 3 (M ν ) ⊆ S 3 3ν (M) allows us to extend the results from Corollary 1 also to the spaces S 3 3ν (M), replacing h by h/ν. This is equivalent to an estimate with h being replaced byĥ, whereĥ is the maximum diameter of the polynomial pieces. Similar bounds can be obtained for spaces of the form S 3 3ν+1 (M) and S 3 3ν+2 (M), by allowing also 3 × 4, 3 × 5, 4 × 4 and 5 × 5 macro-elements. Thus, the approximation error bounds are optimal also for the spaces S 3 κ (M) obtained by spline refinement, i.e., for 0 ≤ ≤ 2 and 4 ≤ m ≤ p + 1 we have for all ϕ ∈ H m (Ω) for some constant C > 0 depending on ρ.
Similarly, one can also show bounds for spline spaces of higher degrees p ≥ 4 and regularity 1.

Basis construction
In the following, we describe how to compute the basis functions corresponding to one quadrilateral Q in the mesh. We define for every vertex six basis functions to interpolate the C 2 data. For every edge we define one basis function to interpolate the normal derivative at the edge midpoint. The remaining basis functions inside the element (with vanishing traces and derivatives on the element boundary) are selected to be standard Bernstein polynomials (for p = 5) or standard B-splines (for p ∈ {3, 4}). See [46,50] for basics on B-splines. In case of C 1 quadrilateral macro-elements, as in Definition 6, the patches Q are split into 2 × 2 (for p = 4) or 3 × 3 (for p = 3) polynomial sub-elements. Since we need κ ≥ max(1, 6 − p), they represent the spline elements with the least number of inner knots, allowing a separation of degrees of freedom at the vertices.
To simplify the construction, we build a basis with respect to a slightly modified dual basis. Instead of point evaluations at the interior, we use integral-based functionals that are dual to the Bernstein polynomials (or B-splines).
Before we go into the details, we discuss the Bernstein-Bézier (or B-spline) representation of basis functions. Letb 5 j be the Bernstein polynomials of degree 5, i.e., for 0 ≤ j ≤ 5 and ξ ∈ [0, 1],b 5 j (ξ ) = 5 j ξ j (1 − ξ) 5−j and letμ 5 i be the corresponding dual functionals, e.g., as in [25], i.e.,μ 5 i (b 5 j ) = δ j i . For p = 4 we consider the spline space with one inner knot at 1 2 with multiplicity two, denoted by S 4,2 2 , having the basisb 4 j and corresponding dual basisμ 4 j for 0 ≤ j ≤ 6. For p = 3, we consider the spline space with two inner knots at 1 3 and 2 3 with multiplicity two, denoted by S 3,1 3 , having the basisb 3 j and dual basiŝ μ 3 j for 0 ≤ j ≤ 7, see [46,50]. Let n = 11 − p. For each basis function β ∈ P p Q , the pull-back β = β • F Q possesses a tensor-product representation, having coefficients d j 1 ,j 2 ∈ R,  To simplify the construction, we replace the point evaluation functionals Λ p 2,Q by the dual functionals of mapped tensor-product Bernstein polynomials or B-splines Here, μ p j 1 ,j 2 form a dual basis for b in such a way that it is dual to

Patch interior basis functions
We have (n − 4) 2 basis functions

Edge basis functions
The edge basis functions β (k) ε , corresponding to edge ε (k) , are given by if n ε (k) points inwards. The To construct a global C 1 basis, the edge function must be positive on one side and negative on the other, which is achieved by selecting D[β

Vertex basis functions
The vertex basis functions β with vanishing C 2 -data at all other vertices. All representations of basis functions are verifiable via symbolic computation, e.g., by using Mathematica. We visualize the non-zero entries of the matrices V

Numerical examples
The goal is to demonstrate the potential of using the proposed C 1 spaces over quadrilateral meshes M for solving fourth order PDEs over domains Ω with piecewise linear boundary. This is done on the basis of a particular example, namely for the biharmonic equation is the resulting refined quadrilateral mesh obtained from M with corresponding sets of quadrilaterals Q h , edges E h and vertices V h . Note that in the refinement process, each quadrilateral of the current mesh is split regularly into four sub-quadrilaterals. Moreover, in all examples below, the functions g, g 1 and g 2 from problem (9.1) are computed from an exact solution u, and the resulting Dirichlet boundary data g 1 and g 2 are L 2 projected and strongly imposed to the numerical solution u h ∈ S p (M h ).
In more detail, we solve the biharmonic equation (9.1) in the following way: we split the such that and S p ∂ (M h ) spans all traces and normal derivatives at the boundary. Let The decomposition (9.2) allows to first solve the quadratic minimization problem (9.3), and then the variational problem (9.4). A detailed isogeometric formulation can be found, e.g., in [3,31]. Fig. 6 (top row, left), we solve the biharmonic equation (9.1) over the corresponding bilinear multi-patch domain by using the Brenner-Sung quadrilateral and macro-element spaces S p (M h ) for polynomial degrees p ∈ {3, 4, 5}. The considered exact solution is given by 5) and is shown in Fig. 6 (top row, right). The resulting L ∞ -errors as well as the relative L 2 -, H 1 -and H 2 -errors with respect to the number of degrees of freedom (NDOF) are shown in the middle and bottom rows of Fig. 6, and decrease with optimal orders

Example 2
We compare the C 1 quadrilateral spaces S p (M h ) for polynomial degrees p ∈ {3, 4, 5} as constructed in this paper with the C 1 isogeometric spaces A h for the cases (p, r) = (3, 1), (p, r) = (4, 2) and (p, r) = (5, 3) as generated in [30]. The spaces S p (M h ) and A h coincide at the coarsest level of discretization; however, A h is refined by standard spline knot-insertion, yielding C p−2 regularity at the new mesh lines. Therefore, when refining, it holds S p (M h ) ⊂ A h for (p, r) = (3, 1), and A h ⊂ S p (M h ) for (p, r) = (4, 2) and (p, r) = (5, 3). For this comparison, we solve the biharmonic equation (9.1) for the exact solution (9.5) see Fig. 7 (top row, right), on the bilinearly parametrized multi-patch domain Ω determined by the mesh M shown in Fig. 7 (top row, left). The resulting L ∞ -errors as well as the  , respectively. While the spaces S p (M h ) perform slightly better than the spaces A h for the case p = 3, it is in the opposite way around for the case p = 5, as expected due to the inclusions between the spaces S p (M h ) and A h recalled above.

Example 3
In the previous examples, the meshes were always nested, even though the spaces were not. Thus, when refining using a regular split, the elements tend to become closer in shape to parallelograms. This is not necessary for optimal convergence, as can be seen in the present example. Here, we reproduce the mesh refinement presented in [2, Figure 1], for the meshes M h as depicted in Fig. 8 (top row), and solve the biharmonic equation for the exact solution introduce local spaces for every element; hence, no uniform space on the parameter domain exists. Thus, we are guaranteed to reproduce polynomials of degree p, even though not all bivariate polynomials of maximum degree p are present on the parameter domain Q. This is different from [2], where a fixed space of polynomials on the parameter domain (e.g., for the serendipity element) yields optimal convergence rates on the presented mesh if and only if the space contains all bivariate polynomials of maximum degree p.

Example 4
The goal is to compare the Brenner-Sung quadrilateral with the Argyris triangle of degree p = 5, comparing the spaces S 5 (M h ) and S 5 (T h ), respectively, where T h is the resulting refined triangular mesh obtained via splitting each triangle in a regular way into four sub-triangles. For this, we solve the biharmonic for the computational domain from for the computational domain from Fig. 10 (top row, right), and fulfills for both cases homogeneous boundary conditions of order 1. While in Fig. 9 the more regular configuration is used for the quadrilateral mesh compared to the triangular one, it is in the opposite way around for the meshes in Fig. 10. The numerical results, which are shown in the middle and bottom rows of Figs. 9 and 10, and which are compared with respect to the number of degrees of freedom (NDOF), indicate that the Brenner-Sung quadrilateral spaces S 5 (M h ) perform significantly better than the Argyris triangle spaces S 5 (T h ) for the more "quad-regular" case (cf. Fig. 9) and just slightly worse for the more "triangle-regular" case (cf. Fig. 10). However, the rates are not affected, as in all considered instances, the resulting L ∞ -errors as well as the relative L 2 -,

Extension to isoparametric/isogeometric elements
It is possible to extend the construction from bilinearly mapped quadrilaterals Q, with F Q ∈ (Q 1 ) 2 , to fully isoparametric elements, with F Q ∈ (P p Q ) 2 . This leads to the C 1 multi-patch isogeometric space proposed in [30] and based on the previous works [14,27,28,31]. There are however two significant differences between this work and [30]. First, a local (i.e., element-wise) definition of the space and degrees of freedom is not known in the isogeometric case 1 . Furthermore, higher-degree element parametrizations need to fulfill a condition (named analysis-suitable G 1 in the papers mentioned above) that requires a suitable parametrization refitting, as in [28].
Modifications of elements near curved boundaries were also discussed and resolved successfully in [5] for a C 1 space of degree 4 and 5 over a quadrilateral mesh. There, the authors presented the construction of a minimal determining set (similar to a dual basis), without giving explicit degrees of freedom or a basis representation, and without a convergence analysis. Note that the proofs and techniques in the convergence analysis presented in this paper do not directly carry over to the case of elements with curved boundaries.
However, it is possible to obtain more general planar and surface parametrizations by template mapping. Given a planar template domain Ω T , which allows a quadrangulation with mesh M T , one can construct the C 1 space S p (M T ) over that quadrilateral mesh. One can then parameterize a domain of interest Ω = G(Ω T ) represented by a regular mapping G ∈ (S p (M T )) d , which could be planar with curved, Lipschitz boundaries (d = 2) or a surface (d = 3). On such a surface parametrization, one can construct a basis using the isoparametric principle. The approximation properties of S p (M T ) on Ω T then also carry over to the mapped space on Ω. Such a discretization may be used for solving fourth order PDEs on surfaces such as a Kirchhoff-Love shell problem. 1 In the definition of the space P p Q and the associated degrees of freedom as developed in [30], the normal derivative (∂ nε ϕ • F Q )|ε, is replaced by a suitable directional derivative (d · ∇ϕ • F Q )|ε: in the bilinear case d is the unitary normal n up to rescaling, while in general d depends on the adjacent patches.

Conclusion
We have described the construction of a novel family of C 1 quadrilateral finite elements, extending the Brenner-Sung quadrilateral construction from [9], possessing similar degrees of freedom as the classical Argyris triangle [1]. The presented method allows the simple design of polynomial as well as of spline elements. Among others, we have introduced a simple and local basis for the C 1 quadrilateral space, and have stated for particular cases explicit formulas for the Bézier or spline coefficients of the basis functions. We have also studied several properties of the C 1 quadrilateral space such as the optimal approximation properties of the space. Furthermore, the C 1 quadrilateral spaces are perfectly suited for solving fourth order PDEs, which has been demonstrated on the basis of several numerical examples solving the biharmonic equation on different quadrilateral meshes.
The classical Argyris triangle spaces and the Brenner-Sung quadrilateral spaces of degree p ≥ 5 not only possess similar degrees of freedom, but also reproduce the same functions as traces and normal derivatives along edges. Thus, they are compatible within a mixed triangle and quadrilateral mesh, cf. [22]. However, extensions to lower degrees p ∈ {3, 4} are not straightforward, as the macro-elements presented in this paper are based on B-spline refinement (knot insertion), which are not directly compatible with the most common triangle constructions of low degree based on Clough-Tocher or Powell-Sabin splits. The extension of our approach to quadrilateral meshes with curved boundaries and surfaces, going beyond the approach discussed in Section 10, is also of interest for future studies.

Appendix A: Pull-backs of derivatives
In the following, we derive the pull-backs of the partial derivatives ∂ 1 ϕ, ∂ 2 ϕ, ∂ 2 1 ϕ, ∂ 1 ∂ 2 ϕ and ∂ 2 2 ϕ, i.e., their representations in parameters (ξ 1 , ξ 2 ). We have that and, after applying the chain rule, obtain In addition, we define the correction terms We then have (i,j ) , i.e., b p (ξ 1 ) T R (k) (X)b p (ξ 2 ) , are constructed in such a way that they satisfy the C 2 interpolation conditions at the vertex v (k) and have vanishing derivatives up to second order at all others. Moreover, they have vanishing function values and gradients along the edges ε (k+1) and ε (k+2) . However, their normal derivatives along the edges ε (k) and ε (k−1) are (in general) not polynomials of degree four, i.e., they do not satisfy the condition (∂ n ε ( ) ϕ•F Q )|ε( ) ∈ P 4 in Eq. 3.1 for ∈ {k−1, k}, and their normal derivatives do not vanish at the edge midpoints, i.e., ∂ n ε ( ) ϕ(m ε ( ) ) = 0 for ∈ {k − 1, k}. The correction terms in Eq. B.1 are defined in such a way that both conditions are satisfied for the functions β (k) v,i , with i ∈ {0, 1, 2, 3, 4, 5}, defined in Eq. 8.1. The functions for degrees p = 4 and p = 3 are constructed analogously.

B.2 Precomputations for degree p = 4
For macro-element constructions of degree p = 3 and p = 4 the involved matrices are different from those for degree p = 5. For p = 4 we have and the correction terms   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/.