A nodal type polynomial finite element exact sequence over quadrilaterals
aa r X i v : . [ m a t h . NA ] O c t A nodal type polynomial finite element exact sequence overquadrilaterals ∗ Xinchen Zhou † a , Zhaoliang Meng b , Xin Fan c , and Zhongxuan Luo b,ca Faculty of Electronic Information and Electrical Engineering, Dalian University of Technology, Dalian 116024,China b School of Mathematical Sciences, Dalian University of Technology, Dalian, 116024, China c School of Software, Dalian University of Technology, Dalian, 116620, China
October 16, 2018
Abstract
This work proposes two nodal type nonconforming finite elements over convex quadrilaterals,which are parts of a finite element exact sequence. Both elements are of 12 degrees of freedom(DoFs) with polynomial shape function spaces selected. The first one is designed for fourth orderelliptic singular perturbation problems, and the other works for Brinkman problems. Numericalexamples are also provided.
Keywords:
Nodal type; polynomial finite element; exact sequence; quadrilateral meshes.
Let Ω ⊂ R be a simply connected Lipschitz domain. The de Rham complex determined by the followingexact sequence 0 H (Ω) (cid:2) H (Ω) (cid:3) L (Ω) 0 , curl div (1.1)also known as the Stokes complex, is well understood and widely applied in the analysis for manyproblems in solid and fluid mechanics. Typical model problems are biharmonic and Stokes problems,whose solutions can be efficiently approximated by suitable finite element methods. In particular, adivergence-free Stokes element often plays a role in a certain discretization of (1.1) with some biharmonicelement. A comprehensive review on this topic can be found in [15]. Generally speaking, there are threetypes of finite element exact sequences approximating (1.1). The first type are completely conforming,namely, all their components are subspaces of the corresponding forms in (1.1). Typical examples includethe sequences derived from the Argyris element [10], the singular Zienkiewicz element [14], the Bogner-Fox-Schmit element [18] and the family from spline or macroelements [7, 11, 1], etc. The second type aresemi-conforming, that is, their 0-forms are H -nonconforming but H -conforming, and their 1-forms are H -nonconforming but H (div)-conforming. The sequence constructed via the modified Morley element[20, 16] and its higher order extension [13] are of this type. The rectangular Adini element was also ∗ This project is supported by NNSFC (Nos. 61733002, 61572096, 61432003, 61720106005, 61502107) and “the Funda-mental Research Funds for the Central Universities”. † Corresponding author: [email protected] H (Ω) H (div; Ω) L (Ω) 0 . curl div (1.2)Indeed, the modification [20, 16] for the Morley-Crouzeix-Raviart sequence is a compromise for thisdilemma.Note that all the aforementioned examples are designed for triangular or rectangular meshes. However,there are fewer researches on the approximation for (1.1) and (1.2) over general convex quadrilaterals,on which we will give a brief review. For the first type approximation for (1.1), the H -conformingFraijes de Veubeke-Sander element [6, 23] is a successful candidate for biharmonic problems. Owing toa normal aggregation trick, a subspace method, namely, the reduced Fraijes de Veubeke-Sander elementwas designed [5]. For H -conforming approximation of the incompressible flow, Neilan and Sap [19]introduced a divergence-free Stokes element from the Fraijes de Veubeke-Sander element. As far as thesecond type approximation is concerned, Bao et al. [2] proposed a H -conforming element for fourth ordersingular perturbation problems by enriching a spline element space by bubble functions. Note that allthese elements are spline-based, and so a cell-refinement procedure cannot be avoided. Comparatively,polynomial shape functions are simple to represent and easy to compute, and therefore they are oftenmore preferred. This has been taken into consideration for the third type approximation. Utilizing thePark-Sheen biharmonic element [21], Zhang [27] generalized the Morley-Crouzeix-Raviart sequence togeneral quadrilaterals, but again there is no evidence of its ability to approximate (1.2). Recently, apolynomial modification was proposed by Zhou et al. [30], which works for both (1.1) and (1.2). We mustpoint out that, the Adini complex [12] and the rectangular Morley complex [25, 26] are also successfulfor the discretization of both (1.1) and (1.2), but their convergence severely relies on the regularity andsymmetry of the rectangular cell, and therefore cannot be generalized to arbitrary convex quadrilateralsin a obvious manner. Moreover, we discover that the number of global DoFs of the reduced Fraijes deVeubeke-Sander element [5] is significantly less than those of the semi-conforming [2] and completelynonconforming counterparts [30] benefitting from the nodal type structure.This work devotes to the construction of a nonconforming finite element exact sequence on generalconvex quadrilateral meshes for approximating both (1.1) and (1.2), enjoying the advantages that theelements therein are of nodal type structure, and their shape functions are polynomials. In fact, the0-form dealing with fourth order elliptic singular perturbation problems is, in a pseudo H -conformingmanner with respect to (1.2), a direct generalization of the modified nonconforming Zienkiewicz element[24] due to Wang, Shi and Xu. The DoFs are values and gradients of at vertices. For the 1-form designedfor Brinkman problems, we select vertex values and edge normal means as the DoFs. Optimal anduniform error estimates are also given for both elements with respect to their associated model problems.From some numerical tests, one can observe that the performances of both elements are consistent withour theoretical findings.The rest of this work is arranged as follows. In Section 2, the nonconforming finite element workingfor fourth order elliptic singular perturbation problems is defined on quadrilateral meshes. Section 3introduces the element for Brinkman problems, and shows that both the two elements are parts of a finiteelement exact sequence. Numerical examples are given in Section 4 to verify the theoretical analysis.Throughout the work, standard notations in Sobolev spaces are adopted. For a domain D ⊂ R , n and t will be the unit outward normal and tangent vectors on ∂D , respectively. The notation P k ( D )2enotes the usual polynomial space over D of degree no more than k . The norms and semi-norms oforder m in the Sobolev spaces H m ( D ) are indicated by k · k m,D and | · | m,D , respectively. The space H m ( D ) is the closure in H m ( D ) of C ∞ ( D ). We also adopt the convention that L ( D ) := H ( D ), wherethe inner-product is denoted by ( · , · ) D . These notations of norms, semi-norms and inner-products alsowork for vector- and matrix-valued Sobolev spaces, where the subscript Ω will be omitted if the domain D = Ω. Moreover, the positive constant C independent of the mesh size h and parameters ε , ν and α inthe model problems might be different in different places. Let K be an arbitrary convex quadrilateral. The four vertices of K are given by V , V , V , V in acounterclockwise order, and the i th edge of K is denoted by E i = V i V i +1 , whose equation is written as l i ( x, y ) = 0, i = 1 , , ,
4. Here and throughout the paper, the index i is taken modulo four. For each E i , M i denotes its midpoint, and n i and t i will be its unit normal and tangential directions. The equationsof lines through M M , M M , V V and V V read as m ( x, y ) = 0, m ( x, y ) = 0, l ( x, y ) = 0 and l ( x, y ) = 0, respectively. Moreover, we assume that all the aforementioned line equations are uniquelydetermined by l ( M ) = l ( M ) = l ( M ) = l ( M ) = m ( M ) = m ( M ) = l ( V ) = l ( V ) = 1 . (2.1)In order to describe the construction, we recall an auxiliary affine transformation for each K generatedby decomposing the standard bilinear mapping (see also [21, 9, 29]). The reference square b K = [ − , is determined by its vertices b V = ( − , − T , b V = (1 , − T , b V = (1 , T and b V = ( − , T . The bilinearmapping F K : b K → K such that b V i is mapped into V i for each i can be decomposed as F K = A K ◦ S K with A K : e K → K and S K : b K → e K defined by A K ( e x ) = A e x + b , S K ( b x ) = b x + b x b y s , e x = ( e x, e y ) T ∈ e K, b x = ( b x, b y ) T ∈ b K, where A is a 2 × b , d and s are two-dimensional vectors given by A = 14 ( V − V − V + V , V + V − V − V ) , d = 14 ( V − V + V − V ) , b = 14 ( V + V + V + V ) , s = ( s , s ) T = A − d . (2.2)Figure 1 gives an example of the intermediate reference element e K and the auxiliary affine transformation A K . We shall denote a point on e K by e V = ( e x, e y ) T if it equals A − K ( V ) for a point V = ( x, y ) T on K , andan edge in e K by e E if it equals A − K ( E ) for an edge E in K . Note that f M i is also the midpoint of b V i b V i +1 and e V i = b V i + ( − ( i +1) s . Furthermore, since K is convex, one must have | s | + | s | < . (2.3)Similarly, we write the function e f = f ◦ A K defined over e K for a function f over K . A simple calculationwill derive e l = 12 (cid:18) − s s − e x + e y + 1 (cid:19) , e l = 12 (cid:18) − e x + s s + 1 e y + 1 (cid:19) , e l = 12 (cid:18) s s + 1 e x − e y + 1 (cid:19) , e l = 12 (cid:18)e x − s s − e y + 1 (cid:19) , e l = − e x + e y + s − s s − s + 1) , e l = e x + e y + s + s s + s + 1) , e m = e x, e m = e y. (2.4)30 , e x e y b V b V b V b V e V e V e V e V ss e m = 0 e m = 0 e l = 0 e l = 0 e K A K KV V V V m = 0 m = 0 l = 0 l = 0Figure 1: Affine mapping A K from the intermediate reference quadrilateral e K to a general K . To design the element for fourth order singular perturbation problem, we first introduce an auxiliaryelement, which extends the rectangular Adini element to general convex quadrilaterals in a pseudo- C manner. Definition 2.1.
The quadrilateral finite element ( K, W − K , T − K ) is defined as follows: • K is a convex quadrilateral; • W − K = P ( K ) ⊕ span { φ , φ } is the shape function space where φ = ( s − s + 1) l l m m − s s l l m + s l l m m ,φ = ( s − s + 1) l l m m − s s l l m + s l l m m . The parameters s and s are defined via (2.2). • T − K = { τ j , j = 1 , , . . . , } is the DoF set where τ j ( w ) = w ( V j ) , ( τ j +4 ( w ) , τ j +8 ( w )) T = ∇ w ( V j ) , j = 1 , . . . , . Write p = l l l , p = l l l , p = l l l , p = φ and q = l l l , q = l l l , q = l l l , q = φ .We also define the nodal functionals λ ( w ) = | E | ∂w∂ t ( V ) , λ ( w ) = | E | ∂w∂ t ( V ) , λ ( w ) = | E | ∂w∂ t ( V ) , λ ( w ) = | E | ∂w∂ t ( V ) ,µ ( w ) = | E | ∂w∂ t ( V ) , µ ( w ) = | E | ∂w∂ t ( V ) , µ ( w ) = | E | ∂w∂ t ( V ) , µ ( w ) = | E | ∂w∂ t ( V )and the 4 × M and N by setting M i,j = λ i ( p j ), N i,j = µ i ( q j ), i, j = 1 , , ,
4. The followinglemma is helpful to verify the unisolvency of (
K, W − K , T − K ). Lemma 2.2.
The matrices M and N are nonsingular.Proof. For i, j = 1 , , ,
4, note that functionals | E i | ∂w∂ t i ( V j ) = | e E i | ∂ e w∂ e t i ( e V j ), therefore we can calculate theentries of M and N on e K rather than the physical K in variables s and s . Using (2.4) we set f ( s , s ) = ( s + s − s − s − s − s + 1) , f ( s , s ) = ( s − s + 1)( s + s + 1)( s − s + 1) ,f ( s , s ) = ( s + s + 1)( s − s − s − s + 1) , f ( s , s ) = ( s − s + 1)( s + s − s − s + 1) (2.5)4nd a direct computation gives M = (cid:0) M T , M T (cid:1) T and N = (cid:0) N T , N T (cid:1) T , where M = f (cid:18) ( s + s − / ( s −
1) 0 ( s − s − / ( s − s + 1) − ( s + 1) ( s + s − s − s + 1) / ( s −
1) 0 0 − ( s + 1) ( s + s − (cid:19) , M = f (cid:18) s + s + 1) / ( s + 1) 1 − ( s − ( s − s − s − s − / ( s + 1) 0 − ( s − ( s − s − (cid:19) , N = f (cid:18) − ( s − s + 1) / ( s + 1) − ( s + s − / ( s + s + 1) ( s − ( − s + s + 1)0 ( s + s + 1) / ( s + 1) 0 ( s − ( − s + s + 1) (cid:19) , N = f (cid:18) − ( s − s − / ( s −
1) 0 − − ( s + 1) ( s + s − s + s − / ( s −
1) 0 0 − ( s + 1) ( s + s − (cid:19) . Hence, by a symbolic computation,det M = 4 f f ( s − s − s + s − , det N = 4 f f ( s + s − s + s − . It then follows from (2.3) that det M = 0 and det N = 0, which is the desired result. Lemma 2.3.
The element ( K, W − K , T − K ) is well-defined.Proof. Set r = l l , r = l l , r = l l , r = l l , then our goal is to show W − K = span { p i , q i , r i , i =1 , , , } and that τ j ( w ) = 0 for w ∈ W − K will derive w = 0. First we see all p i , q i and r i are linearlyindependent. Indeed, if w = X i =1 ( α i p i + β i q i + γ i r i ) = 0 , (2.6)for α i , β i , γ i ∈ R , then τ i ( w ) = 0 for i = 1 , , ,
4. Noting that τ i ( p j ) = τ i ( q j ) = 0 , τ i ( r j ) = r i ( V i ) δ ij , r i ( V i ) = 0 , i, j = 1 , , , , we obtain γ i = 0, i = 1 , , ,
4. Moreover, the constructions of all q i imply λ i ( q j ) = 0 , µ i ( p j ) = 0 , i, j = 1 , , , , and therefore it follows from λ i ( w ) = 0 and µ i ( w ) = 0 that X j =1 α j λ i ( p j ) = 0 , X j =1 β j µ i ( q j ) = 0 . Taking i = 1 , , ,
4, we find M ( α , α , α , α ) T = N ( β , β , β , β ) T = , which implies α = α = α = α = β = β = β = β = 0 according to Lemma 2.2. Hence, all p i , q i and r i are linearly independent, namely, dim span { p i , q i , r i , i = 1 , , , } = dim W − K = 12. Moreover,all p i , q i and r i are members of W − K , and thus these two sets are equal. Repeat the same process abovefrom (2.6) apart from the assumption w = 0, the unisolvency is derived since τ j ( w ) = 0 for j = 5 , . . . , λ i ( w ) = µ i ( w ) = 0, i = 1 , , , φ and φ are necessary. 5 emma 2.4. For all w ∈ W − K , it holds that | E i | Z E i w d s = 12 ( w ( V i ) + w ( V i +1 )) − | E i | (cid:18) ∂w∂ t i ( V i +1 ) − ∂w∂ t i ( V i ) (cid:19) , i = 1 , , , . (2.7) Proof.
For each i , let ξ i ∈ P ( E i ) be taken such that ξ i ( V i ) = − ξ i ( V i +1 ) = 1. Then by integratingby parts, one must have1 | E i | Z E i w d s = 12 (cid:18) w ( V i ) + w ( V i +1 ) − Z E i ∂w∂ t i ξ i d s (cid:19) = 0 , i = 1 , , , , ∀ w ∈ C ( K ) . (2.8)If w ∈ P ( K ), the Simpson quadrature rule implies1 | E i | Z E i ∂w∂ t i ξ i d s = 16 (cid:18) ∂w∂ t i ( V i +1 ) − ∂w∂ t i ( V i ) (cid:19) , which along with (2.8) leads to (2.7). For φ and φ , we can directly calculate on e K that1 | e E i | Z e E i e φ j d e s = 12 ( e φ j ( e V i ) + e φ j ( e V i +1 )) − | e E i | ∂ e φ j ∂ e t i ( e V i +1 ) − ∂ e φ j ∂ e t i ( e V i ) ! = 0 , i = 1 , , , , j = 1 , , and so the assertion is verified. Remark 2.5. If K is a rectangle, s = s = 0 . This element degenerates to the traditional Adini element. Let us now introduce the following H ( K )-bubble function space W bK over K : W bK = span { b , b x, b y, b l l } , b = l l l l . (2.9)The DoFs in T bK with respect to these bubble functions are τ bi ( w ) = 1 | E i | Z E i ∂w∂ n i d s, i = 1 , , , . Lemma 2.6. If w ∈ W bK such that all τ bi ( w ) = 0 , i = 1 , , , , then w = 0 .Proof. Over K we define the following polynomials b = 1, b = m , b = m , b = l l and b i,j = b b j /l i for i, j = 1 , , ,
4, then we turn to the intermediate reference quadrilateral e K . Considerthe 4 × B − determined by B − i,j = 1 | E i | Z E i b i,j d s = 1 | e E i | Z e E i e b i,j d e s, i, j = 1 , , , . By (2.4) and a direct computation, we find B − = (cid:0) ( B − ) T , ( B − ) T , ( B − ) T , ( B − ) T (cid:1) T , where B − = 16 f (cid:18) − , s ( s − s + 1) , s + 5 s + 55( s + 1) , − ( s + s − s − s − s + s + 1)( s − s + 1) (cid:19) , B − = 16 f (cid:18) , − s − s + 55( s − , − s ( s + 1)5( s − , s − s − s − s + 1) (cid:19) , B − = 16 f (cid:18) − , s ( s + 1)5( s − , − − s + 5 s − s − , − (cid:19) , B − = 16 f (cid:18) , − s + 5 s + 55( s + 1) , − s ( s − s + 1) , s + s − s + s + 1) (cid:19) , f i , i = 1 , , , B − = f f f f (( s + s ) − s s ( s + s ) + 9( s + s ) − s s + 15( s + s ) − s − s + 1)( s − s + 1)( s + s + 1)( s − s + 1) . Note from (2.3) that the factor of the numerator( s + s ) − s s ( s + s ) + 9( s + s ) − s s + 15( s + s ) − ≤ ( s + s ) + 9( s + s ) + 15( s + s ) − < B − is nonsingular.Next, we turn to the matrix B with B i,j = τ bi ( b b j ). Our aim is to show B is nonsingular as W bK = span { b b j , j = 1 , , , } . For each i and j , since B i,j = 1 | E i | (cid:18)Z E i ∂l i ∂ n i b i,j d s + Z E i ∂ ( b i,j ) ∂ n i l i d s (cid:19) = ∂l i ∂ n i B − i,j , then det B = Y i =1 ∂l i ∂ n i ! det B − = 0 , which completes the proof.We are in a position to propose the element for fourth order elliptic singular perturbation problemsby the normal aggregation strategy. The DoFs are at vertices over a quadrilateral, which is a nodal typeconstruction. Definition 2.7.
The quadrilateral finite element ( K, W K , T K ) is defined by: • K is a convex quadrilateral; • W K is the shape function space: W K = (cid:26) w ∈ W − K ⊕ W bK : 1 | E i | Z E i ∂w∂ n i d s = 12 (cid:18) ∂w∂ n i ( V i ) + ∂w∂ n i ( V i +1 ) (cid:19) , i = 1 , , , (cid:27) ; (2.10) • T K = T − K is the DoF set.Here W bK and ( K, W − K , T − K ) have been given in (2.9) and Definition 2.1, respectively. Theorem 2.8.
The element ( K, W K , T K ) is well-defined. Moreover, (2.7) holds for all w ∈ W K and P ( K ) ⊂ W K .Proof. The four relations in (2.10) hint that dim W K ≥
12. It suffices to show if w ∈ W K fulfilling τ j ( w ) = 0 for j = 1 , , . . . ,
12 then w = 0. Write w = w − + w b with w − ∈ W − K and w b ∈ W bK .The assumption above gives τ j ( w − ) = 0 as τ j ( w b ) = 0 for all j , and by Lemma 2.3 we find w − = 0.Moreover, the relations in (2.10) will lead to τ bi ( w ) = τ bi ( w b ) = 0, i = 1 , , ,
4, and by Lemma 2.6 we get w b = 0, which means w = 0 and the unisolvency has been derived. The next assertion is trivial as w b arebubble functions. By noting that the relations in (2.10) hold for P ( K ), the last assertion immediatelyfollows.The nodal basis representation can be obtained as follows. Let ψ − j ∈ W − K , j = 1 , , . . . ,
12 be thenodal basis of (
K, W − K , T − K ) and ψ bj ∈ W bK , j = 1 , , , τ bi ( ψ bj ) = δ ij , i = 1 , , ,
4. Then ψ j = ψ − j + X i =1 c i,j ψ bi , j = 1 , , . . . , K, W K , T K ), where the coefficients c i,j are determinedthrough c i,j = 12 ∂ψ − j ∂ n i ( V i ) + ∂ψ − j ∂ n i ( V i +1 ) ! − | E i | Z E i ∂ψ − j ∂ n i d s. Let Ω ⊂ R be a polygonal domain and ∂ Ω be its boundary. For a given f ∈ L (Ω), the fourth orderelliptic singular perturbation problem appears as: Find u such that ε ∆ u − ∆ u = f in Ω ,u = ∂u∂ n = 0 on ∂ Ω , (2.11)where ε is the singular perturbation parameter tending to zero. A weak formulation is to find u ∈ H (Ω)such that ε ( ∇ u, ∇ v ) + ( ∇ u, ∇ v ) = ( f, v ) , ∀ v ∈ H (Ω) . (2.12)Let {T h } be a family of quasi-uniform and shape-regular partitions of Ω consisting of convex quadri-laterals. For a cell K ∈ T h , h K denotes the diameter of K , and so the parameter h := max K ∈T h h K .The sets of all vertices, interior vertices, boundary vertices, edges, interior edges and boundary edgesare correspondingly denoted by V h , V ih , V bh , E h , E ih and E bh . For each E ∈ E h , n E is a fixed unit vectorperpendicular to E and t E is a vector obtained by rotating n E by ninety degree counterclockwisely.Moreover, for E ∈ E ih , the jump of a function v across E is defined as [ v ] E = v | K − v | K , where K and K are the cells sharing E as a common edge, and n E points from K to K . For E ∈ E bh , we set[ v ] E = v | K if E is an edge of K .We now define the finite element space W h by setting W h = n w ∈ L (Ω) : w | K ∈ W K , ∀ K ∈ T h , w and ∇ w are continuous at all V ∈ V ih and vanishes at all V ∈ V bh o . Then the finite element approximation of (2.12) is: Find u h ∈ W h fulfilling ε a h ( u h , v h ) + b h ( u h , v h ) = ( f, v h ) , ∀ v h ∈ W h , (2.13)where a h ( u h , v h ) = X K ∈T h ( ∇ u h , ∇ v h ) K , b h ( u h , v h ) = X K ∈T h ( ∇ u h , ∇ v h ) K . Moreover, we define a discrete semi-norm by setting ||| v ||| ε,h = ε | v | ,h + | v | ,h with | v | m,h = X K ∈T h | v | m,K , m = 1 , . Clearly, owing to the definition of W h , we observe that |||·||| ε,h is a norm on V h . Thus, by the Lax-Milgramlemma, the problem (2.13) has unique solution.For each K ∈ T h and s >
0, we define the interpolation operator I K : H s ( K ) → W K according toTheorem 2.8 such that τ j ( I K w ) = τ j ( w ), j = 1 , , . . . ,
12. Since I K v = v for all v ∈ P ( K ) and {T h } isquasi-uniform and shape regular, we find | w − I h w | j,K ≤ Ch k − j | w | k,K , ∀ w ∈ H k ( K ) ∩ H s ( K ) , j = 0 , , , k = 2 , . (2.14)8hen the global interpolation operator I h : H (Ω) ∩ H s (Ω) → W h is set as I h | K = I K . Moreover,owing to Theorem 2 .
8, one has Z E [ w h ] E d s = Z E (cid:20) ∂w h ∂ n E (cid:21) E d s = 0 , ∀ E ∈ E h (2.15)via (2.7) and (2.10). The Strang lemma says ||| u − u h ||| ε,h ≤ C inf v h ∈ W h ||| u − v h ||| ε,h + sup w h ∈ W h E ε,h ( u, w h ) ||| w h ||| ε,h ! with the consistency error E ε,h ( u, w h ) = ε a h ( u, w h ) + b h ( u, w h ) − ( f, w h ) . Hence, applying the proof of Theorem 1 in Chen et al.’s work [4] and invoking (2.14), (2.15), we getinf v h ∈ W h ||| u − v h ||| ε,h ≤ Ch ||| u − I h u ||| ε,h ≤ Ch ( ε | u | + | u | ) ,E ε,h ( u, w h ) ≤ Ch ( ε | u | + | u | + k f k ) ||| w h ||| ε,h , ∀ w h ∈ W h , which leads to the following convergence result. Theorem 2.9.
Let u ∈ H (Ω) and u h ∈ W h be the solutions of (2.12) and (2.13), respectively. Then ||| u − u h ||| ε,h ≤ Ch ( ε | u | + | u | + k f k ) . From Theorem 2.9, this element ensures a linear convergence order with respect to h , uniformly in ε ,provided that ε | u | , | u | are uniformly bounded. However these terms might blow up when ε tends tozero. The next result, following a similar line of Theorem 4.3 in [25], guarantees a uniform convergencerate under the impact of such boundary layers. Theorem 2.10.
Assume Ω is a convex domain. Let u ∈ H (Ω) and u h ∈ W h be the solutions of (2.12)and (2.13), respectively. Then it holds the uniform error estimate ||| u − u h ||| ε,h ≤ Ch k f k . Let us turn to the construction of a vector-valued nodal type finite element. As in the scalar case, weshall prove the unisolvency of two auxiliary elements. For a general convex quadrilateral K , the element( K, V − K , Σ − K ) and the H (div; K )-bubble element ( K, V bK , Σ bK ) are defined through V − K = [ P ( K )] ⊕ span { curl x , curl x y, curl xy , curl y , curl φ , curl φ } , Σ − K = { σ j , j = 1 , , . . . , } , V bK = curl W bK , Σ bK = { σ bj , j = 1 , , , } , where φ , φ have been given in Definition 2.1, and the DoFs are σ j ( v ) = Z E j v · n j d s, ( σ j +4 ( v ) , σ j +8 ( v )) T = v ( V j ) , σ bj ( v ) = Z E j v · t j d s, j = 1 , , , . emma 3.1. Both ( K, V − K , Σ − K ) and ( K, V bK , Σ bK ) are well-defined.Proof. We only deal with ( K, V − K , Σ − K ) as the latter is much simpler. Define v = ( x, y ) T , then V − K = span { v } ⊕ span { curl w : w ∈ W − K } . (3.1)Suppose v = c v + curl w ∈ V − K for some w ∈ W − K such that σ j ( v ) = 0 for all j , then cσ j ( v ) + σ j ( curl w ) = 0 , j = 1 , , . . . , . (3.2)However, we notice that div v = 2, and by Green’s formula X i =1 σ i ( curl w ) = Z K div curl w d x = 0 , X i =1 σ i ( v ) = Z K div v d x = 2 | K | 6 = 0 . Summing over (3.2) for j = 1 , , , c = 0, and therefore σ j ( curl w ) = 0, j = 1 , , . . . ,
12. Hence, itsuffices to show curl w = . Indeed, we can select w such that w ( V ) = 0 without changing the value ofeach σ j ( curl w ). Then w ( V i +1 ) = w ( V i ) + Z E i ∂w∂ t d s = w ( V i ) + σ i ( curl w ) = 0 , i = 1 , , . (3.3)Moreover, ∇ w ( V j ) = ( − σ j +8 ( curl w ) , σ j +4 ( curl w )) T = , j = 1 , , , . (3.4)As a consequence of (3.3),(3.4) and Lemma 2.3, we find w = 0, and so curl w = , which implies v = .The proof is done.Parallel to Lemma 2.7, the following fact is crucial for the convergence in the Darcy limit. Lemma 3.2.
For all v ∈ V − K , it holds that | E i | Z E i v · n ξ i d s = 16 ( v ( V i +1 ) − v ( V i )) · n , i = 1 , , , , (3.5) where ξ i ∈ P ( E i ) has been defined in the proof of Lemma 2.4.Proof. According to (3.1), we shall verify this relation for v and all curl w , w ∈ W − K . Since v ∈ [ P ( K )] ,the Simpson quadrature rule ensures (3.5). On the other hand, if w ∈ W − K , substituting (2.7) into (2.8)will derive (3.5) for v = curl w , which completes the proof.Now we introduce the nodal type vector-valued element for Brinkman problems. Definition 3.3.
The finite element ( K, V K , Σ K ) is determined through: • K is a convex quadrilateral; • V K is the shape function space: V K = (cid:26) v ∈ V − K ⊕ V bK : 1 | E i | Z E i v · t i d s = 12 ( v ( V i ) + v ( V i +1 )) · t i , i = 1 , , , (cid:27) ; (3.6) • Σ K = Σ − K is the DoF set. Theorem 3.4.
The element ( K, V K , Σ K ) is well-defined. Moreover, (3.5) holds for all v ∈ V K and [ P ( K )] ⊂ V K .Proof. The proof is very similar to that of Theorem 2.8 and thus omitted.10 .2 Applied to Brinkman problems
Consider the following Brinkman problem of porous media flow over Ω: For given f ∈ [ L (Ω)] and g ∈ L (Ω), find the velocity u and the pressure p satisfying − div ( ν ∇ u ) + α u + ∇ p = f in Ω , div u = g in Ω , u = on ∂ Ω . (3.7)Here we assume that parameters ν, α ≥ να = 0. A weak formulation of (3.7) is tofind ( u , p ) ∈ [ H (Ω)] × L (Ω) satisfying a ( u , v ) − b ( v , p ) = ( f , v ) , ∀ v ∈ [ H (Ω)] ,b ( u , q ) = ( g, q ) , ∀ q ∈ L (Ω) , (3.8)with the bilinear forms a ( u , v ) = ν ( ∇ u , ∇ v ) + α ( u , v ) , b ( v , q ) = (div v , q ) . This problem has a unique solution due to the following inf-sup conditionsup v ∈ [ H (Ω)] b ( v , q ) k v k ≥ C k q k , ∀ q ∈ L (Ω) (3.9)according to [3] for all possible ν and α .Let {T h } be given as in Subsection 2.4. We select the following finite element spaces V h and P h : V h = n v ∈ [ L (Ω)] : v | K ∈ V K , ∀ K ∈ T h , Z E [ v · n E ] E d s = 0 for all E ∈ E h , and v is continuous at all V ∈ V ih and vanishes at all V ∈ V bh o .P h = (cid:8) q ∈ L (Ω) : q | K ∈ P ( K ) , ∀ K ∈ T h (cid:9) . If we write div h | K = div on K , then we have the divergence-free condition div h V h ⊂ P h . A discreteformulation of (3.8) will be given as: Find ( u h , p h ) ∈ V h × P h , such that a h ( u h , v h ) − b h ( v h , p h ) = ( f , v h ) , ∀ v h ∈ V h ,b h ( u h , q h ) = ( g, q h ) , ∀ q h ∈ P h , (3.10)where a h ( · , · ) and b h ( · , · ) are discrete versions of a ( · , · ) and b ( · , · ), respectively: a h ( u h , v h ) = ν X K ∈T h ( ∇ u h , ∇ v h ) K + α ( u h , v h ) , b h ( v h , q h ) = (div h v h , q h ) . Moreover, the norm k · k ,h and the semi-norms | · | ,h , k · k a h are equipped by k v h k ,h = X K ∈T h k v h k ,K , | v h | ,h = X K ∈T h | v h | ,K , k v h k a h = a h ( v h , v h ) . Clearly, k · k a h is a norm on V h .For each K ∈ T h and s >
0, the nodal interpolation operator Π K : [ H s ( K )] → V K is defined via σ j ( Π K v ) = σ j ( v ), j = 1 , , . . . ,
12. Like the scalar case, we have from Theorem 3.4 that | v − Π K v | j,K ≤ Ch k − j | v | k,K , ∀ v ∈ [ H k ( K ) ∩ H s ( K )] , j = 0 , , k = 1 , . (3.11)11he global interpolation operator Π h : [ H (Ω) ∩ H s (Ω)] → V h is naturally set as Π h | K = Π K . Since Π K preserves normal integral on E ⊂ ∂K for all K , we find through integrating by parts that b h ( Π h v , q h ) = b ( v , q h ) , ∀ v ∈ [ H (Ω) ∩ H s (Ω)] , ∀ q h ∈ P h . (3.12)Owing to the Scott-Zhang smoothing strategy [22], Π h can be modified into Π h interpolating con-tinuously from [ H (Ω)] to V h . Meanwhile, (3.12) holds for all v ∈ [ H (Ω)] if Π h is replaced by Π h .Hence, by Fortin’s trick and (3.9), the following discrete inf-sup condition is derived:sup v h ∈ V h b h ( v h , q h ) k v h k ,h ≥ sup v ∈ [ H (Ω)] b h ( Π h v , q h ) k Π h v k ,h ≥ sup v ∈ [ H (Ω)] b ( v , q h ) C k v k ≥ C k q h k , ∀ q h ∈ P h . (3.13)Then by Theorem 3.1 in [30], (3.10) has a unique solution ( u h , p h ) ∈ V h × P h , and k u − u h k a h ≤ C (cid:18) inf v h ∈ Z h ( g ) k u − v h k a h + sup w h ∈ V h E h ( u , p, w h ) k w h k a h (cid:19) , k p − p h k ≤ C (cid:20) k p − P h p k + M / (cid:18) inf v h ∈ Z h ( g ) k u − v h k a h + sup w h ∈ V h E h ( u , p, w h ) k w h k a h (cid:19)(cid:21) , (3.14)where P h is the L -projection operator from L (Ω) to P h , M = max { ν, α } and Z h ( g ) = { v h ∈ V h : b h ( v h , q h ) = ( g, q h ) , ∀ q h ∈ P h } ,E h ( u , p, w h ) = X K ∈T h (cid:18) − ν Z ∂K ∂ u ∂ n · w h d s + Z ∂K p w h · n d s (cid:19) . Now we are in a position to estimate each term in (3.14). To this end, let u ∈ [ H (Ω) ∩ H (Ω)] bethe weak velocity solution of (3.8). It follows from (3.12) that Π h u ∈ Z h ( g ), and therefore by (3.11)inf v h ∈ Z h ( g ) k u − v h k a h ≤ k u − Π h u k a h ≤ Ch ( ν / + α / h ) | u | . (3.15)On the other hand, by (3.5) and (3.6), Theorem 3.4 ensures Z E q [ v · n E ] E d s = 0 , ∀ q ∈ P ( E ) , Z E [ v · t E ] E d s = 0 , ∀ E ∈ E h . If p ∈ H (Ω), then following the spirit of the consistency error analysis in [30], we have E h ( u , p, w h ) ≤ (cid:26) Ch ( ν / | u | + ν − / h | p | ) , if ν = 0; Ch ( ν / | u | + α − / | p | ) , if α = 0 . (3.16)Substituting (3.15) and (3.16) into (3.14), we will obtain the following convergence result. Theorem 3.5.
Let ( u , p ) ∈ (cid:0) [ H (Ω) ∩ H (Ω)] (cid:1) × ( L (Ω) ∩ H (Ω)) be the weak solution of (3.8). Thediscrete solution of (3.10) is given by ( u h , p h ) ∈ V h × P h . Then the following error estimates hold: k u − u h k a h ≤ Ch h ( ν / + α / h ) | u | + min { C ν − / h, C α − / }| p | i , k p − p h k ≤ Ch n | p | + M / h ( ν / + α / h ) | u | + min { C ν − / h, C α − / }| p | io , (3.17) where we set α − / = + ∞ if α = 0 , and ν − / = + ∞ if ν = 0 .
12s the scalar case, boundary layers might appear if ν →
0. In such a Darcy limit, | u | , | p | and | p | might explode. We need a uniform convergence result instead of Theorem 3.5. To this end, Ω is assumedto be a convex polygonal domain with vertices x j , j = 1 , . . . , N on ∂ Ω. We also introduce the space H (Ω) = (cid:26) q ∈ H (Ω) ∩ L (Ω) : Z Ω | q ( x ) | | x − x j | d x < ∞ , j = 1 , . . . , N (cid:27) with the norm k q k , + = k q k + N X j =1 Z Ω | q ( x ) | | x − x j | d x . The following result is an analogue counterpart of Theorem 3.3 in [30], whose proof will be omitted.
Theorem 3.6.
Assume that Ω is convex, and α = 1 , ν ≤ in (3.8). Moreover, the known terms f ∈ [ H (Ω)] and g ∈ H (Ω) . Let ( u , p ) ∈ [ H (Ω)] × L (Ω) be the weak solution of (3.8). The discretesolution of (3.10) is given by ( u h , p h ) ∈ V h × P h . Then we have the following uniform error estimate k u − u h k a h + k p − p h k ≤ Ch / ( k f k + k g k , + ) . In this section, we will see that the finite element spaces W h , V h and P h constitute a discrete de Rhamcomplex. Theorem 3.7.
The following finite element sequence is exact. W h V h P h , curl h div h (3.18) where curl h | K = curl on K .Proof. We have shown div h V h ⊂ P h . Moreover, we know that div h is surjective due to the discreteinf-sup condition (3.13). Next, we show curl h W h ⊂ V h . On one hand, for all K ∈ T h , owing to thedefinitions of V − K and V bK , one has curl ( W − K ⊕ W bK ) ⊂ V − K ⊕ V bK . Furthermore, by the proof of Lemma3.2, the relation (3.6) holds for v = curl w , ∀ w ∈ W K . Thus we find curl W K ⊂ V K , ∀ K ∈ T h . On theother hand, ∀ w h ∈ W h , the definition of W h ensures the continuous conditions in the definition of V h for v h = curl h w h , which gives curl h W h ⊂ V h . To verify curl h W h = Z h := { v h ∈ V h : div h v h = 0 } , itsuffices to show the dimensions of this two spaces are the same since curl h W h ⊂ Z h . Let N i V , N i E and N K be the numbers of interior vertices, interior edges and cells in T h , respectively. Then by using Euler’sformula N i V − N i E + N K = 1, we havedim Z h = dim V h − dim (div h V h ) = dim V h − dim P h = (2 N i V + N i E ) − ( N K −
1) = 3 N i V = dim W h = dim( curl h W h ) , which implies the exactness of the sequence. Remark 3.8.
With a Scott-Zhang smoothing trick [22] acting on ∇ H (Ω) , the interpolation operator I h in (2.14) can be modified into I h to work on the whole H (Ω) rather than H (Ω) ∩ H s (Ω) (seee.g. [14]). As a consequence, we have the following commutative diagram: H (Ω) (cid:2) H (Ω) (cid:3) L (Ω) 00 W h V h P h . curl I h div Π h P h curl h div h (3.19)13 emark 3.9. We end this section by remarking that, the exact sequence (3.18) and commuting diagram(3.19) can be adapted to a more general mesh type, namely, mixed meshes consisting of both triangles andquadrilaterals, in light of the pseudo- C property of W h and the pseudo- H (div) property of V h . In fact, fora quadrilateral cell K , we can still select W K as in (2.10) and V K as in (3.6). But if K is a triangle, themodified nonconforming Zienkiewicz element space due to Wang et al. [24] will be a successful candidatefor W K , and V K can also be obtained in a similar manner as in Definition 3.3 from W K . Then W h , V h are formulated as before, and analogous counterparts of the error estimates Theorems 2.9, 2.10, 3.5 and3.6 are also appropriate. Some numerical examples are provided in this section. Let the solution domain Ω be the unit square[0 , , where three types of convex quadrilateral meshes are considered. As for the first type, eachmesh T h is generated by an n × n uniform rectangular partition. Figure 2(a) provides an example.Meshes of the second type consists of uniform trapezoids, see Figure 2(b). As shown in Figure 2(c),the random partitions are demonstrated as well, which are generated by stochastically deforming thefirst-type partitions with at most 20%. In the following examples, the 16-node Gauss quadrature rule isadopted when the entries of stiffness matrices are accumulated for all meshes. (a) A mesh with uniform rectan-gular partition (b) A mesh with uniform trape-zoidal partition (c) A mesh with nonuniform ran-domly perturbed partition Figure 2: Three types of quadrilateral partitions of Ω.Before the numerical experiment, we give a brief analysis in terms the computational cost in com-parison with some other elements working for the same problems over the same meshes reviewed inthe introduction part. For fourth order elliptic singular perturbation problems, we have reviewed two H -nonconforming elements in literature, see [2] for the H -conforming construction and [30] for the H -nonconforming one. For our tested n × n meshes, both elements are edge-based and the numbersof global DoFs are about 5 n . As far as W h in this work is concerned, this number will be about 3 n ,reducing the computational costs in some degree benefiting from the nodal type structure. The reduced H -conforming Fraijes de Veubeke-Sander element [5] has the same DoFs as ours, but the shape functionspace is spline-based, which is less preferred in practical applications than our polynomial selection. ForBrinkman problems, we investigate the H -conforming construction in [19] and a nonconforming one in[30]. Again, the number of global DoFs of V h in this work is about 4 n , much less than those of the twoaforementioned examples: 8 n for the element in [19] and 6 n for the other.We now check the performance of the finite element space W h applied to fourth order elliptic singularperturbation problems. The exact solution of (2.11) is arranged as u = sin (2 πx ) sin (2 πy ) . (4.1)14n Table 1, we list the errors in the energy norm ||| u − u h ||| ε,h with different values of ε and h . The resultsfor the biharmonic equation ∆ u = f as well as the Poisson problem − ∆ u = f with pure Dirichletboundary conditions are also presented. As predicted in Theorem 2.9, the first order convergence rate isobserved for all possible ε . ε n = 4 n = 8 n = 16 n = 32 n = 64 orderrectangular meshes:biharmonic 2.909E0 1.315E0 5.913E-1 2.804E-1 1.368E-1 1.041 2.913E0 1.315E0 5.914E-1 2.804E-1 1.368E-1 1.042 − − − − − − ||| u − u h ||| ε,h produced by W h applied to the given fourth order elliptic singularperturbation problem through (4.1) over three kinds of meshes.We then turn to the performance of the mixed finite element pair V h × P h applied to the Brinkmanproblem. We fix the parameter α = 1 in (3.7) and test different ν ∈ (0 , ν = 0, α = 1) and Stokes problem ( ν = 1, α = 0) are also investigated. The exact solution of(3.7) is determined by u = curl (sin ( πx ) sin ( πy )) , p = sin( πx ) − /π. (4.2)Tables 2 and 3 show the velocity and pressure errors, respectively. The optimal convergence rate is alsoachieved for all possible parameters. References [1] D. N. Arnold, J. Qin. Quadratic velocity/linear pressure Stokes elements, in Advances in Computer Methodsfor Partial Differential Equations VII, eds. R. Vichnevetsky and R.S. Steplemen, 1992.[2] Y. Bao, Z. Meng, Z. Luo. A C -Nonconforming Quadrilateral Finite Element for the Fourth-Order EllipticSingular Perturbation Problem. ESAIM Math. Model. Numer. Anal., 2018, to appear.[3] D. Boffi, F. Brezzi, M. Fortin. Mixed Finite Element Methods and Applications. Springer, 2013.[4] S. Chen, Y. Zhao, D. Shi. Non C nonconforming elements for elliptic fourth order singular perturbationproblem. J. Comput. Math., 23(2): 185–198, 2005.[5] P. G. Ciarlet. The Finite Element Method for Elliptic Problems, Studies in Mathematics and Its Applications,Vol. 4. North-Holland Publishing Company, Amsterdam-New York-Oxford, 1978. / n = 4 n = 8 n = 16 n = 32 n = 64 orderrectangular meshes:Stokes 3.186E0 1.503E0 6.926E-1 3.324E-1 1.631E-1 1.031 3.190E0 1.503E0 6.927E-1 3.324E-1 1.631E-1 1.032 − − − − − − k u − u h k a h produced by V h × P h applied to the Brinkman problem through(4.2) over three kinds of meshes. [6] J. F. Ciavaldini, J. C. N´ed´elec. Sur l’´el´ement de Fraeijs de Veubeke et Sander. RAIRO Analyse num´erique,1974, 8: 29–45.[7] S. H. Christiansen, K. Hu. Generalized Finite Element Systems for smooth differential forms and Stokesproblem. Numer. Math., 2018, 2(140): 327–371.[8] M. Crouzeix, P. -A. Raviart. Conforming and nonconforming finite element methods for solving the stationaryStokes equations. RAIRO 7, 1973, 3: 33–76.[9] E. Dubach, R. Luce, J. M. Thomas. Pseudo-conforming polynomial finite elements on quadrilaterals. Int. J.Comput. Math., 2009, 86(10-11): 1798–1816.[10] R. Falk, M. Neilan. Stokes complexes and the construction of stable finite elements with pointwise massconservation. SIAM J. Numer. Anal., 2013, 51(2): 1308–1326.[11] G. Fu, J. Guzm´an, M. Neilan. Exact smooth piecewise polynomial sequences on Alfeld splits. arXiv preprintarXiv:1807.05883, 2018.[12] A. Gillette, K. Hu, S. Zhang. Nonstandard finite element de Rham complexes on cubical meshes. arXivpreprint arXiv:1804.04390, 2018.[13] J. Guzman, M. Neilan. A family of nonconforming elements for the Brinkman problem. IMA J. Numer.Anal., 2012, 32: 1484–1508.[14] J. Guzm´an, M. Neilan. Conforming and divergence-free Stokes elements on general triangular meshes. Math.Comput., 2014, 83(285): 15–36.[15] V. John, A. Linke, C. Merdon, M. Neilan, L. G. Rebholz. On the divergence constraint in mixed finiteelement methods for incompressible flows. SIAM Review, 2017, 59(3): 492–544.[16] K. A. Mardal, X.-C Tai, R. Winther. A robust finite element method for Darcy-Stokes flow. SIAM J. Numer.Anal., 2002, 40: 1605–1631.[17] L.S.D. Morley. The triangular equilibrium element in the solution of plate bending problems. Aeronaut.,1968, 19: 149–169. / n = 4 n = 8 n = 16 n = 32 n = 64 orderrectangular meshes:Stokes 4.593E-1 0.201E-1 5.810E-2 2.223E-2 1.027E-2 1.111 4.616E-1 0.202E-1 5.827E-2 2.225E-2 1.027E-2 1.112 − − − − − − k p − p h k produced by V h × P h applied to the Brinkman problem through(4.2) over three kinds of meshes. [18] M. Neilan, D. Sap. Stokes elements on cubic meshes yielding divergence-free approximations. Calcolo, 2016,53(3): 263–283.[19] M. Neilan, D. Sap. Macro Stokes elements on quadrilaterals. Inter J. of Numer. Anal. Model., 2018, 15(4-5):729–745.[20] T. Nilssen, X.-C. Tai, R. Winther. A robust nonconforming H -element. Math. Comput., 2001, 70(234):489–505.[21] C. Park, D. Sheen. A quadrilateral Morley element for biharmonic equations. Numer. Math., 2013, 124:395–413.[22] L. R. Scott, S. Zhang. Finite element interpolation of nonsmooth functions satisfying boundary conditions.Math. Comput., 1990, 54(190): 483–493.[23] F. de Verbeke. A conforming finite element for plate bending. J. Solids Structure, 1968, 108: 4–95.[24] M. Wang, Z. Shi, J. Xu. A new class of Zienkiewicz-type non-conforming element in any dimensions. Numer.Math., 2007, 106(2): 335–347.[25] L. Wang, Y. Wu, X. Xie. Uniformly stable rectangular elements for fourth order elliptic singular perturbationproblems. Numer. Methods for Partial Differential Equations, 29(3): 721–737, 2013.[26] S. Zhang, X. Xie, Y. Chen. Low order nonconforming rectangular finite element methods for Dary-Stokesproblem. J. Comput. Math., 2009, 27: 400–424.[27] S. Zhang. Stable finite element pair for Stokes problem and discrete Stokes complex on quadrilateral grids.Numer. Math., 2016, 133: 371–408.[28] S. Zhang. On optimal finite element schemes for biharmonic equation. arXiv preprint arXiv:1805.03851, 2018.[29] X. Zhou, Z. Meng, Z. Luo. New nonconforming finite elements on arbitrary convex quadrilateral meshes. J.Comput. Appl. Math., 2016, 296: 798–814.[30] X. Zhou, Z. Meng, X. Fan, Z. Luo. Nonconforming polynomial mixed finite element for the Brinkman problemover quadrilateral meshes. Comput. Math. Appl., 2018, 76(4): 877–892.-element. Math. Comput., 2001, 70(234):489–505.[21] C. Park, D. Sheen. A quadrilateral Morley element for biharmonic equations. Numer. Math., 2013, 124:395–413.[22] L. R. Scott, S. Zhang. Finite element interpolation of nonsmooth functions satisfying boundary conditions.Math. Comput., 1990, 54(190): 483–493.[23] F. de Verbeke. A conforming finite element for plate bending. J. Solids Structure, 1968, 108: 4–95.[24] M. Wang, Z. Shi, J. Xu. A new class of Zienkiewicz-type non-conforming element in any dimensions. Numer.Math., 2007, 106(2): 335–347.[25] L. Wang, Y. Wu, X. Xie. Uniformly stable rectangular elements for fourth order elliptic singular perturbationproblems. Numer. Methods for Partial Differential Equations, 29(3): 721–737, 2013.[26] S. Zhang, X. Xie, Y. Chen. Low order nonconforming rectangular finite element methods for Dary-Stokesproblem. J. Comput. Math., 2009, 27: 400–424.[27] S. Zhang. Stable finite element pair for Stokes problem and discrete Stokes complex on quadrilateral grids.Numer. Math., 2016, 133: 371–408.[28] S. Zhang. On optimal finite element schemes for biharmonic equation. arXiv preprint arXiv:1805.03851, 2018.[29] X. Zhou, Z. Meng, Z. Luo. New nonconforming finite elements on arbitrary convex quadrilateral meshes. J.Comput. Appl. Math., 2016, 296: 798–814.[30] X. Zhou, Z. Meng, X. Fan, Z. Luo. Nonconforming polynomial mixed finite element for the Brinkman problemover quadrilateral meshes. Comput. Math. Appl., 2018, 76(4): 877–892.