¹

1 Introduction↩︎

In this work, we develop superconvergent and divergence-free finite element methods for the Stokes equation on a bounded domain \(\Omega \subset \mathbb{R}^d\) (\(d \geq 2\)) with external force \(\boldsymbol{f}\in L^{2}(\Omega;\mathbb{R}^d)\): \[\label{stokes} \begin{cases} -\Delta \boldsymbol{u}- \nabla p=\boldsymbol{f} & \text{in } \Omega, \\ \operatorname{div}\boldsymbol{u} =0 & \text{in } \Omega, \\ \boldsymbol{u} =0 & \text{on } \partial \Omega, \end{cases}\tag{1}\] where \(\boldsymbol{u}\) and \(p\) denote the velocity and pressure.

Classical stable finite element pairs for 1 include the Taylor–Hood element [1], the MINI element [2], and the nonconforming \(P_1\)–\(P_0\) element [3]. More examples and a detailed discussion are given in [4]. However, these classical pairs are not pointwise divergence-free, which weakens mass conservation and reduces pressure robustness in error estimates. A review of the importance of divergence-free constraints and their numerical implications can be found in [5].

Divergence-free finite element pairs include the Scott–Vogelius element [6], which requires special mesh conditions; smooth finite element pairs [7]–[9] with supersmooth degrees of freedom; conforming pairs on split meshes [10]–[15]; conforming rational pairs [16], [17]; and low-order, normal-continuous nonconforming pairs [18]–[20]. We also refer to [21]–[23] for discontinuous Galerkin methods, and to [24]–[27] for virtual element methods.

One can impose the divergence-free condition by using an \(H(\operatorname{div})\)-conforming finite element space for the velocity. Specifically, we use the Raviart–Thomas (RT) element space [28], [29] or the Brezzi–Douglas–Marini (BDM) element space [30], [31]: \[\mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}:=\{\boldsymbol{v}_h\in H_{0}(\operatorname{div}, \Omega): \boldsymbol{v}_h|_T\in \mathbb{P}_{k}(T;\mathbb{R}^d) + \mathbb{H}_{\ell}(T)\boldsymbol{x} \;\text{ for } T\in\mathcal{T}_h\}\] for the velocity, with integers \(\ell=k\) (\({\rm RT}_k\)) or \(\ell=k-1\) (\({\rm BDM}_{k}\)), and the discontinuous piecewise polynomial space \(\mathbb{P}_{\ell}(\mathcal{T}_h)\) for the pressure.

Since \(\mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\not\subset H^1(\Omega;\mathbb{R}^d)\), the discretization of \(-\Delta \boldsymbol{u}\) is nontrivial. One approach is to use discontinuous Galerkin (DG) methods for the Laplacian operator, as in [22], [32]. The DG method in [22] introduces a penalty term to ensure stability, while the hybridizable DG (HDG) method in [32] employs \(k\)th-order polynomial spaces for all variables. In fact, it suffices to take \(p_h\in\mathbb{P}_{k-1}(\mathcal{T}_h)\) when the velocity \(\boldsymbol{u}_h\in {\rm BDM}_{k}\).

In the vorticity–velocity–pressure formulation [33], [34], the vorticity \(\omega = {\rm curl\,}\boldsymbol{u}\) is introduced, and \[-\Delta \boldsymbol{u} = {\rm curl\,}\omega -{\rm grad\,}\operatorname{div}\boldsymbol{u} = {\rm curl\,}{\rm curl\,}\boldsymbol{u} - {\rm grad\,}\operatorname{div}\boldsymbol{u}\] is interpreted as the Hodge Laplacian. However, as observed in [33], [34], for no-slip Dirichlet boundary conditions (\(\boldsymbol{u}=0\) on \(\partial \Omega\)), there are stability issues and reduced convergence rates. This was later explained in [35]: the suboptimal convergence rate comes from the lack of a proper Hilbert complex in this setting.

The lost convergence rate can be partly recovered by mesh symmetry. In the triangular MAC (TMAC) schemes for the Stokes equation [36], mass lumping is applied to the vorticity so that it can be eliminated, and the velocity and pressure convergence rates are improved by mesh symmetry. Using the element pair \(\textrm{RT}_0\)–\(\textrm{P}_0\) for velocity and pressure (i.e., \(k=\ell=0\)) yields \[|I_{0,0}^{\operatorname{div}}\boldsymbol{u} - \boldsymbol{u}_h|_{1,h} + \|Q_{0} p - p_h\| \lesssim h^{\min\{1,s\}}|\ln h|^{1/2},\] where \(s \geq \tfrac{1}{2}\) measures mesh symmetry and \(I^{\rm \operatorname{div}}_{k,\ell}\) is the interpolation operator into \(\mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\). With the pair \(\textrm{BDM}_1\)–\(\textrm{P}_0\) (i.e., \(k=1\), \(\ell=0\)), the TMAC scheme achieves \[|I_{1,0}^{\operatorname{div}}\boldsymbol{u} - \boldsymbol{u}_h|_{1,h} + \|Q_{0} p - p_h\| \lesssim h^{1+\min\{1/2,s\}}.\]

In [37], [38], the mass-conserving mixed formulation with stresses (MCS) finite element method was developed for the pseudostress–velocity–pressure formulation, also employing \(H(\operatorname{div})\)-conforming velocity spaces. The stress space in [38] is slightly nonconforming, leading to the standard \(h^k\) convergence instead of the superconvergent \(h^{k+1}\) rate. In [37], a weak symmetry constraint was introduced and the stress space was enriched, yielding \(h^{k+1}\) convergence for the \(\mathrm{RT}_k\)–\(\mathrm{P}_k\) pair with \(k \ge 1\).

A divergence free weak virtual element method was developed in [25], which requires stabilization and does not yield superconvergence. The main idea is to introduce a piecewise vector-valued polynomial \(\boldsymbol{\lambda}_h \in \mathring{\Lambda}_k:=\mathbb{P}_{k}(\mathring{\mathcal{F}}_h;\mathbb{R}^{d-1})\) defined on faces to enforce tangential continuity. We follow this approach and improve it to a stabilization-free scheme with superconvergence. Define the weak deviatoric gradient operator element-wise by \[(\operatorname{dev}{\rm grad\,}_w(\boldsymbol{v}, \boldsymbol{\mu}), \boldsymbol{\tau})_T =-(\boldsymbol{v},\operatorname{div}\boldsymbol{\tau})_T +(\boldsymbol{n}\cdot\boldsymbol{v},\boldsymbol{n}^{\intercal}\boldsymbol{\tau}\boldsymbol{n})_{\partial T} +(\boldsymbol{\mu}, \Pi_F\boldsymbol{\tau}\boldsymbol{n})_{\partial T}\] for any \(\boldsymbol{\tau}\in \mathbb{P}_{k}(T;\mathbb{T})\) and \(T\in\mathcal{T}_h\). The tensor space is reduced to the traceless tensor \(\mathbb{T}\), since for \(\boldsymbol{\sigma} := {\rm grad\,}\boldsymbol{u}\) with \(\operatorname{div}\boldsymbol{u} = 0\), we have \(\operatorname{tr}\boldsymbol{\sigma} = \operatorname{tr}{\rm grad\,}\boldsymbol{u} = \operatorname{div}\boldsymbol{u} = 0\).

We propose the following mixed finite element method: find \(\boldsymbol{u}_h\in \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\), \(\boldsymbol{\lambda}_h\in \mathring{\Lambda}_k\), and \(p_h\in\mathbb{P}_{\ell}(\mathcal{T}_h)/\mathbb{R}\) with integers \(k \geq 0\) and \(\ell\in\{k, k-1\}\), such that \[\tag{2} \begin{align} \tag{3} (\operatorname{dev}{\rm grad\,}_w(\boldsymbol{u}_h,\boldsymbol{\lambda}_h), \operatorname{dev}{\rm grad\,}_w(\boldsymbol{v}_h,\boldsymbol{\mu}_h)) + (\operatorname{div}\boldsymbol{v}_h, p_h) &= (\boldsymbol{f}, \boldsymbol{v}_h), \\ \tag{4} (\operatorname{div}\boldsymbol{u}_h, q_h)&=0 \end{align}\] for all \(\boldsymbol{v}_h\in \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\), \(\boldsymbol{\mu}_h\in\mathring{\Lambda}_k\), and \(q_h\in\mathbb{P}_{\ell}(\mathcal{T}_h)/\mathbb{R}\). Notably, the discrete method 2 requires no additional stabilization, despite relying on the weak differential operator \(\operatorname{dev}{\rm grad\,}_w\).

We derive superconvergent estimates for the velocity and pressure: \[\begin{align} \|\boldsymbol{\sigma} - \boldsymbol{\sigma}_h\| + \|\operatorname{dev}{\rm grad\,}_w (I_{k,k}^{\operatorname{div}}\boldsymbol{u} - \boldsymbol{u}_h) \| + \|Q_{\ell} p - p_h\| \lesssim h^{k+1} |\boldsymbol{\sigma}|_{k+1}, \end{align}\] where \(\boldsymbol{\sigma} = \operatorname{dev}{\rm grad\,}\boldsymbol{u},\) and \(\boldsymbol{\sigma}_h:= \operatorname{dev}{\rm grad\,}_w(\boldsymbol{u}_h, \boldsymbol{\lambda}_h)\). Unlike the TMAC scheme [36], such superconvergence does not require the symmetry of the triangulation. It is further utilized to construct a new superconvergent approximation of \(\boldsymbol{u}\).

Our method allows \(\ell = 0\) an important case not covered by MCS [37], [38]. We provide an illustration of the case \((k, \ell) = (0, 0)\), which can be treated as a generalization of the classical MAC scheme [39] on rectangular grids to triangulations and achieve first order \(h\) convergence, no loss of convergence compared with the vorticity-velocity-pressure formulation. One can use \(\textrm{BDM}_1\)–\(\textrm{P}_0\) pair, \((k, \ell) = (1, 0)\), to get a second order \(h^2\) scheme.

We propose a hybridization of 2 by relaxing the normal continuity of the velocity, which enables the choice \((k,\ell)=(0,-1)\). The velocity is approximated by the piecewise constant vector space, and the normal continuity is enforced through a nonconforming linear element space for the pressure \(p_h\). Even in this low-order setting, the method achieves at least first-order convergence, which is a notable feature.

The remainder of this paper is organized as follows. Section 2 discusses the weak div stability and the distributional discretization of the vector Laplacian. Section 3 develops mixed finite element methods for the Stokes equation, and Section 4 presents their equivalent discrete formulations. Finally, Section 5 reports numerical experiments that confirm and illustrate the theoretical results.

2 Distributional Discretization of the Vector Laplacian↩︎

In this section, we discretize the weak divergence operator \[\operatorname{div}: L^2(\Omega; \mathbb{T}) \to H^{-1}(\Omega; \mathbb{R}^d),\] and apply it to the vector Laplacian problem. We introduce weak divergence and deviatoric gradient operators, establish their inf-sup stability, and derive a superconvergent error estimate for divergence-free velocity spaces.

2.1 Notation↩︎

Let \(\Omega \subset \mathbb{R}^d\) (\(d \ge 2\)) be a bounded polytope with boundary \(\partial\Omega\). For a bounded domain \(D\) and integer \(m \ge 0\), let \(H^m(D)\) be the Sobolev space with norm \(\|\cdot\|_{m,D}\) and semi-norm \(|\cdot|_{m,D}\), and let \(H_0^m(D)\) be the closure of \(C_0^\infty(D)\) under \(\|\cdot\|_{m,D}\). Set \(L^2(D)=H^0(D)\) with inner product \((\cdot,\cdot)_D\), and write \(\|\cdot\|_{0,D}\) as \(\|\cdot\|_D\). When \(D=\Omega\), the subscript \(D\) is omitted. Let \(h_D\) denote the diameter of \(D\), and \(\boldsymbol{n}_{\partial D}\) the unit outward normal, abbreviated as \(\boldsymbol{n}\).

We use \[\begin{align} H(\operatorname{div},D) &=\{\boldsymbol{v}\in L^2(D;\mathbb{R}^d):\operatorname{div}\boldsymbol{v}\in L^2(D)\},\\ H_0(\operatorname{div},D)&=\{\boldsymbol{v}\in H(\operatorname{div},D):\boldsymbol{v}\cdot\boldsymbol{n}=0 \text{ on }\partial D\}. \end{align}\] Let \(L_0^2(D)=\{q\in L^2(D):\int_D q\,\,{\rm d}x=0\}.\) For any \(B\subseteq L^2(D)\), denote \(B/\mathbb{R}=B\cap L_0^2(D)\).

Let \(\{\mathcal{T}_h\}_{h>0}\) be a regular family of simplicial meshes of \(\Omega\), where \(h = \max_{T \in \mathcal{T}_h} h_T\) and \(h_T = \operatorname{diam}(T)\). Denote by \(\mathcal{F}_h\) and \(\mathring{\mathcal{F}}_h\) the sets of all faces and interior faces, respectively. For each face \(F\in\mathcal{F}_h\), fix a unit normal \(\boldsymbol{n}_F\). For \(T\in\mathcal{T}_h\), let \(\mathcal{F}(T)\) be the set of all \((d-1)\)-dimensional faces of \(T\).

For a face \(F\) and vector \(\boldsymbol{v}\in\mathbb{R}^d\), define its tangential projection \[\Pi_F\boldsymbol{v} = (\boldsymbol{I} - \boldsymbol{n}_F\boldsymbol{n}_F^{\intercal})\boldsymbol{v}.\] For two adjacent elements \(T^{\pm}\) sharing \(F\), define the jump \[\llbracket v \rrbracket|_F = v^+\boldsymbol{n}_F\cdot\boldsymbol{n}_{\partial T^{+}} + v^-\boldsymbol{n}_F\cdot\boldsymbol{n}_{\partial T^{-}},\] and on a boundary face \(F\subset\partial\Omega\), set \(\llbracket v \rrbracket = v|_F\,\boldsymbol{n}_F\cdot\boldsymbol{n}\).

For a domain \(D\) and integer \(k\ge0\), let \(\mathbb{P}_k(D)\) be the space of polynomials on \(D\) of degree \(\le k\), \(\mathbb{H}_k(D) = \mathbb{P}_k(D)\setminus\mathbb{P}_{k-1}(D)\), and \(Q_{k,D}\) the \(L^2\)-projection onto \(\mathbb{P}_k(D)\). Set \(\mathbb{M} = \mathbb{R}^{d\times d}\), \(\mathbb{K}\) for skew-symmetric, and \(\mathbb{T}\) for traceless matrices. For \(\boldsymbol{\tau}\in\mathbb{M}\), \[\operatorname{dev}\boldsymbol{\tau} = \boldsymbol{\tau} - \tfrac{1}{d}(\operatorname{tr}\boldsymbol{\tau})\boldsymbol{I}, \qquad \operatorname{skw}\boldsymbol{\tau} = \tfrac{1}{2}(\boldsymbol{\tau}-\boldsymbol{\tau}^{\intercal}).\]

Define \[\begin{align} &H^s(\mathcal{T}_h) = \{v\in L^2(\Omega): v|_T\in H^s(T)\} = \prod_{T\in \mathcal{T}_h} H^s(T),\\ &\mathbb{P}_k(\mathcal{T}_h) = \{v\in L^2(\Omega): v|_T\in\mathbb{P}_k(T)\} = \prod_{T\in \mathcal{T}_h} \mathbb{P}_{k}(T),\\ &\mathbb{P}_k(\mathcal{F}_h) = \{v\in L^2(\mathcal{F}_h): v|_F\in\mathbb{P}_k(F)\} = \prod_{F\in \mathcal{F}_h} \mathbb{P}_{k}(F), \end{align}\] and similarly define \(\mathbb{P}_k(\mathring{\mathcal{F}}_h)\) for interior faces only. For a space \(B(D)\) or \(B(\mathcal{T}_h)\), define \[B(D; \mathbb{X}) = B(D) \otimes \mathbb{X}, \quad B(\mathcal{T}_h; \mathbb{X}) = B(\mathcal{T}_h) \otimes \mathbb{X},\] where \(\mathbb{X} \in \{\mathbb{R}^d, \mathbb{M}, \mathbb{T}, \mathbb{K}\}\). Let \(Q_{k}\) denote the \(L^2\)-projection operator from \(L^2(\Omega)\) onto \(\mathbb{P}_k(\mathcal{T}_h)\); its vector- or tensor-valued version is also denoted by \(Q_{k}\). We use \({\rm grad\,}_h\) (\(\nabla_h\) for scalar function) and \(\operatorname{div}_h\) to represent the element-wise gradient and divergence operators, respectively.

2.2 Weak div stability for the scalar Laplacian↩︎

We first consider the surjective map \[\operatorname{div}: L^{2}(\Omega;\mathbb{R}^d) \to H^{-1}(\Omega),\] and discuss its distributional finite element discretization. This operator is part of the distributional de Rham complex introduced in [40].

The \(L^2\) space is discretized by the discontinuous polynomial space \[\Sigma_k^{-1}(\mathcal{T}_h;\mathbb{R}^d) = \mathbb{P}_k(\mathcal{T}_h; \mathbb{R}^d),\] which we abbreviate as \(\Sigma_k^{-1}\). A natural norm on \(\Sigma_k^{-1}\) is the standard \(L^2\) norm \(\|\cdot\|\).

Next, we introduce the product space of discontinuous polynomial spaces on \(\mathcal{T}_h\times \mathcal{F}_h\): \[M_{k-1,k}^{-1} := \mathbb{P}_k(\mathcal{T}_h) \times \mathbb{P}_k(\mathcal{F}_h).\] An element of \(M_{k-1,k}^{-1}\) is denoted by \(u = (u_0, u_b)\). A natural product-type inner product is defined by \[((u_0, u_b), (v_0, v_b))_{0,h} = (u_0, v_0) + \sum_{F\in\mathcal{F}_h} h_F (u_b, v_b)_{F},\] where the scaling factor \(h_F\) ensures that the two terms are of comparable magnitude. The corresponding norm on \(M_{k-1,k}^{-1}\) is \[\|u\|_{0,h}^2 = \|u_0\|^2 + \sum_{F\in \mathcal{F}_h} \|h_F^{1/2} u_b\|_{F}^2.\] The subspace \[\mathring{M}_{k-1,k}^{-1} = \mathbb{P}_k(\mathcal{T}_h) \times \mathbb{P}_k(\mathring{\mathcal{F}}_h)\] can be viewed as a subspace of \({M}_{k-1,k}^{-1}\) by setting the values on boundary faces to zero.

For \(\sigma \in \Sigma_k^{-1}\), define the operator \(\operatorname{div}_w: \Sigma_k^{-1}\to \mathring{M}_{k-1,k}^{-1}\) by \[\operatorname{div}_w \sigma = \big(\operatorname{div}_T \sigma, - h_F^{-1}[\sigma \cdot \boldsymbol{n}]|_F\big)_{T\in \mathcal{T}_h,~F\in \mathring{\mathcal{F}}_h} \in \mathring{M}_{k-1,k}^{-1}.\] When \(\sigma \in \Sigma_k^{-1}\cap H(\operatorname{div},\Omega)\), we have \(\operatorname{div}_w\sigma = (\operatorname{div}\sigma,0)\) since \([\sigma \cdot \boldsymbol{n}]|_F = 0\) for all interior faces \(F\); in this case we simply write \(\operatorname{div}_w\sigma = \operatorname{div}\sigma\).

A graph norm associated with \(\operatorname{div}_w\) is defined as \[\|\sigma\|_{\operatorname{div}_w}^2 := \|\sigma\|^2 + \|\operatorname{div}_w\sigma\|_{0,h}^2 = \|\sigma\|^2 + \|\operatorname{div}_h\sigma\|^2 + \|h^{-1/2}[\sigma \cdot \boldsymbol{n}]\|_{\mathring{\mathcal{F}}_h}^2,\] where the scaling is chosen so that the last two terms are of comparable order. For \(\sigma \in \Sigma_k^{-1}\cap H(\operatorname{div},\Omega)\), we recover \(\|\sigma\|_{\operatorname{div}_w} = \|\sigma\|_{\operatorname{div}}\).

The following inf-sup condition can be established in a manner similar to Theorem 4, and hence the proof is omitted.

Consider the mixed formulation of the Poisson equation \(-\Delta u = f\) with Dirichlet boundary condition \(u|_{\partial \Omega} = 0\). Find \(\sigma \in \Sigma_k^{-1}\) and \(u\in \mathring{M}_{k-1,k}^{-1}\) such that \[\label{eq:hy} \begin{align} (\sigma, \tau) &= (\nabla_w u, \tau) \quad \forall\,\tau \in \Sigma_k^{-1},\\ (\operatorname{div}_w \sigma, v) &= -(f, v_0) \quad\; \forall\,v\in \mathring{M}_{k-1,k}^{-1}. \end{align}\tag{5}\] This formulation is precisely the hybridized mixed finite element method using the \({\rm BDM}_k\)–\(\mathbb{P}_{k-1}\) pair [41]. The problem 5 is well-posed by the inf-sup condition ?? .

2.3 Finite element spaces↩︎

In the Stokes problem, we use the weak divergence operator to discretize the vector Laplacian. Let \(\boldsymbol{\sigma} = {\rm grad\,}\boldsymbol{u}\) and \(-\Delta \boldsymbol{u} = -\operatorname{div}\boldsymbol{\sigma}\). To distinguish the vector formulation from the scalar case, we use boldface symbols such as \(\boldsymbol{\sigma}\) and \(\boldsymbol{u}\).

By the tensor-product construction, we obtain the weak divergence operator \[\operatorname{div}_w: \Sigma_k^{-1}(\mathbb{M}) \to \mathring{M}_{k-1,k}^{-1}(\mathbb{R}^d),\] where \[\Sigma_k^{-1}(\mathbb{X}):=\prod_{T\in \mathcal{T}_h} \mathbb{P}_{k}(T; \mathbb{X})\quad\text{for }\mathbb{X}=\mathbb{M} \text{ or } \mathbb{T};\quad \mathring{M}_{k-1,k}^{-1}(\mathbb{R}^d):=\mathbb{R}^d\otimes\mathring{M}_{k-1,k}^{-1}.\] The tensor space is restricted to the traceless subspace \(\mathbb{T}\), since for \(\boldsymbol{\sigma} = {\rm grad\,}\boldsymbol{u}\) with \(\operatorname{div}\boldsymbol{u} = 0\), we have \(\operatorname{tr}\boldsymbol{\sigma} = \operatorname{div}\boldsymbol{u} = 0\). We then impose suitable continuity conditions on the range space and identify the subspaces of \(\Sigma_k^{-1}(\mathbb{T})\) that lead to weakly divergence-stable pairs.

Spaces for the velocity↩︎

We employ \(H(\operatorname{div})\)-conforming finite elements to discretize the velocity field \(\boldsymbol{u}\). For a simplex \(T\) and integers \(k \ge 0\) and \(\ell \ge -1\), the local space of shape functions is defined by \[\mathbb{V}_{k,\ell}^{\operatorname{div}}(T) := \mathbb{P}_{k}(T;\mathbb{R}^d) + \mathbb{H}_{\ell}(T)\boldsymbol{x}, \quad \ell = k \text{ or } k - 1.\] For \(\ell \ge 0\), the degrees of freedom (DoFs) are given by (cf. [42]) \[\tag{6} \begin{align} (\boldsymbol{v}\cdot\boldsymbol{n}, q)_F, & \quad q\in\mathbb{P}_{k}(F), \quad F\in\mathcal{F}(T), \tag{7}\\ (\boldsymbol{v}, \boldsymbol{q})_T, & \quad \boldsymbol{q}\in {\rm grad\,}\mathbb{P}_{\ell}(T)\oplus \{\boldsymbol{q}\in \mathbb{P}_{k-1}(T;\mathbb{R}^d): \boldsymbol{q}\cdot\boldsymbol{x}=0\}. \tag{8} \end{align}\] When \(\ell = k\), the space corresponds to the Raviart–Thomas (RT) element [28], [29], where the interior moments 8 span \(\mathbb{P}_{k-1}(T;\mathbb{R}^d)\). For \(\ell = k-1\), the Brezzi–Douglas–Marini (BDM) element [30], [31], the interior moments 8 forms a strict subspace.

The corresponding global finite element spaces are defined as \[\begin{align} \mathbb{V}_{k,\ell}^{\operatorname{div}}(\mathcal{T}_h) &:= \{\boldsymbol{v}_h\in H(\operatorname{div},\Omega): \boldsymbol{v}_h|_T \in \mathbb{V}_{k,\ell}^{\operatorname{div}}(T) \text{ for all } T\in\mathcal{T}_h\},\\ \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}(\mathcal{T}_h) &:= \mathbb{V}_{k,\ell}^{\operatorname{div}}(\mathcal{T}_h)\cap H_0(\operatorname{div},\Omega). \end{align}\] For simplicity, we omit \(\mathcal{T}_h\) in the notation afterwards.

To study the weak div stability, we embed \(\mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\) into \(\mathring{M}_{k-1,k}^{-1}(\mathbb{R}^d)\).

The embedding \(E\) can be extended to \(H^1(\Omega;\mathbb{R}^d)\), but the injectivity no longer holds.

Spaces for the pseudo-stress↩︎

For an integer \(k \ge 0\), the local space of shape functions is \(\mathbb{P}_{k}(T;\mathbb{T})\), with degrees of freedom (DoFs) defined by \[\tag{9} \begin{align}\tag{10} \int_{F} \boldsymbol{t}_i^{\intercal}\boldsymbol{\tau}\boldsymbol{n}\, q \, \mathrm{d}S, & \quad q \in \mathbb{P}_{k}(F), \;i = 1,2,\ldots, d-1,\;F\in\mathcal{F}(T), \\ \tag{11} \int_{T} \boldsymbol{\tau}: \boldsymbol{q} \, \mathrm{d}x, & \quad \boldsymbol{q} \in \mathbb{P}_{k-1}(T; \mathbb{T}). \end{align}\] The unisolvence of these DoFs is established in [37], [43]. The corresponding global finite element space is \[\begin{align} \Sigma_{k}^{\rm tn}(\mathcal{T}_h) := \{\boldsymbol{\tau}_h \in \Sigma_k^{-1}(\mathbb{T}) : \text{ all DoFs in~\eqref{eq:curldivdof} are single-valued}\}. \end{align}\] Functions in \(\Sigma_{k}^{\rm tn}\) are tangential-normal continuous, i.e., \([\Pi_F\boldsymbol{\sigma}\boldsymbol{n}]_F = 0\) for all \(F \in \mathring{\mathcal{F}}_h\).

The weak divergence operator \(\operatorname{div}_{w}\) defines a bilinear form on \(\Sigma_k^{\rm tn} \times \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\): \[\begin{align} (\operatorname{div}_{w}\boldsymbol{\sigma}, \boldsymbol{v})_{0,h} := (\operatorname{div}_{w}\boldsymbol{\sigma}, E\boldsymbol{v})_{0,h} &= \sum_{T\in\mathcal{T}_h}(\operatorname{div}\boldsymbol{\sigma}, \boldsymbol{v})_T - \sum_{F\in\mathring{\mathcal{F}}_h}([\boldsymbol{n}^{\intercal}\boldsymbol{\sigma}\boldsymbol{n}], \boldsymbol{n}\cdot\boldsymbol{v})_F\\ &= - \sum_{T\in\mathcal{T}_h}(\boldsymbol{\sigma}, \operatorname{dev}{\rm grad\,}\boldsymbol{v})_T + \sum_{F\in\mathcal{F}_h}(\Pi_F\boldsymbol{\sigma}\boldsymbol{n}, [\Pi_F\boldsymbol{v}])_F. \end{align}\]

2.4 Weak div stability↩︎

For each \(T \in \mathcal{T}_h\), let \(I_T^{\rm tn}: H^1(T;\mathbb{T}) \to \mathbb{P}_{k}(T; \mathbb{T})\) denote the local interpolation operator defined by DoFs 9 . The corresponding global interpolation operator \(I_k^{\rm tn}: H^1(\mathcal{T}_h;\mathbb{T}) \to \Sigma_k^{-1}(\mathbb{T})\) is given by \[(I_k^{\rm tn}\boldsymbol{\tau})|_T := I_T^{\rm tn}(\boldsymbol{\tau}|_T), \quad \forall\, T \in \mathcal{T}_h, \;\boldsymbol{\tau} \in H^1(\mathcal{T}_h; \mathbb{T}).\] Then \(I_k^{\rm tn}\boldsymbol{\tau}\in \Sigma_{k}^{\rm tn}\) for \(\boldsymbol{\tau}\in H^1(\Omega;\mathbb{T})\). For any \(T\in\mathcal{T}_h\) and \(\boldsymbol{\tau}\in H^m(T;\mathbb{T})\) with \(1\le m\le k+1\), the standard interpolation estimate holds: \[\label{eq:Ihtnprop2} \|\boldsymbol{\tau}-I_T^{\rm tn}\boldsymbol{\tau}\|_{T} + h_T|\boldsymbol{\tau}-I_T^{\rm tn}\boldsymbol{\tau}|_{1,T} \lesssim h_T^m|\boldsymbol{\tau}|_{m,T}.\tag{12}\]

We now prove the following important commutative property.

This leads to the following inf-sup condition.

If the velocity space is enriched to \({\rm RT}_k\), i.e., \(\boldsymbol{v}_h \in \mathring{\mathbb{V}}_{k,k}^{\operatorname{div}}\), then \(\operatorname{dev}{\rm grad\,}_h \boldsymbol{v}_h \in \Sigma_k^{-1}(\mathbb{T})\). Since DoF 11 to define \(I_k^{\rm tn}\) is for moments of degree \(k-1\), the commutative property ?? no longer holds. Consequently, the inf-sup condition may fail for the pair \(\Sigma_k^{\rm tn}\)-\(\mathring{\mathbb{V}}_{k,k}^{\operatorname{div}}\).

Define the weak deviatoric gradient \(\operatorname{dev}{\rm grad\,}_w: \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\to \Sigma_k^{\rm tn}\) as the adjoint of \(\operatorname{div}_w\). For \(\boldsymbol{v}\in \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\), \(\operatorname{dev}{\rm grad\,}_w \boldsymbol{v}\in \Sigma_k^{\rm tn}\) is defined by \[(\operatorname{dev}{\rm grad\,}_w \boldsymbol{v}, \boldsymbol{\tau}) = -(\boldsymbol{v}, \operatorname{div}_w \boldsymbol{\tau})_{0,h}, \quad \forall\, \boldsymbol{\tau}\in \Sigma_k^{\rm tn}.\] Since \(\operatorname{div}_{w}\Sigma_k^{\rm tn}=\mathring{\mathbb{V}}_{k,k-1}^{\operatorname{div}}\), the operator \(\operatorname{dev}{\rm grad\,}_w: \mathring{\mathbb{V}}_{k,k-1}^{\operatorname{div}}\to \Sigma_k^{\rm tn}\) is injective, and thus \(\|\operatorname{dev}{\rm grad\,}_w(\cdot)\|\) defines a norm on \(\mathring{\mathbb{V}}_{k,k-1}^{\operatorname{div}}\).

The inf-sup condition ?? can then be written equivalently as \[\label{eq:infsupdevgrad} \inf_{\boldsymbol{v}_h\in \mathring{\mathbb{V}}_{k,k-1}^{\operatorname{div}}} \sup_{\boldsymbol{\tau}_h\in \Sigma_{k}^{\rm tn}} \frac{(\operatorname{div}_w \boldsymbol{\tau}_h, \boldsymbol{v}_h)_{0,h}}{\|\boldsymbol{\tau}_h\|\,\|\operatorname{dev}{\rm grad\,}_w \boldsymbol{v}_h\|} = 1 > 0.\tag{13}\]

Consider the mixed finite element method for the vector Laplacian \[-\Delta \boldsymbol{u} = \boldsymbol{f}, \quad \text{ with } \; \boldsymbol{u}\in H_0^1(\Omega;\mathbb{R}^d).\] Find \(\boldsymbol{\sigma}_h\in \Sigma_k^{\rm tn}\) and \(\boldsymbol{u}_h\in \mathring{\mathbb{V}}_{k,k-1}^{\operatorname{div}}\) such that \[\label{eq:vecLapBDM} \begin{align} (\boldsymbol{\sigma}_h, \boldsymbol{\tau}_h) + (\boldsymbol{u}_h, \operatorname{div}_w \boldsymbol{\tau}_h)_{0,h} &= 0, \qquad\qquad\; \forall\,\boldsymbol{\tau}_h\in \Sigma_k^{\rm tn},\\ (\operatorname{div}_w\boldsymbol{\sigma}_h, \boldsymbol{v}_h)_{0,h} &= -(\boldsymbol{f}, \boldsymbol{v}_h), \quad \forall\,\boldsymbol{v}_h\in \mathring{\mathbb{V}}_{k,k-1}^{\operatorname{div}}. \end{align}\tag{14}\] The formulation 14 is well posed by the inf-sup conditions ?? or 13 . However, the errors \(\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_h\|\) and \(\|\operatorname{dev}{\rm grad\,}\boldsymbol{u}-\operatorname{dev}{\rm grad\,}_w I^{\operatorname{div}}_{k,k-1}\boldsymbol{u}\|\) are coupled. As a result, the dominant contribution is the \(H^1\)-type error of \(\boldsymbol{u}\) of order \(h^k\), rather than the \(L^2\) error of \(\boldsymbol{\sigma}\) of order \(h^{k+1}\). This suboptimality comes from the loss of exact divergence-free conformity: \[(\operatorname{div}_w\boldsymbol{\sigma}_h, \boldsymbol{v}_h)_{0,h} = 0 \;\;\forall\,\boldsymbol{v}_h\in \mathring{\mathbb{V}}_{k,k-1}^{\operatorname{div}} \;\;\nRightarrow\;\; (\operatorname{div}_w\boldsymbol{\sigma}_h, E\boldsymbol{v})_{0,h} = 0 \;\;\forall\,\boldsymbol{v}\in H^1(\Omega;\mathbb{R}^d).\]

2.5 Commuting property for the weak dev grad operator↩︎

For each \(T \in \mathcal{T}_h\), let \(I_{(k,\ell),T}^{\rm div}: H^1(T;\mathbb{R}^d) \to \mathbb{V}_{k,\ell}^{\operatorname{div}}(T)\) denote the local interpolation operator defined by the DoFs in 6 . The corresponding global interpolation operator \(I_{k,\ell}^{\operatorname{div}}: H^1(\mathcal{T}_h;\mathbb{R}^d)\to L^2(\Omega;\mathbb{R}^d)\) is given by \[(I_{k,\ell}^{\operatorname{div}}\boldsymbol{v})|_T := I_{(k,\ell),T}^{\rm div}(\boldsymbol{v}|_T), \quad \forall\, T \in \mathcal{T}_h, \;\boldsymbol{v} \in H^1(\mathcal{T}_h; \mathbb{R}^d).\] Then \(I_{k,\ell}^{\operatorname{div}}\boldsymbol{v}\in \mathbb{V}_{k,\ell}^{\operatorname{div}}\) for any \(\boldsymbol{v}\in H^1(\Omega;\mathbb{R}^d)\). For \(T\in\mathcal{T}_h\) and \(\boldsymbol{v}\in H^m(T;\mathbb{R}^d)\) with \(1\le m\le k+1\), the interpolation satisfies \[\label{eq:Ihdivprop2} \|\boldsymbol{v}-I_{(k,\ell),T}^{\rm div}\boldsymbol{v}\|_{T} + h_T|\boldsymbol{v}-I_{(k,\ell),T}^{\rm div}\boldsymbol{v}|_{1,T} \lesssim h_T^m|\boldsymbol{v}|_{m,T}.\tag{15}\]

We are going to prove the following commutative property: \[\label{eq:Qkgrad} Q_{k}^{\rm tn}\operatorname{dev}{\rm grad\,}= \operatorname{dev}{\rm grad\,}_w I_{k,k}^{\operatorname{div}},\tag{16}\] where \(Q_k^{\rm tn}: L^2(\Omega; \mathbb{T})\to \Sigma_{k}^{\rm tn}\) is the \(L^2\)-projection.

The \(\mathrm{RT}_k\) space is required to ensure the interior moments 8 cover \(\mathbb{P}_{k-1}(T;\mathbb{R}^d)\) for \(\operatorname{div}\boldsymbol{\tau}_h\), while for the \(\mathrm{BDM}_k\) space \(\mathbb{V}_{k,k-1}^{\operatorname{div}}\), the interior moments are not sufficient. Hence, the commutative property ?? does not hold for \(I_{k,k-1}^{\operatorname{div}}\).

The weak operator \(\operatorname{div}_w\) can be regarded as a discretization of \(\operatorname{div}\) in the distributional sense. We can preserve the exact divergence condition in the following sense.

2.6 Restriction to divergence-free subspaces↩︎

The conditions ?? and ?? are incompatible. Weak div stability holds only on the smaller space \(\mathbb{V}_{k,k-1}^{\operatorname{div}}\), while the exact divergence-freeness requires the larger space \(\mathbb{V}_{k,k}^{\operatorname{div}}\). We resolve this conflict by considering the divergence-free subspaces.

We have the following commuting property (cf. [44]): \[\label{eq:Ihdivprop1} \operatorname{div}\!\big(I_{k,\ell}^{\operatorname{div}}\boldsymbol{v}\big) = Q_{\ell}\big(\operatorname{div}\boldsymbol{v}\big), \qquad \forall\,\boldsymbol{v}\in H^1(\Omega;\mathbb{R}^d),\tag{17}\] which implies \(\operatorname{div}\!\big(I_{k,\ell}^{\operatorname{div}}\boldsymbol{v}\big)=0\) for all \(\boldsymbol{v}\in H^1(\Omega;\mathbb{R}^d)\cap\ker(\operatorname{div})\). The additional interior moments in \({\rm RT}_k\) compared with \({\rm BDM}_k\) serve only to enrich the range of the divergence operator. Hence, for any \(\boldsymbol{v}\in H^1(\Omega;\mathbb{R}^d)\cap\ker(\operatorname{div})\) and \(k\ge 1\), \[I_{k,k}^{\operatorname{div}}\boldsymbol{v} = I_{k,k-1}^{\operatorname{div}}\boldsymbol{v} \in \mathbb{V}_{k,k-1}^{\operatorname{div}}\cap \ker(\operatorname{div}).\]

Consider the mixed finite element method for the vector Laplacian \[\label{eq:vectorLapdivfree} -\Delta \boldsymbol{u} = \boldsymbol{f}, \quad \text{with } \boldsymbol{u}\in H_0^1(\Omega; \mathbb{R}^d)\cap \ker(\operatorname{div}).\tag{18}\] Find \(\boldsymbol{\sigma}_h\in \Sigma_k^{\rm tn}\) and \(\boldsymbol{u}_h\in \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\cap \ker(\operatorname{div})\) such that \[\tag{19} \begin{align} \tag{20} (\boldsymbol{\sigma}_h, \boldsymbol{\tau}_h) + (\boldsymbol{u}_h, \operatorname{div}_w \boldsymbol{\tau}_h)_{0,h} &= 0, \qquad\qquad\; \forall\, \boldsymbol{\tau}_h\in \Sigma_k^{\rm tn},\\ \tag{21} (\operatorname{div}_w \boldsymbol{\sigma}_h, \boldsymbol{v}_h)_{0,h} &= -(\boldsymbol{f}, \boldsymbol{v}_h), \quad \forall\, \boldsymbol{v}_h\in \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\cap \ker(\operatorname{div}). \end{align}\]

The velocity space \(\boldsymbol{v}_h|_T\in \mathbb{V}_{k,\ell}^{\operatorname{div}}(T)\) and \(\mathbb{P}_k(T;\mathbb{R}^d)\subseteq \mathbb{V}_{k,\ell}^{\operatorname{div}}(T)\), and the optimal order of convergence for \(H^1\) type error is \(h^k\). So ?? is a superconvergent error estimate for \(\boldsymbol{u}\).

3 Mixed Finite Element Methods for Stokes Equation↩︎

In this section, we introduce the pressure as a Lagrange multiplier to enforce the divergence-free condition, and another multiplier to relax the tangential-normal continuity of the stress space. The resulting mixed finite element methods for the Stokes equation are pointwise divergence-free and superconvergent.

3.1 Spaces of Lagrange multipliers↩︎

The subspace \(\mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\cap \ker(\operatorname{div})\) is not constructed explicitly. Instead, we enforce the constraint \(\operatorname{div}\boldsymbol{u} = 0\) by introducing a Lagrange multiplier, which has the physical interpretation of pressure.

In the discretization, we use the piecewise polynomial space \(\mathbb{P}_{\ell}(\mathcal{T}_h)\) to approximate the pressure \(p\). By the commutative property 17 , it follows that \[\operatorname{div}\mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}=\mathbb{P}_{\ell}(\mathcal{T}_h)/\mathbb{R}.\]

We introduce the Lagrange multiplier spaces \[\Lambda_k = \mathbb{P}_k(\mathcal{F}_h; \mathbb{R}^{d-1}) \quad\textrm{and}\quad \mathring{\Lambda}_k = \mathbb{P}_k(\mathring{\mathcal{F}}_h; \mathbb{R}^{d-1})\] to relax the tangential-normal continuity condition in the space \(\Sigma_k^{\rm tn}\). The embedding \(E: \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}} \times \mathring{\Lambda}_k \to \mathring{M}_{k-1,k}^{-1}(\mathbb{R}^d)\) is extended as \[(\boldsymbol{u}, \boldsymbol{\lambda}) \mapsto \big(Q_{k-1,T}\boldsymbol{u},\, (\boldsymbol{\lambda},\, \boldsymbol{n}_F\cdot \boldsymbol{u})_F\big)_{T\in\mathcal{T}_h,\, F\in\mathring{\mathcal{F}}_h},\] which remains injective.

The stress space is relaxed to \(\Sigma_k^{-1}(\mathbb{T}) = \prod_{T\in \mathcal{T}_h} \mathbb{P}_k(T; \mathbb{T})\). The weak divergence operator \(\operatorname{div}_w\) then induces the bilinear form on \(\Sigma_k^{-1}(\mathbb{T})\times (\mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\times \mathring{\Lambda}_k)\): \[\label{eq:divwT} (\operatorname{div}_w\boldsymbol{\sigma}, (\boldsymbol{v}, \boldsymbol{\mu}))_{0,h} = \sum_{T\in \mathcal{T}_h}(\operatorname{div}\boldsymbol{\sigma}, \boldsymbol{v})_T - \sum_{F\in\mathring{\mathcal{F}}_h}\!\Big( ([\boldsymbol{n}^{\intercal}\boldsymbol{\sigma}\boldsymbol{n}],\, \boldsymbol{n}\!\cdot\!\boldsymbol{v})_F + ([\Pi_F\boldsymbol{\sigma}\boldsymbol{n}],\, \boldsymbol{\mu})_F\Big),\tag{22}\] which induces a mapping, still denoted by \(\operatorname{div}_w: \Sigma_k^{-1}(\mathbb{T})\to (\mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\times \mathring{\Lambda}_k)'\cong (\mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\times \mathring{\Lambda}_k)\) by the \((\cdot,\cdot)_{0,h}\) inner product through the embedding \(E\).

Its adjoint, still denoted by \(\operatorname{dev}{\rm grad\,}_w: \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}} \times \mathring{\Lambda}_k \to \Sigma_k^{-1}(\mathbb{T})\), is defined by \[(\boldsymbol{\sigma}, \operatorname{dev}{\rm grad\,}_w (\boldsymbol{v}, \boldsymbol{\mu})) := -(\operatorname{div}_w\boldsymbol{\sigma}, (\boldsymbol{v}, \boldsymbol{\mu}))_{0,h},\] for \(\boldsymbol{\sigma}\in \Sigma_k^{-1}(\mathbb{T})\), \(\boldsymbol{v}\in \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\), and \(\boldsymbol{\mu}\in \mathring{\Lambda}_k\).

The operator \(\operatorname{dev}{\rm grad\,}_w\) can be extended to \[\operatorname{dev}{\rm grad\,}_w: H^1(\mathcal{T}_h;\mathbb{R}^d)\times L^2(\mathcal{F}_h;\mathbb{R}^{d-1}) \to \Sigma_k^{-1}(\mathbb{T})\] as follows: for \((\boldsymbol{v}, \boldsymbol{\mu})\in H^1(\mathcal{T}_h;\mathbb{R}^d)\times L^2(\mathcal{F}_h;\mathbb{R}^{d-1})\), let \(\operatorname{dev}{\rm grad\,}_w(\boldsymbol{v}, \boldsymbol{\mu})\in \Sigma_k^{-1}(\mathbb{T})\) be defined elementwise by \[(\operatorname{dev}{\rm grad\,}_w(\boldsymbol{v}, \boldsymbol{\mu}), \boldsymbol{\tau})_T = -(\boldsymbol{v},\operatorname{div}\boldsymbol{\tau})_T + (\boldsymbol{n}\cdot\boldsymbol{v},\, \boldsymbol{n}^{\intercal}\boldsymbol{\tau}\boldsymbol{n})_{\partial T} + (\boldsymbol{\mu},\, \Pi_F\boldsymbol{\tau}\boldsymbol{n})_{\partial T},\] for all \(\boldsymbol{\tau}\in \mathbb{P}_{k}(T;\mathbb{T})\) and \(T\in\mathcal{T}_h\).

3.2 Weak div stability↩︎

We first establish the following inf-sup condition.

We use \((\cdot,\cdot)_{0,h}\) as a Riesz map to identify \(E(\mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}} \times \mathring{\Lambda}_k)\) with its dual. The bilinear form 22 thus induces a mapping \[\operatorname{div}_w: \Sigma_k^{-1}(\mathbb{T})\to \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}} \times \mathring{\Lambda}_k.\] The inf-sup condition ?? implies \(\operatorname{div}_w\) is surjective. Consequently the adjoint \(\operatorname{dev}{\rm grad\,}_w\) is injective and \[\| \operatorname{dev}{\rm grad\,}_w (\cdot, \cdot)\| \text{ defines a norm on } \mathring{\mathbb{V}}_{k,k-1}^{\operatorname{div}}\times \mathring{\Lambda}_k.\] Notice that the velocity space cannot be enlarged to \(\mathring{\mathbb{V}}_{k,k}^{\operatorname{div}}\). But later on we will use the coercivity on the divergence free subspaces.

3.3 Mixed finite element methods↩︎

With the operator \(\operatorname{dev}{\rm grad\,}_w\), we formulate the following mixed finite element method for the Stokes equation 1 : Find \(\boldsymbol{u}_h\in \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\), \(\boldsymbol{\lambda}_h\in\mathring{\Lambda}_k\), and \(p_h\in\mathbb{P}_{\ell}(\mathcal{T}_h)/\mathbb{R}\), for integers \(k \ge 0\) and \(\ell\in\{k,\,k-1\}\), such that \[\tag{23} \begin{align} \tag{24} (\operatorname{dev}{\rm grad\,}_w(\boldsymbol{u}_h,\boldsymbol{\lambda}_h),\, \operatorname{dev}{\rm grad\,}_w(\boldsymbol{v}_h,\boldsymbol{\mu}_h)) + (\operatorname{div}\boldsymbol{v}_h,\, p_h) &= (\boldsymbol{f},\, \boldsymbol{v}_h), \\ \tag{25} (\operatorname{div}\boldsymbol{u}_h,\, q_h) &= 0, \end{align}\] for all \(\boldsymbol{v}_h\in \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\), \(\boldsymbol{\mu}_h\in\mathring{\Lambda}_k\), and \(q_h\in\mathbb{P}_{\ell}(\mathcal{T}_h)/\mathbb{R}\).

We circumvent the inf-sup condition, as the error analysis will be derived through the equivalence to 19 without relying on the inf-sup constant.

3.4 Error analysis↩︎

Hereafter, let \(\boldsymbol{u} \in H_0^1(\Omega;\mathbb{R}^d)\) denote the solution of the Stokes equation 1 , and set \(\boldsymbol{\sigma} = {\rm grad\,}\boldsymbol{u} = \operatorname{dev}{\rm grad\,}\boldsymbol{u}\). Then the first equation in 1 becomes \[\label{stokes1} -\operatorname{div}\boldsymbol{\sigma}- \nabla p=\boldsymbol{f} \quad \text{ in } \Omega.\tag{27}\] Moreover, let \(\boldsymbol{\lambda} \in L^2(\mathcal{F}_h; \mathbb{R}^{d-1})\) be defined facewise by \[\boldsymbol{\lambda}|_F = \Pi_F \boldsymbol{u}, \quad \forall\,F \in \mathcal{F}_h.\] This notation will be used throughout the subsequent analysis.

The estimate ?? demonstrates that the mixed finite element method 23 is pressure-robust.

3.5 \(L^2\) error estimate↩︎

The superconvergence of \(\|I_{k,k}^{\operatorname{div}}\boldsymbol{u}-\boldsymbol{u}_h\|\) can be obtained through the duality argument. Introduce the dual problem: Find \(\tilde{\boldsymbol{u}}\in H_0^1(\Omega;\mathbb{R}^d)\) and \(\tilde{p}\in L^2(\Omega)/\mathbb{R}\) satisfy \[\label{dualstokes} \begin{cases} - \Delta \tilde{\boldsymbol{u}} + \nabla \tilde{p}= I_{k,k}^{\operatorname{div}}\boldsymbol{u}-\boldsymbol{u}_h& \text{ in } \Omega, \\ \qquad\;\;\operatorname{div}\tilde{\boldsymbol{u}} =0 & \text{ in } \Omega. \end{cases}\tag{28}\] We assume that the dual problem 28 has the \(H^2\)-regularity \[\label{eq:regularity} \|\tilde{\boldsymbol{u}}\|_2 \lesssim \|I_{k,k}^{\operatorname{div}}\boldsymbol{u}-\boldsymbol{u}_h\|.\tag{29}\] We refer to [45] and [46] for the \(H^2\)-regularity 29 in convex domains.

3.6 Postprocessing↩︎

We can construct a superconvergent approximation of \(\boldsymbol{u}\) using the estimates ?? and ?? . For each \(T \in \mathcal{T}_h\), find \(\boldsymbol{u}_h^{*}|_T\in \mathbb{P}_{k+1}(T;\mathbb{R}^d)\) and \(p_h^{*}\in\mathbb{P}_k(T)\cap L_0^2(T)\) such that \[\tag{31} \begin{align} (\boldsymbol{u}_h^{*}\cdot\boldsymbol{n},q)_F &= (\boldsymbol{u}_h\cdot \boldsymbol{n},q)_F \quad\;\;\; \forall\,q\in\mathbb{P}_0(F),\;F\in\mathcal{F}(T), \tag{32}\\ ({\rm grad\,}\boldsymbol{u}_h^{*},{\rm grad\,}\boldsymbol{v})_T +(\operatorname{div}\boldsymbol{v},p_h^{*})_T &= (\boldsymbol{\sigma}_h,{\rm grad\,}\boldsymbol{v})_T \quad \forall\,\boldsymbol{v}\in\widetilde{\mathbb{P}}_{k+1}(T;\mathbb{R}^d), \tag{33}\\ (\operatorname{div}\boldsymbol{u}_h^{*},q)_T &= 0 \qquad\qquad\qquad\; \forall\,q\in{\mathbb{P}}_{k}(T)\cap L_0^2(T), \tag{34} \end{align}\] where \(\widetilde{\mathbb{P}}_{k+1}(T;\mathbb{R}^d) := \{\boldsymbol{v}_h\in\mathbb{P}_{k+1}(T;\mathbb{R}^d): Q_{0,F}(\boldsymbol{v}_h\cdot\boldsymbol{n})=0\;\;\forall\,F\subset\partial T\}\).

The local problems 31 are well-posed. From 32 we get \[Q_{0,T}(\operatorname{div}\boldsymbol{u}_h^{*}) =Q_{0,T}(\operatorname{div}\boldsymbol{u}_h)=0;\] combined with 34 , this yields \(\operatorname{div}\boldsymbol{u}_h^{*}=0\) on each \(T\in\mathcal{T}_h\). Thus \(\boldsymbol{u}_h^{*}\in \mathbb{P}_{k+1}(\mathcal{T}_h; \mathbb{R}^d)\) is pointwise divergence-free elementwise.

4 Equivalent Discrete Methods↩︎

In this section, we explore several discrete methods that are equivalent to the mixed finite element scheme 23 .

4.1 Hybridization of the velocity↩︎

Recall the \(H^1\)-nonconforming virtual element in [50], [51]. For integers \(k\geq0\) and \(\ell=k,k-1\), the shape function space is defined as \[V_{k+1,\ell+2}^{\rm VE}(T)=\{v\in H^1(T):\Delta v\in\mathbb{P}_{\ell}(T), \partial_n v|_F\in\mathbb{P}_{k}(F) \textrm{ for } F\in\mathcal{F}(T)\}.\] The DoFs for \(V_{k+1,\ell+2}^{\rm VE}(T)\) are given by \[\tag{36} \begin{align} \tag{37} (v,q)_F, & \quad q\in \mathbb{P}_{k}(F), \, F\in\mathcal{F}(T),\\ \tag{38} (v,q)_T, & \quad q\in \mathbb{P}_{\ell}(T). \end{align}\] The global virtual element space \(V_{k+1,\ell+2}^{\rm VE} = V_{k+1,\ell+2}^{\rm VE}(\mathcal{T}_h)\) by \[\begin{align} V_{k+1,\ell+2}^{\rm VE}(\mathcal{T}_h) &= \{v\in L^2(\Omega):v|_T\in V_{k+1,\ell+2}^{\rm VE}(T) \textrm{ for each }T\in\mathcal{T}_h; \textrm{ DoF \eqref{ve-dof1} is } \\ &\qquad\qquad\qquad\qquad\qquad\quad\quad\;\;\textrm{ single-valued across each face in \mathring{\mathcal{F}}_h}\}. \end{align}\] The virtual element space \(V_{k+1,\ell+2}^{\rm VE}(\mathcal{T}_h)\) satisfies the weak continuity condition \[([\![v]\!], q)_F = 0 \quad \forall\, v\in V_{k+1,\ell+2}^{\rm VE}(\mathcal{T}_h), \; q\in\mathbb{P}_{k}(F), \; F\in\mathring{\mathcal{F}}_h.\] When \((k,\ell) = (0,-1)\), the space \(V_{1,1}^{\rm VE}\) reduces to the Crouzeix–Raviart (CR) element [3]. For \((k,\ell) = (0,0)\), the space \(V_{1,2}^{\rm VE}\) corresponds to the enriched Crouzeix–Raviart element [52], [53].

Find \(\boldsymbol{u}_h\in \prod_{T\in \mathcal{T}_h}\mathbb{V}_{k,\ell}^{\operatorname{div}}(T)\), \(\boldsymbol{\lambda}_h\in\mathbb{P}_{k}(\mathring{\mathcal{F}}_h;\mathbb{R}^{d-1})\) and \(p_h\in V_{k+1,\ell+2}^{\rm VE}(\mathcal{T}_h)/\mathbb{R}\) such that \[\tag{39} \begin{align} \tag{40} (\operatorname{dev}{\rm grad\,}_w(\boldsymbol{u}_h,\boldsymbol{\lambda}_h), \operatorname{dev}{\rm grad\,}_w(\boldsymbol{v}_h,\boldsymbol{\mu}_h)) - (\boldsymbol{v}_h, \nabla_hp_h) &= (\boldsymbol{f}, \boldsymbol{v}_h), \\ \tag{41} (\boldsymbol{u}_h, \nabla_hq_h)&=0 \end{align}\] for all \(\boldsymbol{v}_h\in \prod_{T\in \mathcal{T}_h}\mathbb{V}_{k,\ell}^{\operatorname{div}}(T)\), \(\boldsymbol{\mu}_h\in\mathbb{P}_{k}(\mathring{\mathcal{F}}_h;\mathbb{R}^{d-1})\) and \(q_h\in V_{k+1,\ell+2}^{\rm VE}(\mathcal{T}_h)/\mathbb{R}\).

The normal continuity for the velocity space is relaxed in 39 which is particularly useful when \(\ell = -1\). As a local basis for the space \(\mathring{\mathbb{V}}_{0,-1}^{\operatorname{div}}\) is not readily available, we can implement the equivalent mixed method 39 instead of the mixed method 23 .

Even for \((k, \ell)=(0, -1)\), our method has the following estimates: \[\begin{align} &\|\boldsymbol{\sigma} - \boldsymbol{\sigma}_h\|_{0,h} + \|{\rm grad\,}_h(\boldsymbol{u}-\boldsymbol{u}_h^{*})\| + \|\operatorname{dev}{\rm grad\,}_w(I_{0,0}^{\operatorname{div}}\boldsymbol{u} - \boldsymbol{u}_h,Q_{0,\mathcal{F}_h}\boldsymbol{\lambda}-\boldsymbol{\lambda}_h)\| \lesssim h |\boldsymbol{u}|_{2},\\ &\|\boldsymbol{u} - \boldsymbol{u}_h\|_0+h\|{\rm grad\,}_h(\boldsymbol{u} - \boldsymbol{u}_h)\| \lesssim h( |\boldsymbol{u}|_{2} + |\boldsymbol{u}|_{1}),\\ &\|\boldsymbol{u}-\boldsymbol{u}_h^{*}\|+\|I_{0,0}^{\operatorname{div}}\boldsymbol{u} - \boldsymbol{u}_h\| \lesssim h^{2}(|\boldsymbol{u}|_{2} + \|\operatorname{skw}{\rm grad\,}\boldsymbol{f}\|). \end{align}\]

4.2 Stabilization-free virtual element methods↩︎

For \(k\geq0\), define the space of shape functions for an \(H^1\)-nonconforming virtual element as \[\begin{align} \mathbb{V}_{k+1}^{\rm VE}(T):=\{\boldsymbol{v}\in H^1(T;\mathbb{R}^d): &\operatorname{div}\boldsymbol{v}\in\mathbb{P}_{k}(T), \textrm{ there exists some } s\in L^2(T) \\ & \textrm{ such that } \Delta\boldsymbol{v}+\nabla s\in \big(\mathbb{P}_{k-1}(T;\mathbb{R}^d)\cap\ker(\cdot\boldsymbol{x})\big), \\ &\textrm{ and } (\partial_n\boldsymbol{v}+s\boldsymbol{n})|_F\in\mathbb{P}_{k}(F;\mathbb{R}^d) \;\forall~F\in\mathcal{F}(T)\}. \end{align}\] Following the argument in [26], the local virtual element space \(\mathbb{V}_{k+1}^{\rm VE}(T)\) is uniquely determined by the following DoFs \[\tag{43} \begin{align} (\boldsymbol{v}, \boldsymbol{q})_F, & \quad \boldsymbol{q}\in\mathbb{P}_{k}(F; \mathbb{R}^d), F\in\mathcal{F}(T), \tag{44}\\ (\boldsymbol{v}, \boldsymbol{q})_T, & \quad \boldsymbol{q}\in\mathbb{P}_{k-1}(T;\mathbb{R}^d). \tag{45} \end{align}\] Clearly we have \(\mathbb{P}_{k+1}(T;\mathbb{R}^d)\subseteq\mathbb{V}_{k+1}^{\rm VE}(T)\), and \(\mathbb{V}_{1}^{\rm VE}(T)=\mathbb{P}_1(T;\mathbb{R}^d)\). Then define the global virtual element space by \[\begin{align} \mathring{\mathbb{V}}_{h}^{\rm VE} &:= \{\boldsymbol{v}\in L^2(\Omega;\mathbb{R}^d):\boldsymbol{v}|_T\in\mathbb{V}_{k+1}^{\rm VE}(T) \textrm{ for each }T\in\mathcal{T}_h; \textrm{ DoF \eqref{vev-dof1} is } \\ &\qquad\qquad\qquad\quad\quad\;\;\textrm{ single-valued across each face in \mathring{\mathcal{F}}_h, and vanishes on \partial\Omega}\}. \end{align}\] A mixed virtual element method for the Stokes equation 1 is to find \(\boldsymbol{u}_h\in \mathring{\mathbb{V}}_{h}^{\rm VE}\) and \(p_h\in\mathbb{P}_{k}(\mathcal{T}_h)/\mathbb{R}\) with integer \(k\geq 0\), such that \[\tag{46} \begin{align} \tag{47} (Q_{k}(\operatorname{dev}{\rm grad\,}_h\boldsymbol{u}_h), Q_{k}(\operatorname{dev}{\rm grad\,}_h\boldsymbol{v}_h)) + (\operatorname{div}_h\boldsymbol{v}_h, p_h) &= (\boldsymbol{f}, I_{k,k}^{\operatorname{div}}\boldsymbol{v}_h), \\ \tag{48} (\operatorname{div}_h\boldsymbol{u}_h, q_h)&=0 \end{align}\] for all \(\boldsymbol{v}_h\in \mathring{\mathbb{V}}_{h}^{\rm VE}\) and \(q_h\in\mathbb{P}_{k}(\mathcal{T}_h)/\mathbb{R}\). When \(k=0\), the mixed method 46 is exactly the modified Crouzeix-Raviart element method in [55].

As a consequence of ?? , \(\|Q_{k}(\operatorname{dev}{\rm grad\,}_h\boldsymbol{v}_h)\|\) defines a norm on the space \(\mathring{\mathbb{V}}_{h}^{\rm VE}\cap\ker(\operatorname{div}_h)\).

4.3 Discrete methods based on the pseudostress-velocity-pressure formulation↩︎

By introducing \(\boldsymbol{\sigma}_h:= \operatorname{dev}{\rm grad\,}_w\boldsymbol{u}_h\in\Sigma_{k}^{\rm tn}\), we can rewrite the mixed finite element method 23 as the MCS method [38]: find \(\boldsymbol{\sigma}_h\in \Sigma_{k}^{\rm tn}\), \(\boldsymbol{u}_h\in\mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}\) and \(p_h\in\mathbb{P}_{\ell}(\mathcal{T}_h)/\mathbb{R}\) such that \[\begin{align} (\boldsymbol{\sigma}_h,\boldsymbol{\tau}_h)+b_h(\boldsymbol{\tau}_h,q_h;\boldsymbol{u}_h)&=0 & & \forall~\boldsymbol{\tau}_h \in \Sigma_{k}^{\rm tn}, q_h\in\mathbb{P}_{\ell}(\mathcal{T}_h)/\mathbb{R}, \\ b_h(\boldsymbol{\sigma}_h,p_h;\boldsymbol{v}_h) &=-(\boldsymbol{f}, \boldsymbol{v}_h) & &\forall~\boldsymbol{v}_h \in \mathring{\mathbb{V}}_{k,\ell}^{\operatorname{div}}, \end{align}\] where the bilinear form \[\begin{align} b_h(\boldsymbol{\tau}_h,q_h;\boldsymbol{v}_h)&:=\sum_{T\in \mathcal{T}_{h}}(\operatorname{div}\boldsymbol{\tau}_h, \boldsymbol{v}_h)_T-\sum_{T\in \mathcal{T}_{h}}(\boldsymbol{n}^\intercal\boldsymbol{\tau}_h\boldsymbol{n}, \boldsymbol{n}\cdot\boldsymbol{v}_h)_{\partial T} - (\operatorname{div}\boldsymbol{v}_h, q_h) \\ &\;=-\sum_{T\in \mathcal{T}_{h}}(\boldsymbol{\tau}_h, {\rm grad\,}\boldsymbol{v}_h)_T+\sum_{T\in \mathcal{T}_{h}}(\Pi_F\boldsymbol{\tau}_h\boldsymbol{n}, \Pi_F\boldsymbol{v}_h)_{\partial T} - (\operatorname{div}\boldsymbol{v}_h, q_h). \end{align}\]

The subspace \[\left\{ \boldsymbol{\tau}_h \in \Sigma_{k}^{\rm tn} : \Pi_F \boldsymbol{\tau}_h \boldsymbol{n} \in \mathbb{P}_{k-1}(F; \mathbb{R}^{d-1}) \;\text{for all } F \in \mathcal{F}_h \right\}\] was employed in the MCS method [38] to discretize the traceless stress tensor \(\boldsymbol{\sigma}\). However, the constraint \(\Pi_F \boldsymbol{\tau}_h \boldsymbol{n} \in \mathbb{P}_{k-1}(F; \mathbb{R}^{d-1})\) prevents the discrete scheme in [38] from achieving superconvergence properties such as the estimates \(\|\boldsymbol{\sigma} - \boldsymbol{\sigma}_h\|\) and \(\|Q_{\ell} p - p_h\|\) in ?? , and \(\|I_{k,k}^{\operatorname{div}}\boldsymbol{u} - \boldsymbol{u}_h\|\) in ?? .

Superconvergent results were later obtained in [37] by enriching the stress space, but their approach requires \(k = \ell \ge 1\), whereas our formulation also covers low-order cases: \((k,\ell) = (0,-1)\), \((0,0)\), and \((1,0)\). Our approach also leads to a conceptually simpler and technically more transparent error analysis.

5 Numerical examples↩︎

In this section, we present numerical examples to validate the error estimates of the mixed finite element method 23 and the post-processing procedure 31 . The test problems are adapted from [37], which primarily considered higher-order cases with \(k=\ell\ge1\). Here, we focus on low-order settings with \(\ell=0\) and \(k=0,1\).

The computational domain is \(\Omega=(0,1)^d\) for \(d=2,3\), uniformly partitioned into simplicial meshes. All numerical experiments were developed based on the MATLAB package \(i\)FEM [56].

The numerical errors \(\|\boldsymbol{u}-\boldsymbol{u}_h\|\), \(\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_h\|\), and \(\|p - p_h\|\) are reported in Table 1 for Example 5.1 and Table 3 for Example 5.2, while the post-processed errors \(\|\boldsymbol{u} - \boldsymbol{u}_h^{*}\|\) and \(\|{\rm grad\,}_h(\boldsymbol{u}-\boldsymbol{u}_h^{*})\|\) are given in Table 2 for Example 5.1 and Table 4 for Example 5.2. The results in these tables demonstrate that \[\|\boldsymbol{u}-\boldsymbol{u}_h\|\eqsim\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_h\|=\mathcal{O}(h^{k+1}),\quad \|p - p_h\|=\mathcal{O}(h),\] and for the post-processed solution \(\boldsymbol{u}_h^{*}\), \[\|\boldsymbol{u} - \boldsymbol{u}_h^{*}\|=\mathcal{O}(h^{k+2}),\quad \|{\rm grad\,}_h(\boldsymbol{u}-\boldsymbol{u}_h^{*})\|=\mathcal{O}(h^{k+1}),\] which are consistent with the theoretical estimates ?? and ?? –?? .

References↩︎

---

1. The first author was supported by NSF DMS-2309777 and DMS-2309785. The second author was supported by the National Natural Science Foundation of China Project 12171300.↩︎