An optimal piecewise cubic nonconforming finite element scheme for the planar biharmonic equation on general triangulations

TL;DR

This paper introduces a new nonconforming finite element scheme using piecewise cubic polynomials for the planar biharmonic equation, achieving optimal convergence on general triangulations with a novel basis construction.

Contribution

It develops a piecewise cubic nonconforming finite element scheme with proven optimal convergence and a novel basis construction for the biharmonic equation on general triangulations.

Findings

Achieves $ ext{O}(h^2)$ convergence rate in energy norm.

Establishes an equivalence between discrete Hessian and Laplacian forms.

Constructs a locally supported basis without Ciarlet's triple.

Abstract

This paper presents a nonconforming finite element scheme for the planar biharmonic equation which applis piecewise cubic polynomials ($P_3$) and possesses $\mathcal{O}(h^2)$ convergence rate in energy norm on general shape-regular triangulations. Both Dirichlet and Navier type boundary value problems are studied. The basis for the scheme is a piecewise cubic polynomial space, which can approximate the $H^4$ functions with $\mathcal{O}(h^2)$ accuracy in broken $H^2$ norm. Besides, an equivalence $(\nabla_h^2\ \cdot,\nabla_h^2\ \cdot)=(\Delta_h\ \cdot,\Delta_h\ \cdot)$, which is usually not true for nonconforming finite element spaces, is proved on the newly designed spaces. The finite element space does not correspond to a finite element defined with Ciarlet's triple; however, a set of locally supported basis functions of the finite element space is still figured out. The notion of the finite element Stokes complex plays an important role in the analysis and also the construction of the basis functions.

Equations

$$\begin{array}[]{ccccccccc}0&\xrightarrow{\boldsymbol{inclusion}}&H^{2}_{0}&\xrightarrow{\boldsymbol{\mathrm{\nabla}}}&(H^{1}_{0})^{2}&\xrightarrow{\mathrm{rot}}&L^{2}_{0}&\xrightarrow{\boldsymbol{\int\cdot}}&0.\end{array}$$

$$\begin{array}[]{ccccccccc}0&\xrightarrow{\boldsymbol{inclusion}}&H^{2}\cap H^{1}_{0}&\xrightarrow{\boldsymbol{\mathrm{\nabla}}}&(H^{1})^{2}\cap H_{0}({\rm rot})&\xrightarrow{\mathrm{rot}}&L^{2}_{0}&\xrightarrow{\boldsymbol{\int\cdot}}&0.\end{array}$$

$${\mathcal{X}}_h^{b,+1}:=\{a\in {\mathcal{X}}_h^i,a\text{mbox}isconnectedto{\mathcal{X}}_h^{b}bye\in {\mathcal{E}}_h^i\},\text{mbox}and{\mathcal{X}}_h^{i,-1}:={\mathcal{X}}_h^i\setminus {\mathcal{X}}_h^{b,+1};$$

$${\mathcal{X}}_h^{b,+k}:=\{a\in {\mathcal{X}}_h^{i,-(k-1)},a\text{mbox}isconnectedto{\mathcal{X}}_h^{b,+(k-1)}bye\in {\mathcal{E}}_h^i\},\text{mbox}and{\mathcal{X}}_h^{i,-k}:={\mathcal{X}}_h^{i,-(k-1)}\setminus {\mathcal{X}}_h^{b,+k}.$$

$$\underset{\text{undertilde}v{}_h\in \text{undertilde}G{}_{h0}^2,\Vert {\nabla }_h\text{undertilde}v{}_h{\Vert }_{0,\Omega }=1}{\sup }({\mathrm{div}}_h\text{undertilde}v{}_h,q_h)⩾C\Vert q_h{\Vert }_{0,\Omega },\forall q_h\in {\mathbb{P}}_{h0}^1.$$

$$({\nabla }_h\text{undertilde}w{}_h,{\nabla }_h\text{undertilde}v{}_h)=(\mathrm{div}{}_h\text{undertilde}w{}_h,{\mathrm{div}}_h\text{undertilde}v{}_h)+({\mathrm{rot}}_h\text{undertilde}w{}_h,{\mathrm{rot}}_h\text{undertilde}v{}_h).$$

$$\begin{gathered} {\int }_K\nabla \text{undertilde}w{}_h:\nabla \text{undertilde}v{}_h=-{\int }_K\Delta \text{undertilde}w{}_h\text{undertilde}v{}_h+{\int }_{\partial K}{\partial }_{\mathbf{n}}\text{undertilde}w{}_h\text{undertilde}v{}_h=-{\int }_K\Delta \text{undertilde}w{}_h\text{undertilde}v{}_h=-{\int }_K(\nabla \mathrm{div}+\mathrm{curl}\mathrm{rot})\text{undertilde}w{}_h\text{undertilde}v{}_h \\ =-{\int }_K(\nabla \mathrm{div}+\mathrm{curl}\mathrm{rot})\text{undertilde}w{}_h\text{undertilde}v{}_h+{\int }_{\partial K}(\mathrm{div}\text{undertilde}w{}_h\text{undertilde}v{}_h\cdot \mathbf{n}+\mathrm{rot}\text{undertilde}w{}_h\text{undertilde}v{}_h\cdot \mathbf{t})={\int }_K\mathrm{div}\text{undertilde}w{}_h\mathrm{div}\text{undertilde}v{}_h+{\int }_K\mathrm{rot}\text{undertilde}w{}_h\mathrm{rot}\text{undertilde}v{}_h; \end{gathered}$$

$$\text{undertilde}w{}_h=\text{undertilde}w{}_h^1+\text{undertilde}w{}_h^2,\text{mbox}and\text{undertilde}v{}_h=\text{undertilde}v{}_h^1+\text{undertilde}v{}_h^2.$$

$$({\nabla }_h\text{undertilde}w{}_h,{\nabla }_h\text{undertilde}v{}_h)=({\nabla }_h\text{undertilde}w{}_h^1,{\nabla }_h\text{undertilde}v{}_h^1)+({\nabla }_h\text{undertilde}w{}_h^1,{\nabla }_h\text{undertilde}v{}_h^2)+({\nabla }_h\text{undertilde}w{}_h^2,{\nabla }_h\text{undertilde}v{}_h^1)+({\nabla }_h\text{undertilde}w{}_h^2,{\nabla }_h\text{undertilde}v{}_h^2),$$

$$\begin{array}[]{ccccccccc}0&\xrightarrow{\boldsymbol{inclusion}}&A^{3}_{h0}&\xrightarrow{\boldsymbol{\mathrm{\nabla}}_{h}}&\undertilde{G}{}_{h0}^{\rm 2,r}&\xrightarrow{\mathrm{rot}_{h}}&\mathbb{P}_{h0}^{1}&\xrightarrow{\int\cdot}&0.\end{array}$$

$$B_h^3:=\{v\in {\mathbb{P}}_h^3:v\text{mbox}iscontinuousata\in {\mathcal{X}}_h;{\int }_e[[v]]=0,\text{mbox}and{\int }_e{p}_e[[{\partial }_{\mathbf{n}}v]]=0,\forall p_e\in P_1(e),\forall e\in {\mathcal{E}}_h^i\},$$

$$B_{ht}^3:=\{v\in B_h^3:v(a)=0,a\in {\mathcal{X}}_h^b;{\int }_{e}v=0,\forall e\in {\mathcal{E}}_h^b\},$$

$$B_{h0}^3:=\{v\in B_{ht}^3:{\int }_e{p}_e{\partial }_{\mathbf{n}}v=0,\forall p_e\in P_1(e),\forall e\in {\mathcal{E}}_h^b\}.$$

$$({\nabla }_h^2{w}_h,{\nabla }_h^2{v}_h)=({\Delta }_h{w}_h,{\Delta }_h{v}_h).$$

$$({\nabla }_h^2{w}_h,{\nabla }_h^2{v}_h)=({\mathrm{div}}_h{\nabla }_h{w}_h,{\mathrm{div}}_h{\nabla }_h{v}_h)+({\mathrm{rot}}_h{\nabla }_h{w}_h,{\mathrm{rot}}_h{\nabla }_h{v}_h)=({\Delta }_h{w}_h,{\Delta }_h{v}_h),\forall w_h,v_h\in B_{ht}^3.$$

$$\begin{array}[]{ccccccccc}0&\xrightarrow{\boldsymbol{inclusion}}&B_{h0}^{3}&\xrightarrow{\boldsymbol{\mathrm{\nabla}}_{h}}&\undertilde{G}{}_{h0}^{2}&\xrightarrow{\mathrm{rot}_{h}}&\mathbb{P}_{h0}^{1}&\xrightarrow{\int\cdot}&0.\end{array}$$

$$\begin{array}[]{ccccccccc}0&\xrightarrow{\boldsymbol{inclusion}}&B_{ht}^{3}&\xrightarrow{\boldsymbol{\mathrm{\nabla}}_{h}}&\undertilde{G}{}_{ht}^{2}&\xrightarrow{\mathrm{rot}_{h}}&\mathbb{P}_{h0}^{1}&\xrightarrow{\int\cdot}&0.\end{array}$$

$$\mbox{Dirichlet\ type}\quad\left\{\begin{array}[]{rl}\Delta^{2}u=f&\mbox{in}\,\Omega,\\ u=\partial_{\mathbf{n}}u=0,&\mbox{on}\,\partial\Omega,\end{array}\right.$$

$$\mbox{Navier\ type}\quad\left\{\begin{array}[]{rl}\Delta^{2}z=f&\mbox{in}\,\Omega;\\ z=\Delta z=0,&\mbox{on}\,\partial\Omega.\end{array}\right.$$

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

$$({\nabla }^{2}z,{\nabla }^{2}v)=(f,v),\forall v\in H^2\cap H_0^1(\Omega ).$$

$$a_h(u_h,v_h):=({\nabla }_h^2{u}_h,{\nabla }_h^2{v}_h)=(f,v_h),\forall v_h\in B_{h0}^3;$$

$$a_h(z_h,v_h)=(f,v_h),\forall v_h\in B_{ht}^3.$$

$$\Vert {\nabla }_h^2(u-u_h){\Vert }_{0,\Omega }⩽C(h^{m-2}|u{|}_{m,\Omega }+h^2\Vert f{\Vert }_{0,\Omega }).$$

$$\Vert {\nabla }_h^2(z-z_h){\Vert }_{0,\Omega }⩽C(h^{m-2}|z{|}_{m,\Omega }+h^2\Vert f{\Vert }_{0,\Omega }).$$

$$\Vert {\nabla }_h(u-u_h){\Vert }_{0,\Omega }⩽C(h^{m-1}|u{|}_{m,\Omega }+h^3\Vert f{\Vert }_{0,\Omega }),$$

$$\Vert {\nabla }_h(u-u_h){\Vert }_{0,\Omega }⩽C(h^{m-2+\delta }|u{|}_{m,\Omega }+h^3\Vert f{\Vert }_{0,\Omega }),1/2<\delta ⩽1.$$

$$\left\{\begin{array}[]{lll}(\nabla\undertilde{\varphi},\nabla\undertilde{\psi})+({\rm rot}\undertilde{\psi},p)&=(-\Delta\undertilde{\varphi},\undertilde{\psi}),&\forall\,\undertilde{\psi}\in\undertilde{H}{}^{1}_{0}(\Omega);\\ ({\rm rot}\undertilde{\varphi},q)&=0,&\forall\,q\in L^{2}_{0}(\Omega).\end{array}\right.$$

$$\left\{\begin{array}[]{lll}(\nabla_{h}\undertilde{\varphi}{}_{h},\nabla_{h}\undertilde{\psi}{}_{h})+({\rm rot}_{h}\undertilde{\psi}{}_{h},p_{h})&=(-\Delta\undertilde{\varphi},\undertilde{\psi}{}_{h}),&\forall\,\undertilde{\psi}{}_{h}\in\undertilde{G}{}^{2}_{h0};\\ ({\rm rot}_{h}\undertilde{\varphi}{}_{h},q_{h})&=0,&\forall\,q_{h}\in\mathbb{P}_{h0}^{1}.\end{array}\right.$$

$${\mathbb{I}}_{h0}^{B}w:=w_h.$$

Full text

The finite element space does not correspond to a finite element defined with Ciarlet’s triple; however, a set of locally supported basis functions of the finite element space is still figured out. The notion of the finite element Stokes complex plays an important role in the analysis and also the construction of the basis functions.

The author is supported by NCMIS of CAS and NSFC Grants Nos 11471026 and 11871465.

1. Introduction

In order to obtain a simpler interior structure, in the study of the numerical analysis of partial differential equations, lower-degree polynomials are often expected to be used with respect to the same convergence rate. Finite element schemes with polynomials of degrees not higher than $k$ for $H^{\bf m}$ problem that possess convergence rates of $\mathcal{O}(h^{\bf k+1-m})$ in energy norm for solutions in $H^{\bf k+1}$ are called optimal. According to [27], this illustrates both the highest accuracy with respect to certain degree of polynomials and the smallest shape function space with respect to certain convergence rate, and is a critical characteristic for the finite element methodology. Motivated by the fundamental problem aforementioned, this paper concerns the optimal finite element scheme for the biharmonic equation with piecewise cubic polynomials on general triangulations.

A brief review of relevant works

Papers on optimal schemes can be found focusing mainly on low-order problems. For the lowest-differentiation-order ($H^{1}$) elliptic problems, the standard Lagrangian elements can yield optimal approximation on the simplicial grids of an arbitrary dimension. Further, the optimal nonconforming element spaces of $k$-th degrees are also constructed, c.f., e.g., [12], [19], and [11] for the cases $k=1$, $k=2$, and $k=3$, respectively, and [5] for general $k$. For higher-differentiation-order ($H^{m}$, $m>1$) elliptic problems, minimal-degree approximations have been studied with the lowest accuracy order. Specifically, when the subdivision comprises simplexes, a systematic family of nonconforming finite elements has been proposed by [39] for $H^{m}$ elliptic partial differential equations in $\mathbb{R}^{n}$ for any $n\geqslant m$ with polynomials with degree $m$. Besides, the constructions of finite element functions that do not depend on cell-by-cell definitions can be found in [32, 25, 47], wherein minimal-degree finite element spaces are defined on general quadrilateral grids for $H^{1}$ and $H^{2}$ problems. In contrast to these existing lowest order researches, the construction of higher-accuracy-order optimal schemes for higher-differentiation-order problems is complicated, even for the planar biharmonic problem.

Conforming finite elements for biharmonic equation requires the $C^{1}$ continuity assumption. It is well-known that with polynomials of degrees $k\geqslant 5$, spaces of $\mathcal{C}^{1}$ continuous piecewise polynomials can be constructed with local basis. Moreover, these spaces perform optimal approximations of $H^{2}$ functions with sufficient smoothness[2, 45, 28, 44, 15]. With polynomials of degrees $2\leqslant k\leqslant 4$, spaces of $\mathcal{C}^{1}$ continuous piecewise polynomials can be shown to provide optimal approximation when the triangulation is of some special structures, such as the Powell–Sabin and Powell–Sabin–Heindl triangulations[34, 23, 33], criss-cross triangulations [46], Hsieh–Clough–Tocher triangulation[9], and Sander–Veubeke triangulation[36, 17]. The conditions on the grids can be relaxed, but they are generally required on at least some part of the triangulation[31, 30, 8]. On general triangulations, as is shown in [16], optimal approximation cannot be obtained with $\mathcal{C}^{1}$ continuous piecewise polynomials of degree $k<5$. It is illustrated in [1] that not all the basis functions can be determined locally on general grids. We would particularly recall a counterexample that, as studied in [13, 14, 4], the $\mathcal{C}^{1}-P_{3}$ scheme is only $\mathcal{O}(h)$ order convergent in energy norm on a triangulation obtained by subdividing a rectangular domain with three groups of parallel lines (cf. Figure 1), which is even though one of the simplest and most regular triangulations.

In contrast, a nonconforming finite element methodology, namely, the Morley element [29], which uses piecewise quadratic polynomials with a convergence rate of $\mathcal{O}(h)$, was shown to perform optimally for $k=2$. However, to the best of our knowledge, optimal piecewise cubic or quartic finite element schemes (either conforming or nonconforming) for a planar biharmonic equation with $\mathcal{O}(h^{2})$ or $\mathcal{O}(h^{3})$ convergence rate have not been discovered. We remark that several $\mathcal{O}(h^{2})$ ordered finite element methods are designed with piecewise cubic polynomials enriched with higher-degree bubbles (e.g., [21, 38]). As the degrees of the functions exceed three, these methods are not considered optimal here. For a biharmonic problem in higher dimensions and other problems with higher orders, bigger difficulties can be expected.

Main results of the present paper

In this paper, a space $B^{3}_{h}$ is constructed with piecewise cubic polynomials, whose subspaces $B^{3}_{ht}$ and $B^{3}_{h0}$ are proved to provide optimal approximation of $H^{2}\cap H^{1}_{0}$ and $H^{2}_{0}$, respectively. Finite element schemes that apply the two subspaces to the biharmonic equation with Navier and Dirichlet boundary conditions, respectively, are nonconforming, but the consistency errors are both of $\mathcal{O}(h^{2})$ order. Thus the finite element schemes are optimal, and the optimality can be proved on any shape regular grids on both convex and nonconvex polygonal domains.

Further, for any two functions $w_{h},v_{h}\in B^{3}_{ht}$, it can be proved that $(\nabla_{h}^{2}w_{h},\nabla_{h}^{2}v_{h})=(\Delta_{h}w_{h},\Delta_{h}v_{h})$, which is seldom true for nonconforming finite elements. This property makes the finite element spaces fit for the discretization of biLaplacian operator $\Delta\mathcal{A}\Delta$ with varying coefficient $\mathcal{A}$.

Two approaches of implementing the schemes are suggested. One is to figure out their local basis functions: the finite element scheme does not correspond to a finite element in Ciarlet’s triple; but the finite element spaces do possess local basis functions that each is supported in the patch of a vertex or the patch of an edge. The other is to decompose the finite element scheme to three decoupled subproblems, which are either a Poisson system or a Stokes system, to solve sequentially. Note that the optimal solvers for discrete Poisson system and Stokes problem have been very well developed, and the latter approach suggests indeed a method to solve the finite element problem with optimal cost.

Main technical ingredient of the present paper

For the nonconforming finite element space $B_{h}^{3}$, to control the consistency error, sufficient restrictions on the interfacial continuity have to be imposed across the edges of the cells. However, the constraints on the continuity are overdetermined in comparison to local shape functions; hence, the global finite element space do not correspond to a local finite element defined with Ciarlet’s triple. The functions can be viewed as nonsmooth splines. Consequently, several challenges arise in both theoretical analysis and practical implementation, even on counting the dimension of the space. To avoid these challenges, in this paper, indirect methods are adopted; namely, the construction and utilization of discretized Stokes complexes constitute the bulk of the task in the construction of the space and schemes. This indirect approach is viewed as the main ingredient of the paper.

Discretized Stokes complexes are finite element analogs of the 2D Stokes complexes (or the de Rham complex with enhanced regularity), which read corresponding to the boundary condition:

[TABLE]

and

[TABLE]

In the complex, the combination of the successive two operators vanish, and the kernel of the latter one is exactly the range of the former one. The finite element complexes have been widely used for stability analysis (c.f.[3]), and, in this paper, the important role they play is four-folded:

1. (1)

It is used for approximation analysis. We construct two discretized Stokes complexes that start with finite element spaces $B^{3}_{h0}$ and $B^{3}_{ht}$, respectively, for $H^{2}$ and estimates the approximation error of $B^{3}_{h0}$($B^{3}_{ht}$) by estimating the discretization error of the auxiliary finite element discretization of the Stokes problem. This way, we prove the $\mathcal{O}(h^{2})$ approximation accuracy of $B^{3}_{h0}$($B^{3}_{ht}$) in energy norm for $H^{4}$ functions. Moreover, the proof does not require a convexity assumption on the domain. 2. (2)

Different from existing nonconforming finite elements such as the Morley element, for $w_{h},v_{h}\in B^{3}_{ht}$, $(\nabla_{h}^{2}w_{h},\nabla_{h}^{2}v_{h})=(\Delta_{h}w_{h},\Delta_{h}v_{h})$, the operations done cell by cell. This makes the finite element suitable for, e.g., $\Delta\mathcal{A}\Delta$ with varying coefficient $\mathcal{A}$; see [41] for a practical application. This property is, once again, proved by the aid of the discretized Stokes complex. 3. (3)

Further, though the finite element space does not correspond to a finite element defined in Ciarlet’s triple, the finite element spaces do admit a set of basis functions, each of which is supported in a patch of a vertex or a patch of an edge. Again, the discretized Stokes complexes play crucial roles in proving the existence of the locally supported basis functions. 4. (4)

Finally, we remark, beyond bringing ease in constructing and analyzing the schemes, the discretized Stokes complex is also helpful to the implementation and numerical solution of the systems by the aid of the discretized Poisson and discretized Stokes systems; we also refer to [42, 24, 43, 50, 35, 20, 18] for relevant discussions.

Bibliographic remark

This paper collects some original results from the unpublished arXiv preprints 1805.03851([48]) authored by the same author as the present paper.

Organization of the paper

The remaining of the paper is organized as follows. Section 2 presents some finite element spaces and finite element complexes. Section 3 presents two optimal nonconforming finite element schemes, including the construction, theoretical analysis, for the two kinds of boundary value problems, respectively. Two approaches of implementation are given in Section 4. Finally, in Section 5, some conclusions and further discussions are given.

2. Finite element spaces and finite element complexes

2.1. Preliminaries

In what follows, we use $\Omega$ to denote a simply connected polygonal domain, and $\nabla$, ${\rm curl}$, ${\rm div}$, ${\rm rot}$, and $\nabla^{2}$ to denote the gradient operator, curl operator, divergence operator, rot operator, and Hessian operator, respectively. As usual, we use $H^{2}(\Omega)$, $H^{2}_{0}(\Omega)$, $H^{1}(\Omega)$, $H^{1}_{0}(\Omega)$, $H({\rm rot},\Omega)$, $H_{0}({\rm rot},\Omega)$, and $L^{2}(\Omega)$ to denote certain Sobolev spaces, and specifically, denote $\displaystyle L^{2}_{0}(\Omega):=\{w\in L^{2}(\Omega):\int_{\Omega}wdx=0\}$, $\undertilde{H}{}^{1}_{0}(\Omega):=(H^{1}_{0}(\Omega))^{2}$, and $\undertilde{H}{}^{1}_{t}(\Omega)=(H^{1}(\Omega))^{2}\cap H_{0}({\rm rot},\Omega)$. Furthermore, we denoted vector-valued quantities by $``\undertilde{~{}}"$, while $\undertilde{v}{}^{1}$ and $\undertilde{v}{}^{2}$ denote the two components of the function $\undertilde{v}$. We use $(\cdot,\cdot)$ to represent $L^{2}$ inner product, and $\langle\cdot,\cdot\rangle$ to denote the duality between a space and its dual. Without ambiguity, we use the same notation $\langle\cdot,\cdot\rangle$ for different dualities, and it can occasionally be treated as $L^{2}$ inner product for certain functions. We use the subscript $``\cdot_{h}"$ to denote the dependence on triangulation. In particular, an operator with the subscript $``\cdot_{h}"$ indicates that the operation is performed cell-by-cell. Finally, $\;\stackrel{{\scriptstyle\raisebox{-11.09654pt}{=}}}{{{\sim}}}\;$ denotes equality up to a constant. The hidden constants depend on the domain, and when triangulation is involved, they also depend on the shape regularity of the triangulation, but they do not depend on $h$ or any other mesh parameter.

Let $\mathcal{T}_{h}$ be a shape-regular triangular subdivision of $\Omega$ with mesh size $h$, such that $\overline{\Omega}=\cup_{T\in\mathcal{T}_{h}}\overline{T}$. Denote by $\mathcal{E}_{h}$, $\mathcal{E}_{h}^{i}$, $\mathcal{E}_{h}^{b}$, $\mathcal{X}_{h}$, $\mathcal{X}_{h}^{i}$, $\mathcal{X}_{h}^{b}$ and $\mathcal{X}_{h}^{c}$ the set of edges, interior edges, boundary edges, vertices, interior vertices, boundary vertices and corners, respectively. For any edge $e\in\mathcal{E}_{h}$, denote by $\mathbf{n}_{e}$ and $\mathbf{t}_{e}$ the unit normal and tangential vectors of $e$, respectively, and denote by $\llbracket\cdot\rrbracket_{e}$ the jump of a given function across $e$; if particularly $e\in\mathcal{E}_{h}^{b}$, $\llbracket\cdot\rrbracket_{e}$ stands for the evaluation of the function on $e$. The subscript ${\cdot}_{e}$ can be dropped when there is no ambiguity brought in.

Denote

[TABLE]

further, denote with $\mathcal{X}^{i,-(k-1)}_{h}\neq\emptyset$,

[TABLE]

The smallest $k$ such that $\mathcal{X}_{h}^{i,-(k-1)}=\mathcal{X}_{h}^{b,+k}$ is called the number of levels of the triangulation.

For a triangle $T$, we use $P_{k}(T)$ to denote the set of polynomials on $K$ of degrees not higher than $k$. In a similar manner, $P_{k}(e)$ is defined on an edge $e$. We define $\undertilde{P}{}_{k}(T)=P_{k}(T)^{2}$ and similarly is $\undertilde{P}{}_{k}(e)$ defined. We use $a_{i}$, $i=1,2,3$ for the vertices of $T$ in an anticlockwise order, $e_{i}$, $i=1,2,3$ for the edges opposite to $a_{i}$, respectively, and $\lambda_{i}$, $i=1,2,3$ the barycentric coordinates.

Also, we denote basic finite element spaces by

• •

$\mathcal{L}^{k}_{h}:=\{w\in H^{1}(\Omega):w|_{T}\in P_{k}(T),\ \forall\,T\in\mathcal{T}_{h}\}$, $\mathcal{L}^{k}_{h0}:=\mathcal{L}^{k}\cap H^{1}_{0}(\Omega)$, $k\geqslant 1$;

• •

$\mathbb{P}^{k}_{h}:=\{w\in L^{2}(\Omega):w|_{T}\in P_{k}(T)\}$, $\mathbb{P}^{k}_{h0}:=\mathbb{P}^{k}_{h}\cap L^{2}_{0}(\Omega)$, $k\geqslant 0$;

• •

$\undertilde{S}{}_{h}^{k}:=(\mathbb{P}^{k}_{h})^{2}\cap\undertilde{H}{}^{1}(\Omega)$, $k\geqslant 1$, $\undertilde{S}{}^{k}_{ht}:=\undertilde{S}{}^{k}_{h}\cap H_{0}({\rm rot},\Omega)$ and $\undertilde{S}{}^{k}_{h0}:=\undertilde{S}{}^{k}_{h}\cap\undertilde{H}{}^{1}_{0}(\Omega)$;

• •

$\undertilde{G}{}_{h}^{k}:=\{\undertilde{v}\in(\mathbb{P}^{k}_{h})^{2}:\int_{e}p_{e}\llbracket\undertilde{v}^{j}\rrbracket=0,\ \forall\,p_{e}\in P_{k-1}(e),\ \forall\,e\in\mathcal{E}_{h}^{i},\ j=1,2\}$, $k\geqslant 1$, $\undertilde{G}{}^{k}_{ht}:=\{\undertilde{v}\in\undertilde{G}{}^{k}_{h}:\int_{e}p_{e}\undertilde{v}\cdot\mathbf{t}_{e}=0,\ \forall\,e\in\mathcal{E}_{h}^{b}\ \mbox{and}\ p_{e}\in P_{k-1}(e)\}$, and $\undertilde{G}{}^{k}_{h0}:=\{\undertilde{v}\in\undertilde{G}{}^{k}_{h}:\int_{e}p_{e}\undertilde{v}^{j}=0,\ \forall\,e\in\mathcal{E}_{h}^{b}\ \mbox{and}\ p_{e}\in P_{k-1}(e),\ j=1,2\}.$

Namely, $\undertilde{S}{}^{k}_{h}$ consists of continuous functions, and $\undertilde{G}{}^{k}_{h}$ consists of $(k-1)^{\rm th}$ order moment-continuous functions. Particularly, the space $\undertilde{G}{}^{2}_{h}$ corresponds to the famous Fortin-Soulie element [19]. The following stability result is well-known.

Lemma 1**.**

[19]** There exists a generic constant $C$ depending on the domain and the regularity of the grid, such that

[TABLE]

Remark 2**.**

By the symmetry between the two components of $\undertilde{H}{}^{1}(\Omega)$, Lemma 1 remains true when “${\rm div}{}_{h}$” is replaced by “${\rm rot}_{h}$.”

Denote $\undertilde{\mathcal{B}}{}^{2}_{h0}:=\{\undertilde{\phi}{}_{h}:(\undertilde{\phi}{}_{h}|_{T})^{j}\in{\rm span}\{(\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2})-2/3\},\ j=1,2,\ \forall\,T\in\mathcal{T}_{h}\}$ and evidently the first order moments of $\undertilde{\phi}{}_{h}$ vanish along any edge of $\mathcal{T}_{h}$. Then $\undertilde{G}{}^{2}_{h0}=\undertilde{S}{}^{2}_{h0}\oplus\undertilde{\mathcal{B}}{}^{2}_{h0}$ (c.f. [19]). Also $\undertilde{G}{}^{2}_{ht}=\undertilde{S}{}^{2}_{ht}\oplus\undertilde{\mathcal{B}}{}^{2}_{h0}$. Note that, as is known, $\undertilde{G}{}^{2}_{h}=\undertilde{S}{}^{2}_{h}+\undertilde{\mathcal{B}}{}^{2}_{h0}$ is not a direct sum. The decomposition can be generalized to even $k$ (c.f. [5]).

Lemma 3**.**

For any $\undertilde{w}{}_{h},\undertilde{v}{}_{h}\in\undertilde{G}{}^{2}_{ht}$, it holds that

[TABLE]

Proof.

Firstly, (4) holds for any $\undertilde{w}{}_{h},\undertilde{v}{}_{h}\in\undertilde{S}{}^{2}_{ht}\subset\undertilde{H}{}^{1}_{t}(\Omega)$. Secondly, (4) holds for any $\undertilde{w}{}_{h}\in\undertilde{G}{}^{2}_{ht}$ and $\undertilde{v}{}_{h}\in\mathcal{B}_{h0}$; actually, for any $K\in\mathcal{T}_{h}$,

[TABLE]

here we have used the fact that $\partial_{\mathbf{n}}\undertilde{w}{}_{h}$, ${\rm div}\undertilde{w}{}_{h}$ and ${\rm rot}\undertilde{w}{}_{h}$ are all linear polynomials along the edges of $K$ and that the first order moments of $\undertilde{v}{}_{h}$ vanish along the edges of $K$.

Now, given $\undertilde{w}{}_{h},\undertilde{v}{}_{h}\in\undertilde{G}{}^{2}_{ht}$, there exist uniquely $\undertilde{w}{}_{h}^{1},\undertilde{v}{}_{h}^{1}\in\undertilde{S}{}^{2}_{ht}$ and $\undertilde{w}{}_{h}^{2},\undertilde{v}{}_{h}^{2}\in\undertilde{\mathcal{B}}{}^{2}_{h0}$, such that

[TABLE]

Thus

[TABLE]

and $({\rm div}_{h}\undertilde{w}{}_{h},{\rm div}_{h}\undertilde{v}{}_{h})+({\rm rot}_{h}\undertilde{w}{}_{h},{\rm rot}_{h}\undertilde{v}{}_{h})$ can be decomposed to four corresponding parts. Then (4) can be established for every pair of the parts, and the proof is completed. ∎

Remark 4**.**

It is known that (4) holds for $\undertilde{H}{}^{1}_{t}$ functions but in general not for nonconforming finite element functions (such as the Crouzeix-Raviart element functions). This lemma reveals that the nonconforming space $\undertilde{G}{}^{2}_{ht}$ is in some sense like a conforming one.

2.2. An auxiliary finite element Stokes complex

Given a grid $\mathcal{T}_{h}$, define

• •

$A^{3}_{h}:=\{w_{h}\in L^{2}(\Omega):w_{h}|_{T}\in P_{3}(T);w_{h}(a)\ \mbox{is\ continuous\ at}\ a\in\mathcal{X}_{h}\}$;

• •

$A^{3}_{h0}:=\{w_{h}\in A^{3}_{h}:w_{h}(a)=0\ at\ a\in\mathcal{X}_{h}^{b}\};$

• •

$\undertilde{G}{}_{h}^{\rm 2,r}:=\{\undertilde{v}\in(\mathbb{P}^{2}_{h})^{2};\ \int_{e}\llbracket\undertilde{v}\cdot\mathbf{t}_{e}\rrbracket=0,\ \forall\,e\in\mathcal{E}_{h}^{i}\}$;

• •

$\undertilde{G}{}_{h0}^{\rm 2,r}:=\{\undertilde{v}\in\undertilde{G}{}_{h}^{\rm 2,r},\ \int_{e}\undertilde{v}\cdot\mathbf{t}_{e}=0,\ \forall\,e\in\mathcal{E}_{h}^{b}\}.$

Lemma 5**.**

A finite element complex is given by

[TABLE]

Proof.

We adopt the standard counting technique.

Firstly, by Lemma 1, $\mathbb{P}^{1}_{h0}={\rm rot}_{h}\undertilde{G}{}^{2}_{h0}\subset{\rm rot}_{h}\undertilde{G}{}^{\rm 2,r}_{h0}\subset\mathbb{P}^{1}_{h0}$. Secondly, $\nabla_{h}A^{3}_{h0}\subset\{\undertilde{v}{}_{h}\in\undertilde{G}{}^{\rm 2,r}_{h0}:{\rm rot}_{h}\undertilde{v}{}_{h}=0\}$. Thus we only have to check if $\dim(\nabla_{h}A^{3}_{h0})+\dim(\mathbb{P}^{1}_{h0})=\dim(\undertilde{G}{}^{\rm 2,r}_{h0})$, which can be verified by observing that $\dim(A^{3}_{h0})=\#(\mathcal{X}_{h}^{i})+7\#(\mathcal{T}_{h})$, $\dim(\undertilde{G}{}^{\rm 2,r}_{h0})=\#(\mathcal{E}_{h}^{i})+9\#(\mathcal{T}_{h})$ and $\dim(\undertilde{P}{}^{1}_{h0})=3\#(\mathcal{T}_{h})-1$, and by the Euler formula. The proof is completed. ∎

2.3. Finite element spaces for $H^{2}$ and discretized Stokes complexes

Define

[TABLE]

[TABLE]

and

[TABLE]

According to the boundary conditions on $B^{3}_{h0}$, we can recognize them as for $H^{2}$ problems.

Remark 6**.**

Note that, given $v_{h}\in B^{3}_{h}$, on every cell, $v_{h}$ is embedded in 12 restrictions. We can not expect $B^{3}_{h}$ correspond to a finite element defined with Ciarlet’s triple.

Lemma 7**.**

$B^{3}_{ht}=\{w_{h}\in A^{3}_{h0}:\nabla_{h}w_{h}\in\undertilde{G}{}^{2}_{ht}\}$, and $B^{3}_{h0}=\{w_{h}\in A^{3}_{h0}:\nabla_{h}w_{h}\in\undertilde{G}{}^{2}_{h0}\}$.

Proof.

Firstly, by an elementary calculus, the continuity restriction of $B^{3}_{h}$ implies that $\int_{e}p_{e}\llbracket\partial_{\mathbf{t}}v_{h}\rrbracket=0$ for any $p_{e}\in P_{1}(e)$, any $e\in\mathcal{E}_{h}^{i}$ and any $v_{h}\in B^{3}_{h}$. Also, $\int_{e}p_{e}\partial_{\mathbf{t}}v_{h}=0$ for any $p_{3}\in P_{1}(e)$, any $e\in\mathcal{E}_{h}^{b}$ and any $v_{h}\in B^{3}_{ht}$.

By the definitions of $B^{3}_{h0}$ and $A^{3}_{h0}$, $B^{3}_{h0}\subset\{w_{h}\in A^{3}_{h0}:\nabla w_{h}\in\undertilde{G}{}^{2}_{h0}\}$. On the other hand, given $w_{h}\in A^{3}_{h0}$ such that $\nabla_{h}w_{h}\in\undertilde{G}{}^{2}_{h0}$, then $\int_{e}\llbracket\partial_{\bf n_{e}}w_{h}\rrbracket p_{e}=\int_{e}\llbracket\partial_{\bf t_{e}}w_{h}\rrbracket p_{e}=0$ for any $e\in\mathcal{E}_{h}$ and $p_{e}\in P_{1}(e)$. This implies $w_{h}\in B^{3}_{h0}$. Namely $B^{3}_{h0}=\{w_{h}\in A^{3}_{h0}:\nabla_{h}w_{h}\in\undertilde{G}{}^{2}_{h0}\}$. Similarly can $B^{3}_{ht}=\{w_{h}\in A^{3}_{h0}:\nabla_{h}w_{h}\in\undertilde{G}{}^{2}_{ht}\}$ be proved, and all the proof is completed. ∎

Lemma 8**.**

It holds for $w_{h},v_{h}\in B^{3}_{ht}$ that

[TABLE]

Proof.

By Lemma 3, as $\nabla_{h}B^{3}_{ht}\subset\undertilde{G}{}^{2}_{ht}$,

[TABLE]

The proof is completed. ∎

Remark 9**.**

The lemma reveals that the functions in $B^{3}_{ht}$ possess some property like the $H^{2}$ conforming functions.

Theorem 10**.**

Two discretized Stokes complexex are given by

[TABLE]

and

[TABLE]

Proof.

By Lemmas 1, given $p_{h}\in\mathbb{P}^{1}_{h0}$, there exists $\undertilde{\sigma}{}_{h}\in\undertilde{G}{}^{2}_{h0}$, such that ${\rm rot}_{h}\undertilde{\sigma}{}_{h}=p_{h}$. Further, given $\undertilde{\tau}{}_{h}\in\undertilde{G}{}^{2}_{h0}$, such that ${\rm rot}_{h}\undertilde{\tau}{}_{h}=0$, by Lemma 5, there exists $w_{h}\in A^{3}_{h0}$, such that $\undertilde{\tau}{}_{h}=\nabla_{h}w_{h}$. Further, by Lemma 7, $w_{h}\in B^{3}_{h0}$. Therefore, (7) is proved. Similarly can (8) be proved. ∎

Remark 11**.**

A key feature for the proof of Theorem 10 is to construct a bigger finite element complex to cover, e.g., (7); this is accomplished by Lemma 5, where a finite element complex is constructed where the same piecewise polynomial space with lower regularity is used corresponding to (7). A dual way can be to use bigger piecewise polynomial space with the same regularity. A different proof of (7) can be found along this line in [48].

3. Optimal nonconforming finite element schemes for biharmonic equation

We consider the biharmonic equation with $f\in L^{2}(\Omega)$:

[TABLE]

and

[TABLE]

The variational problems are respectively

• •

to find $u\in H^{2}_{0}(\Omega)$ such that

[TABLE]

• •

to find $z\in H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$, such that

[TABLE]

In this section, we consider the nonconforming finite element discretization for them:

• •

find $u_{h}\in B^{3}_{h0}$ such that

[TABLE]

• •

find $z_{h}\in B^{3}_{ht}$ such that

[TABLE]

By the weak continuity of $B^{3}_{ht}$, $|\cdot|_{2,h}$ (namely, $\|\nabla_{h}^{2}\cdot\|_{0,\Omega}$) is a norm on $B^{3}_{ht}$, and (13) and (14) are well-posed.

The main result of this section is contained in the theorem below.

Theorem 12**.**

Let $u$, $u_{h}$, $z$ and $z_{h}$ be solutions of (11) and (13), (12), and (14), respectively. Then, with a generic constant $C$ depending on $\Omega$ and the regularity of the grid only, it holds for $u,z\in H^{m}(\Omega)$, $m=3,4$, that

[TABLE]

and

[TABLE]

Moreover, when $\Omega$ is convex,

[TABLE]

and

[TABLE]

When $\Omega$ is specifically a rectangle, $\delta=1$.

We postpone the proof of Theorem 12 after some technical lemmas.

3.1. Approximation property of $B^{3}_{h0}$

First of all, we define an interpolator to $B^{3}_{h0}$. Given $w\in H^{3}(\Omega)\cap H^{2}_{0}(\Omega)$, set $\undertilde{\varphi}:=\nabla w$, then $\undertilde{\varphi}\in\undertilde{H}^{2}(\Omega)\cap\undertilde{H}{}^{1}_{0}(\Omega)$ and ${\rm rot}\,\undertilde{\varphi}=0$. Indeed, $(\undertilde{\varphi},p\equiv 0)$ solves the incompressible Stokes equation:

[TABLE]

Now, choose $(\undertilde{\varphi}{}_{h},p_{h})\in\undertilde{G}{}^{2}_{h0}\times\mathbb{P}_{h0}^{1}$ such that

[TABLE]

Then, by Theorem 10, there exists a unique $w_{h}\in B^{3}_{h0}$ such that $\nabla_{h}w_{h}=\undertilde{\varphi}{}_{h}$. This way, we define an interpolation operator $\mathbb{I}_{h0}^{B}:H^{3}(\Omega)\cap H^{2}_{0}(\Omega)\to B^{3}_{h0}$ by

[TABLE]

Lemma 13**.**

There exists a constant $C$ such that for any $w\in H^{2}_{0}(\Omega)\cap H^{m}(\Omega)$, $m=3,4$, it holds for $k=2$ that

[TABLE]

If $\Omega$ is convex, then (22) holds for $k=1,2$.

Proof.

By definition, the interpolation error of $\mathbb{I}^{B}_{h0}$ is the discretization error of (20), and (22) can be obtained by standard technique (with $\Omega$ either convex or nonconvex). ∎

3.2. Approximation of $B^{3}_{ht}$

Again, we firstly define an interpolator to $B^{3}_{ht}$. Given $w\in H^{3}(\Omega)\cap H^{1}_{0}(\Omega)$ such that $\Delta w|_{\Gamma}=0$, set $\undertilde{\varphi}:=\nabla w$, then $\undertilde{\varphi}\in\undertilde{H}^{2}(\Omega)\cap H_{0}({\rm rot},\Omega)$, ${\rm rot}\,\undertilde{\varphi}=0$ and $({\rm div}\undertilde{\varphi})|_{\Gamma}=0$. Indeed, $(\undertilde{\varphi},p\equiv 0)$ solves the incompressible Stokes equation:

[TABLE]

Now, choose $(\undertilde{\varphi}{}_{h},p_{h})\in\undertilde{G}{}^{2}_{ht}\times\mathbb{P}_{h0}^{1}$ such that

[TABLE]

Then, by Theorem 10, there exists a unique $w_{h}\in B^{3}_{ht}$ such that $\nabla_{h}w_{h}=\undertilde{\varphi}{}_{h}$. This way, we define an interpolation operator $\mathbb{I}_{ht}^{B}:H^{3}(\Omega)\cap H^{1}_{0}(\Omega)\to B^{3}_{ht}$ by

[TABLE]

Lemma 14**.**

There exists a constant $C$ such that for any $w\in H^{1}_{0}(\Omega)\cap H^{m}(\Omega)$ such that $\Delta w|_{\Gamma}=0$, $m=3,4$, it holds that

[TABLE]

If $\Omega$ is convex, then

[TABLE]

If specifically $\Omega$ is rectangle, $\kappa=2$.

Proof.

By definition, the interpolation error of $\mathbb{I}^{B}_{ht}$ is the discretization error of (24), and (26) and (27) can be obtained by standard technique (with $\Omega$ either convex or nonconvex). We only have to note that the regularity of the auxiliary Stokes problem on convexs domain can be affected under the boundary condition of this kind. Specifically, we refer to [6] for the full regularity of (10) and thus of the auxiliary Stokes problem (23) on rectangles. ∎

3.3. Convergence analysis of the nonconforming scheme

For suitable $\varphi$ and $\psi$, define the bilinear forms

[TABLE]

[TABLE]

and

[TABLE]

Lemma 15**.**

There exists a constant $C$ such that it holds for any $\varphi\in H^{2}_{0}(\Omega)\cap H^{k}(\Omega)$, $w_{h}\in B^{3}_{h0}+H^{2}_{0}(\Omega)$, and $k=3,4$ that,

[TABLE]

[TABLE]

Proof.

Given $e\in\mathcal{E}_{h}$, by the definition of $B_{h0}^{3}$, $\fint_{e}p_{e}\llbracket\partial_{\mathbf{n}_{e}}w_{h}\rrbracket_{e}=0,\ p_{e}\in P_{1}(e)$; for the tangential direction, $\fint_{e}p_{e}\llbracket\partial_{\mathbf{t}_{e}}w_{h}\rrbracket_{e}=(p_{e}(L_{e})\llbracket w_{h}\rrbracket_{e}(L_{e})-p_{e}(R_{e})\llbracket w_{h}\rrbracket_{e}(R_{e}))-\fint_{e}\partial_{\mathbf{t}_{e}}p_{e}\llbracket w_{h}\rrbracket_{e}=0$. Hence,

[TABLE]

Therefore, (31) follows by standard techniques.

Now, define $\Pi_{h}^{2}$ the nodal interpolation to $\mathcal{L}_{h0}^{2}$ by

[TABLE]

It is easy to verify that the operator is well-defined. Moreover,

[TABLE]

By Green’s formula,

[TABLE]

Therefore,

[TABLE]

By (34),

[TABLE]

Further,

[TABLE]

Summing all above proves (37). ∎

Similarly, we have the lemma below.

Lemma 16**.**

There exists a constant $C$ such that it holds for any $\varphi\in H^{1}_{0}(\Omega)\cap H^{k}(\Omega)$ so that $(\Delta\varphi)|_{\Gamma}=0$, $w_{h}\in B^{3}_{ht}+H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$ and $k=3,4$ that,

[TABLE]

[TABLE]

Proof of Theorem 12

The proof follows a similar approach as the one in [37], with some technical modifications. By Strang lemma,

[TABLE]

The approximation error estimate follows by Lemma 13. By Lemma 15,

[TABLE]

which completes the proof of (15).

Now, we turn our attention to the proof of (17) for convex $\Omega$. Denote $u^{\Pi}_{h}=\mathbb{I}^{B}_{h0}u$. Then, by Lemma 13, $\|\nabla_{h}^{j}(u-u^{\Pi}_{h})\|_{0,\Omega}\leqslant Ch^{4-j}|u|_{4,\Omega}$, $j=1,2$. Denote by $\Pi_{h}^{1}$ the nodal interpolation onto $\mathcal{L}^{1}_{h0}$, then $\Pi_{h}^{1}(u^{\Pi}_{h}-u_{h})\in H^{1}_{0}(\Omega)$. Set $\varphi\in H^{3}(\Omega)\cap H^{2}_{0}(\Omega)$ such that

[TABLE]

then when $\Omega$ is convex, \|\varphi\|_{3,\Omega}\raisebox{-4.2679pt}{;\stackrel{{\scriptstyle\raisebox{-11.09654pt}{=}}}{{{\sim}}};}\|\Pi_{h}^{1}(u^{\Pi}_{h}-u_{h})\|_{1,\Omega}. By Green’s formula,

[TABLE]

Further, set $\varphi_{h}^{\Pi}=\mathbb{I}^{B}_{h0}\varphi$, and

[TABLE]

Therefore, $\|\nabla\Pi_{h}^{1}(u^{\Pi}_{h}-u_{h})\|_{0,\Omega}^{2}\leqslant C|\varphi|_{3,\Omega}(h^{m-1}|u|_{m,\Omega}+h^{3}\|\Delta^{2}u\|_{0,\Omega}),$ and $\|\nabla\Pi_{h}^{1}(u_{h}^{\Pi}-u_{h})\|_{0,\Omega}\leqslant C(h^{m-1}|u|_{m,\Omega}+h^{3}\|\Delta^{2}u\|_{0,\Omega})$. Finally,

[TABLE]

The proof of (18) and (16) is basically the same. The convergence rate for the $H^{1}$ norm of the error is slightly lost due to the lost of the regularity of the model problem (9) on general convex polygons. The proof is completed.∎

3.4. A variant formulation for bi-Laplacian equation with varying coefficient

The bi-Laplacian equation $\Delta(\mathcal{A}\Delta u)=f$, where $\mathcal{A}$ is a non-constant coefficient with positive lower and upper bounds, is frequently dealt with in applications. The equation arises in, e.g., the Helmholtz transmission eigenvalue problem in acoustics (c.f., e.g., [10, 26, 40]). The variational problem is to find $u\in H^{2}_{0}(\Omega)$ such that

[TABLE]

Correspondingly, we consider the nonconforming finite element discretization:

find $u_{h}\in B^{3}_{h0}$ such that

[TABLE]

Lemma 17**.**

The finite element problem (39) admits a unique solution.

Proof.

By Lemma 8, the bilinear form $\tilde{a}_{h}(\cdot,\cdot)$ is coercive on $B^{3}_{h0}$ with respect to the norm $|\cdot|_{2,h}$. The well-posedness of (39) follows by Lax-Milgrem lemma. The proof is completed. ∎

Similar to Theorem 12, we can establish and prove the theorem below.

Theorem 18**.**

Let $u$ and $u_{h}$ be solutions of (38) and (39), respectively. Then, with a generic constant $C$ depending on $\mathcal{A}$, $\Omega$ and the regularity of the grid only, it holds for $u\in H^{m}(\Omega)$, $m=3,4$, that

[TABLE]

Moreover, when $\Omega$ is convex,

[TABLE]

Remark 19**.**

For the bi-Laplacian equation with non-constant coefficient $\mathcal{A}$, the finite element scheme of the formulation (39) is a natural alternative. When the formulation (39) is used on, e.g., the Morley element, however, the scheme is not well-posed without extra stabilisations. Higher regularity of $B^{3}_{h0}$ here makes it fit for the formulation (39).

Remark 20**.**

Similarly, by Lemma 8, a bilinear form induced by $(\mathcal{A}\Delta_{h}\cdot,\Delta_{h}\cdot)$ can be used for $\Delta(\mathcal{A}\Delta u)=f$ with Navier type boundary condition.

4. On the implementation of the schemes

In this section, we present two approaches to implement the schemes. One is to figure out the locally supported basis functions of $B^{3}_{h0}$ and $B^{3}_{ht}$, and the other is to decompose the finite element system to three sub-problems to be solved sequentially. The former approach makes the scheme fit for the general finite element programing procedure, and the latter approach, as the sub-problems are Poisson systems and a Stokes system, makes the finite element problems optimally solvable.

4.1. Locally supported basis functions of the finite element spaces

4.1.1. Structure of weakly rot-free space

Let $T$ be a triangle with $a_{i}$ and $e_{i}$, $i=1,2,3$, being its vertices and edges. Define $\undertilde{P}{}_{2}({\rm rot},w0):=\{\undertilde{p}\in\undertilde{P}{}_{2}(T):\int_{T}{\rm rot}\undertilde{p}=0\}$. Then $\dim(\undertilde{P}{}_{2}({\rm rot},w0))=11$.

Denote

$\undertilde{\eta}{}_{a_{i}}^{x}$**: **

such that $\undertilde{\eta}{}_{a_{i}}^{x}(a_{j})=(\delta_{ij},0)^{\top}$; $\fint_{e_{k}}\undertilde{\eta}{}_{a_{i}}^{x}=\undertilde{0}$; $i,j,k=1,2,3$;

$\undertilde{\eta}{}_{a_{i}}^{y}$**: **

such that $\undertilde{\eta}{}_{a_{i}}^{y}(a_{j})=(0,\delta_{ij})^{\top}$; $\fint_{e_{k}}\undertilde{\eta}{}^{y}_{a_{i}}=\undertilde{0}$; $i,j,k=1,2,3$.

$\undertilde{\eta}{}_{a_{i}}$**: **

such that $\undertilde{\eta}{}_{a_{i}}(a_{j})=\undertilde{0}$; $\fint_{e_{k}}\undertilde{\eta}{}_{a_{i}}\cdot\mathbf{t}_{e_{k},a_{i}}=1-\delta_{ik}$ and $\int_{e_{k}}\undertilde{\eta}{}_{a_{i}}\cdot\mathbf{n}_{e_{k},a_{i}}=0$, where $\mathbf{t}_{e_{k},a_{i}}$ is the unit tangential vector along $e_{k}$ starting from $a_{i}$ and $\mathbf{n}_{e_{k},a_{i}}$ is the normal direction of $e_{k}$; $i,j,k=1,2,3$;

$\undertilde{\eta}{}_{e_{i}}$**: **

such that $\undertilde{\eta}{}_{e_{i}}(a_{j})=\undertilde{0}$; $\fint_{e_{k}}\undertilde{\eta}{}_{e_{i}}\cdot\undertilde{\tau}{}_{e_{k}}=0$, $\fint_{e_{k}}\undertilde{\eta}{}_{e_{i}}\cdot\mathbf{n}_{e_{k}}=\delta_{ik}$; $i,j,k=1,2,3$.

The functions form a frame of $\undertilde{P}{}_{2}({\rm rot},w0)$. Indeed, $\undertilde{\eta}{}_{e_{1}}+\undertilde{\eta}{}_{e_{2}}+\undertilde{\eta}{}_{e_{3}}=0$, while we have the lemma below.

Lemma 21**.**

All $\undertilde{\eta}{}^{x}_{a_{i}}$, $\undertilde{\eta}{}^{y}_{a_{i}}$ and $\undertilde{\eta}{}_{e_{i}}$ for $i=1,2,3$ and any two of $\undertilde{\eta}{}_{a_{i}}$ among $i=1,2,3$ form a basis of $\undertilde{P}{}_{2}({\rm rot},w0)$.

Analogically, denote $\undertilde{S}{}^{2}_{h}({\rm rot},w0):=\{\undertilde{v}\in\undertilde{S}{}^{2}_{h}:({\rm rot}\undertilde{v},q)=0,\ \forall\,q\in\mathbb{P}^{0}_{h}\}$, $\undertilde{S}{}^{2}_{h0}({\rm rot},w0):=\{\undertilde{v}\in\undertilde{S}{}^{2}_{h0}:({\rm rot}\undertilde{v},q)=0,\ \forall\,q\in\mathbb{P}^{0}_{h0}\}$ and $\undertilde{S}{}^{2}_{ht}({\rm rot},w0):=\{\undertilde{v}\in\undertilde{S}{}^{2}_{ht}:({\rm rot}\undertilde{v},q)=0,\ \forall\,q\in\mathbb{P}^{0}_{h0}\}$.

Meanwhile, for $a\in\mathcal{X}_{h}$, denote by $P_{a}$ the union of triangles of which $a$ is a vertex, namely the patch associated with $a$; for $e\in\mathcal{E}_{h}$, denote by $P_{e}$ the patch associated with $e$. Denote, with respect to $a\in\mathcal{X}_{h}$ and $e\in\mathcal{E}_{h}$, functions in $\undertilde{S}{}^{2}_{h}$ as,

$\undertilde{\varphi}{}_{a}^{x}$**: **

such that $\undertilde{\varphi}{}_{a}^{x}(a)=(1,0)^{\top}$; $\undertilde{\varphi}{}_{a}^{x}(a^{\prime})=\undertilde{0}$ on $a\neq a^{\prime}\in\mathcal{X}_{h}$; $\fint_{e^{\prime}}\undertilde{\varphi}{}_{a}^{x}=\undertilde{0}$ on $e^{\prime}\in\mathcal{E}_{h}$;

$\undertilde{\varphi}{}_{a}^{y}$**: **

such that $\undertilde{\varphi}{}_{a}^{y}(a)=(0,1)^{\top}$; $\undertilde{\varphi}{}_{a}^{y}(a^{\prime})=\undertilde{0}$ on $a\neq a^{\prime}\in\mathcal{X}_{h}$; $\fint_{e^{\prime}}\undertilde{\varphi}{}^{y}_{a}=\undertilde{0}$ on $e^{\prime}\in\mathcal{E}_{h}^{i}$;

$\undertilde{\varphi}{}_{P_{a}}$**: **

such that $\undertilde{\varphi}{}_{P_{a}}(a^{\prime})=\undertilde{0}$ on $a^{\prime}\in\mathcal{X}_{h}$; $\fint_{e}\undertilde{\varphi}{}_{P_{a}}=\undertilde{0}$ on $e\in\mathcal{E}_{h}$ and $a\not\in e$; $\fint_{e}\undertilde{\varphi}{}_{P_{a}}\cdot\mathbf{t}_{e,P_{a}}=1$ and $\int_{e}\undertilde{\varphi}{}_{P_{a}}\cdot\mathbf{n}_{e,P_{a}}=0$ on $e\subset P_{a}$ and $a\in e$, where $\mathbf{t}_{e,P_{a}}$ is the unit tangential vector along $e$ starting from $a$ and $\mathbf{n}_{e,P_{a}}$ is the anticlockwise normal direction of $e$ with respect to $P_{a}$;

$\undertilde{\varphi}{}_{e}$**: **

such that $\fint_{e}\undertilde{\varphi}{}_{e}\cdot\undertilde{\tau}{}_{e}=0$, $\fint_{e}\undertilde{\varphi}{}_{e}\cdot\mathbf{n}_{e}=1$, and $\undertilde{\varphi}{}_{e}$ vanishes on $\Omega\setminus\mathring{P_{e}}$.

Lemma 22**.**

The set $\{\undertilde{\varphi}{}^{x}_{a},\undertilde{\varphi}{}^{y}_{a},\undertilde{\varphi}{}_{P_{a}},\undertilde{\varphi}{}_{e}\}_{a\in\mathcal{X}^{i}_{h},\ e\in\mathcal{E}^{i}_{h}}$ forms a basis of $\undertilde{S}{}^{2}_{h0}({\rm rot},w0)$; namely

[TABLE]

Proof.

By direct calculation, the functions $\undertilde{\varphi}{}_{a}^{x}$, $\undertilde{\varphi}{}_{a}^{y}$, $\undertilde{\varphi}{}_{P_{a}}$ and $\undertilde{\varphi}{}_{e}$ all belong to $\undertilde{S}{}^{2}_{h0}({\rm rot},w0)$. By their definitions, the functions $\{\undertilde{\varphi}{}^{x}_{a},\undertilde{\varphi}{}^{y}_{a},\undertilde{\varphi}{}_{e}\}_{a\in\mathcal{X}_{h}^{i},\ e\in\mathcal{E}_{h}^{i}}$ are linearly independent, and the summation ${\rm span}\{\undertilde{\varphi}{}_{P_{a}}\}_{a\in\mathcal{X}_{h}^{i}}+{\rm span}\{\undertilde{\varphi}{}^{x}_{a},\undertilde{\varphi}{}^{y}_{a},\undertilde{\varphi}{}_{e}\}_{a\in\mathcal{X}_{h}^{i},\ e\in\mathcal{E}_{h}^{i}}$ is direct. Since $\dim(\undertilde{S}{}^{2}_{h0}({\rm rot},w0))=\dim(\undertilde{S}{}^{2}_{h0})-\dim(\mathbb{P}^{0}_{h0})=\#(\mathcal{(}N)_{h0}^{i})+3\#(\mathcal{X}_{h0}^{i})$, it remains for us to show $\{\undertilde{\varphi}{}_{e}\}_{e\in\mathcal{E}_{h}^{i}}$ are linearly independent.

Assume there exist $\{\alpha_{a}\}_{a\in\mathcal{X}}\subset\mathbb{R}$ with $\alpha_{a}=0$ for $a\in\mathcal{X}_{h}^{b}$, such that $\undertilde{\psi}=\sum_{a\in\mathcal{X}_{h}}\alpha_{a}\undertilde{\varphi}{}_{P_{a}}\equiv 0$. By the definition of $\undertilde{\varphi}{}_{P_{a}}$, for any $e\in\mathcal{E}_{h}^{i}$, $|\fint_{e}\undertilde{\psi}\cdot\mathbf{n}_{e}|=|\alpha_{a_{e}^{L}}-\alpha_{a_{e}^{R}}|$, where $a_{e}^{L}$ and $a_{e}^{R}$ are the two ends of $e$; thus $\alpha_{a_{e}^{L}}=\alpha_{a_{e}^{R}}$ for every $e\in\mathcal{E}_{h}^{i}$. Since $\alpha_{a}=0$ for $a\in\mathcal{X}_{h}^{b}$, $\alpha_{a}=0$ for $a\in\mathcal{X}_{h}^{b,+1}$; recursively, we obtain $\alpha_{a}=0$ for $a\in\mathcal{X}_{h}^{b,+j}$ level by level, and finally $\alpha_{a}=0$ for $a\in\mathcal{X}_{h}$.

The proof is completed by noting the two sides of (42) have the same dimension. ∎

For $a\in\mathcal{X}_{h}^{b}\setminus\mathcal{X}^{c}_{h}$, denote by $\mathbf{n}^{b}_{a}$ the outward unit normal vector of $\partial\Omega$ at $a$. Thus, for $a\in\mathcal{X}_{h}^{b}\setminus\mathcal{X}^{c}_{h}$, denote

[TABLE]

Lemma 23**.**

$\undertilde{S}{}^{2}_{ht}({\rm rot},w0)=\undertilde{S}{}^{2}_{h0}({\rm rot},w0)\oplus{\rm span}\{\undertilde{\varphi}{}_{a}^{b}:a\in\mathcal{X}_{h}^{b}\setminus\mathcal{X}^{c}_{h}\}\oplus{\rm span}\{\undertilde{\varphi}{}_{e}:e\in\mathcal{E}_{h}^{b}\}$.

Proof.

By Lemma 21, $\undertilde{\varphi}{}^{b}_{a}$ with $a\in\mathcal{X}_{h}^{b}\setminus\mathcal{X}^{c}_{h}$ are linearly independent, and the right hand side is a direct sum included in the left hand side. On the other hand,

[TABLE]

This proves the assertion. ∎

4.1.2. Structure of piecewise rot-free space

Denote $\undertilde{P}{}_{2}({\rm rot},0)=\{\undertilde{p}\in\undertilde{P}{}_{2}(T):{\rm rot}\undertilde{p}=0\}$, and $\dim(\undertilde{P}{}_{2}({\rm rot},0))=9$. Denote by $\phi_{T}$ the bubble function $(\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2})-2/3$, and define a mapping $\mathcal{F}_{T}$ from $\undertilde{P}{}_{2}({\rm rot},w0)$ to $\undertilde{P}{}_{2}({\rm rot},0)$ by

[TABLE]

Since $\int_{T}{\rm rot}\undertilde{\eta}=0$ in the formula above, the mapping is well defined. It can be verified that $\mathcal{F}_{T}(\undertilde{\eta}{}^{x}_{a_{1}}+\undertilde{\eta}{}^{x}_{a_{2}}+\undertilde{\eta}{}^{x}_{a_{3}})=0$ and $\mathcal{F}_{T}(\undertilde{\eta}{}^{y}_{a_{1}}+\undertilde{\eta}{}^{y}_{a_{2}}+\undertilde{\eta}{}^{y}_{a_{3}})=0$. A frame of $\undertilde{P}{}_{2}({\rm rot},0)$ is presented in the lemma below.

Lemma 24**.**

Any two $\mathcal{F}_{T}\undertilde{\eta}{}^{x}_{a_{i}}$ among $i=1,2,3$, any two $\mathcal{F}_{T}\undertilde{\eta}{}^{y}_{a_{i}}$ among $i=1,2,3$, any two $\mathcal{F}_{T}\undertilde{\eta}{}_{a_{i}}$ among $i=1,2,3$, and all $\mathcal{F}_{T}\undertilde{\eta}{}_{e_{i}}$ for $i=1,2,3$ form a basis of $\undertilde{P}{}_{2}({\rm rot},0)$.

Denote $\undertilde{G}{}^{2}_{h}({\rm rot},0):=\{\undertilde{v}\in\undertilde{G}{}^{2}_{h}:{\rm rot}_{h}\undertilde{v}=0\}$, $\undertilde{G}{}^{2}_{h0}({\rm rot},0):=\{\undertilde{v}\in\undertilde{G}{}^{2}_{h0}:{\rm rot}_{h}\undertilde{v}=0\}$ and $\undertilde{G}{}^{2}_{ht}({\rm rot},0):=\{\undertilde{v}\in\undertilde{G}{}^{2}_{ht}:{\rm rot}_{h}\undertilde{v}=0\}$. Define an operator $\mathcal{F}_{h}:\undertilde{S}{}^{2}_{h}({\rm rot},w0)\to\undertilde{G}{}^{2}_{h}({\rm rot},0)$ by

[TABLE]

Since $\int_{T}{\rm rot}\undertilde{\varphi}{}_{h}=0$ on any $T$ for $\undertilde{\varphi}{}_{h}\in\undertilde{S}{}^{2}_{h}({\rm rot},w0)$, $\mathcal{F}_{h}$ is well defined. Indeed, $(\mathcal{F}_{h}\undertilde{\varphi}{}_{h})|_{T}=\mathcal{F}_{T}(\undertilde{\varphi}{}_{h}|_{T})$.

Lemma 25**.**

$\mathcal{F}_{h}$* is a bijection between $\undertilde{S}{}^{2}_{h0}({\rm rot},w0):=\{\undertilde{v}\in\undertilde{S}{}^{2}_{h0}:({\rm rot}\undertilde{v},q)=0,\ \forall\,q\in\mathbb{P}^{0}_{h0}\}$ and $\undertilde{G}{}^{2}_{h0}({\rm rot},0):=\{\undertilde{v}\in\undertilde{G}{}^{2}_{h0}:{\rm rot}_{h}\undertilde{v}=0\}$, and a bijection between $\undertilde{S}{}^{2}_{ht}({\rm rot},w0):=\{\undertilde{v}\in\undertilde{S}{}^{2}_{ht}:({\rm rot}\undertilde{v},q)=0,\ \forall\,q\in\mathbb{P}^{0}_{h0}\}$ and $\undertilde{G}{}^{2}_{ht}({\rm rot},0):=\{\undertilde{v}\in\undertilde{G}{}^{2}_{ht}:{\rm rot}_{h}\undertilde{v}=0\}$.*

Proof.

Since $\undertilde{S}{}^{2}_{ht}\cap\undertilde{\mathcal{B}}{}_{h0}=\{0\}$, $\mathcal{F}_{h}$ is an injection on $\undertilde{S}{}^{2}_{ht}({\rm rot},w0)$.

Given $\undertilde{\gamma}{}_{h}\in\undertilde{G}{}^{2}_{ht}({\rm rot},0)$, decompose it to $\undertilde{\gamma}{}_{h}=\undertilde{\gamma}{}_{h}^{1}+\undertilde{\gamma}{}_{h}^{2}$ such that $\undertilde{\gamma}{}_{h}^{1}\in\undertilde{S}{}^{2}_{ht}$ and $\undertilde{\gamma}{}_{h}^{2}\in\undertilde{\mathcal{B}}{}_{h0}$. As ${\rm rot}(\undertilde{\gamma}{}_{h}|_{T})=0$ and $\int_{T}{\rm rot}(\phi_{T},0)=\int_{T}{\rm rot}(0,\phi_{T})=0$ on every cell $T$, $\int_{T}{\rm rot}(\undertilde{\gamma}{}_{h}^{1}|_{T})=0$. Namely $\undertilde{\gamma}{}_{h}^{1}\in\undertilde{S}{}^{2}_{ht}({\rm rot},w0)$. This way $\mathcal{F}_{h}$ is a bijection between $\undertilde{S}{}^{2}_{ht}({\rm rot},w0)$ and $\undertilde{G}{}^{2}_{ht}({\rm rot},0)$.

Similarly we can prove $\mathcal{F}_{h}\undertilde{S}{}^{2}_{h0}({\rm rot},w0)=\undertilde{G}{}^{2}_{h0}({\rm rot},0)$, and the proof is completed. ∎

By Lemmas 22, 22 and 25, we can prove the lemmas below.

Lemma 26**.**

The set $\{\mathcal{F}_{h}\undertilde{\varphi}{}^{x}_{a},\mathcal{F}_{h}\undertilde{\varphi}{}^{y}_{a},\mathcal{F}_{h}\undertilde{\varphi}{}_{P_{a}},\mathcal{F}_{h}\undertilde{\varphi}{}_{e}\}_{a\in\mathcal{X}^{i}_{h},\ e\in\mathcal{E}^{i}_{h}}$ forms a basis of $\undertilde{G}{}^{2}_{h0}({\rm rot},0)$.

Lemma 27**.**

$\undertilde{G}{}^{2}_{ht}({\rm rot},0)=\undertilde{G}{}^{2}_{h0}({\rm rot},0)\oplus{\rm span}\{\mathcal{F}_{h}\undertilde{\varphi}{}_{a}^{b}:a\in\mathcal{X}_{h}^{b}\setminus\mathcal{X}^{c}_{h}\}\oplus{\rm span}\{\mathcal{F}_{h}\undertilde{\varphi}{}_{e}:e\in\mathcal{E}_{h}^{b}\}$.

4.1.3. Locally supported basis functions of $B^{3}_{h0}$ and $B^{3}_{ht}$

Now we are going to show that $B^{3}_{h0}$ admits a set of basis functions with vertex-patch-based supports.

Theorem 28**.**

The space $B^{3}_{h0}$ admits a set of basis functions each is supported in a patch of some vertex.

Proof.

By the exact sequence (7), we got to know that the piecewise gradient $\nabla_{h}$ is a bijection between $B^{3}_{h0}$ and $\undertilde{G}{}^{2}_{h0}({\rm rot},0)$. Further, by Lemma 26, the set

[TABLE]

form a basis of $B^{3}_{h0}$. Note again that both $(\nabla_{h})^{-1}$ and $\mathcal{F}_{h}$ preserve the locality of the support; this is verified by viewing the patch as a specific triangulation. Namely (45) is a basis each supported in the patch of a vertex. The proof is completed. ∎

Similar to Theorem 28, we have the description below.

Theorem 29**.**

The space $B^{3}_{ht}$ admits a set of basis functions each is supported in a patch of some vertex.

Proof.

By the complex (8), again, $\nabla_{h}$ is a bijection between $B^{3}_{ht}$ and $\undertilde{G}{}^{2}_{ht}({\rm rot},0)$. By Lemma 27,

[TABLE]

and a locally supported basis of $B^{3}_{ht}$ follows. The proof is completed. ∎

We use the notation below for convenience:

[TABLE]

[TABLE]

[TABLE]

We remark here all these $w^{\prime}s$ can be obtained by straightforward calculation, as, again, both $(\nabla_{h})^{-1}$ and $\mathcal{F}_{h}$ preserve the locality of the supports and can be done cell by cell. Though the space $B^{3}_{h}$ does not correspond to a finite element defined by Ciarlet’s triple, these $w^{\prime}s$ play the same role as that by the usual nodal basis functions. Substituting these functions into the common routine generates finite element codes of the schemes (13) and (14) in a standard way.

4.2. Implementation by decomposition

In this subsection, alternatively, we suggest a decomposition procedure, and the schemes (13) and (14) can be implemented without the explicit construction of the basis functions.

Lemma 30**.**

Let $u_{h}^{*}$ be obtained by the following procedure:

1. (1)

find $r_{h}\in A^{3}_{h0}$ such that

[TABLE] 2. (2)

with $r_{h}$ obtained, find $(\undertilde{\varphi}{}_{h},p_{h})\in\undertilde{G}{}^{2}_{h0}\times\mathbb{P}^{1}_{h0}$ such that

[TABLE] 3. (3)

with $\undertilde{\varphi}{}_{h}$ obtained, find $u_{h}^{*}\in A^{3}_{h0}$ such that

[TABLE]

Let $u_{h}$ be the solution of (13). Then, $u_{h}^{*}=u_{h}$.

Lemma 31**.**

Let $z_{h}^{*}$ be obtained by the following procedure:

1. (1)

find $r_{h}\in A^{3}_{h0}$ such that

[TABLE] 2. (2)

with $r_{h}$ obtained, find $(\undertilde{\varphi}{}_{h},p_{h})\in\undertilde{G}{}^{2}_{ht}\times\mathbb{P}^{1}_{h0}$ such that

[TABLE] 3. (3)

with $\undertilde{\varphi}{}_{h}$ obtained, find $z_{h}^{*}\in A^{3}_{h0}$ such that

[TABLE]

Let $z_{h}$ be the solution of (14). Then, $z_{h}^{*}=z_{h}$.

Lemmas 30 and 31 follows from Theorem 10 and Lemma 7. The scheme (49) is not a convergent one for the Poisson equation, but it is well-posed based on the continuity of $A^{3}_{h0}$ on vertices. With the formulations presented in Lemmas 30 and 31, the spaces used for Poisson equations and Stokes problems only are easy to formulate; to solve the system only needs solving two Poisson systems and one Stokes systems one by one, each of which can be solved with various optimal solvers in a friendly way.

Remark 32**.**

The decompositions as in Lemmas 30 and 31 can also be established for (39) and the one in Remark 20, respectively.

5. Conclusion and discussion

In this paper, based on theoretical analysis by an indirect approach, a constructive answer is given to the question if an optimal scheme can be designed for the biharmonic equation with piecewise cubic polynomials on general triangulations; the schemes work optimally, e.g., on triangulations shown in Figure 1. Beside the theoretical meaning, the scheme can find its application onto practical problems. For example, a high order scheme has been implemented based on $B^{3}_{h}$ for the Helmholtz transmission eigenvalue problem from inverse scattering and accoustics in [41]; we refer there for many numerical experiments about schemes with $B^{3}_{h}$. The practical usage of the scheme can be thus illustrated.

This paper relies on construction and utilization of discretized Stokes complexes based on the $\undertilde{G}{}^{2}_{h0}-\mathbb{P}^{1}_{h0}$ pair. The space $\undertilde{G}{}^{k}_{h}$ with $k=3$ corresponds to the Crouzeix–Falk pair studied in [11]. In that paper, the authors proved that the pair $\undertilde{G}{}^{3}_{h0}-\mathbb{P}^{2}_{h0}$ is stable “for most reasonable meshes.” Moreover, they presented a conjecture that the pair is stable “for any triangulation of a convex polygon satisfying the minimal angle condition and containing an interior vertex.” Recently, some triangulations where $\undertilde{G}{}^{3}_{h0}-\mathbb{P}_{h0}^{2}$ is stable or at least ${\rm div}\undertilde{G}{}^{3}_{h0}=\mathbb{P}_{h0}^{2}$ are introduced in [22]. This hints the possibility to generalize the concept for optimal quartic element schemes (see [48] for details).

It is worthwhile pointing out, in this paper, we focus on the primal schemes only. There have been various kinds of schemes that considered new variables and/or conduct the second order differentiation in a dual way, such as the mixed element method, local DG method, hybridized DG method, CDG method, weak Galerkin method, and so forth. We remark that the literature on related works in this context is vast, but we will not discuss them in this paper. Moreover, based on the space $B^{3}_{h0}$($B^{3}_{ht}$), DG schemes can be designed. One may be able to construct, for example, a weakly over-penalized IP method (like [7]) or IPDG method with optimal convergence rate robust with respect to the penalization paremeter ([49]) with piecewise cubic polynomials.

The spaces $A^{3}_{h}$ and $B^{3}_{h}$ each belongs to a systematic family which reads:

[TABLE]

and

[TABLE]

The spaces $A^{k}_{h0}$ and $B_{h0}^{k}$ can be defined corresponding to the boundary conditions of $H^{1}_{0}(\Omega)$ and $H^{2}_{0}(\Omega)$, respectively. It is now known that $B_{h(0)}^{k}$ is an optimally consistent finite element space for biharmonic equation ($k=2,3$) for arbitrary triangulations. For $k=4$, as discussed above, the assertion holds on most “reasonable” triangulations. Can the family work optimally with arbitrary $k\geqslant 2$ and can it be generalized to a higher dimension and even higher-order problems? This question could be of interest in future research.