# Iterative solution to the biharmonic equation in mixed form discretized by
the Hybrid High-Order method111Funding: This research has been partially
funded by the European Union (ERC, NEMESIS, project number 101115663). Views
and opinions expressed are however those of the author(s) only and do not
necessarily reflect those of the European Union or the European Research
Council Executive Agency. PFA and MV have been partially funded by PRIN2020 n.
20204LN5N5 _“Advanced polyhedral discretisations of heterogeneous PDEs for
multiphysics problems”_ research grant, funded by the Italian Ministry of
Universities and Research (MUR). The present research is part of the
activities of “Dipartimento di Eccelllenza 2023-2027”. PFA and MV are members
of INdAM-GNCS.

P. F. Antonietti, P. Matalon22footnotemark: 2 222Corresponding author, M.
Verani22footnotemark: 2 MOX, Department of Mathematics, Politecnico di Milano,
Italy<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider the solution to the biharmonic equation in mixed form discretized
by the Hybrid High-Order (HHO) methods. The two resulting second-order
elliptic problems can be decoupled via the introduction of a new unknown,
corresponding to the boundary value of the solution of the first Laplacian
problem. This technique yields a global linear problem that can be solved
iteratively via a Krylov-type method. More precisely, at each iteration of the
scheme, two second-order elliptic problems have to be solved, and a normal
derivative on the boundary has to be computed. In this work, we specialize
this scheme for the HHO discretization. To this aim, an explicit technique to
compute the discrete normal derivative of an HHO solution of a Laplacian
problem is proposed. Moreover, we show that the resulting discrete scheme is
well-posed. Finally, a new preconditioner is designed to speed up the
convergence of the Krylov method. Numerical experiments assessing the
performance of the proposed iterative algorithm on both two- and three-
dimensional test cases are presented.

Keywords: partial differential equations, biharmonic equation, hybrid high-
order.

##  1 Introduction

Let $\Omega\subset\mathbb{R}^{d}$, $d\in\\{2,3\\}$, be an open, bounded,
polyhedral domain with smooth boundary $\partial\Omega$. In this work, we
address the numerical approximation of the solution to the biharmonic equation

$\displaystyle\Delta^{2}\psi=\mathsf{f}\quad$ $\displaystyle\text{in }\Omega,$
(1a)
$\displaystyle\psi=\mathsf{g}_{\mathrm{D}},\;\partial_{\mathrm{n}}\psi=\mathsf{g}_{\mathrm{N}}\quad$
$\displaystyle\text{on }\partial\Omega,$ (1b)

where the load function $\mathsf{f}\in L^{2}(\Omega)$,
$\mathsf{g}_{\mathrm{D}}\in H^{\nicefrac{{1}}{{2}}}(\partial\Omega)$ and
$\mathsf{g}_{\mathrm{N}}\in H^{\nicefrac{{-1}}{{2}}}(\partial\Omega)$ are
prescribed. As usual, for $X\subset\overline{\Omega}$ and $s\in\mathbb{R}$, we
denote by $H^{s}(X)$ the standard Sobolev space of index $s$. For all
$\mathsf{g}\colon\partial\Omega\to\mathbb{R}$, we define the subspace
$H_{\mathsf{g}}^{s}(\Omega)\mathrel{\mathop{:}}=\\{v\in H^{s}(\Omega)\mid
v_{|\partial\Omega}=\mathsf{g}\\}$. We denote by $(\cdot,\cdot)$ (resp.
$\langle\cdot,\cdot\rangle$) the standard $L^{2}$-inner product in $\Omega$
(resp. on $\partial\Omega$). Moreover, $\partial_{\mathrm{n}}\cdot$ denotes
the outer normal derivative on $\partial\Omega$.

Equation (1) typically models the bending of a clamped plate supporting a
load. By introducing the unknown $\omega:=-\Delta\psi$, (1a) can be rewritten
into two second-order elliptic equations, yielding the mixed formulation

$\displaystyle-\Delta\omega$ $\displaystyle=\mathsf{f}$ $\displaystyle\text{in
}\Omega,$ (2) $\displaystyle-\Delta\psi$ $\displaystyle=\omega$
$\displaystyle\text{in }\Omega,$
$\displaystyle\psi=\mathsf{g}_{\mathrm{D}},\;\partial_{\mathrm{n}}\psi$
$\displaystyle=\mathsf{g}_{\mathrm{N}}$ $\displaystyle\text{on
}\partial\Omega.$

This mixed form naturally arises in fluid dynamics, where $\psi$ represents
the stream and $\omega$ the vorticity. In the plate bending model, $\psi$
represents the deflection and $\omega$ the bending moment or shear resultant
force. There are several advantages in employing the mixed formulation (2)
over the primal form (1). First, while the weak solution of the primal form
(1) is to be found in $H^{2}(\Omega)$, that of the mixed one lies in
$H^{1}(\Omega)$, for which approximation spaces are easier to be constructed.
Second, the splitting (2) allows, after introducing an additional unknown, the
use of fast, scalable solvers available for second-order elliptic equations.
In particular, we will consider here the Hybrid High-Order (HHO) method [8,
10, 9], a non-conforming polyhedral discretization allowing arbitrary
polynomial degrees of approximation, and exhibiting optimal convergence rates.
In this context, the mixed form (2) allows to take advantage of the recent
works on multigrid methods designed for the diffusion equation discretized by
HHO schemes; see [12, 13, 11, 14].

Recently, ad-hoc HHO discretizations of the primal formulation (1) have been
designed in [17], and [16] generalizes the methods to other boundary
conditions. See also [15] for HHO discretizations of singularly perturbed
fourth-order problems. Previously, related models have also been tackled with
HHO methods, such as the Kirchhoff–Love plate bending model problem in [2],
and the Cahn–Hilliard equation in [4]. The mixed form (2) leads to a saddle-
point algebraic system, for which special solvers and preconditioners have
been proposed, e.g., in [20, 26, 22]. Nonetheless, in this form, the Laplacian
equations are coupled, preventing the use of fast and scalable solvers
specifically designed for symmetric positive-definite (SPD) matrices. To
address this issue, techniques introducing new variables have been designed in
order to transform (2) into a series of decoupled problems, whose associated
algebraic formulations lead to symmetric and positive-definite algebraic
problems. Among them, we mention the method proposed by Glowinski, Ciarlet,
Raviart and Pironneau [5, 6, 19], where the authors introduce the unknown
$\lambda\mathrel{\mathop{:}}=\omega_{|\partial\Omega}$. In that setting, one
solves a sequence of Dirichlet problems and iteratively improves the solution
as to enforce the prescribed Neumann condition. In a symmetric fashion, the
technique proposed by Falk [18] introduces
$\lambda\mathrel{\mathop{:}}=\partial_{\mathrm{n}}\omega$, and one solves a
sequence of Neumann problems while iteratively improving the solution as to
enforce the prescribed Dirichlet condition. In both schemes, $\lambda$ is the
solution of a linear, symmetric, elliptic equation of the form
$\mathcal{L}(\lambda)=b$, in which the evaluation of $\mathcal{L}$ involves
the solution of two Laplacian problems. In the discrete setting, solving this
equation then corresponds to solving an _implicit_ linear system, i.e. whose
matrix is not explicitly known but one can compute its action on a vector.
This specific configuration is well suited for the use of iterative methods
where the operator is applied to a vector at each iteration, without requiring
explicit knowledge of the matrix coefficients. Gradient descent algorithms
and, more specifically, Krylov methods such as conjugate gradients are ideal
candidates in this setting.

In this work, we focus on the approach of Glowinski _et al._ [5, 6, 19], that
we recall in Section 2. The use of an iterative method for the solution of the
global linear problem yields an iterative scheme where each iteration consists
of three steps: (i) the solution of a Laplace problem; (ii) the solution of a
subsequent Laplacian problem, using the solution of step (i) as a source term;
(iii) the computation of the normal derivative on the boundary of the solution
at step (ii). In Section 3, we recall the HHO discretization of the Laplacian
problem with Dirichlet boundary conditions. Then, Section 4 proposes a
computable algorithm for the discrete normal derivative of the HHO solution of
a diffusion problem. This step is indeed the main ingredient to define the
global discrete scheme.

The work of Glowinski _et al._ was applied to the standard Finite Element
Method (FEM). While that setting does not raise any issue regarding the well-
posedness of the discrete problem, as we will show in Section 5, in the
context of HHO, the resulting problem requires to be stabilized. A
stabilization method is proposed and the well-posedness of the discrete
problem is proved.

In the context of two-dimensional FEMs, a preconditioner ensuring a
convergence rate independent of the mesh size, was proposed in [23]. It is,
however, restricted to two-dimensional problems. In Section 6, we therefore
propose a novel preconditioner, applicable to general polytopal meshes. The
main idea consists in building an approximate, sparse matrix of the problem,
where each column $j$ is computed by solving the two Laplacian problems and
evaluating the normal derivative only in a restricted neighbourhood of the
$j^{th}$ degree of freedom (DoF). Finally, numerical experiments are reported
in Section 7. Various types of two- and three-dimensional meshes, including
polygonal meshes, are used. The scheme exhibits a convergence rate scaling as
$\mathcal{O}(h^{k+2})$ in $L^{2}$-norm, where $k$ denotes the polynomial
degree corresponding to the face unknowns of the HHO method.

##  2 The continuous splitting of the biharmonic problem

Following [5, 6, 19], we start from equation (2) and we introduce the new
unknown $\lambda:=\omega_{|\partial\Omega}\in
H^{\nicefrac{{1}}{{2}}}(\partial\Omega)$. Supposing $\lambda$ known, $\omega$
and $\psi$ are successively recovered by solving the following Dirichlet
problems:

$\left\\{\begin{aligned} -\Delta\omega&=\mathsf{f}&\quad&\textrm{ in
}\Omega,\\\ \omega&=\lambda&\quad&\textrm{ on
}\partial\Omega,\end{aligned}\right.$ (3a) $\left\\{\begin{aligned}
-\Delta\psi&=\omega&\quad&\textrm{ in }\Omega,\\\
\psi&=\mathsf{g}_{\mathrm{D}}&\quad&\textrm{ on
}\partial\Omega.\end{aligned}\right.$ (3b)

While the Dirichlet datum $\mathsf{g}_{\mathrm{D}}$ is explicitly enforced on
the solution in (3b), the enforcement of the Neumann condition
$\partial_{\mathrm{n}}\psi=\mathsf{g}_{\mathrm{N}}$ defines a problem in the
unknown $\lambda$, which is derived in the following way. For all $\mu\in
H^{\nicefrac{{1}}{{2}}}(\partial\Omega)$, we denote by
$(\omega(\mu),\psi(\mu))$ the solution of

$\left\\{\begin{aligned} -\Delta\omega(\mu)&=\mathsf{f}&\quad&\textrm{ in
}\Omega,\\\ \omega(\mu)&=\mu&\quad&\textrm{ on
}\partial\Omega,\end{aligned}\right.\qquad\qquad\left\\{\begin{aligned}
-\Delta\psi(\mu)&=\omega(\mu)&\quad&\textrm{ in }\Omega,\\\
\psi(\mu)&=\mathsf{g}_{\mathrm{D}}&\quad&\textrm{ on
}\partial\Omega.\end{aligned}\right.$

By construction, the solution $(\omega,\psi)$ of (3) corresponds to
$(\omega(\lambda),\psi(\lambda))$. Solving (3) then boils down to finding
$\lambda\in H^{\nicefrac{{1}}{{2}}}(\partial\Omega)$ such that

$\partial_{\mathrm{n}}\psi(\lambda)=\mathsf{g}_{\mathrm{N}}.$ (4)

In order to derive a linear problem from (4), the constraints related to
$\mathsf{f}$ and $\mathsf{g}_{\mathrm{D}}$ are eliminated through the
introduction of
$(\omega_{0},\psi_{0})\mathrel{\mathop{:}}=(\omega(0),\psi(0))$ (the choice of
$\mu=0$ is arbitrary). For all $\mu\in
H^{\nicefrac{{1}}{{2}}}(\partial\Omega)$, we then denote by
$(\mathring{\omega}({\mu}),\mathring{\psi}({\mu}))$ the solution of the
following sequence of problems, with vanishing load and Dirichlet function:

$\left\\{\begin{aligned} -\Delta\mathring{\omega}({\mu})&=0&\quad&\textrm{ in
}\Omega,\\\ \mathring{\omega}({\mu})&=\mu&\quad&\textrm{ on
}\partial\Omega,\end{aligned}\right.\qquad\qquad\left\\{\begin{aligned}
-\Delta\mathring{\psi}({\mu})&=\mathring{\omega}({\mu})&\quad&\textrm{ in
}\Omega,\\\ \mathring{\psi}({\mu})&=0&\quad&\textrm{ on
}\partial\Omega,\end{aligned}\right.$ (5)

respectively. Equation (4) can now be reformulated as the linear problem

$\mathcal{L}(\lambda)=b,\qquad\qquad
b:=\partial_{\mathrm{n}}\psi_{0}-\mathsf{g}_{\mathrm{N}},$ (6)

where $\mathcal{L}\colon H^{\nicefrac{{1}}{{2}}}(\partial\Omega)\to
H^{\nicefrac{{-1}}{{2}}}(\partial\Omega)$ is the linear operator defined such
that for all $\mu\in H^{\nicefrac{{1}}{{2}}}(\partial\Omega)$,

$\mathcal{L}(\mu):=-\partial_{\mathrm{n}}\mathring{\psi}({\mu}).$ (7)

The operator $\mathcal{L}$ is proved to be continuous, symmetric and positive-
definite in [19, Lem. 2.1].

##  3 HHO discretization of the Laplacian problem

In this section, we briefly recall (see, e.g., [8, Chap. 2] for extended
details) the discrete HHO formulation of the following problem: find
$u\colon\Omega\to\mathbb{R}$ such that

$\left\\{\begin{aligned} -\Delta u&=f&\quad&\textrm{ in }\Omega,\\\
u&=g_{\mathrm{D}}&\quad&\textrm{ on }\partial\Omega,\end{aligned}\right.$ (8)

where $f\colon\Omega\to\mathbb{R}$ and $g_{\mathrm{D}}\in
H^{\nicefrac{{1}}{{2}}}(\partial\Omega)$. A weak solution of (8) is obtained
via the variational formulation: find $u\in H^{1}_{g_{\mathrm{D}}}(\Omega)$
such that

$a(u,v)=(f,v)\qquad\forall v\in H_{0}^{1}(\Omega),$ (9)

where the bilinear form $a$ is such that $a(v,w)=(\nabla v,\nabla w)$, for all
$v,w\in H^{1}(\Omega)$.

###  3.1 Mesh definition and notation

Let the couple $(\mathcal{T}_{h},\mathcal{F}_{h})$ define a mesh of the domain
$\Omega\subset\mathbb{R}^{d}$, $d\in\\{1,2,3\\}$: $\mathcal{T}_{h}$ is a set
of disjoint, open, polyhedral elements such that
$\bigcup_{T\in\mathcal{T}_{h}}\overline{T}=\overline{\Omega}$;
$\mathcal{F}_{h}$ is the set of element faces;
$h\mathrel{\mathop{:}}=\max_{T\in\mathcal{T}_{h}}h_{T}$ with $h_{T}$ denoting
the diameter of $T\in\mathcal{T}_{h}$. The mesh is assumed to match the
geometrical requirements of [8, Def. 1.4] and, when asymptotic behaviours are
studied, polytopal regular mesh sequences in the sense defined by [8, Def.
1.9] are considered. Let $\mathcal{T}_{h}^{\mathrm{B}}$ the subset of
$\mathcal{T}_{h}$ collecting the elements located at the boundary of the
domain. We also define the following subsets of $\mathcal{F}_{h}$:

  * •

$\mathcal{F}_{h}^{\mathrm{I}}$, collecting the interior faces;

  * •

$\mathcal{F}_{h}^{\mathrm{B}}$, collecting the boundary faces;

  * •

$\mathcal{F}_{T}$, collecting the faces of $T$, for all $T\in\mathcal{T}_{h}$.

We denote by $\mathbf{n}_{\partial T}$ the unit normal vector to $\partial T$
pointing outward of $T$. For all $T\in\mathcal{T}_{h}$ (resp.
$F\in\mathcal{F}_{h}$), we denote by $(\cdot,\cdot)_{T}$ (resp.
$\langle\cdot,\cdot\rangle_{F}$) the standard inner product of $L^{2}(T)$
(resp. $L^{2}(F)$) or $L^{2}(T)^{d}$. We also denote by
$\langle\cdot,\cdot\rangle_{\mathcal{F}_{T}}\mathrel{\mathop{:}}=\sum_{F\in\mathcal{F}_{T}}\langle\cdot,\cdot\rangle_{F}$

###  3.2 Local and broken polynomial spaces

The HHO method hinges on discrete unknowns representing polynomial functions
local to elements and faces. So, for all $m\in\mathbb{N}_{0}$ and all
$T\in\mathcal{T}_{h}$ (resp. $F\in\mathcal{F}_{h}$), we denote by
$\mathbb{P}^{m}(T)$ (resp. $\mathbb{P}^{m}(F)$) the space spanned by the
restriction to $T$ (resp. $F$) of $d$-variate polynomials of total degree
$\leq m$. From these local polynomial spaces, we can construct the following
broken polynomial spaces supported by the mesh and its skeleton:

$\displaystyle\mathbb{P}^{m}(\mathcal{T}_{h})$
$\displaystyle\mathrel{\mathop{:}}=\left\\{v_{\mathcal{T}_{h}}\mathrel{\mathop{:}}=(v_{T})_{T\in\mathcal{T}_{h}}\mid
v_{T}\in\mathbb{P}^{m}(T)\;\;\forall T\in\mathcal{T}_{h}\right\\},$
$\displaystyle\mathbb{P}^{m}(\mathcal{F}_{h})$
$\displaystyle\mathrel{\mathop{:}}=\left\\{v_{\mathcal{F}_{h}}\mathrel{\mathop{:}}=(v_{F})_{F\in\mathcal{F}_{h}}\mid
v_{F}\in\mathbb{P}^{m}(F)\;\;\forall F\in\mathcal{F}_{h}\right\\},$

respectively. The local space $\mathbb{P}^{m}(\mathcal{F}_{T})$ is defined
analogously for all $T\in\mathcal{T}_{h}$. For all cell or face $X$, we denote
by $\pi_{X}^{m}\colon L^{2}(X)\to\mathbb{P}^{m}(X)$ the local
$L^{2}$-orthogonal projector onto the space $\mathbb{P}^{m}(X)$. By patching
up those local projectors, we denote by $\pi_{\mathcal{T}_{h}}^{m}\colon
L^{2}(\Omega)\to\mathbb{P}^{m}(\mathcal{T}_{h})$ the piecewise
$L^{2}$-orthogonal projector onto $\mathbb{P}^{m}(\mathcal{T}_{h})$, and by
$\pi_{\mathcal{F}_{h}^{\mathrm{B}}}^{m}\colon
L^{2}(\partial\Omega)\to\mathbb{P}^{m}(\mathcal{F}_{h}^{\mathrm{B}})$ the
piecewise $L^{2}$-orthogonal projector onto
$\mathbb{P}^{m}(\mathcal{F}_{h}^{\mathrm{B}})$.

###  3.3 Discrete hybrid formulation

Given the polynomial degrees $k\in\mathbb{N}_{0}$ and $l\in\\{k,k+1\\}$. The
global and local spaces of _hybrid_ variables are defined as

$\displaystyle\underline{U}_{h}$
$\displaystyle\mathrel{\mathop{:}}=\left\\{\underline{v}_{h}\mathrel{\mathop{:}}=\left(v_{\mathcal{T}_{h}},v_{\mathcal{F}_{h}}\right)\in\mathbb{P}^{l}(\mathcal{T}_{h})\times\mathbb{P}^{k}(\mathcal{F}_{h})\right\\},$
$\displaystyle\underline{U}_{T}$
$\displaystyle\mathrel{\mathop{:}}=\left\\{\underline{v}_{T}\mathrel{\mathop{:}}=\left(v_{T},v_{\mathcal{F}_{T}}\right)\in\mathbb{P}^{l}(T)\times\mathbb{P}^{k}(\mathcal{F}_{T})\right\\}\quad\forall
T\in\mathcal{T}_{h},$

respectively. For any $\underline{v}_{h}\in\underline{U}_{h}$, we denote by
$\underline{v}_{T}\in\underline{U}_{T}$ its restriction to
$T\in\mathcal{T}_{h}$. Boundary data are strongly accounted for in the
following subspaces:

$\mathbb{P}^{k,g_{\mathrm{D}}}(\mathcal{F}_{h})\mathrel{\mathop{:}}=\left\\{v_{\mathcal{F}_{h}}\in\mathbb{P}^{k}(\mathcal{F}_{h})\mid
v_{F}=\pi_{F}^{k}\,g_{\mathrm{D}}\;\;\forall
F\in\mathcal{F}_{h}^{\mathrm{B}}\right\\},\qquad\underline{U}_{h,g_{\mathrm{D}}}\mathrel{\mathop{:}}=\mathbb{P}^{l}(\mathcal{T}_{h})\times\mathbb{P}^{k,g_{\mathrm{D}}}(\mathcal{F}_{h}).$

In particular, homogeneous Dirichlet conditions are strongly enforced in
$\underline{U}_{h,0}$. The global HHO bilinear form associated to the
variational formulation of problem (8) is defined as
$a_{h}\colon\underline{U}_{h}\times\underline{U}_{h}\to\mathbb{R}$ such that
$a_{h}(\underline{v}_{h},\underline{w}_{h})\mathrel{\mathop{:}}=\sum_{T\in\mathcal{T}_{h}}a_{T}(\underline{v}_{T},\underline{w}_{T})$
where the local bilinear form
$a_{T}\colon\underline{U}_{T}\times\underline{U}_{T}\to\mathbb{R}$ is defined
as

$a_{T}(\underline{v}_{T},\underline{w}_{T})\mathrel{\mathop{:}}=(\nabla
p_{T}^{k+1}\underline{v}_{T},\nabla
p_{T}^{k+1}\underline{w}_{T})_{T}+s_{T}(\underline{v}_{T},\underline{w}_{T})\qquad\forall
T\in\mathcal{T}_{h},$ (10)

In this expression, the first term is responsible for consistency while the
second is required to ensure stability of the scheme. The consistency term
involves the _local potential reconstruction operator_
$p_{T}^{k+1}\colon\underline{U}_{T}\to\mathbb{P}^{k+1}(T)$ defined such that,
for all $\underline{v}_{T}\in\underline{U}_{T}$, it satisfies

$\displaystyle(\nabla p_{T}^{k+1}\underline{v}_{T},\nabla
w)_{T}=-(v_{T},\Delta w)_{T}+\langle v_{\mathcal{F}_{T}},\nabla
w\cdot\mathbf{n}_{\partial T}\rangle_{\mathcal{F}_{T}}\qquad\quad\forall
w\in\mathbb{P}^{k+1}(T),$ (11a)
$\displaystyle(p_{T}^{k+1}\underline{v}_{T},1)_{T}=(v_{T},1)_{T}.$ (11b)

Given the local interpolant $\underline{v}_{T}\in\underline{U}_{T}$ of a
function $v\in L^{2}(\Omega)$, $p_{T}^{k+1}$ reconstructs an approximation of
$v$ of degree $k+1$. By patching together local contributions, we define the
global operator
$p_{h}^{k+1}\colon\underline{U}_{h}\to\mathbb{P}^{k+1}(\mathcal{T}_{h})$ such
that
$(p_{h}^{k+1}\underline{v}_{h})_{|T}\mathrel{\mathop{:}}=p_{T}^{k+1}\underline{v}_{T}$
for all $T\in\mathcal{T}_{h}$. The stabilization bilinear form $s_{T}$ depends
on its arguments through the _difference operators_
$\delta_{T}:\underline{U}_{T}\to\mathbb{P}^{k}(T)$ and
$\delta_{TF}:\underline{U}_{T}\to\mathbb{P}^{k}(F)$ for all
$F\in\mathcal{F}_{T}$, defined such that, for all
$\underline{v}_{T}\in\underline{U}_{T}$,

$\delta_{T}\underline{v}_{T}\mathrel{\mathop{:}}=\pi_{T}^{l}(p_{T}^{k+1}\underline{v}_{T}-v_{T})$
and
$\delta_{TF}\underline{v}_{T}\mathrel{\mathop{:}}=\pi_{F}^{k}(p_{T}^{k+1}\underline{v}_{T}-v_{F})$
for all $F\in\mathcal{F}_{T}$.  (12)

These operators capture the higher-order correction that the reconstruction
$p_{T}^{k+1}$ adds to the element and face unknowns, respectively. A classical
expression for
$s_{T}\colon\underline{U}_{T}\times\underline{U}_{T}\to\mathbb{R}$ is

$s_{T}(\underline{v}_{T},\underline{w}_{T})\mathrel{\mathop{:}}=\sum_{F\in\mathcal{F}_{T}}h_{F}^{-1}\;\mathfrak{s}_{TF}(\underline{v}_{T},\underline{w}_{T}),\qquad\qquad\mathfrak{s}_{TF}(\underline{v}_{T},\underline{w}_{T})\mathrel{\mathop{:}}=\langle(\delta_{TF}-\delta_{T})\underline{v}_{T},(\delta_{TF}-\delta_{T})\underline{w}_{T}\rangle_{F}.$
(13)

If $l=k+1$, we can also use the simpler formula (used, e.g., in [4])

$\mathfrak{s}_{TF}(\underline{v}_{T},\underline{w}_{T})\mathrel{\mathop{:}}=\langle\pi_{F}^{k}(v_{T}-v_{F}),\pi_{F}^{k}(w_{T}-w_{F})\rangle_{F}.$
(14)

The global HHO problem reads: find
$\underline{u}_{h}\in\underline{U}_{h,g_{\mathrm{D}}}$ such that

$a_{h}(\underline{u}_{h},\underline{v}_{h})=(f,v_{\mathcal{T}_{h}})\qquad\forall\underline{v}_{h}\in\underline{U}_{h,0}.$
(15)

The final approximation is obtained through the post-processing step

$u_{h}\mathrel{\mathop{:}}=p_{h}^{k+1}\underline{u}_{h}\in\mathbb{P}^{k+1}(\mathcal{T}_{h}).$
(16)

###  3.4 Local problems and cell unknown recovery operator

It follows from this local construction that the cell unknowns are only
locally coupled. For all $T\in\mathcal{T}_{h}$, we define the linear operator
$\Theta_{T}\colon\mathbb{P}^{k}(\mathcal{F}_{T})\to\mathbb{P}^{l}(T)$ such
that for all $v_{\mathcal{F}_{T}}\in\mathbb{P}^{k}(\mathcal{F}_{T})$,
$\Theta_{T}v_{\mathcal{F}_{T}}$ is the unique solution of the local problem

$a_{T}((\Theta_{T}v_{\mathcal{F}_{T}},0),(w_{T},0))=-a_{T}((0,v_{\mathcal{F}_{T}}),(w_{T},0))\qquad\forall
w_{T}\in\mathbb{P}^{l}(T).$ (17)

We denote by $\Theta_{T}$ the _cell unknown recovery operator_. In order to
shorten the notation of the hybrid couple
$\underline{v}_{T}:=(\Theta_{T}v_{\mathcal{F}_{T}},v_{\mathcal{F}_{T}})$, we
define the operator
$\underline{\Theta}_{T}\colon\mathbb{P}^{k}(\mathcal{F}_{T})\to\underline{U}_{T}$
such that for all $v_{\mathcal{F}_{T}}$,

$\underline{\Theta}_{T}v_{\mathcal{F}_{T}}:=(\Theta_{T}v_{\mathcal{F}_{T}},v_{\mathcal{F}_{T}}).$
(18)

Finally, the associated global operators
$\Theta_{\mathcal{T}_{h}}\colon\mathbb{P}^{k}(\mathcal{F}_{h})\to\mathbb{P}^{l}(\mathcal{T}_{h})$
and
$\underline{\Theta}_{h}\colon\mathbb{P}^{k}(\mathcal{F}_{h})\to\underline{U}_{h}$
are defined locally such that for all $T\in\mathcal{T}_{h}$:

$\displaystyle(\Theta_{\mathcal{T}_{h}}v_{\mathcal{F}_{h}})_{|T}\mathrel{\mathop{:}}=\Theta_{T}v_{\mathcal{F}_{T}}\qquad\forall
v_{\mathcal{F}_{h}}\in\mathbb{P}^{k}(\mathcal{F}_{h}),$ (19)
$\displaystyle(\underline{\Theta}_{h}v_{\mathcal{F}_{h}})_{|T\times\partial
T}\mathrel{\mathop{:}}=\underline{\Theta}_{T}v_{\mathcal{F}_{T}}\qquad\forall
v_{\mathcal{F}_{h}}\in\mathbb{P}^{k}(\mathcal{F}_{h}).$

$\underline{\Theta}_{h}$ verifies the following useful property:

######  Lemma 1.

For all $v_{\mathcal{F}_{h}}\in\mathbb{P}^{k}(\mathcal{F}_{h})$ and
$\underline{w}_{h}\in\underline{U}_{h}$, it holds that

$\displaystyle
a_{h}(\underline{\Theta}_{h}v_{\mathcal{F}_{h}},\underline{w}_{h})$
$\displaystyle=a_{h}(\underline{\Theta}_{h}v_{\mathcal{F}_{h}},\underline{\Theta}_{h}w_{\mathcal{F}_{h}}).$
(20)

###### Proof.

See Appendix A. ∎

######  Remark 1.

(Static condensation) In practice, problem (15) is solved through the
equivalent, condensed formulation: find
$u_{\mathcal{F}_{h}}\in\mathbb{P}^{k,g_{\mathrm{D}}}(\mathcal{F}_{h})$ such
that

$\widehat{a}_{h}(u_{\mathcal{F}_{h}},v_{\mathcal{F}_{h}})=(f,\Theta_{\mathcal{T}_{h}}v_{\mathcal{F}_{h}})\qquad\forall
v_{\mathcal{F}_{h}}\in\mathbb{P}^{k,0}(\mathcal{F}_{h}),$ (21)

where
$\widehat{a}_{h}\colon\mathbb{P}^{k}(\mathcal{F}_{h})\times\mathbb{P}^{k}(\mathcal{F}_{h})\to\mathbb{R}$
is such that for all
$v_{\mathcal{F}_{h}},w_{\mathcal{F}_{h}}\in\mathbb{P}^{k}(\mathcal{F}_{h})$,

$\widehat{a}_{h}(v_{\mathcal{F}_{h}},w_{\mathcal{F}_{h}}):=a_{h}(\underline{\Theta}_{h}v_{\mathcal{F}_{h}},\underline{\Theta}_{h}w_{\mathcal{F}_{h}}).$
(22)

Refer to [7, Prop. 4] for the proof.

##  4 Discrete normal derivative

In this section, we derive an approximation, in the HHO context, of the normal
derivative on $\partial\Omega$. This formula will play a crucial role in the
approximation scheme of the biharmonic problem.

Let $u$ be the solution of the boundary value problem (8) and let
$\mathcal{H}\colon H^{\nicefrac{{1}}{{2}}}(\partial\Omega)\to H^{1}(\Omega)$
be a linear operator such that for all $v\in
H^{\nicefrac{{1}}{{2}}}(\partial\Omega)$ defined on the boundary,
$\mathcal{H}v$ extends $v$ in the interior of $\Omega$. As the trace operator
is surjective, such an operator exists. By Green’s formula, it holds that

$\langle\partial_{\mathrm{n}}u,v\rangle=(\nabla u,\nabla\mathcal{H}v)+(\Delta
u,\mathcal{H}v)\qquad\forall v\in H^{\nicefrac{{1}}{{2}}}(\partial\Omega),$

where the $L^{2}(\partial\Omega)$-inner product notation
$\langle\cdot,\cdot\rangle$ is employed here to denote a duality pairing.
Using the bilinear form $a(\cdot,\cdot)$, and given that $-\Delta u=f$, the
above equation becomes

$\langle\partial_{\mathrm{n}}u,v\rangle=a(u,\mathcal{H}v)-(f,\mathcal{H}v)\qquad\forall
v\in H^{\nicefrac{{1}}{{2}}}(\partial\Omega).$ (23)

In the literature, the operator $g_{\mathrm{D}}\mapsto\partial_{\mathrm{n}}u$
is called Poincaré-Steklov operator, or Dirichlet-to-Neumann map; see, e.g.
[25]. Equation (23) is the well-known variational formula for the computation
of the normal derivative. Now, given the solution $\underline{u}_{h}$ of the
discrete problem (15), we define its normal derivative on the boundary faces,
denoted by
$\partial_{\mathrm{n},h}(\underline{u}_{h})\in\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})$,
through a discrete counterpart of (23) in the HHO setting. Namely,
$\partial_{\mathrm{n},h}(\underline{u}_{h})$ verifies

$\langle\partial_{\mathrm{n},h}(\underline{u}_{h}),v_{\mathcal{F}_{h}^{\mathrm{B}}}\rangle=a_{h}(\underline{u}_{h},\underline{\mathcal{H}}_{h}v_{\mathcal{F}_{h}^{\mathrm{B}}})-(f,\mathcal{H}_{\mathcal{T}_{h}}v_{\mathcal{F}_{h}^{\mathrm{B}}})\qquad\forall
v_{\mathcal{F}_{h}^{\mathrm{B}}}\in\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}}),$
(24)

where $\underline{\mathcal{H}}_{h}$ is an extension/lifting operator expressed
in the hybrid setting, i.e.

$\underline{\mathcal{H}}_{h}\mathrel{\mathop{:}}=(\mathcal{H}_{\mathcal{T}_{h}},\mathcal{H}_{\mathcal{F}_{h}}),\qquad\mathcal{H}_{\mathcal{T}_{h}}\colon\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})\to\mathbb{P}^{l}(\mathcal{T}_{h}),\qquad\mathcal{H}_{\mathcal{F}_{h}}\colon\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})\to\mathbb{P}^{k}(\mathcal{F}_{h}).$

For all
$v_{\mathcal{F}_{h}^{\mathrm{B}}}\in\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})$
and $F\in\mathcal{F}_{h}$, we define

$(\mathcal{H}_{\mathcal{F}_{h}}v_{\mathcal{F}_{h}^{\mathrm{B}}})_{|F}\mathrel{\mathop{:}}=\begin{cases}(v_{\mathcal{F}_{h}^{\mathrm{B}}})_{|F}&\text{
if }F\in\mathcal{F}_{h}^{\mathrm{B}},\\\ 0&\text{ otherwise}.\end{cases}$

Then, we set

$\mathcal{H}_{\mathcal{T}_{h}}\mathrel{\mathop{:}}=\Theta_{\mathcal{T}_{h}}\mathcal{H}_{\mathcal{F}_{h}},$

where $\Theta_{\mathcal{T}_{h}}$ is defined as in (19).

######  Remark 2.

(Condensed formula) In order to facilitate the practical computation of
$\partial_{\mathrm{n},h}(\underline{u}_{h})$, equation (24) can be rewritten
in terms of the _condensed_ bilinear form $\widehat{a}_{h}$ (cf. (22)). First
of all, it follows from the definition of $\underline{\Theta}_{h}$ that
$\underline{\mathcal{H}}_{h}=\underline{\Theta}_{h}\mathcal{H}_{\mathcal{F}_{h}}$.
Then, by using the symmetry of $a_{h}$ with property (20) and the definition
(22) of $\widehat{a}_{h}$, we have

$a_{h}(\underline{u}_{h},\underline{\mathcal{H}}_{h}v_{\mathcal{F}_{h}^{\mathrm{B}}})=a_{h}(\underline{u}_{h},\underline{\Theta}_{h}\mathcal{H}_{\mathcal{F}_{h}}v_{\mathcal{F}_{h}^{\mathrm{B}}})=a_{h}(\underline{\Theta}_{h}u_{\mathcal{F}_{h}},\underline{\Theta}_{h}\mathcal{H}_{\mathcal{F}_{h}}v_{\mathcal{F}_{h}^{\mathrm{B}}})=\widehat{a}_{h}(u_{\mathcal{F}_{h}},\mathcal{H}_{\mathcal{F}_{h}}v_{\mathcal{F}_{h}^{\mathrm{B}}}).$

Equation (24) then becomes

$\langle\partial_{\mathrm{n},h}(\underline{u}_{h}),v_{\mathcal{F}_{h}^{\mathrm{B}}}\rangle=\widehat{a}_{h}(u_{\mathcal{F}_{h}},\mathcal{H}_{\mathcal{F}_{h}}v_{\mathcal{F}_{h}^{\mathrm{B}}})-(f,\mathcal{H}_{\mathcal{T}_{h}}v_{\mathcal{F}_{h}^{\mathrm{B}}})\qquad\forall
v_{\mathcal{F}_{h}^{\mathrm{B}}}\in\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}}).$
(25)

Notice that compared to (24), which requires the knowledge of
$\underline{u}_{h}$, formula (25) only involves $u_{\mathcal{F}_{h}}$.
Consequently, the following hold: (i) $u_{\mathcal{T}_{h}}$ does not have to
be computed, (ii) the matrix used to compute the first term is smaller, as the
matrix representation of $\widehat{a}_{h}$ is a Schur complement where the
cell unknowns have been eliminated.

The computation of the discrete normal derivative as proposed here is
numerically validated in Section 7.2. The experiments show a convergence in
$\mathcal{O}(h^{k+1})$ in $L^{2}$-norm.

######  Remark 3.

(Choice of discrete normal derivative) In this scheme, the normal derivative
on the boundary is defined by a discrete version of (23), i.e. through an
integration by parts involving a lifting operator, which is the way employed
by Glowinski _et al._ in their original scheme [19] with standard FEM.
However, contrary to FEM, HHO methods natively include a flux formulation (cf.
[8, Lem. 2.25]), which we could lean on instead of (23). However, the
structure of formula (23) directly interplays the two Laplacian problems in
variational form, allowing a natural reformulation of the problem through a
symmetric positive-definite bilinear form (see Lemma 3).

##  5 Discrete HHO problem

In this section, the global, discrete realization of problem (6) is
formulated.

###  5.1 Bilinear form

In order to devise a stable problem, we will make use of the following
bilinear form in the hybrid space $\underline{U}_{h}$. For all
$\underline{v}_{h},\underline{w}_{h}\in\underline{U}_{h}$, we introduce

$(\underline{v}_{h},\underline{w}_{h})^{\star}\mathrel{\mathop{:}}=(v_{\mathcal{T}_{h}},w_{\mathcal{T}_{h}})+\sum_{T\in\mathcal{T}_{h}^{\mathrm{B}}}\sum_{F\in\mathcal{F}_{T}}h_{F}\;\mathfrak{s}_{TF}^{\star}(\underline{v}_{T},\underline{w}_{T}),\qquad\quad\mathfrak{s}_{TF}^{\star}(\underline{v}_{T},\underline{w}_{T})\mathrel{\mathop{:}}=\langle\pi_{F}^{k}(v_{T}-v_{F}),\pi_{F}^{k}(w_{T}-w_{F})\rangle_{F}.$
(26)

Formula (26) defines an inner product-like bilinear form based on the
$L^{2}$-inner product of the cell unknowns, to which a stabilizing term has
been added. The stabilizing term is inspired from [4, Eq. (26)]. Remark that
$\mathfrak{s}_{TF}^{\star}$ is the same as (14), that one can employ to
stabilize the Laplacian bilinear form if $l=k+1$. The scaling factor $h_{F}$
is selected to ensure dimensional homogeneity with the consistency term it is
added to. Note that $(\cdot,\cdot)^{\star}$ does not define an inner product
in $\underline{U}_{h}$, due to the fact that only the boundary cells are
involved in the stabilizing term. In spite of that, we allow ourselves the use
of an inner product notation, since $(\cdot,\cdot)^{\star}$ shall undertake
the role of the $L^{2}$-inner product in the discrete problem and, thus,
carries the same semantics.

$(\cdot,\cdot)^{\star}$ remains symmetric and positive semi-definite. If
definiteness is not globally ensured, it is ensured “locally” in the cells
where the stabilization term is applied. One might say that
$(\cdot,\cdot)^{\star}$ enjoys a property of _local stability_ in the boundary
cells, which is formalized by

######  Lemma 2.

(Local stability on the boundary)

$(\underline{v}_{h},\underline{v}_{h})^{\star}=0\qquad\Longrightarrow\qquad\underline{v}_{T}=0\quad\forall
T\in\mathcal{T}_{h}^{\mathrm{B}}.$ (27)

###### Proof.

Let $T\in\mathcal{T}_{h}^{\mathrm{B}}$. Assume
$\underline{v}_{h}\in\underline{U}_{h}$ such that
$(\underline{v}_{h},\underline{v}_{h})^{\star}=0$. Then, by definition (26),

$(v_{T},v_{T})_{T}+\sum_{F\in\mathcal{F}_{T}}h_{F}\;\mathfrak{s}_{TF}^{\star}(\underline{v}_{T},\underline{v}_{T})=0.$

As $\mathfrak{s}_{TF}^{\star}$ is positive semi-definite, all terms are non-
negative. Consequently, $v_{T}=0$ and
$\mathfrak{s}_{TF}^{\star}(\underline{v}_{T},\underline{v}_{T})=0$ for all
$F\in\mathcal{F}_{T}$. We then have

$0=\mathfrak{s}_{TF}^{\star}(\underline{v}_{T},\underline{v}_{T})=\|\pi_{F}^{k}(v_{T}-v_{F})\|_{F}^{2}=\|v_{F}\|_{F}^{2},$

where we have used the facts that $v_{T}=0$ and $\pi_{F}^{k}v_{F}=v_{F}$. We
conclude that $v_{F}=0$. ∎

Let us now describe the discrete, variational form of the continuous operator
(7). For all $\mu\in\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})$, we denote
by
$(\mathring{\underline{\omega}}_{h}^{\mu},\mathring{\underline{\psi}}_{h}^{\mu})\in\underline{U}_{h,\mu}\times\underline{U}_{h,0}$
the solution of the discrete problems

$\displaystyle
a_{h}(\mathring{\underline{\omega}}_{h}^{\mu},\underline{v}_{h})$
$\displaystyle=0$
$\displaystyle\forall\underline{v}_{h}\in\underline{U}_{h,0},$ (28a)
$\displaystyle a_{h}(\mathring{\underline{\psi}}_{h}^{\mu},\underline{v}_{h})$
$\displaystyle=(\mathring{\underline{\omega}}_{h}^{\mu},\underline{v}_{h})^{\star}$
$\displaystyle\forall\underline{v}_{h}\in\underline{U}_{h,0}.$ (28b)

For $\mu,\eta\in\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})$, and
$(\mathring{\underline{\omega}}_{h}^{\mu},\mathring{\underline{\psi}}_{h}^{\mu})$
the solution of (28) associated with $\mu$, we define the bilinear form
$\ell_{h}\colon\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})\times\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})\to\mathbb{R}$
as the computation of
$-\partial_{\mathrm{n},h}(\mathring{\underline{\psi}}_{h}^{\mu})$ via formula
(24), such that

$\ell_{h}(\mu,\eta)\mathrel{\mathop{:}}=-a_{h}(\mathring{\underline{\psi}}_{h}^{\mu},\underline{\mathcal{H}}_{h}\eta)+(\mathring{\underline{\omega}}_{h}^{\mu},\underline{\mathcal{H}}_{h}\eta)^{\star}.$
(29)

Remark that we have used $(\cdot,\cdot)^{\star}$ in (28b) and (29) instead of
$(\cdot,\cdot)$ in the reference formulas (15) and (24), respectively. Given
that formula (29) does not explicitly exhibit symmetry or positive-
definiteness, we prove an equivalent reformulation, easier to analyze.

######  Lemma 3.

(Reformulation of $\ell_{h}$) For $\mu$ and
$\eta\in\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})$, denote by
$\mathring{\underline{\omega}}_{h}^{\mu}$ and
$\mathring{\underline{\omega}}_{h}^{\eta}$ their associated solutions of
(28a), respectively. Then

$\ell_{h}(\mu,\eta)=(\mathring{\underline{\omega}}_{h}^{\mu},\mathring{\underline{\omega}}_{h}^{\eta})^{\star}.$
(30)

###### Proof.

For $\mu$ and $\eta\in\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})$, let
$(\mathring{\underline{\omega}}_{h}^{\mu},\mathring{\underline{\psi}}_{h}^{\mu})$
and
$(\mathring{\underline{\omega}}_{h}^{\eta},\mathring{\underline{\psi}}_{h}^{\eta})$
be their respective, associated solutions of (28). Let us start with the
definition (29) of $\ell_{h}$:

$\ell_{h}(\mu,\eta)\mathrel{\mathop{:}}=-a_{h}(\mathring{\underline{\psi}}_{h}^{\mu},\underline{\mathcal{H}}_{h}\eta)+(\mathring{\underline{\omega}}_{h}^{\mu},\underline{\mathcal{H}}_{h}\eta)^{\star}.$
(31)

By writing
$\underline{\mathcal{H}}_{h}\eta=(\underline{\mathcal{H}}_{h}\eta-\mathring{\underline{\omega}}_{h}^{\eta})+\mathring{\underline{\omega}}_{h}^{\eta}$,
the first term becomes

$a_{h}(\mathring{\underline{\psi}}_{h}^{\mu},\underline{\mathcal{H}}_{h}\eta)=a_{h}(\mathring{\underline{\psi}}_{h}^{\mu},\underline{\mathcal{H}}_{h}\eta-\mathring{\underline{\omega}}_{h}^{\eta})+a_{h}(\mathring{\underline{\psi}}_{h}^{\mu},\mathring{\underline{\omega}}_{h}^{\eta}).$
(32)

As $\underline{\mathcal{H}}_{h}$ does not modify the boundary face unknowns,
we have $\underline{\mathcal{H}}_{h}\eta\in\underline{U}_{h,\eta}$.
Additionally,
$\mathring{\underline{\omega}}_{h}^{\eta}\in\underline{U}_{h,\eta}$ by
definition, so the difference
$(\underline{\mathcal{H}}_{h}\eta-\mathring{\underline{\omega}}_{h}^{\eta})\in\underline{U}_{h,0}$.
Consequently, as $\mathring{\underline{\psi}}_{h}^{\mu}$ verifies (28b), it
holds that

$a_{h}(\mathring{\underline{\psi}}_{h}^{\mu},\underline{\mathcal{H}}_{h}\eta-\mathring{\underline{\omega}}_{h}^{\eta})=(\mathring{\underline{\omega}}_{h}^{\mu},\underline{\mathcal{H}}_{h}\eta-\mathring{\underline{\omega}}_{h}^{\eta})^{\star}=(\mathring{\underline{\omega}}_{h}^{\mu},\underline{\mathcal{H}}_{h}\eta)^{\star}-(\mathring{\underline{\omega}}_{h}^{\mu},\mathring{\underline{\omega}}_{h}^{\eta})^{\star}.$
(33)

For the second term of (32), given that
$\mathring{\underline{\omega}}_{h}^{\eta}$ is solution of (28a) and
$\mathring{\underline{\psi}}_{h}^{\mu}\in\underline{U}_{h,0}$, we can use
$\mathring{\underline{\psi}}_{h}^{\mu}$ as a test function in (28a) to show,
since $a_{h}$ is symmetric, that

$a_{h}(\mathring{\underline{\psi}}_{h}^{\mu},\mathring{\underline{\omega}}_{h}^{\eta})=0.$
(34)

Now, plugging (33) and (34) into (32) gives

$a_{h}(\mathring{\underline{\psi}}_{h}^{\mu},\underline{\mathcal{H}}_{h}\eta)=(\mathring{\underline{\omega}}_{h}^{\mu},\underline{\mathcal{H}}_{h}\eta)^{\star}-(\mathring{\underline{\omega}}_{h}^{\mu},\mathring{\underline{\omega}}_{h}^{\eta})^{\star}.$
(35)

Finally, plugging (35) into (31) yields (30). ∎

Using formulation (30), we clearly have

######  Theorem 1.

$\ell_{h}$, defined by (29), is symmetric and positive-definite.

###### Proof.

Formulation (30) clearly shows symmetry. Additionally, for all
$\mu\in\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})$, we have

$\ell_{h}(\mu,\mu)=(\mathring{\underline{\omega}}_{h}^{\mu},\mathring{\underline{\omega}}_{h}^{\mu})^{\star}.$
(36)

Since $(\cdot,\cdot)^{\star}$ is positive semi-definite, so is $\ell_{h}$.
Now, assume $\ell_{h}(\mu,\mu)=0$, which, by (36), means
$(\mathring{\underline{\omega}}_{h}^{\mu},\mathring{\underline{\omega}}_{h}^{\mu})^{\star}=0$.
The property of local stability (27) implies that, for all
$T\in\mathcal{T}_{h}^{\mathrm{B}}$, we have
$\mathring{\underline{\omega}}_{T}^{\mu}=0$, i.e.
$\mathring{\omega}_{T}^{\mu}=0$ and
$\mathring{\omega}_{\mathcal{F}_{T}}^{\mu}=0$. In particular, for all
$F\in\mathcal{F}_{h}^{\mathrm{B}}$, $\mu_{F}=\mathring{\omega}_{F}^{\mu}=0$.
This being true for all $F\in\mathcal{F}_{h}^{\mathrm{B}}$, we have $\mu=0$,
which proves that $\ell_{h}$ is positive-definite. ∎

###  5.2 Discrete scheme

Let
$(\underline{\omega}_{h}^{0},\underline{\psi}_{h}^{0})\in\underline{U}_{h,0}\times\underline{U}_{h,\mathsf{g}_{\mathrm{D}}}$
be the solution of

$\displaystyle a_{h}(\underline{\omega}_{h}^{0},\underline{v}_{h})$
$\displaystyle=(\mathsf{f},v_{\mathcal{T}_{h}})$
$\displaystyle\forall\underline{v}_{h}\in\underline{U}_{h,0},$ (37a)
$\displaystyle a_{h}(\underline{\psi}_{h}^{0},\underline{v}_{h})$
$\displaystyle=(\omega_{\mathcal{T}_{h}}^{0},v_{\mathcal{T}_{h}})$
$\displaystyle\forall\underline{v}_{h}\in\underline{U}_{h,0}.$ (37b)

The discrete formulation of problem (6) reads: find
$\lambda_{\mathcal{F}_{h}^{\mathrm{B}}}\in\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})$
such that

$\ell_{h}(\lambda_{\mathcal{F}_{h}^{\mathrm{B}}},\mu)=\langle\partial_{\mathrm{n},h}(\underline{\psi}_{h}^{0})-\mathsf{g}_{\mathrm{N}},\;\mu\rangle\qquad\forall\mu\in\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}}).$
(38)

We also define the linear operator
$\mathcal{L}_{h}\colon\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})\to\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})$
associated to the bilinear form $\ell_{h}$ such that, for all
$\mu\in\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})$,

$\langle\mathcal{L}_{h}\mu,\eta\rangle\mathrel{\mathop{:}}=\ell_{h}(\mu,\eta)\qquad\forall\eta\in\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}}).$
(39)

Since $\ell_{h}$ is positive-definite, problem (38) is well-posed. Its
solution exists and is unique. Once it is solved, we compute
$(\underline{\omega}_{h},\underline{\psi}_{h})\in\underline{U}_{h,\lambda_{\mathcal{F}_{h}^{\mathrm{B}}}}\times\underline{U}_{h,\mathsf{g}_{\mathrm{D}}}$,
solution of

$\displaystyle a_{h}(\underline{\omega}_{h},\underline{v}_{h})$
$\displaystyle=(\mathsf{f},v_{\mathcal{T}_{h}})$
$\displaystyle\forall\underline{v}_{h}\in\underline{U}_{h,0},$ (40a)
$\displaystyle a_{h}(\underline{\psi}_{h},\underline{v}_{h})$
$\displaystyle=(\omega_{\mathcal{T}_{h}},v_{\mathcal{T}_{h}})$
$\displaystyle\forall\underline{v}_{h}\in\underline{U}_{h,0}.$ (40b)

Finally, the discrete approximation of $(\omega,\psi)$ is given by

$(\omega_{h},\psi_{h})\mathrel{\mathop{:}}=(p_{h}^{k+1}\underline{\omega}_{h},p_{h}^{k+1}\underline{\psi}_{h}).$

Note that while $(\cdot,\cdot)^{\star}$ is required within $\ell_{h}$ for
problem (38) to be well-posed, it has to be used neither for the enforcement
of the source functions in (37) and (40), nor for the computation of
$\partial_{\mathrm{n},h}(\underline{\psi}_{h}^{0})$ (via formula (24)) in the
right-hand side of (38).

##  6 Preconditioner

In this section, we introduce the preconditioned iterative strategy for the
solution of the algebraic realization of (38).

###  6.1 Algebraic setting

For a fixed basis for the discrete HHO spaces, let
$N_{\mathrm{B}}\mathrel{\mathop{:}}=\operatorname{dim}\left(\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})\right)$.
Introduce the algebraic counterpart
$\mathbf{L}_{h}\colon\mathbb{R}^{N_{\mathrm{B}}}\to\mathbb{R}^{N_{\mathrm{B}}}$
of the linear operator $\mathcal{L}_{h}$ defined by (39). We recall that
$\mathbf{L}_{h}$ can be viewed as a matrix in
$\mathbb{R}^{N_{\mathrm{B}}\times N_{\mathrm{B}}}$, whose coefficients are not
explicitly known, but whose application to a vector can be computed. The
algebraic realization of problem (38) reads: find
$\boldsymbol{\lambda}\in\mathbb{R}^{N_{\mathrm{B}}}$ such that

$\mathbf{L}_{h}(\boldsymbol{\lambda})=\mathbf{b},$ (41)

where $\mathbf{b}\in\mathbb{R}^{N_{\mathrm{B}}}$ is an algebraic
representation of
$\partial_{\mathrm{n},h}(\underline{\psi}_{h}^{0})-\mathsf{g}_{\mathrm{N}}$ in
the chosen basis.

###  6.2 An approximate, sparse matrix

The implicit system (41) is solved using a preconditioned, flexible conjugate
gradient (PFCG) algorithm. This section describes the construction of the
preconditioner. The use of a flexible version of the Krylov method is
justified by the preconditioner being non-symmetric [1, 3]. Refer to Section
6.4 for a discussion about the symmetric property of the preconditioner.

Introducing $\left(\mathbf{e}_{j}\right)_{j=1,\dots,N_{\mathrm{B}}}$ the
canonical basis of $\mathbb{R}^{N_{\mathrm{B}}}$, the $j^{\text{th}}$ column
of $\mathbf{L}_{h}$ is given by the evaluation of
$\mathbf{L}_{h}(\mathbf{e}_{j})$. Note that $\mathbf{L}_{h}(\mathbf{e}_{j})$
is a dense vector. Consequently, the explicit computation of $\mathbf{L}_{h}$
would result in a dense matrix. The preconditioner proposed here consists in
the explicit computation of a _sparse_ matrix $\widetilde{\mathbf{L}}_{h}$
approximating $\mathbf{L}_{h}$. In a nutshell, while the computation of each
column $\mathbf{L}_{h}(\mathbf{e}_{j})$ involves solving two Laplacian
problems in the whole domain, our approximation
$\widetilde{\mathbf{L}}_{h}(\mathbf{e}_{j})$ consists in solving those
problems in a restricted neighbourhood of the $j^{th}$ DoF.

Let $j\in\\{1,\dots,N_{\mathrm{B}}\\}$ a fixed unknown index. Let
$F^{j}\in\mathcal{F}_{h}^{\mathrm{B}}$ the face supporting the DoF associated
to $j$, and $T^{j}\in\mathcal{T}_{h}$ the unique element owning $F^{j}$.
Define $\mathcal{T}_{h}^{j}\subset\mathcal{T}_{h}$ a neighbourhood of $T^{j}$:
roughly, a set of elements around $T^{j}$, including $T^{j}$. Let $\Omega^{j}$
be the open, connected set such that
$\overline{\Omega^{j}}\mathrel{\mathop{:}}=\cup_{T\in\mathcal{T}_{h}^{j}}\overline{T}$.
Figure 1(a) illustrates (in grey) such a neighbourhood, where $F^{j}$ is
represented by a thick line and $T^{j}$ by a darker triangle. Relatively to
$\Omega^{j}$, one defines the sets of interior and boundary faces by
$(\mathcal{F}_{h}^{\mathrm{I}})^{j}\mathrel{\mathop{:}}=\\{F\in\mathcal{F}_{h}^{\mathrm{I}}\mid
F\subset\Omega^{j}\\}$ and
$(\mathcal{F}_{h}^{\mathrm{B}})^{j}\mathrel{\mathop{:}}=\\{F\in\mathcal{F}_{h}\mid
F\subset\partial\Omega^{j}\\}$, respectively.

Consider the operator $\widetilde{\mathcal{L}}_{h}^{j}$, counterpart of
$\mathcal{L}_{h}$ in $\Omega^{j}$, in which
$(\mathcal{T}_{h},\mathcal{F}_{h}^{\mathrm{I}},\mathcal{F}_{h}^{\mathrm{B}})$
is replaced with
$(\mathcal{T}_{h}^{j},(\mathcal{F}_{h}^{\mathrm{I}})^{j},(\mathcal{F}_{h}^{\mathrm{B}})^{j})$
for the solution of the Laplacian subproblems, and
$\mathcal{F}_{h}^{\mathrm{B}}$ is replaced with
$\mathcal{F}_{h}^{\mathrm{B}}\cap(\mathcal{F}_{h}^{\mathrm{B}})^{j}$ for the
computation of the normal derivative. This latter point indicates that the
normal derivative is computed on $\partial\Omega^{j}\cap\partial\Omega$ only.
For both subproblems, homogeneous Dirichlet conditions are enforced on the
part of $\partial\Omega^{j}$ interior to $\Omega$. Define
$N_{\mathrm{B}}^{j}\mathrel{\mathop{:}}=\operatorname{dim}\mathbb{P}^{k}((\mathcal{F}_{h}^{\mathrm{B}})^{j})$,
$\left(\widetilde{\mathbf{e}}_{j}\right)_{j=1,\dots,N_{\mathrm{B}}^{j}}$ the
canonical basis of $\mathbb{R}^{N_{\mathrm{B}}^{j}}$, and
$\widetilde{\mathbf{L}}_{h}^{j}\colon\mathbb{R}^{N_{\mathrm{B}}^{j}}\to\mathbb{R}^{N_{\mathrm{B}}^{j}}$
the algebraic representation of $\widetilde{\mathcal{L}}_{h}^{j}$.

Supposing that the $j^{\text{th}}$ basis function of
$\mathbb{P}^{k}(\mathcal{F}_{h}^{\mathrm{B}})$ corresponds to the
$m^{\text{th}}$ basis function of
$\mathbb{P}^{k}((\mathcal{F}_{h}^{\mathrm{B}})^{j})$, we compute
$\widetilde{\mathbf{L}}_{h}^{j}(\widetilde{\mathbf{e}}_{m})$. As
$\widetilde{\mathbf{L}}_{h}^{j}(\widetilde{\mathbf{e}}_{m})$ yields the normal
derivative only on $\partial\Omega^{j}\cap\partial\Omega$, we build the final,
approximate column $\widetilde{\mathbf{L}}_{h}(\mathbf{e}_{j})$ from
$\widetilde{\mathbf{L}}_{h}^{j}(\widetilde{\mathbf{e}}_{m})$ by setting zero
coefficients where the domain boundary is not covered by the neighbourhood
boundary.

We claim that $\widetilde{\mathbf{L}}_{h}(\mathbf{e}_{j})$ yields a sparse
approximation of $\mathbf{L}_{h}(\mathbf{e}_{j})$. To support our claim, in
Figures 1(b) and 1(c) we plot, in iso-value curves, the solutions to the first
and second Laplacian problems computed in $\Omega^{j}$. Moreover, Figures 1(d)
and 1(e) plot the solutions of the Laplacian subproblems computed in the whole
domain. We observe that, owing to the homogeneity of the first equation and
the homogeneous Dirichlet condition everywhere except on one face, the
solution culminates on that face and decreases towards zero as it goes further
away. Roughly speaking, our strategy consists in approximating the solution
locally in $\Omega^{j}$ while imposing it to be zero in the rest of the
domain. Then, this first truncated solution is used as the source function in
the second problem, also solved in $\Omega^{j}$. Here, note that the quantity
of interest is not the solution itself, but its normal derivative on the
boundary. In particular, the neighbourhood of the considered face concentrates
the most information, insofar as the solution flattens as it goes away from
it. In that sense, $\widetilde{\mathbf{L}}_{h}(\mathbf{e}_{j})$ yields a
reasonable approximation of $\mathbf{L}_{h}(\mathbf{e}_{j})$.

(a) Edge $F^{j}$ (bold segment), element $T^{j}$ (dark grey) and neighbourhood
$\Omega^{j}$ (light grey)

(b) Solution of the first Laplacian problem computed in the neighbourhood

(c) Solution of the second Laplacian problem computed in the neighbourhood

(d) Solution of the first Laplacian problem computed in the whole domain

(e) Solution of the second Laplacian problem computed in the whole domain

Figure 1: Illustration of the action of the preconditioner.

###  6.3 Computation in practice

Regarding the actual computation of $\widetilde{\mathbf{L}}_{h}$, one can note
that the HHO matrix blocks used to solve the Laplacian problems in
$\Omega^{j}$ and to compute the normal derivative can be extracted from the
global ones. The sparsity of $\widetilde{\mathbf{L}}_{h}$ is controlled by the
number of boundary faces in the chosen neighbourhood, while the accuracy of
the approximation depends on the size of the neighbourhood. In particular,
remark that choosing $\mathcal{T}_{h}^{j}\mathrel{\mathop{:}}=\mathcal{T}_{h}$
for all $j$ yields the actual matrix $\mathbf{L}_{h}$. In practice, we
construct $\mathcal{T}_{h}^{j}$ by adding successive layers of neighbours
around $T^{j}$. The number of layers is defined by the parameter
$\alpha\in\mathbb{N}_{0}$. $\mathcal{T}_{h}^{j}$ is defined as
$\mathcal{T}_{h}^{j}(\alpha)$, where
$\mathcal{T}_{h}^{j}(0)\mathrel{\mathop{:}}=\\{T^{j}\\}$ and
$\mathcal{T}_{h}^{j}(\alpha)\mathrel{\mathop{:}}=\mathcal{T}_{h}^{j}(\alpha-1)\cup\\{T\in\mathcal{T}_{h}\mid\exists
T^{\prime}\in\mathcal{T}_{h}^{j}(\alpha-1)\text{ s.t.
}\overline{T}\cap\overline{T^{\prime}}\not=\emptyset\\}$ for $\alpha\geq 1$.
Note that this definition understands the neighbouring relationship as having
at least one vertex in common. The neighbourhood represented in Figure 1(a)
corresponds to $\alpha=3$.

###  6.4 Discussion on symmetry

Choosing a constant $\alpha$ for all $j$ allows to preserve a symmetric
sparsity pattern. Indeed, considering
$(i,j)\in\\{1,\dots,N_{\mathrm{B}}\\}^{2}$,
$(\widetilde{\mathbf{L}}_{h})_{ij}\not=0$ implies that
$F^{i}\subset\partial\Omega^{j}$. Consequently,
$(\widetilde{\mathbf{L}}_{h})_{ji}\not=0$ implies that, reciprocally,
$F^{j}\subset\partial\Omega^{i}$. However, symmetry itself is not preserved,
inasmuch as $(\widetilde{\mathbf{L}}_{h})_{ij}$ results from computations in
the neighbourhood $\Omega^{j}$ whereas $(\widetilde{\mathbf{L}}_{h})_{ji}$
results from computations in $\Omega^{i}\not=\Omega^{j}$.

To enforce the symmetry of the matrix $\mathbf{L}_{h}$, one must ensure that
for all couple $(i,j)\in\\{1,\dots,N_{\mathrm{B}}\\}^{2}$,
$\Omega^{i}=\Omega^{j}$. This requirement leads to a patch-based method: the
set of boundary faces must be split into _patches_ , i.e. connected subsets of
boundary faces. For each patch, one single neighbourhood must be defined to
process all DoFs included in that patch. This yields a block-diagonal matrix
$\widetilde{\mathbf{L}}_{h}$, one block corresponding to each patch.
Algebraically, such a preconditioner is in fact an approximate block Jacobi
preconditioner, insofar as each diagonal block approximates the corresponding
diagonal block of the exact matrix $\mathbf{L}_{h}$. This symmetric
preconditioner was numerically tested. Even with the facts that a non-flexible
Krylov method can be used and that the matrix $\widetilde{\mathbf{L}}_{h}$ is
easier to factorize, the convergence of the method is significantly slower
than with the unsymmetric preconditioner defined in Section 6.3. Therefore,
only the latter has been used in the numerical experiments of Section 7.

##  7 Numerical experiments

###  7.1 Experimental setup

The implicit system (41) is solved iteratively by the PFCG method, using the
preconditioning technique described in Section 6. Given the iterate
$\widetilde{\boldsymbol{\lambda}}\in\mathbb{R}^{N_{\mathrm{B}}}$, the
corresponding residual vector is defined as
$\mathbf{r}\mathrel{\mathop{:}}=\mathbf{b}-\mathbf{L}_{h}(\widetilde{\boldsymbol{\lambda}})$.
The PFCG algorithm stops when
$||\mathbf{r}||_{2}/||\mathbf{b}||_{2}<\varepsilon$, where $\varepsilon>0$ is
a fixed tolerance and $||\cdot||_{2}$ denotes the Euclidean norm on
$\mathbb{R}^{N_{\mathrm{B}}}$. All experiments are conducted choosing $l=k$,
but we stress that the same qualitative results are obtained with $l=k+1$. In
2D, the linear systems corresponding to the Laplacian subproblems are solved
by Cholesky factorization; in 3D, we use a $p$-multigrid algorithm on top of
the algebraic multigrid method designed in [11]. In that latter case, the same
tolerance $\varepsilon$ must be used to stop the Laplacian solvers and the
global PFCG algorithm. To apply the preconditioner, we solve
$\widetilde{\mathbf{L}}_{h}$ using the BiCGSTAB algorithm with tolerance
$\varepsilon$. The computations are run on an 8-core processor (AMR M1 Pro)
clocked at 3228 MHz.

###  7.2 Preliminary: convergence of the discrete normal derivative

We first assert the validity of the discrete normal derivative proposed in
Section 4 by evaluating its order of convergence. Let
$\Omega\mathrel{\mathop{:}}=(0,1)^{2}$ and $u\colon(x,y)\mapsto\sin(4\pi
x)\sin(4\pi y)$ the manufactured solution of the boundary value problem (8),
where $f$ and $g_{\mathrm{D}}$ are defined accordingly. The problem is
discretized on a sequence of successively refined meshes, and we consider
$\partial_{\mathrm{n},h}(\underline{u}_{h})$ computed by (25). The exact
normal derivative $\partial_{\mathrm{n}}u$ is known, and we assess the
relative $L^{2}$-error
$\|\partial_{\mathrm{n}}u-\partial_{\mathrm{n},h}(\underline{u}_{h})\|_{L^{2}(\partial\Omega)}/\|\partial_{\mathrm{n}}u\|_{L^{2}(\partial\Omega)}$
for $k\in\\{0,1,2,3\\}$. Figures 2(a) and 2(b) present the experimental
results on uniform Cartesian meshes and on unstructured, triangular Delaunay
meshes, respectively. These experiments assess a convergence in
$\mathcal{O}(h^{k+1})$ (also observed for $l=k+1$), which is in line with the
literature on finite elements [21, Cor. 6.1]. Indeed, the convergence order of
the normal derivative relies on that of the Laplacian scheme: if the discrete
solution of the Laplacian scheme converges in $\mathcal{O}(h^{p})$, with
$p>1$, then the discrete normal derivative on the boundary converges in
$\mathcal{O}(h^{p-1})$. Transposed to HHO, since
$u_{h}\mathrel{\mathop{:}}=p_{h}^{k+1}\underline{u}_{h}$ converges in
$\mathcal{O}(h^{k+2})$, it is then expected that
$\partial_{\mathrm{n},h}(\underline{u}_{h})$ converge in
$\mathcal{O}(h^{k+1})$. Note that in the case of a non-convex domain with re-
entrant corners, this convergence order is expected to be reduced due to the
corner singularities [24]. One way to mitigate this issue is the use of graded
meshes (cf. [24]).

$10^{-2}$$10^{-1}$$10^{-11}$$10^{-6}$$10^{-1}$1234$h$$L^{2}$-error$k=0$$k=1$$k=2$$k=3$
(a) Uniform, Cartesian meshes.

$10^{-2}$$10^{-1}$$10^{-11}$$10^{-6}$$10^{-1}$1234$h$ (b) Delaunay meshes.

Figure 2: Convergence of the discrete normal derivative
$\partial_{\mathrm{n},h}(\underline{u}_{h})$ for a smooth solution.

###  7.3 Convergence rate of the scheme

On the unit square, we consider $\psi(x,y)=x\sin(\pi y)e^{-xy}$ and its
corresponding $\omega=-\Delta\psi$, the manufactured solution of (2), with
$\mathsf{f},\mathsf{g}_{\mathrm{D}}$ and $\mathsf{g}_{\mathrm{N}}$ determined
accordingly. The square is decomposed into a sequence of successively finer
meshes, and the numerical scheme is applied for $k\in\\{0,1,2,3\\}$. Figure 3
presents the evolution of the relative $L^{2}$-error achieved by the
approximation $(\omega_{h},\psi_{h})$ with respect to $(\omega,\psi)$, on a
sequence of Cartesian meshes. Similarly, Figure 4 shows the analogous results
obtained with a sequence of polygonal meshes, each constructed from an initial
Cartesian mesh by agglomeration and face collapsing. One such mesh is
represented in Figure 6. We empirically interpret these results in light of
Glowinski _et al._ ’s theoretical findings in [19] for the continuous high-
order FEM. According to [19, Section 3.3], since $\psi$ is regular enough, the
FEM solution $(\omega_{h}^{\textrm{FEM}},\psi_{h}^{\textrm{FEM}})$ of order
$p\geq 3$ verifies

$\|\psi_{h}^{\textrm{FEM}}-\psi\|_{L^{2}(\Omega)}+h^{2}\|\omega_{h}^{\textrm{FEM}}-\omega\|_{L^{2}(\Omega)}\leq
Ch^{p+1}\|\psi\|_{H^{p+1}(\Omega)},$ (42)

where $C$ is a constant independent of $h$ and $\psi$. In our HHO setting, the
above error estimate would translate in the following one: if $k\geq 2$, then

$\|\psi_{h}-\psi\|_{L^{2}(\Omega)}+h^{2}\|\omega_{h}-\omega\|_{L^{2}(\Omega)}\leq
Ch^{k+2}\|\psi\|_{H^{k+2}(\Omega)}.$ (43)

The experimental results of Figure 3 validate the estimate (43), not only for
$k\geq 2$, but also for $k\in\\{0,1\\}$. More precisely, while Figure 3(b)
shows that the estimate is sharp for $\psi$ (i.e. $\mathcal{O}(h^{k+2})$),
Figure 3(a) shows, for $\omega$, a faster convergence than indicated by (43).
Namely, the convergence seems to be $\mathcal{O}(h^{k+1})$ for $k=0$ and
$\mathcal{O}(h^{k+\nicefrac{{1}}{{2}}})$ for $k\geq 1$. The experiments of
Figure 5 on polygonal meshes exhibit the same convergence orders, except for
$\omega_{h}$ with $k=0$, where the error estimate reduces to
$\mathcal{O}(h^{k+\nicefrac{{1}}{{2}}})$, in accordance with the other values
of $k$. Based on these experiments, we conjecture that the following error
estimate is sharp for all $k\geq 0$:

$\|\psi_{h}-\psi\|_{L^{2}(\Omega)}+h^{\nicefrac{{3}}{{2}}}\|\omega_{h}-\omega\|_{L^{2}(\Omega)}\leq
Ch^{k+2}\|\psi\|_{H^{k+2}(\Omega)}.$ (44)

Nonetheless, we keep in mind that $k=0$ seems to be a special case for
$\omega_{h}$, where superconvergence may be observed: besides the order 1
obtained on Cartesian meshes (Figure 3(a)), we observe the order 2 if we
choose as an exact solution the polynomial function
$\psi(x,y)=x^{4}(x-1)^{2}y^{4}(y-1)^{2}$, as illustrated in Figure 5.

######  Remark 4.

(Comparison with [17]) The HHO discretizations of [17] rely on the primal
formulation (1) of the equation. They hinge on the local approximation space
$\mathbb{P}^{k+2}(T)\times\mathbb{P}^{k+1}(\mathcal{F}_{T})\times\mathbb{P}^{k}(\mathcal{F}_{T})$
in 2D (resp.
$\mathbb{P}^{k+2}(T)\times\mathbb{P}^{k+2}(\mathcal{F}_{T})\times\mathbb{P}^{k}(\mathcal{F}_{T})$
in 3D), where the associated set of DoFs aims at approximating
$(\psi_{|T},\psi_{|\partial T},\nabla\psi\cdot\mathbf{n}_{\partial T})$. With
this setting, the methods achieve a convergence order of $k+3$ in
$L^{2}$-norm. In comparison, the present method requires the space
$\mathbb{P}^{k+1}(T)\times\mathbb{P}^{k+1}(\mathcal{F}_{T})$ to obtain the
same convergence order, i.e. one less polynomial degree for the cell unknowns
and no DoF dedicated to the approximation of
$\nabla\psi\cdot\mathbf{n}_{\partial T}$. However, if our method is
structurally lighter in DoFs than [17], one cannot conclude regarding the
computational cost and overall efficiency. Attempting a fair comparison
necessitates to account for the fact that, in spite of the lower number of
unknowns and the existence of fast solvers for second-order elliptic systems,
two such systems have to be solved at each iteration of the iterative process.
In comparison, the monolithic system arising from [17], although larger and
exhibiting a more challenging structure, has to be solved only once. If an
iterative solver is used, notwithstanding the preconditioner, no inner system
has to be solved at each iteration, which suggests a lower cost per iteration
than the present method. However, for large-scale problems, scalability
becomes essential: on top of the cost per iteration, convergence speed and
massive parallelization become crucial criteria. In this context, the
efficiency of a method strongly depends on the availability of an efficient
preconditioner.

$10^{-2}$$10^{-1}$$10^{-3}$$10^{-6}$$10^{-9}$1.01.52.53.5$h$$L^{2}$-error$k=0$$k=1$$k=2$$k=3$
(a) $\omega$

$10^{-2}$$10^{-1}$$10^{-3}$$10^{-6}$$10^{-9}$$10^{-12}$2345$h$ (b) $\psi$

Figure 3: Convergence in $L^{2}$-norm of the discrete solution
$(\omega_{h},\psi_{h})$ with respect to the exact solution
$\psi(x,y)=x\sin(\pi y)e^{-xy}$. Square domain discretized by Cartesian
meshes.

$10^{-2}$$10^{-1}$$10^{-3}$$10^{-6}$$10^{-9}$0.51.52.53.5$h$$L^{2}$-error$k=0$$k=1$$k=2$$k=3$
(a) $\omega$

$10^{-2}$$10^{-1}$$10^{-3}$$10^{-6}$$10^{-9}$$10^{-12}$2345$h$ (b) $\psi$

Figure 4: Convergence in $L^{2}$-norm of the discrete solution
$(\omega_{h},\psi_{h})$ with respect to the exact solution
$\psi(x,y)=x\sin(\pi y)e^{-xy}$. Square domain discretized by polygonal
meshes.

$10^{-2}$$10^{-1}$$10^{-2}$$10^{-5}$$10^{-8}$2.01.52.53.5$h$$L^{2}$-error$k=0$$k=1$$k=2$$k=3$
(a) $\omega$

$10^{-2}$$10^{-1}$$10^{0}$$10^{-3}$$10^{-6}$$10^{-9}$$10^{-12}$2345$h$ (b)
$\psi$

Figure 5: Convergence in $L^{2}$-norm of the discrete solution
$(\omega_{h},\psi_{h})$ with respect to the exact solution
$\psi(x,y)=x^{4}(x-1)^{2}y^{4}(y-1)^{2}$. Square domain discretized by
Cartesian meshes. Figure 6: Example of a polygonal mesh of the unit square
$\Omega=(0,1)^{2}$.

###  7.4 Preconditioner

We evaluate in this section the efficiency of the preconditioner designed in
Section 6. We recall that the parameter $\alpha\in\mathbb{N}_{0}$ is an
indicator of the size of the neighbourhoods used for the approximation of the
matrix.

####  7.4.1 Convergence vs. cost: the parameter $\alpha$

The considered test case is the square domain partitioned into $256\times 256$
Cartesian elements, with $k=1$. The tolerance is set to
$\varepsilon=10^{-14}$. Let us assess how the preconditioner improves the
speed of convergence. Figure 7(a) plots the decay of the algebraic residual
$||\mathbf{r}||_{2}/||\mathbf{b}||_{2}$ without preconditioner (in red) and
with preconditioner (in blue), considering increasing values of $\alpha$. One
can confirm that, as expected, the higher the value of $\alpha$, the better
the convergence rate. However, increasing $\alpha$ increases the setup cost of
the preconditioner. Indeed, (i) the Laplacian problems are solved in larger
subdomains, (ii) the approximate matrix $\widetilde{\mathbf{L}}_{h}$ is
denser, which increases the cost of its factorization. Figure 7(b) presents
the computational cost, measured in CPU time, to reach the algebraic solution.
One can see that a trade-off between cost and convergence rate must be made to
achieve optimal performance. If raising $\alpha$ is a good strategy up to
$\alpha=6$, above that threshold, the better convergence rate does not make up
for the cost of the setup phase, hence making the overall CPU time slightly
increase. We stress that other choices of $k$ yield the same qualitative
results.

$0$$20$$40$$60$$80$$100$$10^{-2}$$10^{-6}$$10^{-10}$$10^{-14}$Iterationsresidualunprec.
FCGPFCG($\alpha=2$)PFCG($\alpha=4$)PFCG($\alpha=6$)PFCG($\alpha=8$)PFCG($\alpha=10$)
(a) Decay of the algebraic residual

$2$$4$$6$$8$$10$$0$$50$$100$$150$$\alpha$CPU time (s)iterationssetup (b) CPU
time vs. $\alpha$

Figure 7: (a) compares the speed of convergence of the (P)FCG w.r.t. the
normalized algebraic residual $||\mathbf{r}||_{2}/||\mathbf{b}||_{2}$. (b)
compares the CPU time consumed to reach a tolerance of $10^{-14}$, in function
of $\alpha$. The test case is the square domain partitioned into $256\times
256$ Cartesian elements, $k=1$.

####  7.4.2 PFCG convergence rate: $k$-independence and $h$-dependency
assessment

Problem (2) on the unit square with exact solution $\psi(x,y)=x\sin(\pi
y)e^{-xy}$ is solved on a sequence of 2D grids composed of $N^{2}$ Cartesian
elements, $N\in\\{32,64,128,256,512\\}$. Table 1 presents, for
$k\in\\{0,1,2,3\\}$, the number of PFCG iterations executed to reach the
convergence criterion with $\varepsilon=10^{-8}$. The parameter of the
preconditioner is set to $\alpha=8$. For comparison, the number of iterations
without preconditioning is reported in brackets. Firstly, one can remark that
the use of the preconditioner allows a convergence rate independent of $k$,
while a dependency is observed without preconditioner. Secondly, we
approximate the asymptotic convergence rate with respect to the problem size
by using the number of iterations measured on the two finer meshes. This
yields a dependency in $\mathcal{O}(h^{-0.41})$. Additionally, Table 2 reports
the CPU times of the setup and iteration phases. The setup includes the
Cholesky factorization of the Laplacian matrix and the assembly of the
preconditioning matrix. In brackets, the same information without
preconditioning is displayed. Based on these results, one can see that the
compensation of setup cost by the gain in iteration time is more and more
efficient as the problem grows larger.

The equivalent experiment is performed in 3D (exact solution
$\psi(x,y)=xz\sin(\pi y)e^{-xy}$) with unstructured tetrahedral meshes. As the
preconditioner is costlier in 3D, we choose $\alpha=2$. The number of PFCG
iterations (independent of $k$) is presented in Table 3, which exhibits a
dependency in $\mathcal{O}(h^{-0.45})$. The CPU times of setup and iteration
phases are reported for $k=0$. In brackets, the same information without
preconditioning is displayed.

$N=$ | 32 | 64 | 128 | 256 | 512
---|---|---|---|---|---
$k=0$ | 13 (19) | 19 (25) | 26 (33) | 33 (42) | 44 (56)
$k=1$ | 13 (30) | 19 (36) | 26 (47) | 33 (58) | 44 (75)
$k=2$ | 13 (37) | 19 (47) | 26 (62) | 33 (77) | 44 (101)
$k=3$ | 13 (44) | 19 (58) | 26 (75) | 33 (96) | 44 (122)
Table 1: Test case: $N^{2}$ Cartesian elements, $\varepsilon=10^{-8}$. Number of PFCG($\alpha=8$) iterations. In brackets, number of unpreconditioned FCG iterations. $N=$ | 32 | 64 | 128 | 256 | 512
---|---|---|---|---|---
$k=0$ | setup | 0.4 | (0.0) | 0.9 | (0.0) | 1.9 | (0.0) | 3.9 | (0.0) | 10.7 | (2.0)
iter. | 0.1 | (0.1) | 0.3 | (0.4) | 1.8 | (2.3) | 9.4 | (11.8) | 50.8 | (66.7)
$k=1$ | setup | 1.2 | (0.0) | 2.8 | (0.0) | 5.7 | (0.2) | 13.9 | (1.7) | 41.0 | (17.2)
iter. | 0.2 | (0.5) | 1.1 | (2.2) | 6.3 | (11.5) | 34.1 | (60.1) | 194.9 | (336.3)
$k=2$ | setup | 4.8 | (0.0) | 10.1 | (0.0) | 22.0 | (0.6) | 51.0 | (5.6) | 147.8 | (57.0)
iter. | 0.6 | (1.7) | 3.5 | (8.7) | 20.1 | (48.6) | 109.1 | (251.9) | 646.2 | (1450.4)
$k=3$ | setup | 18.2 | (0.0) | 40.5 | (0.1) | 86.8 | (1.2) | 188.4 | (13.4) | 508.0 | (155.7)
iter. | 1.8 | (6.3) | 10.8 | (33.2) | 61.4 | (177.3) | 319.4 | (1028.9) | 2028.4 | (6271.7)
Table 2: Test case: $N^{2}$ Cartesian elements, $\varepsilon=10^{-8}$. CPU times (in seconds) of the setup and iteration phases. In brackets: analogous quantity without preconditioning. $h\approx$ | $h_{0}$ | $h_{0}/2$ | $h_{0}/4$ | $h_{0}/8$
---|---|---|---|---
Elements | $3373$ | $22869$ | $162167$ | $1224468$
iter. counts | 14 (29) | 20 (40) | 27 (50) | 37 (68)
setup time | 5.0 (0.1) | 21.8 (1.1) | 82.8 (8.6) | 581.3 (64.2)
iter. time | 5.3 (0.8) | 24.4 (7.1) | 41.1 (74.9) | 759.4 (1362.1)
Table 3: Test case: Cubic domain, unstructured tetrahedral mesh,
$\varepsilon=10^{-8}$. Number of PFCG($\alpha=2$) iterations, along with the
setup and iteration CPU times (in seconds) for $k=0$. In brackets: analogous
quantity without preconditioning.

##  8 Conclusion

In this work, we extended the scheme of Glowinski _et al._ to HHO
discretizations, yielding an iterative method for the mixed solution of the
biharmonic equation. Its main advantage lies in the fact that it can be
implemented from an existing diffusion code with limited development costs.
Namely, notwithstanding the PFCG solver, one requires only the implementation
of (29) and the right-hand side of (28b). Additionally, for large problems, it
allows the use of fast, elliptic solvers for an enhanced time to solution.
Future work will focus on the a priori error analysis of the proposed method
and on improving the preconditioner, in order to obtain a convergence rate of
the PFCG that is independent of the mesh size, and achieve a scalable
behaviour.

## References

  * [1] R. Blaheta.  GPCG-generalized preconditioned CG method and its use with non-linear and non-symmetric displacement decomposition preconditioners.  Numer. Linear Algebra Appl., 9(6-7):527–550, 2002.  doi:10.1002/nla.295.
  * [2] 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(2):393–421, 2018.  doi:10.1051/m2an/2017065.
  * [3] H. Bouwmeester, A. Dougherty, and A. V. Knyazev.  Nonsymmetric Preconditioning for Conjugate Gradient and Steepest Descent Methods 1.  Procedia Computer Science, 51:276–285, 2015.  doi:10.1016/j.procs.2015.05.241.
  * [4] F. Chave, D. A. Di Pietro, F. Marche, and F. Pigeonneau.  A Hybrid High-Order method for the Cahn-Hilliard problem in mixed form.  SIAM J. Numer. Anal., 54(3):1873–1898, 2016.  doi:10.1137/15M1041055.
  * [5] P. G. Ciarlet and P.-A. Raviart.  A mixed finite element method for the biharmonic equation.  In Mathematical aspects of finite elements in partial differential equations (Proc. Sympos., Math. Res. Center, Univ. Wisconsin, Madison, Wis., 1974), pages 125–145. Publication No. 33, 1974\.
  * [6] P.G. Ciarlet and R. Glowinski.  Dual iterative techniques for solving a finite element approximation of the biharmonic equation.  Comput. Methods Appl. Mech. Eng., 5(3):277–295, 1975.  doi:10.1016/0045-7825(75)90002-X.
  * [7] 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(3):635–650, 2016.  doi:10.1051/m2an/2015051.
  * [8] D. A. Di Pietro and J. Droniou.  The Hybrid High-Order method for polytopal meshes.  Number 19 in Modeling, Simulation and Application. Springer International Publishing, 2020.  doi:10.1007/978-3-030-37203-3.
  * [9] D. A. Di Pietro and A. Ern.  A hybrid high-order locking-free method for linear elasticity on general meshes.  Comput. Methods Appl. Mech. Eng., 283:1–21, 2015.  doi:10.1016/j.cma.2014.09.009.
  * [10] 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(4):461–472, 2014.  doi:10.1515/cmam-2014-0018.
  * [11] D. A. Di Pietro, F. Hülsemann, P. Matalon, P. Mycek, and U. Rüde.  Algebraic multigrid preconditioner for statically condensed systems arising from lowest-order hybrid discretizations.  SIAM J. Sci. Comput., 45(4):S329–S350, 2023.  doi:10.1137/21M1429849.
  * [12] D. A. Di Pietro, F. Hülsemann, P. Matalon, P. Mycek, U. Rüde, and D. Ruiz.  An $h$-multigrid method for Hybrid High-Order discretizations.  SIAM J. Sci. Comput., 43(5):S839–S861, 2021.  doi:10.1137/20M1342471.
  * [13] D. A. Di Pietro, F. Hülsemann, P. Matalon, P. Mycek, U. Rüde, and D. Ruiz.  Towards robust, fast solutions of elliptic equations on complex domains through hybrid high-order discretizations and non-nested multigrid methods.  Int. J. Numer. Methods Eng., 122(22):6576–6595, 2021.  doi:10.1002/nme.6803.
  * [14] D. A. Di Pietro, P. Matalon, P. Mycek, and U. Rüde.  High-order multigrid strategies for hho discretizations of elliptic equations.  Numer. Linear Algebra Appl., 30(1):e2456, 2023.  doi:10.1002/nla.2456.
  * [15] Z. Dong and A. Ern.  Hybrid high-order method for singularly perturbed fourth-order problems on curved domains.  ESAIM: Math. Model. Numer. Anal., 55(6):3091–3114, 2021.  doi:10.1051/m2an/2021081.
  * [16] Z. Dong and A. Ern.  $C^{0}$-hybrid high-order methods for biharmonic problems, 2022.  doi:10.48550/ARXIV.2206.07074.
  * [17] Z. Dong and A. Ern.  Hybrid High-Order and Weak Galerkin Methods for the Biharmonic Problem.  SIAM J. Numer. Anal., 60(5):2626–2656, 2022.  doi:10.1137/21M1408555.
  * [18] Richard S. Falk.  Approximation of the Biharmonic Equation by a Mixed Finite Element Method.  SIAM J. Numer. Anal., 15(3):556–567, 1978.  doi:10.1137/0715036.
  * [19] R. Glowinski and O. Pironneau.  Numerical Methods for the First Biharmonic Equation and for the Two-Dimensional Stokes Problem.  SIAM Review, 21(2):167–212, 1979.  doi:10.1137/1021028.
  * [20] I. Gustafsson.  A Preconditioned Iterative Method for the Solution of the Biharmonic Problem.  IMA Journal of Numerical Analysis, 4(1):55–67, 1984.  doi:10.1093/imanum/4.1.55.
  * [21] J. M. Melenk and B. Wohlmuth.  Quasi-optimal approximation of surface based lagrange multipliers in finite element methods.  SIAM J. Numer. Anal., 50(4):2064–2087, 2012.  doi:10.1137/110832999.
  * [22] Milan D. Mihajlović and David J. Silvester.  Efficient parallel solvers for the biharmonic equation.  Parallel Computing, 30(1):35–55, 2004.  doi:10.1016/j.parco.2003.07.002.
  * [23] P. Peisker.  On the numerical solution of the first biharmonic equation.  ESAIM: Math. Model. Numer. Anal., 22(4):655–676, 1988.  doi:10.1051/m2an/1988220406551.
  * [24] J. Pfefferer and M. Winkler.  Finite element error estimates for normal derivatives on boundary concentrated meshes.  SIAM J. Numer. Anal., 57(5):2043–2073, 2019.  doi:10.1137/18M1181341.
  * [25] A. Quarteroni and A. Valli.  Theory and Application of Steklov-Poincaré Operators for Boundary-Value Problems.  In Renato Spigler, editor, Applied and Industrial Mathematics, pages 179–203. Springer Netherlands, Dordrecht, 1991.  doi:10.1007/978-94-009-1908-2_14.
  * [26] M. B. van Gijzen.  Conjugate Gradient-like solution algorithms for the mixed finite element approximation of the biharmonic equation, applied to plate bending problems.  Comput. Methods Appl. Mech. Eng., 121(1):121–136, 1995.  doi:10.1016/0045-7825(94)00737-8.

##  Appendix A Proof of Lemma 1

Let $v_{\mathcal{F}_{h}}\in\mathbb{P}^{k}(\mathcal{F}_{h})$ and
$\underline{w}_{h}\in\underline{U}_{h}$. For all $T\in\mathcal{T}_{h}$,
$\Theta_{T}v_{\mathcal{F}_{T}}$ is by definition solution of (17). By
linearity of $a_{h}$, combining both terms of (17) yields

$a_{T}(\underline{\Theta}_{T}v_{\mathcal{F}_{T}},(w_{T},0))=0.$

That being true for all $T\in\mathcal{T}_{h}$ proves the global counterpart

$a_{h}(\underline{\Theta}_{h}v_{\mathcal{F}_{h}},(w_{\mathcal{T}_{h}},0))=0.$
(45)

By splitting $\underline{w}_{h}$ into
$(w_{\mathcal{T}_{h}},0)+(0,w_{\mathcal{F}_{h}})$ and by linearity of $a_{h}$,
we have

$a_{h}(\underline{\Theta}_{h}v_{\mathcal{F}_{h}},\underline{w}_{h})=a_{h}(\underline{\Theta}_{h}v_{\mathcal{F}_{h}},(w_{\mathcal{T}_{h}},0))+a_{h}(\underline{\Theta}_{h}v_{\mathcal{F}_{h}},(0,w_{\mathcal{F}_{h}})).$

Property (45) allows to cancel the first term of the right-hand side. Then, by
writing
$(0,w_{\mathcal{F}_{h}})=(\Theta_{\mathcal{T}_{h}}w_{\mathcal{F}_{h}},w_{\mathcal{F}_{h}})-(\Theta_{\mathcal{T}_{h}}w_{\mathcal{F}_{h}},0)=\underline{\Theta}_{h}w_{\mathcal{F}_{h}}-(\Theta_{\mathcal{T}_{h}}w_{\mathcal{F}_{h}},0)$
and using the linearity of $a_{h}$, the second term rewrites

$a_{h}(\underline{\Theta}_{h}v_{\mathcal{F}_{h}},(0,w_{\mathcal{F}_{h}}))=a_{h}(\underline{\Theta}_{h}v_{\mathcal{F}_{h}},\underline{\Theta}_{h}w_{\mathcal{F}_{h}})-a_{h}(\underline{\Theta}_{h}v_{\mathcal{F}_{h}},(\Theta_{\mathcal{T}_{h}}w_{\mathcal{F}_{h}},0)).$
(46)

Using (45) again to cancel the last term of (46) finishes the proof of (20).