High accuracy analysis of a nonconforming discrete Stokes complex over rectangular meshes

TL;DR

This paper develops a high-accuracy analysis of a nonconforming discrete Stokes complex on rectangular meshes, achieving superconvergence rates under certain conditions and verifying results with numerical tests.

Contribution

It introduces a simple, low-degree nonconforming element for biharmonic problems that attains higher convergence rates on uniform rectangular meshes, with rigorous proofs and numerical validation.

Findings

Achieves $O(h^2)$ convergence in discrete $H^2$-norm on uniform meshes.

Attains $O(h^3)$ convergence in discrete $H^1$-norm for convex domains.

Numerical tests confirm theoretical superconvergence results.

Abstract

This work is devoted to the high accuracy analysis of a discrete Stokes complex over rectangular meshes with a simple structure. The 0-form in the complex is a non $C^0$ nonconforming element space for biharmonic problems. This plate element contains only 12 degrees of freedom (DoFs) over a rectangular cell with a zeroth order weak continuity for the normal derivative, therefore only the lowest convergence order can be obtained by a standard consistency error analysis. Nevertheless, we prove that, if the rectangular mesh is uniform, an $O(h^2)$ convergence rate in discrete $H^2$-norm will be achieved. Moreover, based on a duality argument, it has an $O(h^3)$ convergence order in discrete $H^1$-norm if the solution region is convex. The 1-form and 2-form constitute a divergence-free pair for incompressible flow. We also show its higher accuracy than that derived from a usual error estimate under uniform partitions, which explains the phenomenon observed in our previous work. Numerical tests verify our theoretical results.

Tables

Table 1. Table 1 : The discrete H 2 superscript 𝐻 2 H^{2} , H 1 superscript 𝐻 1 H^{1} and L 2 superscript 𝐿 2 L^{2} errors produced by W h subscript 𝑊 ℎ W_{h} applied to the biharmonic problem determined by ( 5.1 ) for different n 𝑛 n .

$n$	${\left|u-u_h\right|}_{2,h}$	order	${\left|u-u_h\right|}_{1,h}$	order	${\Vert u-u_h\Vert }_0$	order
$4$	1.209E0		7.041E-2		9.209E-3
$8$	3.528E-1	1.78	9.173E-3	2.94	4.139E-4	4.48
$16$	8.880E-2	1.99	1.063E-3	3.11	2.060E-5	4.33
$32$	2.140E-2	2.05	1.242E-4	3.10	1.202E-6	4.10
$64$	5.198E-3	2.04	1.487E-5	3.06	7.428E-8	4.02

Table 2. Table 2 : The discrete H 2 superscript 𝐻 2 H^{2} , H 1 superscript 𝐻 1 H^{1} and L 2 superscript 𝐿 2 L^{2} errors produced by the Adini element for the biharmonic problem determined by ( 5.1 ) for different n 𝑛 n .

$n$	${\left|u-u_h\right|}_{2,h}$	order	${\left|u-u_h\right|}_{1,h}$	order	${\Vert u-u_h\Vert }_0$	order
$4$	1.112E0		1.270E-1		2.195E-2
$8$	3.113E-1	1.84	3.060E-2	2.05	6.283E-3	1.80
$16$	7.681E-2	2.02	7.642E-3	2.00	1.636E-3	1.94
$32$	1.901E-2	2.01	1.915E-3	2.00	4.137E-4	1.98
$64$	4.738E-3	2.00	4.791E-4	2.00	1.037E-4	2.00

Table 3. Table 3 : The errors for the velocity and the pressure produced by 𝐕 h × P h subscript 𝐕 ℎ subscript 𝑃 ℎ \bm{V}_{h}\times P_{h} for the Stokes problem determined by ( 5.2 ) for different n 𝑛 n .

$n$	${\left|𝒖-{𝒖}_h\right|}_{1,h}$	order	${\Vert 𝒖-{𝒖}_h\Vert }_0$	order	${\Vert p-p_h\Vert }_0$	order	${\Vert p-p_h^{\ast }\Vert }_0$	order
$4$	2.473E0		1.100E-1		6.569E-1		1.045E0
$8$	7.505E-1	1.72	1.657E-2	2.73	3.485E-1	0.91	3.187E-1	1.71
$16$	1.968E-1	1.93	2.305E-3	3.03	1.772E-1	0.98	6.173E-2	2.37
$32$	4.846E-2	2.02	2.399E-4	3.08	8.932E-2	0.99	9.938E-3	2.64
$64$	1.188E-2	2.03	2.878E-5	3.06	4.477E-2	1.00	1.869E-3	2.41

Equations

$$h_x=|E_{1,K}|=|E_{3,K}|=x_K^{\prime \prime }-x_K^{\prime },h_y=|E_{2,K}|=|E_{4,K}|=y_K^{\prime \prime }-y_K^{\prime }.$$

$${\tau }_{i,K}(w)$$

$${\tau }_{i+8,K}(w)$$

$${\mathbf{V}}_K=[P_1(K){]}^2\oplus \mathrm{span}\{\mathrm{curl}{x}^3,\mathrm{curl}{x}^{2}y,\mathrm{curl}x{y}^2,\mathrm{curl}{y}^3,\mathrm{curl}{x}^4,\mathrm{curl}{y}^4\};$$

$${\sigma }_{i,K}(\mathbf{v})$$

$${\sigma }_{i+8,K}(\mathbf{v})$$

$${\xi }_{1,K}$$

$${\xi }_{2,K}$$

$$\displaystyle W_{h}=\Bigg{\{}w$$

$$\displaystyle\mbox{vanishes at all $V\in\mathcal{V}_{h}^{b}$},~{}\int_{E}\left[w\right]_{E}\,\mathrm{d}s=\int_{E}\left[\frac{\partial w}{\partial\bm{n}_{E}}\right]_{E}\,\mathrm{d}s=0~{}\mbox{for all $E\in\mathcal{E}_{h}$}\Bigg{\}};$$

$$\displaystyle\bm{V}_{h}=\Bigg{\{}\bm{v}$$

$$\displaystyle\int_{E}[\bm{v}\cdot\bm{n}_{E}]_{E}\xi\,\mathrm{d}s=0,~{}\forall\xi\in P_{1}(E),~{}\int_{E}[\bm{v}\cdot\bm{t}_{E}]_{E}\,\mathrm{d}s=0,~{}\forall E\in\mathcal{E}_{h}\Bigg{\}};$$

$$\displaystyle P_{h}=\Big{\{}q$$

$$0$$

$$|v{|}_{m,h}={\left(\sum _{K\in {\mathcal{T}}_h}|v{|}_{m,K}^2\right)}^{1/2},m=0,1,2,3,$$

$$|w-{\mathcal{I}}_{h}w{|}_{j,h}\le C{h}^{k-j}|w{|}_k,k=3,4,j=0,1,2,3,\forall w\in H_0^2(\Omega )\cap H^4(\Omega ).$$

$${\left|\mathbf{v}-{\mathbf{\Pi }}_h\mathbf{v}\right|}_{1,h}={\left|\mathrm{curl}\varphi -{\mathrm{curl}}_h{\mathcal{I}}_h\varphi \right|}_{1,h}\le C|\varphi -{\mathcal{I}}_h\varphi {|}_{2,h}\le C{h}^2|\varphi {|}_4.$$

$${\Delta }^{2}u=f$$

$$u=\frac{\partial u}{\partial \mathbf{n}}=0$$

$$({\nabla }^{2}u,{\nabla }^{2}v)=(f,v),\forall v\in H_0^2(\Omega ),$$

$$\sum _{K\in {\mathcal{T}}_h}({\nabla }^2{u}_h,{\nabla }^2{v}_h{)}_K=(f,v_h),\forall v_h\in W_h.$$

$${\int }_E{\left[\frac{\partial w_h}{\partial {\mathbf{n}}_E}\right]}_E\mathrm{d}s=0,\forall E\in {\mathcal{E}}_h,\forall w_h\in W_h.$$

$$|u-u_h{|}_{2,h}\le C{h}^2|u{|}_4.$$

$$|u-u_h{|}_{2,h}\le C\left(\underset{v_h\in W_h}{\inf }|u-v_h{|}_{2,h}+\underset{w_h\in W_h}{\sup }\frac{|E_{2,h}(u,w_h)|}{|w_h{|}_{2,h}}\right)$$

$$E_{2,h}(u,w_h)=\sum _{K\in {\mathcal{T}}_h}({\nabla }^{2}u,{\nabla }^2{w}_h{)}_K-(f,w_h).$$

$$\underset{v_h\in W_h}{\inf }|u-v_h{|}_{2,h}\le |u-{\mathcal{I}}_{h}u{|}_{2,h}\le C{h}^2|u{|}_4,$$

$${\mathcal{I}}_K^{S_2}(w_h{|}_K)(V_{i,K})=(w_h{|}_K)(V_{i,K}),{\int }_{E_{i,K}}{\mathcal{I}}_K^{S_2}(w_h{|}_K)\mathrm{d}s={\int }_{E_{i,K}}{w}_h{|}_K\mathrm{d}s,i=1,2,3,4,$$

$$E_{2,h}(u,w_h)=$$

$$=$$

$$+\sum _{K\in {\mathcal{T}}_h}{\int }_{\partial K}\frac{{\partial }^{2}u}{\partial {\mathbf{n}}^2}\frac{\partial w_h}{\partial \mathbf{n}}\mathrm{d}s+\sum _{K\in {\mathcal{T}}_h}{\int }_{\partial K}\frac{{\partial }^{2}u}{\partial \mathbf{n}\partial \mathbf{t}}\frac{\partial w_h}{\partial \mathbf{t}}\mathrm{d}s$$

$$:=$$

$${\int }_K\nabla (w_h-{\mathcal{I}}_h^{S_2}{w}_h)\mathrm{d}\mathbf{x}={\int }_{\partial K}(w_h-{\mathcal{I}}_h^{S_2}{w}_h)\mathbf{n}\mathrm{d}s=0$$

$$|I_1(u,w_h)|$$

$$\{\begin{array}{rl}& \le \sum _{K\in {\mathcal{T}}_h}Ch|\nabla \Delta u{|}_{1,K}Ch|w_h{|}_{2,K}\le C{h}^2|u{|}_4|w_h{|}_{2,h};\\ & \le \sum _{K\in {\mathcal{T}}_h}Ch|\nabla \Delta u{|}_{1,K}C{h}^2|w_h{|}_{3,K}\le C{h}^3|u{|}_4|w_h{|}_{3,h}.\end{array}$$

$$|I_2(u,w_h)|\le \Vert f{\Vert }_0\Vert w_h-{\Pi }_h{w}_h{\Vert }_0\{\begin{array}{rl}& \le C{h}^2\Vert f{\Vert }_0|w_h{|}_{2,h}\le C{h}^2|u{|}_4|w_h{|}_{2,h};\\ & \le C{h}^3\Vert f{\Vert }_0|w_h{|}_{3,h}\le C{h}^3|u{|}_4|w_h{|}_{3,h}.\end{array}$$

$${\int }_E{\left[\frac{\partial w_h}{\partial {\mathbf{t}}_E}\right]}_E{\xi }_E\mathrm{d}s=0,\forall {\xi }_E\in P_1(E),\forall E\in {\mathcal{E}}_h,\forall w_h\in W_h,$$

Full text

Keywords: High accuracy analysis; discrete Stokes complex; uniform rectangular meshes; dual error estimate*

1 Introduction

It is generally not an easy task to construct finite elements for biharmonic problems with high accuracy but a simple structure. Although there are several conforming elements achieving high order convergence rate, such as the Argyris element **[1]**, the Bogner-Fox-Schmit (BFS) element **[3]**, the Hsieh-Clough-Tocher (HCT) element **[6]** and the Fraijes de Veubeke-Sander (FVS) element **[5, 28]**, they suffer from the drawbacks caused by the strong $C^{1}$ continuity requirement. For example, many degrees of freedom (DoFs) have to be utilized, the degrees of the polynomials in shape function spaces must be very high, or the structures of shape function spaces are too complicated. An alternative way is to adopt nonconforming elements. A successful construction is the Adini element on rectangular meshes, whose DoFs are of the nodal type at the vertices of a mesh, resulting in fewer total DoFs. In **[20]**, it was shown that a second order convergence rate can be guaranteed in energy norm if the rectangular cells are of the same size, which is based on the interior symmetry of a rectangle. Thereafter, Luo and Lin **[21]** pointed out that the result above is still valid without such a uniformity assumption. This fact was recently extended to higher-dimensional cases in **[17]**. However, one can only expect the same convergence order in discrete $H^{1}$- and $L^{2}$-norms due to the lower bound estimate given by Hu et al. **[16, 17]**. In recent years, the design for nonconforming elements draw rapidly increasing attention. The bubble function technique has become a standard tool to achieve a high convergence order if the desired element is $C^{0}$-continuous **[12, 29, 9, 10]**. There are also some completely nonconforming constructions utilizing a high order inter-element continuity **[30]**. In addition, we remark that the purpose of high accuracy can also be achieved by nonstandard methods such as the double set parameter (DSP) method **[27]** or superconvergence analysis and postprocessing for low order elements **[22, 18]**.

Biharmonic elements are often used to design divergence-free Stokes elements with discrete Stokes complexes. A divergence-free Stokes element is usually preferred than a non divergence-free one, as the former benefits from many advantages, say, the conservation law is preserved numerically for incompressible flow, the discrete scheme seems more robust and accurate with respect to the time discretization. A comprehensive review of this topic can be found in **[19]**. Conforming examples include the Stokes elements and the discrete Stokes complexes derived from the Argyris element **[11]**, the BFS element **[23]**, the HCT element **[4]** and the FVS element **[24]**, etc. High order nonconforming complexes such as in **[15]** and **[30]** are also successful constructions. Recently, the de Rham complex from the Adini element was also shown by Gillette et al. **[13]** as a member of a nonstandard family.

In 1996, in order to construct a family of biharmonic elements by the DSP method, Chen and Shi **[7]** designed an intermediate rectangular element of 12 local DoFs, which are the four point values, the four edge integrals of the shape functions and the edge integrals of their normal derivatives. The shape function space is selected as $P_{3}$ plus a space spanned by two monomials of degree four. It was then used as a non $C^{0}$ nonconforming element for fourth order elliptic singular perturbation problems with a uniform convergence property **[8]**. Very recently, we extended this method to general convex quadrilateral meshes with some slight modifications to the shape function space, and induced a divergence-free Stokes element for the Brinkman problem **[31]**. The associated discrete Stokes complex was also provided. Based on a standard consistency error estimate, only the lowest order convergence rates of this complex can be derived. Surprisingly, in our numerical tests, we found that the velocity has an $O(h^{2})$ approximation order in energy norm if the mesh consists of uniform rectangles, one order higher than that given by the theoretical analysis therein. Moreover, as far as we are aware, although there have been many researches on the DSP method motivated by **[7]**, the high accuracy analysis for the 12-DoF intermediate rectangular element introduced in **[7]** is still unknown.

The aim of this work is to give a high accuracy analysis of the intermediate rectangular element given in **[7]** and the induced discrete Stokes complex. We prove, by the special properties of the shape function space, that if the rectangular mesh is uniform, this plate element has actually an $O(h^{2})$ convergence rate in discrete $H^{2}$-norm if it works for the biharmonic problem. In addition, with the aid of an auxiliary biharmonic element and through a duality argument, an $O(h^{3})$ convergence order in discrete $H^{1}$-norm can be achieved, provided that the solution region is convex. In comparison with the well-known Adini element, the numbers of local DoFs are the same. Although there are more global DoFs, this element space is highly nonconforming with many edge-oriented basis functions, enjoying a low brand width when the stiffness matrix is accumulated. This might be more convenient in implementation, especially for parallel computing. Moreover, the convergence rate in discrete $H^{1}$-norm is one order faster than the Adini counterpart. From this element, a nonconforming discrete complex is presented using the strategy in **[31]**. The 1-form and 2-form constitute a divergence-free Stokes pair for incompressible flow. Under uniform partitions, we also show its higher accuracy than that derived from a usual error estimate, in both energy norm and $L^{2}$-norm for the velocity. The convergence order of the pressure can be improved by a simple postprocessing. This complex has the same key features as that provided in our previous work **[31]** over uniform rectangular meshes, and so a rigorous explanation for the phenomenon observed in **[31]** is obtained.

The rest of this work is arranged as follows. In Section 2, we give an overview of the investigated discrete Stokes complex. Section 3 provides the high accuracy analysis of the 0-form for biharmonic problems. The error estimate for the Stokes pair determined by the 1-form and 2-form will be presented in Section 4. Finally, numerical examples are given in Section 5 to verify our theoretical analysis.

Standard notations in Sobolev spaces are used throughout this work. For a domain $D\subset\mathbb{R}^{2}$, $\bm{n}$ and $\bm{t}$ will be the unit outward normal and tangent vectors on $\partial D$, respectively. The notation $P_{k}(D)$ denotes the usual polynomial space over $D$ of degree no more than $k$. The norms and semi-norms of order $m$ in the Sobolev spaces $H^{m}(D)$ are indicated by $\|\cdot\|_{m,D}$ and $|\cdot|_{m,D}$, respectively. The space $H_{0}^{m}(D)$ is the closure in $H^{m}(D)$ of $C_{0}^{\infty}(D)$. The relation between $\bm{H}_{0}(\mathrm{div};\Omega)$ and $\bm{H}(\mathrm{div};\Omega)$ is in a similar manner. We also adopt the convention that $L^{2}(D):=H^{0}(D)$, where the inner-product is denoted by $(\cdot,\cdot)_{D}$. The functions in its subspace $L_{0}^{2}(D)$ are of zero integral. These notations of norms, semi-norms and inner-products also work for vector- and matrix-valued Sobolev spaces, where the subscript $\Omega$ will be omitted if the domain $D=\Omega$. Moreover, $C$ is a positive constant independent of the mesh size $h$ and may be different in different places.

2 Preliminaries: A nonconforming discrete Stokes complex

Let $\Omega\subset\mathbb{R}^{2}$ be a polygonal domain and $\partial\Omega$ be its boundary. We assume that $\Omega$ can be uniformly partitioned into rectangular cells, denoted by $\mathcal{T}_{h}$, with $h$ being the length of the diagonal of each rectangle. For a cell $K\in\mathcal{T}_{h}$, the vertices are designated as $V_{i,K}$, $i=1,2,3,4$. We select the coordinate system fulfilling that $V_{1,K}=(x^{\prime}_{K},y^{\prime}_{K})$, $V_{2,K}=(x^{\prime\prime}_{K},y^{\prime}_{K})$, $V_{3,K}=(x^{\prime\prime}_{K},y^{\prime\prime}_{K})$, $V_{4,K}=(x^{\prime}_{K},y^{\prime\prime}_{K})$. The edges $E_{i,K}=V_{i,K}V_{i+1,K}$ are parallel to the coordinate axes, whose lengths are naturally given by

[TABLE]

Note that here $i$ is taken modulo four. We can now give the following two elements over $K$. See Figure 1 as an example.

Definition 2.1**.**

(See also [7]) The element $(K,W_{K},T_{K})$ is defined through

• •

$K\in\mathcal{T}_{h}$* is a rectangle;*

• •

$W_{K}=P_{3}(K)\oplus\mathrm{span}\{x^{4},y^{4}\}$* is the shape function space;*

• •

$T_{K}=\{\tau_{i,K},~{}i=1,2,\ldots,12\}$* is the DoF set where*

[TABLE]

Definition 2.2**.**

The element $(K,\bm{V}_{K},\Sigma_{K})$ is defined through

• •

$K\in\mathcal{T}_{h}$* is a rectangle;*

• •

$\bm{V}_{K}$* is the shape function space, where*

[TABLE]

• •

$\Sigma_{K}=\{\sigma_{i,K},~{}i=1,2,\ldots,12\}$* is the DoF set:*

[TABLE]

where $\xi_{i,K}\in P_{1}(K)$ are monomials in variable $x$ or $y$ such that

[TABLE]

Theorem 2.3**.**

These two elements are well-defined.

Proof.

The unisolvency of $W_{K}$ with respect to $T_{K}$ can be found in [7]. The assertion for the second element can be directly obtained by the technique provided in Theorem 2.5 in [31] using the first element. ∎

For a given partition $\mathcal{T}_{h}$, the sets of all vertices, interior vertices, boundary vertices, edges, interior edges and boundary edges are correspondingly denoted by $\mathcal{V}_{h}$, $\mathcal{V}_{h}^{i}$, $\mathcal{V}_{h}^{b}$, $\mathcal{E}_{h}$, $\mathcal{E}_{h}^{i}$ and $\mathcal{E}_{h}^{b}$. For each $E\in\mathcal{E}_{h}$, $\bm{n}_{E}$ is a fixed unit vector perpendicular to $E$ and $\bm{t}_{E}$ is a vector obtained by rotating $\bm{n}_{E}$ by ninety degree counterclockwisely. Moreover, for $E\in\mathcal{E}_{h}^{i}$, the jump of a function $v$ across $E$ is defined as $[v]_{E}=v|_{K_{1}}-v|_{K_{2}}$, where $K_{1}$ and $K_{2}$ are the cells sharing $E$ as a common edge, and $\bm{n}_{E}$ points from $K_{1}$ to $K_{2}$. For $E\in\mathcal{E}_{h}^{b}$, we set $[v]_{E}=v|_{K}$ if $E$ is an edge of $K$. Under these notations, we introduce the following finite element spaces over $\Omega$:

[TABLE]

For each $K\in\mathcal{T}_{h}$, we define interpolation operators $\mathcal{I}_{K}:H^{2}(K)\rightarrow W_{K}$ and $\bm{\Pi}_{K}:[H^{1}(K)]^{2}\rightarrow\bm{V}_{K}$, such that $\tau_{i,K}(\mathcal{I}_{K}w)=\tau_{i,K}(w)$, $\sigma_{i,K}(\bm{\Pi}_{K}\bm{v})=\sigma_{i,K}(\bm{v})$, $i=1,2,\ldots,12$. The global interpolation operators $\mathcal{I}_{h}:H_{0}^{2}(\Omega)\rightarrow W_{h}$ and $\bm{\Pi}_{h}:[H_{0}^{1}(\Omega)]^{2}\rightarrow\bm{V}_{h}$ are naturally set as $\mathcal{I}_{h}|_{K}=\mathcal{I}_{K}$ and $\bm{\Pi}_{h}|_{K}=\bm{\Pi}_{K}$, $\forall K\in\mathcal{T}_{h}$. Similarly, we write $\mathcal{P}_{h}$ as the $L^{2}$-projection operator from $L^{2}(\Omega)$ to $P_{h}$. Moreover, the differential operators $\bm{\mathrm{curl}}$ and $\mathrm{div}$ have their discrete versions $\bm{\mathrm{curl}}_{h}$ and $\mathrm{div}_{h}$, respectively, determined by $\bm{\mathrm{curl}}_{h}|_{K}=\bm{\mathrm{curl}}$ on $K$, and $\mathrm{div}_{h}|_{K}=\mathrm{div}$ on $K$, $\forall K\in\mathcal{T}_{h}$. The following commutative diagram provides a discrete Stokes complex.

Theorem 2.4**.**

It holds that

[TABLE]

Moreover, the sequence in the second row is exact.

Proof.

Please see Theorem 4.1 in [31] for quite a similar argument. ∎

Define the semi-norms

[TABLE]

then $|\cdot|_{2,h}$ and $|\cdot|_{1,h}$ are norms over $W_{h}$ and $\bm{V}_{h}$, respectively. Since $\mathcal{I}_{K}$ preserves $P_{3}(K)$ for each $K\in\mathcal{T}_{h}$, the Bramble-Hilbert lemma gives

[TABLE]

Moreover, if $\bm{v}\in[H_{0}^{1}(\Omega)\cap H^{3}(\Omega)]^{2}$ and $\mathrm{div}\,\bm{v}=0$, we can find $\bm{v}=\bm{\mathrm{curl}}\,\phi$ for some $\phi\in H_{0}^{2}(\Omega)\cap H^{4}(\Omega)$. The commutative diagram (2.1) implies

[TABLE]

For a cell $K\in\mathcal{T}_{h}$, we write $\mathcal{P}_{k}^{K}$ as the $L^{2}$-projection operator from $L^{2}(K)$ to $P_{k}(K)$. If $E$ is an edge of $K$, we also define $\mathcal{P}_{k}^{E}$ that $L^{2}$-projects the traces of the functions in $H^{1}(K)$ on $E$ to $P_{k}(E)$.

3 High accuracy analysis of $W_{h}$ for biharmonic problems

Let us turn to the application of $W_{h}$ to the biharmonic problem. Given $f\in H^{-2}(\Omega)$, the biharmonic problem is to find $u$ such that

[TABLE]

A weak form can be represented as: Find $u\in H_{0}^{2}(\Omega)$ such that

[TABLE]

where $\nabla^{2}v$ is the Hessian matrix of $v$. Note that we have fixed the Poisson ratio $\sigma=0$. For other $\sigma\in(0,1/2]$, the coming analysis has no intrinsic difference. The discrete form of (3.2) reads: Find $u_{h}\in W_{h}$ such that

[TABLE]

This problem has a unique solution due to the Lax-Milgram theorem. Note that

[TABLE]

Hence, based on a classical consistency error analysis, one can only predict an $O(h)$ convergence order in discrete $H^{2}$-norm. Nevertheless, the result below hints a higher accuracy.

Theorem 3.1**.**

Let $u\in H_{0}^{2}(\Omega)\cap H^{4}(\Omega)$ and $u_{h}\in W_{h}$ be the solutions of (3.2) and (3.3), respectively. If $\mathcal{T}_{h}$ is uniform, then

[TABLE]

Proof.

Let us begin with the Strang lemma

[TABLE]

with

[TABLE]

It follows from (2.2) that

[TABLE]

therefore it suffices to consider the consistency error $E_{2,h}(u,w_{h})$.

If $u\in H^{4}(\Omega)$ then $f\in L^{2}(\Omega)$. Let $\mathcal{I}^{S_{2}}_{K}$ be the nodal interpolation operator from $W_{K}$ to $S_{2}(K)$, the second order serendipity element space $P_{2}(K)\oplus\mathrm{span}\{x^{2}y,xy^{2}\}$, satisfying

[TABLE]

and then set $\mathcal{I}^{S_{2}}_{h}$ via $\mathcal{I}^{S_{2}}_{h}|_{K}=\mathcal{I}^{S_{2}}_{K}$, $\forall K\in\mathcal{T}_{h}$. Owing to the weak continuity of $w_{h}$, we find $\mathcal{I}^{S_{2}}_{h}w_{h}\in H_{0}^{1}(\Omega)$. An integration by parts using (3.1) gives

[TABLE]

As far as $I_{1}(u,w_{h})$ is concerned, noting that

[TABLE]

and $\mathcal{I}^{S_{2}}_{K}$ preserves $P_{2}(K)$, we derive

[TABLE]

One can also estimate $|I_{2}(u,w_{h})|$ as

[TABLE]

Moreover, since the vertex values and the edge integrals of $w_{h}$ are continuous, integrating by parts results in

[TABLE]

which asserts that

[TABLE]

Hence, it remains to estimate $|I_{3}(u,w_{h})|$.

To this end, notice from the weak continuity of $w_{h}$ that

[TABLE]

Let us consider the summation involving $I_{1,K}(u,w_{h})$. Clearly,

[TABLE]

On the other hand, introduce the bilinear form $I^{\prime}_{1,K}(u,w_{h})$:

[TABLE]

It follows from the inverse inequality that

[TABLE]

Moreover, we find

[TABLE]

Next, an easy verification shows

[TABLE]

therefore by the Simpson quadrature rule,

[TABLE]

where $i^{\prime\prime}=3$ and $i^{\prime}=4$ for $i=3$, and $i^{\prime\prime}=2$ and $i^{\prime}=1$ for $i=1$. This gives

[TABLE]

and thus

[TABLE]

It follows from (3.16) and (3.19) that

[TABLE]

Owing to (3.13), (3.15), (3.20) and the Bramble-Hilbert lemma, we derive

[TABLE]

Substituting (3.21) into (3.12) results in

[TABLE]

Now we are in the position to estimate the right-hand side in (3.22). Indeed,

[TABLE]

where we have used the integration by parts formula, the weak continuity of $w_{h}$ and the uniformity of the partition $\mathcal{T}_{h}$. For each $K\in\mathcal{T}_{h}$,

[TABLE]

Note that $x^{3}$ and $x^{4}$ are independent of $y$, therefore a closer observation for $\mathcal{I}^{S_{2}}_{h}$ gives

[TABLE]

Hence, it follows from a standard argument and the inverse inequality that

[TABLE]

Combining this estimate with (3.22) we get

[TABLE]

Similarly, it holds that

[TABLE]

which along with (3.8), (3.9), (3.11) and (3.12) gives

[TABLE]

and the proof is completed. ∎

Thanks to the zeroth order weak continuity for the gradients of the functions in $W_{h}$, we can obtain the discrete $H^{1}$ error estimate using the duality argument provided in **[26, 25]**. To proceed, we need to construct an auxiliary $C^{0}$ nonconforming element with the first order orthogonality for normal derivatives, which will be carried out in two steps. In the first step, define the local shape function space

[TABLE]

where $b_{K}$ is a 4-degree bubble polynomial fulfilling $b_{K}|_{\partial K}=0$. The functions $\phi_{i,K}$ are selected such that $\widetilde{W}_{K}^{*}$ is unisolvent to the DoF set

[TABLE]

According to **[9]**, a successful selection of these $\phi_{i,\widehat{K}}$ over the reference square $\widehat{K}=[-1,1]^{2}$ can be given as

[TABLE]

then the physical $\phi_{i,K}$ will be obtained by an affine equivalent technique, $i=1,2,\ldots,8$. The second step is to take a subspace of $\widetilde{W}_{K}^{*}$, named as $W_{K}^{*}$, satisfying the following four linear relations of the DoFs in $\widetilde{T}_{K}^{*}$:

[TABLE]

where $i^{\prime\prime}=3$ and $i^{\prime}=4$ for $i=3$, and $i^{\prime\prime}=2$ and $i^{\prime}=1$ for $i=1$; $j^{\prime\prime}=3$ and $j^{\prime}=2$ for $j=2$, and $j^{\prime\prime}=4$ and $j^{\prime}=1$ for $j=4$. By the Simpson quadrature rule, we find $P_{3}(K)\subset W_{K}^{*}$. The actual DoF set is then set as

[TABLE]

and the global finite element space $W_{h}^{*}$ reads

[TABLE]

Clearly, $W_{h}^{*}\subset H_{0}^{1}(\Omega)$. For a fixed $K$, the interpolation $\mathcal{I}_{K}^{*}$ from $H^{3}(K)$ to $W_{K}^{*}$ is defined such that $\tau(\mathcal{I}_{K}^{*}v)=\tau(v)$, $\forall v\in H^{3}(K)$, $\forall\tau\in T_{K}^{*}$. The global version $\mathcal{I}_{h}^{*}$ from $H^{3}(\Omega)\cap H_{0}^{2}(\Omega)$ to $W_{h}^{*}$ is naturally determined by $\mathcal{I}_{h}^{*}|_{K}=\mathcal{I}_{K}^{*}$, $\forall K\in\mathcal{T}_{h}$. This element is indeed a rectangular extension of the triangular $C^{0}$ nonconforming element designed by Gao et al. **[12]**.

Theorem 3.2**.**

Let the solution domain $\Omega$ be convex. In addition, under the same assumptions in Theorem 3.1, we have

[TABLE]

Proof.

We first set $e_{h}=u-u_{h}$ and consider the following dual problem: Find $\phi\in H_{0}^{2}(\Omega)\cap H^{3}(\Omega)$, such that

[TABLE]

For the convex domain $\Omega$, it follows from [14] that $\Delta^{2}$ is an isomorphism from $H_{0}^{2}(\Omega)\cap H^{3}(\Omega)$ onto $H^{-1}(\Omega)$, and therefore (3.26) has a unique solution with the bound

[TABLE]

We will use the decomposition

[TABLE]

The first summation is estimated with the aid of (3.5) as

[TABLE]

As far as the second summation in (3.28) is concerned, note that $\mathcal{I}^{S_{2}}_{h}e_{h}\in H_{0}^{1}(\Omega)$, and a density argument using (3.26) with an integration by parts shows that

[TABLE]

where we have used the weak continuity of $e_{h}$, see e.g. (3.4) and (3.10). The first two terms in the last line of (3.30) are bounded by

[TABLE]

according to (3.5) and (3.27). Hence it remains to estimate the last term in the last line of (3.30).

To this end, we select $\phi^{*}=\mathcal{I}_{h}^{*}\phi\in W_{h}^{*}$, and then $\mathcal{I}_{h}\phi^{*}\in W_{h}$ due to the definitions of $W_{h}^{*}$ and $W_{h}$. As a result, it follows from (3.3) and (3.7) that

[TABLE]

Invoking (3.8), (3.9) and (3.11), one sees

[TABLE]

Furthermore, by the triangle inequality,

[TABLE]

Hence, we shall estimate $I_{3}(u,\mathcal{I}_{h}\phi^{*})$. Indeed, we write $I_{3}(u,\mathcal{I}_{h}\phi^{*})$ as

[TABLE]

Note that the linear relations in (3.25) and the weak continuity of $\phi^{*}$ ensure that

[TABLE]

which gives

[TABLE]

Next we consider $I_{3}(u,\mathcal{I}_{h}\phi^{*}-\phi^{*})$. Just as in (3.12), by the weak continuity of both $W_{h}$ and $W_{h}^{*}$, we can write

[TABLE]

Let us investigate the relation between $I_{1,K}(u,\phi^{*})$ and $I^{\prime}_{1,K}(u,\phi^{*})$ defined in (3.12) and (3.14) for each $K\in\mathcal{T}_{h}$. It is clear that (3.16) still holds if $w_{h}$ is replaced by $\phi^{*}$. Moreover, a careful comparison between (3.17) and the first line in (3.25) shows that they are precisely the same thing if $w_{h}$ in (3.17) is replaced by $v$ in (3.25). Hence, one can follow the same argument as in (3.17) and (3.18) to derive an analogue of (3.19). As a consequence,

[TABLE]

which asserts by (3.13), (3.15), (3.20) and the Bramble-Hilbert lemma that

[TABLE]

Substituting (3.37) into (3.36) results in

[TABLE]

The right-hand side of (3.38) can be estimated by an integration by parts:

[TABLE]

The last summation immediately vanishes as $\phi^{*}\in H_{0}^{1}(\Omega)$. Therefore using the argument in (3.24), one derives

[TABLE]

which along with (3.38) gives

[TABLE]

In a similar fashion, it holds that

[TABLE]

Substituting (3.39), (3.40) into (3.36) and combining (3.32), (3.33) and (3.35) with (3.31), we have

[TABLE]

Note that we have again used the regularity condition (3.27). Finally, it follows from (3.28), (3.29) and (3.30) that

[TABLE]

The desired conclusion is obtained by dividing both sides by $|e_{h}|_{1,h}$. ∎

Remark 3.3**.**

In our previous work [31] over general quadrilateral meshes, the shape function space $W_{K}$ is of degree six, so that the edge integrals can be replaced by edge midpoint values in the DoF set, which seems cheaper in practical applications, especially for nonhomogeneous boundary conditions. Nevertheless, when applied to uniform rectangular meshes, the key features (3.20) and (3.23) remain valid, and so the error estimates in Theorems 3.1 and 3.2 still hold.

4 High accuracy analysis of $\bm{V}_{h}\times P_{h}$ for Stokes problems

In a similar fashion, we will give the high accuracy analysis of the Stokes element $\bm{V}_{h}\times P_{h}$. The Stokes problem for incompressible flow reads: For a given force field $\bm{f}$, find the velocity $\bm{u}$ and the pressure $p$ satisfying

[TABLE]

A weak formulation of (4.1) is to find a pair $(\bm{u},p)\in[H_{0}^{1}(\Omega)]^{2}\times L_{0}^{2}(\Omega)$ such that

[TABLE]

We will approximate (4.2) by seeking $(\bm{u}_{h},p_{h})\in\bm{V}_{h}\times P_{h}$ yielding

[TABLE]

with bilinear forms

[TABLE]

Following a classical consistency error analysis, only an $O(h)$ convergence order can be derived for $|\bm{u}-\bm{u}_{h}|_{1,h}$, since one only has the following weak continuity for the tangential component:

[TABLE]

As before, we have the following higher order error estimate.

Theorem 4.1**.**

The discrete problem (4.3) has a unique solution $(\bm{u}_{h},p_{h})\in\bm{V}_{h}\times P_{h}$. Let $(\bm{u},p)\in\left([H_{0}^{1}(\Omega)\cap H^{3}(\Omega)]^{2}\right)\times(L_{0}^{2}(\Omega)\cap H^{2}(\Omega))$ be the weak solution of (4.2) fulfilling $\bm{u}=\bm{\mathrm{curl}}\,\phi$, $\phi\in H_{0}^{2}(\Omega)\cap H^{4}(\Omega)$. If $\mathcal{T}_{h}$ is uniform, the following estimates hold:

[TABLE]

Proof.

From the commutative diagram (2.1) we see $\mathrm{div}_{h}\bm{V}_{h}\subset P_{h}$ and

[TABLE]

Hence, using Fortin’s trick and the continuous inf-sup condition, we obtain the following discrete version

[TABLE]

By a standard argument of the mixed finite element method, e.g., Theorem 3.1 in [31] or Theorem 5.2.6 in [2], (4.3) has a unique solution $(\bm{u}_{h},p_{h})\in\bm{V}_{h}\times P_{h}$. Moreover,

[TABLE]

where

[TABLE]

and the consistency term $E_{1,h}(\bm{u},p,\bm{w}_{h})$ is given by

[TABLE]

It follows from (2.3) that

[TABLE]

as $\bm{\Pi}_{h}\bm{u}\in\bm{Z}_{h}$. We shall therefore consider the consistency error. To this end, we rewrite $E_{1,h}(\bm{u},p,\bm{w}_{h})$ as

[TABLE]

In fact, the weak continuity of $\bm{w}_{h}$ implies

[TABLE]

Using a similar strategy as in (3.11) one sees

[TABLE]

thus it remains to estimate $J_{2}(\bm{u},\bm{w}_{h})$.

If we write $\bm{u}=(u_{1},u_{2})^{T}$ and $\bm{w}_{h}=(w_{h1},w_{h2})^{T}$, notice from the weak continuity of $\bm{w}_{h}$ that

[TABLE]

Let us again focus only on $J_{1,K}(u_{1},w_{h1})$. On one hand,

[TABLE]

On the other hand, define

[TABLE]

then

[TABLE]

Moreover, it follows from the definition of $\bm{V}_{K}$ that $w_{h1}|_{K}\in P_{2}(K)\oplus\mathrm{span}\{y^{3}\}$, and hence,

[TABLE]

The Simpson quadrature rule ensures

[TABLE]

where $i^{\prime\prime}=3$ and $i^{\prime}=4$ for $i=3$, and $i^{\prime\prime}=2$ and $i^{\prime}=1$ for $i=1$. An argument just like in (3.17) and (3.18) leads to

[TABLE]

Collecting (4.11), (4.12) and (4.14) with the Bramble-Hilbert lemma gives

[TABLE]

Substituting (4.15) into (4.10) results in

[TABLE]

Let us turn to the right-hand side in (4.16). We need to define the following interpolation operator $\Pi_{E_{i,K}}$ from $P_{2}(K)\oplus\mathrm{span}\{y^{3}\}$ to $P_{1}(E_{i,K})$ by setting

[TABLE]

A simple calculation shows that this interpolation is well-defined. Moreover,

[TABLE]

Hence, according to the weak continuity of the normal components of $\bm{w}_{h}$ over $E_{2,K}$ and $E_{4,K}$, it holds that

[TABLE]

Hence,

[TABLE]

and therefore by (4.16),

[TABLE]

An almost identical argument applies to $J_{2,K}(u_{2},w_{h2})$ along with (4.10) gives

[TABLE]

and then

[TABLE]

which establishes the error estimate for $|\bm{u}-\bm{u}_{h}|_{1,h}$ due to (4.6). Again owing to (4.6), the estimate for the pressure is simply based on the fact above and the approximation order $O(h)$ of the projection $\mathcal{P}_{h}$. ∎

From Theorem 4.1, the velocity error has a higher accuracy of order $O(h^{2})$, while the pressure can only achieve an $O(h)$ convergence rate, as the approximation error is only $O(h)$. Nevertheless, this can be improved by a simple postprocessing based on the idea from **[11]** or **[24]**. For any $K\in\mathcal{T}_{h}$, we define $p_{K}^{*}\in P_{1}(K)$ by setting

[TABLE]

Then the postprocessed solution $p_{h}^{*}$ is discontinuous and piecewise linear, defined via $p_{h}^{*}|_{K}=p_{K}^{*}$, $\forall K\in\mathcal{T}_{h}$. The following result can be derived in a similar manner as Theorem 4.4 in **[24]**, and so the proof is omitted.

Theorem 4.2**.**

Under the same assumptions in Theorem 4.1, we have

[TABLE]

In what follows, we shall provide the $L^{2}$-error estimate for the velocity by a duality argument. To begin with, an $H(\mathrm{div})$-conforming but $H^{1}$-nonconforming element will be constructed using $W_{h}^{*}$. Over each $K\in\mathcal{T}_{h}$, the shape function space $\bm{V}_{K}^{*}$ is defined by

[TABLE]

and the DoF set is given as

[TABLE]

Using the technique provided in Theorem 2.5 in **[31]** and the unisolvency of $W_{K}^{*}$ with respect to $T_{K}^{*}$, the element $(K,\bm{V}_{K}^{*},\Sigma_{K}^{*})$ is well-defined. By the Simpson quadrature rule applied to $(x,y)^{T}$ and (3.25), the following relations hold:

[TABLE]

with $i^{\prime\prime}=3$ and $i^{\prime}=4$ for $i=3$, and $i^{\prime\prime}=2$ and $i^{\prime}=1$ for $i=1$; $j^{\prime\prime}=3$ and $j^{\prime}=2$ for $j=2$, and $j^{\prime\prime}=4$ and $j^{\prime}=1$ for $j=4$. The global finite element space $\bm{V}_{h}^{*}$ is set as

[TABLE]

Clearly, $\bm{V}_{h}^{*}\in\bm{H}_{0}(\mathrm{div};\Omega)$. For a fixed $K$, the interpolation $\bm{\Pi}_{K}^{*}$ from $[H^{2}(K)]^{2}$ to $\bm{V}_{K}^{*}$ is defined such that $\sigma(\bm{\Pi}_{K}^{*}\bm{v})=\sigma(\bm{v})$, $\forall\bm{v}\in[H^{2}(K)]^{2}$, $\forall\sigma\in\Sigma_{K}^{*}$. The global version $\bm{\Pi}_{h}^{*}$ from $[H^{2}(\Omega)\cap H_{0}^{1}(\Omega)]^{2}$ to $\bm{V}_{h}^{*}$ is naturally determined by $\bm{\Pi}_{h}^{*}|_{K}=\bm{\Pi}_{K}^{*}$, $\forall K\in\mathcal{T}_{h}$.

Remark 4.3**.**

Using a standard argument (see e.g. Theorem 4.1 in [31]), we have constructed the following discrete Stokes complex:

[TABLE]

Theorem 4.4**.**

Let the solution domain $\Omega$ be convex. In addition, under the same assumptions in Theorem 4.1, we have

[TABLE]

Proof.

Set $\bm{e}_{h}=\bm{u}-\bm{u}_{h}$ and consider the following dual problem: Find $\bm{\psi}$ and $\chi$ satisfying

[TABLE]

with the weak form

[TABLE]

where $(\bm{\psi},\chi)\in[H_{0}^{1}(\Omega)]^{2}\times L_{0}^{2}(\Omega)$. Owing to the convexity of $\Omega$, it follows from [14] that $(\bm{\psi},\chi)\in[H^{2}(\Omega)]^{2}\times H^{1}(\Omega)$ yielding the regularity

[TABLE]

Multiplying $\bm{e}_{h}$ on both sides of (4.20) and integrating by parts show that

[TABLE]

On the other hand, we define $\bm{\psi}^{*}=\bm{\Pi}_{h}^{*}\bm{\psi}\in\bm{V}_{h}^{*}$ and then $\bm{\Pi}_{h}\bm{\psi}^{*}\in\bm{V}_{h}$ in light of the definitions of $\bm{V}_{h}^{*}$ and $\bm{V}_{h}$ and the fact that $\bm{\psi}^{*}\in\bm{H}_{0}(\mathrm{div};\Omega)$. Moreover, on each $K$, by Green’s formula,

[TABLE]

therefore $\bm{\Pi}_{h}\bm{\psi}^{*}\in\bm{Z}_{h}$. Take $\bm{v}_{h}=\bm{\Pi}_{h}\bm{\psi}^{*}$ in (4.3), multiply $\bm{\Pi}_{h}\bm{\psi}^{*}$ on both sides of (4.1) and integrate by parts, then

[TABLE]

The difference of (4.22) and (4.23) gives

[TABLE]

The first term is bounded by the triangle inequality, the Bramble-Hilbert lemma and (4.5):

[TABLE]

The second term is estimated using the lowest order inter-element orthogonality (4.4) and (4.8):

[TABLE]

As far as the last term is concerned, we see from (4.9) that

[TABLE]

Hence, it suffices to investigate $J_{2}(\bm{u},\bm{\Pi}_{h}\bm{\psi}^{*})$.

To proceed, we decompose $J_{2}(\bm{u},\bm{\Pi}_{h}\bm{\psi}^{*})$ by

[TABLE]

It follows from (4.19) and the continuity of $\bm{\psi}^{*}$ at the two endpoints of each edge that

[TABLE]

therefore

[TABLE]

If we write $\bm{\Pi}_{h}\bm{\psi}^{*}-\bm{\psi}^{*}=(\psi_{h1}^{-},\psi_{h2}^{-})^{T}$ for short, then

[TABLE]

From the first relation in (4.19), (4.14) and the derivation of (4.14), we have

[TABLE]

and thus by the Bramble-Hilbert lemma as in (4.15),

[TABLE]

Substituting (4.31) into (4.28) results in

[TABLE]

Moreover, since $\bm{\psi}^{*}\in\bm{H}_{0}(\mathrm{div};\Omega)$, we find by considering (4.18) that

[TABLE]

and therefore by (4.32),

[TABLE]

The analysis for $J_{2,K}(u_{2},\psi_{h2}^{-})$ is similar, which along with (4.28), (4.29) (4.30) gives

[TABLE]

Collecting this result and (4.25), (4.26), (4.25) with (4.24), one gets

[TABLE]

where we have used the regularity (4.21). Finally, this theorem is established by dividing both sides by $\|\bm{e}_{h}\|_{0}$. ∎

Remark 4.5**.**

Again, the mixed finite element designed for Brinkman problems in [31] has similar properties as the pair $\bm{V}_{h}\times P_{h}$ in this work, and so the arguments in Theorems 4.1 and 4.4 are also valid. This explains the high accuracy phenomenon observed in the numerical tests in [31] under uniform rectangular partitions.

5 Numerical examples

Numerical tests are given in this section. The solution domain is set as $[0,2]\times[0,1]$, with uniform $n\times n$ rectangular partitions $\{\mathcal{T}_{h}\}$ constructed. It is clear that

[TABLE]

We first test $W_{h}$ for the biharmonic problem (3.1), where the exact solution is given by

[TABLE]

The errors in various norms are illustrated in Table 1. One sees that the convergence orders in discrete $H^{2}$- and $H^{1}$-norms are precisely $O(h^{2})$ and $O(h^{3})$, respectively, as predicted in Theorems 3.1 and 3.2. The $L^{2}$ error seems to be of order four, which will be studied in our future work.

As a comparison, we also check the performance of the well-known Adini element **[20]** for the same problem. The numbers of local DoFs are both 12 for both elements, and the global DoFs of the Adini element are fewer than those of $W_{h}$. However, the space $W_{h}$ is highly nonconforming, enjoying a cheap local communication when the method is implemented, especially for parallel computing. As the error analysis in **[20, 21, 17]**, the Adini element can also achieve a second order convergence rate in discrete $H^{2}$-norm. Moreover, the absolute errors are sightly lower than the those produced by $W_{h}$ due to the strong continuity in the tangential direction. However, the errors in discrete $H^{1}$- and $L^{2}$-norms are also only $O(h^{2})$, lower than those of $W_{h}$, which is consistent with the lower bound estimate given in **[16, 17]**.

We end this work by testing the divergence-free Stokes element $\bm{V}_{h}\times P_{h}$ for the model problem (4.1) determined by

[TABLE]

The numerical results are provided in Table 3. One can observe that the convergence orders for the velocity, the pressure and the postprocessed pressure are achieved, as predicted in Theorems 4.1 and 4.2. Moreover, the $L^{2}$ errors of the velocity have an $O(h^{3})$ convergence order, which agrees with the assertion in Theorem 4.4.