1 Introduction↩︎

We shall develop a symmetric-tensor distributional mixed method for the fourth-order elliptic singular perturbation problem: \[\label{eq:fourthorderequation} \begin{cases} \varepsilon^2\Delta^2 u-\Delta u=f & \text{in }\Omega,\\ u=\partial_n u=0 & \text{on }\partial\Omega, \end{cases}\tag{1}\] where \(\Omega\subset\mathbb{R}^d\) (\(d\ge2\)) is a bounded polytope and \(\varepsilon>0\) is the perturbation parameter. As \(\varepsilon\to0\), the operator formally degenerates to the Poisson operator, whereas the normal-derivative condition \(\partial_n u=0\) has no counterpart in the reduced problem. This incompatibility may generate boundary layers near \(\partial\Omega\) and is the main difficulty in deriving optimal error estimates uniform in \(\varepsilon\).

Several scalar primal discretizations have been proposed for 1 . Conforming \(H^2\) methods discretize the fourth- and second-order terms directly, but require \(C^1\) finite elements and are difficult to implement in general [1]. Many \(C^1\)-free primal methods are robust in suitable norms, including nonconforming \(H^2\) methods [2], \(C^0\) interior-penalty methods [3], Morley–Wang–Xu type methods [4]–[7], and virtual element methods [8]–[10]. However, parameter-uniform estimates for such methods are often suboptimal in the boundary-layer regime.

Optimal or high-order parameter-uniform results within scalar primal formulations have been obtained by restoring compatibility with the reduced Poisson problem or by using layer-adapted meshes, including Nitsche-type methods [11], [12], hybrid high-order methods [13], Nitsche-modified Morley–Wang–Xu methods [7], and layer-adapted \(C^0\) interior-penalty methods [14]. These approaches typically involve penalty or stabilization terms. The interpolation-based method of [15] also gives optimal parameter-uniform estimates by introducing an interpolation from an \(H^2\) finite element space to an \(H^1\) finite element space in both the load term and the Laplacian bilinear form.

Mixed formulations provide another approach by exposing the fourth-order structure through a moment tensor. For 1 , the Hellan–Herrmann–Johnson (HHJ) method [16]–[18] was used in two dimensions in [19], but without a parameter-robust error analysis. Robust and optimal mixed methods based on row-wise \(H(\operatorname{div})\)-conforming tensor discretizations were developed in [20]; in that setting the tensor variable is not symmetric. In the present problem, however, the moment variable is naturally symmetric because it is proportional to the Hessian. In Kirchhoff–Love plate models [21], this variable may be interpreted as a scaled bending-moment tensor. The symmetry of such stress or moment tensors is consistent with the symmetry of the Cauchy stress tensor in continuum mechanics, which follows from angular-momentum balance; see, e.g., [22].

The central idea of this paper is to introduce \(\boldsymbol{\sigma}=\varepsilon^2\nabla^2u\) as a symmetric tensor field. Then 1 can be written as \[\label{eq:intro-mixed-system} \begin{cases} \varepsilon^{-2}\boldsymbol{\sigma}=\nabla^2u & \text{in }\Omega,\\ \operatorname{div}\operatorname{div}\boldsymbol{\sigma}-\Delta u=f & \text{in }\Omega,\\ u=\partial_nu=0 & \text{on }\partial\Omega. \end{cases}\tag{2}\] The resulting distributional mixed formulation 4 seeks \(\boldsymbol{\sigma}\in H^{-1}(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})\) and \(u\in H_0^1(\Omega)\), where \[H^{-1}(\operatorname{div}\operatorname{div},\Omega;\mathbb{S}) :=\{\boldsymbol{\tau}\in L^2(\Omega;\mathbb{S}): \operatorname{div}\operatorname{div}\boldsymbol{\tau}\in H^{-1}(\Omega)\},\] and \(\mathbb{S}\) denotes the space of symmetric tensors. The well-posedness of this formulation is tied to the exactness of the terminal part of the distributional \(\operatorname{div}\operatorname{div}\) complex \[\label{intro:shortdistritdivdiv} H^{-1}(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})\xrightarrow{\operatorname{div}\operatorname{div}}H^{-1}(\Omega)\to0.\tag{3}\]

At the discrete level, we use the distributional finite element \(\operatorname{div}\operatorname{div}\) complexes in [23]. The moment variable is approximated by normal–normal continuous symmetric tensor elements [23]–[25]. The scalar approximation uses an \(H^1\)-nonconforming virtual element space and a polynomial multiplier on interior codimension-two subsimplices. The resulting distributional mixed method 33 possesses optimal parameter-uniform error estimates ?? –?? , independent of the presence of boundary layers. We also give a tangential–normal reconstruction of the normal–normal continuous symmetric tensor element in which all boundary degrees of freedom are located on codimension-one faces. The normal–normal continuous symmetric tensor elements have applications in the tangential-displacement and normal-normal-stress (TDNNS) method for elasticity [26] and Reissner–Mindlin plates [27], as well as low-order mixed methods for elasticity [28], [29].

In two dimensions, we establish a precise connection between the distributional mixed method 33 and the classical HHJ method. Specifically, the \(H^1\)-nonconforming virtual element space together with the codimension-two multiplier can be identified with a Lagrange finite element space; see Figure 1. Under this identification, the discrete distributional bilinear form agrees with the HHJ bilinear form. The proposed method can therefore be viewed as an arbitrary-dimensional extension of the two-dimensional HHJ method. The standard HHJ coupling with a Lagrange scalar space is essentially two-dimensional and does not directly yield a discretization in arbitrary dimension. We overcome this obstruction by coupling the moment space with the \(H^1\)-nonconforming virtual element space and the polynomial multiplier.

As a further consequence, relaxing the normal–normal continuity gives a hybridized formulation. After local elimination of the tensor variable, the hybridized method is equivalent to stabilization-free weak Galerkin and \(H^2\)-nonconforming virtual element methods.

The rest of the paper is organized as follows. Section 2 presents the continuous distributional mixed formulation. Section 3 constructs the discrete spaces and establishes the norm equivalence. Section 4 introduces the mixed method and proves the robust error estimates. Section 5 derives the hybridized formulation and its equivalent weak Galerkin and \(H^2\)-nonconforming virtual element formulations. Section 6 reports numerical experiments.

2 Symmetric-Tensor Distributional Mixed Formulation↩︎

In this section, we present a symmetric-tensor distributional mixed formulation of 1 , establish its well-posedness and equivalence with the primal weak formulation, and state the regularity assumptions used in the error analysis.

2.1 Notation↩︎

Throughout the paper, \(\Omega\subset\mathbb{R}^d\) (\(d\ge2\)) denotes a bounded polytope. For a bounded domain \(D\subset\mathbb{R}^d\) and \(m\ge0\), we denote by \(H^m(D)\) the standard Sobolev space on \(D\) with norm \(\|\cdot\|_{m,D}\) and seminorm \(|\cdot|_{m,D}\). We write \(H_0^1(D)\) and \(H_0^2(D)\) for the closures of \(C_0^\infty(D)\) in \(H^1(D)\) and \(H^2(D)\), respectively, and use \((\cdot,\cdot)_D\) to denote the \(L^2(D)\) inner product. When \(D=\Omega\), the subscript is omitted. The duality pairing between a Banach space \(V\) and its dual \(V'\) is denoted by \(\langle\cdot,\cdot\rangle_{V'\times V}\), or simply by \(\langle\cdot,\cdot\rangle\) if no confusion arises.

For an integer \(k\ge 0\), let \(\mathbb{P}_k(D)\) denote the space of polynomials on \(D\) of total degree at most \(k\). If \(D\) is a vertex, we identify \(\mathbb{P}_k(D)\) with \(\mathbb{R}\). We adopt the convention \(\mathbb{P}_k(D)=\{0\}\) for \(k<0\). Let \[\mathbb{S}:=\{\boldsymbol{\tau}\in\mathbb{R}^{d\times d}: \boldsymbol{\tau}^{\intercal}=\boldsymbol{\tau}\}\] be the space of symmetric tensors. For \(\boldsymbol{\tau}\in\mathbb{R}^{d\times d}\), set \[\operatorname{sym}\boldsymbol{\tau} :=\frac{1}{2}(\boldsymbol{\tau}+\boldsymbol{\tau}^{\intercal}).\] For a finite-dimensional space \(\mathbb{X}\), we use \(L^2(D;\mathbb{X})\) and \(H^m(D;\mathbb{X})\) for the corresponding \(\mathbb{X}\)-valued Lebesgue and Sobolev spaces, and set \(\mathbb{P}_k(D;\mathbb{X}):=\mathbb{P}_k(D)\otimes\mathbb{X}\). We denote by \(Q_{k,D}\) the \(L^2(D;\mathbb{X})\)-orthogonal projection onto \(\mathbb{P}_k(D;\mathbb{X})\), with the convention \(Q_{-1,D}=0\).

For a scalar function \(v\), we use \(\nabla v\) and \(\nabla^2 v\) for its gradient and Hessian, respectively. For a vector-valued function \(\boldsymbol{v}\), we write \(\nabla\boldsymbol{v}:=\nabla\otimes\boldsymbol{v}\). For a tensor-valued function \(\boldsymbol{\tau}=(\tau_{ij})_{i,j=1}^d\), the divergence is taken row-wise, namely, \((\operatorname{div}\boldsymbol{\tau})_i:=\sum_{j=1}^d\partial_j\tau_{ij}\) for \(i=1,\dots,d\), and \(\operatorname{div}\operatorname{div}\boldsymbol{\tau}\) is understood in the distributional sense whenever needed.

Let \(\mathcal{T}_h\) be a conforming shape-regular simplicial triangulation of \(\Omega\). For each \(T\in\mathcal{T}_h\), we write \(h_T:=\operatorname{diam}(T)\) and set \(h:=\max_{T\in\mathcal{T}_h}h_T\). We denote by \(\Delta_j(T)\) the set of all \(j\)-dimensional subsimplices of \(T\). The sets of all faces and all interior faces are denoted by \(\mathcal{F}_h\) and \(\mathring{\mathcal{F}}_h\), respectively, while \(\mathcal{E}_h\) and \(\mathring{\mathcal{E}}_h\) denote the sets of all \((d-2)\)-dimensional subsimplices and the interior ones. For \(e\in\mathcal{E}_h\), we define the patch of elements around \(e\) by \(\omega_e:=\bigcup\{T\in\mathcal{T}_h:\;e\subset\overline{T}\}\).

For each simplex \(T\in\mathcal{T}_h\) with vertices \(\{\texttt{v}_i\}_{i=0}^d\), we denote by \(F_i\) the face opposite to \(\texttt{v}_i\), by \(\{\lambda_i\}_{i=0}^d\) the corresponding barycentric coordinates, and by \(\boldsymbol{n}_{\partial T}\) the unit outward normal on \(\partial T\). On each local face \(F_i\subset\partial T\), we write \(\boldsymbol{n}_i:=\boldsymbol{n}_{\partial T}|_{F_i}\). For a global face \(F\in\mathcal{F}_h\), \(\boldsymbol{n}_F\) denotes a fixed unit normal on \(F\), chosen as the outward unit normal to \(\partial\Omega\) when \(F\subset\partial\Omega\). If \(F\in\mathring{\mathcal{F}}_h\), we denote by \(T^+\) and \(T^-\) the two elements sharing \(F\), with \(\boldsymbol{n}_F\) taken outward on \(T^+\). For \(e\subset \partial F\), let \(\boldsymbol{n}_{F,e}\) denote the outward unit normal to \(e\) within the hyperplane containing \(F\). We also set the edge vector \(\boldsymbol{t}_{i,j}:=\texttt{v}_j-\texttt{v}_i\) for \(i,j=0,1,\dots,d\). For a face \(F\in\mathcal{F}_h\), \(\operatorname{div}_F\) denotes the surface divergence on \(F\). For \(\boldsymbol{\tau}\in\mathbb{S}\), its normal–normal component on a local face \(F\subset\partial T\) is defined by \(\boldsymbol{\tau}_{nn}:=\boldsymbol{n}_{\partial T}^{\intercal} \boldsymbol{\tau}\boldsymbol{n}_{\partial T}\).

For \(\mathcal{G}_h\in\{\mathcal{T}_h,\mathcal{F}_h,\mathring{\mathcal{F}}_h, \mathcal{E}_h,\mathring{\mathcal{E}}_h\}\), we define the broken Sobolev space and the piecewise polynomial space by \[H^s(\mathcal{G}_h):=\prod_{G\in\mathcal{G}_h} H^s(G), \qquad \mathbb{P}_k(\mathcal{G}_h):=\prod_{G\in\mathcal{G}_h}\mathbb{P}_k(G).\] For a finite-dimensional space \(\mathbb{X}\), we set \[H^s(\mathcal{G}_h;\mathbb{X}):=H^s(\mathcal{G}_h)\otimes \mathbb{X}, \qquad \mathbb{P}_k(\mathcal{G}_h;\mathbb{X}):= \mathbb{P}_k(\mathcal{G}_h)\otimes \mathbb{X} .\] We denote by \(Q_{k,\mathcal{G}_h}\) the elementwise \(L^2\)-orthogonal projection onto \(\mathbb{P}_k(\mathcal{G}_h;\mathbb{X})\), and set \(Q_{k,h}:=Q_{k,\mathcal{T}_h}\). Here and throughout, functions defined on interior faces or interior \((d-2)\)-subsimplices are understood to be extended by zero to the corresponding boundary entities whenever needed.

For piecewise smooth scalar-, vector-, or tensor-valued functions, we use \(\nabla_h\), \(\nabla_h^2\), \(\operatorname{div}_h\), and \(\operatorname{div}\operatorname{div}_h\) to denote the elementwise gradient, Hessian, divergence, and double divergence, respectively. We also use the broken norm and seminorm \[\|v\|_{m,h}^2:=\sum_{T\in\mathcal{T}_h}\|v\|_{m,T}^2, \qquad |v|_{m,h}^2:=\sum_{T\in\mathcal{T}_h}|v|_{m,T}^2,\] with the same notation understood componentwise for vector- and tensor-valued functions.

For a face \(F\in\mathcal{F}_h\) and a piecewise smooth scalar-, vector-, or tensor-valued function \(\phi\), we define the jump by \[[\![\phi]\!]|_F := \begin{cases} \phi^+\boldsymbol{n}_F\cdot\boldsymbol{n}_{\partial T^+} + \phi^-\boldsymbol{n}_F\cdot\boldsymbol{n}_{\partial T^-}, & \text{if } F=\partial T^+\cap\partial T^-\in\mathring{\mathcal{F}}_h,\\ \phi|_F, & \text{if } F\subset\partial\Omega, \end{cases}\] where \(\phi^\pm\) denote the traces from \(T^\pm\).

Throughout this paper, \(C\) denotes a generic positive constant independent of the mesh size \(h\) and the perturbation parameter \(\varepsilon\). We write \(a\lesssim b\) if \(a\le Cb\), and \(a\eqsim b\) if \(a\lesssim b\) and \(b\lesssim a\).

We shall use the following Green identity for the \(\operatorname{div}\operatorname{div}\) operator; see [30] and [31].

2.2 A distributional mixed formulation↩︎

We formulate 2 in a symmetric-tensor distributional setting. The symmetric tensor variable \(\boldsymbol{\sigma}\) is kept as a primary unknown, with only the natural regularity dictated by the second equation in 2 . Since \(\operatorname{div}\operatorname{div}\boldsymbol{\sigma}\) is tested against functions in \(H_0^1(\Omega)\), it is natural to introduce \[H^{-1}(\operatorname{div}\operatorname{div},\Omega;\mathbb{S}) := \left\{ \boldsymbol{\tau}\in L^2(\Omega;\mathbb{S}): \operatorname{div}\operatorname{div}\boldsymbol{\tau}\in H^{-1}(\Omega) \right\},\] equipped with the norm \[\|\boldsymbol{\tau}\|_{H^{-1}(\operatorname{div}\operatorname{div})}^2 := \|\boldsymbol{\tau}\|_0^2+\|\operatorname{div}\operatorname{div}\boldsymbol{\tau}\|_{-1}^2,\; \|\operatorname{div}\operatorname{div}\boldsymbol{\tau}\|_{-1} := \sup_{v\in H_0^1(\Omega),\,v\neq0} \frac{\langle \operatorname{div}\operatorname{div}\boldsymbol{\tau},v\rangle}{|v|_1}.\] For the singularly perturbed problem, we also use the parameter-dependent norm \[\|\boldsymbol{\tau}\|_{\varepsilon^{-1}L^2\cap H^{-1}(\operatorname{div}\operatorname{div})}^2 := \varepsilon^{-2}\|\boldsymbol{\tau}\|_0^2 +\|\operatorname{div}\operatorname{div}\boldsymbol{\tau}\|_{-1}^2.\] The distributional mixed formulation is to find \((\boldsymbol{\sigma},u)\in H^{-1}(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})\times H_0^1(\Omega)\) such that \[\tag{4} \begin{align} a(\boldsymbol{\sigma}, \boldsymbol{\tau}) + b(\boldsymbol{\tau}, u) &= 0, && \forall\, \boldsymbol{\tau} \in H^{-1}(\operatorname{div}\operatorname{div},\Omega;\mathbb{S}), \tag{5}\\ b(\boldsymbol{\sigma}, v) - c(u, v) &= -(f, v), && \forall\, v \in H_0^1(\Omega), \tag{6} \end{align}\] where \[a(\boldsymbol{\sigma}, \boldsymbol{\tau}) := \varepsilon^{-2}(\boldsymbol{\sigma}, \boldsymbol{\tau}), \qquad b(\boldsymbol{\tau}, v) := -\langle \operatorname{div}\operatorname{div}\boldsymbol{\tau}, v \rangle, \qquad c(u, v) := (\nabla u, \nabla v).\] The well-posedness of the distributional mixed formulation 4 is governed by the exactness of the terminal part of the distributional \(\operatorname{div}\operatorname{div}\) complex, \[H^{-1}(\operatorname{div}\operatorname{div},\Omega;\mathbb{S})\xrightarrow{\operatorname{div}\operatorname{div}}H^{-1}(\Omega)\to0.\]

For comparison, we recall the standard primal weak formulation of 1 : find \(u\in H_0^2(\Omega)\) such that \[\label{eq:primal-weak} \varepsilon^2(\nabla^2u,\nabla^2v)+(\nabla u,\nabla v)=(f,v) \qquad \forall\,v\in H_0^2(\Omega).\tag{7}\] Since the bilinear form on the left-hand side of 7 is continuous and coercive on \(H_0^2(\Omega)\), problem 7 admits a unique solution by the Lax–Milgram theorem.

2.3 Regularity assumptions↩︎

To state the regularity assumptions used in the robust error analysis, we let \(\bar u\) denote the solution of the limiting Poisson problem \[\label{eq:poisson} \begin{cases} -\Delta \bar{u} = f, & \text{in } \Omega, \\ \bar{u} = 0, & \text{on } \partial\Omega. \end{cases}\tag{8}\] We regard \(\bar u\) as the regular part of the solution, and \(u-\bar u\) as the singularly perturbed remainder.

We assume the elliptic regularity estimate: for some \(s\ge2\), \[\label{eq:regularity-assumption} \|\bar{u}\|_{s} \lesssim \|f\|_{s-2}.\tag{9}\] If \(\Omega\) is semiconvex, or if the closure of \(\Omega\) has uniformly positive reach, estimate 9 for \(s=2\) holds; see [33]–[37]. In particular, every convex domain is semiconvex.

We further assume that \[\label{eq:perturbation-regularity} |u-\bar{u}|_{1} + \varepsilon |u|_{2} + \varepsilon^{2} |u|_{3} \lesssim \varepsilon^{1/2}\|f\|_{0}.\tag{10}\] Estimate 10 is available on convex domains in two and three dimensions; see, for example, [2], [12]. More recently, [38] established a related arbitrary-dimensional counterpart of 10 .

3 Discrete Spaces, Weak Hessian and Norm Equivalences↩︎

This section introduces the discrete spaces and interpolation operators used in the mixed method, and then proves the weak-Hessian norm equivalence on which the stability analysis rests.

3.1 Normal–normal continuous finite element spaces for symmetric tensors↩︎

Finite elements for symmetric tensors with normal–normal continuity have been constructed in [16]–[18], [23]–[26]. We use here a tangential–normal construction whose boundary degrees of freedom (DoFs) are supported only on faces.

For each simplex \(T\in\mathcal{T}_h\) and each integer \(k\ge 1\), we take \(\mathbb{P}_k(T;\mathbb{S})\) as the local shape function space. The DoFs are defined as follows: \[\tag{11} \begin{align} \bigl(\boldsymbol{\tau}_{nn},q\bigr)_F, &\qquad q\in\mathbb{P}_k(F),\quad F\in\Delta_{d-1}(T), \tag{12}\\ \bigl(\boldsymbol{n}_i^{\intercal}\boldsymbol{\tau}\,\boldsymbol{t}_{0,j},q\bigr)_{F_i}, &\qquad q\in\mathbb{P}_k(F_i),\quad i=1,\ldots,d-2,\quad j=i+1,\ldots,d, \tag{13}\\ (\boldsymbol{\tau},\boldsymbol{q})_T, &\qquad \boldsymbol{q}\in\mathbb{P}_{k-1}(T;\mathbb{S}). \tag{14} \end{align}\] For \(d=2\), the DoFs 13 are absent. The boundary DoFs 12 –13 are associated only with the \((d-1)\)-dimensional faces of \(T\). The interior moments in 14 are taken against the full symmetric tensor-valued polynomial space of degree at most \(k-1\).

To prove unisolvence of the DoFs 11 for \(\mathbb{P}_k(T;\mathbb{S})\), we use the following direct-sum decomposition of \(\mathbb{S}\), adapted to the tangential–normal decomposition of [39], [40].

Let \(\{N_m\}_{m=1}^{\dim\mathbb{S}}\) be the basis of \(\mathbb{S}\) given by ?? , and let \(\{N_m'\}_{m=1}^{\dim\mathbb{S}}\) be the dual basis, defined by \(N_m':N_\ell=\delta_{m\ell}\). Here \(\delta_{m\ell}\) denotes the Kronecker delta.

Accordingly, we define the global finite element space \[\Sigma_{k,h}^{nn} := \{ \boldsymbol{\tau}\in \Sigma^{nn}(\mathcal{T}_h;\mathbb{S}): \boldsymbol{\tau}|_T\in\mathbb{P}_k(T;\mathbb{S})\;\; \forall\,T\in\mathcal{T}_h \},\] where \[\Sigma^{nn}(\mathcal{T}_h;\mathbb{S}) := \{ \boldsymbol{\tau}\in H^1(\mathcal{T}_h;\mathbb{S}) : [\![\boldsymbol{\tau}_{nn}]\!]|_F=0 \;\; \forall\, F\in\mathring{\mathcal{F}}_h \}.\] Thus \(\Sigma_{k,h}^{nn}\) imposes only normal–normal continuity across interelement faces, with no additional edge or vertex continuity constraints.

For \(\boldsymbol{\tau}\in\Sigma^{nn}(\mathcal{T}_h;\mathbb{S})\), let \(I_h^{nn}\boldsymbol{\tau}\in\Sigma_{k,h}^{nn}\) be the canonical interpolant associated with the DoFs 11 . By a standard scaling argument, for any \(1\le s\le k+1\) and any \(T\in\mathcal{T}_h\), we have \[\label{eq:intererrorSigma} \|\boldsymbol{\tau}-I_h^{nn}\boldsymbol{\tau}\|_{0,T} \lesssim h_T^{s}\,|\boldsymbol{\tau}|_{s,T}, \qquad \forall\,\boldsymbol{\tau}\in \Sigma^{nn}(\mathcal{T}_h;\mathbb{S}) \cap H^{s}(\mathcal{T}_h;\mathbb{S}), T\in\mathcal{T}_h.\tag{17}\]

3.2 \(H^1\)-nonconforming virtual element space↩︎

For the scalar variable, we use the \(H^1\)-nonconforming virtual element space introduced in [41], [42]. This space is compatible with the weak continuity structure of \(\Sigma_{k,h}^{nn}\) and the weak Hessian framework developed below. For each integer \(k\ge1\), the local shape function space is defined by \[\label{eq:vem-space-def} V_{k}^{\mathrm{VE}}(T) := \bigl\{ v \in H^1(T): \Delta v \in \mathbb{P}_{k-2}(T), \partial_n v|_F \in \mathbb{P}_{k-1}(F) \; \forall\, F \in \Delta_{d-1}(T) \bigr\}.\tag{18}\] The space \(V_k^{\mathrm{VE}}(T)\) contains \(\mathbb{P}_k(T)\) and coincides with \(\mathbb{P}_1(T)\) when \(k=1\). The associated DoFs are \[\tag{19} \begin{align} \frac{1}{|F|}(v,q)_F, &\qquad q \in \mathbb{P}_{k-1}(F),\; F \in \Delta_{d-1}(T), \tag{20} \\ \frac{1}{|T|}(v,q)_T, &\qquad q \in \mathbb{P}_{k-2}(T). \tag{21} \end{align}\] By a scaling argument, we obtain the \(L^2\)-norm equivalence \[\label{eq:vem-norm-equiv} \|v\|_{0,T}^2\eqsim \|Q_{k-2,T}v\|_{0,T}^2 + \sum_{F\in\Delta_{d-1}(T)}h_F\|Q_{k-1,F}v\|_{0,F}^2 \quad\forall\,v\in V_k^{\mathrm{VE}}(T).\tag{22}\]

The associated global \(H^1\)-nonconforming virtual element space is \[\label{eq:vem-global-space} \mathring V_{k,h}^{\mathrm{VE}} := \{ v\in \mathring{H}_1^{\rm nc}(\mathcal{T}_h) :\; v|_T\in V_k^{\mathrm{VE}}(T)\;\;\forall\,T\in\mathcal{T}_h \},\tag{23}\] where \[\mathring{H}_1^{\rm nc}(\mathcal{T}_h) := \{ v\in H^1(\mathcal{T}_h) : [\![Q_{k-1,F}v]\!]|_F=0 \;\; \forall\,F\in\mathcal{F}_h \}.\] For \(k=1\), the space \(\mathring V_{1,h}^{\mathrm{VE}}\) coincides with the classical Crouzeix–Raviart finite element space [43]. This space satisfies the following weak continuity property: \[\label{eq:vem-weak-continuity} ([\![v]\!], q)_F = 0 \quad \forall\, v \in \mathring V_{k,h}^{\mathrm{VE}},\; q \in \mathbb{P}_{k-1}(F),\; F \in \mathcal{F}_h.\tag{24}\] For each \(T\in\mathcal{T}_h\), the following norm equivalence holds [20]: \[\label{eq:vem-grad-projection-equivalence} \|Q_{k-1,T}\nabla v\|_{0,T} \eqsim \|\nabla v\|_{0,T}, \qquad \forall\, v\in V_k^{\mathrm{VE}}(T).\tag{25}\]

Let \(I_{k,h}^{\mathrm{VE}} : \mathring{H}_1^{\mathrm{nc}}(\mathcal{T}_h) \to \mathring V_{k,h}^{\mathrm{VE}}\) denote the global interpolation operator associated with the DoFs 19 . When no ambiguity arises, we abbreviate \(I_{k,h}^{\mathrm{VE}}\) as \(I_h^{\mathrm{VE}}\). The following estimate holds (see [42]): for any \(1 \le s \le k+1\) and any \(T \in \mathcal{T}_h\), \[\label{eq:intererrorVE} \|v-I_h^{\mathrm{VE}}v\|_{0,T} + h_T|v-I_h^{\mathrm{VE}}v|_{1,T} \lesssim h_T^{s} |v|_{s,T}, \qquad \forall\, v\in H_0^1(\Omega)\cap H^{s}(\mathcal{T}_h).\tag{26}\] For \(v\in\mathring{H}_1^{\mathrm{nc}}(\mathcal{T}_h)\), let \(v^{\mathrm{CR}}:=I_{1,h}^{\mathrm{VE}}v\in \mathring V_{1,h}^{\mathrm{VE}}\) be the nonconforming linear interpolant of \(v\).

We also need the following estimate relating face jumps to the multiplier.

By the definition of \(I_h^{\mathrm{VE}}\) together with \(Q_{k-1,h}\), the following commuting property holds [20]: \[\label{eq:vem-hdiv-commuting} Q_{k-1,h}\,\nabla_h(I_h^{\mathrm{VE}}v) = Q_{k-1,h}\,\nabla v, \qquad \forall\, v\in H_0^1(\Omega)+\mathring V_{k,h}^{\mathrm{VE}}.\tag{27}\]

We also employ the standard \(H(\operatorname{div})\)-conforming interpolation operator \(I_h^{\operatorname{div}}\) into the BDM space of degree \(k-1\); see, for example, [30], [39], [44]–[46]. For \(\boldsymbol{w}\in H^1(\Omega;\mathbb{R}^d)\), \(v_h\in\mathring V_{k,h}^{\mathrm{VE}}\), and \(k\ge2\), we have (cf. [20]) \[\begin{gather} \label{eq:hdiv-projection-difference} \bigl( Q_{k-1,h}\boldsymbol{w}-I_h^{\operatorname{div}}\boldsymbol{w},\, \nabla_h v_h \bigr) \\ = \sum_{T\in\mathcal{T}_h}\sum_{F\in\Delta_{d-1}(T)} \Bigl( \bigl(Q_{k-1,T}\boldsymbol{w}-\boldsymbol{w}\bigr)\cdot\boldsymbol{n},\, Q_{k-1,F}v_h - Q_{k-2,T}v_h \Bigr)_F . \end{gather}\tag{28}\] In addition, \(I_h^{\operatorname{div}}\) satisfies the commuting relation [47] \[\label{eq:commutationIhdiv} Q_{k-2,h}\operatorname{div}\boldsymbol{w} = \operatorname{div}I_h^{\operatorname{div}}\boldsymbol{w}.\tag{29}\]

3.3 Weak Hessian and norm equivalence↩︎

For \((v,\mu)\in \mathring H_1^{\mathrm{nc}}(\mathcal{T}_h) \times L^2(\mathring{\mathcal{E}}_h)\), we define the weak Hessian \(\nabla_w^2(v,\mu)\in \Sigma_{k,h}^{nn}\) as follows: for all \(\boldsymbol{\tau}\in\Sigma_{k,h}^{nn}\), \[(\nabla_w^2(v,\mu),\boldsymbol{\tau}) = \sum_{T \in \mathcal{T}_h}\Big((v,\operatorname{div}\operatorname{div}\boldsymbol{\tau})_T -(v,\operatorname{tr}_2(\boldsymbol{\tau}))_{\partial T} +\sum_{e\in\Delta_{d-2}(T)} (\mu,\operatorname{tr}_e(\boldsymbol{\tau}))_e\Big).\]

For later use, we introduce two discrete \(H^2\)-seminorms on \(\mathring H_1^{\mathrm{nc}}(\mathcal{T}_h)\times L^2(\mathring{\mathcal{E}}_h)\): \[\begin{align} |\!|\!|(v,\mu)|\!|\!|_{2,\mathrm{CR}}^2 :=& \sum_{T\in\mathcal{T}_h} h_T^{-4}\|Q_{k-2,T}(v-v^{\mathrm{CR}})\|_{0,T}^2 \\ & +\sum_{T\in\mathcal{T}_h}\sum_{F\in\Delta_{d-1}(T)} h_T^{-3}\|Q_{k-1,F}(v-v^{\mathrm{CR}})\|_{0,F}^2 \\ & +\sum_{T\in\mathcal{T}_h}\sum_{e\in\Delta_{d-2}(T)} h_T^{-2}\|Q_{k,e}(v^{\mathrm{CR}}-\mu)\|_{0,e}^2 +\sum_{F\in\mathcal{F}_h} h_F^{-1}\|[\![\partial_n v^{\mathrm{CR}}]\!]\|_{0,F}^2, \end{align}\] and \[\begin{align} |\!|\!|(v,\mu)|\!|\!|_{2,h}^2 :={}& |Q_{k-1,h}\nabla_h v|_{1,h}^2 + \sum_{F\in\mathcal{F}_h} h_F^{-1} \|[\![Q_{k-1,h}\nabla_h v]\!]\|_{0,F}^2 \\ &+ \sum_{T\in\mathcal{T}_h} \sum_{e\in\Delta_{d-2}(T)} h_T^{-2} \|Q_{k,e}(v^{\mathrm{CR}}-\mu)\|_{0,e}^2 . \end{align}\] By the broken Poincaré inequality in [48] and the norm equivalence 25 , the seminorm \(|\!|\!|\cdot|\!|\!|_{2,h}\) is indeed a norm on \(\mathring V_{k,h}^{\mathrm{VE}}\times \mathbb{P}_k(\mathring{\mathcal{E}}_h)\).

4 Symmetric-Tensor Distributional Mixed Methods↩︎

This section presents a distributional mixed method for 1 . The method couples normal–normal continuous symmetric tensor finite elements with an \(H^1\)-nonconforming virtual element space and a codimension-two multiplier. We identify its two-dimensional specialization with an HHJ-type formulation and prove optimal-order and parameter-robust a priori error estimates.

4.1 Distributional mixed method↩︎

For \(k\ge1\), recall from [23] the weak div div operator \((\operatorname{div}\operatorname{div})_w\) from \(\Sigma_{k,h}^{-1}\) to \(\mathring{M}^{-1}_{k-2,k-1,k,k}\): \[(\operatorname{div}\operatorname{div})_w\boldsymbol{\sigma} := ((\operatorname{div}\operatorname{div})_T\boldsymbol{\sigma}, -h_F^{-1}[\operatorname{tr}_2(\boldsymbol{\sigma})]|_F, h_F^{-3}[\boldsymbol{n}^{\intercal}\boldsymbol{\sigma}\boldsymbol{n}]|_F, h_e^{-2}[ \operatorname{tr}_e(\boldsymbol{\sigma})]|_e),\] where \(\operatorname{tr}_2(\boldsymbol{\sigma})\) and \(\operatorname{tr}_e(\boldsymbol{\sigma})\) are defined in ?? , and the broken spaces are \[\begin{align} \Sigma_{k,h}^{-1}&:=\mathbb{P}_k(\mathcal{T}_h;\mathbb{S}), \quad \mathring{M}_{k-2,k-1,k,k}^{-1} := \mathbb{P}_{k-2}(\mathcal{T}_h)\times \mathbb{P}_{k-1}(\mathring{\mathcal{F}}_h)\times \mathbb{P}_k(\mathring{\mathcal{F}}_h) \times \mathbb{P}_k(\mathring{\mathcal{E}}_h). \end{align}\] For the space \(\mathring{M}_{k-2,k-1,k,k}^{-1}\), define the weighted inner product \[\begin{align} ((u_0, u_b, u_n, u_e), (v_0, v_b, v_n, v_e))_{0,h} &:= \sum_{T\in \mathcal{T}_h}(u_0, v_0)_{T} + \sum_{F\in\mathcal{F}_h}h_F(u_b, v_b)_{F} \\ &\quad+ \sum_{F\in\mathcal{F}_h}h_F^3(u_n, v_n)_{F}+ \sum_{e\in\mathcal{E}_h}h_e^2(u_e, v_e)_{e}. \end{align}\] The scaling makes all components dimensionally consistent with the cellwise \(L^2\) inner product \((u_0,v_0)\). Thus, for \(\boldsymbol{\sigma}\in \Sigma_{k,h}^{-1}\) and \(v=(v_0,v_b,v_n,v_e)\in \mathring{M}_{k-2,k-1,k,k}^{-1}\), \[\begin{align} ((\operatorname{div}\operatorname{div})_w\boldsymbol{\sigma}, v)_{0,h} = & \sum_{T\in \mathcal{T}_h}\big((\operatorname{div}\operatorname{div}\boldsymbol{\sigma}, v_0)_T - (\operatorname{tr}_2(\boldsymbol{\sigma}), v_b)_{\partial T}\big) \\ & + \sum_{T\in \mathcal{T}_h}(\boldsymbol{n}^{\intercal} \boldsymbol{\sigma}\boldsymbol{n},v_n\boldsymbol{n}_F\cdot\boldsymbol{n})_{\partial T}+ \sum_{T\in \mathcal{T}_h}\sum_{e\in \Delta_{d-2}(T)}(\operatorname{tr}_e(\boldsymbol{\sigma}), v_e)_e. \end{align}\] If \(\boldsymbol{\sigma}\) is normal–normal continuous across interelement faces, i.e., \(\boldsymbol{\sigma} \in \Sigma^{nn}(\mathcal{T}_h;\mathbb{S})\), this reduces to \[\begin{align} ((\operatorname{div}\operatorname{div})_w\boldsymbol{\sigma}, v)_{0,h} = & \sum_{T\in \mathcal{T}_h}\big((\operatorname{div}\operatorname{div}\boldsymbol{\sigma}, v_0)_T - (\operatorname{tr}_2(\boldsymbol{\sigma}), v_b)_{\partial T}\big) \\ & + \sum_{T\in \mathcal{T}_h}\sum_{e\in \Delta_{d-2}(T)}(\operatorname{tr}_e(\boldsymbol{\sigma}), v_e)_e. \end{align}\]

By the unisolvence of the DoFs 19 for \(V_k^{\mathrm{VE}}(T)\), the map \((Q_{k-2,T}, Q_{k-1,F})_{T,F}\) is an isomorphism from \(\mathring{V}_{k,h}^{\mathrm{VE}}\) onto \(\mathbb{P}_{k-2}(\mathcal{T}_h)\times \mathbb{P}_{k-1}(\mathring{\mathcal{F}}_h)\). This motivates the following distributional mixed method for solving the fourth-order elliptic singular perturbation problem 1 , i.e. a discretization of the mixed formulation 4 : find \(\boldsymbol{\sigma}_h \in \Sigma_{k,h}^{nn}\), \(u_h \in \mathring{V}_{k,h}^{\mathrm{VE}}\), and \(\lambda_h \in \mathbb{P}_k(\mathring{\mathcal{E}}_h)\) such that \[\tag{33} \begin{align} a(\boldsymbol{\sigma}_h, \boldsymbol{\tau}_h) + b_h(\boldsymbol{\tau}_h;u_h,\lambda_h) &= 0, \qquad\qquad\forall\,\boldsymbol{\tau}_h \in \Sigma_{k,h}^{nn}, \tag{34}\\ b_h(\boldsymbol{\sigma}_h; v_h,\mu_h) - c_h(u_h, v_h) &= -(f, \tilde{v}_h),\quad\forall\, v_h \in \mathring{V}_{k,h}^{\mathrm{VE}}, \mu_h\in\mathbb{P}_{k}(\mathring{\mathcal{E}}_h), \tag{35} \end{align}\] where \(\tilde{v}_h := Q_{k-2,h} v_h + (I-Q_{k-2,h}) v_h^{\mathrm{CR}}\), and \[\begin{align} b_h(\boldsymbol{\tau}_h;v_h,\mu_h) &:= -\sum_{T\in\mathcal{T}_h}\Bigl[ (\operatorname{div}\operatorname{div}\boldsymbol{\tau}_h, Q_{k-2,T} v_h)_T -(\operatorname{tr}_2(\boldsymbol{\tau}_h), Q_{k-1,F}v_h)_{\partial T} \\ &\qquad\qquad\quad +\sum_{e\in\Delta_{d-2}(T)}(\operatorname{tr}_e(\boldsymbol{\tau}_h), \mu_h)_e \Bigr], \\ c_h(u_h,v_h) &:= (Q_{k-1,h}\nabla_h u_h,Q_{k-1,h}\nabla_h v_h). \end{align}\]

By the definition of the weak Hessian \(\nabla_w^2\), \[b_h(\boldsymbol{\tau};v,\mu) = -(\boldsymbol{\tau}, \nabla_w^2 (v, \mu)),\quad\forall\,\boldsymbol{\tau} \in \Sigma_{k,h}^{nn}, v \in \mathring H_1^{\mathrm{nc}}(\mathcal{T}_h), \mu \in L^2(\mathring{\mathcal{E}}_h).\] As an immediate consequence of ?? , the following commutative property holds: \[\label{eq:weakdivdivinterproperty} b_h(\boldsymbol{\tau};v - I_h^{\mathrm{VE}} v, \mu - Q_{k,\mathcal{E}_h} \mu) = 0, \quad \forall\,\boldsymbol{\tau} \in \Sigma_{k,h}^{nn}.\tag{36}\]

We next prove the well-posedness of the distributional mixed method 33 .

4.2 Relation to the Hellan–Herrmann–Johnson method↩︎

Now we discuss the relation between the distributional mixed method 33 and the classical Hellan–Herrmann–Johnson (HHJ) method [16]–[18] for solving the fourth-order elliptic singular perturbation problem 1 in two dimensions.

Let the two-dimensional Lagrange finite element space of degree \(k+1\) with homogeneous boundary conditions be \[\mathring{V}_{k+1,h}^{L} := \{ w_h\in H_0^1(\Omega): w_h|_T\in \mathbb{P}_{k+1}(T)\quad \forall\,T\in\mathcal{T}_h\}.\] The space \(\mathring{V}_{k,h}^{\mathrm{VE}}\times \mathbb{P}_k(\mathring{\mathcal{E}}_h)\) is naturally isomorphic to \(\mathring{V}_{k+1,h}^{L}\). This isomorphism is the key step in identifying 33 with a two-dimensional HHJ-type method.

Thus, under the isomorphism \(\mathcal{I}_h\), the bilinear form \(b_h(\boldsymbol{\tau};v,\mu)\) coincides in two dimensions with the HHJ bilinear form \(b_h^{\mathrm{HHJ}}(\boldsymbol{\tau},\mathcal{I}_h(v,\mu))\).

For comparison, the HHJ method for the two-dimensional problem 1 considered in [19] reads as follows: find \((\boldsymbol{\sigma}_h,w_h)\in \Sigma_{k,h}^{nn}\times \mathring{V}_{k+1,h}^{L}\) such that \[\tag{37} \begin{align} \tag{38} \varepsilon^{-2}(\boldsymbol{\sigma}_h,\boldsymbol{\tau}_h) + b_h^{\mathrm{HHJ}}(\boldsymbol{\tau}_h,w_h) &= 0, &&\forall\,\boldsymbol{\tau}_h\in \Sigma_{k,h}^{nn},\\ \tag{39} b_h^{\mathrm{HHJ}}(\boldsymbol{\sigma}_h,v_h) - (\nabla w_h,\nabla v_h) &= -(f,v_h), &&\forall\,v_h\in \mathring{V}_{k+1,h}^{L}. \end{align}\]

Compared with the HHJ method 37 , the HHJ-type formulation ?? uses the projected gradient term \[(Q_{k-1,h}\nabla w_h,Q_{k-1,h}\nabla \chi_h)=(Q_{k-1,h}\nabla_hu_h,Q_{k-1,h}\nabla_hv_h)\] instead of \((\nabla w_h,\nabla v_h)\), and it uses the modified load term \((f,\tilde{\chi}_h)\) associated with the nonconforming scalar variable. The fourth-order HHJ bilinear form is unchanged after identifying \(w_h\) with \((v_h,\mu_h)\) through \(\mathcal{I}_h\). Thus the difference lies in the discretization of the second-order term: ?? , equivalently 33 , uses an \(H^1\)-nonconforming scalar space, whereas 37 uses an \(H^1\)-conforming Lagrange space.

The classical HHJ method with a Lagrange scalar space is intrinsically two-dimensional. By replacing that scalar space with the pair \(\mathring{V}_{k,h}^{\mathrm{VE}}\times\mathbb{P}_k(\mathring{\mathcal{E}}_h)\) and retaining the normal–normal continuous tensor space \(\Sigma_{k,h}^{nn}\), the distributional mixed method 33 extends the HHJ framework to arbitrary spatial dimension.

4.3 Error analysis↩︎

We equip \(\mathring{V}_{k,h}^{\mathrm{VE}}\times \mathbb{P}_{k}(\mathring{\mathcal{E}}_h)\) with the parameter-dependent norm \[|\! |\!|(v_h,\mu_h)|\!|\!|_{\varepsilon,h}^2 := \varepsilon^2|\! |\!|(v_h,\mu_h)|\!|\!|_{2,h}^2 + | v_h |_{1,h}^2.\]

Let \((\boldsymbol{\sigma},u)\in H^{-1}(\operatorname{div}\operatorname{div},\Omega;\mathbb{S}) \times H_0^1(\Omega)\) be the solution of 4 . Whenever the codimension-two traces of \(u\) are well-defined, we define \(\lambda\in L^2(\mathring{\mathcal{E}}_h)\) by \[\lambda|_e:=u|_e\qquad \forall\, e\in\mathring{\mathcal{E}}_h,\] and set \(\boldsymbol{\sigma}_I:=I_h^{nn}\boldsymbol{\sigma}\) whenever \(I_h^{nn}\boldsymbol{\sigma}\) is well-defined. Define the residual functional \[\mathcal{R}_h(v_h,\mu_h) := b_h(\boldsymbol{\sigma}_I;v_h,\mu_h) - c_h(u,v_h) + (f,\widetilde{v}_h).\]

We next estimate the consistency residual associated with 35 .

For the discrete solution \((\boldsymbol{\sigma}_h,u_h,\lambda_h)\) of 33 , define the errors \[e_{\boldsymbol{\sigma}}:=\boldsymbol{\sigma}_I-\boldsymbol{\sigma}_h, \qquad e_u:=I_h^{\mathrm{VE}}u-u_h, \qquad e_\lambda:=Q_{k,\mathcal{E}_h}\lambda-\lambda_h .\]

We next derive a parameter-robust error estimate for the distributional method 33 .

5 Hybridization and Equivalent Formulations↩︎

This section hybridizes the normal–normal continuity of the tensor variable in 33 . After local elimination of the stress variable, the hybridized method yields a stabilization-free weak Galerkin formulation. We then prove that this weak Galerkin formulation is equivalent to an \(H^2\)-nonconforming virtual element method.

5.1 Hybridized normal–normal continuity↩︎

To hybridize the normal–normal continuity of the stress variable, we introduce the multiplier space \(\mathbb{P}_k(\mathring{\mathcal{F}}_h)\) and seek the stress variable in the broken space \[\Sigma_{k,h}^{-1}:=\mathbb{P}_k(\mathcal{T}_h;\mathbb{S}).\] The multiplier imposes normal–normal continuity weakly. The hybridized mixed formulation of 33 reads as follows: find \((\boldsymbol{\sigma}_h,u_h,\gamma_h,\lambda_h) \in \Sigma_{k,h}^{-1}\times \mathring V_{k,h}^{\mathrm{VE}}\times \mathbb{P}_k(\mathring{\mathcal{F}}_h)\times \mathbb{P}_k(\mathring{\mathcal{E}}_h)\) such that \[\tag{51} \begin{align} a(\boldsymbol{\sigma}_h,\boldsymbol{\tau}_h) + b_h^{\mathrm{hyb}}(\boldsymbol{\tau}_h;u_h,\gamma_h,\lambda_h) &=0, \tag{52}\\ b_h^{\mathrm{hyb}}(\boldsymbol{\sigma}_h;v_h,\chi_h,\mu_h) - c_h(u_h,v_h) &=-(f,\widetilde{v}_h), \tag{53} \end{align}\] for all \((\boldsymbol{\tau}_h,v_h,\chi_h,\mu_h) \in \Sigma_{k,h}^{-1}\times \mathring V_{k,h}^{\mathrm{VE}}\times \mathbb{P}_k(\mathring{\mathcal{F}}_h)\times \mathbb{P}_k(\mathring{\mathcal{E}}_h),\) where \[\begin{align} b_h^{\mathrm{hyb}}(\boldsymbol{\tau}_h; v_h,\chi_h,\mu_h) &:= -\sum_{T\in\mathcal{T}_h}\Bigl[ (\operatorname{div}\operatorname{div}\boldsymbol{\tau}_h, Q_{k-2,T} v_h)_T +(\boldsymbol{n}_F^{\intercal}\boldsymbol{\tau}_h\boldsymbol{n}_{\partial T}, \chi_h)_{\partial T} \\ &\qquad\quad -(\operatorname{tr}_2(\boldsymbol{\tau}_h), Q_{k-1,F}v_h)_{\partial T}+\sum_{e\in\Delta_{d-2}(T)}(\operatorname{tr}_e(\boldsymbol{\tau}_h), \mu_h)_e \Bigr]. \end{align}\] If \(\boldsymbol{\tau}_h\in\Sigma_{k,h}^{nn}\), then \[\label{eq:hyb-reduction} b_h^{\mathrm{hyb}}(\boldsymbol{\tau}_h; v_h,\chi_h,\mu_h)=b_h(\boldsymbol{\tau}_h;v_h,\mu_h)\;\;\forall\, v_h \in \mathring{V}_{k,h}^{\mathrm{VE}}, \chi_h\in\mathbb{P}_k(\mathring{\mathcal{F}}_h), \mu_h \in \mathbb{P}_k(\mathring{\mathcal{E}}_h).\tag{54}\]

5.2 Stabilization-free weak Galerkin method↩︎

For \((v,\chi,\mu)\in \mathring H_1^{\mathrm{nc}}(\mathcal{T}_h) \times L^2(\mathring{\mathcal{F}}_h) \times L^2(\mathring{\mathcal{E}}_h)\), we define the hybrid weak Hessian \(\nabla_{w,\mathrm{hyb}}^2(v,\chi,\mu)\in\Sigma_{k,h}^{-1}\) elementwise by \[\begin{align} (\nabla_{w,\mathrm{hyb}}^2(v,\chi,\mu),\boldsymbol{\tau})_T :={}& (Q_{k-2,T}v,\operatorname{div}\operatorname{div}\boldsymbol{\tau})_T + (\chi,\boldsymbol{n}_F^{\intercal}\boldsymbol{\tau} \boldsymbol{n}_{\partial T})_{\partial T} \\ &- (Q_{k-1,F}v,\operatorname{tr}_2(\boldsymbol{\tau}))_{\partial T} + \sum_{e\in\Delta_{d-2}(T)} (\mu,\operatorname{tr}_e(\boldsymbol{\tau}))_e \end{align}\] for all \(\boldsymbol{\tau}\in\mathbb{P}_k(T;\mathbb{S})\) and \(T\in\mathcal{T}_h\).

By definition, for \(\boldsymbol{\tau} \in \Sigma_{k,h}^{-1}\), \(v \in \mathring H_1^{\mathrm{nc}}(\mathcal{T}_h)\), \(\chi\in L^2(\mathring{\mathcal{F}}_h)\) and \(\mu \in L^2(\mathring{\mathcal{E}}_h)\) \[b_h^{\mathrm{hyb}}(\boldsymbol{\tau};v,\chi,\mu) = -(\boldsymbol{\tau},\nabla_{w,\mathrm{hyb}}^2(v,\chi,\mu)).\] Consequently, the first equation 52 is equivalent to \[\boldsymbol{\sigma}_h = \varepsilon^2 \nabla_{w,\mathrm{hyb}}^2(u_h,\gamma_h,\lambda_h).\] Substituting this identity into 53 , we obtain the following stabilization-free weak Galerkin formulation: find \((u_h,\gamma_h,\lambda_h)\in\mathring M_{k,h}^{\mathrm{hyb}}\) such that \[\label{eq:WG} \varepsilon^2 (\nabla_{w,\mathrm{hyb}}^2(u_h,\gamma_h,\lambda_h), \nabla_{w,\mathrm{hyb}}^2(v,\chi,\mu)) + (Q_{k-1,h}\nabla_h u_h,Q_{k-1,h}\nabla_h v) = (f,\widetilde{v})\tag{56}\] for all \((v,\chi,\mu)\in\mathring M_{k,h}^{\mathrm{hyb}}\), where \[\mathring M_{k,h}^{\mathrm{hyb}} := \mathring V_{k,h}^{\mathrm{VE}} \times \mathbb{P}_k(\mathring{\mathcal{F}}_h) \times \mathbb{P}_k(\mathring{\mathcal{E}}_h).\]

5.3 Stabilization-free \(H^2\)-nonconforming virtual element formulation↩︎

We recall the \(H^2\)-nonconforming virtual element space from [42]. For \(k\ge1\), let \[\begin{align} \mathring{W}_{k+2,h}^{\mathrm{VE}} := \{u\in L^2(\Omega):\;& u|_T\in W_{k+2}^{\mathrm{VE}}(T)\; \textrm{ for } T\in\mathcal{T}_h;\;\textrm{ all the DoFs in \eqref{eq:H2NCVEM-DOFs} are } \\ & \text{single-valued on } \mathring{\mathcal{F}}_h \textrm{ and } \mathring{\mathcal{E}}_h,\;\text{and vanish on }\partial\Omega\}, \end{align}\] where the local shape function space is \[\begin{align} W_{k+2}^{\mathrm{VE}}(T):=&\{u\in H^2(T): \Delta^2u\in\mathbb{P}_{k-2}(T);\;\operatorname{tr}_e(\nabla^2u)\in\mathbb{P}_k(e)\, \textrm{ for } e\in\Delta_{d-2}(T); \\ &\quad\; \operatorname{tr}_1(\nabla^2u)|_F\in\mathbb{P}_k(F),\, \operatorname{tr}_2(\nabla^2u)|_F\in\mathbb{P}_{k-1}(F) \, \textrm{ for } F\in\Delta_{d-1}(T)\}. \end{align}\] A unisolvent set of DoFs is given by \[\label{eq:H2NCVEM-DOFs} (Q_{k-2,T}u,\; Q_{k-1,F}u,\; Q_{k,F}(\partial_{n_F}u),\; Q_{k,e}u).\tag{57}\]

Since the first two groups of DoFs in 57 coincide with those of \(\mathring V_{k,h}^{\mathrm{VE}}\), the DoFs in 57 induce the natural isomorphism \(Q_M:\mathring{W}_{k+2,h}^{\mathrm{VE}}\to\mathring M_{k,h}^{\mathrm{hyb}}\) defined by \[Q_Mv_h:= \bigl(I_h^{\mathrm{VE}}v_h,\;Q_{k,F}(\partial_{n_F}v_h),\;Q_{k,e}v_h\bigr),\quad F\in\mathring{\mathcal{F}}_h,\,e\in\mathring{\mathcal{E}}_h.\]

Under the isomorphism \(Q_M\), the weak Galerkin scheme 56 can be rewritten as the following stabilization-free \(H^2\)-nonconforming virtual element method: find \(u_h\in\mathring{W}_{k+2,h}^{\mathrm{VE}}\) such that \[\label{eq:VEMWG} \varepsilon^2 (Q_{k,h}\nabla_h^2 u_h,Q_{k,h}\nabla_h^2 v_h) + (Q_{k-1,h}\nabla_h I_h^{\mathrm{VE}}u_h, Q_{k-1,h}\nabla_h I_h^{\mathrm{VE}}v_h) = (f,\widetilde{v}_h)\tag{58}\] for all \(v_h\in\mathring{W}_{k+2,h}^{\mathrm{VE}}\). By the norm equivalence in [23], \(\|Q_{k,h}\nabla_h^2 v_h\|_0\) is a norm on \(\mathring{W}_{k+2,h}^{\mathrm{VE}}\). Hence 58 is well posed.

6 Numerical Results↩︎

This section reports numerical experiments for the method 33 . The convergence tests are performed on the unit cube \(\Omega=(0,1)^3\) using uniformly refined simplicial meshes. We also include a boundary-layer visualization that illustrates stress concentration near the boundary. All computations are carried out in MATLAB using \(i\)FEM [50].

6.1 Error Estimates for a Smooth Exact Solution↩︎

We first use a smooth exact solution to test the robust and optimal error estimates. For the robust estimate, define \[\mathrm{Err}_1 := \varepsilon^{-1}\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_h\|_{0} + |\!|\!|(u-u_h,\lambda-\lambda_h)|\!|\!|_{\varepsilon,h}.\] For the optimal estimate, we use \[\mathrm{Err}_u := |\!|\!|(I_h^{\mathrm{VE}}u-u_h, Q_{k,\mathcal{E}_h}\lambda-\lambda_h)|\!|\!|_{\varepsilon,h}, \qquad \mathrm{Err}_\sigma := \varepsilon^{-1}\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_h\|_{0}.\]

The results in Table 1 are consistent with the robust estimate ?? . After the initial coarse-mesh regime, the errors decay at the expected first-order rate, \[\mathrm{Err}_1 = O(h).\] Moreover, neither the magnitude of the errors nor the observed rates deteriorate as \(\varepsilon\) decreases.

4pt

The results in Table 2 agree with ?? . In both reported cases, the observed rates on refined meshes confirm the optimal second-order convergence \[\mathrm{Err}_u = O(h^{2}), \qquad \mathrm{Err}_\sigma = O(h^{2}).\]

4pt

6.2 The Reduced-Limit Problem and Boundary Layer↩︎

We next test the reduced-limit estimate ?? using the following computable error: \[\mathrm{Err}_2 := \varepsilon \Bigl( \sum_{T \in \mathcal{T}_h} h_T^{-4} \big\| Q_{k-2,T}(\bar u-u_h-\bar u^{\mathrm{CR}}+u_h^{\mathrm{CR}}) \big\|_{0,T}^2 \Bigr)^{1/2} + \left\|\nabla \bar u-Q_{k-1,h}\nabla_h u_h \right\|_0 .\]

The results in Table 3 agree with ?? . For the reported lowest-order case \(k=1\), the errors exhibit the expected behavior \[\mathrm{Err}_2 = O(h).\] The observed rates are stable with respect to \(\varepsilon\).

4pt

Figure 2 shows that, as \(\varepsilon\) decreases, the large values of this quantity become increasingly localized near the boundary, consistent with the expected boundary-layer behavior.

References↩︎