# Hybrid high-order and weak Galerkin methods for the biharmonic problem

Zhaonan Dong  Inria, 2 rue Simone Iff, 75589 Paris, France, and CERMICS, Ecole
des Ponts, 77455 Marne-la-Vallée 2, France<EMAIL_ADDRESS>Alexandre
Ern  CERMICS, Ecole des Ponts, 77455 Marne-la-Vallée 2, France, and Inria, 2
rue Simone Iff, 75589 Paris, France<EMAIL_ADDRESS>
###### Abstract

We devise and analyze two hybrid high-order (HHO) methods for the numerical
approximation of the biharmonic problem. The methods support polyhedral
meshes, rely on the primal formulation of the problem, and deliver
$O(h^{k+1})$ $H^{2}$-error estimates when using polynomials of order $k\geq 0$
to approximate the normal derivative on the mesh (inter)faces. Both HHO
methods hinge on a stabilization in the spirit of Lehrenfeld–Schöberl for
second-order PDEs. The cell unknowns are polynomials of order $(k+2)$ that can
be eliminated locally by means of static condensation. The face unknowns
approximating the trace of the solution on the mesh (inter)faces are
polynomials of order $(k+1)$ in the first HHO method which is valid in
dimension two and uses an original stabilization involving the canonical
hybrid finite element, and they are of order $(k+2)$ for the second HHO method
which is valid in arbitrary dimension and uses only $L^{2}$-orthogonal
projections in the stabilization. A comparative discussion with the weak
Galerkin methods from the literature is provided, highlighting the close
connections and the improvements proposed herein. Additionally, we show how
the two HHO methods can be combined with a Nitsche-like boundary-penalty
technique to weakly enforce the boundary conditions. An originality in the
devised Nitsche’s technique is to avoid any penalty parameter that must be
large enough. Finally, numerical results showcase the efficiency of the
proposed methods, and indicate that the HHO methods can generally outperform
discontinuous Galerkin methods and even be competitive with $C^{0}$-interior
penalty methods on triangular meshes.

##  1 Introduction

Fourth-order PDEs are encountered in the modeling of various physical
phenomena, such as plate bending, thin-plate elasticity,
microelectromechanical systems, and the Cahn–Hilliard phase-field model. In
the present work, we are concerned with the following model problem:

$\displaystyle\Delta^{2}u$ $\displaystyle=f$ $\displaystyle\text{in
$\Omega$},$ (1) $\displaystyle u$ $\displaystyle=0$ $\displaystyle\text{on
$\partial\Omega$},$ $\displaystyle\partial_{n}u$ $\displaystyle=0$
$\displaystyle\text{on $\partial\Omega$},$

where $\Omega$ is an open, bounded, polytopal, Lipschitz set in
$\mathbb{R}^{d}$, $d\geq 2$, with boundary $\partial\Omega$, the load $f$ is
in $L^{2}(\Omega)$, and $\partial_{n}$ denotes the normal derivative on
$\partial\Omega$. Non-homogeneous boundary conditions and a boundary condition
on the second-normal derivative can be readily incorporated. Instead,
considering more singular loads is nontrivial for the present purpose. We also
emphasize that the present developments hinge on the weak formulation of (1)
involving the Hessian.

The main goal of this work is to devise and analyze a discretization method
for (1) offering two main features: (i) it supports polyhedral meshes (the
mesh cells can be polyhedra as such or have a simple shape but contain hanging
nodes); (ii) it hinges on the primal formulation of the problem, thereby
leading to a symmetric positive definite system matrix. There are already some
methods in the literature achieving these goals. These methods can be loosely
classified into three groups, depending on the dimension of the smallest
geometric object to which discrete unknowns are attached. This criterion is
relevant since it influences the stencil of the method, and it also influences
the level of conformity that can be achieved in the approximation of the
solution. The methods in the first group were developed in the practically
important case where $d=2$. They attach discrete unknowns to the mesh
vertices, edges, and cells and can achieve $C^{1}$-conformity. Salient
examples are the $C^{1}$-conforming virtual element methods (VEM) from [6, 13]
and the $C^{0}$-conforming VEM from [50]. Another example of method in this
group is the nonconforming VEM from [3, 51] where the approximation is,
however, (fully) nonconforming. The methods in the second group attach
discrete unknowns only to the mesh faces and cells for $d\geq 2$. They are
amenable to static condensation (meaning that the cell unknowns can be
eliminated locally leading to a global problem coupling only the face
unknowns), and they provide a nonconforming approximation to the solution. The
two salient examples are the weak Galerkin (WG) methods from [38, 49, 48] and
the hybrid high-order (HHO) method from [4]. Finally, the methods in the third
group attach discrete unknowns only to the mesh cells for $d\geq 2$ and belong
to the class of interior penalty discontinuous Galerkin (IPDG) methods. These
are also nonconforming methods, and they were developed for the model problem
(1) in [37, 42, 31]. We mention that on specific meshes composed of simplices
or cuboids, there are variants of the above methods achieving
$C^{0}$-conformity, such as the $C^{0}$-WG method from [40, 12] and the
$C^{0}$-IPDG from [25, 5]. Furthermore, important examples of nonconforming
finite elements on simplicial meshes are the Morley element [36, 47] and the
Hsieh–Clough–Tocher (HCT) element (see, e.g., [15, Chap. 6]).

In the present work, we focus on HHO methods. HHO methods were introduced in
[21] for locking-free linear elasticity and in [22] for linear diffusion. The
two ingredients of HHO methods are a local reconstruction operator and a local
stabilization operator in each mesh cell. For second-order PDEs, the aim of
the first operator is to reconstruct locally a gradient from the cell and the
face unknowns, and the aim of the second operator is to penalize in a least-
squares sense the difference between the trace of the cell unknown and the
face unknown on every mesh face. HHO methods have undergone a vigorous
development in the last few years; to cite a few examples, we mention
Navier–Stokes flows [23], elastoplastic problems [2], Tresca friction problems
[14], spectral problems [9], and magnetostatics [11]. HHO methods were
embedded into the broad framework of hybridizable dG (HDG) methods in [17] by
reformulating the HHO equations as local balance equations with equilibrated
numerical fluxes. Moreover, HHO methods are closely related to WG methods,
which were also embedded into the HDG framework in [16, Sec. 6.6]. The
reconstruction operator in the HHO method corresponds to the weak gradient in
WG methods. Hence, HHO and WG methods differ only in the choice of the
discrete unknowns and in the design of stabilization. Although the close
connections between HHO and WG methods should be mutually beneficial, these
connections are, in the authors’ opinion, not sufficiently explicit in the
literature, and the title of the present work is also meant to draw the
community’s attention on this opportunity.

The design of the stabilization turns out to be a key ingredient so that the
method leads to _optimal_ error estimates. By this, we mean, in the case of a
second-order elliptic PDE, that the method delivers an $O(h^{k+1})$
$H^{1}$-error estimate, where $k\geq 0$ is the degree of the face unknowns.
Notice that this criterion is consistent with the classical properties of
hybridized mixed finite element methods. The point we want to make here is
that optimality cannot be reached on general meshes if one uses plain least-
squares stabilization, i.e., a more subtle design of the stabilization is
required. If the cell unknowns are of degree $k$, optimality is achieved in
[21, 22] by means of a stabilization that uses the reconstruction operator
(this is the first occurrence of this idea in the broad framework of HDG
methods). Alternatively, if the cell unknowns are of degree $(k+1)$, one can
use the Lehrenfeld–Schöberl (LS) stabilization [35], as in [17] for HHO
methods and in [39] for WG methods. Although the LS stabilization does not use
the reconstruction operator, it is not a plain least-squares stabilization,
since an orthogonal projection is applied to the trace of the cell unknowns.
We mention that it is also possible to achieve optimality _without
stabilization_ for second-order PDEs if one uses Raviart–Thomas functions of
degree $k$ to reconstruct the gradient (see [1, 19]). However, optimality is
lost if one reconstructs the gradient in larger polynomial spaces (the
convergence rate is in general $O(h^{k})$), since the normal component of the
reconstructed gradient on the mesh faces is too rich to be captured by the
face unknowns. Another possibility is to enrich the space for the gradient
reconstruction by suitable bubble functions based on the notion of
$M$-decomposition devised for HDG methods [18].

unknowns | cell | face | grad | $k$ | ref.
---|---|---|---|---|---
WG | $k+2$ | $k+2$ | $[k+1]^{d}$ | $k\geq 0$ | [38]
| $k+2$ | $k+2$ | $k+1$ | $k\geq 0$ | [38]
| $k+2$ | $k+1$ | $k+1$ | $k\geq 0$ | [49]
| 1 | 1 | $[1]^{d}$ | $k=0$ | [48]
HHO | $k$ | $k$ | $[k]^{d}$ | $k\geq 1$ | [4]
HHO(A) | $k+2$ | $k+1$ | $k$ | $k\geq 0$ | present ($d=2$)
HHO(B) | $k+2$ | $k+2$ | $k$ | $k\geq 0$ | present ($d\geq 2$)
Table 1: Discrete unknowns in HHO and WG methods from the literature and the
present work. In the column ‘grad’, the notation $[\cdot]^{d}$ means that the
full gradient is approximated; otherwise, only the normal derivative is
approximated. For all the methods, the integer $k$ is fixed by the fact that
the method delivers an $O(h^{k+1})$ $H^{2}$-error estimate.

To discretize fourth-order PDEs, HHO and WG methods use cell unknowns that are
meant to approximate the solution in each mesh cell, face unknowns that are
meant to approximate its trace on each mesh (inter)face, and additional face
unknowns that are meant to approximate either its full gradient trace or only
its normal derivative on each mesh (inter)face. The HHO and WG methods from
the literature and the present HHO methods are described in Table 1 in terms
of their discrete unknowns. To put all the methods on the same ground and
allow for a fair comparison, the polynomial degree $k$ is such that all the
methods in the table deliver an $O(h^{k+1})$ $H^{2}$-error estimate.
Consistently with the terminology adopted above for second-order elliptic
PDEs, the method can be viewed as _optimal_ if the order of the face unknowns
approximating the trace of the gradient (or of the normal derivative) is of
degree $k$. As seen from Table 1, the WG methods from the literature do not
meet this criterion. For instance, the HHO method from [4] with $k=1$
converges with one order higher than the WG method from [48] while using the
same discrete unknowns. The lack of optimality is related to the use of a
plain least-squares stabilization. Instead, the HHO method from [4] and the
present HHO methods are optimal, and this is reflected by a more elaborate
design of the stabilization. Notice that for fourth-order PDEs, this also
means that the Hessian (and not only the Laplacian) has to be reconstructed
locally. In [4], the stabilization design follows the spirit of [21, 22] in
that it uses a Hessian-based deflection reconstruction operator. In the
present methods, the design is performed in the spirit of the LS
stabilization. Another difference with [4] is that the present methods only
introduce face unknowns approximating the normal derivative of the solution on
the mesh (inter)faces (and not the full gradient trace). As a result, and
despite the slight increase in the degree of its face unknowns approximating
the solution trace, the present HHO methods involve less globally coupled
unknowns than in [4]; see the discussion in Remarks 3.1 and 5.2. Moreover, we
allow here for $k\geq 0$, whereas [4] requires $k\geq 1$. We also mention that
the increase of cell unknowns compared to [4] has a moderate impact on
computational costs owing to static condensation. This slight overhead is
actually compensated by the simplification in the stabilization term (see
Section 7 for further discussion).

Let us briefly summarize the main novelties and results of the present work:
(i) Two novel and computationally effective HHO methods leading to optimal
$O(h^{k+1})$ $H^{2}$-error estimates with polynomials of order $k\geq 0$ to
approximate the normal derivative; (ii) An original design in 2D using, for
the first time in HHO methods, the canonical hybrid finite element in the
stabilization; (iii) The HHO methods do not feature stabilization parameters
that must be large enough (only positive), in contrast with dG and
$C^{0}$-IPDG methods. (iv) A numerical study showing the attractive
performances of the proposed methods, which in particular can outperform dG
methods (except for low polynomial orders and Voronoi-like meshes where the
number of faces is quite large) and even be competitive with $C^{0}$-IPDG and
HCT methods on simplicial meshes; (v) A variant of the HHO methods using a
Nitsche-type boundary-penalty technique to weakly enforce the boundary
conditions. We notice in particular that the development of Nitsche’s
boundary-penalty technique is instrumental to deal with domains with curved
boundary (in the wake of [8, 7] for elliptic interface problems) and to derive
a robust approximation method in the case of singularly perturbed regimes.
These results are explored in our recent work [24]. We also emphasize that our
Nitsche technique does not need the penalty parameter to be large enough. This
is the first time this property is met for fourth-order PDEs, and to this
purpose, we adapt ideas from [34, 7] derived for second-order PDEs.
Heuristically, the reason for circumventing the constraint on having a large
enough penalty parameter is that the reconstruction operator in HHO methods
avoids the need to introduce an additional consistency term as in the standard
Nitsche method.

As a final remark, we mention that our main error estimates are established
for an exact solution that belongs to the broken Sobolev space
$H^{k+3}(\mathcal{T}_{h})$ (where $\mathcal{T}_{h}$ denotes the underlying
mesh) and to the Sobolev space $H^{2+s}(\Omega)$ with $s>\frac{3}{2}$. This
latter assumption follows the rather classical paradigm in the analysis of
nonconforming methods and is invoked when bounding the consistency error. As
discussed in Remark 4.7, the regularity gap can be lowered to $s>1$ by
adapting the techniques developed in [29] and [28, Chap. 40&41] in the context
of second-order elliptic PDEs. We also notice that quasi-optimal error
estimates for general loads in $H^{-2}(\Omega)$ are derived in [45, 44] for
the Morley element and the $C^{0}$-IPDG method (see also [10] for further
results in the case of various lowest-order methods). The techniques in [45,
44] require to modify the right-hand side of the discrete problem by means of
bubble functions and a $C^{1}$-smoother. These ideas have been adapted to HHO
methods for second-order elliptic PDEs with loads in $H^{-1}(\Omega)$ in [30].
We expect that the extension to the biharmonic problem could follow a similar
path for $d=2$, whereas for $d=3$, one difficulty is related, irrespective of
the considered discretization method, to the lack of a well-established and
computable $C^{1}$-smoother of arbitrary order.

The rest of this work is organized as follows. We introduce some basic
notation, the mesh assumptions, and some analysis tools in Section 2. In
Section 3, we introduce the HHO method in the 2D setting employing the
canonical hybrid finite element to design the stabilization. In Section 4, we
present the stability and error analysis of the method introduced in Section
3. In Section 5, we present the second HHO method, this time valid in
arbitrary dimension, and we outline the main changes in the stability and
error analysis from Section 4. Our numerical results indicate that in two
dimensions, the first HHO method from Section 3 is more effective than the
second method from Section 5. In Section 6, we combine the above HHO methods
with Nitsche’s boundary-penalty technique. Finally, numerical results
showcasing the computational advantages of the proposed HHO methods are
presented in Section 7.

##  2 Model problem and discrete setting

In this section, we introduce some basic notation, the weak formulation of the
model problem, and the discrete setting to formulate and analyze the HHO
discretization.

###  2.1 Basic notation and weak formulation

We use standard notation for the Lebesgue and Sobolev spaces and, in
particular, for the fractional-order Sobolev spaces, we consider the
Sobolev–Slobodeckij seminorm based on the double integral. For an open,
bounded, Lipschitz set $S$ in $\mathbb{R}^{d}$, $d\in\\{1,2,3\\}$, we denote
by $(v,w)_{S}$ the $L^{2}(S)$-inner product, and we employ the same notation
when $v$ and $w$ are vector- or matrix-valued fields. We denote by $\nabla w$
the (weak) gradient of $w$ and by $\nabla^{2}w$ its (weak) Hessian. Let
${\bm{n}}$ be the unit outward normal vector on the boundary $\partial S$ of
$S$. Assuming that the functions $v$ and $w$ are smooth enough, we have the
following integration by parts formula:

$(\Delta^{2}v,w)_{S}=(\nabla^{2}v,\nabla^{2}w)_{S}+(\nabla\Delta
v,{\bm{n}}w)_{\partial S}-(\nabla^{2}v{\bm{n}},\nabla w)_{\partial S}.$ (2)

Whenever the context is unambiguous, we denote by $\partial_{n}$ the (scalar-
valued) normal derivative on $\partial S$ and by $\partial_{t}$ the
($\mathbb{R}^{d-1}$-valued) tangential derivative. We also denote by
$\partial_{nn}v$ the (scalar-valued) normal-normal second-order derivative and
by $\partial_{nt}v$ the ($\mathbb{R}^{d-1}$-valued) normal-tangential second-
order derivative. The integration by parts formula (2) can then be rewritten
as

$(\Delta^{2}v,w)_{S}=(\nabla^{2}v,\nabla^{2}w)_{S}+(\partial_{n}\Delta
v,w)_{\partial S}-(\partial_{nn}v,\partial_{n}w)_{\partial
S}-(\partial_{nt}v,\partial_{t}w)_{\partial S}.$ (3)

In what follows, the set $S$ is always a polytope so that its boundary can be
decomposed into a finite union of planar faces with disjoint interiors.
Expressions involving the tangential derivative on $\partial S$ are then
implicitly understood to be evaluated as a summation over the faces composing
$\partial S$.

Using the above integration by parts formula, the following weak formulation
of the model problem (1) is classically derived: Find $u\in H^{2}_{0}(\Omega)$
such that

$(\nabla^{2}u,\nabla^{2}v)_{\Omega}=(f,v)_{\Omega},\qquad\forall v\in
H^{2}_{0}(\Omega),$ (4)

The well-posedness of (4) is proven, e.g., in [32, Section 1.5].

######  Remark 2.1 (Non-homogeneous conditions)

Since the domain $\Omega$ is a polytope, its boundary can be split into
$\\{\partial\Omega_{i}\\}_{i=1}^{N}$ $(d-1)$-dimensional planar faces with
disjoint interiors. Let $g_{D}$ and $g_{N}$ be boundary data such that
$g_{D}|_{\partial\Omega_{i}}\in H^{\frac{3}{2}}(\partial\Omega_{i})$ and
$g_{N}|_{\partial\Omega_{i}}\in H^{\frac{1}{2}}(\partial\Omega_{i})$ for all
$i\in\\{1,\ldots,N\\}$, as well as $g_{D}\in C^{0}(\partial\Omega)$. Then, one
can enforce the non-homogeneous boundary conditions $u=g_{D}$ and
$\partial_{n}u=g_{N}$ on all the faces $\partial\Omega_{i}$; see [33, Sections
1.5 &1.6].

###  2.2 Inverse, trace, and Poincaré inequalities

Let $\mathcal{T}_{h}$ be a mesh covering $\Omega$ exactly. The mesh
$\mathcal{T}_{h}$ can have cells that are disjoint open polytopes in
$\mathbb{R}^{d}$ (with planar faces), and hanging nodes are possible. A
generic mesh cell is denoted by $K\in\mathcal{T}_{h}$, its diameter by
$h_{K}$, and its unit outward normal by ${\bm{n}}_{K}$. We assume that the
mesh belongs to a shape-regular mesh sequence $(\mathcal{T}_{h})_{h>0}$ in the
sense of [21]. In a nutshell, any mesh $\mathcal{T}_{h}$ admits a matching
simplicial submesh $\mathcal{T}_{h}^{\prime}$ such that any cell (or face) of
$\mathcal{T}_{h}^{\prime}$ is a subset of exactly one cell (or at most one
face) of $\mathcal{T}_{h}$. Moreover, there exists a mesh-regularity parameter
$\rho>0$ such that for all $h>0$, all $K\in\mathcal{T}_{h}$, and all
$S\in\mathcal{T}_{h}^{\prime}$ such that $S\subset K$, we have $\rho h_{S}\leq
r_{S}$ and $\rho h_{K}\leq h_{S}$, where $r_{S}$ denotes the inradius of the
simplex $S$. The mesh faces are collected in the set $\mathcal{F}_{h}$, which
is split as the set $\mathcal{F}_{h}^{\mathrm{b}}$ containing the mesh
boundary faces and the set $\mathcal{F}_{h}^{\mathrm{i}}$ containing the mesh
interfaces. In this work, we make the mild additional assumption that the mesh
faces are connected; the reason for this is that we will consider an
approximation operator on the mesh faces that is only $H^{1}$-stable, and not
$L^{2}$-stable, so that we will need to invoke some polynomial approximation
properties directly on the mesh faces (see (14)). Let ${\bm{n}}_{F}$ denote
the unit normal vector orienting any mesh face $F\in\mathcal{F}_{h}$. The
direction of ${\bm{n}}_{F}$ is arbitrary, but fixed once and for all, for all
$F\in\mathcal{F}_{h}^{\mathrm{i}}$, and ${\bm{n}}_{F}:={\bm{n}}$ for all
$F\in\mathcal{F}_{h}^{\mathrm{b}}$. For any mesh cell $K\in\mathcal{T}_{h}$,
the mesh faces composing its boundary $\partial K$ are collected in the set
$\mathcal{F}_{\partial K}$. The shape-regularity of the mesh sequence implies
that for all $K\in\mathcal{T}_{h}$ and all $F\in\mathcal{F}_{\partial K}$ with
diameter $h_{F}$, the length scales $h_{K}$ and $h_{F}$ are uniformly
equivalent and that $\\#(\mathcal{F}_{\partial K})$ is uniformly bounded.

Let us recall some important analysis tools. We refer the reader, e.g., to
[20, Sec. 1.4] for the proofs of Lemma 2.2 and Lemma 2.3 and to [46] for the
derivation of the Poincaré inequality in $H^{2}$ from the corresponding
inequality in $H^{1}$. For all $K\in\mathcal{T}_{h}$ and all $k\geq 0$,
$\mathbb{P}^{k}(K)$ denotes the linear space composed of the restriction to
$K$ of polynomials of total degree at most $k$.

######  Lemma 2.2 (Discrete inverse and trace inequalities)

Let $\mathcal{T}_{h}$ belong to a shape-regular mesh sequence and let $k\geq
0$. There are constants $C_{\mathrm{inv}}^{\mathrm{tr}}$ and
$C_{\mathrm{inv}}$, only depending on the mesh shape-regularity, the
polynomial degree $k$, and the space dimension $d$, such that for all
$v_{h}\in\mathbb{P}^{k}(K)$ and all $K\in\mathcal{T}_{h}$,

$\displaystyle\|{v}_{h}\|_{\partial K}$ $\displaystyle\leq
C_{\mathrm{inv}}^{\mathrm{tr}}h_{K}^{-\frac{1}{2}}\|{v}_{h}\|_{K},$ (5)
$\displaystyle\|\nabla{v}_{h}\|_{K}$ $\displaystyle\leq
C_{\mathrm{inv}}h_{K}^{-1}\|{v}_{h}\|_{K},$ (6)

######  Lemma 2.3 (Multiplicative trace inequality)

Let $\mathcal{T}_{h}$ belong to a shape-regular mesh sequence. There is a
constant $C_{\mathrm{mt}}$, only depending on the mesh shape-regularity and
the space dimension $d$, such that for all $v\in H^{1}(K)$ and all
$K\in\mathcal{T}_{h}$,

$\|{v}\|_{\partial K}\leq
C_{\mathrm{mt}}\big{(}h_{K}^{-\frac{1}{2}}\|{v}\|_{K}+\|{v}\|_{K}^{\frac{1}{2}}\|\nabla{v}\|_{K}^{\frac{1}{2}}\big{)}.$
(7)

######  Lemma 2.4 (Poincaré inequality)

Let $\mathcal{T}_{h}$ belong to a shape-regular mesh sequence. There is a
constant $C_{\mathrm{P}}$, only depending on the mesh shape-regularity and the
space dimension $d$, such that for all $v\in H^{2}(K)^{\perp}:=\\{v\in
H^{2}(K)\;|\;(v,\xi)_{K}=0,\;\forall\xi\in\mathbb{P}^{1}(K)\\}$ and all
$K\in\mathcal{T}_{h}$,

$h_{K}^{-2}\|{v}\|_{K}+h_{K}^{-1}\|\nabla{v}\|_{K}\leq
C_{\mathrm{P}}\|\nabla^{2}{v}\|_{K}.$ (8)

######  Remark 2.5 (Discrete inverse inequality on faces)

Similarly to (6) and recalling that the diameter of any face
$F\in\mathcal{F}_{\partial K}$ is uniformly equivalent to $h_{K}$, one can
prove that there is a constant $\widetilde{C}_{\mathrm{inv}}$, only depending
on the mesh shape-regularity, the polynomial degree $k$, and the space
dimension $d$, such that

$\|\partial_{t}{v}_{h}\|_{F}\leq\widetilde{C}_{\mathrm{inv}}h_{K}^{-1}\|{v}_{h}\|_{F},$
(9)

for all $v_{h}\in\mathbb{P}^{k}(K)$, all $K\in\mathcal{T}_{h}$, and all
$F\in\mathcal{F}_{\partial K}$.

######  Remark 2.6 (Fractional multiplicative trace inequality)

Let $\mathcal{T}_{h}$ belong to a shape-regular mesh sequence. Let
$s\in(\frac{1}{2},1]$. There is a constant $C_{\mathrm{fmt}}$, only depending
on the mesh shape-regularity and the space dimension $d$, such that for all
$v\in H^{s}(K)$ and all $K\in\mathcal{T}_{h}$,

$\|{v}\|_{\partial K}\leq
C_{\mathrm{fmt}}\big{(}h_{K}^{-\frac{1}{2}}\|{v}\|_{K}+h_{K}^{s-\frac{1}{2}}|v|_{H^{s}(K)}\big{)}.$
(10)

The proof when $K$ is a simplex can be found in [26, Lem. 7.2]. In the general
case, for every subface of a face in $\mathcal{F}_{\partial K}$, one carves a
subsimplex inside $K$ whose height is uniformly equivalent to $h_{K}$. Notice
that for $s=1$, (10) is a simple consequence of (7) and Young’s inequality
since $|v|_{H^{1}(K)}=\|\nabla v\|_{K}$.

###  2.3 Polynomial approximation in cells and on faces

Let $k\geq 0$ and let $\Pi^{k+2}_{K}$ be the $L^{2}$-orthogonal projection
onto $\mathbb{P}^{k+2}(K)$. Since the mesh cells can be decomposed into a
finite number of subsimplices, the approximation properties of $\Pi^{k+2}_{K}$
can be established by proceeding as in [26, Lem. 5.4].

######  Lemma 2.7 (Polynomial approximation in $K$)

Let $\mathcal{T}_{h}$ belong to a shape-regular mesh sequence. Let $k\geq 0$.
There is a constant $C_{\mathrm{app}}$, only depending on the mesh shape-
regularity, the polynomial degree $k$, and the space dimension $d$, such that
for all $t\in[0,k+3]$, all $m\in\\{0,\ldots,\lfloor t\rfloor\\}$, all $v\in
H^{t}(K)$, and all $K\in\mathcal{T}_{h}$,

$|v-\Pi^{k+2}_{K}(v)|_{H^{m}(K)}\leq
C_{\mathrm{app}}h^{t-m}_{K}|{v}|_{H^{t}(K)}.$ (11)

Another useful property of $\Pi^{k+2}_{K}$ results from the multiplicative
trace inequality (7) and the Poincaré inequality (8). Indeed, we infer that
there is a constant $C_{\Pi}$, only depending on the mesh shape-regularity,
the polynomial degree $k$, and the space dimension $d$, such that for all
$v\in H^{2}(K)$ and all $K\in\mathcal{T}_{h}$,

$h_{K}^{-\frac{3}{2}}\|v-\Pi^{k+2}_{K}(v)\|_{\partial
K}+h_{K}^{-\frac{1}{2}}\|\nabla(v-\Pi^{k+2}_{K}(v))\|_{\partial K}\leq
C_{\Pi}\|\nabla^{2}(v-\Pi^{k+2}_{K}(v))\|_{K}.$ (12)

We will use two operators for the polynomial approximation on the mesh faces.
The first one is an $L^{2}$-orthogonal projection. Specifically, letting
$\mathbb{P}^{k}(\mathcal{F}_{\partial K}):=\times_{F\in\mathcal{F}_{\partial
K}}\mathbb{P}^{k}(F)$ for all $k\geq 0$ and all $K\in\mathcal{T}_{h}$, we
denote by $\Pi^{k}_{\partial K}$ the $L^{2}$-orthogonal projection onto
$\mathbb{P}^{k}(\mathcal{F}_{\partial K})$. Notice that $\Pi^{k}_{\partial
K}(v)$ can be computed independently for each face $F\in\mathcal{F}_{\partial
K}$. The second operator is specific to the 2D setting where the mesh faces
are straight segments. On the reference interval $\widehat{I}:=(-1,1)$, the
canonical hybrid finite element of degree $(k+1)$ has for its degrees of
freedom the value at the two endpoints and, for $k\geq 1$, the integrals on
$\widehat{I}$ weighted by a chosen set of basis functions in
$\mathbb{P}^{k-1}(\widehat{I})$ (see, e.g., [27, Sec. 6.3.3 & 7.6] or [41,
Thm. 3.14]). For all $F\in\mathcal{F}_{h}$, let
$J_{F}^{k+1}:H^{1}(F)\rightarrow\mathbb{P}^{k+1}(F)$ be the corresponding
interpolation operator generated using geometric affine mappings. Then, the
two key identities satisfied by $J_{F}^{k+1}$ are for all $v\in H^{1}(F)$,

$\big{(}\partial_{t}(v-J_{F}^{k+1}(v)),\xi\big{)}_{F}=0,\;\forall\xi\in\mathbb{P}^{k}(F),\qquad(v-J_{F}^{k+1}(v),\theta)_{F}=0,\;\forall\theta\in\mathbb{P}^{k-1}(F),$
(13)

or, in more compact form,
$\Pi_{F}^{k}(\partial_{t}v)=\partial_{t}(J_{F}^{k+1}(v))$ and
$\Pi^{k-1}_{F}\circ J_{F}^{k+1}=\Pi^{k-1}_{F}$. Moreover, $J_{F}^{k+1}$
satisfies the following approximation properties: There is $C_{J}$, only
depending on the mesh shape-regularity, the polynomial degree $k$, and the
space dimension $d$, such that

$\|v-J_{F}^{k+1}(v)\|_{F}\leq
C_{J}h_{F}\|\partial_{t}v\|_{F},\qquad\|v-J_{F}^{k+1}(v)\|_{F}\leq
C_{J}h_{F}^{2}\|\partial_{tt}v\|_{F},$ (14)

for all $v\in H^{1}(F)$ and all $v\in H^{2}(F)$, respectively. Notice that for
$k=0$, $J_{F}^{1}$ coincides with the Lagrange interpolate on $F$ based on its
two endpoints.

In what follows, it is convenient to rewrite (13) and (14) on the whole
boundary of every mesh cell $K\in\mathcal{T}_{h}$. Letting
$H^{l}(\mathcal{F}_{\partial K}):=\\{v\in L^{2}(\partial K)\;|\;v|_{F}\in
H^{l}(F),\;\forall F\in\mathcal{F}_{\partial K}\\}$ with $l\in\\{1,2\\}$,
$J_{\partial K}^{k+1}:H^{1}(\mathcal{F}_{\partial
K})\rightarrow\mathbb{P}^{k+1}(\mathcal{F}_{\partial K})$ is defined facewise
by setting $J_{\partial K}^{k+1}(v)|_{F}:=J_{F}^{k+1}(v|_{F})$ for all $v\in
H^{1}(\mathcal{F}_{\partial K})$. Recalling that the tangential derivative is
understood to act facewise, we obtain for all $v\in
H^{1}(\mathcal{F}_{\partial K})$,

$\big{(}\partial_{t}(v-J_{\partial K}^{k+1}(v)),\xi\big{)}_{\partial
K}=0,\;\forall\xi\in\mathbb{P}^{k}(\mathcal{F}_{\partial
K}),\qquad(v-J_{\partial K}^{k+1}(v),\theta)_{\partial
K}=0,\;\forall\theta\in\mathbb{P}^{k-1}(\mathcal{F}_{\partial K}).$ (15)

Moreover, there is $\widetilde{C}_{J}$ having the same dependencies as $C_{J}$
such that

$\|v-J_{\partial K}^{k+1}(v)\|_{\partial
K}\leq\widetilde{C}_{J}h_{K}\|\partial_{t}v\|_{\partial
K},\qquad\|v-J_{\partial K}^{k+1}(v)\|_{\partial
K}\leq\widetilde{C}_{J}h_{K}^{2}\|\partial_{tt}v\|_{\partial K},$ (16)

for all $v\in H^{1}(\mathcal{F}_{\partial K})$ and all $v\in
H^{2}(\mathcal{F}_{\partial K})$, respectively, where we used that $h_{F}$ and
$h_{K}$ are uniformly equivalent for all $F\in\mathcal{F}_{\partial K}$.

##  3 HHO method for the 2D biharmonic problem

Let $k\geq 0$ be the polynomial degree. For all $K\in\mathcal{T}_{h}$, the
local HHO space is

$\widehat{V}_{K}^{k}:=\mathbb{P}^{k+2}(K)\times\mathbb{P}^{k+1}(\mathcal{F}_{\partial
K})\times\mathbb{P}^{k}(\mathcal{F}_{\partial K}).$ (17)

A generic element in $\widehat{V}_{K}^{k}$ is denoted
$\widehat{v}_{K}:=(v_{K},v_{\partial K},\gamma_{\partial K})$ with
$v_{K}\in\mathbb{P}^{k+2}(K)$, $v_{\partial
K}\in\mathbb{P}^{k+1}(\mathcal{F}_{\partial K})$, and $\gamma_{\partial
K}\in\mathbb{P}^{k}(\mathcal{F}_{\partial K})$. The first component of
$\widehat{v}_{K}$ aims at representing the solution inside the mesh cell, the
second its trace on the cell boundary, and the third its normal derivative on
the cell boundary (along the direction of the outward normal ${\bm{n}}_{K}$).
In what follows, it is implicitly understood that within integrals over
$\partial K$, the symbol $\partial_{n}$ means ${\bm{n}}_{K}{\cdot}\nabla$.

###  3.1 Reconstruction and stabilization

The HHO method is formulated locally by means of a reconstruction and a
stabilization operator. The local reconstruction operator
$R_{K}:\widehat{V}_{K}^{k}\rightarrow\mathbb{P}^{k+2}(K)$ is such that, for
all $\widehat{v}_{K}:=(v_{K},v_{\partial K},\gamma_{\partial
K})\in\widehat{V}_{K}^{k}$, $R_{K}(\widehat{v}_{K})\in\mathbb{P}^{k+2}(K)$ is
determined by solving the following well-posed problem:

$\displaystyle(\nabla^{2}R_{K}(\widehat{v}_{K}),\nabla^{2}w)_{K}={}$
$\displaystyle(\nabla^{2}{v}_{K},\nabla^{2}w)_{K}+(v_{K}-v_{\partial
K},\partial_{n}\Delta w)_{\partial K}-(\partial_{n}v_{K}-\gamma_{\partial
K},\partial_{nn}w)_{\partial K}$ (18)
$\displaystyle-(\partial_{t}(v_{K}-v_{\partial K}),\partial_{nt}w)_{\partial
K},\quad\forall w\in\mathbb{P}^{k+2}(K),$
$\displaystyle(R_{K}(\widehat{v}_{K}),\xi)_{K}={}$
$\displaystyle(v_{K},\xi)_{K},\quad\forall\xi\in\mathbb{P}^{1}(K).$

When computing $R_{K}(\widehat{v}_{K})$, one actually takes
$w\in\mathbb{P}^{k+2}(K)^{\perp}:=\\{w\in\mathbb{P}^{k+2}(K)\;|\;(w,\xi)_{K}=0,\;\forall\xi\in\mathbb{P}^{1}(K)\\}$
since the equation is trivial whenever $w\in\mathbb{P}^{1}(K)$. Moreover,
owing to the integration by parts formula (3), we infer that

$(\nabla^{2}R_{K}(\widehat{v}_{K}),\nabla^{2}w)_{K}=(v_{K},\Delta^{2}w)_{K}-(v_{\partial
K},\partial_{n}\Delta w)_{\partial K}+(\gamma_{\partial
K},\partial_{nn}w)_{\partial K}+(\partial_{t}v_{\partial
K},\partial_{nt}w)_{\partial K}.$ (19)

This expression shows that in the rightmost term on the right-hand side, we
take advantage of the face component $v_{\partial K}$ to represent the
tangential derivative at the boundary of $K$. Notice also that for $k=0$, the
second term on the right-hand side vanishes.

The local stabilization bilinear form $S_{\partial K}$ is defined such that,
for all
$(\widehat{v}_{K},\widehat{w}_{K})\in\widehat{V}_{K}^{k}\times\widehat{V}_{K}^{k}$
with $\widehat{v}_{K}:=(v_{K},v_{\partial K},\gamma_{\partial K})$ and
$\widehat{w}_{K}:=(w_{K},w_{\partial K},\chi_{\partial K})$,

$\displaystyle S_{\partial K}(\widehat{v}_{K},\widehat{w}_{K}):={}$
$\displaystyle h_{K}^{-3}\big{(}J^{k+1}_{\partial K}(v_{\partial
K}-v_{K}),J^{k+1}_{\partial K}(w_{\partial K}-{w}_{K})\big{)}_{\partial K}$
(20) $\displaystyle+h_{K}^{-1}\big{(}\Pi^{k}_{\partial K}(\gamma_{\partial
K}-\partial_{n}{v}_{K}),\Pi^{k}_{\partial K}(\chi_{\partial
K}-\partial_{n}{w}_{K})\big{)}_{\partial K}.$

Notice the use of the interpolation operator $J^{k+1}_{\partial K}$ for the
first term on the right-hand side. The reconstruction and stabilization
operators are combined together to build the local bilinear form $a_{K}$ on
$\widehat{V}_{K}^{k}\times\widehat{V}_{K}^{k}$ such that

$a_{K}(\widehat{v}_{K},\widehat{w}_{K}):=(\nabla^{2}R_{K}(\widehat{v}_{K}),\nabla^{2}R_{K}(\widehat{w}_{K}))_{K}+S_{\partial
K}(\widehat{v}_{K},\widehat{w}_{K}).$ (21)

###  3.2 The global discrete problem

We define the global HHO space as

$\widehat{V}_{h}^{k}:=\mathbb{P}^{k+2}(\mathcal{T}_{h})\times\mathbb{P}^{k+1}(\mathcal{F}_{h})\times\mathbb{P}^{k}(\mathcal{F}_{h}).$
(22)

A generic element in $\widehat{V}_{h}^{k}$ is denoted
$\widehat{v}_{h}:=(v_{\mathcal{T}_{h}},v_{\mathcal{F}_{h}},\gamma_{\mathcal{F}_{h}})$
with $v_{\mathcal{T}_{h}}:=(v_{K})_{K\in\mathcal{T}_{h}}$,
$v_{\mathcal{F}_{h}}:=(v_{F})_{F\in\mathcal{F}_{h}}$, and
$\gamma_{\mathcal{F}_{h}}:=(\gamma_{F})_{F\in\mathcal{F}_{h}}$, where
$\gamma_{F}$ is meant to approximate the normal derivative in the direction of
the unit normal vector ${\bm{n}}_{F}$ orienting $F$. For all
$K\in\mathcal{T}_{h}$, the local components of $\widehat{v}_{h}$ are collected
in the triple $\widehat{v}_{K}:=(v_{K},v_{\partial K},\gamma_{\partial
K})\in\widehat{V}_{K}^{k}$ with $v_{\partial K}|_{F}:=v_{F}$ and
$\gamma_{\partial K}|_{F}:=({\bm{n}}_{F}{\cdot}{\bm{n}}_{K})\gamma_{F}$ for
all $F\in\mathcal{F}_{\partial K}$. Notice that the way the face components of
$v_{\partial K}$ are assigned follows the usual way of HHO methods for second-
order elliptic PDEs, whereas the definition of the face components of
$\gamma_{\partial K}$ takes into account the orientation of the faces in
$\mathcal{F}_{\partial K}$. Furthermore, we enforce the homogeneous boundary
conditions strongly by considering the subspace

$\widehat{V}_{h0}^{k}:=\\{\widehat{v}_{h}\in\widehat{V}_{h}^{k}\;|\;v_{F}=\gamma_{F}=0,\;\forall
F\in\mathcal{F}_{h}^{\mathrm{b}}\\}.$ (23)

The discrete HHO problem for the 2D biharmonic problem is as follows: Find
$\widehat{u}_{h}\in\widehat{V}_{h0}^{k}$ such that

$a_{h}(\widehat{u}_{h},\widehat{w}_{h})=\ell(w_{\mathcal{T}_{h}}),\qquad\forall\widehat{w}_{h}\in\widehat{V}_{h0}^{k},$
(24)

where the global discrete bilinear form $a_{h}$ and the global linear form
$\ell$ are assembled cellwise as

$a_{h}(\widehat{v}_{h},\widehat{w}_{h}):=\sum_{K\in\mathcal{T}_{h}}a_{K}(\widehat{v}_{K},\widehat{w}_{K}),\qquad\ell(w_{\mathcal{T}_{h}}):=(f,w_{\mathcal{T}_{h}})_{\Omega}=\sum_{K\in\mathcal{T}_{h}}(f,{w}_{K})_{K}.$
(25)

Notice that only the first component of the triple $\widehat{w}_{h}$ is used
to evaluate the right-hand side in (24). An important observation is that the
discrete problem (24) is amenable to static condensation. Indeed, the cell
unknowns can be eliminated locally in every mesh cell, leading to a global
problem where the only remaining unknowns are those attached to the mesh
faces, i.e., those in
$\mathbb{P}^{k+1}(\mathcal{F}_{h})\times\mathbb{P}^{k}(\mathcal{F}_{h})$.

######  Remark 3.1 (Comparison with [4])

The present HHO method is cheaper than the one from [4] where the globally
coupled unknowns are in
$\mathbb{P}^{k}(\mathcal{F}_{h})\times\mathbb{P}^{k}(\mathcal{F}_{h};\mathbb{R}^{2})$
with $k\geq 1$. Indeed, in the present method, there are $(2k+3)$, $k\geq 0$,
unknowns per mesh interface, whereas this number is $(3k+3)$, $k\geq 1$, in
[4]. On the other hand, the present method is slightly more expensive
regarding static condensation since the number of cell unknowns is
$\frac{1}{2}(k+3)(k+4)$ vs. $\frac{1}{2}(k+2)(k+3)$ in [4]. However, the
slight overhead incurred in the static condensation is compensated by the
simpler form of the stabilization; see Section 7 for more insight into the
computational costs.

######  Remark 3.2 (Variant)

It is also possible to consider the slightly cheaper choice
$\widehat{V}_{K}^{k}:=\mathbb{P}^{k+1}(K)\times\mathbb{P}^{k+1}({\partial
K})\times\mathbb{P}^{k}({\partial K})$ with $k\geq 0$. With this choice, the
number of cell unknowns to be statically condensed is slightly reduced, but
the size of the global problem coupling all the face unknowns is unchanged.
Notice that this choice requires to modify the stabilization bilinear form by
setting

$\displaystyle S_{\partial K}(\widehat{v}_{K},\widehat{w}_{K}):={}$
$\displaystyle
h_{K}^{-4}\big{(}\Pi^{k+1}_{K}(v_{K}-R_{K}(\widehat{v}_{K})),\Pi^{k+1}_{K}(w_{K}-R_{K}(\widehat{w}_{K}))\big{)}_{K}$
$\displaystyle+h_{K}^{-3}\big{(}J^{k+1}_{\partial K}(v_{\partial
K}-R_{K}(\widehat{v}_{K})),J^{k+1}_{\partial K}(w_{\partial
K}-R_{K}(\widehat{w}_{K}))\big{)}_{\partial K}$
$\displaystyle+h_{K}^{-1}\big{(}\Pi^{k}_{\partial K}(\gamma_{\partial
K}-\partial_{n}R_{K}(\widehat{v}_{K})),\Pi^{k}_{\partial K}(\chi_{\partial
K}-\partial_{n}R_{K}(\widehat{w}_{K}))\big{)}_{\partial K}.$ (26)

The analysis of this variant will not be detailed herein, but this variant
will be included in the numerical investigations presented in Section 7.

######  Remark 3.3 (Boundary conditions)

In the non-homogeneous case, the HHO solution is sought in the space
$\widehat{V}_{h}^{k}$, whereas the test functions are still taken in the space
$\widehat{V}_{h0}^{k}$. The value of the components of the HHO solution
attached to the mesh boundary faces is then assigned by means of the
projections $J^{k+1}_{\partial K}({g_{D}}|_{\partial K})$ and
$\Pi^{k}_{\partial K}({g_{N}}|_{\partial K})$. The convergence analysis
proceeds as for homogeneous boundary conditions, up to straightforward
adaptations when bounding the consistency error. Details are omitted for
brevity. It is also possible to consider the boundary conditions
$u=\partial_{nn}u=0$ on $\partial\Omega$. The discrete HHO bilinear form is
still defined as above, but the discrete problem (24) now involves the
subspace
$\widehat{V}_{h0}^{k}:=\\{\widehat{v}_{h}\in\widehat{V}_{h}^{k}\;|\;v_{F}=0,\;\forall
F\in\mathcal{F}_{h}^{\mathrm{b}}\\}$, i.e., the boundary condition $u=0$ is
still strongly enforced, whereas the boundary condition $\partial_{nn}u=0$ is
weakly enforced.

##  4 Stability and error analysis

In this section, we perform the stability and error analysis of the HHO method
devised in the previous section for the 2D biharmonic problem. We first
establish a local stability property for the bilinear form $a_{K}$ together
with the well-posedness of the discrete problem (24). Then, we introduce a
suitable reduction operator leading to optimal approximation properties, we
bound the corresponding consistency error, and finally we derive the error
estimate.

In what follows, the symbol $C$ denotes a generic positive constant whose
value can change at each occurrence, provided this value only depends on the
mesh shape-regularity, the polynomial degree $k$, and the space dimension $d$.

###  4.1 Stability and well-posedness

We equip the local HHO space $\widehat{V}_{K}^{k}$ with the $H^{2}$-like
seminorm such that for all $\widehat{v}_{K}:=(v_{K},v_{\partial
K},\gamma_{\partial K})\in\widehat{V}_{K}^{k}$,

$|\widehat{v}_{K}|^{2}_{\widehat{V}_{K}^{k}}:=\|\nabla^{2}v_{K}\|_{K}^{2}+h_{K}^{-3}\|v_{\partial
K}-v_{K}\|_{\partial K}^{2}+h_{K}^{-1}\|\gamma_{\partial
K}-\partial_{n}v_{K}\|_{\partial K}^{2}.$ (27)

######  Lemma 4.1 (Local stability and boundedness)

There is a real number $\alpha>0$, depending only on the mesh shape-
regularity, the polynomial degree $k$, and the space dimension $d$, such that
for all $\widehat{v}_{K}\in\widehat{V}_{K}^{k}$ and all $K\in\mathcal{T}_{h}$,

$\alpha|\widehat{v}_{K}|^{2}_{\widehat{V}_{K}^{k}}\leq\|\nabla^{2}R_{K}(\widehat{v}_{K})\|_{K}^{2}+S_{\partial
K}(\widehat{v}_{K},\widehat{v}_{K})\leq\alpha^{-1}|\widehat{v}_{K}|^{2}_{\widehat{V}_{K}^{k}}.$
(28)

###### Proof.

(1) Lower bound. Using the reconstruction defined in (18) with
$w:=v_{K}\in\mathbb{P}^{k+2}(K)$, we have

$\displaystyle(\nabla^{2}R_{K}(\widehat{v}_{K}),\nabla^{2}v_{K})_{K}={}$
$\displaystyle\|\nabla^{2}{v}_{K}\|^{2}_{K}+(v_{K}-v_{\partial
K},\partial_{n}\Delta v_{K})_{\partial K}$ (29)
$\displaystyle-(\partial_{n}v_{K}-\gamma_{\partial
K},\partial_{nn}v_{K})_{\partial K}-(\partial_{t}(v_{K}-v_{\partial
K}),\partial_{nt}v_{K})_{\partial K}$ $\displaystyle={}$
$\displaystyle\|\nabla^{2}{v}_{K}\|^{2}_{K}+(J_{\partial
K}^{k+1}(v_{K}-v_{\partial K}),\partial_{n}\Delta v_{K})_{\partial K}$
$\displaystyle-(\Pi_{\partial K}^{k}(\partial_{n}v_{K}-\gamma_{\partial
K}),\partial_{nn}v_{K})_{\partial K}-(\partial_{t}J_{\partial
K}^{k+1}(v_{K}-v_{\partial K}),\partial_{nt}v_{K})_{\partial K},$

where we used the two identities from (15) together with $\partial_{n}\Delta
v_{K}\in\mathbb{P}^{k-1}(\mathcal{F}_{\partial K})$ (if $k\geq 1$, otherwise
this term vanishes) and
$\partial_{nn}v_{K},\partial_{nt}v_{K}\in\mathbb{P}^{k}(\mathcal{F}_{\partial
K})$. Using the Cauchy–Schwarz inequality together with the inverse
inequalities (6), (5), and (9), we infer that

$\displaystyle\|\nabla^{2}{v}_{K}\|^{2}_{K}\leq{}$
$\displaystyle\|\nabla^{2}R_{K}(\widehat{v}_{K})\|_{K}\|\nabla^{2}{v}_{K}\|_{K}+\|J_{\partial
K}^{k+1}(v_{K}-v_{\partial K})\|_{\partial K}\|\partial_{n}\Delta
v_{K}\|_{\partial K}$ $\displaystyle+\|\Pi_{\partial
K}^{k}(\partial_{n}v_{K}-\gamma_{\partial K})\|_{\partial
K}\|\partial_{nn}v_{K}\|_{\partial K}+\|\partial_{t}J_{\partial
K}^{k+1}(v_{K}-v_{\partial K})\|_{\partial K}\|\partial_{nt}v_{K}\|_{\partial
K}$ $\displaystyle\leq{}$
$\displaystyle\|\nabla^{2}R_{K}(\widehat{v}_{K})\|_{K}\|\nabla^{2}{v}_{K}\|_{K}+CS_{\partial
K}(\widehat{v}_{K},\widehat{v}_{K})^{\frac{1}{2}}\|\nabla^{2}{v}_{K}\|_{K},$

which shows that

$\|\nabla^{2}{v}_{K}\|_{K}\leq\|\nabla^{2}R_{K}(\widehat{v}_{K})\|_{K}+CS_{\partial
K}(\widehat{v}_{K},\widehat{v}_{K})^{\frac{1}{2}}.$ (30)

Moreover, since $v_{\partial K}-v_{K}=J_{\partial K}^{k+1}(v_{\partial
K}-v_{K})-(v_{K}-J_{\partial K}^{k+1}(v_{K}))$, the triangle inequality
implies that

$\displaystyle h_{K}^{-\frac{3}{2}}\|v_{\partial K}-v_{K}\|_{\partial
K}\leq{}$ $\displaystyle h_{K}^{-\frac{3}{2}}\|J_{\partial
K}^{k+1}(v_{\partial K}-v_{K})\|_{\partial
K}+h_{K}^{-\frac{3}{2}}\|v_{K}-J_{\partial K}^{k+1}(v_{K})\|_{\partial K}$
$\displaystyle\leq{}$ $\displaystyle S_{\partial
K}(\widehat{v}_{K},\widehat{v}_{K})^{\frac{1}{2}}+Ch_{K}^{\frac{1}{2}}\|\partial_{tt}v_{K}\|_{\partial
K}$ $\displaystyle\leq{}$ $\displaystyle S_{\partial
K}(\widehat{v}_{K},\widehat{v}_{K})^{\frac{1}{2}}+C\|\nabla^{2}v_{K}\|_{K},$

where we used the approximation property (16) and the discrete trace
inequality (5). Combining this bound with (30) proves that

$h_{K}^{-\frac{3}{2}}\|v_{\partial K}-v_{K}\|_{\partial K}\leq
C\big{(}\|\nabla^{2}R_{K}(\widehat{v}_{K})\|_{K}+S_{\partial
K}(\widehat{v}_{K},\widehat{v}_{K})^{\frac{1}{2}}\big{)}.$ (31)

Proceeding similarly shows that

$\displaystyle h_{K}^{-\frac{1}{2}}\|\gamma_{\partial
K}-\partial_{n}v_{K}\|_{\partial K}\leq{}$ $\displaystyle
h_{K}^{-\frac{1}{2}}\|\Pi_{\partial K}^{k}(\gamma_{\partial
K}-\partial_{n}v_{K})\|_{\partial
K}+h_{K}^{-\frac{1}{2}}\|\partial_{n}v_{K}-\Pi_{\partial
K}^{k}(\partial_{n}v_{K})\|_{\partial K}$ $\displaystyle\leq{}$ $\displaystyle
S_{\partial
K}(\widehat{v}_{K},\widehat{v}_{K})^{\frac{1}{2}}+Ch_{K}^{\frac{1}{2}}\|\partial_{nt}v_{K}\|_{\partial
K}$ $\displaystyle\leq{}$ $\displaystyle S_{\partial
K}(\widehat{v}_{K},\widehat{v}_{K})^{\frac{1}{2}}+C\|\nabla^{2}v_{K}\|_{K},$

where we used the approximation properties of $\Pi_{\partial K}^{k}$ and the
discrete trace inequality (5). Combining this bound with (30) proves that

$h_{K}^{-\frac{1}{2}}\|\gamma_{\partial K}-\partial_{n}v_{K}\|_{\partial
K}\leq C\big{(}\|\nabla^{2}R_{K}(\widehat{v}_{K})\|_{K}+S_{\partial
K}(\widehat{v}_{K},\widehat{v}_{K})^{\frac{1}{2}}\big{)}.$ (32)

Finally, combining the bounds (30), (31), and (32) proves the lower bound in
(28).
(2) Upper bound. Using this time
$w:=R_{K}(\widehat{v}_{K})\in\mathbb{P}^{k+2}(K)$ in the reconstruction
defined in (18) and proceeding as above shows that

$\|\nabla^{2}R_{K}(\widehat{v}_{K})\|_{K}\leq\|\nabla^{2}{v}_{K}\|_{K}+CS_{\partial
K}(\widehat{v}_{K},\widehat{v}_{K})^{\frac{1}{2}}.$ (33)

Moreover, still proceeding as above, we infer that

$\displaystyle h_{K}^{-\frac{3}{2}}\|J_{\partial K}^{k+1}(v_{\partial
K}-v_{K})\|_{\partial K}\leq{}$ $\displaystyle
h_{K}^{-\frac{3}{2}}\|v_{\partial K}-v_{K}\|_{\partial
K}+h_{K}^{-\frac{3}{2}}\|v_{K}-J_{\partial K}^{k+1}(v_{K})\|_{\partial K}$
$\displaystyle\leq{}$ $\displaystyle h_{K}^{-\frac{3}{2}}\|v_{\partial
K}-v_{K}\|_{\partial K}+C\|\nabla^{2}v_{K}\|_{K},$

and

$h_{K}^{-\frac{1}{2}}\|\Pi_{\partial K}^{k}(\gamma_{\partial
K}-\partial_{n}v_{K})\|_{\partial K}\leq
h_{K}^{-\frac{1}{2}}\|\gamma_{\partial K}-\partial_{n}v_{K}\|_{\partial
K}+C\|\nabla^{2}v_{K}\|_{K}.$

Putting the above two bounds together shows that

$S_{\partial K}(\widehat{v}_{K},\widehat{v}_{K})^{\frac{1}{2}}\leq
C|\widehat{v}_{K}|_{\widehat{V}_{K}^{k}}.$ (34)

Finally, the combination of (33) and (34) proves the upper bound in (28). ∎

We equip the space $\widehat{V}_{h0}^{k}$ with the norm

$\|\widehat{v}_{h}\|_{\widehat{V}_{h0}^{k}}:=\sum_{K\in\mathcal{T}_{h}}|\widehat{v}_{K}|^{2}_{\widehat{V}_{K}^{k}},\qquad\forall\widehat{v}_{h}\in\widehat{V}_{h0}^{k}.$
(35)

To show that this indeed defines a norm, consider
$\widehat{v}_{h}\in\widehat{V}_{h0}^{k}$ such that
$\|\widehat{v}_{h}\|_{\widehat{V}_{h0}^{k}}=0$. Then, for all
$K\in\mathcal{T}_{h}$, $v_{K}\in\mathbb{P}^{1}(K)$ and $v_{\partial K}=v_{K}$
and $\gamma_{\partial K}=\partial_{n}v_{K}$ on $\partial K$. For any cell
$K\in\mathcal{T}_{h}$ having at least one boundary face, say
$F\in\mathcal{F}_{\partial K}\cap\mathcal{F}_{h}^{\mathrm{b}}$, we have
$v_{F}=\gamma_{F}=0$ by definition of $\widehat{V}_{h0}^{k}$. Since $v_{K}$ is
affine and its gradient vanishes identically on $F$, $\nabla v_{K}$ vanishes
in $K$, and since $v_{K}$ vanishes on $F$, we infer that $v_{K}$ vanishes
identically in $K$. This implies that $v_{F}=\gamma_{F}=0$ for all
$F\in\mathcal{F}_{\partial K}$. We can then propagate the reasoning one layer
of cells further inside the domain, and by repeating the process, we reach all
the cells composing the mesh. Thus, the three components of the triple
$\widehat{v}_{h}$ vanish identically everywhere.

######  Corollary 4.2 (Coercivity and well-posedness)

The discrete bilinear form $a_{h}$ is coercive on $\widehat{V}_{h0}^{k}$, and
the discrete problem (24) is well-posed.

###### Proof.

Summing the lower bound in (28) over all the mesh cells shows the following
coercivity property:

$a_{h}(\widehat{v}_{h},\widehat{v}_{h})\geq\alpha\|\widehat{v}_{h}\|_{\widehat{V}_{h0}^{k}}^{2},\qquad\forall\widehat{v}_{h}\in\widehat{V}_{h0}^{k}.$
(36)

The well-posedness of (24) then follows from the Lax–Milgram lemma. ∎

###  4.2 Local reduction operator and polynomial approximation

For all $K\in\mathcal{T}_{h}$, we define the local reduction operator
$\mathcal{\widehat{I}}^{k}_{K}:H^{2}(K)\rightarrow\widehat{V}_{K}^{k}$ such
that for all $v\in H^{2}(K)$,

$\mathcal{\widehat{I}}^{k}_{K}(v):=\big{(}\Pi_{K}^{k+2}(v),J_{\partial
K}^{k+1}(v),\Pi_{\partial K}^{k}({\bm{n}}_{K}{\cdot}\nabla
v)\big{)}\in\widehat{V}_{K}^{k}.$ (37)

Moreover, the $H^{2}$-elliptic projection
$\mathcal{E}_{K}(v):H^{2}(K)\rightarrow\mathbb{P}^{k+2}(K)$ is defined such
that

$\displaystyle(\nabla^{2}(\mathcal{E}_{K}(v)-v),\nabla^{2}w)_{K}$
$\displaystyle=0,$ $\displaystyle\forall w\in\mathbb{P}^{k+2}(K),$ (38)
$\displaystyle(\mathcal{E}_{K}(v)-v,\xi)_{K}$ $\displaystyle=0,$
$\displaystyle\forall\xi\in\mathbb{P}^{1}(K).$

The following lemma states the two key properties of the local reduction
operator defined in (37).

######  Lemma 4.3 (Local reduction operator)

We have $R_{K}\circ\mathcal{\widehat{I}}^{k}_{K}=\mathcal{E}_{K}$ for all
$K\in\mathcal{T}_{h}$. Moreover, for all $K\in\mathcal{T}_{h}$ and all $v\in
H^{2}(K)$, we have

$S_{\partial
K}(\mathcal{\widehat{I}}^{k}_{K}(v),\mathcal{\widehat{I}}^{k}_{K}(v))^{\frac{1}{2}}\leq
C\|\nabla^{2}(v-\Pi_{K}^{k+2}(v))\|_{K}.$ (39)

###### Proof.

Let $K\in\mathcal{T}_{h}$ and let $v\in H^{2}(K)$.
(1) Using the definition (19) of the reconstruction operator, we infer that
for all $w\in\mathbb{P}^{k+2}(K)$, we have

$\displaystyle(\nabla^{2}R_{K}(\mathcal{\widehat{I}}^{k}_{K}(v)),\nabla^{2}w)_{K}={}$
$\displaystyle(\Pi_{K}^{k+2}(v),\Delta^{2}w)_{K}$ $\displaystyle-(J_{\partial
K}^{k+1}(v),\partial_{n}\Delta w)_{\partial K}+(\Pi_{\partial
K}^{k}(\partial_{n}v),\partial_{nn}w)_{\partial K}+(\partial_{t}(J_{\partial
K}^{k+1}(v)),\partial_{nt}w)_{\partial K}.$

Since $\Delta^{2}w\in\mathbb{P}^{k-2}(K)$ for $k\geq 2$ (and vanishes
otherwise), $\partial_{n}\Delta w\in\mathbb{P}^{k-1}(K)$ for $k\geq 1$ (and
vanishes otherwise), $\partial_{nn}w\in\mathbb{P}^{k}(K)$, and
$\partial_{nt}w\in\mathbb{P}^{k}(K)$, the $L^{2}$-orthogonality properties of
$\Pi_{K}^{k+2}$ and $\Pi_{\partial K}^{k}$, together with the identities (15)
satisfied by $J_{\partial K}^{k+1}$ imply that

$\displaystyle(\nabla^{2}R_{K}(\mathcal{\widehat{I}}^{k}_{K}(v)),\nabla^{2}w)_{K}=(v,\Delta^{2}w)_{K}-(v,\partial_{n}\Delta
w)_{\partial K}+(\partial_{n}v,\partial_{nn}w)_{\partial
K}+(\partial_{t}v,\partial_{nt}w)_{\partial K}=(\nabla^{2}v,\nabla^{2}w)_{K}.$

Moreover, for all $\xi\in\mathbb{P}^{1}(K)$, we have
$(R_{K}(\mathcal{\widehat{I}}^{k}_{K}(v)),\xi)_{K}=(\Pi_{K}^{k+2}(v),\xi)_{K}=(v,\xi)_{K}$
for all $\xi\in\mathbb{P}_{1}(K)$. The above two identities prove that
$R_{K}(\mathcal{\widehat{I}}^{k}_{K}(v))=\mathcal{E}_{K}(v)$ for all $v\in
H^{2}(K)$. Thus, $R_{K}\circ\mathcal{\widehat{I}}^{k}_{K}=\mathcal{E}_{K}$.
(2) Let us now prove (39). Recalling the definitions (20) and (37), we have

$S_{\partial
K}(\mathcal{\widehat{I}}^{k}_{K}(v),\mathcal{\widehat{I}}^{k}_{K}(v))^{\frac{1}{2}}\leq
h_{K}^{-\frac{3}{2}}\|J^{k+1}_{\partial K}(v-\Pi^{k+2}_{K}(v))\|_{\partial
K}+h_{K}^{-\frac{1}{2}}\|\Pi^{k}_{\partial
K}(\partial_{n}v-\partial_{n}\Pi^{k+2}_{K}(v))\|_{\partial K},$ (40)

where we used that $J^{k+1}_{\partial K}\circ J^{k+1}_{\partial
K}=J^{k+1}_{\partial K}$ and $\Pi^{k}_{\partial K}\circ\Pi^{k}_{\partial
K}=\Pi^{k}_{\partial K}$. We start with the first term in (40) where we set
$\phi:=v-\Pi^{k+2}_{K}(v)$. Notice that $\phi\in H^{2}(K)^{\perp}$ and that
$\phi|_{\partial K}\in H^{1}(\mathcal{F}_{\partial K})$. Invoking the triangle
inequality, the approximation property (16), and the trace inequality (12)
shows that

$\displaystyle h_{K}^{-\frac{3}{2}}\|J^{k+1}_{\partial K}(\phi)\|_{\partial
K}$ $\displaystyle\leq h_{K}^{-\frac{3}{2}}\|\phi\|_{\partial
K}+h_{K}^{-\frac{3}{2}}\|\phi-J^{k+1}_{\partial K}(\phi)\|_{\partial K}$
$\displaystyle\leq h_{K}^{-\frac{3}{2}}\|\phi\|_{\partial
K}+Ch_{K}^{-\frac{1}{2}}\|\partial_{t}\phi\|_{\partial K}\leq
C\|\nabla^{2}\phi\|_{K}=C\|\nabla^{2}(v-\Pi^{k+2}_{K}(v))\|_{K}.$

Moreover, for the second term in (40), we invoke the $L^{2}(\partial
K)$-stability of $\Pi_{\partial K}^{k}$ and the trace inequality (12) to show
that

$\displaystyle h_{K}^{-\frac{1}{2}}\|\Pi^{k}_{\partial
K}(\partial_{n}v-\partial_{n}(\Pi^{k+2}_{K}(v))\|_{\partial K}$
$\displaystyle\leq
h_{K}^{-\frac{1}{2}}\|\partial_{n}(v-\Pi^{k+2}_{K}(v))\|_{\partial K}\leq
C\|\nabla^{2}(v-\Pi^{k+2}_{K}(v))\|_{K}.$

Combining the above two bounds with (40) proves (39). ∎

To bound the consistency error in the next section, we will consider a norm
that is stronger than the $H^{2}$-norm. For all $K\in\mathcal{T}_{h}$ and all
$v\in H^{2+s}(K)$, $s>\frac{3}{2}$, we consider the following norm:

$\displaystyle\|v\|^{2}_{\sharp,K}:=\|\nabla^{2}v\|_{K}^{2}+h^{3}_{K}\|\partial_{n}\Delta
v\|_{\partial K}^{2}+h_{K}\|\partial_{nn}v\|_{\partial
K}^{2}+h_{K}\|\partial_{nt}v\|_{\partial K}^{2}.$ (41)

######  Lemma 4.4 (Approximation)

The following holds true for all $K\in\mathcal{T}_{h}$ and all $v\in
H^{2+s}(K)$, $s>\frac{3}{2}$:

$\|v-\mathcal{E}_{K}(v)\|_{\sharp,K}\leq C\|v-\Pi_{K}^{k+2}(v)\|_{\sharp,K}.$
(42)

###### Proof.

Using the triangle inequality, we have

$\|v-\mathcal{E}_{K}(v)\|_{\sharp,K}\leq\|v-\Pi_{K}^{k+2}(v)\|_{\sharp,K}+\|\mathcal{E}_{K}(v)-\Pi_{K}^{k+2}(v)\|_{\sharp,K},$

so that we only need to bound the second term on the right-hand side. Owing to
the discrete inverse and trace inequalities (5) and (6), we readily infer that

$\|\mathcal{E}_{K}(v)-\Pi_{K}^{k+2}(v)\|_{\sharp,K}\leq
C\|\nabla^{2}(\mathcal{E}_{K}(v)-\Pi_{K}^{k+2}(v))\|_{K},$

so that it remains to bound
$\|\nabla^{2}(\mathcal{E}_{K}(v)-\Pi_{K}^{k+2}(v))\|_{K}$. Recalling that
$\mathcal{E}_{K}=R_{K}\circ\mathcal{\widehat{I}}^{k}_{K}$, using the
definition (18) of the reconstruction operator, and reasoning as in the proof
of Lemma 4.3 to remove the various projection operators, we infer that

$\displaystyle(\nabla^{2}\mathcal{E}_{K}(v),\nabla^{2}w)_{K}={}$
$\displaystyle(\nabla^{2}\Pi_{K}^{k+2}(v),\nabla^{2}w)_{K}+(\Pi_{K}^{k+2}(v)-v,\partial_{n}\Delta
w)_{\partial K}$
$\displaystyle-(\partial_{n}(\Pi_{K}^{k+2}(v)-v),\partial_{nn}w)_{\partial
K}-(\partial_{t}(\Pi_{K}^{k+2}(v)-v),\partial_{nt}w)_{\partial K},$

for all $w\in\mathbb{P}^{k+2}(K)$. Taking
$w:=\mathcal{E}_{K}(v)-\Pi_{K}^{k+2}(v)$, and invoking the Cauchy–Schwarz
inequality together with the discrete trace and inverse inequalities (5) and
(6), we infer that

$\displaystyle\|\nabla^{2}(\mathcal{E}_{K}(v)-\Pi_{K}^{k+2}(v))\|_{K}$
$\displaystyle\leq
C\big{(}h^{-\frac{3}{2}}_{K}\|v-\Pi_{K}^{k+2}(v)\|_{\partial
K}+h^{-\frac{1}{2}}_{K}\|\nabla(v-\Pi_{K}^{k+2}(v))\|_{\partial K}\big{)}$
$\displaystyle\leq C\|\nabla^{2}(v-\Pi_{K}^{k+2}(v))\|_{K},$

where the last bound follows from the trace inequality (12). This completes
the proof. ∎

###  4.3 Bound on consistency error

The global reduction operator
$\mathcal{\widehat{I}}_{h}^{k}:H^{2}(\Omega)\rightarrow\widehat{V}_{h}^{k}$ is
defined such that for all $v\in H^{2}(\Omega)$,

$\mathcal{\widehat{I}}_{h}^{k}(v):=\big{(}(\Pi_{K}^{k+2}(v))_{K\in\mathcal{T}_{h}},(J^{k+1}_{F}(v))_{F\in\mathcal{F}_{h}},(\Pi_{F}^{k}({\bm{n}}_{F}{\cdot}\nabla
v))_{F\in\mathcal{F}_{h}}\big{)}\in\widehat{V}_{h}^{k},$ (43)

recalling that $v$ and $\nabla v$ are single-valued on every
$F\in\mathcal{F}_{h}^{\mathrm{i}}$ for all $v\in H^{2}(\Omega)$. Importantly,
we notice that for all $K\in\mathcal{T}_{h}$, the local components of
$\mathcal{\widehat{I}}_{h}^{k}(v)$ attached to $K$ and the faces composing its
boundary are $\mathcal{\widehat{I}}^{k}_{K}(v|_{K})$. Moreover, for the exact
solution $u$ of (4), we have
$\mathcal{\widehat{I}}_{h}^{k}(u)\in\widehat{V}_{h0}^{k}$. We define the
consistency error $\delta_{h}\in(\widehat{V}_{h0}^{k})^{\prime}$ such that

$\langle\delta_{h},\widehat{w}_{h}\rangle:=\ell(w_{\mathcal{T}_{h}})-a_{h}(\mathcal{\widehat{I}}_{h}^{k}(u),\widehat{w}_{h}),\qquad\forall\widehat{w}_{h}\in\widehat{V}_{h0}^{k},$
(44)

where $\langle\cdot,\cdot\rangle$ denotes the duality pairing between
$(\widehat{V}_{h0}^{k})^{\prime}$ and $\widehat{V}_{h0}^{k}$.

######  Lemma 4.5 (Consistency)

Assume that $u\in H^{2+s}(\Omega)$ with $s>\frac{3}{2}$. The following holds
true:

$\|\delta_{h}\|_{(\widehat{V}_{h0}^{k})^{\prime}}:=\sup_{\widehat{w}_{h}\in\widehat{V}_{h0}^{k}}\frac{|\langle\delta_{h},\widehat{w}_{h}\rangle|}{\|\widehat{w}_{h}\|_{\widehat{V}_{h0}^{k}}}\leq
C\left(\sum_{K\in\mathcal{T}_{h}}\|u-\Pi_{K}^{k+2}(u)\|^{2}_{\sharp,K}\right)^{\frac{1}{2}}.$
(45)

###### Proof.

Let $\widehat{w}_{h}\in\widehat{V}_{h0}^{k}$. Using the definition of $\ell$
in (25), the PDE and the boundary conditions satisfied by the exact solution
$u$, and integrating by parts cellwise, we infer that

$\ell(w_{\mathcal{T}_{h}})=\sum_{K\in\mathcal{T}_{h}}\Big{\\{}(\nabla^{2}u,\nabla^{2}{w}_{K})_{K}+(\partial_{n}\Delta
u,w_{K})_{\partial K}-(\partial_{nn}u,\partial_{n}w_{K})_{\partial
K}-(\partial_{nt}u,\partial_{t}w_{K})_{\partial K}\Big{\\}}.$

The assumption $u\in H^{2+s}(\Omega)$ with $s>\frac{3}{2}$ implies that
$(\partial_{n}\Delta u)|_{\partial K}$, $(\partial_{nn}u)|_{\partial K}$, and
$(\partial_{nt}u)|_{\partial K}$ are meaningful in $L^{2}(\partial K)$ and
single-valued at every mesh interface. Moreover, since $w_{\partial K}$,
$\partial_{t}w_{\partial K}$, and $\chi_{\partial K}$ are single-valued at
every mesh interface $\mathcal{F}_{h}^{\mathrm{i}}$ and vanish at each mesh
boundary face since $\widehat{w}_{h}\in\widehat{V}_{h0}^{k}$, we have

$\displaystyle\ell(w_{\mathcal{T}_{h}})={}$
$\displaystyle\sum_{K\in\mathcal{T}_{h}}\Big{\\{}(\nabla^{2}u,\nabla^{2}{w}_{K})_{K}+(\partial_{n}\Delta
u,w_{K}-w_{\partial K})_{\partial K}$
$\displaystyle-(\partial_{nn}u,\partial_{n}w_{K}-\chi_{\partial K})_{\partial
K}-(\partial_{nt}u,\partial_{t}(w_{K}-w_{\partial K}))_{\partial K}\Big{\\}}.$

Since $a_{h}$ is assembled cellwise (see (25)) and the local components of
$\mathcal{\widehat{I}}_{h}^{k}(u)$ are $\mathcal{\widehat{I}}^{k}_{K}(u|_{K})$
for all $K\in\mathcal{T}_{h}$, we infer that
$a_{h}(\mathcal{\widehat{I}}_{h}^{k}(u),\widehat{w}_{h})=\sum_{K\in\mathcal{T}_{h}}a_{K}(\mathcal{\widehat{I}}^{k}_{K}(u|_{K}),\widehat{w}_{K})$.
Using the definition (21) of $a_{K}$, the definition (18) of
$R_{K}(\widehat{w}_{K})$, and the identity
$R_{K}\circ\mathcal{\widehat{I}}^{k}_{K}=\mathcal{E}_{K}$ from Lemma 4.3 leads
to

$\displaystyle a_{h}(\mathcal{\widehat{I}}_{h}^{k}(u),\widehat{w}_{h})={}$
$\displaystyle\sum_{K\in\mathcal{T}_{h}}\Big{\\{}(\nabla^{2}\mathcal{E}_{K}(u),\nabla^{2}{w}_{K})_{K}+(\partial_{n}\Delta\mathcal{E}_{K}(u),w_{K}-w_{\partial
K})_{\partial K}$
$\displaystyle-(\partial_{nn}\mathcal{E}_{K}(u),\partial_{n}w_{K}-\chi_{\partial
K})_{\partial
K}-(\partial_{nt}\mathcal{E}_{K}(u),\partial_{t}(w_{K}-w_{\partial
K}))_{\partial K}+S_{\partial
K}(\mathcal{\widehat{I}}^{k}_{K}(u),\widehat{w}_{K})\Big{\\}}.$

Defining the function $\eta$ cellwise as
$\eta|_{K}:=u|_{K}-\mathcal{E}_{K}(u|_{K})$ for all $K\in\mathcal{T}_{h}$, we
infer that

$\displaystyle\langle\delta_{h},\widehat{w}_{h}\rangle={}$
$\displaystyle\sum_{K\in\mathcal{T}_{h}}\Big{\\{}(\nabla^{2}\eta,\nabla^{2}{w}_{K})_{K}+(\partial_{n}\Delta\eta,w_{K}-w_{\partial
K})_{\partial K}$ (46)
$\displaystyle-(\partial_{nn}\eta,\partial_{n}w_{K}-\chi_{\partial
K})_{\partial K}-(\partial_{nt}\eta,\partial_{t}(w_{K}-w_{\partial
K}))_{\partial K}-S_{\partial
K}(\mathcal{\widehat{I}}^{k}_{K}(u),\widehat{w}_{K})\Big{\\}}.$

(Notice that $(\nabla^{2}\eta,\nabla^{2}{w}_{K})_{K}=0$, but we keep this term
since it can be bounded as the other ones.) Let us denote by
$\mathcal{T}_{1,K}$ the first four addends on the right-hand side and by
$\mathcal{T}_{2,K}$ the fifth addend. We bound $\mathcal{T}_{1,K}$ by the
Cauchy–Schwarz inequality and also invoke the inverse inequality (9).
Recalling the definition (41) of the $\|{\cdot}\|_{\sharp,K}$-norm, this
yields

$|\mathcal{T}_{1,K}|\leq
C\|\eta\|_{\sharp,K}|\widehat{w}_{K}|_{\widehat{V}_{K}^{k}}.$

Moreover, owing to (39) and the upper bound in (28), we have

$|\mathcal{T}_{2,K}|\leq S_{\partial
K}(\mathcal{\widehat{I}}^{k}_{K}(u),\mathcal{\widehat{I}}^{k}_{K}(u))^{\frac{1}{2}}S_{\partial
K}(\widehat{w}_{K},\widehat{w}_{K})^{\frac{1}{2}}\leq
C\|\nabla^{2}(u-\Pi_{K}^{k+2}(u))\|_{K}|\widehat{w}_{K}|_{\widehat{V}_{K}^{k}}.$

Altogether, this implies that

$|\langle\delta_{h},\widehat{w}_{h}\rangle|\leq
C\left(\sum_{K\in\mathcal{T}_{h}}\|\eta\|^{2}_{\sharp,K}+\|\nabla^{2}(u-\Pi_{K}^{k+2}(u))\|_{K}^{2}\right)^{\frac{1}{2}}\|\widehat{w}_{h}\|_{\widehat{V}_{h0}^{k}}.$

Invoking Lemma 4.4, this completes the proof. ∎

###  4.4 Error estimate

We are now ready to establish the main result concerning the error analysis.

######  Theorem 4.6 ($H^{2}$-error estimate)

Assume that $u\in H^{2+s}(\Omega)$ with $s>\frac{3}{2}$. The following holds
true:

$\sum_{K\in\mathcal{T}_{h}}\|\nabla^{2}(u-R_{K}(\widehat{u}_{K}))\|_{K}^{2}\leq
C\sum_{K\in\mathcal{T}_{h}}\|u-\Pi_{K}^{k+2}(u)\|_{\sharp,K}^{2}.$ (47)

Consequently, if $k\geq 1$, assuming $u|_{K}\in H^{k+3}(K)$ for all
$K\in\mathcal{T}_{h}$, we have

$\sum_{K\in\mathcal{T}_{h}}\|\nabla^{2}(u-R_{K}(\widehat{u}_{K}))\|_{K}^{2}\leq
C\sum_{K\in\mathcal{T}_{h}}\big{(}h_{K}^{k+1}|u|_{H^{k+3}(K)}\big{)}^{2},$
(48)

and if $k=0$, letting $\sigma:=\min(s-1,1)\in(\frac{1}{2},1]$, we have

$\sum_{K\in\mathcal{T}_{h}}\|\nabla^{2}(u-R_{K}(\widehat{u}_{K}))\|_{K}^{2}\leq
C\sum_{K\in\mathcal{T}_{h}}\big{(}h_{K}(|u|_{H^{3}(K)}+h_{K}^{\sigma}|u|_{H^{3+\sigma}(K)})\big{)}^{2}.$
(49)

###### Proof.

Set
$\widehat{e}_{h}:=\mathcal{\widehat{I}}_{h}^{k}(u)-\widehat{u}_{h}\in\widehat{V}_{h0}^{k}$,
so that
$a_{h}(\widehat{e}_{h},\widehat{e}_{h})=-\langle\delta_{h},\widehat{e}_{h}\rangle$.
The coercivity property (36) implies that

$\alpha\|\widehat{e}_{h}\|_{\widehat{V}_{h0}^{k}}^{2}\leq
a_{h}(\widehat{e}_{h},\widehat{e}_{h})=-\langle\delta_{h},\widehat{e}_{h}\rangle\leq\|\delta_{h}\|_{(\widehat{V}_{h0}^{k})^{\prime}}\|\widehat{e}_{h}\|_{\widehat{V}_{h0}^{k}},$

so that
$\|\widehat{e}_{h}\|_{\widehat{V}_{h0}^{k}}\leq\frac{1}{\alpha}\|\delta_{h}\|_{(\widehat{V}_{h0}^{k})^{\prime}}$.
Since
$\sum_{K\in\mathcal{T}_{h}}\|\nabla^{2}R_{K}(\widehat{e}_{K})\|_{K}^{2}\leq
C\|\widehat{e}_{h}\|_{\widehat{V}_{h0}^{k}}^{2}$, we infer from Lemma 4.5 that

$\sum_{K\in\mathcal{T}_{h}}\|\nabla^{2}R_{K}(\widehat{e}_{K})\|_{K}^{2}\leq
C\sum_{K\in\mathcal{T}_{h}}\|u-\Pi_{K}^{k+2}(u)\|^{2}_{\sharp,K}.$

Since
$u-R_{K}(\widehat{u}_{K})=(u-\mathcal{E}_{K}(u))+R_{K}(\widehat{e}_{K})$, the
triangle inequality combined with Lemma 4.4 and the above bound proves (47).
Furthermore, (48) results from (47) and the approximation properties of
$\Pi_{K}^{k+2}$ (using Lemma 2.7 and the multiplicative trace inequality (7)).
Finally, (49) is proved similarly to (48), but this time invoking the
fractional multiplicative trace inequality (10), in particular to bound
$h_{K}^{\frac{3}{2}}\|\partial_{n}\Delta(u-\Pi_{K}^{2}(u))\|_{\partial
K}=h_{K}^{\frac{3}{2}}\|\partial_{n}\Delta u\|_{\partial K}$. ∎

######  Remark 4.7 (Regularity gap)

The error estimates in Theorem 4.6 require $u\in H^{2+s}(\Omega)$ with
$s>\frac{3}{2}$. This global regularity requirement on the exact solution can
be lowered to $s>1$ by using the techniques developed in [29] and [28, Chap.
40&41] in the context of second-order elliptic PDEs. Indeed, the crucial point
is to give a meaning to $\partial_{n}\Delta u$ on each mesh face, and this can
be done by applying the tools from [29, 28] to the field $\nabla\Delta u$.
Notice that the requirement $u\in H^{2+s}(\Omega)$ with $s>\frac{3}{2}$ is,
however, less stringent than the one resulting from achieving optimal decay
rates as soon as $k\geq 1$ (see (48)).

##  5 HHO method in arbitrary dimension

In this section, we adapt the material from the above two sections to devise
and analyze an HHO method to approximate the biharmonic problem in arbitrary
dimension $d\geq 2$. The main difference with the previous section is that the
interpolation operator $J_{\partial K}^{k+1}$ is no longer available if $d\geq
3$. The idea in this section is to raise the degree of the face unknowns
representing the solution trace to $(k+2)$, and to consider $L^{2}$-orthogonal
projections to lead the analysis. Thus, letting $k\geq 0$ be the polynomial
degree, the local HHO space considered in this section is such that for all
$K\in\mathcal{T}_{h}$,

$\widehat{V}_{K}^{k}:=\mathbb{P}^{k+2}(K)\times\mathbb{P}^{k+2}(\mathcal{F}_{\partial
K})\times\mathbb{P}^{k}(\mathcal{F}_{\partial K}).$ (50)

######  Remark 5.1 ($d=3$)

In 3D, on tetrahedral meshes, one can also generalize the HHO method from the
previous section by considering the canonical hybrid finite element of degree
$(k+2)$ on the mesh faces.

###  5.1 Reconstruction, stabilization, discrete problem, and stability

The local reconstruction operator is still defined by (18) (or, equivalently,
(19)). Instead, the local stabilization bilinear form $S_{\partial K}$ has to
be slightly modified and is now such that for all
$(\widehat{v}_{K},\widehat{w}_{K})\in\widehat{V}_{K}^{k}\times\widehat{V}_{K}^{k}$,

$S_{\partial K}(\widehat{v}_{K},\widehat{w}_{K}):=h_{K}^{-3}\big{(}v_{\partial
K}-v_{K},w_{\partial K}-{w}_{K}\big{)}_{\partial
K}+h_{K}^{-1}\big{(}\Pi^{k}_{\partial K}(\gamma_{\partial
K}-\partial_{n}{v}_{K}),\Pi^{k}_{\partial K}(\chi_{\partial
K}-\partial_{n}{w}_{K})\big{)}_{\partial K}.$ (51)

Notice that only $L^{2}$-orthogonal projections are considered. The local
bilinear form $a_{K}$ is defined on
$\widehat{V}_{K}^{k}\times\widehat{V}_{K}^{k}$ as in (21).

The global HHO space is now defined as

$\widehat{V}_{h}^{k}:=\mathbb{P}^{k+2}(\mathcal{T}_{h})\times\mathbb{P}^{k+2}(\mathcal{F}_{h})\times\mathbb{P}^{k}(\mathcal{F}_{h}).$
(52)

Focusing for simplicity on homogeneous boundary conditions, we consider the
subspace $\widehat{V}_{h0}^{k}$ obtained by zeroing out all the components
attached to the mesh boundary faces. The discrete HHO problem is as follows:
Find $\widehat{u}_{h}\in\widehat{V}_{h0}^{k}$ such that

$a_{h}(\widehat{u}_{h},\widehat{w}_{h})={\ell}(w_{\mathcal{T}_{h}}),\qquad\forall
w_{h}\in\widehat{V}_{h0}^{k},$ (53)

where $a_{h}$ and ${\ell}$ are still defined as in (25). Moreover, as in the
2D setting, the discrete problem (53) is amenable to static condensation,
whereby the cell unknowns are eliminated locally in every mesh cell, leading
to a global problem where the only remaining unknowns are those attached to
the mesh faces, i.e., those in
$\mathbb{P}^{k+2}(\mathcal{F}_{h})\times\mathbb{P}^{k}(\mathcal{F}_{h})$.

Finally, it is readily seen that the local stability and boundedness property
stated in Lemma 4.1 still holds true. Therefore, the discrete bilinear form
$a_{h}$ is coercive on $\widehat{V}_{h0}^{k}$, so that the discrete problem
(53) is well-posed owing to the Lax–Milgram lemma.

######  Remark 5.2 (Comparison with [4])

In the present HHO method, the global problem after static condensation
features $(2k+4)$, $k\geq 0$, unknowns per mesh interface, whereas this number
is $(4k+4)$, $k\geq 1$, for [4].

######  Remark 5.3 (Boundary conditions)

In the non-homogeneous case, similarly to Remark 3.3, the value of the
components of the HHO solution attached to the mesh boundary faces is assigned
by means of the projections $\Pi^{k+2}_{\partial K}(g_{D}|_{\partial K})$ and
$\Pi^{k}_{\partial K}(g_{N}|_{\partial K})$. It is also possible to enforce
the boundary conditions $u=\partial_{nn}u=0$ on $\partial\Omega$ by proceeding
as in Remark 3.3.

###  5.2 Polynomial approximation, consistency and error estimate

For all $K\in\mathcal{T}_{h}$, the local reduction operator
$\mathcal{\widehat{I}}^{k}_{K}:H^{2}(K)\rightarrow\widehat{V}_{K}^{k}$ is now
defined such that for all $v\in H^{2}(K)$,

$\mathcal{\widehat{I}}^{k}_{K}(v):=(\Pi_{K}^{k+2}(v),\Pi_{\partial
K}^{k+2}(v),\Pi_{\partial K}^{k}({\bm{n}}_{K}{\cdot}\nabla
v))\in\widehat{V}_{K}^{k}.$ (54)

We also define the operator
$\widetilde{\mathcal{E}}_{K}:=R_{K}\circ\mathcal{\widehat{I}}^{k}_{K}:H^{2}(K)\rightarrow\mathbb{P}^{k+2}(K)$.
Although this operator is no longer the $H^{2}$-elliptic projection, we can
show that it still enjoys the same approximation properties as those derived
in Lemma 4.4. Recall that the $\|{\cdot}\|_{\sharp,K}$-norm is defined in
(41).

######  Lemma 5.4 (Polynomial approximation)

The following holds true for all $K\in\mathcal{T}_{h}$ and all $v\in
H^{2+s}(K)$ with $s>\frac{3}{2}$:

$\|v-\widetilde{\mathcal{E}}_{K}(v)\|_{\sharp,K}\leq
C\|v-\Pi_{K}^{k+2}(v)\|_{\sharp,K}.$ (55)

Moreover, for all $K\in\mathcal{T}_{h}$ and all $v\in H^{2}(K)$, we have

$S_{\partial
K}(\mathcal{\widehat{I}}^{k}_{K}(v),\mathcal{\widehat{I}}^{k}_{K}(v))^{\frac{1}{2}}\leq
C\|v-\Pi_{K}^{k+2}(v)\|_{\sharp,K}.$ (56)

###### Proof.

(1) Using the definition (18) of $R_{K}$, the definition (54) of
$\mathcal{\widehat{I}}^{k}_{K}$, and the orthogonality property of the
$L^{2}$-projections $\Pi_{K}^{k+2}$ and $\Pi_{\partial K}^{k+2}$, we infer
that for all $w\in\mathbb{P}^{k+2}(K)$,

$\displaystyle(\nabla^{2}\widetilde{\mathcal{E}}_{K}(v),\nabla^{2}w)_{K}={}$
$\displaystyle(\nabla^{2}\Pi_{K}^{k+2}(v),\nabla^{2}w)_{K}-(\Pi_{\partial
K}^{k+2}(v)-\Pi_{K}^{k+2}(v),\partial_{n}\Delta w)_{\partial K}$
$\displaystyle+(\Pi_{\partial
K}^{k}(\partial_{n}v)-\partial_{n}\Pi_{K}^{k+2}(v),\partial_{nn}w)_{\partial
K}+(\partial_{t}(\Pi_{\partial
K}^{k+2}(v)-\Pi_{K}^{k+2}(v)),\partial_{nt}w)_{\partial K}$ $\displaystyle={}$
$\displaystyle(\nabla^{2}\Pi_{K}^{k+2}(v),\nabla^{2}w)_{K}-(v-\Pi_{K}^{k+2}(v),\partial_{n}\Delta
w)_{\partial K}$
$\displaystyle+(\partial_{n}(v-\Pi_{K}^{k+2}(v)),\partial_{nn}w)_{\partial
K}+(\partial_{t}(\Pi_{\partial
K}^{k+2}(v)-\Pi_{K}^{k+2}(v)),\partial_{nt}w)_{\partial K},$

Taking $w:=\widetilde{\mathcal{E}}_{K}(v)-\Pi_{K}^{k+2}(v)$, rearranging the
terms, and invoking the Cauchy–Schwarz inequality together with the inverse
inequalities (5), (6), (9) leads to

$\|\nabla^{2}(\widetilde{\mathcal{E}}_{K}(v)-\Pi_{K}^{k+2}(v))\|_{K}\\\ \leq
C\big{(}h_{K}^{-\frac{3}{2}}\|v-\Pi_{K}^{k+2}(v)\|_{\partial
K}+h_{K}^{-\frac{1}{2}}\|\partial_{n}(v-\Pi_{K}^{k+2}(v))\|_{\partial
K}+h_{K}^{-\frac{3}{2}}\|\Pi_{\partial
K}^{k+2}(v)-\Pi_{K}^{k+2}(v)\|_{\partial K}\big{)}.$

Concerning the rightmost term, we observe that $\Pi_{\partial
K}^{k+2}(v)-\Pi_{K}^{k+2}(v)=\Pi_{\partial K}^{k+2}(v-\Pi_{K}^{k+2}(v))$, so
that using the $L^{2}$-stability of $\Pi^{k+2}_{\partial K}$, we obtain

$\|\nabla^{2}(\widetilde{\mathcal{E}}_{K}(v)-\Pi_{K}^{k+2}(v))\|_{K}\leq
C\big{(}h_{K}^{-\frac{3}{2}}\|v-\Pi_{K}^{k+2}(v)\|_{\partial
K}+h_{K}^{-\frac{1}{2}}\|\partial_{n}(v-\Pi_{K}^{k+2}(v))\|_{\partial
K}\big{)}.$

The trace inequality (12) then shows that

$\|\nabla^{2}(\widetilde{\mathcal{E}}_{K}(v)-\Pi_{K}^{k+2}(v))\|_{K}\leq
C\|\nabla^{2}(v-\Pi_{K}^{k+2}(v))\|_{K}.$

The proof of (55) can now be completed by invoking the triangle inequality.
(2) Let us now prove (56). We have

$S_{\partial
K}(\mathcal{\widehat{I}}^{k}_{K}(v),\mathcal{\widehat{I}}^{k}_{K}(v))^{\frac{1}{2}}\leq
h_{K}^{-\frac{3}{2}}\|\Pi^{k+2}_{\partial K}(v)-\Pi^{k+2}_{K}(v)\|_{\partial
K}+h_{K}^{-\frac{1}{2}}\|\Pi^{k}_{\partial
K}(\partial_{n}v-\partial_{n}\Pi^{k+2}_{K}(v))\|_{\partial K},$

where we used that $\Pi^{k}_{\partial K}\circ\Pi^{k}_{\partial
K}=\Pi^{k}_{\partial K}$ for the second term on the right-hand side. To bound
the first term on the right-hand side, we invoke the same arguments as in the
first step of this proof leading to

$h_{K}^{-\frac{3}{2}}\|\Pi^{k+2}_{\partial K}(v)-\Pi^{k+2}_{K}(v)\|_{\partial
K}\leq C\|\nabla^{2}(v-\Pi_{K}^{k+2}(v))\|_{K}.$

Furthermore, the second term has already been bounded in the proof of Lemma
4.3. This completes the proof. ∎

The global reduction operator
$\mathcal{\widehat{I}}_{h}^{k}:H^{2}(\Omega)\rightarrow\widehat{V}_{h}^{k}$ is
defined such that for all $v\in H^{2}(\Omega)$,

$\mathcal{\widehat{I}}_{h}^{k}(v):=\big{(}(\Pi_{K}^{k+2}(v))_{K\in\mathcal{T}_{h}},(\Pi^{k+2}_{F}(v))_{F\in\mathcal{F}_{h}},(\Pi_{F}^{k}({\bm{n}}_{F}{\cdot}\nabla
v))_{F\in\mathcal{F}_{h}}\big{)}\in\widehat{V}_{h}^{k},$ (57)

so that the local components of $\mathcal{\widehat{I}}_{h}^{k}(v)$ are
$\mathcal{\widehat{I}}^{k}_{K}(v|_{K})$ for all $K\in\mathcal{T}_{h}$. The
consistency error $\delta_{h}\in(\widehat{V}_{h0}^{k})^{\prime}$ can now be
defined as in (44) and it can be bounded as in Lemma 4.5. Finally, the error
estimate and its proof are the same as those from Theorem 4.6 (and are not
repeated for brevity).

##  6 HHO method with Nitsche’s boundary penalty

In this section, we combine the HHO methods devised in the previous sections
with Nitsche’s boundary-penalty technique to enforce the boundary conditions
in a weak manner. For brevity, we only discuss the HHO method presented in
Section 3, but the following developments can be readily applied to the HHO
method from Section 5. To allow for a bit more generality, we detail here the
case of non-homogeneous boundary conditions. Thus, the model problem is as
follows:

$\Delta^{2}u=f\;\text{in $\Omega$},\qquad
u=g_{D},\;\partial_{n}u=g_{N}\;\text{in $\partial\Omega$},$ (58)

where the assumptions on the boundary data $g_{D}$ and $g_{N}$ are given in
Remark 2.1. We set ${\bm{G}}:=(\nabla u)|_{\partial\Omega}$ and notice that
${\bm{G}}$ is explicitly known in terms of the boundary data $g_{D}$ and
$g_{N}$ since ${\bm{G}}=g_{N}{\bm{n}}+(\partial_{t}g_{D}){\bm{t}}$.

Hinging on the ideas from [8, 7] for second-order elliptic PDEs, the HHO-
Nitsche (HHO-N) method devised in this section does not place any discrete
unknown on the mesh boundary faces, but only in the mesh cells and the mesh
interfaces. Thus, for every mesh cell $K\in\mathcal{T}_{h}$, we define the
subsets

$\partial K^{\mathrm{i}}:=\overline{\partial K\cap\Omega},\qquad\partial
K^{\mathrm{b}}:=\partial K\cap\partial\Omega,$ (59)

as well as $\mathcal{F}_{\partial K^{\mathrm{i}}}:=\mathcal{F}_{\partial
K}\cap\mathcal{F}_{h}^{\mathrm{i}}$ and $\mathcal{F}_{\partial
K^{\mathrm{b}}}:=\mathcal{F}_{\partial K}\cap\mathcal{F}_{h}^{\mathrm{b}}$.
The mesh cells having at least one boundary face are collected in the subset
$\mathcal{T}_{h}^{\mathrm{b}}:=\\{K\in\mathcal{T}_{h}\;|\;\mathcal{F}_{\partial
K^{\mathrm{b}}}\neq\emptyset\\}$, and we set
$\mathcal{T}_{h}^{\mathrm{i}}:=\mathcal{T}_{h}\setminus\mathcal{T}_{h}^{\mathrm{b}}$.

Letting $k\geq 0$ be the polynomial degree, the local HHO-N space is such that
for all $K\in\mathcal{T}_{h}$,

$\widehat{V}_{K}^{k}:=\mathbb{P}^{k+2}(K)\times\mathbb{P}^{k+1}(\mathcal{F}_{\partial
K^{\mathrm{i}}})\times\mathbb{P}^{k}(\mathcal{F}_{\partial K^{\mathrm{i}}}),$
(60)

and the corresponding global HHO-N space is now defined as

$\widehat{V}_{h}^{k}:=\mathbb{P}^{k+2}(\mathcal{T}_{h})\times\mathbb{P}^{k+1}(\mathcal{F}_{h}^{\mathrm{i}})\times\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{i}}).$
(61)

###  6.1 Reconstruction, stabilization, discrete problem, and stability

The definition of the local reconstruction operator is slightly modified with
respect to (18). Indeed,
$R_{K}^{\mathrm{i}}:\widehat{V}_{K}^{k}\rightarrow\mathbb{P}^{k+2}(K)$ is now
such that for all $\widehat{v}_{K}\in\widehat{V}_{K}^{k}$,

$\displaystyle(\nabla^{2}R_{K}^{\mathrm{i}}(\widehat{v}_{K}),\nabla^{2}w)_{K}={}$
$\displaystyle(\nabla^{2}{v}_{K},\nabla^{2}w)_{K}+(v_{K}-v_{\partial
K},\partial_{n}\Delta w)_{\partial
K^{\mathrm{i}}}-(\partial_{n}v_{K}-\gamma_{\partial
K},\partial_{nn}w)_{\partial K^{\mathrm{i}}}$ (62)
$\displaystyle-(\partial_{t}(v_{K}-v_{\partial K}),\partial_{nt}w)_{\partial
K^{\mathrm{i}}}+(v_{K},\partial_{n}\Delta w)_{\partial K^{\mathrm{b}}}-(\nabla
v_{K},\nabla\partial_{n}w)_{\partial K^{\mathrm{b}}},$

for all $w\in\mathbb{P}^{k+2}(K)^{\perp}$ together with the condition
$(R_{K}^{\mathrm{i}}(\widehat{v}_{K}),\xi)_{K}=(v_{K},\xi)_{K}$ for all
$\xi\in\mathbb{P}^{1}(K)$. Equivalently, owing to the integration by parts
formula (3), we have

$(\nabla^{2}R_{K}^{\mathrm{i}}(\widehat{v}_{K}),\nabla^{2}w)_{K}=(v_{K},\Delta^{2}w)_{K}-(v_{\partial
K},\partial_{n}\Delta w)_{\partial K^{\mathrm{i}}}+(\gamma_{\partial
K},\partial_{nn}w)_{\partial K^{\mathrm{i}}}+(\partial_{t}v_{\partial
K},\partial_{nt}w)_{\partial K^{\mathrm{i}}}.$ (63)

Dropping the integral over $\partial K^{\mathrm{b}}$ for the three rightmost
terms in (63) is, loosely speaking, a consistent operation in the case of
homogeneous boundary conditions. In the general case, we need to lift the
boundary data in every mesh cell $K\in\mathcal{T}_{h}^{\mathrm{b}}$ by means
of the lifting operator
$\mathcal{L}_{K}:H^{2}(K)\rightarrow\mathbb{P}^{k+2}(K)$ such that for all
$v\in H^{2}(K)$,

$(\nabla^{2}\mathcal{L}_{K}(v),\nabla^{2}w)_{K}=-(v,\partial_{n}\Delta
w)_{\partial K^{\mathrm{b}}}+(\nabla v,\nabla{\partial_{n}w})_{\partial
K^{\mathrm{b}}},$ (64)

for all $w\in\mathbb{P}^{k+2}(K)^{\perp}$, together with the condition
$(\mathcal{L}_{K}(v),\xi)_{K}=0$ for all $\xi\in\mathbb{P}^{1}(K)$. Notice
that $\mathcal{L}_{K}(u|_{K})$ is fully computable from the boundary data
$g_{D}$ and $g_{N}$. For convenience, we set $\mathcal{L}_{K}(v):=0$ for all
$K\in\mathcal{T}_{h}^{\mathrm{i}}$.

The local stabilization bilinear form $S_{\partial K}$ is also slightly
modified and is now such that for all
$(\widehat{v}_{K},\widehat{w}_{K})\in\widehat{V}_{K}^{k}\times\widehat{V}_{K}^{k}$,
we have $S_{\partial
K}(\widehat{v}_{K},\widehat{w}_{K}):=S^{\mathrm{i}}_{\partial
K}(\widehat{v}_{K},\widehat{w}_{K})+S^{\mathrm{b}}_{\partial K}(v_{K},w_{K})$
with

$\displaystyle S^{\mathrm{i}}_{\partial
K}(\widehat{v}_{K},\widehat{w}_{K}):={}$ $\displaystyle
h_{K}^{-3}\big{(}J^{k+1}_{\partial K}(v_{\partial K}-v_{K}),J^{k+1}_{\partial
K}(w_{\partial K}-{w}_{K})\big{)}_{\partial K^{\mathrm{i}}}$
$\displaystyle+h_{K}^{-1}\big{(}\Pi^{k}_{\partial K}(\gamma_{\partial
K}-\partial_{n}{v}_{K}),\Pi^{k}_{\partial K}(\chi_{\partial
K}-\partial_{n}{w}_{K})\big{)}_{\partial K^{\mathrm{i}}},$ (65) $\displaystyle
S^{\mathrm{b}}_{\partial K}(v_{K},w_{K}):={}$
$\displaystyle{}h_{K}^{-3}\big{(}v_{K},{w}_{K}\big{)}_{\partial
K^{\mathrm{b}}}+h_{K}^{-1}\big{(}\nabla{v}_{K},\nabla{w}_{K}\big{)}_{\partial
K^{\mathrm{b}}},$ (66)

where $S^{\mathrm{b}}_{\partial K}$ represents the boundary-penalty
contribution and acts only on the cell components. We emphasize that
$S^{\mathrm{b}}_{\partial K}$ does not need to be scaled by a weighting
coefficient to be taken large enough. Finally, the local bilinear form $a_{K}$
is defined on $\widehat{V}_{K}^{k}\times\widehat{V}_{K}^{k}$ as in (21).

The discrete HHO-N problem is as follows: Find
$\widehat{u}_{h}\in\widehat{V}_{h}^{k}$ such that

$a_{h}(\widehat{u}_{h},\widehat{w}_{h})={\ell}_{h}(\widehat{w}_{h}),\qquad\forall\widehat{w}_{h}\in\widehat{V}_{h}^{k},$
(67)

where $a_{h}$ is still assembled cellwise as in (25) yielding

$a_{h}(\widehat{v}_{h},\widehat{w}_{h}):=\sum_{K\in\mathcal{T}_{h}}(\nabla^{2}R_{K}^{\mathrm{i}}(\widehat{v}_{K}),\nabla^{2}R_{K}^{\mathrm{i}}(\widehat{w}_{K}))_{K}+\sum_{K\in\mathcal{T}_{h}}S^{\mathrm{i}}_{\partial
K}(\widehat{v}_{K},\widehat{w}_{K})+\sum_{K\in\mathcal{T}_{h}^{\mathrm{b}}}S^{\mathrm{b}}_{\partial
K}(v_{K},w_{K}),$ (68)

whereas the linear form ${\ell}_{h}$ now acts as follows:

$\displaystyle\ell_{h}(\widehat{w}_{h}):=$
$\displaystyle\sum_{K\in\mathcal{T}_{h}}(f,{w}_{K})_{K}+\sum_{K\in\mathcal{T}_{h}^{\mathrm{b}}}\Big{\\{}h_{K}^{-3}(g_{D},w_{K})_{\partial
K^{\mathrm{b}}}+h_{K}^{-1}({\bm{G}},\nabla w_{K})_{\partial K^{\mathrm{b}}}$
(69) $\displaystyle+\big{(}g_{D},\partial_{n}\Delta
R_{K}^{\mathrm{i}}(\widehat{w}_{K})\big{)}_{\partial
K^{\mathrm{b}}}-\big{(}{\bm{G}},\nabla\partial_{n}R_{K}^{\mathrm{i}}(\widehat{w}_{K})\big{)}_{\partial
K^{\mathrm{b}}}\Big{\\}}.$

Notice that

$\ell_{h}(\widehat{w}_{h})=\sum_{K\in\mathcal{T}_{h}}(f,{w}_{K})_{K}+\sum_{K\in\mathcal{T}_{h}^{\mathrm{b}}}\Big{\\{}S^{\mathrm{b}}_{\partial
K}(u|_{K},w_{K})-(\nabla^{2}\mathcal{L}_{K}(u|_{K}),\nabla^{2}R_{K}^{\mathrm{i}}(\widehat{w}_{K}))_{K}\Big{\\}}.$
(70)

Notice also that
$\ell_{h}(\widehat{w}_{h})=\sum_{K\in\mathcal{T}_{h}}(f,{w}_{K})_{K}=(f,w_{\mathcal{T}_{h}})_{\Omega}$
if the boundary conditions are homogeneous (so that only the cell component of
$\widehat{w}_{h}$ is needed to assemble $\ell_{h}$). As in the previous
sections, the discrete problem (67) is amenable to static condensation,
whereby all the cell unknowns are eliminated locally in every mesh cell,
leading to a global problem where the only remaining unknowns are those
attached to the mesh interfaces.

It is easy to see that the local stability and boundedness property stated in
Lemma 4.1 still holds true in the updated $H^{2}$-seminorm

$|\widehat{v}_{K}|^{2}_{\widehat{V}_{K}^{k}}:=\|\nabla^{2}v_{K}\|_{K}^{2}+h_{K}^{-3}\|v_{\partial
K}-v_{K}\|_{\partial K^{\mathrm{i}}}^{2}+h_{K}^{-1}\|\gamma_{\partial
K}-\partial_{n}v_{K}\|_{\partial K^{\mathrm{i}}}^{2}+S^{\mathrm{b}}_{\partial
K}(v_{K},v_{K}),$ (71)

recalling that $S^{\mathrm{b}}_{\partial
K}(v_{K},v_{K})=h_{K}^{-3}\|v_{K}\|^{2}_{\partial
K^{\mathrm{b}}}+h_{K}^{-1}\|\nabla v_{K}\|^{2}_{\partial K^{\mathrm{b}}}$.
Since
$\|\widehat{v}_{h}\|_{\widehat{V}_{h}^{k}}^{2}:=\sum_{K\in\mathcal{T}_{h}}|\widehat{v}_{K}|^{2}_{\widehat{V}_{K}^{k}}$
defines a norm on the global HHO space $\widehat{V}_{h}^{k}$ defined in (61),
the discrete bilinear form $a_{h}$ is coercive on $\widehat{V}_{h}^{k}$, and
the discrete problem (67) is well-posed owing to the Lax–Milgram lemma.

###  6.2 Polynomial approximation, consistency and error estimate

For all $K\in\mathcal{T}_{h}$, we define the local reduction operator
$\mathcal{\widehat{I}}^{k}_{K}:H^{2}(K)\rightarrow\widehat{V}_{K}^{k}$ such
that for all $v\in H^{2}(K)$,

$\mathcal{\widehat{I}}^{k}_{K}(v):=\big{(}\Pi_{K}^{k+2}(v),J_{\partial
K^{\mathrm{i}}}^{k+1}(v),\Pi_{\partial
K^{\mathrm{i}}}^{k}({\bm{n}}_{K}{\cdot}\nabla
v)\big{)}\in\widehat{V}_{K}^{k},$ (72)

with obvious notation regarding the operators $J_{\partial
K^{\mathrm{i}}}^{k+1}$ and $\Pi_{\partial K^{\mathrm{i}}}^{k}$. Let us set

$\mathcal{E}_{K}^{\mathrm{i}}:=R_{K}^{\mathrm{i}}\circ\mathcal{\widehat{I}}^{k}_{K}:H^{2}(K)\rightarrow\mathbb{P}^{k+2}(K).$
(73)

A straightforward verification (omitted for brevity) shows that the operator
$\mathcal{E}_{K}:=\mathcal{E}_{K}^{\mathrm{i}}+\mathcal{L}_{K}:H^{2}(K)\rightarrow\mathbb{P}^{k+2}(K)$
coincides indeed with the $H^{2}$-elliptic projection defined in (38).
Therefore, owing to Lemma 4.4, there is $C$ such that for all
$K\in\mathcal{T}_{h}$ and all $v\in H^{2+s}(K)$, $s>\frac{3}{2}$,

$\|v-(\mathcal{E}^{\mathrm{i}}_{K}(v)+\mathcal{L}_{K}(v))\|_{\sharp,K}\leq
C\|v-\Pi_{K}^{k+2}(v)\|_{\sharp,K},$ (74)

where the $\|{\cdot}\|_{\sharp,K}$-norm is defined in (41). Moreover, by
restricting the arguments to the mesh interfaces in the proof of Lemma 4.3, we
infer that there is $C$ such that for all $K\in\mathcal{T}_{h}$ and all $v\in
H^{2}(K)$,

$S^{\mathrm{i}}_{\partial
K}(\mathcal{\widehat{I}}^{k}_{K}(v),\mathcal{\widehat{I}}^{k}_{K}(v))^{\frac{1}{2}}\leq
C\|\nabla^{2}(v-\Pi_{K}^{k+2}(v))\|_{K}.$ (75)

The global reduction operator
$\mathcal{\widehat{I}}_{h}^{k}:H^{2}(\Omega)\rightarrow\widehat{V}_{h}^{k}$ is
defined such that for all $v\in H^{2}(\Omega)$,

$\mathcal{\widehat{I}}_{h}^{k}(v):=\big{(}(\Pi_{K}^{k+2}(v))_{K\in\mathcal{T}_{h}},(J^{k+1}_{F}(v))_{F\in\mathcal{F}_{h}^{\mathrm{i}}},(\Pi_{F}^{k}({\bm{n}}_{F}{\cdot}\nabla
v))_{F\in\mathcal{F}_{h}^{\mathrm{i}}}\big{)}\in\widehat{V}_{h}^{k},$ (76)

recalling that $v$ and $\nabla v$ are single-valued on every
$F\in\mathcal{F}_{h}^{\mathrm{i}}$ for all $v\in H^{2}(\Omega)$. As above, the
local components of $\mathcal{\widehat{I}}_{h}^{k}(v)$ attached to $K$ and its
faces in $\mathcal{F}_{\partial K^{\mathrm{i}}}$ are
$\mathcal{\widehat{I}}^{k}_{K}(v|_{K})$ for all $K\in\mathcal{T}_{h}$. We
define the consistency error $\delta_{h}\in(\widehat{V}_{h}^{k})^{\prime}$
such that
$\langle\delta_{h},\widehat{w}_{h}\rangle:=\ell(w_{\mathcal{T}_{h}})-a_{h}(\mathcal{\widehat{I}}_{h}^{k}(u),\widehat{w}_{h})$,
for all $\widehat{w}_{h}\in\widehat{V}_{h}^{k}$, where
$\langle\cdot,\cdot\rangle$ now denotes the duality pairing between
$(\widehat{V}_{h}^{k})^{\prime}$ and $\widehat{V}_{h}^{k}$.

######  Lemma 6.1 (Consistency)

Assume that $u\in H^{2+s}(\Omega)$ with $s>\frac{3}{2}$. The following holds
true:

$\|\delta_{h}\|_{(\widehat{V}_{h}^{k})^{\prime}}:=\sup_{\widehat{w}_{h}\in\widehat{V}_{h}^{k}}\frac{|\langle\delta_{h},\widehat{w}_{h}\rangle|}{\|\widehat{w}_{h}\|_{\widehat{V}_{h}^{k}}}\leq
C\left(\sum_{K\in\mathcal{T}_{h}}\|u-\Pi_{K}^{k+2}(u)\|^{2}_{\sharp,K}\right)^{\frac{1}{2}}.$
(77)

###### Proof.

The proof is similar to that of Lemma 4.5, so we only sketch it. Let
$\widehat{w}_{h}\in\widehat{V}_{h}^{k}$ having local components
$(w_{K},w_{\partial K},\chi_{\partial K})$ for all $K\in\mathcal{T}_{h}$. On
the one hand, we have

$\displaystyle\ell_{h}(\widehat{w}_{h})+\sum_{K\in\mathcal{T}_{h}^{\mathrm{b}}}(\nabla^{2}\mathcal{L}_{K}(u|_{K}),\nabla^{2}(R_{K}^{\mathrm{i}}(\widehat{w}_{K})))_{K}$
$\displaystyle=\sum_{K\in\mathcal{T}_{h}}\Big{\\{}(\nabla^{2}u,\nabla^{2}{w}_{K})_{K}+(\partial_{n}\Delta
u,w_{K})_{\partial K^{\mathrm{b}}}-(\nabla\partial_{n}u,\nabla
w_{K})_{\partial K^{\mathrm{b}}}$ $\displaystyle\quad+(\partial_{n}\Delta
u,w_{K}-w_{\partial K})_{\partial
K^{\mathrm{i}}}-(\partial_{nn}u,\partial_{n}w_{K}-\chi_{\partial K})_{\partial
K^{\mathrm{i}}}-(\partial_{nt}u,\partial_{t}(w_{K}-w_{\partial K}))_{\partial
K^{\mathrm{i}}}\Big{\\}}$
$\displaystyle\quad+\sum_{K\in\mathcal{T}_{h}^{\mathrm{b}}}\Big{\\{}h_{K}^{-3}\left(u,w_{K}\right)_{\partial
K^{\mathrm{b}}}+h_{K}^{-1}\left(\nabla u,\nabla w_{K}\right)_{\partial
K^{\mathrm{b}}}\Big{\\}}.$

On the other hand, recalling that
$\mathcal{E}_{K}=\mathcal{E}_{K}^{\mathrm{i}}+\mathcal{L}_{K}$, we have

$\displaystyle
a_{h}(\mathcal{\widehat{I}}_{h}^{k}(u),\widehat{w}_{h})+\sum_{K\in\mathcal{T}_{h}^{\mathrm{b}}}(\nabla^{2}\mathcal{L}_{K}(u|_{K}),\nabla^{2}(R_{K}^{\mathrm{i}}(\widehat{w}_{K})))_{K}$
$\displaystyle=\sum_{K\in\mathcal{T}_{h}}\Big{\\{}(\nabla^{2}\mathcal{E}_{K}(u),\nabla^{2}R^{\mathrm{i}}_{K}(\widehat{w}_{K}))_{K}+S^{\mathrm{i}}_{\partial
K}(\mathcal{\widehat{I}}^{k}_{K}(u),\widehat{w}_{K})\Big{\\}}+\sum_{K\in\mathcal{T}_{h}^{\mathrm{b}}}S^{\mathrm{b}}_{\partial
K}(\Pi_{K}^{k+2}(u),w_{K}),$

so that we have

$\displaystyle
a_{h}(\mathcal{\widehat{I}}_{h}^{k}(u),\widehat{w}_{h})+\sum_{K\in\mathcal{T}_{h}^{\mathrm{b}}}(\nabla^{2}\mathcal{L}_{K}(u|_{K}),\nabla^{2}(R_{K}^{\mathrm{i}}(\widehat{w}_{K})))_{K}$
$\displaystyle=\sum_{K\in\mathcal{T}_{h}}\Big{\\{}(\nabla^{2}\mathcal{E}_{K}(u),\nabla^{2}{w}_{K})_{K}+(\partial_{n}\Delta\mathcal{E}_{K}(u),w_{K})_{\partial
K^{\mathrm{b}}}-(\nabla\partial_{n}\mathcal{E}_{K}(u),\nabla w_{K})_{\partial
K^{\mathrm{b}}}$
$\displaystyle\quad+(\partial_{n}\Delta\mathcal{E}_{K}(u),w_{K}-w_{\partial
K})_{\partial
K^{\mathrm{i}}}-(\partial_{nn}\mathcal{E}_{K}(u),\partial_{n}w_{K}-\chi_{\partial
K})_{\partial
K^{\mathrm{i}}}-(\partial_{nt}\mathcal{E}_{K}(u),\partial_{t}(w_{K}-w_{\partial
K}))_{\partial K^{\mathrm{i}}}$ $\displaystyle\quad+S^{\mathrm{i}}_{\partial
K}(\mathcal{\widehat{I}}^{k}_{K}(u),\widehat{w}_{K})\Big{\\}}+\sum_{K\in\mathcal{T}_{h}^{\mathrm{b}}}\Big{\\{}h_{K}^{-3}(\Pi_{K}^{k+2}(u),w_{K})_{\partial
K^{\mathrm{b}}}+h_{K}^{-1}(\nabla\Pi_{K}^{k+2}(u),\nabla w_{K})_{\partial
K^{\mathrm{b}}}\Big{\\}}.$

Defining the function $\eta$ cellwise as
$\eta|_{K}:=u|_{K}-\mathcal{E}_{K}(u|_{K})$ for all $K\in\mathcal{T}_{h}$, we
infer that

$\displaystyle\langle\delta_{h},\widehat{w}_{h}\rangle=\sum_{K\in\mathcal{T}_{h}}\Big{\\{}(\nabla^{2}\eta,\nabla^{2}{w}_{K})_{K}+(\partial_{n}\Delta\eta,w_{K})_{\partial
K^{\mathrm{b}}}-(\nabla\partial_{n}\eta,\nabla w_{K})_{\partial
K^{\mathrm{b}}}$ $\displaystyle+(\partial_{n}\Delta\eta,w_{K}-w_{\partial
K})_{\partial
K^{\mathrm{i}}}-(\partial_{nn}\eta,\partial_{n}w_{K}-\chi_{\partial
K})_{\partial
K^{\mathrm{i}}}-(\partial_{nt}\eta,\partial_{t}(w_{K}-w_{\partial
K}))_{\partial K^{\mathrm{i}}}-S^{\mathrm{i}}_{\partial
K}(\mathcal{\widehat{I}}^{k}_{K}(u),\widehat{w}_{K})\Big{\\}}$
$\displaystyle+\sum_{K\in\mathcal{T}_{h}^{\mathrm{b}}}\Big{\\{}h_{K}^{-3}(u-\Pi_{K}^{k+2}(u),w_{K})_{\partial
K^{\mathrm{b}}}+h_{K}^{-1}(\nabla(u-\Pi_{K}^{k+2}(u)),\nabla w_{K})_{\partial
K^{\mathrm{b}}}\Big{\\}}.$

All the terms on the right-hand side can now be bounded by means of the
Cauchy–Schwarz inequality. For the first, fourth, fifth, and sixth terms, we
use (74), for the seventh term (involving $S^{\mathrm{i}}_{\partial K}$), we
use (75), and for the eighth and ninth terms, we invoke the trace inequality
(12). ∎

We are now ready to establish our main error estimate.

######  Theorem 6.2 ($H^{2}$-error estimate)

Assume that $u\in H^{2+s}(\Omega)$ with $s>\frac{3}{2}$. The following holds
true:

$\sum_{K\in\mathcal{T}_{h}}\|\nabla^{2}(u-R^{\mathrm{i}}_{K}(\widehat{u}_{K})-\mathcal{L}_{K}(u))\|_{K}^{2}\leq
C\sum_{K\in\mathcal{T}_{h}}\|u-\Pi_{K}^{k+2}(u)\|_{\sharp,K}^{2}.$ (78)

Consequently, if $k\geq 1$, assuming $u|_{K}\in H^{k+3}(K)$ for all
$K\in\mathcal{T}_{h}$, we have

$\sum_{K\in\mathcal{T}_{h}}\|\nabla^{2}(u-R^{\mathrm{i}}_{K}(\widehat{u}_{K})-\mathcal{L}_{K}(u))\|_{K}^{2}\leq
C\sum_{K\in\mathcal{T}_{h}}\big{(}h_{K}^{k+1}|u|_{H^{k+3}(K)}\big{)}^{2},$
(79)

and if $k=0$, letting $\sigma:=\min(s-1,1)\in(\frac{1}{2},1]$, we have

$\sum_{K\in\mathcal{T}_{h}}\|\nabla^{2}(u-R^{\mathrm{i}}_{K}(\widehat{u}_{K})-\mathcal{L}_{K}(u))\|_{K}^{2}\leq
C\sum_{K\in\mathcal{T}_{h}}\big{(}h_{K}(|u|_{H^{3}(K)}+h_{K}^{\sigma}|u|_{H^{3+\sigma}(K)})\big{)}^{2}.$
(80)

###### Proof.

As in the proof of Theorem 4.6, one shows that

$\sum_{K\in\mathcal{T}_{h}}\|\nabla^{2}R^{\mathrm{i}}_{K}(\widehat{e}_{K})\|_{K}^{2}\leq
C\sum_{K\in\mathcal{T}_{h}}\|u-\Pi_{K}^{k+2}(u)\|^{2}_{\sharp,K},$

where
$\widehat{e}_{h}:=\mathcal{\widehat{I}}_{h}^{k}(u)-\widehat{u}_{h}\in\widehat{V}_{h}^{k}$
is the discrete error. Since
$u-R^{\mathrm{i}}_{K}(\widehat{u}_{K})-\mathcal{L}_{K}(u)=(u-\mathcal{E}^{\mathrm{i}}_{K}(u)-\mathcal{L}_{K}(u))+R^{\mathrm{i}}_{K}(\widehat{e}_{K})$
for all $K\in\mathcal{T}_{h}$, the triangle inequality combined with (74) and
the above bound on the discrete error proves (47). Finally, (79) and (80) are
established by invoking the same arguments as above. ∎

##  7 Numerical examples

In this section, we present numerical examples to illustrate the theoretical
results on the present HHO methods and also to compare their numerical
performance with respect to other methods from the literature.

###  7.1 Convergence rates and computational performance of HHO methods

We select $f$ on $\Omega:=(0,1)^{2}$ so that the exact solution to (1) is
$u(x,y)=\sin(\pi x)^{2}\sin(\pi y)^{2}$ with homogeneous boundary conditions.
We consider the two HHO methods analyzed above. For clarity, we term “HHO(A)”
the method introduced in Section 3 with discrete unknowns in
$\mathbb{P}^{k+2}(\mathcal{T}_{h})\times\mathbb{P}^{k+1}(\mathcal{F}_{h}^{\mathrm{i}})\times\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{i}})$
and “HHO(B)” the method introduced in Section 5 with discrete unknowns in
$\mathbb{P}^{k+2}(\mathcal{T}_{h})\times\mathbb{P}^{k+2}(\mathcal{F}_{h}^{\mathrm{i}})\times\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{i}})$.
Additionally, we consider the method termed “HHO(C)” mentioned in Remark 3.2
where the discrete unknowns are in
$\mathbb{P}^{k+1}(\mathcal{T}_{h})\times\mathbb{P}^{k+1}(\mathcal{F}_{h}^{\mathrm{i}})\times\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{i}})$.
We employ polynomial degrees $k\in\\{0,\ldots,5\\}$. Since we consider various
polynomial degrees, and despite an $hp$-analysis falls beyond the present
scope, we implement the stabilization terms in (20), (51), and (26) with
$h_{K}^{-1}$ replaced by $(k+1)^{2}h_{K}^{-1}$ for all $K\in\mathcal{T}_{h}$.
All the computations were run with Matlab R2018a on the NEF platform at INRIA
Sophia Antipolis Méditerranée using 12 cores, and all the linear systems after
static condensation are solved using the `backslash` function. The algorithm
for solving the symmetric positive definite linear systems is the Cholesky
factorization.

Figure 1: Two examples of polygonal (Voronoi-like) meshes with $64$ (left) and
$1{,}024$ polygons

.

Figure 2:  Convergence of HHO(A) method in $H^{2}$\- and $L^{2}$-(semi)norms
on polygonal and rectangular meshes.

Let us first verify the convergence rates obtained with the HHO(A) method with
$k\in\\{0,1,2,3\\}$. We consider a sequence of successively refined
rectangular meshes and a sequence of successively refined polygonal (Voronoi-
like) meshes (generated through the PolyMesher Matlab library [43]). Two
examples of polygonal meshes are shown in Figure 1 (in general, the cells do
not contain more than $8$ edges). We measure relative errors in the (broken)
$H^{2}$-seminorm and in the $L^{2}$-norm, both quantities being evaluated
using the reconstruction of the HHO solution cellwise. The errors are reported
in Figure 2 as a function of $\mathrm{DoFs}^{1/2}$, where $\mathrm{DoFs}$
denotes the total number of globally coupled discrete unknowns (that is, the
face unknowns). We observe that the $H^{2}$-error converges at the optimal
rate $O(h^{k+1})$, as predicted in Theorem 4.6. The $L^{2}$-error converges at
the optimal rate $O(h^{k+3})$, except for $k=0$ where the rate is only
$O(h^{2})$; all these rates are consistent with what can be expected from a
duality argument (not detailed herein for brevity; see [4, 38] for examples of
this argument for HHO and WG methods).

Figure 3:  Convergence of HHO(A,B,C) methods in $H^{2}$\- and
$L^{2}$-(semi)norms on polygonal meshes.

Figure 4:  Comparison of computational times (in seconds) spent in
reconstruction, stabilization, and static condensation for the three
HHO(A,B,C) methods on a polygonal mesh with $16{,}384$ cells and polynomial
degrees $k\in\\{0,\ldots,5\\}$.

Let us now compare the three HHO(A,B,C) methods. The same relative errors as
in Figure 2 are reported in Figure 3. The results show that the three HHO
methods converge with the same rates, and that the accuracy reached on a given
mesh with a given polynomial degree is quite close for the three methods. We
mention that the three methods are sensitive to conditioning issues that arise
for high polynomial degree when the error is already quite low (typically
below $10^{-8}$ in the $H^{2}$-seminorm), and the HHO(C) method is somewhat
more sensitive. It is instructive to have a closer look at how the
computational costs related to the assembling of the system matrix are spent
between the tasks of reconstruction, stabilization, and static condensation.
The results are reported in Figure 4 on a polygonal mesh with $16{,}384$ cells
(and $49{,}014$ edges) and polynomial degrees $k\in\\{0,\ldots,5\\}$. Quite
importantly, the local reconstruction operator is computed based on equation
(19). Indeed, using (18) instead results in a more intricate assembling of the
right-hand side, increasing by a factor ranging from 2.5 (for $k=0$) to 3.5
(for $k=5$) the time spent in reconstruction. Figure 4 shows that the time
spent on static condensation is always marginal. Moreover, we can see that the
somewhat more elaborate design of the stabilization in the HHO(C) method is
reflected by a somewhat larger computational cost than for the HHO(A,B)
methods. The (perhaps a bit unexpected) consequence is that the HHO(A,B)
methods require altogether less assembly time than the HHO(C) method although
their number of discrete unknowns is larger. Finally, we notice that the
reconstruction time is always larger than the stabilization time, and this
trend gets more pronounced for larger $k$. To sum up, the most computationally
effective method based on these results is HHO(A). In what follows, we only
consider this method and simply call it “HHO” method.

###  7.2 Comparison with DG, $C^{0}$-IPDG, and FEM

In this section, we compare the computational performance of HHO with the
fully nonconforming dG method on polygonal and simpler meshes, and with the
$C^{0}$-IPDG, Morley, and HCT methods on triangular meshes.

HHO, triangular mesh | dG, triangular mesh
---|---
order |  $\\#$ DoFs | assembling | solving | order |  $\\#$ DoFs | assembling | solving
$k=0$ | $146{,}688$ | $275.5$ | $7.0$ | $\ell=2$ | $196{,}608$ | $472.0$ | $16.0$
$k=1$ | $244{,}480$ | $882.0$ | $18.3$ | $\ell=3$ | $327{,}680$ | $1300.9$ | $41.1$
$k=2$ | $342{,}272$ | $2076.3$ | $33.2$ | $\ell=4$ | $491{,}520$ | $2965.4$ | $96.0$
$k=3$ | $440{,}064$ | $4062.0$ | $53.0$ | $\ell=5$ | $688{,}128$ | $5940.6$ | $195.1$

HHO, polygonal mesh | dG, polygonal mesh
---|---
order |  $\\#$ DoFs | assembling | solving | order |  $\\#$ DoFs | assembling | solving
$k=0$ | $145{,}554$ | $251.1$ | $17.3$ | $\ell=2$ | $98{,}304$ | $420.5$ | $12.7$
$k=1$ | $242{,}590$ | $770.2$ | $44.7$ | $\ell=3$ | $163{,}840$ | $1160.9$ | $33.1$
$k=2$ | $339{,}626$ | $1784.3$ | $86.9$ | $\ell=4$ | $245{,}760$ | $2647.3$ | $78.4$
$k=3$ | $436{,}662$ | $3496.7$ | $149.9$ | $\ell=5$ | $344{,}064$ | $5304.7$ | $155.8$

Table 2: Comparison of total DoFs, assembling time, and solving time for the
HHO and dG methods. The polynomial degree is chosen so that both methods
deliver the same decay rates on the $H^{2}$-error. Upper table: triangular
mesh composed of $32{,}768$ cells; lower table: polygonal mesh composed of
$16{,}384$ cells.

Let us consider first the dG method. To put HHO and dG on a fair comparison
basis, we compare the HHO method with face polynomial degree $k\geq 0$ to the
dG method with cell polynomial degree $\ell=k+2$, so that both methods deliver
the same decay rates on the $H^{2}$-error. A comparison of total DoFs,
assembling time (including static condensation if applicable), and solving
time for both methods is provided in Table 2. We consider a triangular mesh
and a polygonal mesh (with $32{,}768$ and $16{,}384$ cells, respectively). The
first observation is that HHO always leads to less DoFs, and to smaller times
spent on assembling. The main reason is that the HHO DoFs are attached to the
mesh faces rather than the mesh cells. Although there are more faces than
cells in a given mesh (the more so as the cells are polygons with many faces),
the polynomial spaces in cells are richer than those on faces. Moreover, the
degree of the cell polynomials in the dG method is larger than the degree of
the face polynomials in the HHO method ($(k+2)$ vs. $\\{k,k+1\\}$). Another
reason for the lower assembling times with HHO is that the evaluation of
numerical fluxes in dG methods actually leads to a more expensive evaluation
of face-related quantities. The conclusions are, however, slightly different
if one considers the solving time (since the assembling stage can be fully
parallelized, the solving time becomes dominant in highly parallel
architectures). The results in Table 2 show that on triangular meshes (where
cells have a moderate number of faces), the solving time for HHO is always
smaller than that for dG. Instead, on polygonal (Voronoi-like) meshes, the
solving time for dG is smaller for low polynomial degrees (up to $2$), whereas
the solving time for HHO becomes again smaller for higher polynomial degrees.
The observation on polygonal meshes and low polynomial degrees indicates that
although the stencil of HHO methods is quite compact, it is still less compact
than that of dG methods. In particular, all the discrete unknowns attached to
the faces sharing a given mesh cell are coupled. Figure 5 provides a more
thorough viewpoint on the above results by highlighting the relative
efficiency of both methods measured as the time needed to reach a certain
error threshold in the $H^{2}$-seminorm. The time is either the assembling
time (which is more representative of a serial implementation) or the solving
time (which is more representative of a parallel implementation). We can see
that on triangular and rectangular meshes, for all polynomial orders, the HHO
method reaches an error threshold with less assembling or solving time than
the dG method. The same conclusion is reached on polygonal meshes for the
polynomial degree $k=3$ and both times as well as for $k\in\\{0,1,2\\}$ and
assembling time, whereas for $k\in\\{0,1,2\\}$ and solving time, the
efficiency of both methods is comparable.

Figure 5:  Comparison of HHO and dG methods: relative $H^{2}$-seminorm error as a function of assembling time (upper row) and solving time (lower row) on a sequence of polygonal (left), rectangular (center), and triangular (right) meshes and polynomial order $k\in\\{0,1,2,3\\}$ for HHO and $\ell=k+2$ for dG. $k=0$ |  $\\#$ DoFs | assembling | solving | $k=1$ |  $\\#$ DoFs | assembling | solving
---|---|---|---|---|---|---|---
Morley | $65{,}025$ | 22.9 | 4.3 | HCT | 97,283 | 169.3 | 19.8
HHO | $146{,}688$ | 275.5 | 7.0 | HHO | 244,480 | 882.0 | 18.3
$C^{0}$-IPDG | $65{,}025$ | 369.5 | 9.3 |  $C^{0}$-IPDG | 130,560 | 1318.8 | 27.0
Table 3: Comparison of total DoFs, assembling time, and solving time for the
HHO, $C^{0}$-IPDG, Morley, and HCT methods. The polynomial degree is chosen so
that all the methods in the same column deliver the same decay rates on the
$H^{2}$-error. Triangular mesh composed of $32{,}768$ cells, $49{,}408$ edges,
and $16{,}641$ vertices.

Let us now compare the efficiency of the HHO method to the $C^{0}$-IPDG,
Morley, and HCT methods on a sequence of successively refined triangulations
with $32$, $128$, $512$, $2{,}048$, $8{,}192$, and $32{,}768$ cells. As above,
the comparison is made between methods delivering the same decay rates on the
$H^{2}$-error. This means that the HHO method with polynomial degree $k\geq 0$
is compared with the $C^{0}$-IPDG with degree $\ell=k+2$. Moreover, the
HHO($k=0$) and the $C^{0}$-IPDG($\ell=2$) methods are compared with the Morley
element, and the HHO($k=1$) and the $C^{0}$-IPDG($\ell=3$) methods are
compared with the HCT element. Table 3 reports the total number of DoFs, the
assembling time, and the solving time for all the methods on the finest
triangular mesh. We can see that in the lowest-order case, both the assembling
and solving times for the Morley element are (much) smaller than those for the
HHO($k=0$) method, which are, in turn, smaller than those for the
$C^{0}$-IPDG($\ell=2$) method. The conclusion for the higher-order case is the
same concerning the lower times for HHO($k=1$) with respect to
$C^{0}$-IPDG($\ell=3$), whereas only the assembling time for HCT is (much)
smaller than that for HHO($k=1$), the solving time being instead comparable.
One reason for this good performance of HHO compared with HCT can be that the
stencil of HCT leads to a more dense system matrix, as a result of the method
attaching DoFs to the mesh vertices. Figure 6 reports the error measured in
the $H^{2}$-seminorm as a function of assembling and solving time, thereby
providing a comparison of the efficiency of the various methods on all the
considered triangulations. We notice that the Morley element is the most
efficient among the lowest-order methods, whereas the efficiency of the HHO
method is better than that of $C^{0}$-IPDG, and it is better than that of the
HCT element if the solving time is considered, whereas the conclusion is
reverted if the assembling time is considered.

Figure 6:  Comparison of HHO, $C^{0}$-IPDG, Morley, and HCT methods: relative
$H^{2}$-seminorm error as a function of assembling time (upper row) and
solving time (lower row) on a sequence of triangular meshes.

###  7.3 Tests on the HHO-Nitsche method

To conclude, let us briefly illustrate that the proposed HHO-Nitsche (HHO-N)
method with a weak enforcement of the boundary conditions performs as well as
the HHO method with a strong enforcement of the boundary conditions. We select
$f$ and the non-homogeneous boundary data $g_{D}$ and $g_{N}$ such that on
$\Omega:=(0,1)^{2}$, the exact solution is $u(x,y)=\sin(\pi x)^{2}\sin(\pi
y)^{2}+\exp{(-(x-0.5)^{2}-(y-0.5)^{2})}$. We consider the same sequence of
polygonal meshes and the same polynomial degrees as in Section 7.1. Figure 7
presents the relative errors measured in the $H^{2}$-seminorm and the
$L^{2}$-norm using cellwise the reconstruction operator for their evaluation.
We compare the HHO and HHO-N methods. Both methods employ the same number of
globally coupled DoFs. We can see from Figure 7 that the errors produced by
both methods are quite close in all cases.

Figure 7:  Convergence of HHO and HHO-N methods in $H^{2}$\- and
$L^{2}$-(semi)norms on polygonal meshes.

### Acknowledgment

The use of the NEF computing platform at Inria Sophia Antipolis Méditerranée
is gratefully acknowledged.

## References

  * [1] M. Abbas, A. Ern, and N. Pignet, Hybrid High-Order methods for finite deformations of hyperelastic materials, Comput. Mech., 62 (2018), pp. 909–928.
  * [2] , A Hybrid High-Order method for incremental associative plasticity with small deformations, Comput. Methods Appl. Mech. Engrg., 346 (2019), pp. 891–912.
  * [3] P. F. Antonietti, G. Manzini, and M. Verani, The fully nonconforming virtual element method for biharmonic problems, Math. Models Methods Appl. Sci., 28 (2018), pp. 387–407.
  * [4] F. Bonaldi, D. A. Di Pietro, G. Geymonat, and F. Krasucki, A hybrid high-order method for Kirchhoff-Love plate bending problems, ESAIM Math. Model. Numer. Anal., 52 (2018), pp. 393–421.
  * [5] S. C. Brenner and L.-Y. Sung, $C^{0}$ interior penalty methods for fourth order elliptic boundary value problems on polygonal domains, J. Sci. Comput., 22/23 (2005), pp. 83–118.
  * [6] F. Brezzi and L. D. Marini, Virtual element methods for plate bending problems, Comput. Methods Appl. Mech. Engrg., 253 (2013), pp. 455–462.
  * [7] E. Burman, M. Cicuttin, G. Delay, and A. Ern, An unfitted hybrid high-order method with cell agglomeration for elliptic interface problems, SIAM J. Sci. Comput., 43 (2021), pp. A859–A882.
  * [8] E. Burman and A. Ern, An unfitted hybrid high-order method for elliptic interface problems, SIAM J. Numer. Anal., 56 (2018), pp. 1525–1546.
  * [9] V. Calo, M. Cicuttin, Q. Deng, and A. Ern, Spectral approximation of elliptic operators by the hybrid high-order method, Math. Comp., 88 (2019), pp. 1559–1586.
  * [10] C. Carstensen and N. Nataraj, Lowest-order equivalent nonstandard finite element methods for biharmonic plates, ESAIM Math. Model. Numer. Anal., 56 (2022), pp. 41–78.
  * [11] F. Chave, D. A. Di Pietro, and S. Lemaire, A discrete Weber inequality on three-dimensional hybrid spaces with application to the HHO approximation of magnetostatics, Math. Models Methods Appl. Sci., 32 (2022), pp. 175–207.
  * [12] G. Chen and M. Feng, A $C^{0}$-weak Galerkin finite element method for fourth-order elliptic problems, Numer. Methods Partial Differential Equations, 32 (2016), pp. 1090–1104.
  * [13] C. Chinosi and L. D. Marini, Virtual element method for fourth order problems: $L^{2}$-estimates, Comput. Math. Appl., 72 (2016), pp. 1959–1967.
  * [14] F. Chouly, A. Ern, and N. Pignet, A hybrid high-order discretization combined with Nitsche’s method for contact and Tresca friction in small strain elasticity, SIAM J. Sci. Comput., 42 (2020), pp. A2300–A2324.
  * [15] P. G. Ciarlet, The finite element method for elliptic problems, vol. 40 of Classics in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002.
  * [16] B. Cockburn, Static condensation, hybridization, and the devising of the HDG methods, in Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations, G. R. Barrenechea, F. Brezzi, A. Cangiani, and E. H. Georgoulis, eds., vol. 114 of Lecture Notes in Computational Science and Engineering, Switzerland, 2016, Springer, pp. 129–178.
  * [17] B. Cockburn, D. A. Di Pietro, and A. Ern, Bridging the Hybrid High-Order and hybridizable discontinuous Galerkin methods, ESAIM Math. Model. Numer. Anal., 50 (2016), pp. 635–650.
  * [18] B. Cockburn, G. Fu, and F. J. Sayas, Superconvergence by $M$-decompositions. Part I: General theory for HDG methods for diffusion, Math. Comp., 86 (2017), pp. 1609–1641.
  * [19] D. A. Di Pietro, J. Droniou, and G. Manzini, Discontinuous Skeletal Gradient Discretisation Methods on polytopal meshes, J. Comput. Phys., 355 (2018), pp. 397–425.
  * [20] D. A. Di Pietro and A. Ern, Mathematical aspects of discontinuous Galerkin methods, vol. 69 of Mathématiques & Applications (Berlin) [Mathematics & Applications], Springer, Heidelberg, 2012.
  * [21] , A Hybrid High-Order locking-free method for linear elasticity on general meshes, Comput. Meth. Appl. Mech. Engrg., 283 (2015), pp. 1–21.
  * [22] D. A. Di Pietro, A. Ern, and S. Lemaire, An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators, Comput. Meth. Appl. Math., 14 (2014), pp. 461–472.
  * [23] D. A. Di Pietro and S. Krell, A Hybrid High-Order method for the steady incompressible Navier–Stokes problem, J. Sci. Comput., 74 (2018), pp. 1677–1705.
  * [24] Z. Dong and A. Ern, Hybrid high-order method for singularly perturbed fourth-order problems on curved domains, ESAIM Math. Model. Numer. Anal., 55 (2021), pp. 3091–3114.
  * [25] G. Engel, K. Garikipati, T. J. R. Hughes, M. G. Larson, L. Mazzei, and R. L. Taylor, Continuous/discontinuous finite element approximations of fourth-order elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity, Comput. Methods Appl. Mech. Engrg., 191 (2002), pp. 3669–3750.
  * [26] A. Ern and J.-L. Guermond, Finite element quasi-interpolation and best approximation, ESAIM Math. Model. Numer. Anal. (M2AN), 51 (2017), pp. 1367–1385.
  * [27] , Finite Elements I: Approximation and Interpolation, vol. 72 of Texts in Applied Mathematics, Springer Nature, Cham, Switzerland, 2021.
  * [28] , Finite Elements II: Galerkin Approximation, Elliptic and Mixed PDEs, vol. 73 of Texts in Applied Mathematics, Springer Nature, Cham, Switzerland, 2021.
  * [29] , Quasi-optimal nonconforming approximation of elliptic PDEs with contrasted coefficients and $H^{1+r}$, $r>0$, regularity, Found. Comput. Math. (Published online), (2021).
  * [30] A. Ern and P. Zanotti, A quasi-optimal variant of the hybrid high-order method for elliptic partial differential equations with $H^{-1}$ loads, IMA J. Numer. Anal., 40 (2020), pp. 2163–2188.
  * [31] E. H. Georgoulis and P. Houston, Discontinuous Galerkin methods for the biharmonic problem, IMA J. Numer. Anal., 29 (2009), pp. 573–594.
  * [32] V. Girault and P.-A. Raviart, Finite element methods for Navier-Stokes equations, vol. 5 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, 1986.  Theory and algorithms.
  * [33] P. Grisvard, Elliptic problems in nonsmooth domains, vol. 69 of Classics in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2011.
  * [34] C. Lehrenfeld, Removing the stabilization parameter in fitted and unfitted symmetric Nitsche formulations, in Proc. of ECCOMAS 2016, 2016.
  * [35] C. Lehrenfeld and J. Schöberl, High order exactly divergence-free hybrid discontinuous Galerkin methods for unsteady incompressible flows, Comput. Methods Appl. Mech. Engrg., 307 (2016), pp. 339–361.
  * [36] L. Morley, The triangular equilibrium element in the solution of plate bending problems, Aero. Quart., 19 (1968), pp. 149–169.
  * [37] I. Mozolevski and E. Süli, A priori error analysis for the $hp$-version of the discontinuous Galerkin finite element method for the biharmonic equation, Comput. Methods Appl. Math., 3 (2003), pp. 596–607.
  * [38] L. Mu, J. Wang, and X. Ye, Weak Galerkin finite element methods for the biharmonic equation on polytopal meshes, Numer. Methods Partial Differential Equations, 30 (2014), pp. 1003–1029.
  * [39] , A weak Galerkin finite element method with polynomial reduction, J. Comput. Appl. Math., 285 (2015), pp. 45–58.
  * [40] L. Mu, J. Wang, X. Ye, and S. Zhang, A $C^{0}$-weak Galerkin finite element method for the biharmonic equation, J. Sci. Comput., 59 (2014), pp. 473–495.
  * [41] C. Schwab, $p$– and $hp$–Finite element methods: Theory and applications in solid and fluid mechanics, Oxford University Press: Numerical mathematics and scientific computation, 1998.
  * [42] E. Süli and I. Mozolevski, $hp$-version interior penalty DGFEMs for the biharmonic equation, Comput. Methods Appl. Mech. Engrg., 196 (2007), pp. 1851–1863.
  * [43] C. Talischi, G. H. Paulino, A. Pereira, and I. F. M. Menezes, Polymesher: A general-purpose mesh generator for polygonal elements written in Matlab, Struct. Multidisc. Optim., 45 (2012), pp. 309–328.
  * [44] A. Veeser and P. Zanotti, Quasi-optimal nonconforming methods for symmetric elliptic problems. III—Discontinuous Galerkin and other interior penalty methods, SIAM J. Numer. Anal., 56 (2018), pp. 2871–2894.
  * [45] , Quasi-optimal nonconforming methods for symmetric elliptic problems. II—Overconsistency and classical nonconforming elements, SIAM J. Numer. Anal., 57 (2019), pp. 266–292.
  * [46] R. Verfürth, Error estimates for some quasi-interpolation operators, M2AN Math. Model. Numer. Anal., 33 (1999), pp. 695–713.
  * [47] M. Wang and J. Xu, The Morley element for fourth order elliptic equations in any dimensions, Numer. Math., 103 (2006), pp. 155–169.
  * [48] X. Ye, S. Zhang, and Z. Zhang, A new $P_{1}$ weak Galerkin method for the biharmonic equation, J. Comput. Appl. Math., 364 (2020), pp. 12337, 10\.
  * [49] R. Zhang and Q. Zhai, A weak Galerkin finite element scheme for the biharmonic equations by using polynomials of reduced order, J. Sci. Comput., 64 (2015), pp. 559–585.
  * [50] J. Zhao, S. Chen, and B. Zhang, The nonconforming virtual element method for plate bending problems, Math. Models Methods Appl. Sci., 26 (2016), pp. 1671–1687.
  * [51] J. Zhao, B. Zhang, S. Chen, and S. Mao, The Morley-type virtual element for plate bending problems, J. Sci. Comput., 76 (2018), pp. 610–629.