¹

1 Introduction↩︎

Let \(\Omega\) be a bounded polyhedral domain in \(\mathbb{R}^d\), \(d\ge2\). In this paper, we develop hybridizable staggered discontinuous Galerkin (HSDG) methods for the \(2m\)-th order polyharmonic problem, for \(m\ge1\), \[\label{intro:polyharmonic} \begin{cases} (-1)^m\Delta^m u=f, & \textrm{in }\Omega,\\ u=\partial_n u=\cdots=\partial_n^{m-1}u=0, & \textrm{on }\partial\Omega. \end{cases}\tag{1}\] The method is designed for general polytopal meshes, arbitrary order \(m\geq 1\), arbitrary dimension \(d\geq 2\), and all polynomial degrees \(k\ge0\). Its main ingredients are a staggered primal–dual mesh, a mixed tensor formulation, local tensor enrichments for low-order spaces, and a hybridization procedure that reduces the global system to scalar and face unknowns.

The starting point is the mixed formulation \[\label{intro:mixfem} \sigma=\nabla^m u,\qquad (-1)^m\operatorname{div}^m\sigma=f .\tag{2}\] Formally, the relation \(\sigma=\nabla^m u\) can be written as \[(\sigma,\tau)=(-1)^m(\operatorname{div}^m\tau,u),\] so the scalar unknown only needs to be in \(L^2(\Omega)\), while the tensor unknown carries the \(H(\operatorname{div}^m)\)-type regularity. For \(m=1\), this is the classical mixed formulation based on \(H(\operatorname{div})\)-conforming spaces, including Raviart–Thomas, Brezzi–Douglas–Marini, Nédélec, and related elements [1]–[5]. For \(m=2\), several \(H(\operatorname{div}\operatorname{div})\)-conforming finite element spaces have recently been developed [6]–[10]. For general \(m\), constructing explicit \(H(\operatorname{div}^m)\)-conforming tensor finite elements is difficult, especially on general polytopes.

The staggered discontinuous Galerkin (SDG) framework provides a different way to impose the mixed regularity. SDG methods, initiated by Chung and Engquist [11], [12], use staggered primal and dual grids to obtain stable and locally conservative discretizations. They have been extended from triangular meshes to polygonal meshes in [13]–[15], and a primal SDG method on polytopal meshes was recently developed in [16]. The present work extends [16] from second-order problems to arbitrary-order polyharmonic equations.

Given a polyhedral primal mesh \(\mathcal{K}_h\), each element \(K\) is split by connecting an interior point of \(K\) to its faces, giving a refined mesh \(\mathcal{K}_h^{\rm R}\). The dual mesh \(\mathcal{K}_h^*\) is formed by grouping refined cells that share a primal face. The key SDG idea is to impose complementary continuity on the primal and dual meshes. The tensor variable \(\sigma_h\) satisfies the required \(H(\operatorname{div}^m)\)-type conformity on dual elements, while the scalar variable \(u_h\) is locally \(H^m\)-conforming on primal elements. Thus both variables are globally discontinuous, but their continuity properties are staggered and complementary. This gives a stable mixed discretization of 2 without constructing globally \(H(\operatorname{div}^m)\)-conforming tensor elements on the primal polytopal mesh.

The scalar space is chosen elementwise as \(\mathbb{P}_{k+m}(K)\) on each primal element. The main difficulty is to construct tensor spaces \(\Sigma_{h,k}\) for an arbitrary number of face-normal directions, order \(m\), dimension \(d\), and polynomial degree \(k\). Two local enrichments are used. The trace enrichment recovers missing normal trace layers when \(k<m-1\). The bubble enrichment provides enough tensor directions for the local stability construction. If \(\nu_K\) is the largest number of pairwise non-parallel face normals of \(K\), then tensor layers up to order \(\lfloor m/\nu_K\rfloor\) are enough.

The method is naturally hybridizable. Following the hybridization framework for mixed methods [17], [18], we relax the primal-face continuity of the tensor traces, use a fully broken tensor space on the refined mesh, and introduce scalar trace unknowns on primal faces. The weak \(m\)-th gradient \(\nabla^m_w\) is defined so that the tensor unknown can be eliminated element by element. After static condensation, the global problem involves only the scalar variable and its face traces and is equivalent to a weak Galerkin formulation without penalty stabilization. Stability comes from the staggered mixed structure and the local tensor trace construction, which yield coercivity of the discrete weak gradient after hybridization.

We now place the proposed method in context. The primal formulation of 1 uses the bilinear form \((\nabla^m u,\nabla^m v)\) and therefore requires an \(H^m\)-conforming discrete space. On simplicial meshes, this leads to highly smooth finite elements, such as those in [19]–[23]. A geometric decomposition of the simplicial lattice was introduced in [24], [25] and used for explicit basis construction and implementation in [26]. These conforming finite element spaces are explicit, but they require high polynomial degree \(2^d (m-1)+1\).

On general polytopes, conforming virtual element methods provide a flexible alternative. The \(H^m\)-conforming virtual elements of arbitrary degree \(k\ge m\) on polytopal meshes were constructed in [27], generalizing earlier two- and three-dimensional constructions [28]–[30]. These methods work on arbitrary polytopes and arbitrary degree, but the local shape functions are virtual and a nontrivial stabilization term is needed.

Nonconforming methods reduce the smoothness requirement. On simplicial meshes, minimal \(H^m\)-nonconforming elements were constructed for \(m\le d\) [31], [32] and later extended to \(m=d+1\) [33]; arbitrary \(m\) and \(d\) can also be handled by interior penalty techniques [34]. The minimal nonconforming elements correspond to the lowest-order case \(k=0\) in our notation. A related minimal bubble-enriched construction was proposed in [35]. On general polytopes, \(H^m\)-nonconforming virtual elements of arbitrary degree were developed in [36], [37]. These VEMs, however, require nontrivial stabilization. The present paper follows a different route: the primary object is a staggered mixed DG method, and its stability is obtained from complementary primal–dual continuity and tensor trace liftings, not from a primal nonconforming space with a stabilization term.

After local tensor elimination, the hybridized SDG method can be interpreted as a weak Galerkin or VEM-like method without penalty stabilization. In this interpretation, the scalar unknown and its face traces determine a local weak \(m\)-th gradient. Unlike standard VEM projections on a single polytopal element, the projection here is enriched from a polynomial space on \(K\) to a piecewise polynomial space on the barycentric split \(K^{\rm R}\). This larger local projection space replaces penalty stabilization.

The analysis is carried out under standard shape-regularity assumptions on the polytopal mesh and on the associated refined subcells. The constants in the stability and error estimates are uniform with respect to the mesh size \(h\) and may depend only on the fixed polynomial degree, the order \(m\), the dimension \(d\), and the mesh regularity parameters. These assumptions are stated precisely in Section 4.1.

The remainder of this article is organized as follows. Section 2 introduces the notation, tensor calculus, and staggered mesh structure. Section 3 proves the tensor polynomial decompositions and local trace constructions. Section 4 constructs the SDG finite element spaces and proves the well-posedness of the mixed method. Section 5 develops the hybridized SDG formulation and proves the error estimates. Section 6 reports numerical experiments supporting the theoretical results.

2 Symmetric Tensors and Decompositions↩︎

In this section, we recall basic notation for tensor spaces, contractions, tensor differential operators, and symmetric tensors. We then use simplicial lattices to index bases of symmetric tensor spaces and to describe polynomial layers relative to a face. Finally, we derive a geometric \(t\)-\(n\) decomposition of symmetric tensors and a face-based spanning decomposition that will be used in the construction of the local tensor spaces.

2.1 Tensors↩︎

Denote the space of \(d\)-dimensional tensors of order \(m\) by \(\mathbb{R}^{d,m}:=(\mathbb{R}^d)^{\otimes m}.\) Let \(\{\boldsymbol{e}_1,\ldots,\boldsymbol{e}_d\}\) be an orthonormal basis of \(\mathbb{R}^d\). Then \[\left\{ \boldsymbol{e}_{i_1}\otimes\boldsymbol{e}_{i_2}\otimes\cdots\otimes\boldsymbol{e}_{i_m} : i_\ell\in\{1,\ldots,d\}, \;\ell=1,\ldots,m \right\}\] is a basis of \(\mathbb{R}^{d,m}\). Hence any \(\tau\in\mathbb{R}^{d,m}\) can be written as \[\tau = \tau_{i_1,\ldots,i_m} \boldsymbol{e}_{i_1}\otimes\boldsymbol{e}_{i_2}\otimes\cdots\otimes\boldsymbol{e}_{i_m},\] where repeated indices are summed over \(1,\ldots,d\). Therefore \(\dim\mathbb{R}^{d,m}=d^m.\) Equivalently, the basis can be represented by the integer lattice \(\{1,\ldots,d\}^m\); see Fig. 1 (a).

Let \(\tau\in\mathbb{R}^{d,m}\) and \(\gamma\in\mathbb{R}^{d,n}\) with \(0\le n\le m\). The contraction \[\tau\mathbin{\lrcorner}\gamma\in\mathbb{R}^{d,m-n}, \quad (\tau\mathbin{\lrcorner}\gamma)_{i_1,\ldots,i_{m-n}} := \tau_{i_1,\ldots,i_{m-n},j_1,\ldots,j_n} \gamma_{j_1,\ldots,j_n}.\] In particular, if \(n=1\), then \(\tau\mathbin{\lrcorner}\gamma\) is the tensor-vector contraction, extending the matrix-vector product. If \(n=m\), then \[\tau:\gamma := \tau\mathbin{\lrcorner}\gamma = \tau_{j_1,\ldots,j_m}\gamma_{j_1,\ldots,j_m}\] defines the Euclidean inner product on \(\mathbb{R}^{d,m}\), and the induced norm is the Frobenius norm.

For \(\tau\in\mathbb{R}^{d,m}\) and \(\boldsymbol{v}\in\mathbb{R}^d\), the outer product \(\tau\otimes\boldsymbol{v}\in\mathbb{R}^{d,m+1}\) is defined by \[(\tau\otimes\boldsymbol{v})_{i_1,\ldots,i_m,j} = \tau_{i_1,\ldots,i_m}v_j.\]

Let \(\Omega\subset\mathbb{R}^d\) be a domain. A tensor function on \(\Omega\) is a tensor whose coefficients are functions on \(\Omega\). For \(s\ge0\), define \(H^s(\Omega;\mathbb{R}^{d,m})\) as the space of tensor functions whose coefficients belong to \(H^s(\Omega)\). Let \(\mathbb{P}_k(T)\) be the space of polynomials of degree at most \(k\) on \(T\), with the convention \(\mathbb{P}_k(T)=\{0\}\) for \(k<0.\) For a finite-dimensional linear space \(V\), define \[\mathbb{P}_k(T;V):=\mathbb{P}_k(T)\otimes V.\] Equivalently, \(\mathbb{P}_k(T;V)\) is the space of \(V\)-valued polynomials of degree at most \(k\) on \(T\).

For \(\tau\in H^1(\Omega;\mathbb{R}^{d,m})\), the divergence is the tensor in \(L^2(\Omega;\mathbb{R}^{d,m-1})\) defined by \[(\operatorname{div}\tau)_{i_1,\ldots,i_{m-1}} = \partial_j\tau_{i_1,\ldots,i_{m-1},j}.\] For \(0\le s\le m\), let \(\operatorname{div}^s\) denote the operator obtained by applying divergence \(s\) times. With \(\nabla=(\partial_1,\ldots,\partial_d)^{\intercal}\), this can be written symbolically as \[\label{eq:tensor-div-symbolic} \operatorname{div}^{s}\tau = \tau\mathbin{\lrcorner}\nabla^s .\tag{3}\]

2.2 Symmetric tensors↩︎

Denote by \(\mathbb{S}^{d,m}\) the space of \(d\)-dimensional symmetric tensors of order \(m\). It is the subspace of \(\mathbb{R}^{d,m}\) consisting of tensors \(\tau\) such that \[\tau_{i_1,\ldots,i_m} = \tau_{i_{\sigma(1)},\ldots,i_{\sigma(m)}} \qquad \text{for all }\sigma\in\mathfrak S_m,\] where \(\mathfrak S_m\) is the symmetric group on \(\{1,\ldots,m\}\).

The symmetrization operator \(\operatorname{sym}:\mathbb{R}^{d,m}\to\mathbb{S}^{d,m}\) is defined by \[\operatorname{sym}(\boldsymbol{v}_1\otimes\cdots\otimes\boldsymbol{v}_m) = \frac{1}{m!} \sum_{\sigma\in\mathfrak S_m} \boldsymbol{v}_{\sigma(1)} \otimes\cdots\otimes \boldsymbol{v}_{\sigma(m)}, \qquad \boldsymbol{v}_1,\ldots,\boldsymbol{v}_m\in\mathbb{R}^d,\] and extended by linearity.

For \(v\in H^m(\Omega)\), define the \(m\)-th gradient by \[\nabla^m v := \nabla^{\otimes m}v = (\underbrace{\nabla\otimes\cdots\otimes\nabla}_{m})v, \quad (\nabla^m v)_{i_1,\ldots,i_m} = \partial_{i_1}\cdots\partial_{i_m}v.\] If \(v\) is sufficiently smooth, mixed partial derivatives commute. Hence \(\nabla^m v\in\mathbb{S}^{d,m}.\)

Let \(\alpha=(\alpha_1,\ldots,\alpha_d)\in\mathbb{N}^d\) be a multi-index, with \(|\alpha|=\sum_{i=1}^d\alpha_i\). For \(v\in C^m\), the \(m\)-th derivatives can also be indexed by \[D^\alpha v := \frac{\partial^m v}{\partial x_1^{\alpha_1}\cdots\partial x_d^{\alpha_d}}, \qquad |\alpha|=m.\] This multi-index notation gives a simplicial lattice representation of symmetric tensors.

2.3 Simplicial lattice representation↩︎

A simplicial lattice [24], [38] of degree \(r\) and dimension \(d\) is \[\mathbb{T}^{d}_r = \left\{ \alpha=(\alpha_0,\alpha_1,\ldots,\alpha_d)\in\mathbb{N}^{d+1} :\alpha_0+\alpha_1+\cdots+\alpha_d=r \right\}.\] Each \(\alpha\in\mathbb{T}^d_r\) is called a node of the lattice. Figure 1 (b) shows \(\mathbb{T}^2_3\).

For the symmetric tensor \(\operatorname{sym}(\boldsymbol{e}_{i_1}\otimes\cdots\otimes\boldsymbol{e}_{i_m})\), the order of the factors in \(\boldsymbol{e}_{i_1}\otimes\cdots\otimes\boldsymbol{e}_{i_m}\) is irrelevant up to permutation. Only the number of occurrences of each basis vector is relevant. Let \(\alpha_k\) be the number of factors equal to \(\boldsymbol{e}_k\) for \(k=1,\ldots,d\). Then \(\alpha=(\alpha_1,\ldots,\alpha_d)\) satisfies \(|\alpha|=m\). This set of multi-indices is identified with the lattice \(\mathbb{T}^{d-1}_m\). Therefore \[\left\{ \operatorname{sym}\!\left(\vcenter{\scalebox{1.5}{\otimes}}_{k=1}^d \boldsymbol{e}_k^{\alpha_k}\right) :\alpha\in\mathbb{T}^{d-1}_m \right\},\] where \(\boldsymbol{e}_k^{\alpha_k}:=\boldsymbol{e}_k^{\otimes\alpha_k}\), is a basis of \(\mathbb{S}^{d,m}\). Hence \[\dim\mathbb{S}^{d,m} = |\mathbb{T}^{d-1}_m| = \binom{m+d-1}{m},\] which is much smaller than \(\dim\mathbb{R}^{d,m}=d^m\).

2.4 Layer decomposition of the simplicial lattice↩︎

Let \(F\) be a face of a \(d\)-simplex \(T\). We label the barycentric coordinate opposite to \(F\) by index \(0\). Then each lattice point \(\alpha\in\mathbb{T}^d_k\) can be written as \[\alpha = (\alpha_0,\alpha_1,\ldots,\alpha_d) = \alpha_0\boldsymbol{e}_0+(0,\alpha_F),\] where \(\boldsymbol{e}_0=(1,0,\ldots,0)\) and \(\alpha_F=(\alpha_1,\ldots,\alpha_d)\in\mathbb{T}^{d-1}_{k-\alpha_0}.\) Define the distance of \(\alpha\) to \(F\) by \[\operatorname{dist}(\alpha,F):=\alpha_0, \qquad 0\le \alpha_0\le k.\] For \(0\le j\le k\), define the \(j\)-th layer by \[L_k^d(F,j) := \left\{ j\boldsymbol{e}_0+(0,\alpha_F): \alpha_F\in\mathbb{T}^{d-1}_{k-j} \right\}.\] Thus \(L_k^d(F,j)\) is the set of lattice points in \(\mathbb{T}^d_k\) at distance \(j\) from \(F\). The union of the layers \(L_k^d(F,j)\) for \(j=\ell,\ldots,k\) is a translated copy of the smaller lattice \(\mathbb{T}^d_{k-\ell}\); see Fig. 2 (a).

2.5 Layer decomposition of polynomial spaces↩︎

Let \(T\) be a \(d\)-simplex with barycentric coordinates \(\{\lambda_i, i=0,\ldots, d\}\). We will apply the layer decomposition ?? to the Bernstein basis of \(\mathbb{P}_k(T)\), \[\mathbb{P}_k(T) = \operatorname{span}\{\lambda^\alpha:= \lambda_0^{\alpha_0}\lambda_1^{\alpha_1}\cdots\lambda_d^{\alpha_d} :\alpha\in\mathbb{T}^d_k\}.\] Since \(\lambda_0\) is the barycentric coordinate opposite to \(F\), the \(j\)-th lattice layer gives \[\operatorname{span}\{\lambda^\alpha:\alpha\in L_k^d(F,j)\} = \lambda_0^j\, \mathbb{P}_{k-j}(F).\]

2.6 Layer decomposition of symmetric tensors↩︎

We use the layer decomposition of the lattice \(\mathbb{T}^{d-1}_m\) to decompose the symmetric tensor space \(\mathbb{S}^{d,m}\) with respect to a face \(F\).

Choose a local \(t\)-\(n\) coordinate system associated with \(F\), \[(\boldsymbol{t}_1,\ldots,\boldsymbol{t}_{d-1},\boldsymbol{n}_F),\] where \(\boldsymbol{t}_1,\ldots,\boldsymbol{t}_{d-1}\) form an orthonormal basis of the tangent space of \(F\) and \(\boldsymbol{n}_F\) is a unit normal vector.

For a tangential multi-index \(\alpha_F=(\alpha_1,\ldots,\alpha_{d-1})\), define \[\boldsymbol{t}^{\alpha_F} := \vcenter{\scalebox{1.5}{\otimes}}_{i=1}^{d-1}\boldsymbol{t}_i^{\alpha_i}.\] Using the \(t\)-\(n\) basis, the symmetric tensor space is spanned by \[\mathbb{S}^{d,m} = \operatorname{span} \left\{ \operatorname{sym}\big( \boldsymbol{n}_F^\ell\otimes\boldsymbol{t}^{\alpha_F} \big) : 0\le\ell\le m,\; \alpha_F\in\mathbb{T}^{d-2}_{m-\ell} \right\}.\] Define the \(\ell\)-th normal layer by \[\mathbb{S}^{(\ell)}_{F,m} := \operatorname{span} \left\{ \operatorname{sym}\big( \boldsymbol{n}_F^\ell\otimes\boldsymbol{t}^{\alpha_F} \big) : \alpha_F\in\mathbb{T}^{d-2}_{m-\ell} \right\}, \quad 0\le\ell\le m.\]

2.7 Geometric decomposition using all faces↩︎

If \(K\) is a bounded full-dimensional polytope, then \[d \le \nu_K \le \#\{F:\;F\subset\partial K\}.\] Indeed, if the face normals had fewer than \(d\) pairwise non-parallel directions, then they would span a proper subspace of \(\mathbb{R}^d\). A nonzero vector orthogonal to this span would be tangent to all supporting hyperplanes of \(K\), contradicting the boundedness of \(K\). The upper bound follows from the definition. The lower bound \(\nu_K=d\) is achieved by parallelepipeds, while for a \(d\)-simplex one has \(\nu_K=d+1\).

For a fixed \(K\), write \(\nu:=\nu_K\). Then \[\nu(\left\lfloor\frac{m}{\nu}\right\rfloor+1)>m.\] This elementary inequality will be used in the proof below.

3 Finite Elements for Symmetric Tensors↩︎

In this section, we construct local finite element spaces for symmetric tensor fields. We first introduce a scalar single-face finite element and its normal Taylor decomposition. We then combine this scalar decomposition with the \(t\)-\(n\) decomposition of symmetric tensors to describe normal trace spaces and bubble spaces. Finally, we assemble these single-face spaces on a polytopal element.

3.1 Single-face finite element↩︎

Let \(T\subset\mathbb{R}^d\) be a single-face cell with distinguished face \(F\) and vertex \(c\) opposite to \(F\). Let \(\lambda_c\) be the affine function satisfying \[\lambda_c(c)=1, \qquad \lambda_c|_F=0.\] Then \(\nabla\lambda_c\perp F\) and \[\nabla\lambda_c=(\partial_n\lambda_c)\boldsymbol{n}_F, \qquad \partial_n\lambda_c=\frac{1}{\operatorname{dist}(c,F)}.\]

Although \(T=\operatorname{conv}(c,F)\) need not be a simplex, we may use simplex coordinates as auxiliary affine coordinates. Choose a nondegenerate \((d-1)\)-simplex \(S_F\subset F\) and let \(\widehat T_F:=\operatorname{conv}(c,S_F).\) Let \(\lambda_c,\lambda_1,\ldots,\lambda_d\) be the barycentric coordinates of \(\widehat T_F\). These functions are affine on the whole space \(\mathbb{R}^d\). Therefore their restrictions to \(T\) may be used to represent polynomials on \(T\). In particular, \[\mathbb{P}_k(T) = \operatorname{span}\{\lambda^\alpha:\alpha\in\mathbb{T}^d_k\},\] where \[\lambda^\alpha = \lambda_c^{\alpha_0}\lambda_F^{\alpha_F}, \qquad \lambda_F^{\alpha_F} = \lambda_1^{\alpha_1}\cdots\lambda_d^{\alpha_d}.\] We identify a polynomial on \(F\) with its extension to \(T\) obtained by using the same monomials \(\lambda_F^{\alpha_F}\), now viewed as affine functions on \(T\). With this convention, for \(k\ge0\) and \(i\ge0\), the polynomial space has the normal-layer decomposition \[\label{eq:Pkgeodecomp1} \mathbb{P}_k(T) = \left( \vcenter{\scalebox{1.5}{\oplus}}_{j=0}^{i-1}\lambda_c^j\mathbb{P}_{k-j}(F) \right) \vcenter{\scalebox{1.5}{\oplus}} \lambda_c^i\mathbb{P}_{k-i}(T),\tag{4}\] with the convention \(\mathbb{P}_s=\{0\}\) if \(s<0\).

Introduce the single-face trace operator \[\operatorname{tr}_F^{H^m}v = (v,\partial_n v,\ldots,\partial_n^{m-1}v)|_F.\] Then 4 gives \[\mathbb{B}_k^{H^m}(T) := \ker(\operatorname{tr}_F^{H^m})\cap\mathbb{P}_k(T) = \lambda_c^m\mathbb{P}_{k-m}(T).\] When \(k<m\), \(\mathbb{B}_k^{H^m}(T)=\{0\}\).

3.2 Polynomial symmetric tensors and normal traces↩︎

Combining the polynomial decomposition 4 with the \(t\)-\(n\) decomposition of \(\mathbb{S}^{d,m}\), we obtain \[\label{eq:PkSdecomp} \mathbb{P}_k(T;\mathbb{S}^{d,m}) = \vcenter{\scalebox{1.5}{\oplus}}_{i=0}^{k} \vcenter{\scalebox{1.5}{\oplus}}_{\ell=0}^{m} \lambda_c^i\mathbb{P}_{k-i}(F;\mathbb{S}^{(\ell)}_{F,m}).\tag{5}\] This decomposition is indexed by pairs \((i,\ell)\), where \(i\) is the polynomial layer in the normal direction and \(\ell\) is the number of normal tensor factors; see Fig. 3.

Define the normal trace component \[\gamma_s(\tau) := ((\operatorname{div}^s\tau)\mathbin{\lrcorner}\boldsymbol{n}_F)|_F, \qquad 0\le s\le m-1.\] The full normal trace is \[\operatorname{tr}_F^{\operatorname{div}^m}\tau := (\gamma_0(\tau),\gamma_1(\tau),\ldots,\gamma_{m-1}(\tau)).\] For \(\tau=a(x)\eta\) with a scalar function \(a\) and a constant tensor \(\eta\in\mathbb{R}^{d,m}\), we have \[(\operatorname{div}^s\tau)_{i_1\cdots i_{m-s}} = \partial_{j_1}\cdots\partial_{j_s} \big(a\,\eta_{i_1\cdots i_{m-s}j_1\cdots j_s}\big) = (\partial_{j_1}\cdots\partial_{j_s}a) \eta_{i_1\cdots i_{m-s}j_1\cdots j_s}.\] Therefore, \[\operatorname{div}^s(a\eta) = \eta\mathbin{\lrcorner}\nabla^s a.\]

Define \[\mathbb{B}_k^{\operatorname{div}^m}(T;\mathbb{S}^{d,m}) := \ker(\operatorname{tr}_F^{\operatorname{div}^m}) \cap \mathbb{P}_k(T;\mathbb{S}^{d,m}).\] Lemma 6 gives the following characterization of the polynomial trace and bubble spaces.

Figure 3 shows the \((i,\ell)\)-indexed block layout. The lower-triangular cells \(i\ge\ell\) form the bubble space, while the upper-triangular cells \(i<\ell\) are detected by the normal traces.

3.3 Single-face enriched tensor space↩︎

The natural local polynomial tensor space is \(\mathbb{P}_k(T;\mathbb{S}^{d,m})\). To obtain stability for all polynomial degrees, we enrich this space by trace and bubble components. The trace enrichment is needed only when \(k<m-1\), while the bubble enrichment depends on the geometry of \(K\).

3.3.1 Trace cells↩︎

When \(k\ge m-1\), the polynomial tensor space already contains all trace layers. When \(k<m-1\), the trace components \(\gamma_i\), \(k+1\le i\le m-1\), are missing because \(\operatorname{div}^i\sigma=0\) for \(\sigma\in\mathbb{P}_k(T;\mathbb{S}^{d,m})\); see the white cells in Fig. 3 (b).

We will add more trace cells row by row and define the full trace-cell space by \[\label{eq:full-trace-space} \operatorname{Tr}(\Sigma_k;T) := \vcenter{\scalebox{1.5}{\oplus}}_{i=0}^{m-1} \lambda_c^i \mathbb{P}_{(k-i)^+}(F)\otimes \operatorname{sym}\big(\boldsymbol{n}_F^{i+1}\otimes\mathbb{S}^{d,m-i-1}\big),\tag{6}\] where \((k-i)^+ = \max\{k-i,0\}\).

3.3.2 Bubble cells↩︎

To prove the inf-sup condition, the local tensor space must contain enough bubble tensors to recover all tensors in \(\mathbb{P}_k(K;\mathbb{S}^{d,m})\). This requirement depends on the geometry of \(K\) through the number \(\nu_K\) of pairwise non-parallel face normals.

For fixed \(K\), write \(\nu:=\nu_K\). If \(\lfloor m/\nu\rfloor>0\), we add extra bubble cells so that all tensor layers up to order \(\lfloor m/\nu\rfloor\) carry degree \(k\). Define \[\label{eq:full-bubble-space} \begin{align} \mathbb{B}(\Sigma_k;T) :={}& \underbrace{ \vcenter{\scalebox{1.5}{\oplus}}_{\ell=0}^{\lfloor m/\nu\rfloor} \lambda_c^{\ell}\mathbb{P}_k(T;\mathbb{S}^{(\ell)}_{F,m}) }_{\mathbb{B}_{\mathrm{ess}}(T)} \vcenter{\scalebox{1.5}{\oplus}} \underbrace{ \vcenter{\scalebox{1.5}{\oplus}}_{\ell=\lfloor m/\nu\rfloor+1}^{m} \lambda_c^{\ell}\mathbb{P}_{k-\ell}(T;\mathbb{S}^{(\ell)}_{F,m}) }_{\mathbb{B}_{\mathrm{rem}}(T)} . \end{align}\tag{7}\]

3.3.3 Stress space and unisolvence↩︎

We define the local single-face stress space by the trace–bubble decomposition \[\label{eq:decomposition-tr-bubble} \Sigma_k(T) := \operatorname{Tr}(\Sigma_k;T) \vcenter{\scalebox{1.5}{\oplus}} \mathbb{B}(\Sigma_k;T).\tag{8}\]

The decomposition 8 gives the following degrees of freedom.

3.4 Finite element for symmetric tensors on a polytope↩︎

We now fix one polytope \(K\) and construct the local tensor space on \(K\). Let \(K^{\rm R}:=\{T_i\}_{i=1}^{n_K}\) be the local refinement associated with the faces \(F_i\subset\partial K\). We assume throughout that such a point \(c_K\) exists and that each \(T_F\) is contained in \(K\). Define the broken local stress space by \[\label{eq:broken-local-stress-space} \Sigma^{-1}_{k}(K^{\rm R}) := \prod_{T\in K^{\rm R}}\Sigma_k(T) = \operatorname{Tr}(\Sigma_k;K) \vcenter{\scalebox{1.5}{\oplus}} \mathbb{B}_k(\Sigma_k;K),\tag{9}\] where \[\operatorname{Tr}(\Sigma_k;K) := \prod_{T_F\in K^{\rm R}}\operatorname{Tr}(\Sigma_k;T_F), \qquad \mathbb{B}_k(\Sigma_k;K) := \prod_{T\in K^{\rm R}}\mathbb{B}(\Sigma_k;T).\] For \(\tau\in\mathbb{B}_k(\Sigma_k;K)\), the exterior traces satisfy \[\gamma_i(\tau)|_{\partial K}=0,\qquad i=0,\ldots,m-1.\] The trace condition is imposed only on the exterior faces \(F_i\subset\partial K\), not on the internal faces of the local refinement.

By Lemma 7 on each refined cell \(T\), we obtain the following DoFs on \(K\).

To prove the discrete inf-sup condition, we use an equivalent set of interior DoFs. Let \[N_m:=\dim\mathbb{S}^{d,m}=\binom{m+d-1}{m}.\] By Lemma 4, we can choose faces \(F_\alpha\subset\partial K\), integers \(0\le r_\alpha\le \lfloor m/\nu\rfloor\), and tensors \[\{\xi_\alpha\in\mathbb{S}^{(r_\alpha)}_{F_\alpha,m}, \, \alpha=1,\ldots,N_m,\} \text{ is a basis of } \mathbb{S}^{d,m}.\] Let \(T_\alpha=T_{F_\alpha}\) and \(\lambda_\alpha\) be the affine function on \(T_\alpha\) satisfying \(\lambda_\alpha|_{F_\alpha}=0\) and \(\lambda_\alpha(c_K)=1\).

Define the essential bubble space on \(K\) as \[\label{eq:essential-bubble-space} \mathbb{B}_{\rm ess}(\Sigma_k;K) := \vcenter{\scalebox{1.5}{\oplus}}_{\alpha=1}^{N_m} \left\{ \lambda_\alpha^{r_\alpha} q_\alpha \xi_\alpha\,\chi_{T_\alpha}: q_\alpha\in\mathbb{P}_k(T_\alpha) \right\} \subset \mathbb{B}_k(\Sigma_k;K),\tag{10}\] where \(\chi_{T_\alpha}\) is the characteristic function of \(T_\alpha\). Let \(\mathbb{B}_{\rm rem}(\Sigma_k;K)\) be the \(L^2(K)\)-orthogonal complement of \(\mathbb{B}_{\rm ess}(\Sigma_k;K)\) in \(\mathbb{B}_k(\Sigma_k;K)\). Then \[\label{eq:bubble-essential-rem-split} \mathbb{B}_k(\Sigma_k;K) = \mathbb{B}_{\rm rem}(\Sigma_k;K) \vcenter{\scalebox{1.5}{\oplus}}^{\bot_{L^2(K)}} \mathbb{B}_{\rm ess}(\Sigma_k;K).\tag{11}\] We will change the DoF corresponding to \(\mathbb{B}_{\rm rem}(\Sigma_k;K)\) to \(\mathbb{P}_k(K;\mathbb{S}^{d,m})\).

4 Staggered Discontinuous Galerkin Methods↩︎

In this section, we define the staggered discontinuous Galerkin method for the polyharmonic problem. The tensor variable is normal-trace continuous across primal faces, while the scalar variable is locally \(H^m\) on each primal element.

4.1 Primal and dual meshes↩︎

Let \(\mathcal{K}_h\) be a shape-regular polytopal partition of \(\Omega\subset\mathbb{R}^d\). Let \(\mathcal{F}_h\) be the set of faces, \(\mathcal{F}_h^\partial\) the set of boundary faces, and \(\mathring{\mathcal{F}}_h:=\mathcal{F}_h\setminus\mathcal{F}_h^\partial\) the set of interior faces. We call \(\mathcal{K}_h\) the primal mesh.

For each \(K\in\mathcal{K}_h\), choose a point \(c_K\in K\). For each face \(F\subset\partial K\), define \[T_F:=\operatorname{conv}(c_K,F), \qquad K^{\rm R}:=\{T_F:\;F\subset\partial K\}.\] The collection of all such cells gives a refinement \(\mathcal{K}_h^{\rm R}\) of \(\mathcal{K}_h\). We assume that the points \(c_K\) are chosen so that \(\mathcal{K}_h^{\rm R}\) is uniformly shape regular, and that the number of faces of each element is uniformly bounded.

Let \(\mathcal{F}_h^{\rm R}\) be the set of faces of the refined mesh \(\mathcal{K}_h^{\rm R}\) and define \(\mathring{\mathcal{F}}_h^* := \mathcal{F}_h^{\rm R}\setminus\mathcal{F}_h.\) The faces in \(\mathring{\mathcal{F}}_h^*\) are the dual faces; they lie inside primal elements. For each primal face \(F\in\mathcal{F}_h\), let \(\omega_F\) be the union of refined cells in \(\mathcal{K}_h^{\rm R}\) that contain \(F\). The dual mesh is \[\mathcal{K}_h^* := \{\omega_F:\;F\in\mathcal{F}_h\}.\]

For each face \(F\), fix a unit normal \(\boldsymbol{n}_F\). When no confusion can arise, we write \(\boldsymbol{n}\) for the relevant unit normal. For an interior face \(F\), the jump \([\cdot]\) is taken with respect to \(\boldsymbol{n}_F\). On a boundary face, we use the convention \([w]=w\).

4.2 Broken spaces and weak operators↩︎

Define the broken scalar space \[H^m(\mathcal{K}_h) := \{v\in L^2(\Omega):\;v|_K\in H^m(K)\quad\text{for all }K\in\mathcal{K}_h\}.\] For the tensor variable, define \[H(\operatorname{div}^m,\mathcal{K}_h^*;\mathbb{S}^{d,m}) := \left\{ \tau\in L^2(\Omega;\mathbb{S}^{d,m}): \tau|_{K^*}\in H(\operatorname{div}^m,K^*;\mathbb{S}^{d,m}) \;\text{for all }K^*\in\mathcal{K}_h^* \right\}.\] Equivalently, \[H(\operatorname{div}^m,\mathcal{K}_h^*;\mathbb{S}^{d,m}) = \left\{ \tau\in H^m(\mathcal{K}_h^{\rm R};\mathbb{S}^{d,m}): [\operatorname{tr}_F^{\operatorname{div}^m}\tau]=0 \quad\text{for all }F\in\mathring{\mathcal{F}}_h \right\},\] where \[\operatorname{tr}_F^{\operatorname{div}^m}\tau = (\gamma_0(\tau),\gamma_1(\tau),\ldots,\gamma_{m-1}(\tau)), \; \gamma_i(\tau) = ((\operatorname{div}^i\tau)\mathbin{\lrcorner}\boldsymbol{n}_F)|_F, \; i = 0, \ldots, m-1.\] Thus the normal trace components are single-valued across primal faces. The tensor variable may jump across dual faces, while the scalar variable may jump across primal faces.

For \(u\in H^m(\mathcal{K}_h)\), define the weak gradient by \[\label{eq:weakgrad} \begin{align} \langle \nabla_w^m u,\tau\rangle :=\;& \sum_{K\in\mathcal{K}_h} (\nabla^m u,\tau)_K\\ &+ \sum_{j=0}^{m-1} (-1)^{m-j} \sum_{F\in\mathcal{F}_h} \left( [\nabla^j u], (\operatorname{div}^{m-1-j}\tau)\mathbin{\lrcorner}\boldsymbol{n} \right)_F, \end{align}\tag{13}\] for all \(\tau\in H(\operatorname{div}^m,\mathcal{K}_h^*;\mathbb{S}^{d,m})\). The jump terms in 13 appear only on primal faces, because \(u\) is single-valued and \(H^m\) inside each primal element.

For \(\sigma\in H(\operatorname{div}^m,\mathcal{K}_h^*;\mathbb{S}^{d,m})\), define the weak divergence by \[\label{eq:weakdiv} \begin{align} \langle \operatorname{div}_w^m\sigma,v\rangle :=\;& \sum_{K^*\in\mathcal{K}_h^*} (\operatorname{div}^m\sigma,v)_{K^*}\\ &+ \sum_{i=0}^{m-1} (-1)^{m-i} \sum_{F^*\in\mathring{\mathcal{F}}_h^*} \left( [(\operatorname{div}^{i}\sigma)\mathbin{\lrcorner}\boldsymbol{n}], \nabla^{m-i-1}v \right)_{F^*}, \end{align}\tag{14}\] for all \(v\in H^m(\mathcal{K}_h)\). The jump terms in 14 appear only on dual faces, because the corresponding normal traces of \(\sigma\) are single-valued on primal faces.

For \(k\ge0\), define the scalar space \[U_{h,k+m} := \{u_h\in L^2(\Omega):\;u_h|_K\in\mathbb{P}_{k+m}(K) \quad\text{for all }K\in\mathcal{K}_h\}.\] Let \(\Sigma_k(T)\) be the local tensor space 8 on each refined cell \(T\in K^{\rm R}\). Define \[\Sigma_k^{-1}(K^{\rm R}) := \prod_{T\in K^{\rm R}}\Sigma_k(T).\] The global tensor space is \[\begin{align} \Sigma_{h,k} := \big\{ \sigma_h\in L^2(\Omega;\mathbb{S}^{d,m}): &\;\sigma_h|_K\in\Sigma_k^{-1}(K^{\rm R}) \quad\text{for all }K\in\mathcal{K}_h, \notag\\ &\;\text{the trace DoFs in \eqref{Ktensordof1} are single-valued across every } F\in\mathring{\mathcal{F}}_h \big\}. \label{eq:global-stress-space} \end{align}\tag{15}\] Then \[\Sigma_{h,k} \subset H(\operatorname{div}^m,\mathcal{K}_h^*;\mathbb{S}^{d,m}).\]

4.3 The Staggered Discontinuous Galerkin method↩︎

For \(\sigma_h,\tau_h\in\Sigma_{h,k}\) and \(v_h\in U_{h,k+m}\), define \[\begin{align} a(\sigma_h,\tau_h) &:= (\sigma_h,\tau_h), \\ b_h(\sigma_h,v_h) &:= (-1)^{m-1} \langle \operatorname{div}_w^m\sigma_h,v_h\rangle = -\langle \nabla_w^m v_h,\sigma_h\rangle. \end{align}\] Equivalently, by 13 , \[\label{eq:bh-expanded} \begin{align} b_h(\sigma_h,v_h) ={}& -\sum_{K\in\mathcal{K}_h} (\sigma_h,\nabla^m v_h)_K \\ &+ \sum_{i=0}^{m-1} (-1)^i \sum_{F\in\mathcal{F}_h} \left( (\operatorname{div}^i\sigma_h)\mathbin{\lrcorner}\boldsymbol{n}, [\nabla^{m-i-1}v_h] \right)_F . \end{align}\tag{16}\]

We reindex \(i\) by \(j\) with relation \[\label{eq:kj-def} j=m-i-1, k_j:=\max\{k-(m-1-j),0\}=(k-i)^+.\tag{17}\] This index \(k_j\) is the effective polynomial degree of the stress trace paired with \([\nabla^j v_h]\) in 16 . Indeed, with \(j=m-i-1\), the jump \([\nabla^j v_h]\) is paired with \[\gamma_{m-1-j}(\sigma_h) = (\operatorname{div}^{m-1-j}\sigma_h)\mathbin{\lrcorner}\boldsymbol{n},\] whose degree is at most \(k_j\). Hence only the projected jump \(Q_{F,k_j}[\nabla^j v_h]\) is seen by \(b_h(\cdot,\cdot)\). Here \(Q_{F,r}\) is the \(L^2(F)\) projection onto \(\mathbb{P}_r(F)\), applied componentwise.

For \(v\in H^m(\mathcal{K}_h)\), define \[\label{eq:U-norm} \|v\|_{m,h}^2 := |v|_{H^m(\mathcal{K}_h)}^2 + \sum_{j=0}^{m-1} \sum_{F\in\mathcal{F}_h} h_F^{2j-2m+1} \left\| Q_{F,k_j} [\nabla^j v] \right\|_F^2.\tag{18}\] By the Poincaré inequality for broken \(H^1\) spaces in [39], \[\label{eq:hybrid-poincare} \|v\|_{H^{m-1}(\mathcal{K}_h)} \lesssim \|v\|_{m,h} \qquad \forall~v\in H^m(\mathcal{K}_h).\tag{19}\] Thus \(\|\cdot\|_{m,h}\) is a norm on \(H^m(\mathcal{K}_h)\).

For \(\tau_h\in\Sigma_{h,k}\), define \[\label{eq:Sigma-norm} \|\tau_h\|_{0,h}^2 := \|\tau_h\|^2 + \sum_{i=0}^{m-1} \sum_{F\in\mathcal{F}_h} h_F^{2i+1} \left\| (\operatorname{div}^i\tau_h)\mathbin{\lrcorner}\boldsymbol{n} \right\|_F^2.\tag{20}\] By inverse trace inequalities on the finite-dimensional local tensor spaces, \[\label{eq:sigmaJ-equivalence} \|\tau_h\|_{0,h} \eqsim \|\tau_h\| \qquad \forall~\tau_h\in\Sigma_{h,k}.\tag{21}\]

The SDG method is: find \((\sigma_h,u_h)\in\Sigma_{h,k}\times U_{h,k+m}\) such that \[\label{SDG} \begin{cases} a(\sigma_h,\tau_h)+b_h(\tau_h,u_h)=0, & \quad \forall~\tau_h\in\Sigma_{h,k},\\ b_h(\sigma_h,v_h)=-(f,v_h), & \quad \forall~v_h\in U_{h,k+m}. \end{cases}\tag{22}\]

By 21 , \(a(\cdot,\cdot)\) is coercive on \(\Sigma_{h,k}\) with respect to \(\|\cdot\|_{0, h}\): \[\label{eq:acoercive} a(\tau_h,\tau_h) = \|\tau_h\|^2 \gtrsim \|\tau_h\|_{0, h}^2.\tag{23}\] Moreover, based on 16 , using the Cauchy-Schwarz inequality with correct scaling, we have \[\label{eq:bcontinuity} |b_h(\sigma_h,v)| \lesssim \|\sigma_h\|_{0, h}\|v\|_{m, h} \qquad \forall~\sigma_h\in\Sigma_{h,k},\quad v\in H^m(\mathcal{K}_h).\tag{24}\]

The key is to verify the following inf-sup condition.

Then by the Babuška–Brezzi theory, we have the following well-posedness.

5 Hybridization↩︎

This section hybridizes the staggered DG method by using a fully broken stress space and scalar trace unknowns only on primal faces. The resulting mixed method is stable, and local elimination of the stress variable gives an equivalent stabilization-free weak Galerkin formulation. An optimal-order error analysis is given in the energy norm.

5.1 Broken spaces and weak operators↩︎

Following the standard hybridization idea [17], we use a fully broken stress space on the refined mesh \(\mathcal{K}_h^{\rm R}\) and scalar trace unknowns on primal faces. For \(j=0,\ldots,m-1\), recall that \(k_j=\max\{k-m+1+j,0\}\) is the degree of the stress trace paired with the \(j\)-th scalar trace. Define the hybrid scalar space by \[\begin{align} M_h := \big\{ u_h={}&(u_0,u_{g,0},\ldots,u_{g,m-1}):\; u_0|_K\in\mathbb{P}_{k+m}(K) \quad \text{for all }K\in\mathcal{K}_h,\\ &\quad u_{g,j}|_F \in \mathbb{P}_{k_j}(F;\mathbb{S}^{d,j}) \quad \text{for all }F\in\mathcal{F}_h,\; j=0,\ldots,m-1 \big\}. \end{align}\] Thus only the traces on primal faces are independent hybrid unknowns. Let \[M_h^0 := \{u_h\in M_h:\;u_{g,j}|_{\partial\Omega}=0,\quad j=0,\ldots,m-1\}.\]

The stress space is fully broken on the refined mesh: \[\label{eq:broken-hybrid-tensor-space} \Sigma_{h,k}^{-1} := \prod_{T\in\mathcal{K}_h^{\rm R}}\Sigma_k(T).\tag{27}\]

For \(u_h\in M_h\), define \(\nabla_w^m u_h\in\Sigma_{h,k}^{-1}\) cellwise. Let \(T=T_F\subset K\) be the refined cell attached to the primal face \(F\subset\partial K\). Then \(\nabla_{w,T}^m u_h\in\Sigma_k(T)\) is defined by \[\label{eq:hybrid-weak-gradient} \begin{align} (\nabla_{w,T}^m u_h,\tau)_T ={}& (\nabla^m u_0,\tau)_T\\ &+ \sum_{j=0}^{m-1} (-1)^{m-1-j} \langle u_{g,j}-\nabla^j u_0, \gamma_{m-1-j}(\tau) \rangle_F , \end{align}\tag{28}\] for all \(\tau\in\Sigma_k(T)\). No independent trace unknowns are introduced on dual faces; equivalently, the traces on dual faces are taken to be the natural traces \(\nabla^j u_0\). Notice that as \(\gamma_{m-1-j}(\tau)\in \mathbb{P}_{k_j}(F;\mathbb{S}^{d,j})\), we can add the projection operator: \[\langle u_{g,j}-\nabla^j u_0, \gamma_{m-1-j}(\tau) \rangle_F = \langle Q_{F, k_j}(u_{g,j}-\nabla^j u_0), \gamma_{m-1-j}(\tau) \rangle_F.\]

5.2 Hybridized variational form↩︎

Define \[a(\sigma_h,\tau_h):=(\sigma_h,\tau_h), \qquad b_h(\tau_h,u_h):=-(\tau_h,\nabla_w^m u_h).\] The hybridized SDG method is: find \((\sigma_h,u_h)\in\Sigma_{h,k}^{-1}\times M_h^0\) such that \[\tag{29} \begin{align} a(\sigma_h,\tau_h)+b_h(\tau_h,u_h)&=0, \quad &&\forall~\tau_h\in\Sigma_{h,k}^{-1},\tag{30}\\ b_h(\sigma_h,v_h)&=-(f,v_0), \quad &&\forall~v_h\in M_h^0.\tag{31} \end{align}\] Since the stress space \(\Sigma_{h,k}^{-1}\) is fully broken on \(\mathcal{K}_h^{\rm R}\), the stress variable can be eliminated cell by cell. Indeed, 30 gives \[(\sigma_h,\tau_h) = (\tau_h,\nabla_w^m u_h) \qquad \forall~\tau_h\in\Sigma_{h,k}^{-1},\] and hence \(\sigma_h=\nabla_w^m u_h.\) Substituting this identity into 31 gives the primal hybridized form: find \(u_h\in M_h^0\) such that \[\label{eq:WG} (\nabla_w^m u_h,\nabla_w^m v_h) = (f,v_0) \qquad \forall~v_h\in M_h^0 .\tag{32}\] Thus the hybridized SDG method is equivalent to a stabilization-free weak Galerkin method. Since \(\sigma_h\) is eliminated locally, the local stress space may be enriched to \(\mathbb{P}_{\bar k}(T;\mathbb{S}^{d,m})\) without adding global unknowns.

For \(u_h\in M_h\), define \[\label{eq:hybrid-u-norm} \begin{align} \|u_h\|_{m,h}^2 := |u_0|_{H^m(\mathcal{K}_h)}^2 + \sum_{j=0}^{m-1} \sum_{K\in\mathcal{K}_h} \sum_{F\subset\partial K} h_F^{2j-2m+1} \|Q_{F,k_j}(u_{g,j}-\nabla^j u_0)\|_F^2 . \end{align}\tag{33}\]

5.3 Error analysis↩︎

All constants below may depend on \(m\), \(d\), \(k\), and shape-regularity parameters, but not on \(h\). Let \[Q_Mu:=(Q_0u,Q_{g,0}u,\ldots,Q_{g,m-1}u)\in M_h^0\] be the componentwise projection. Here \(Q_0u|_K\) is the \(L^2(K)\) projection onto \(\mathbb{P}_{k+m}(K)\), and \(Q_{g,j}u|_F\) is the \(L^2(F)\) projection of \(\nabla^j u\) onto \(\mathbb{P}_{k_j}(F;\mathbb{S}^{d,j})\) on primal faces.

An \(L^2\)-error estimate can be derived by a standard Aubin–Nitsche duality argument. Since this argument is routine, we omit the details.

6 Numerical Experiments↩︎

We present numerical results for the hybridized SDG method 29 applied to the polyharmonic problem 1 . The experiments verify the predicted convergence rates and test robustness on two-dimensional polygonal meshes and three-dimensional tetrahedral meshes.

6.1 Test problems and mesh configurations↩︎

On the unit square \(\Omega=(0,1)^2\), we use the manufactured solutions \[u(x,y)=\sin^m(\pi x)\sin^m(\pi y), \qquad m=2,3.\] The right-hand side is chosen accordingly. The errors are measured by \[\|\sigma-\sigma_h\|_{0,h},\qquad |u-u_0|_{H^m(\mathcal{K}_h)},\qquad \|u-u_0\|.\] Here \(h\) denotes the mesh size and NT denotes the number of elements.

In two dimensions, we test two families of polygonal meshes. The first consists of convex polygonal meshes generated by connecting element centroids of a triangular mesh. The second consists of concave polygonal meshes constructed from a perturbed square grid. Representative meshes are shown in Figure 5.

For the convex polygonal meshes used in Fig 5, the number of pairwise non-parallel face-normal directions is \(\nu=3\). Thus \(\lfloor m/\nu\rfloor=0\) for \(m=2\), while \(\lfloor m/\nu\rfloor=1\) for \(m=3\). For the concave polygonal meshes, we have \(\nu=4\), and hence \(\lfloor m/\nu\rfloor=0\) for both \(m=2\) and \(m=3\). Therefore, bubble enrichment is not needed in these two-dimensional tests, except for the triharmonic experiment on convex polygonal meshes. The low-order cases still require trace enrichment when \(k<m-1\).

6.2 Two-dimensional biharmonic problem↩︎

Tables 1–2 report the errors for the biharmonic problem, corresponding to \(m=2\), on convex and concave polygonal meshes. The energy-type errors converge with order \(k+1\), and the \(L^2\) error \(\|u-u_0\|\) shows higher-order convergence.

6.3 Two-dimensional triharmonic problem↩︎

Tables 3–4 report the errors for the triharmonic problem, corresponding to \(m=3\), on convex and concave polygonal meshes. These tests include trace enrichment for \(k=0,1\), but no bubble enrichment is required because \(\lfloor m/\nu\rfloor=0\) for both mesh families. The energy-type errors converge with order \(k+1\), and the \(L^2\) error \(\|u-u_0\|\) shows higher-order convergence.

6.4 Three-dimensional biharmonic problem↩︎

We also test the biharmonic problem in three dimensions on the unit cube \(\Omega=(0,1)^3\), with exact solution \[u(x,y,z)=\sin^2(\pi x)\sin^2(\pi y)\sin^2(\pi z).\] The domain is discretized by tetrahedral meshes. Although tetrahedra are simplicial elements rather than general polyhedra, they are polytopes and provide a useful three-dimensional test of the dimension-independent construction. For a tetrahedron, \(\nu=4\), and hence \(\left\lfloor\frac{m}{\nu}\right\rfloor = \left\lfloor\frac{2}{4}\right\rfloor =0.\) Thus no bubble enrichment is needed in this test. Table 5 shows that the method retains the predicted convergence behavior in three dimensions: the energy-type errors converge with order \(k+1\), and the \(L^2\) error shows higher-order convergence.

References↩︎

---

1. The first author was partially supported by NSF grant DMS-2309785, the third author was supported by NSFC grant 12431014, and the fourth author was supported by NSFC grant 12401483.↩︎