# A pressure-robust HHO method for the solution of the incompressible
Navier–Stokes equations on general meshes

Daniel Castanon<EMAIL_ADDRESS>Instituto de
Investigaciones en Matemáticas Aplicadas y en Sistemas, Universidad Nacional
Autónoma de México, Circuito Escolar s/n, Ciudad Universitaria C.P. 04510 Cd.
Mx. (México)  Université Côte d’Azur, CNRS, Inria team Coffee, LJAD, Nice,
France  IMAG, Univ Montpellier, CNRS, Montpellier, France  Daniele A. Di
<EMAIL_ADDRESS>IMAG, Univ Montpellier, CNRS,
Montpellier, France

###### Abstract

In a recent work [11], we have introduced a pressure-robust Hybrid High-Order
method for the numerical solution of the incompressible Navier–Stokes
equations on matching simplicial meshes. Pressure-robust methods are
characterized by error estimates for the velocity that are fully independent
of the pressure. A crucial question was left open in that work, namely whether
the proposed construction could be extended to general polytopal meshes. In
this paper we provide a positive answer to this question. Specifically, we
introduce a novel divergence-preserving velocity reconstruction that hinges on
the solution inside each element of a mixed problem on a subtriangulation,
then use it to design discretizations of the body force and convective terms
that lead to pressure robustness. ††margin:  R#1 and R#2 remark: The proposed
method now works for arbitrary $k\geq 0$ and general meshes An in-depth
theoretical study of the properties of this velocity reconstruction, and their
reverberation on the scheme, is carried out for arbitrary polynomial degrees
$k\geq 0$ and meshes composed of general polytopes. The theoretical
convergence estimates and the pressure robustness of the method are confirmed
by an extensive panel of numerical examples.

Key words: Hybrid High-Order methods, incompressible Navier–Stokes equations,
general meshes, pressure robustness
MSC 2010: 65N08, 65N30, 65N12, 35Q30, 76D05

##  1 Introduction

This paper focuses on numerical approximations of the Navier–Stokes equations
robust with respect to large irrotational body forces. Specifically, we
address a nontrivial question left open in the previous work [11], namely
whether robustness can be achieved on general polyhedral meshes such as the
ones supported by the Hybrid High-Order (HHO) method [22, 16].

Let $\Omega\subset\mathbb{R}^{3}$ denote an open, bounded, simply connected
polyhedral domain with Lipschitz boundary $\partial\Omega$. Let $\nu>0$ be the
kinematic viscosity of the fluid and $\boldsymbol{f}\in L^{2}(\Omega)^{3}$ a
given vector field representing a body force. Setting $\boldsymbol{U}\coloneq
H_{0}^{1}(\Omega)^{3}$ and
$P\coloneq{L}^{2}_{0}({\Omega})=\left\\{q\in{L}^{2}({\Omega}):\int_{\Omega}q=0\right\\}$,
we consider the Navier–Stokes problem: Find
$(\boldsymbol{u},p)\in\boldsymbol{U}\times P$ such that

$\displaystyle\nu\int_{\Omega}\nabla\boldsymbol{u}:\nabla\boldsymbol{v}+\int_{\Omega}((\nabla\times\boldsymbol{u})\times\boldsymbol{u})\cdot\boldsymbol{v}-\int_{\Omega}(\nabla\cdot\boldsymbol{v})p$
$\displaystyle=\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{v}$
$\displaystyle\forall\boldsymbol{v}\in\boldsymbol{U},$ (1a)
$\displaystyle\int_{\Omega}(\nabla\cdot\boldsymbol{u})q$ $\displaystyle=0$
$\displaystyle\forall q\in{L}^{2}({\Omega}).$ (1b)

Above, $\nabla\cdot$ and $\nabla\times$ denote, respectively, the divergence
and curl operators, while $\times$ is the cross product of two vectors. The
convective term in (1a) is expressed in rotational form, so $p$ is here the
Bernoulli pressure, which is related to the kinematic pressure
$p_{\text{kin}}$ by the equation
$p=p_{\text{kin}}+\frac{1}{2}|\boldsymbol{u}|^{2}$.

The domain $\Omega$ being simply connected, we have the following Hodge
decomposition of the body force (see, e.g., [2, Section 4.3]):

$\boldsymbol{f}=\boldsymbol{g}+\lambda\nabla\psi,$ (2)

where $\boldsymbol{g}$ is the curl of a function in $\boldsymbol{H}({\bf
curl};\Omega)$ the tangent trace of which vanishes on $\partial\Omega$,
$\psi\in H^{1}(\Omega)$ is such that $\|\nabla\psi\|_{L^{2}(\Omega)^{3}}=1$,
and $\lambda\in\mathbb{R}^{+}$. It is well-known that, at the continuous
level, the velocity field is entirely determined by the first component in the
decomposition (2). This property, however, does not carry out automatically to
the discrete level. The development of numerical methods that possess this
property, and which are sometimes referred to in the literature as _pressure-
robust_ , has been an active field of research over the last few years; see,
e.g., [26, 42, 41, 33, 1] concerning finite element methods on standard
meshes.

Recently, the mathematical community have become interested in the development
of arbitrary-order approximation methods that support more general meshes than
standard finite elements and which can include, e.g., polyhedral elements and
non-matching interfaces. A representative but by far non exhaustive list of
references concerning incompressible flow problems includes [18, 19, 23, 29,
8, 5, 4, 52]; see also the recent works [10, 12] concerning non-Newtonian
fluids. Pressure-robust variations of the HHO method on matching simplicial
meshes for the Stokes and Navier–Stokes problem have been proposed,
respectively, in [20, 11].

The development of pressure-robust methods on polyhedral meshes is, however, a
challenging task. Some of the first genuinely pressure-robust polyhedral
methods for the Stokes equations have been proposed in [43, 51, 53]. These
methods handle the lowest order case using a velocity reconstruction in
$\boldsymbol{H}({\rm div};\Omega)$ introduced in [13] and relying on
Wachspress (generalized barycentric) coordinates. This approach has two
shortcomings: first, the faces of each (convex) polyhedral element must be
either triangles or parallelograms; second, error estimates for the
approximated velocity would require gradient bounds for the Wachspress
coordinates on an arbitrary convex polyhedron, the derivation of which
remains, to the best of our knowledge, an open problem. Regarding arbitrary-
order methods on general meshes, a pressure-robust Virtual Element method has
been recently proposed in [27] for the Stokes equations. The extension of this
method to the Navier–Stokes equations remains, to the best of our knowledge,
an open problem. A pressure-robust discretization scheme for the full
Navier–Stokes equations has been proposed in [35] based on the staggered
Discontinuous Garlekin method. This method solves for three unknowns (the
pressure, the velocity, and its gradient), thus leading to larger algebraic
systems. Recently, a novel HHO method for which pressure-robustness has been
numerically demonstrated has been proposed in [9]. This method uses a larger
pressure space than the one considered in the present work, and the derivation
of rigorous pressure-robust error estimates is still to be done. An entirely
different approach to pressure-robustness on polyhedral meshes has also been
recently pursued in [49], hinging on the compatibility features of Discrete de
Rham [17, 15] and Virtual Element methods. While this approach leads to a
fully pressure-robust, arbitrary-order method, it is based on a curl-curl
formulation of the viscous term, which does not lend itself naturally to the
treatment of certain standard boundary conditions.

In the present work, we propose a novel fully pressure-robust HHO method for
the Navier–Stokes problem (1) that works in space dimension two and three and
supports general meshes composed of polytopal elements. The cornerstone of the
method is a local divergence-preserving reconstruction of the velocity built
inside each mesh element $T$ by solving a mixed problem inspired by [38, 39,
40]††margin:  R#1 major point: checking and citing related work on a
subtriangulation of $T$; see also [50]. The assumptions made in Section 2.1
for each element $T$ enable us to derive the required continuity and
approximation bounds for this reconstruction. Robustness with respect to large
irrotational body forces is achieved by leveraging the divergence-preserving
velocity reconstruction in the discretisation of both the convective term and
the body force, so that similar properties as the ones discussed in [11,
Section 4.3 and Lemma 7] are obtained for these terms.

The rest of the paper is organised as follows. In Section 2 we introduce the
discrete setting, including mesh assumptions, notation, and the novel
divergence-preserving velocity reconstruction. Section 3 contains the discrete
problem and the main results of the analysis, with particular focus on the
definition and properties of the discrete convective trilinear form. A
complete panel of two-dimensional numerical tests on a variety of polygonal
meshes is provided in Section 4, including a comparison with the standard HHO
scheme of [8].

##  2 Discrete setting

The following exposition focuses on the three-dimensional case $d=3$, the two-
dimensional case $d=2$ being a special instance of the latter as detailed in
Remark 13 below.

###  2.1 Mesh

Following [16, Definition 1.4], we consider a polyhedral mesh defined as a
couple $\mathcal{M}_{h}\coloneq(\mathcal{T}_{h},\mathcal{F}_{h})$, where
$\mathcal{T}_{h}$ is a finite collection of polyhedral elements which we
additionally assume to be convex ††margin:  R#1 & R#2 points: Where the
convexity enters and non-convex elements ( see Remark 5 below on how to relax
this assumption), while $\mathcal{F}_{h}$ is a finite collection of planar
faces $F$. For any mesh element or face
$X\in\mathcal{T}_{h}\cup\mathcal{F}_{h}$, we denote by $|X|$ its Hausdorff
measure and by $h_{X}$ its diameter, so that the meshsize satisfies
$h=\max_{T\in\mathcal{T}_{h}}h_{T}$. Boundary faces lying on $\partial\Omega$
and internal faces contained in $\Omega$ are collected in the sets
$\mathcal{F}_{h}^{{\rm b}}$ and $\mathcal{F}_{h}^{{\rm i}}$, respectively. For
each mesh element $T\in\mathcal{T}_{h}$, we denote by $\mathcal{F}_{T}$ the
set collecting the faces that lie on the boundary $\partial T$ of $T$ and, for
all $F\in\mathcal{F}_{T}$, we denote by $\boldsymbol{n}_{TF}$ the (constant)
unit vector normal to $F$ and pointing out of $T$.

It is assumed that $\mathcal{M}_{h}$ belongs to a regular mesh sequence
$(\mathcal{M}_{h})_{h}$ in the sense of [16, Definition 1.9]. This assumption
entails the existence of a matching simplicial submesh
$\mathfrak{M}_{h}\coloneq(\mathfrak{T}_{h},\mathfrak{F}_{h})$ of
$\mathcal{M}_{h}$ with the following properties: $\mathfrak{T}_{h}$ is a
finite collection of simplicial elements; for any simplex
$\tau\in\mathfrak{T}_{h}$ , there is a unique mesh element
$T\in\mathcal{T}_{h}$ such that $\tau\subset T$; for any simplicial face
$\sigma\in\mathfrak{F}_{h}$ and any mesh face $F\in\mathcal{F}_{h}$ , either
$\sigma\cap F=\emptyset$ or $\sigma\subset F$. For $T\in\mathcal{T}_{h}$, we
define $\mathfrak{T}_{T}$ as the set of all simplices of $\mathfrak{T}_{h}$
contained in $T$ (see Figure 1) and $\mathfrak{F}_{T}^{{\rm i}}$ as the set of
faces of $\mathfrak{F}_{h}$ that lie in the interior of $T$. For
$F\in\mathcal{F}_{h}$, $\mathfrak{F}_{F}$ denotes the set of simplicial faces
$\sigma$ for which $\sigma\subset F$, and we let
$\boldsymbol{n}_{\sigma}\coloneq\boldsymbol{n}_{TF}$, and
$\boldsymbol{n}_{\tau\sigma}\coloneq\boldsymbol{n}_{\sigma}$ for the unique
element $\tau\in\mathfrak{T}_{T}$, $T\in\mathcal{T}_{h}$, which contains
$\sigma$. ††margin:  R#1 minor point: Explaining Fig. 1(b), Justification of
(3) Additional notations for mesh elements and faces are introduced at the
beginning of Section 2.5 and illustrated in Figure 1. For future use, we
notice that,  by [16, Lemma 1.12], mesh regularity implies the existence of an
integer $N\geq 0$ depending only on the mesh regularity parameter such that

$\displaystyle\max_{h}\max_{T\in\mathcal{T}_{h}}\operatorname{card}(\mathfrak{T}_{T})\leq
N\qquad\text{and}\qquad\max_{h}\max_{T\in\mathcal{T}_{h}}\operatorname{card}(\mathcal{F}_{T})\leq
N.$ (3)

$F_{1}$
---
$F_{2}$
---
$F_{3}$
---
$F_{4}$
---
$F_{5}$
---
$F_{6}$
---
$\tau_{1}$
---
$\tau_{2}$
---
$\tau_{3}$
---
$\tau_{4}$
---
$\tau_{5}$
---
$\tau_{6}$
---
(a) The elements of $\mathfrak{T}_{T}$ and $\mathcal{F}_{T}$.

$\tau_{1}$
---
$\sigma_{2}$
---
$\sigma_{1}$
---
$\sigma_{3}$
---
$\sigma_{4}$
---
$\sigma_{5}$
---
$F_{1}$
---
$F_{2}$
---
$\mathbf{n}_{\sigma_{2}}$
---
$\tau_{6}$
---
$\mathbf{n}_{\sigma_{3}}$
---
$\mathbf{n}_{\sigma_{1}}$
---
$\mathbf{n}_{TF_{2}}$
---
$\mathbf{n}_{TF_{1}}$
---
(b) A closer look to the bottom part: The faces
$\sigma_{1},\sigma_{2},\sigma_{3}$ are interior faces, i.e.,
$\\{\sigma_{1},\sigma_{2},\sigma_{3}\\}\subset\mathfrak{F}_{T}^{{\rm i}}$. For
the set $\\{\sigma_{4},\sigma_{5}\\}$, we have $\sigma_{4}=F_{1}$ and
$\sigma_{5}=F_{2}$.

Figure 1: An illustration of the sets $\mathfrak{T}_{T},\mathcal{F}_{T}$ and
$\mathfrak{F}_{T}^{{\rm i}}$ for a given element $T\in\mathcal{T}_{h}$ in
$\mathbb{R}^{2}$.

$\boldsymbol{x}_{T}$
---
(a) Pyramidal submesh.

$\boldsymbol{x}_{T}$
---
(b) Non-pyramidal submesh.

Figure 2: Two examples of submeshes $\mathfrak{T}_{T}$ in $\mathbb{R}^{2}$
that satisfy the assumptions of the Section 2.1 . The red dot represents the
vertex $\boldsymbol{x}_{T}$. ††margin:  Referees #1 and #2 remark: For the
general case $k\geq 2$, we need to make an additional assumption

We additionally make the assumption that, for all element
$T\in\mathcal{T}_{h}$, its submesh $\mathfrak{T}_{T}$ is constructed in such
way that all simplices in $\mathfrak{T}_{T}$ have at least one common vertex
(see Remarks 4 and 16 for the technical details of this assumption) . This
vertex will be denoted ${\boldsymbol{x}}_{T}$. In particular, when
$\boldsymbol{x}_{T}$ lies in the interior of $T$, we call $\mathfrak{T}_{T}$ a
pyramidal submesh.  The Figure 2 shows two examples of submeshes that satisfy
the current assumption.

In order to prevent the proliferation of generic constants we write, whenever
possible, $a\lesssim b$ in place of $a\leq Cb$ with $C>0$ independent of
$\nu$, $\lambda$, $h$ and, for local inequalities, also on the mesh element or
face. The dependencies of the hidden constant will be further specified when
relevant. Moreover, we write $a\simeq b$, when both $a\lesssim b$ and
$b\lesssim a$ hold.

###  2.2 Local and broken spaces and projectors

Let $X$ denote a mesh element or face and, for an integer $l\geq 0$, denote by
$\mathcal{P}^{l}(X)$ the space spanned by the restrictions to $X$ of
polynomials in the space variables of total degree $\leq l$. The
$L^{2}$-orthogonal projector
$\pi_{X}^{l}:L^{1}(X)\rightarrow\mathcal{P}^{l}(X)$ is such that, for all
$\zeta\in L^{1}(X)$,

$\int_{X}(\zeta-\pi_{X}^{l}\zeta)w=0\qquad\forall w\in\mathcal{P}^{l}(X).$ (4)

Vector and matrix versions of the $L^{2}$-orthogonal projector are obtained by
applying $\pi_{X}^{l}$ component-wise, and are both denoted with the bold
symbol ${\boldsymbol{\pi}}_{X}^{l}$ in what follows. Optimal approximation
properties for the $L^{2}$-orthogonal projector are proved in [21, Appendix
A.2]; see also [16, Chapter 1], where these estimates are extended to non-star
shaped elements. Specifically, let $s\in\\{0,\dots,l+1\\}$ and
$r\in[1,+\infty]$. Then, it holds, with hidden constant only depending on $l$,
$s$, $r$, and the mesh regularity parameter: For all
$T\in{{\mathcal{T}}_{h}}$, all $\zeta\in W^{s,r}(T)$, and all
$m\in\\{0,\dots,s\\}$,

$|\zeta-\pi_{T}^{l}\zeta|_{W^{m,r}(T)}\lesssim
h_{T}^{s-m}|\zeta|_{W^{s,r}(T)},$ (5a) and, if $s\geq 1$ and $m\leq s-1$,
$h_{T}^{\frac{1}{r}}|\zeta-\pi_{T}^{l}\zeta|_{W^{m,r}({\mathcal{F}}_{T})}\lesssim
h_{T}^{s-m}|\zeta|_{W^{s,r}(T)},$ (5b)

where $W^{m,r}({\mathcal{F}}_{T})$ is the space spanned by functions in
$L^{r}(\partial T)$ that are in $W^{m,r}(F)$ for all $F\in{\mathcal{F}}_{T}$,
endowed with the corresponding broken norm.

At the global level, the space of broken polynomial functions on
${\mathcal{T}}_{h}$ of total degree $\leq l$ is denoted by
$\mathcal{P}^{l}({\mathcal{T}}_{h})$, and $\pi_{h}^{l}$ is the corresponding
$L^{2}$-orthogonal projector such that, for all $\zeta\in L^{1}(\Omega)$,
$(\pi_{h}^{l}\zeta)_{|T}\coloneq\pi_{T}^{l}\zeta_{|T}$ for all
$T\in{\mathcal{T}}_{h}$. Regularity requirements in error estimates will be
expressed in terms of the broken Sobolev spaces $W^{s,r}({\mathcal{T}}_{h})$
spanned by functions in $L^{r}(\Omega)$ the restriction of which to every
$T\in\mathcal{T}_{h}$ is in $W^{s,r}(T)$. We additionally set, as usual,
$H^{s}({\mathcal{T}}_{h})\coloneq W^{s,2}({\mathcal{T}}_{h})$.

###  2.3 Discrete spaces and norms

Let a polynomial degree $k\geq 0$ be fixed. We define the HHO space as usual,
setting

$\displaystyle\uline{\boldsymbol{U}}_{h}^{k}\coloneq\left\\{\uline{\boldsymbol{v}}_{h}=((\boldsymbol{v}_{T})_{T\in{\mathcal{T}}_{h}},(\boldsymbol{v}_{F})_{F\in{\mathcal{F}}_{h}}):\mbox{$\boldsymbol{v}_{T}\in\mathcal{P}^{k}(T)^{3}$
for all $T\in{\mathcal{T}}_{h}$ and
$\boldsymbol{v}_{F}\in\mathcal{P}^{k}(F)^{3}$ for all
$F\in{\mathcal{F}}_{h}$}\right\\}.$

The restrictions of $\uline{\boldsymbol{U}}_{h}^{k}$ and
$\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}_{h}^{k}$ to a generic
mesh element $T\in{\mathcal{T}}_{h}$ are respectively denoted by
$\uline{\boldsymbol{U}}_{T}^{k}$ and
$\uline{\boldsymbol{v}}_{T}=(\boldsymbol{v}_{T},(\boldsymbol{v}_{F})_{F\in{\mathcal{F}}_{T}})$.
The vector of polynomials corresponding to a smooth function over $\Omega$ is
obtained via the global interpolation operator
$\uline{\boldsymbol{I}}_{h}^{k}:H^{1}(\Omega)^{3}\rightarrow\uline{\boldsymbol{U}}_{h}^{k}$
such that, for all $\boldsymbol{v}\in H^{1}(\Omega)^{3}$,

$\uline{\boldsymbol{I}}_{h}^{k}\boldsymbol{v}\coloneq(({\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}_{|T})_{T\in{\mathcal{T}}_{h}},({\boldsymbol{\pi}}_{F}^{k}\boldsymbol{v}_{|F})_{F\in{\mathcal{F}}_{h}}).$
(6)

Its restriction to a generic mesh element $T\in{\mathcal{T}}_{h}$, collecting
the components on $T$ and its faces, is denoted by
$\uline{\boldsymbol{I}}_{T}^{k}$. We furnish $\uline{\boldsymbol{U}}_{h}^{k}$
with the discrete $H^{1}$-like seminorm such that, for all
$\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}_{h}^{k}$,

$\|\underline{\boldsymbol{v}}_{h}\|_{1,h}\coloneq\left(\sum_{T\in{\mathcal{T}}_{h}}\|\underline{\boldsymbol{v}}_{T}\|_{1,T}^{2}\right)^{\nicefrac{{1}}{{2}}},$

where, for all $T\in{\mathcal{T}}_{h}$,

$\|\underline{\boldsymbol{v}}_{T}\|_{1,T}\coloneq\left(\|\nabla{\boldsymbol{v}}_{T}\|_{L^{2}(T)^{3\times
3}}^{2}+|\uline{\boldsymbol{v}}_{T}|_{1,\partial
T}^{2}\right)^{\nicefrac{{1}}{{2}}}\mbox{ with
}|\uline{\boldsymbol{v}}_{T}|_{1,\partial
T}\coloneq\left(\sum_{F\in{\mathcal{F}}_{T}}\\!h_{F}^{-1}\|\boldsymbol{v}_{F}-\boldsymbol{v}_{T}\|^{2}_{L^{2}(F)^{3}}\right)^{\nicefrac{{1}}{{2}}}.$
(7)

The discrete spaces for the velocity and the pressure, respectively accounting
for the wall boundary condition and the zero-average condition, are

$\uline{\boldsymbol{U}}_{h,0}^{k}\coloneq\left\\{\uline{\boldsymbol{v}}_{h}=((\boldsymbol{v}_{T})_{T\in{\mathcal{T}}_{h}},(\boldsymbol{v}_{F})_{F\in{\mathcal{F}}_{h}})\in\uline{\boldsymbol{U}}_{h}^{k}:\boldsymbol{v}_{F}=\boldsymbol{0}\quad\forall
F\in{\mathcal{F}}_{h}^{\rm b}\right\\},\qquad
P_{h}^{k}\coloneq\mathcal{P}^{k}({\cal{T}}_{h})\cap P.$

For all $\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}_{h}^{k}$, we
denote by ${\boldsymbol{v}}_{h}\in\mathcal{P}^{k}({\mathcal{T}}_{h})^{3}$ the
vector-valued broken polynomial function obtained patching element-based
unkowns, that is $(\boldsymbol{v}_{h})_{|T}\coloneq\boldsymbol{v}_{T}$ for all
$T\in{\mathcal{T}}_{h}$. The following discrete Sobolev embeddings in
$\uline{\boldsymbol{U}}_{h,0}^{k}$ have been proved in [21, Proposition 5.4]:
For all $r\in[1,6]$ it holds, for all
$\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}_{h,0}^{k}$,

$\displaystyle\|{\boldsymbol{v}}_{h}\|_{L^{r}(\Omega)^{3}}\lesssim\|\uline{\boldsymbol{v}}_{h}\|_{1,h}.$
(8)

where the hidden constant is independent of both $h$ and
$\uline{\boldsymbol{v}}_{h}$, but possibly depends on $\Omega$, $k$, $r$, and
the mesh regularity parameter. It follows from (8) that the map
$\|{\cdot}\|_{1,h}$ defines a norm on $\uline{\boldsymbol{U}}_{h,0}^{k}$.
Classically, the corresponding dual norm of a linear form
$\mathcal{L}_{h}:\uline{\boldsymbol{U}}_{h,0}^{k}\to\mathbb{R}$ is given by

$\|\mathcal{L}_{h}\|_{1,h,*}\coloneq\sup_{\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}_{h,0}^{k},\|\uline{\boldsymbol{v}}_{h}\|_{1,h}=1}\left|\mathcal{L}_{h}(\uline{\boldsymbol{v}}_{h})\right|.$
(9)

###  2.4 Divergence-preserving local velocity reconstruction

Following [22], for any element $T\in{\mathcal{T}}_{h}$ we define the discrete
divergence operator
$D_{T}^{k}:\uline{\boldsymbol{U}}_{T}^{k}\rightarrow\mathcal{P}^{k}(T)$ such
that, for all $\uline{\boldsymbol{v}}_{T}\in\uline{\boldsymbol{U}}_{T}^{k}$
and all $q\in\mathcal{P}^{k}(T)$,

$\int_{T}D_{T}^{k}\uline{\boldsymbol{v}}_{T}q=-\int_{T}\boldsymbol{v}_{T}\cdot\nabla
q+\sum_{F\in{\mathcal{F}}_{T}}\int_{F}(\boldsymbol{v}_{F}\cdot\boldsymbol{n}_{TF})q.$
(10)

Crucially, the operator $D_{T}^{k}$ satisfies the following commutation
property (see [16, Eq. (8.21)]):

$D_{T}^{k}\uline{\boldsymbol{I}}_{T}^{k}\boldsymbol{v}=\pi_{T}^{k}(\nabla\cdot\boldsymbol{v})\qquad\forall\boldsymbol{v}\in
H^{1}(T)^{3}.$ (11)

To achieve pressure robustness in the sense made precise by Remark 15 below,
we reconstruct divergence-preserving velocity test functions, which are used
for the discretization of the body force and the nonlinear term. Let an
element $T\in{\mathcal{T}}_{h}$ be fixed and, for $\tau\in\mathfrak{T}_{T}$,
denote by
$\boldsymbol{\mathcal{RT}}^{k}(\tau)\coloneq\mathcal{P}^{k}(\tau)^{3}+\boldsymbol{x}\mathcal{P}^{k}(\tau)$
the local Raviart–Thomas–Nédélec space of degree $k$ [47, 45]. We recall that
a function in $\boldsymbol{\mathcal{RT}}^{k}(\tau)$ is uniquely determined by
its polynomial moments of degree up to $(k-1)$ inside $\tau$ and the
polynomial moments of degree $k$ of its normal component on each face
$\sigma\in\mathfrak{F}_{\tau}$ (with $\mathfrak{F}_{\tau}$ denoting the subset
of $\mathfrak{F}_{h}$ collecting the simplicial faces of $\tau$). We
additionally note the following local norm equivalence uniform in $h$:

$\displaystyle\|\boldsymbol{\mathfrak{w}}\|_{{L}^{2}({\tau})^{3}}^{2}\simeq\|{\boldsymbol{\pi}}_{\tau}^{k-1}\boldsymbol{\mathfrak{w}}\|_{{L}^{2}({\tau})^{3}}^{2}+\sum_{\sigma\in\mathfrak{F}_{\tau}}h_{\sigma}\|\boldsymbol{\mathfrak{w}}\cdot\boldsymbol{n}_{\tau\sigma}\|_{{L}^{2}({\sigma})}^{2}\qquad\forall\boldsymbol{\mathfrak{w}}\in\boldsymbol{\mathcal{RT}}^{k}(\tau).$
(12)

We introduce the Raviart–Thomas–Nédélec space of degree $k$ on the matching
simplicial submesh $\mathfrak{T}_{T}$ of $T$ defined as follows:

$\boldsymbol{\mathcal{RT}}^{k}(\mathfrak{T}_{T})\coloneq\left\\{\boldsymbol{\mathfrak{w}}\in{\boldsymbol{H}}_{\text{div}}({T}):\text{$\boldsymbol{\mathfrak{w}}_{|\tau}\in\boldsymbol{\mathcal{RT}}^{k}(\tau)$
for all $\tau\in\mathfrak{T}_{T}$ }\right\\},$

where
${\boldsymbol{H}}_{\text{div}}({T})\coloneq\\{\boldsymbol{\mathfrak{w}}\in
L^{2}(T)^{3}:\nabla\cdot\boldsymbol{\mathfrak{w}}\in{L}^{2}({T})\\}$. We also
introduce the subspace of $\boldsymbol{\mathcal{RT}}^{k}(\mathfrak{T}_{T})$
spanned by functions with zero normal trace on the boundary of $T$:

$\boldsymbol{\mathcal{RT}}^{k}_{0}(\mathfrak{T}_{T})\coloneq\left\\{\boldsymbol{\mathfrak{w}}\in\boldsymbol{\mathcal{RT}}^{k}(\mathfrak{T}_{T}):\boldsymbol{\mathfrak{w}}\cdot\text{$\boldsymbol{n}_{\sigma}=0$
for all $\sigma\in\mathfrak{F}_{F}$ and all $F\in\mathcal{F}_{T}$}\right\\}.$

Recall from Section 2.1 that, for a given element $T\in\mathcal{T}_{h}$, we
denote by $\boldsymbol{x}_{T}$ the common vertex of all simplices in
$\mathfrak{T}_{T}$. With this in mind, we additionally introduce the following
space generated by the Koszul operator ([2, Section 7.2]):

$\displaystyle\boldsymbol{\mathcal{G}}^{{\rm
c},k}(T)\coloneq(\boldsymbol{x}-\boldsymbol{x}_{T})\times\mathcal{P}^{k-1}(T)^{3}\qquad\text{for
}k\geq 1,$

and define $\boldsymbol{\mathcal{G}}^{{\rm
c},-1}(T)\coloneq\boldsymbol{\mathcal{G}}^{{\rm c},0}(T)\coloneq\\{0\\}$.
Observe that we have the following decomposition for $\mathcal{P}^{k}(T)^{3}$
(see [2, Corollary 7.4]):

$\mathcal{P}^{k}(T)^{3}=\nabla\mathcal{P}^{k+1}(T)\oplus\boldsymbol{\mathcal{G}}^{{\rm
c},k}(T),$ (13)

where the direct sum above is not orthogonal in general. Additionally, we
define the $L^{2}$-orthogonal projector on the space
$\boldsymbol{\mathcal{G}}^{{\rm c},k}(T)$ as ${\boldsymbol{\pi}}^{{\rm
c},k}_{\boldsymbol{\mathcal{G}},T}$. Then, the _divergence-preserving velocity
reconstruction_
${\boldsymbol{R}}_{T}^{k}:\uline{\boldsymbol{U}}_{T}^{k}\rightarrow\boldsymbol{\mathcal{RT}}^{k}(\mathfrak{T}_{T})$
is defined, for all
$\uline{\boldsymbol{v}}_{T}\in\uline{\boldsymbol{U}}_{T}^{k}$, as the first
component of the solution of the following mixed problem: Find
$({\boldsymbol{R}}_{T}^{k}\uline{\boldsymbol{v}}_{T},{\mathfrak{\psi}},\boldsymbol{\mathfrak{\zeta}})\in\boldsymbol{\mathcal{RT}}^{k}(\mathfrak{T}_{T})\times\mathcal{P}^{k}(\mathfrak{T}_{T})\times\boldsymbol{\mathcal{G}}^{{\rm
c},k-1}(T)$ such that

$\displaystyle({\boldsymbol{R}}_{T}^{k}\uline{\boldsymbol{v}}_{T})_{|\sigma}\cdot\boldsymbol{n}_{\sigma}$
$\displaystyle=(\boldsymbol{v}_{F}\cdot\boldsymbol{n}_{TF})_{|\sigma}$
$\displaystyle\forall\sigma\in\mathfrak{F}_{F},\,\forall F\in\mathcal{F}_{T},$
(14a)
$\displaystyle\int_{T}(\nabla\cdot{\boldsymbol{R}}_{T}^{k}\uline{\boldsymbol{v}}_{T}){\mathfrak{\phi}}$
$\displaystyle=\int_{T}(D_{T}^{k}\uline{\boldsymbol{v}}_{T}){\mathfrak{\phi}}$
$\displaystyle\forall{\mathfrak{\phi}}\in\mathcal{P}^{k}(\mathfrak{T}_{T}),$
(14b)
$\displaystyle\int_{T}{\boldsymbol{R}}_{T}^{k}\uline{\boldsymbol{v}}_{T}\cdot\boldsymbol{\mathfrak{\xi}}$
$\displaystyle=\int_{T}\boldsymbol{v}_{T}\cdot\boldsymbol{\mathfrak{\xi}}$
$\displaystyle\forall\boldsymbol{\mathfrak{\xi}}\in\boldsymbol{\mathcal{G}}^{{\rm
c},k-1}(T),$ (14c)
$\displaystyle\int_{T}{\boldsymbol{R}}_{T}^{k}\uline{\boldsymbol{v}}_{T}\cdot\boldsymbol{\mathfrak{w}}+\int_{T}(\nabla\cdot\boldsymbol{\mathfrak{w}}){\mathfrak{\psi}}+\int_{T}{\boldsymbol{\mathfrak{w}}}\cdot\boldsymbol{\mathfrak{\zeta}}$
$\displaystyle=\int_{T}\boldsymbol{v}_{T}\cdot\boldsymbol{\mathfrak{w}}$
$\displaystyle\forall\boldsymbol{\mathfrak{w}}\in\boldsymbol{\mathcal{RT}}^{k}_{0}(\mathfrak{T}_{T}).$
(14d)

######  Remark 1 (Allowing more than pyramidal meshes).

A similar divergence-preserving operator has been proposed in ††margin:  R#2
major point: checking and citing related work [40, Section 4.2] in the context
of finite elements pairs with continuous pressures. However, adapting it to
the current HHO framework will restrict the submesh $\mathfrak{T}_{T}$ to be
only a pyramidal submesh (or a vertex patch in the terminology of [40]).
Specifically, using the methodology introduced in [40, Proof of Theorem 12],
to prove the equation (17) below, it will be necessary to construct the
Lagrange hat function of $\boldsymbol{x}_{T}$ (a polynomial function
$q\in\mathcal{P}^{1}(\mathfrak{T}_{T})$ such that $q(\boldsymbol{x}_{T})=1$
and vanishes at the other vertices of $\mathfrak{T}_{T}$) and use its
properties with the crucial restriction that this hat function must vanish at
the boundary of $T$. This is only possible when $\mathfrak{T}_{T}$ is
pyramidal. In the current manuscript, we avoid this restriction using Lemma 3
below.

######  Lemma 2 (Properties of $\boldsymbol{R}_{T}^{k}$).

It holds:

  1. (i)

_Well-posedness._ For a given
$\uline{\boldsymbol{v}}_{T}\in\uline{\boldsymbol{U}}_{T}^{k}$, there exists a
unique solution to problem (14), and it holds that

$\displaystyle\|\boldsymbol{v}_{T}-\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}\|_{{L}^{2}({T})^{3}}\lesssim
h_{T}|\uline{\boldsymbol{v}}_{T}|_{1,\partial T}.$ (15)

  2. (ii)

_Approximation._ For all $\boldsymbol{v}\in\boldsymbol{H}^{k+1}(T)^{3}$, it
holds

$\displaystyle\|\boldsymbol{v}-\boldsymbol{R}_{T}^{k}(\uline{\boldsymbol{I}}_{T}^{k}\boldsymbol{v})\|_{{L}^{2}({T})^{3}}\lesssim
h_{T}^{k+1}|\boldsymbol{v}|_{H^{k+1}(T)^{3}}.$ (16) ††margin:  Referees #1 and
#2 remark: Now we have consistency for arbitrary $k\geq 1$

  3. (iii)

_Consistency._ For a given
$\uline{\boldsymbol{v}}_{T}\in\uline{\boldsymbol{U}}_{T}^{k}$, it holds, for
$k\geq 1$,

$\displaystyle{\boldsymbol{\pi}}^{k-1}_{T}(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T})={\boldsymbol{\pi}}^{k-1}_{T}(\boldsymbol{v}_{T}).$
(17)

The proof makes use of the following Lemma, whose proof is given in Appendix
A.

######  Lemma 3 (Raviart–Thomas lifting of the projection in
$\boldsymbol{\mathcal{G}}^{{\rm c},k-1}(T)$).

Let $T\in\mathcal{T}_{h}$ and a function $\boldsymbol{v}\in L^{2}(T)^{3}$ be
given. Then, for $k\geq 2$, there exists
$\widetilde{\boldsymbol{R}}_{T}^{k}(\boldsymbol{v})\in\boldsymbol{\mathcal{RT}}^{k}_{0}(\mathfrak{T}_{T})$
such that

$\displaystyle{\boldsymbol{\pi}}^{{\rm
c},k-1}_{\boldsymbol{\mathcal{G}},T}\widetilde{\boldsymbol{R}}_{T}^{k}(\boldsymbol{v})$
$\displaystyle={\boldsymbol{\pi}}^{{\rm
c},k-1}_{\boldsymbol{\mathcal{G}},T}\boldsymbol{v},$ (18a)
$\displaystyle\nabla\cdot\widetilde{\boldsymbol{R}}_{T}^{k}(\boldsymbol{v})$
$\displaystyle=0,$ (18b)
$\displaystyle\widetilde{\boldsymbol{R}}_{T}^{k}(\boldsymbol{v})\cdot\boldsymbol{n}_{\sigma}$
$\displaystyle=0\qquad\forall\sigma\in\mathfrak{F}_{T}^{\rm i},$ (18c)
$\displaystyle\|\widetilde{\boldsymbol{R}}_{T}^{k}(\boldsymbol{v})\|_{{L}^{2}({T})^{3}}$
$\displaystyle\lesssim\|\boldsymbol{v}\|_{{L}^{2}({T})^{3}}$ (18d)

######  Remark 4 (The common vertex assumption for $k\in\\{0,1\\}$).

The common vertex assumption described in Section 2.1 is not necessary for
$k\in\\{0,1\\}$ since, for those cases, $\boldsymbol{\mathcal{G}}^{{\rm
c},k-1}$ becomes the trivial space, and Lemma (3) is not needed.

###### Proof.

_(i) Well-posedness._ We prove this item in three parts starting with the
existence and uniqueness of
${\boldsymbol{R}}_{T}^{k}\uline{\boldsymbol{v}}_{T}$.
(i.A) _Existence, uniqueness, and decomposition of
${\boldsymbol{R}}_{T}^{k}\uline{\boldsymbol{v}}_{T}$._ The existence and
uniqueness of a solution to problem (14) follows from the classical theory of
mixed problems given the compatibility of the selected spaces; see, e.g., [14,
Section 14] and Lemma 3 below. In order to prove the a priori estimate (15),
we decompose ${\boldsymbol{R}}_{T}^{k}\uline{\boldsymbol{v}}_{T}$ as follows:

$\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}={\boldsymbol{\mathfrak{v}}}^{\prime}+\boldsymbol{\mathfrak{v}}_{0},$
(19)

where:

  * •

${\boldsymbol{\mathfrak{v}}}^{\prime}\in\boldsymbol{\mathcal{RT}}^{k}(\mathfrak{T}_{T})$
is a lifting of the boundary values defined by prescribing its DOFs as
follows:

$\displaystyle\forall\tau\in\mathfrak{T}_{T},$
$\displaystyle\quad\int_{\tau}\boldsymbol{\mathfrak{v}}^{\prime}\cdot\boldsymbol{\mathfrak{w}}$
$\displaystyle=0,$
$\displaystyle\forall\boldsymbol{\mathfrak{w}}\in\mathcal{P}^{k-1}(\tau)^{3},$
(20a) $\displaystyle\forall\sigma\in\mathfrak{F}_{T}^{\rm i},$
$\displaystyle\quad\int_{\sigma}(\boldsymbol{\mathfrak{v}}^{\prime}\cdot\boldsymbol{n}_{\sigma}){\mathfrak{\phi}}$
$\displaystyle=0$
$\displaystyle\forall{\mathfrak{\phi}}\in\mathcal{P}^{k}(\sigma),$ (20b)
$\displaystyle\forall F\in\mathcal{F}_{T},\,\forall\sigma\in\mathfrak{F}_{F},$
$\displaystyle\quad\boldsymbol{\mathfrak{v}}^{\prime}_{|\sigma}\cdot\boldsymbol{n}_{\sigma}$
$\displaystyle=(\boldsymbol{v}_{F}\cdot\boldsymbol{n}_{TF})_{|\sigma};$ (20c)

  * •

Letting
$\mathcal{P}^{k}_{0}(\mathfrak{T}_{T})\coloneq\left\\{{\mathfrak{\phi}}\in\mathcal{P}^{k}(\mathfrak{T}_{T}):\int_{T}{\mathfrak{\phi}}=0\right\\}$,
$\boldsymbol{\mathfrak{v}}_{0}$ is the first component of the unique solution
to the following mixed problem: Find
$(\boldsymbol{\mathfrak{v}}_{0},{\mathfrak{\psi}},\boldsymbol{\mathfrak{\zeta}})\in\boldsymbol{\mathcal{RT}}^{k}_{0}(\mathfrak{T}_{T})\times\mathcal{P}^{k}_{0}(\mathfrak{T}_{T})\times\boldsymbol{\mathcal{G}}^{{\rm
c},k-1}(T)$ such that

$\displaystyle\int_{T}(\nabla\cdot\boldsymbol{\mathfrak{v}}_{0}){\mathfrak{\phi}}$
$\displaystyle=\int_{T}(D_{T}^{k}\uline{\boldsymbol{v}}_{T}-\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime}){\mathfrak{\phi}}$
$\displaystyle\forall{\mathfrak{\phi}}\in\mathcal{P}^{k}_{0}(\mathfrak{T}_{T}),$
(21a)
$\displaystyle\int_{T}{\boldsymbol{\mathfrak{v}}}_{0}\cdot\boldsymbol{\mathfrak{\xi}}$
$\displaystyle=\int_{T}\boldsymbol{v}_{T}\cdot\boldsymbol{\mathfrak{\xi}}$
$\displaystyle\forall\boldsymbol{\mathfrak{\xi}}\in\boldsymbol{\mathcal{G}}^{{\rm
c},k-1}(T),$ (21b)
$\displaystyle\int_{T}\boldsymbol{\mathfrak{v}}_{0}\cdot\boldsymbol{\mathfrak{w}}+\int_{T}(\nabla\cdot\boldsymbol{\mathfrak{w}}){\mathfrak{\psi}}+\int_{T}{\boldsymbol{\mathfrak{w}}}\cdot\boldsymbol{\mathfrak{\zeta}}$
$\displaystyle=\int_{T}(\boldsymbol{v}_{T}-\boldsymbol{\mathfrak{v}}^{\prime})\cdot\boldsymbol{\mathfrak{w}}$
$\displaystyle\forall\boldsymbol{\mathfrak{w}}\in\boldsymbol{\mathcal{RT}}^{k}_{0}(\mathfrak{T}_{T}),$
(21c)

where we have used equation (20a) along with the fact that
$\boldsymbol{\mathcal{G}}^{{\rm
c},k-1}(T)\subset\mathcal{P}^{k-1}(\mathfrak{T}_{T})^{3}$ in the right hand
side of (14c) to write $\boldsymbol{v}_{T}$ instead of
$\boldsymbol{v}_{T}-\boldsymbol{\mathfrak{v}}^{\prime}$.

(i.B) _Boundedness_. We begin by proving the following estimate:

$\displaystyle\|{\boldsymbol{R}}_{T}^{k}\uline{\boldsymbol{v}}_{T}\|_{{L}^{2}({T})^{3}}$
$\displaystyle\lesssim\|{\boldsymbol{v}}_{T}\|_{{L}^{2}({T})^{3}}+h_{T}\|\uline{\boldsymbol{v}}_{T}\|_{1,T}+\sum_{F\in\mathcal{F}_{T}}h_{F}^{\frac{1}{2}}\|\boldsymbol{v}_{F}\|_{{L}^{2}({F})^{3}}.$
(22)

Let
${L}^{2}_{0}({T})\coloneq\left\\{\xi\in{L}^{2}({T}):\int_{T}\xi=0\right\\}$
and denote by $\theta\in{H}^{1}({T})\cap{L}^{2}_{0}({T})$ the solution of the
equation

$\displaystyle\int_{T}\nabla\theta\cdot\nabla\xi=\int_{T}(D_{T}^{k}\uline{\boldsymbol{v}}_{T}-\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime})\xi\qquad\forall\xi\in{H}^{1}({T})\cap{L}^{2}_{0}({T}).$
(23)

We recall that (23) is the weak form of the following strong Neumann problem

$\displaystyle-\Delta\theta$
$\displaystyle=(D_{T}^{k}\uline{\boldsymbol{v}}_{T}-\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime})$
$\displaystyle\text{in }T,$ (24a) $\displaystyle\frac{\partial\theta}{\partial
n}$ $\displaystyle=0$ $\displaystyle\text{on }\partial T.$ (24b)

Since
$\int_{T}(D_{T}^{k}\uline{\boldsymbol{v}}_{T}-\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime})=0$
by (10) with $q=1$ and using integration by parts for the integral having
$\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime}$ along with (20c), the
compatibility condition for problem (24) is satisfied, yielding existence and
uniqueness of $\theta$. Moreover, since $T$ is a convex polyhedron and the
forcing term is in ${L}^{2}({T})$, then $\theta\in
H^{2}(T)\cap{L}^{2}_{0}({T})$ (see [31, Section 8.2]) with

$\displaystyle|\theta|_{H^{2}(T)}\lesssim\|D_{T}^{k}\uline{\boldsymbol{v}}_{T}-\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})},$
(25)

where it can be checked following the argument in the reference that the
hidden constant does not depend on $T$. Setting $\xi=\theta$ in (23) and using
the Cauchy–Schwarz inequality followed by the Poincaré–Wirtinger inequality
$\|\zeta\|_{{L}^{2}({T})}\leq\frac{h_{T}}{\pi}|\zeta|_{{H}^{1}({T})}$ valid
for all $\zeta\in H^{1}(T)\cap L^{2}_{0}(\Omega)$ (see [46, 3]), we estimate

$\displaystyle|\theta|_{{H}^{1}({T})}\lesssim{h_{T}}\|D_{T}^{k}\uline{\boldsymbol{v}}_{T}-\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})}.$
(26)

Let now $\boldsymbol{z}_{0}\coloneq-\nabla\theta\in{H}^{1}({T})^{3}$. Defining
$\hat{\boldsymbol{z}}_{0}$
$\in\boldsymbol{\mathcal{RT}}^{k}_{0}(\mathfrak{T}_{T})$ as the interpolate of
$\boldsymbol{z}_{0}$ onto $\boldsymbol{\mathcal{RT}}^{k}(\mathfrak{T}_{T})$,
and using the commutation property
$\nabla\cdot\hat{\boldsymbol{z}}_{0}=\pi_{\mathfrak{T}_{T}}^{k}(\nabla\cdot\boldsymbol{z}_{0})$
with $\pi_{\mathfrak{T}_{T}}^{k}$ denoting the $L^{2}$-orthogonal projector
onto $\mathcal{P}^{k}(\mathfrak{T}_{T})$ (see, e.g., [6, Section 2.5.2]), it
is inferred that
$\nabla\cdot\hat{\boldsymbol{z}}_{0}=\pi_{\mathfrak{T}_{T}}^{k}(\nabla\cdot\boldsymbol{z}_{0})=(D_{T}^{k}\uline{\boldsymbol{v}}_{T}-\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime})$,
where we have used (24a) in the last step. Therefore, by (21a),
$\nabla\cdot(\boldsymbol{\mathfrak{v}}_{0}-{\hat{\boldsymbol{z}}_{0}})=0$.
Now, using Lemma 3, let
$\hat{\boldsymbol{z}}_{1}\coloneq\widetilde{\boldsymbol{R}}_{T}^{k}(\boldsymbol{v}_{T}-\hat{\boldsymbol{z}}_{0})\in\boldsymbol{\mathcal{RT}}^{k}_{0}(\mathfrak{T}_{T})$,
and set
$\hat{\boldsymbol{z}}\coloneq\hat{\boldsymbol{z}}_{0}+\hat{\boldsymbol{z}}_{1}$.
Using (18b), we get
$\nabla\cdot(\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}})=0$; moreover,
using (18a) and (21b), we have $\boldsymbol{\pi}^{{\rm
c},k-1}_{\boldsymbol{\cal
G},T}(\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}})=\boldsymbol{0}$.
Taking then
$\boldsymbol{\mathfrak{w}}=\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}}$
as a test function in (21c), we obtain

$\int_{T}\boldsymbol{\mathfrak{v}}_{0}\cdot(\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}})=\int_{T}(\boldsymbol{v}_{T}-\boldsymbol{\mathfrak{v}}^{\prime})\cdot(\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}}).$
(27)

Thus, it is readily seen that

$\displaystyle\|\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}}\|_{{L}^{2}({T})^{3}}^{2}$
$\displaystyle=\int_{T}(\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}})\cdot(\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}})$
$\displaystyle=\int_{T}(\boldsymbol{v}_{T}-\boldsymbol{\mathfrak{v}}^{\prime}-\hat{\boldsymbol{z}})\cdot(\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}})$
$\displaystyle\leq\left(\|\boldsymbol{v}_{T}-\hat{\boldsymbol{z}}\|_{{L}^{2}({T})^{3}}+\|\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})^{3}}\right)\|\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}}\|_{{L}^{2}({T})^{3}}$
$\displaystyle\leq\left(\|\boldsymbol{v}_{T}-\hat{\boldsymbol{z}}_{0}\|_{{L}^{2}({T})^{3}}+\|\hat{\boldsymbol{z}}_{1}\|_{{L}^{2}({T})^{3}}+\|\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})^{3}}\right)\|\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}}\|_{{L}^{2}({T})^{3}}$
$\displaystyle\lesssim\left(\|\boldsymbol{v}_{T}-\hat{\boldsymbol{z}}_{0}\|_{{L}^{2}({T})^{3}}+\|\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})^{3}}\right)\|\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}}\|_{{L}^{2}({T})^{3}}$
$\displaystyle\leq\left(\|\boldsymbol{v}_{T}-\boldsymbol{z}_{0}\|_{{L}^{2}({T})^{3}}+\|\boldsymbol{z}_{0}-\hat{\boldsymbol{z}}_{0}\|_{{L}^{2}({T})^{3}}+\|\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})^{3}}\right)\|\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}}\|_{{L}^{2}({T})^{3}}$
$\displaystyle\eqcolon(\mathfrak{T}_{1}+\mathfrak{T}_{2}+\mathfrak{T}_{3})\|\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}}\|_{{L}^{2}({T})^{3}},$
(28)

where we have used equation (27) in the second step, Cauchy–Schwarz and
triangle inequalities in the third step, the definition of
$\hat{\boldsymbol{z}}$ for the first term in the parentheses and a triangle
inequality in the fourth step, the definition of $\hat{\boldsymbol{z}}_{1}$
along with the bound (18d) in the fifth step, and again a triangle inequality
after inserting $\pm{\boldsymbol{z}}_{0}$ into the first norm in the sixth
step.
To estimate $\mathfrak{T}_{1}$, we begin using the triangle inequality to
obtain

$\displaystyle\mathfrak{T}_{1}$
$\displaystyle\leq\|\boldsymbol{v}_{T}\|_{{L}^{2}({T})^{3}}+\|\boldsymbol{z}_{0}\|_{{L}^{2}({T})^{3}}\lesssim\|\boldsymbol{v}_{T}\|_{{L}^{2}({T})^{3}}+h_{T}\|D_{T}^{k}\uline{\boldsymbol{v}}_{T}-\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})}$
$\displaystyle\lesssim\|\boldsymbol{v}_{T}\|_{{L}^{2}({T})^{3}}+h_{T}\left(\|D_{T}^{k}\uline{\boldsymbol{v}}_{T}\|_{{L}^{2}({T})}+\|\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})}\right),$
(29)

where we have used (26) to bound the second term in the second step, and a
triangle inequality in the last step. For the first term in parentheses,
integrating by parts the right-hand side of (10), applying Cauchy–Schwarz and
discrete trace inequalities, and taking the supremum over
$q\in\mathcal{P}^{k}(T)$, we obtain

$\displaystyle\|D_{T}^{k}\uline{\boldsymbol{v}}_{T}\|_{{L}^{2}({T})}$
$\displaystyle\lesssim\|\uline{\boldsymbol{v}}_{T}\|_{1,T}.$ (30)

To bound the second term in parentheses, we use a discrete inverse inequality
to write

$\displaystyle\|\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})}^{2}\lesssim\sum_{\tau\in\mathfrak{T}_{T}}|\boldsymbol{\mathfrak{v}}^{\prime}|_{H^{1}(\tau)^{3}}^{2}\lesssim
h_{\tau}^{-2}\sum_{\tau\in\mathfrak{T}_{T}}\|\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({\tau})^{3}}^{2}\lesssim
h_{T}^{-2}\|\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})^{3}}^{2},$ (31)

where the fact that $h_{\tau}^{-1}h_{T}\lesssim 1$ for regular mesh sequences
(see [16, Eq. (1.4)]) has been used in the last step. Now, to estimate
$\|\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})^{3}}$, we combine (12)
and (20) to obtain

$\displaystyle\|\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})^{3}}^{2}\simeq\sum_{\sigma\in\mathfrak{F}_{F},F\in\mathcal{F}_{T}}h_{\sigma}\|\boldsymbol{v}_{F}\cdot\boldsymbol{n}_{TF}\|_{{L}^{2}({\sigma})}^{2}\leq\sum_{F\in\mathcal{F}_{T}}h_{F}\|\boldsymbol{v}_{F}\|_{{L}^{2}({F})^{3}}^{2},$
(32)

where we have used the fact that $\sigma\subset F$, the inequality
$h_{\sigma}\leq h_{F}$ valid for any $\sigma\in\mathfrak{F}_{F}$ and any
$F\in\mathcal{F}_{T}$, and the Hölder inequality with exponents $(2,\infty)$
along with $\|\boldsymbol{n}_{TF}\|_{L^{\infty}(F)^{3}}=1$ for the third step.
Therefore plugging (32) into (31), we obtain

$\displaystyle\|\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})}^{2}\lesssim\sum_{F\in\mathcal{F}_{T}}h_{F}^{-1}\|\boldsymbol{v}_{F}\|_{{L}^{2}({F})^{3}}^{2},$

where we have used the equivalence $h_{F}\simeq h_{T}$ valid for regular mesh
sequences (see [16, Eq. (1.6)]). Now, using the following bound valid for non-
negative real numbers $a_{i}$

$\sum_{i}a_{i}^{2}\leq\left(\sum_{i}a_{i}\right)^{2},$ (33)

for the previous inequality, and then plugging the result along with (30) into
(29), it is inferred that

$\mathfrak{T}_{1}\lesssim\|\boldsymbol{v}_{T}\|_{{L}^{2}({T})^{3}}+h_{T}\|\uline{\boldsymbol{v}}_{T}\|_{1,T}+\sum_{F\in\mathcal{F}_{T}}h_{F}^{\frac{1}{2}}\|\boldsymbol{v}_{F}\|_{{L}^{2}({F})^{3}}.$
(34)

To bound the term $\mathfrak{T}_{2}$ in (28), we use standard interpolation
estimates for $\hat{\boldsymbol{z}}_{0}$ (see, e.g., [6, Proposition 2.5.4])
followed by (25) to write

$\mathfrak{T}_{2}\lesssim h_{T}|\theta|_{H^{2}(T)}\lesssim
h_{T}\|D_{T}^{k}\uline{\boldsymbol{v}}_{T}-\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})^{3}}\lesssim
h_{T}\|\uline{\boldsymbol{v}}_{T}\|_{1,T}+\sum_{F\in\mathcal{F}_{T}}h_{F}^{\frac{1}{2}}\|\boldsymbol{v}_{F}\|_{{L}^{2}({F})^{3}},$
(35)

where we have used a triangle inequality followed by (30), (31), and (32) to
conclude. Plugging (34) and (35) into (28), using (32) to estimate
$\mathfrak{T}_{3}$, and simplifying, we obtain

$\displaystyle\|\boldsymbol{\mathfrak{v}}_{0}-\hat{\boldsymbol{z}}\|_{{L}^{2}({T})^{3}}$
$\displaystyle\lesssim\|\boldsymbol{v}_{T}\|_{{L}^{2}({T})^{3}}+h_{T}\|\uline{\boldsymbol{v}}_{T}\|_{1,T}+\sum_{F\in\mathcal{F}_{T}}h_{F}^{\frac{1}{2}}\|\boldsymbol{v}_{F}\|_{{L}^{2}({F})^{3}}.$
(36)

Using the decomposition (19) followed by triangle inequalities, we finally get

$\displaystyle\|{\boldsymbol{v}}_{T}-{\boldsymbol{R}}_{T}^{k}\uline{\boldsymbol{v}}_{T}\|_{{L}^{2}({T})^{3}}$
$\displaystyle\leq\|{\boldsymbol{v}}_{T}-\hat{\boldsymbol{z}}\|_{{L}^{2}({T})^{3}}+\|\hat{\boldsymbol{z}}-\boldsymbol{\mathfrak{v}}_{0}\|_{{L}^{2}({T})^{3}}+\|\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})^{3}}$
$\displaystyle\lesssim\|{\boldsymbol{v}}_{T}-\hat{\boldsymbol{z}}_{0}\|_{{L}^{2}({T})^{3}}+\|\hat{\boldsymbol{z}}_{1}\|_{{L}^{2}({T})^{3}}+\|\hat{\boldsymbol{z}}-\boldsymbol{\mathfrak{v}}_{0}\|_{{L}^{2}({T})^{3}}+\|\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})^{3}}$
$\displaystyle\lesssim\|{\boldsymbol{v}}_{T}-\hat{\boldsymbol{z}}_{0}\|_{{L}^{2}({T})^{3}}+\|\hat{\boldsymbol{z}}-\boldsymbol{\mathfrak{v}}_{0}\|_{{L}^{2}({T})^{3}}+\|\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})^{3}}$
$\displaystyle\lesssim\|{\boldsymbol{v}}_{T}-\boldsymbol{z}_{0}\|_{{L}^{2}({T})^{3}}+\|\boldsymbol{z}_{0}-\hat{\boldsymbol{z}}_{0}\|_{{L}^{2}({T})^{3}}+\|\hat{\boldsymbol{z}}-\boldsymbol{\mathfrak{v}}_{0}\|_{{L}^{2}({T})^{3}}+\|\boldsymbol{\mathfrak{v}}^{\prime}\|_{{L}^{2}({T})^{3}}$
$\displaystyle\lesssim\|\boldsymbol{v}_{T}\|_{{L}^{2}({T})^{3}}+h_{T}\|\uline{\boldsymbol{v}}_{T}\|_{1,T}+\sum_{F\in\mathcal{F}_{T}}h_{F}^{\frac{1}{2}}\|\boldsymbol{v}_{F}\|_{{L}^{2}({F})^{3}},$

where we have used the definition of $\hat{\boldsymbol{z}}_{1}$ along with the
bound (18d) in the third step, a triangle inequality in the fourth step, and
the bounds (32), and (34)–(36) in the last step. Inserting
$\pm{\boldsymbol{v}}_{T}$ into the left-hand side of (22) and using a triangle
inequality followed by the above estimate, (22) follows.
(i.C) _Proof of the bound ( 15)_. Recalling that
$\uline{\boldsymbol{I}}_{T}^{k}$ is obtained restricting the global
interpolator (6) to an element $T$, letting
$\hat{\boldsymbol{\mathfrak{v}}}\coloneq{\boldsymbol{R}}_{T}^{k}(\uline{\boldsymbol{I}}_{T}^{k}\boldsymbol{v}_{T})$,
and using the triangle inequality, we get that

$\|\boldsymbol{v}_{T}-{\boldsymbol{R}}_{T}^{k}\uline{\boldsymbol{v}}_{T}\|_{{L}^{2}({T})^{3}}\leq\|\boldsymbol{v}_{T}-\hat{\boldsymbol{\mathfrak{v}}}\|_{{L}^{2}({T})^{3}}+\|\hat{\boldsymbol{\mathfrak{v}}}-{\boldsymbol{R}}_{T}^{k}\uline{\boldsymbol{v}}_{T}\|_{{L}^{2}({T})^{3}}\eqcolon\mathfrak{T}_{1}+\mathfrak{T}_{2}.$
(37)

By condition (14b), we have that
$\nabla\cdot\hat{\boldsymbol{\mathfrak{v}}}=D_{T}^{k}(\uline{\boldsymbol{I}}_{T}^{k}\boldsymbol{v}_{T})\in\mathcal{P}^{k}(T)\subset\mathcal{P}^{k}(\mathfrak{T}_{T})$.
But, since $\boldsymbol{v}_{T}\in\mathcal{P}^{k}(T)$, the commutation property
(11) gives
$D_{T}^{k}(\uline{\boldsymbol{I}}_{T}^{k}\boldsymbol{v}_{T})=\pi_{T}^{k}(\nabla\cdot\boldsymbol{v}_{T})=\nabla\cdot\boldsymbol{v}_{T}$,
so that $\nabla\cdot(\hat{\boldsymbol{\mathfrak{v}}}-\boldsymbol{v}_{T})=0$.
In addition, by (14c) we have ${\boldsymbol{\pi}}^{{\rm
c},k-1}_{\boldsymbol{\mathcal{G}},T}(\hat{\boldsymbol{\mathfrak{v}}}-\boldsymbol{v}_{T})=0$,
and then observing that
$\hat{\boldsymbol{\mathfrak{v}}}-\boldsymbol{v}_{T}\in\boldsymbol{\mathcal{RT}}^{k}_{0}(\mathfrak{T}_{T})$,
and taking
$\boldsymbol{\mathfrak{w}}=\hat{\boldsymbol{\mathfrak{v}}}-\boldsymbol{v}_{T}$
in (14d), it is inferred that
$\|\hat{\boldsymbol{\mathfrak{v}}}-\boldsymbol{v}_{T}\|_{{L}^{2}({T})^{3}}^{2}=0$,
hence $\mathfrak{T}_{1}=0$.

Let us now estimate the term $\mathfrak{T}_{2}$. By linearity of
${\boldsymbol{R}}_{T}^{k}$, we can write
$\mathfrak{T}_{2}=\|{\boldsymbol{R}}_{T}^{k}(\uline{\boldsymbol{I}}_{T}^{k}\boldsymbol{v}_{T}-\uline{\boldsymbol{v}}_{T})\|_{{L}^{2}({T})^{3}}$.
Hence, using the bound (22), the fact that
$(\uline{\boldsymbol{I}}_{T}^{k}\boldsymbol{v}_{T}-\uline{\boldsymbol{v}}_{T})_{T}=0$
and
$(\uline{\boldsymbol{I}}_{T}^{k}\boldsymbol{v}_{T}-\uline{\boldsymbol{v}}_{T})_{F}=(\boldsymbol{v}_{T}-\boldsymbol{v}_{F})$
for all $F\in\mathcal{F}_{T}$, and recalling the definition (7), we can write

$\displaystyle\mathfrak{T}_{2}$ $\displaystyle\lesssim
h_{T}\|\uline{\boldsymbol{I}}_{T}^{k}\boldsymbol{v}_{T}-\uline{\boldsymbol{v}}_{T}\|_{1,T}+\sum_{F\in\mathcal{F}_{T}}h_{F}^{\frac{1}{2}}\|\boldsymbol{v}_{T}-\boldsymbol{v}_{F}\|_{{L}^{2}({F})^{3}}$
$\displaystyle=h_{T}|\uline{\boldsymbol{v}}_{T}|_{1,\partial
T}+\sum_{F\in\mathcal{F}_{T}}h_{F}h_{F}^{-\frac{1}{2}}\|\boldsymbol{v}_{T}-\boldsymbol{v}_{F}\|_{{L}^{2}({F})^{3}}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\lesssim}\marginpar{\normalcolor\raggedright\scriptsize\textsf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}
R\\#2 point: Typo}}h_{T}|\uline{\boldsymbol{v}}_{T}|_{1,\partial
T},\@add@raggedright$

where we have used the inequality $h_{F}\leq h_{T}$ in the last step. Plugging
this last bound along with $\mathfrak{T}_{1}=0$ into (37), the conclusion
follows.
_(ii) Approximation._ To prove the approximation estimate (16), let
$T\in\mathcal{T}_{h}$ and denote, for the sake of brevity, by
$\hat{\uline{\boldsymbol{v}}}_{T}\coloneq\uline{\boldsymbol{I}}_{T}^{k}\boldsymbol{v}$
the interpolate of $\boldsymbol{v}$ on $\uline{\boldsymbol{U}}_{T}^{k}$. We
begin using the triangle inequality to write

$\displaystyle\|\boldsymbol{v}-{\boldsymbol{R}}_{T}^{k}(\hat{\uline{\boldsymbol{v}}}_{T})\|_{{L}^{2}({T})^{3}}$
$\displaystyle\leq\|\boldsymbol{v}-\hat{\boldsymbol{v}}_{T}\|_{{L}^{2}({T})^{3}}+\|\hat{\boldsymbol{v}}_{T}-{\boldsymbol{R}}_{T}^{k}(\hat{\uline{\boldsymbol{v}}}_{T})\|_{{L}^{2}({T})^{3}}$
$\displaystyle\lesssim\|\boldsymbol{v}-{\boldsymbol{\pi}}_{T}^{k}{\boldsymbol{v}}\|_{{L}^{2}({T})^{3}}+h_{T}|\hat{\uline{\boldsymbol{v}}}_{T}|_{1,\partial
T}\eqcolon\mathfrak{T}_{1}+\mathfrak{T}_{2},$ (38)

where in the last step we have used the definition of
$\uline{\boldsymbol{I}}_{T}^{k}$ and the bound (15) for the first and second
terms, respectively. To bound $\mathfrak{T}_{1}$ we use (5a) with
$(l,m,r,s)=(k,0,2,k+1)$, so we get

$\mathfrak{T}_{1}\lesssim h_{T}^{k+1}|\boldsymbol{v}|_{H^{k+1}(T)^{3}}.$ (39)

Now to bound $\mathfrak{T}_{2}$ we first take the square, use the definition
(7) of the boundary seminorm and the equivalence $h_{T}\simeq h_{F}$ (valid
for regular meshes) to obtain

$\displaystyle(\mathfrak{T}_{2})^{2}$
$\displaystyle=h_{T}^{2}\sum_{F\in{\mathcal{F}}_{T}}\\!h_{F}^{-1}\|\hat{{\boldsymbol{v}}}_{F}-\hat{{\boldsymbol{v}}}_{T}\|^{2}_{L^{2}(F)^{3}}\lesssim\sum_{F\in{\mathcal{F}}_{T}}\\!h_{T}\left(\|{\boldsymbol{\pi}}_{F}^{k}\boldsymbol{v}-\boldsymbol{v}\|^{2}_{L^{2}(F)^{3}}+\|\boldsymbol{v}-{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}\|^{2}_{L^{2}(F)^{3}}\right),$

where in the last step we have used the Young inequality. Now using a triangle
inequality and standard properties of the $L^{2}$-projectors
${\boldsymbol{\pi}}_{F}^{k}$ and ${\boldsymbol{\pi}}_{T}^{k}$ on $F$, we have

$\|{\boldsymbol{\pi}}_{F}^{k}\boldsymbol{v}-\boldsymbol{v}\|^{2}_{L^{2}(F)^{3}}=\inf_{\boldsymbol{w}\in\mathcal{P}^{k}(F)^{3}}\|\boldsymbol{w}-\boldsymbol{v}\|^{2}_{L^{2}(F)^{3}}\leq\|{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}-\boldsymbol{v}\|^{2}_{L^{2}(F)^{3}}.$
(40)

Thus using this for $F\in{\mathcal{F}}_{T}$, the bound (33), and then taking
the square root, it is inferred that

$\mathfrak{T}_{2}\lesssim\sum_{F\in{\mathcal{F}}_{T}}\\!h_{T}^{\frac{1}{2}}\|\boldsymbol{v}-{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}\|_{L^{2}(F)^{3}}.$

Finally, using (5b) with $(l,m,r,s)=(k,0,2,k+1)$ to bound $\mathfrak{T}_{2}$
along with (39), and plugging the result into (38), we conclude.
_(iii) Consistency._ To simplify the notation let us define the space
$\boldsymbol{\mathcal{G}}^{k-1}(T)\coloneq\nabla\mathcal{P}^{k}(T)^{3}$, and
let ${\boldsymbol{\pi}}^{k-1}_{\boldsymbol{\mathcal{G}},T}$ its
$L^{2}$-orthogonal projector. In addition, let
$\boldsymbol{q}\coloneq{\boldsymbol{\pi}}_{T}^{k-1}(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}-\boldsymbol{v}_{T})$.
As mentioned before, the decomposition (13) is not necessarily orthogonal;
nevertheless, by [15, Lemma 1] there exists a _recovery operator_
$\boldsymbol{\mathfrak{R}}_{\boldsymbol{\mathcal{G}},\boldsymbol{\mathcal{G}}^{{\rm
c}}}:\boldsymbol{\mathcal{G}}^{k-1}(T)\times\boldsymbol{\mathcal{G}}^{{\rm
c},k-1}(T)\rightarrow\mathcal{P}^{k-1}(T)^{3}$ such that
$\boldsymbol{q}=\boldsymbol{\mathfrak{R}}_{\boldsymbol{\mathcal{G}},\boldsymbol{\mathcal{G}}^{{\rm
c}}}({\boldsymbol{\pi}}^{k-1}_{\boldsymbol{\mathcal{G}},T}\boldsymbol{q},{{\boldsymbol{\pi}}^{{\rm
c},k-1}_{\boldsymbol{\mathcal{G}}^{,}T}\boldsymbol{q}})$, and
$\|\boldsymbol{q}\|_{L^{2}(T)^{3}}\simeq\|{\boldsymbol{\pi}}^{k-1}_{\boldsymbol{\mathcal{G}},T}\boldsymbol{q}\|_{L^{2}(T)^{3}}+\|{{\boldsymbol{\pi}}^{{\rm
c},k-1}_{\boldsymbol{\mathcal{G}}^{,}T}\boldsymbol{q}}\|_{L^{2}(T)^{3}}$.
Using this last equation and the linearity of ${\boldsymbol{\pi}}_{T}^{k-1}$,
it is enough to show that
$\|{\boldsymbol{\pi}}^{k-1}_{\boldsymbol{\mathcal{G}},T}\boldsymbol{q}\|_{L^{2}(T)^{3}}=\|{{\boldsymbol{\pi}}^{{\rm
c},k-1}_{\boldsymbol{\mathcal{G}}^{,}T}\boldsymbol{q}}\|_{L^{2}(T)^{3}}=0$.
Since ${\boldsymbol{\mathcal{G}}^{{\rm
c},k-1}}\subset\mathcal{P}^{k-1}(T)^{3}$, we have ${{\boldsymbol{\pi}}^{{\rm
c},k-1}_{\boldsymbol{\mathcal{G}}^{,}T}}={{\boldsymbol{\pi}}^{{\rm
c},k-1}_{\boldsymbol{\mathcal{G}}^{,}T}}\circ{\boldsymbol{\pi}}^{k-1}_{T}$,
and by the equation (14c), we infer that $\|{\boldsymbol{\pi}}^{{\rm
c},k-1}_{\boldsymbol{\mathcal{G}}^{,}T}\boldsymbol{q}\|_{L^{2}(T)^{3}}=0$.
Using (14b) along with an integration by parts and the boundary condition
(14a), we infer that, for all $\phi\in\mathcal{P}^{k}(T)$,

$\displaystyle-\int_{T}\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}\cdot\nabla\phi+\cancel{\sum_{F\in\mathcal{T}_{T}}\int_{F}({\boldsymbol{v}}_{F}\cdot\boldsymbol{n}_{TF})\phi}=\int_{T}(D_{T}^{k}\uline{\boldsymbol{v}}_{T})\phi=-\int_{T}\boldsymbol{v}_{T}\cdot\nabla\phi+\cancel{\sum_{F\in{\mathcal{F}}_{T}}\int_{F}(\boldsymbol{v}_{F}\cdot\boldsymbol{n}_{TF})\phi},$

where in the last step we have used the definition (10) of $D_{T}^{k}$. This
shows that
${\boldsymbol{\pi}}^{k-1}_{\boldsymbol{\mathcal{G}},T}(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}-{\boldsymbol{v}}_{T})=0$.
Finally, using
${\boldsymbol{\pi}}^{k-1}_{\boldsymbol{\mathcal{G}},T}={\boldsymbol{\pi}}^{k-1}_{\boldsymbol{\mathcal{G}},T}\circ{\boldsymbol{\pi}}^{k-1}_{T}$,
we obtain
$\|{\boldsymbol{\pi}}^{k-1}_{\boldsymbol{\mathcal{G}},T}\boldsymbol{q}\|_{L^{2}(T)^{3}}=0$,
and (17) follows. ∎

######  Remark 5 (Convexity assumption).

††margin:  Referees #1 & #2 remarks: Where the convexity enters and non-convex
elements

The convexity assumption introduced in Section 2.1 can be relaxed if instead
we make the following assumptions for $T\in\mathcal{T}_{h}$:

  1. 1.

There exist a point $\widetilde{\boldsymbol{x}}_{T}\in T$ such that $T$ is a
star-shaped with respect to it.

  2. 2.

There exists $\tau\in\mathfrak{T}_{T}$ such that contains the ball
$\boldsymbol{\mathcal{B}}(\widetilde{\boldsymbol{x}}_{T},r_{\tau})$ where
$r_{\tau}$ denotes the inradius of $\tau$.

Then instead of solving the problems (23)–(24), invoke the Lemma III.3.1 of
[28] to obtain $z_{0}\in H_{0}^{1}(T)^{3}$ such that
$\nabla\cdot\boldsymbol{z}_{0}=D_{T}^{k}\uline{\boldsymbol{v}}_{T}-\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime}$
with
$\|\boldsymbol{z}_{0}\|_{H^{1}(T)^{3}}\lesssim\|D_{T}^{k}\uline{\boldsymbol{v}}_{T}-\nabla\cdot\boldsymbol{\mathfrak{v}}^{\prime}\|_{L^{2}(T)}$,
and use [28, Eq. (II.5.5)] to get a Poincaré-like inequality
$\|\boldsymbol{z}_{0}\|_{{L}^{2}({T})^{3}}\lesssim
h_{T}|\boldsymbol{z}_{0}|_{H^{1}(T)^{3}}$ and use it in (29) in the proof of
item (i) of Lemma 2.

Let $\boldsymbol{\mathcal{RT}}^{k}(\mathfrak{T}_{h})$ denote the global
(${\boldsymbol{H}}_{\text{div}}({\Omega})$-conforming) Raviart–Thomas–Nédélec
space on $\mathfrak{T}_{h}$. We define the _global velocity reconstruction_
$\boldsymbol{R}_{h}^{k}:\uline{\boldsymbol{U}}_{h}^{k}\rightarrow\boldsymbol{\mathcal{RT}}^{k}(\mathfrak{T}_{h})$
patching the local contributions: For all
$\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}_{h}^{k}$,

$\displaystyle(\boldsymbol{R}_{h}^{k}\uline{\boldsymbol{v}}_{h})_{|T}\coloneq\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}\qquad\forall
T\in{\cal{T}}_{h}.$

Note that $\boldsymbol{R}_{h}^{k}\uline{\boldsymbol{v}}_{h}$ is well-defined,
since its normal components across each mesh interface are continuous as a
consequence of (14a) combined with the single-valuedness of interface
unknowns.

######  Proposition 6 (Sobolev inequalities for the velocity reconstruction).

It holds, for all $r\in[1,6]$ and all
$\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}_{h,0}^{k}$,

$\|\boldsymbol{R}_{h}^{k}\uline{\boldsymbol{v}}_{h}\|_{L^{r}(\Omega)^{3}}\lesssim\|\uline{\boldsymbol{v}}_{h}\|_{1,h},$
(41)

where the hidden constant is independent of both $h$ and
$\uline{\boldsymbol{v}}_{h}$, but possibly depends on $\Omega$, $k$, $r$, and
the mesh regularity parameter.

###### Proof.

Let a mesh element $T\in{\mathcal{T}}_{h}$ be fixed. Inserting
$\pm\boldsymbol{v}_{T}$ into the norm and using a triangle inequality, we can
write

$\begin{split}\|\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}\|_{L^{r}(T)^{3}}&\leq\|\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}-{\boldsymbol{v}}_{T}\|_{L^{r}(T)^{3}}+\|{\boldsymbol{v}}_{T}\|_{L^{r}(T)^{3}}\\\
&=\left(\sum_{\tau\in\mathfrak{T}_{T}}\|\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}-{\boldsymbol{v}}_{T}\|_{L^{r}(\tau)^{3}}^{r}\right)^{\frac{1}{r}}+\|{\boldsymbol{v}}_{T}\|_{L^{r}(T)^{3}}.\end{split}$
(42)

From the discrete Lebesgue embeddings proved in [16, Lemma 1.25], it follows
that, for all $(\alpha,\beta)\in[1,+\infty]$, all
$X\in\mathcal{T}_{h}\cup\mathfrak{T}_{h}$, and all
$\zeta\in\mathcal{P}^{l}(X)$ for $l\geq 0$,

$\|\zeta\|_{L^{\alpha}(X)}\lesssim
h_{X}^{\frac{3}{\alpha}-\frac{3}{\beta}}\|\zeta\|_{L^{\beta}(X)},$ (43)

with hidden constant independent of $h$, $X$, and $\zeta$, but possibly
depending on $l$, $\alpha$, $\beta$, and the mesh regularity parameter. Since
$(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}-{\boldsymbol{v}}_{T})_{|\tau}\in\mathcal{P}^{k+1}(\tau)^{3}$,
we use (43) for $(X,l,\alpha,\beta)=(\tau,k+1,r,2)$ in the term in parentheses
of (42) to write

$\displaystyle\sum_{\tau\in\mathfrak{T}_{T}}\|\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}-{\boldsymbol{v}}_{T}\|_{L^{2}(\tau)^{r}}^{r}$
$\displaystyle\lesssim\sum_{\tau\in\mathfrak{T}_{T}}h_{\tau}^{r\left(\frac{3}{r}-\frac{3}{2}\right)}\|\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}-{\boldsymbol{v}}_{T}\|_{L^{2}(\tau)^{3}}^{r}$
(44) $\displaystyle\lesssim
h_{T}^{r\left(\frac{3}{r}-\frac{3}{2}\right)}\|\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}-{\boldsymbol{v}}_{T}\|_{L^{2}(T)^{3}}^{r}\lesssim
h_{T}^{r\left(\frac{3}{r}-\frac{1}{2}\right)}|\uline{\boldsymbol{v}}_{T}|_{1,\partial
T}^{r},$

where we have used $\tau\subset T$ and $h_{\tau}\leq h_{T}$ for all
$\tau\in\mathfrak{T}_{T}$ along with the uniform bound (3) on
$\operatorname{card}(\mathfrak{T}_{T})$ in the second step, and the estimate
(15) to conclude. Plugging (44) into (42), raising the resulting inequality to
the $r$-th power, using the inequality $(a+b)^{r}\lesssim a^{r}+b^{r}$ valid
for any nonnegative real numbers $a$ and $b$, and summing over
$T\in{\mathcal{T}}_{h}$, we get

$\|\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}\|_{L^{r}(\Omega)^{3}}^{r}\lesssim\sum_{T\in{\mathcal{T}}_{h}}h_{T}^{\frac{6-r}{2}}|\uline{\boldsymbol{v}}_{T}|_{1,\partial
T}^{r}+\|{\boldsymbol{v}}_{h}\|_{L^{r}(\Omega)^{3}}^{r}.$

The proof now continues as that of [11, Proposition 3]. The details are
omitted for the sake of conciseness. ∎

###  2.5 Gradient reconstruction on a submesh

Let an element $T\in{\mathcal{T}}_{h}$ be fixed. For every face
$\sigma\in\mathfrak{F}_{T}^{\rm i}$, we introduce an arbitrary but fixed
ordering of the elements $\tau_{1}$ and $\tau_{2}$ such that
$\sigma\subset\partial\tau_{1}\cap\partial\tau_{2}$, and let
$\boldsymbol{n}_{\sigma}\coloneq\boldsymbol{n}_{\tau_{1}\sigma}=-\boldsymbol{n}_{\tau_{2}\sigma}$,
where $\boldsymbol{n}_{\tau_{i}\sigma},i\in\\{1,2\\}$, denotes the unit vector
normal to $\sigma$ pointing out of $\tau_{i}$ (see Figure 1). With this
convention, for every scalar-valued function $\zeta$ admitting a possibly two-
valued trace on $\sigma$, we define the jump of $\zeta$ across $\sigma$ as

$\llbracket\zeta\rrbracket_{\sigma}\coloneq\zeta_{|\tau_{1}}-\zeta_{|\tau_{2}}.$
(45)

When applied to vector- or tensor-valued functions, the jump operator acts
component-wise.

For any polyomial degree $l\geq 0$, we then define the local gradient
reconstruction
$\boldsymbol{G}_{\mathfrak{T}_{T}}^{l}:\uline{\boldsymbol{U}}_{T}^{k}\rightarrow\mathcal{P}^{l}(\mathfrak{T}_{T})^{3\times
3}$ such that, for all
$\uline{\boldsymbol{v}}_{T}\in\uline{\boldsymbol{U}}_{T}^{k}$ and all
${\boldsymbol{\tau}}\in\mathcal{P}^{l}(\mathfrak{T}_{T})^{3\times 3}$,

$\displaystyle\int_{T}\boldsymbol{G}^{l}_{\mathfrak{T}_{T}}\uline{\boldsymbol{v}}_{T}:{\boldsymbol{\tau}}$
$\displaystyle=\int_{T}\nabla\boldsymbol{v}_{T}:{\boldsymbol{\tau}}+\sum_{F\in{\mathcal{F}}_{T}}\int_{F}(\boldsymbol{v}_{F}-\boldsymbol{v}_{T})\cdot{\boldsymbol{\tau}}\boldsymbol{n}_{TF}$
(46a)
$\displaystyle=-\int_{T}\boldsymbol{v}_{T}\cdot(\nabla\cdot{\boldsymbol{\tau}})+\sum_{\sigma\in\mathfrak{F}_{T}^{\rm
i}}\int_{\sigma}\boldsymbol{v}_{T}\cdot\llbracket{\boldsymbol{\tau}}\rrbracket_{\sigma}\boldsymbol{n}_{\sigma}+\sum_{F\in{\mathcal{F}}_{T}}\int_{F}\boldsymbol{v}_{F}\cdot{\boldsymbol{\tau}}\boldsymbol{n}_{TF},$
(46b)

where we have used an integration by parts to pass to the second line. The
above definition is an extension of the operator
$\boldsymbol{G}_{T}^{l}:\uline{\boldsymbol{U}}_{T}^{k}\rightarrow\mathcal{P}^{l}(T)^{3\times
3}$ introduced in [11, 23], which is defined using
$\mathcal{P}^{l}(T)^{3\times 3}$ instead of
$\mathcal{P}^{l}(\mathfrak{T}_{T})^{3\times 3}$ as a test space, and thus we
have
$\boldsymbol{\pi}_{T}^{l}\boldsymbol{G}_{\mathfrak{T}_{T}}^{l}=\boldsymbol{G}_{T}^{l}$.
The gradient reconstruction $\boldsymbol{G}_{T}^{k}$ will be used in the
viscous term, while the enriched gradient reconstruction
$\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}$ will be used in the convective
term (see Section 3.3).

######  Lemma 7 (Properties of $\boldsymbol{G}_{\mathfrak{T}_{T}}^{l}$).

The operator $\boldsymbol{G}_{\mathfrak{T}_{T}}^{l}$ has the following
properties:

  1. 1.

_Boundedness._ For all
$\uline{\boldsymbol{v}}_{T}\in\uline{\boldsymbol{U}}_{T}^{k}$, it holds

$\|\boldsymbol{G}^{l}_{\mathfrak{T}_{T}}\uline{\boldsymbol{v}}_{T}\|_{L^{2}(T)^{3\times
3}}\lesssim\|\uline{\boldsymbol{v}}_{T}\|_{1,T}.$ (47)

  2. 2.

_Consistency._ For all $\boldsymbol{v}\in H^{k+1}(T)^{3}$ and all $l>k$, it
holds,

$\|\boldsymbol{G}^{l}_{\mathfrak{T}_{T}}\uline{\boldsymbol{I}}_{T}^{k}{\boldsymbol{v}}-\nabla\boldsymbol{v}\|_{L^{2}(T)^{3\times
3}}\lesssim h_{T}^{k}|\boldsymbol{v}|_{H^{k+1}(T)^{3}}.$ (48)

###### Proof.

_(i) Boundedness._ The proof is the same as that of [23, Proposition 1].
_(ii) Consistency._ Let $\boldsymbol{v}\in H^{k+1}(T)^{3}$. For all
$T\in\mathcal{T}_{h}$, using
$\uline{\boldsymbol{I}}_{T}^{k}\boldsymbol{v}=({\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v},({\boldsymbol{\pi}}_{F}^{k}\boldsymbol{v}_{|F})_{F\in\mathcal{F}_{T}})$
into (46b) for the first term and an integration by parts for the second term,
we obtain, for all
${\boldsymbol{\tau}}\in\mathcal{P}^{k}(\mathfrak{T}_{T})^{3\times 3}$,

$\displaystyle\int_{T}(\boldsymbol{G}_{\mathfrak{T}_{T}}^{l}\uline{\boldsymbol{I}}_{T}^{k}\boldsymbol{v}-\nabla\boldsymbol{v}):{\boldsymbol{\tau}}=$
$\displaystyle-\int_{T}({\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}-\boldsymbol{v})\cdot(\nabla\cdot{\boldsymbol{\tau}})+\sum_{F\in{\mathcal{F}}_{T}}\int_{F}({\boldsymbol{\pi}}_{F}^{k}\boldsymbol{v}-\boldsymbol{v})\cdot{\boldsymbol{\tau}}\boldsymbol{n}_{TF}$
(49) $\displaystyle+\sum_{\sigma\in\mathfrak{F}_{T}^{\rm
i}}\int_{\sigma}({\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}-\boldsymbol{v})\cdot\llbracket{\boldsymbol{\tau}}\rrbracket_{\sigma}\boldsymbol{n}_{\sigma}\eqcolon\mathfrak{T}_{1}+\mathfrak{T}_{2}+\mathfrak{T}_{3}.$

We now proceed to bound the terms in the right-hand side. Using Cauchy–Schwarz
and discrete inverse inequalities along with the approximation properties (5a)
of ${\boldsymbol{\pi}}_{T}^{k}$ with ${(l,r,m,s)}={(k,2,0,k+1)}$ we obtain,
for the first term,

$|\mathfrak{T}_{1}|\lesssim
h_{T}^{k}|\boldsymbol{v}|_{H^{k+1}(T)^{3}}\|{\boldsymbol{\tau}}\|_{L^{2}(T)^{3\times
3}}.$ (50)

For the second term, we use a Hölder inequality with exponents $(2,2,\infty)$
along with $\|\boldsymbol{n}_{TF}\|_{L^{\infty}(F)^{3}}=1$ to write

$|\mathfrak{T}_{2}|\leq\sum_{F\in\mathcal{F}_{T}}\|{\boldsymbol{\pi}}_{F}^{k}\boldsymbol{v}-\boldsymbol{v}\|_{L^{2}(F)^{3}}\|{\boldsymbol{\tau}}\|_{L^{2}(F)^{3\times
3}}\\\ \lesssim
h_{T}^{-\nicefrac{{1}}{{2}}}\left(\sum_{F\in\mathcal{F}_{T}}\|{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}-\boldsymbol{v}\|_{L^{2}(F)^{3}}^{2}\right)^{\nicefrac{{1}}{{2}}}\|{\boldsymbol{\tau}}\|_{L^{2}(T)^{3\times
3}}\lesssim
h_{T}^{k}|\boldsymbol{v}|_{H^{k+1}(T)^{3}}\|{\boldsymbol{\tau}}\|_{L^{2}(T)^{3\times
3}},$ (51)

where we have used the inequality (40) together with a discrete trace
inequality in the second step and the trace approximation properties (5b) of
${\boldsymbol{\pi}}_{T}^{k}$ with $(l,r,m,s)=(k,2,0,k+1)$ to conclude.
Let us now consider the third term in (49). Recalling the definition (45) of
the jump operator, we bound each integral over $\sigma\in\mathfrak{F}_{T}^{\rm
i}$ as follows:

$\displaystyle\left|\int_{\sigma}({\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}-\boldsymbol{v})\cdot\llbracket{\boldsymbol{\tau}}\rrbracket_{\sigma}\boldsymbol{n}_{\sigma}\right|$
$\displaystyle\leq\sum_{i=1}^{2}\left|\int_{\sigma}({\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}-\boldsymbol{v})\cdot\left({\boldsymbol{\tau}}_{|\tau_{i}}\boldsymbol{n}_{\tau_{i}\sigma}\right)\right|$
$\displaystyle\leq\sum_{i=1}^{2}\|{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}-\boldsymbol{v}\|_{L^{2}(\sigma)^{3}}\|{\boldsymbol{\tau}}_{|\tau_{i}}\|_{L^{2}(\sigma)^{3\times
3}}$
$\displaystyle\lesssim\sum_{i=1}^{2}\left(h_{\sigma}^{-1}\|{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}-\boldsymbol{v}\|_{L^{2}(\tau_{i})^{3}}+|{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}-\boldsymbol{v}|_{H^{1}(\tau_{i})^{3}}\right)\|{\boldsymbol{\tau}}\|_{L^{2}(\tau_{i})^{3\times
3}}$
$\displaystyle\lesssim\left(h_{T}^{-1}\|{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}-\boldsymbol{v}\|_{L^{2}(T)^{3}}+|{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{v}-\boldsymbol{v}|_{H^{1}(T)^{3}}\right)\|{\boldsymbol{\tau}}\|_{L^{2}(T)^{3\times
3}}$ $\displaystyle\lesssim
h_{T}^{k}|\boldsymbol{v}|_{H^{k+1}(T)^{3}}\|{\boldsymbol{\tau}}\|_{L^{2}(T)^{3\times
3}},$

where we have started with a triangle inequality, used Cauchy–Schwarz and
Hölder inequalities (the latter with exponents $(2,\infty)$) along with
$\|\boldsymbol{n}_{\tau_{i}\sigma}\|_{L^{\infty}(\sigma)^{3}}=1$ in the second
step, local continuous and discrete trace inequalities on the submesh for the
first and second factor, respectively, in the third step, and the fact that
$\tau_{i}\subset T$ for $i\in\\{1,2\\}$ along with the first geometric bound
in (3) and $h_{\sigma}^{-1}\lesssim h_{T}^{-1}$ (consequence of mesh
regularity) in the fourth step. The conclusion follows using the approximation
properties (5a) of ${\boldsymbol{\pi}}_{T}^{k}$ with $(l,r,m,s)=(k,2,0,k+1)$
for the first term in parenthesis and $(l,r,m,s)=(k,2,1,k+1)$ for the second
one. Gathering the above estimates and observing that
$\operatorname{card}(\mathfrak{F}_{T}^{\rm i})\leq
4\operatorname{card}(\mathfrak{T}_{T})\lesssim{1}$ by (3), we obtain

$|\mathfrak{T}_{3}|{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\lesssim}h_{T}^{k}|\boldsymbol{v}|_{H^{k+1}(T)^{3}}\|{\boldsymbol{\tau}}\|_{L^{2}(T)^{3\times
3}}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}.}\marginpar{\normalcolor\raggedright\scriptsize\textsf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}
R\\#1 minor point: Typo}}\@add@raggedright$ (52)

Setting
${\boldsymbol{\tau}}=\boldsymbol{G}_{\mathfrak{T}_{T}}^{l}\uline{\boldsymbol{I}}_{T}^{k}\boldsymbol{v}-{\boldsymbol{\pi}}^{l}_{\mathfrak{T}_{T}}\nabla\boldsymbol{v}$
in (49), using the bounds (50)–(52), and simplifying yields

$\|\boldsymbol{G}_{\mathfrak{T}_{T}}^{l}\uline{\boldsymbol{I}}_{T}^{k}{\boldsymbol{v}}-{\boldsymbol{\pi}}^{l}_{\mathfrak{T}_{T}}\nabla\boldsymbol{v}\|_{L^{2}(T)^{3\times
3}}\lesssim h_{T}^{k}|\boldsymbol{v}|_{H^{k+1}(T)^{3}}.$

Therefore, using a triangle inequality and the approximation properties (5a),
valid as well for ${\boldsymbol{\pi}}_{\tau}^{l}$, with
$(l,r,m,s)=(k,2,0,k+1)$ along with $h_{\tau}\leq h_{T}$, we infer

$\|\boldsymbol{G}_{\mathfrak{T}_{T}}^{l}\uline{\boldsymbol{I}}_{T}^{k}{\boldsymbol{v}}-\nabla\boldsymbol{v}\|_{L^{2}(T)^{3\times
3}}\leq\|\boldsymbol{G}_{\mathfrak{T}_{T}}^{l}\uline{\boldsymbol{I}}_{T}^{k}{\boldsymbol{v}}-{\boldsymbol{\pi}}^{l}_{\mathfrak{T}_{T}}\nabla\boldsymbol{v}\|_{L^{2}(T)^{3\times
3}}+\|{\boldsymbol{\pi}}^{l}_{\mathfrak{T}_{T}}\nabla\boldsymbol{v}-\nabla\boldsymbol{v}\|_{L^{2}(T)^{3\times
3}}\\\ \lesssim h_{T}^{k}|\boldsymbol{v}|_{H^{k+1}(T)^{3}}.\qed$

##  3 Discrete problem

###  3.1 Viscous term and pressure-velocity coupling

The viscous term and the pressure-velocity coupling are the same as in the
standard HHO method; see, e.g., [23, 8]. We briefly recall them here to make
the exposition self-contained.

The viscous bilinear form $a_{h}$:
$\uline{\boldsymbol{U}}_{h}^{k}\times\uline{\boldsymbol{U}}_{h}^{k}\rightarrow\mathbb{R}$
is such that, for all
$\uline{\boldsymbol{w}}_{h},\uline{\boldsymbol{v}}_{h},\in\uline{\boldsymbol{U}}_{h}^{k}$,

$a_{h}(\uline{\boldsymbol{w}}_{h},\uline{\boldsymbol{v}}_{h})\coloneq\sum_{T\in{\mathcal{T}}_{h}}\left(\int_{T}\boldsymbol{G}_{T}^{k}\uline{\boldsymbol{w}}_{T}:\boldsymbol{G}_{T}^{k}\uline{\boldsymbol{v}}_{T}+s_{T}(\uline{\boldsymbol{w}}_{T},\uline{\boldsymbol{v}}_{T})\right),$

where, for any $T\in{\mathcal{T}}_{h}$,
$s_{T}:\uline{\boldsymbol{U}}_{T}^{k}\times\uline{\boldsymbol{U}}_{T}^{k}\rightarrow\mathbb{R}$
denotes a local stabilization bilinear form designed according to the
principles of [16, Assumption 2.4], so that, in particular, there exists
$C_{a}>0$ independent of $h$ (and, clearly, also of $\nu$ and $\lambda$) such
that, for all $\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}_{h}^{k}$,

$C_{a}\|\uline{\boldsymbol{v}}_{h}\|_{1,h}^{2}\leq
a_{h}(\uline{\boldsymbol{v}}_{h},\uline{\boldsymbol{v}}_{h})\leq
C_{a}^{-1}\|\uline{\boldsymbol{v}}_{h}\|_{1,h}^{2}.$ (53)

Recalling the definition (10) of the local divergence $D_{T}^{k}$, the global
pressure-velocity coupling bilinear form
$b_{h}:\uline{\boldsymbol{U}}_{h,0}^{k}\times\mathcal{P}^{k}({\mathcal{T}}_{h})\rightarrow\mathbb{R}$
is such that, for all
$(\uline{\boldsymbol{v}}_{h},q_{h})\in\uline{\boldsymbol{U}}_{h,0}^{k}\times\mathcal{P}^{k}({\mathcal{T}}_{h})$,

$b_{h}(\uline{\boldsymbol{v}}_{h},q_{h}):=-\sum_{T\in{\mathcal{T}}_{h}}\int_{T}D_{T}^{k}\uline{\boldsymbol{v}}_{T}\,q_{T},$

where $q_{T}\coloneq q_{h|T}$. The properties of $b_{h}$ relevant for the
analysis can be found in [16, Lemma 8.12].

###  3.2 Body force

The discretization of the body force leverages the new divergence-preserving
velocity reconstruction introduced in Section 2.4. Specifically, we introduce
the bilinear form
$\ell_{h}:{L}^{2}({\Omega})^{3}\times\uline{\boldsymbol{U}}_{h}^{k}\rightarrow\mathbb{R}$
such that, for any $\boldsymbol{\phi}\in{L}^{2}({\Omega})^{3}$ and any
$\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}_{h}^{k}$,

$\ell_{h}(\boldsymbol{\phi},\uline{\boldsymbol{v}}_{h})\coloneq\int_{\Omega}\boldsymbol{\phi}\cdot\boldsymbol{R}_{h}^{k}\underline{\boldsymbol{v}}_{h}.$

######  Lemma 8 (Properties of $\ell_{h}$).

The bilinear form $\ell_{h}$ has the following properties:

  1. (i)

_Velocity invariance._ Recalling the Hodge decomposition (2) of
$\boldsymbol{f}$, it holds

$\ell_{h}(\boldsymbol{g}+\lambda\nabla\psi,\uline{\boldsymbol{v}}_{h})=\ell_{h}(\boldsymbol{g},\uline{\boldsymbol{v}}_{h})+b_{h}(\uline{\boldsymbol{v}}_{h},\lambda\pi_{h}^{k}\psi)\qquad\forall\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}_{h,0}^{k}.$
(54) ††margin:  Referees #1 and #2 remark: Now we have a body force
consistency error for arbitrary $k\geq 0$

  2. (ii)

_Consistency._ For all $\boldsymbol{\phi}\in{L}^{2}({\Omega})^{3}\cap
H^{k}({\mathcal{T}}_{h})^{3}$,

$\|{\cal{E}}_{\ell,h}(\boldsymbol{\phi};\cdot)\|_{1,h,*}\lesssim
h^{k+1}|\boldsymbol{\phi}|_{H^{k}({\mathcal{T}}_{h})^{3}},$ (55)

where the linear form
${\cal{E}}_{\ell,h}(\boldsymbol{\phi};\cdot):\uline{\boldsymbol{U}}_{h}^{k}\rightarrow\mathbb{R}$,
representing the body force consistency error, is such that

$\displaystyle{\cal{E}}_{\ell,h}(\boldsymbol{\phi};\uline{\boldsymbol{v}}_{h})\coloneq\ell_{h}(\boldsymbol{\phi},\uline{\boldsymbol{v}}_{h})-\int_{\Omega}\boldsymbol{\phi}\cdot\boldsymbol{v}_{h}=\sum_{T\in\mathcal{T}_{h}}\int_{T}{\boldsymbol{\phi}}\cdot(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}-{\boldsymbol{v}}_{T}).$
(56)

###### Proof.

_(i) Velocity invariance._ The proof is the same as in [11, Section 4.3],
using the fact that
$(\nabla\cdot\boldsymbol{R}_{h}^{k}\underline{\boldsymbol{v}}_{h})_{|T}=D_{T}^{k}\uline{\boldsymbol{v}}_{T}\in\mathcal{P}^{k}(T)$,
which is enforced by (14b).
_(ii) Consistency._ We prove the cases $k=0$ and $k\geq 1$ separately.
(ii.A) _The case $k=0$._ Taking absolute values in (56) and using
Cauchy–Schwarz inequalities along with (15) and $h_{T}\leq h$ for all
$T\in\mathcal{T}_{h}$, we can write
$\left|{\cal{E}}_{\ell,h}(\boldsymbol{\phi};\uline{\boldsymbol{v}}_{h})\right|\leq
h\|{\boldsymbol{\phi}}\|_{L^{2}(\Omega)^{3}}\|\uline{\boldsymbol{v}}_{T}\|_{1,h}$.
Passing to the supremum over
$\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}_{h}^{k}$ such that
$\|\uline{\boldsymbol{v}}_{h}\|_{1,h}=1$, we obtain (55).
(ii.B) _The case $k\geq 1$._ Using (17) in (56) and continuing with
Cauchy–Schwarz inequalities, we obtain

$\displaystyle\left|{\cal{E}}_{\ell,h}({\boldsymbol{\phi}};\uline{\boldsymbol{v}}_{h})\right|=\left|\sum_{T\in\mathcal{T}_{h}}\int_{T}({\boldsymbol{\phi}}-{\boldsymbol{\pi}}^{k-1}{\boldsymbol{\phi}})\cdot(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}-{\boldsymbol{v}}_{T})\right|\leq\sum_{T\in\mathcal{T}_{h}}\|{\boldsymbol{\phi}}-{\boldsymbol{\pi}}^{k-1}{\boldsymbol{\phi}}\|_{L^{2}(T)^{3}}\|\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}-{\boldsymbol{v}}_{T}\|_{L^{2}(T)^{3}}.$

Using then the approximation properties (5a) of the $L^{2}$-projector with
$(l,m,r,s)=(k-1,0,2,k)$ for the first factor and the bound (15) for the
second, applying discrete Cauchy–Schwarz inequalities to the sums, and passing
to the supremum over
$\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}_{h}^{k}$ such that
$\|\uline{\boldsymbol{v}}_{h}\|_{1,h}=1$, (55) follows. ∎

###  3.3 Convective term

To discretize the convective term, we introduce the global trilinear form
$t_{h}:\left[\uline{\boldsymbol{U}}_{h}^{k}\right]^{3}\to\mathbb{R}$ such that

$t_{h}(\uline{\boldsymbol{w}}_{h},\uline{\boldsymbol{v}}_{h},\uline{\boldsymbol{z}}_{h})\coloneq\sum_{T\in{\mathcal{T}}_{h}}t_{T}(\uline{\boldsymbol{w}}_{T},\uline{\boldsymbol{v}}_{T},\uline{\boldsymbol{z}}_{T}),$
(57a) where, for any $T\in{\mathcal{T}}_{h}$,
$t_{T}:\left[\uline{\boldsymbol{U}}_{T}^{k}\right]^{3}\rightarrow\mathbb{R}$
is defined as
$t_{T}(\uline{\boldsymbol{w}}_{T},\uline{\boldsymbol{v}}_{T},\uline{\boldsymbol{z}}_{T})\coloneq\int_{T}\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}\uline{\boldsymbol{w}}_{T}\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}\cdot\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}-\int_{T}\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}\uline{\boldsymbol{w}}_{T}\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\cdot\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}.$
(57b)

######  Remark 9 (Reformulation of $t_{h}$).

In practice, it is not necessary to compute the piecewise gradient
reconstruction operators $\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}$ to
evaluate $t_{T}$ and $t_{h}$. As a matter of fact, expanding the piecewise
gradient operator in (57) according to its definition (46a), we have that

$\displaystyle
t_{h}(\uline{\boldsymbol{w}}_{h},\uline{\boldsymbol{v}}_{h},\uline{\boldsymbol{z}}_{h})=$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\left[\int_{T}\nabla{\boldsymbol{w}}_{T}\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}\cdot\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}-\int_{T}\nabla{\boldsymbol{w}}_{T}\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\cdot\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}\right]$
$\displaystyle+\sum_{T\in{\mathcal{T}}_{h}}\sum_{F\in{\mathcal{F}}_{T}}\int_{F}(\boldsymbol{w}_{F}-\boldsymbol{w}_{T})\cdot\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}\cdot\boldsymbol{n}_{TF})$
$\displaystyle-\sum_{T\in{\mathcal{T}}_{h}}\sum_{F\in{\mathcal{F}}_{T}}\int_{F}(\boldsymbol{w}_{F}-\boldsymbol{w}_{T})\cdot\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{v}}_{T}(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\cdot\boldsymbol{n}_{TF}).$

The properties of $t_{h}$ relevant for the analysis are contained in the
following lemma.

######  Lemma 10 (Properties of $t_{h}$).

The trilinear form $t_{h}$ has the following properties:

  1. 1.

_Non-dissipativity._ For all
$\uline{\boldsymbol{w}}_{h},\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}_{h}^{k}$,
it holds that

$t_{h}(\uline{\boldsymbol{w}}_{h},\uline{\boldsymbol{v}}_{h},\uline{\boldsymbol{v}}_{h})=0.$
(58)

  2. 2.

_Boundedness._ There exists a real number $C_{t}>0$ independent of $h$ (and,
clearly, also of $\nu$ and $\lambda$) such that, for all
$\uline{\boldsymbol{w}}_{h},\uline{\boldsymbol{v}}_{h},\uline{\boldsymbol{z}}_{h}\in\uline{\boldsymbol{U}}_{h}^{k}$,

$|t_{h}(\uline{\boldsymbol{w}}_{h},\uline{\boldsymbol{v}}_{h},\uline{\boldsymbol{z}}_{h})|\leq
C_{t}\|\uline{\boldsymbol{w}}_{h}\|_{1,h}\|\uline{\boldsymbol{v}}_{h}\|_{1,h}\|\uline{\boldsymbol{z}}_{h}\|_{1,h}.$
(59)

  3. 3.

_Consistency._ It holds, for all $\boldsymbol{w}\in\boldsymbol{U}\cap
W^{k+1,4}({\mathcal{T}}_{h})^{3}$ such that $\nabla\cdot\boldsymbol{w}=0$ a.e.
in $\Omega$,

$\|{\cal{E}}_{t,h}(\boldsymbol{w};\cdot)\|_{1,h,*}\lesssim
h^{k+1}|\boldsymbol{w}|_{W^{k+1,4}({\mathcal{T}}_{h})^{3}}\|\boldsymbol{w}\|_{W^{1,4}(\Omega)^{3}},$
(60)

where the linear form
${\cal{E}}_{t,h}(\boldsymbol{w};\cdot):\uline{\boldsymbol{U}}_{h}^{k}\rightarrow\mathbb{R}$
representing the consistency error is such that, for all
$\uline{\boldsymbol{z}}_{h}\in\uline{\boldsymbol{U}}_{h}^{k}$,

${\cal{E}}_{t,h}(\boldsymbol{w};\uline{\boldsymbol{z}}_{h})\coloneq\ell_{h}((\nabla\times\boldsymbol{w})\times\boldsymbol{w},\uline{\boldsymbol{z}}_{h})-t_{h}(\uline{\boldsymbol{I}}_{h}^{k}\boldsymbol{w},\uline{\boldsymbol{I}}_{h}^{k}\boldsymbol{w},\uline{\boldsymbol{z}}_{h}).$

###### Proof.

(i) _Non-dissipativity._ Immediate consequence of the definition (57) of
$t_{h}$.
(ii) _Boundedness._ The proof is similar to that of [11, Lemma 7.ii] using the
Hölder inequalities with exponent (2,4,4), the bound
$\eqref{eq:disc:grad:bound}$, and the discrete Sobolev embedding (41) with
$r=4$. The details are omitted for the sake of conciseness.
(iii) _Consistency._ Let
$\uline{\hat{\boldsymbol{w}}}_{h}\coloneq\uline{\boldsymbol{I}}_{h}^{k}\boldsymbol{w}$.
Proceeding as in [11, Lemma 7.iii], we obtain the following decomposition:

$\displaystyle{\cal{E}}_{t,h}(\boldsymbol{w};\uline{\boldsymbol{z}}_{h})$
$\displaystyle=\underbrace{\sum_{T\in{\mathcal{T}}_{h}}\int_{T}(\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}\uline{\hat{\boldsymbol{w}}}_{T}-\nabla\boldsymbol{w})\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\cdot\boldsymbol{w}}_{\mathfrak{T}_{1}}+\underbrace{\sum_{T\in{\mathcal{T}}_{h}}\int_{T}(\nabla\boldsymbol{w}-\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}\uline{\hat{\boldsymbol{w}}}_{T})\boldsymbol{w}\cdot\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}}_{\mathfrak{T}_{2}}$
(61)
$\displaystyle\quad+\underbrace{\sum_{T\in{\mathcal{T}}_{h}}\int_{T}\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}\uline{\hat{\boldsymbol{w}}}_{T}(\boldsymbol{w}-\boldsymbol{R}_{T}^{k}\uline{\hat{\boldsymbol{w}}}_{T})\cdot\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}}_{\mathfrak{T}_{3}}+\underbrace{\sum_{T\in{\mathcal{T}}_{h}}\int_{T}\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}\uline{\hat{\boldsymbol{w}}}_{T}\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\cdot(\boldsymbol{R}_{T}^{k}\uline{\hat{\boldsymbol{w}}}_{T}-\boldsymbol{w}).}_{\mathfrak{T}_{4}}$

We next proceed to estimate the terms
$\mathfrak{T}_{1},\cdots,\mathfrak{T}_{4}$.
(iii.A) _Estimate of $\mathfrak{T}_{1}$._ Following similar steps as in [11,
Lemma 7.iii.A] using the approximation properties (48) of
$\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}$ and its definition (46), we get
that

$|\mathfrak{T}_{1}|{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\lesssim}\marginpar{\normalcolor\raggedright\scriptsize\textsf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}
R\\#2 point:
Typo}}h^{k+1}|\boldsymbol{w}|_{H^{k+1}({\cal{T}}_{h})^{3}}\|\boldsymbol{w}\|_{W^{1,4}(\Omega)^{3}}\|\uline{\boldsymbol{z}}_{h}\|_{1,h}.\@add@raggedright$
(62)

(iii.B) _Estimate of $\mathfrak{T}_{2}$._ For the term $\mathfrak{T}_{2}$ in
(61), inserting $\pm{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w}$ into the second
factor, we get

$\displaystyle\mathfrak{T}_{2}$
$\displaystyle=\sum_{T\in{\mathcal{T}}_{h}}\int_{T}(\nabla\boldsymbol{w}-\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}\uline{\hat{\boldsymbol{w}}}_{T})(\boldsymbol{w}-{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w})\cdot\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}+\sum_{T\in{\mathcal{T}}_{h}}\int_{T}(\nabla\boldsymbol{w}-\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}\uline{\hat{\boldsymbol{w}}}_{T}){\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w}\cdot\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}$
(63) $\displaystyle\eqcolon\mathfrak{T}_{2,1}+\mathfrak{T}_{2,2}.$

We bound $\mathfrak{T}_{2,1}$ using Hölder inequalities with exponents
$(2,4,4)$, then the approximation properties (48) of
$\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}$ and (5a) of
${\boldsymbol{\pi}}_{T}^{0}$ with $(l,m,r,s)=(0,0,4,1)$, and the bound (41)
with $r=4$:

$|\mathfrak{T}_{2,1}|\lesssim
h^{k+1}|\boldsymbol{w}|_{H^{k+1}({\mathcal{T}}_{h})^{3}}|\boldsymbol{w}|_{W^{1,4}(\Omega)^{3}}\|\uline{\boldsymbol{z}}_{T}\|_{1,h}.$
(64)

To estimate $\mathfrak{T}_{2,2}$ in (63), we integrate by parts the term
involving $\nabla\boldsymbol{w}$ and we use, for each element
$T\in{\mathcal{T}}_{h}$, the definition (46b) of
$\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}$ with
$(\uline{\boldsymbol{v}}_{T},{\boldsymbol{\tau}})=(\uline{\hat{\boldsymbol{w}}}_{T},\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\otimes{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w})$
(notice that
$\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\otimes{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w}\in\mathcal{P}^{k+1}(\mathfrak{T}_{T})^{3\times
3}\subset\mathcal{P}^{2(k+1)}(\mathfrak{T}_{T})^{3\times 3}$) to write

$\displaystyle\mathfrak{T}_{2,2}$
$\displaystyle=-\sum_{T\in{\mathcal{T}}_{h}}\sum_{\tau\in\mathfrak{T}_{T}}\int_{\tau}(\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T}^{k}\boldsymbol{w})\cdot\nabla\cdot(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\otimes{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w})$
(65)
$\displaystyle\quad+\sum_{T\in\mathcal{T}_{h}}\sum_{\tau\in\mathfrak{T}_{T}}\sum_{\sigma\in\mathfrak{F}_{T}^{\rm
i}}\int_{\sigma}(\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T}^{k}\boldsymbol{w})\cdot\llbracket\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\otimes{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w}\rrbracket_{\sigma}\boldsymbol{n}_{\sigma}$
$\displaystyle\quad+\sum_{T\in{\mathcal{T}}_{h}}\sum_{F\in{\mathcal{F}}_{T}}\int_{F}(\boldsymbol{w}-{{\boldsymbol{\pi}}}_{F}^{k}\boldsymbol{w})\cdot(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\otimes{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w})\boldsymbol{n}_{TF},$
$\displaystyle\eqcolon\mathfrak{T}_{2,2,1}+\mathfrak{T}_{2,2,2}+\mathfrak{T}_{2,2,3}.$

For $\mathfrak{T}_{2,2,1}$, we first observe that
$\nabla\cdot(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\otimes{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w})=\nabla\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w}+\cancel{\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}(\nabla\cdot{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w})}\in\mathcal{P}^{k}(\mathfrak{T}_{T})^{3}$.
Hence, using Hölder inequalities with exponents $(4,2,4)$, we infer that

$\displaystyle|\mathfrak{T}_{2,2,1}|$
$\displaystyle\leq\sum_{T\in{\mathcal{T}}_{h}}\sum_{\tau\in\mathfrak{T}_{T}}\|\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T}^{k}\boldsymbol{w}\|_{L^{4}(\tau)^{3}}\|\nabla\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\|_{L^{2}(\tau)^{3\times
3}}\|{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w}\|_{L^{4}(\tau)^{3}}$
$\displaystyle\leq\sum_{T\in{\mathcal{T}}_{h}}\|\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T}^{k}\boldsymbol{w}\|_{L^{4}(T)^{3}}\|{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w}\|_{L^{4}(T)^{3}}\sum_{\tau\in\mathfrak{T}_{T}}\|\nabla\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\|_{L^{2}(\tau)^{3\times
3}}$
$\displaystyle\lesssim\sum_{T\in{\mathcal{T}}_{h}}h_{T}^{k+1}|\boldsymbol{w}|_{W^{k+1,4}(T)^{3}}\|\boldsymbol{w}\|_{W^{1,4}(T)^{3}}\sum_{\tau\in\mathfrak{T}_{T}}\|\nabla\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\|_{L^{2}(\tau)^{3\times
3}},$ (66)

where, in the second step, we have used the fact that $\tau\subset T$ for all
$\tau\in\mathfrak{T}_{T}$, while, in the third step, we have used the
approximation properties (5a) of the $L^{2}$-orthogonal projector with
$(l,m,r,s)=(k,0,4,k+1)$ for the first factor and its boundedness for the
second factor. To bound
$\|\nabla\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\|_{L^{2}(\tau)^{3\times
3}}$, we first observe that
$(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T})_{|\tau}$ is in the space
$\mathcal{P}^{k+1}(\tau)^{3}$, and it holds that

$\displaystyle\|\nabla\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}\|_{L^{2}(\tau)^{3\times
3}}$
$\displaystyle\leq\|\nabla(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}-\boldsymbol{z}_{T})\|_{L^{2}(\tau)^{3}}+\|\nabla\boldsymbol{z}_{T}\|_{L^{2}(\tau)^{3}}$
$\displaystyle\lesssim
h_{\tau}^{-1}\|\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}-\boldsymbol{z}_{T}\|_{L^{2}(\tau)^{3}}+\|\nabla\boldsymbol{z}_{T}\|_{L^{2}(\tau)^{3}}$
$\displaystyle\lesssim
h_{\tau}^{-1}\|\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}-\boldsymbol{z}_{T}\|_{L^{2}(T)^{3}}+\|\nabla\boldsymbol{z}_{T}\|_{L^{2}(T)^{3}}$
$\displaystyle\lesssim
h_{\tau}^{-1}h_{T}|\uline{\boldsymbol{z}}_{T}|_{1,\partial
T}+\|\nabla\boldsymbol{z}_{T}\|_{L^{2}(T)^{3}}\lesssim\|\uline{\boldsymbol{z}}_{T}\|_{1,T},$

where we have started with a triangle inequality after inserting
$\pm\nabla\boldsymbol{z}_{T}$, used a local discrete inverse inequality on
$\tau$ in the second step, the fact that $\tau\subset T$ for all
$\tau\in\mathfrak{T}_{T}$ in the third step, the bound (15) in the fourth
step, and the inequality $h_{\tau}^{-1}h_{T}\lesssim 1$ valid for regular mesh
sequences (see [16, Eq. (1.4)]), along with the definition (7) of the
$\|{\cdot}\|_{1,T}$-norm to conclude. Plugging this last inequality into (66)
and using the geometric bound (3) on $\mathfrak{T}_{T}$ along with a discrete
Hölder inequality, we arrive at

$|\mathfrak{T}_{2,2,1}|\lesssim
h^{k+1}|\boldsymbol{w}|_{W^{k+1,4}({\mathcal{T}}_{h})^{3}}\|\boldsymbol{w}\|_{W^{1,4}(\Omega)^{3}}\|\uline{\boldsymbol{z}}_{h}\|_{1,h}.$
(67)

To estimate $\mathfrak{T}_{2,2,2}$ in (65), we insert $\pm\boldsymbol{z}_{T}$
into the first factor inside the jump operator to write

$\displaystyle\mathfrak{T}_{2,2,2}=$
$\displaystyle\sum_{T\in\mathcal{T}_{h}}\sum_{\tau\in\mathfrak{T}_{T}}\sum_{\sigma\in\mathfrak{F}_{T}^{\rm
i}}\int_{\sigma}(\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T}^{k}\boldsymbol{w})\cdot\llbracket(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}-\boldsymbol{z}_{T})\otimes{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w}\rrbracket_{\sigma}\boldsymbol{n}_{\sigma}$
$\displaystyle-\sum_{T\in\mathcal{T}_{h}}\sum_{\tau\in\mathfrak{T}_{T}}\sum_{\sigma\in\mathfrak{F}_{T}^{\rm
i}}\int_{\sigma}\cancel{(\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T}^{k}\boldsymbol{w})\cdot\llbracket\boldsymbol{z}_{T}\otimes{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w}\rrbracket_{\sigma}\boldsymbol{n}_{\sigma}},$

where the second addend cancels since
$\boldsymbol{z}_{T}\otimes{\boldsymbol{\pi}}_{T}^{0}{\boldsymbol{w}}$ is
continuous across the interior faces of $\mathfrak{T}_{T}$. Setting
$\mathfrak{F}_{h}^{\rm i}\coloneq\\{\sigma\in\mathfrak{F}_{T}^{\rm
i}:T\in\mathcal{T}_{h}\\}$ and exchanging the order of the sums, we can now
express $\mathfrak{T}_{2,2,2}$ in the following equivalent form:

$\displaystyle\mathfrak{T}_{2,2,2}=$
$\displaystyle\sum_{\sigma\in\mathfrak{F}_{h}^{\rm
i}}\sum_{i=1}^{2}\int_{\sigma}(\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T_{\sigma}}^{k}\boldsymbol{w})\cdot\left[(\boldsymbol{R}_{T_{\sigma}}^{k}\uline{\boldsymbol{z}}_{T_{\sigma}}-\boldsymbol{z}_{T_{\sigma}})\otimes{\boldsymbol{\pi}}_{T_{\sigma}}^{0}\boldsymbol{w}\right]_{|\tau_{i}}\boldsymbol{n}_{\tau_{i}\sigma},$
(68)

where, for a given $\sigma\in\mathfrak{F}_{h}^{\rm i}$,
$T_{\sigma}\in\mathcal{T}_{h}$ is the element in which $\sigma$ is contained,
while $\tau_{1}$ and $\tau_{2}$ denote the simplices in
$\mathfrak{T}_{T_{\sigma}}$ sharing $\sigma$. To bound the right-hand side of
the above expression, we apply Hölder inequalities with exponents
$(4,2,4,\infty)$ along with
$\|\boldsymbol{n}_{\tau_{i}\sigma}\|_{L^{\infty}(\sigma)^{3}}=1$ to write

$\displaystyle|\mathfrak{T}_{2,2,2}|$
$\displaystyle\leq\sum_{\sigma\in\mathfrak{F}_{h}^{\rm
i}}\sum_{i=1}^{2}\|\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T_{\sigma}}^{k}\boldsymbol{w}\|_{L^{4}(\sigma)^{3}}\|(\boldsymbol{R}_{T_{\sigma}}^{k}\uline{\boldsymbol{z}}_{T_{\sigma}})_{|\tau_{i}}-\boldsymbol{z}_{{T_{\sigma}}}\|_{L^{2}(\sigma)^{3}}\|{\boldsymbol{\pi}}_{T_{\sigma}}^{0}\boldsymbol{w}\|_{L^{4}(\sigma)^{3}}$
(69) $\displaystyle\lesssim\sum_{\sigma\in\mathfrak{F}_{h}^{\rm
i}}\sum_{i=1}^{2}\|\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T_{\sigma}}^{k}\boldsymbol{w}\|_{L^{4}(\sigma)^{3}}h_{\tau_{i}}^{-\frac{1}{2}}\|\boldsymbol{R}_{T_{\sigma}}^{k}\uline{\boldsymbol{z}}_{T_{\sigma}}-\boldsymbol{z}_{T_{\sigma}}\|_{L^{2}(\tau_{i})^{3}}h_{\tau_{i}}^{-\frac{1}{4}}\|{\boldsymbol{\pi}}_{T_{\sigma}}^{0}\boldsymbol{w}\|_{L^{4}(\tau_{i})^{3}}$
$\displaystyle\lesssim\sum_{\sigma\in\mathfrak{F}_{h}^{\rm
i}}\left(h_{T_{\sigma}}^{\frac{1}{4}}\|\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T_{\sigma}}^{k}\boldsymbol{w}\|_{L^{4}(\sigma)^{3}}\right)\|\uline{\boldsymbol{z}}_{T_{\sigma}}\|_{1,T_{\sigma}}\|\boldsymbol{w}\|_{W^{1,4}(T_{\sigma})^{3}},$

where, in the second step, we have used local trace inequalities on the
submesh for the second and third factors while, in the third step, we have
used $\tau_{i}\subset T_{\sigma}$ along with (15), (7), and
$h_{\tau_{i}}^{-1}h_{T_{\sigma}}\lesssim 1$ (consequence of mesh regularity)
for the second factor while, for the third factor, we have used again
$\tau_{i}\subset T_{\sigma}$ along with the boundedness of the
$L^{2}$-orthogonal projector. Using trace inequalities on the submesh along
with the approximation properties of the $L^{2}$-orthogonal projector, we
infer
$h_{T_{\sigma}}^{\frac{1}{4}}\|\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T_{\sigma}}^{k}\boldsymbol{w}\|_{L^{4}(\sigma)^{3}}\lesssim
h_{T_{\sigma}}^{k+1}|\boldsymbol{w}|_{W^{k+1,4}(T_{\sigma})^{3}}$ which,
plugged into (69) and combined with the geometric bound (3), gives

$\displaystyle|\mathfrak{T}_{2,2,2}|\lesssim
h^{k+1}|\boldsymbol{w}|_{W^{k+1,4}({\mathcal{T}}_{h})^{3}}\|\boldsymbol{w}\|_{W^{1,4}(\Omega)^{3}}\|\uline{\boldsymbol{z}}_{h}\|_{1,h}.$
(70)

To bound the term $\mathfrak{T}_{2,2,3}$ in (65), we first insert
$\pm{\boldsymbol{z}}_{T}$ into its second factor to write

$\displaystyle\mathfrak{T}_{2,2,3}=$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\sum_{F\in{\mathcal{F}}_{T}}\int_{F}(\boldsymbol{w}-{{\boldsymbol{\pi}}}_{F}^{k}\boldsymbol{w})\cdot\left[(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}-{\boldsymbol{z}}_{T})\otimes{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w}\right]\boldsymbol{n}_{TF}$
$\displaystyle+\sum_{T\in{\mathcal{T}}_{h}}\sum_{F\in{\mathcal{F}}_{T}}\int_{F}\cancel{(\boldsymbol{w}-{{\boldsymbol{\pi}}}_{F}^{k}\boldsymbol{w})\cdot(\boldsymbol{z}_{T}\otimes{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w})\boldsymbol{n}_{TF}},$

where the second addend cancels by the definition (4) of the
$L^{2}$-orthogonal projector ${\boldsymbol{\pi}}_{F}^{k}$ since
$(\boldsymbol{z}_{T}\otimes{\boldsymbol{\pi}}_{T}^{0})_{|F}\boldsymbol{n}_{TF}\in\mathcal{P}^{k}(F)^{3}$.
Now, we rewrite the equation above as

$\displaystyle\mathfrak{T}_{2,2,3}=$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\sum_{F\in{\mathcal{F}}_{T}}\sum_{\sigma\in\mathfrak{F}_{F}}\int_{\sigma}(\boldsymbol{w}-{{\boldsymbol{\pi}}}_{F}^{k}\boldsymbol{w})\cdot\left[(\boldsymbol{R}_{T}^{k}\uline{\boldsymbol{z}}_{T}-{\boldsymbol{z}}_{T})\otimes{\boldsymbol{\pi}}_{T}^{0}\boldsymbol{w}\right]\boldsymbol{n}_{TF},$

thus, using a similar procedure as for (68)–(69), but for
$\sigma\in\mathfrak{F}_{F}$ and $\tau_{i}=\tau_{\sigma}$ where
$\tau_{\sigma}\in\mathfrak{T}_{T}$ is the simplicial subelement containing
$\sigma$, we infer that

$|\mathfrak{T}_{2,2,3}|\lesssim\sum_{T\in{\mathcal{T}}_{h}}\sum_{F\in{\mathcal{F}}_{T}}\sum_{\sigma\in\mathfrak{F}_{F}}\left(h_{T}^{\frac{1}{4}}\|\boldsymbol{w}-{{\boldsymbol{\pi}}}_{F}^{k}\boldsymbol{w}\|_{L^{4}(\sigma)^{3}}\right)\|\uline{\boldsymbol{z}}_{T}\|_{1,T}\|\boldsymbol{w}\|_{W^{1,4}(T)^{3}},$

in addition, using the fact that

$\displaystyle\|\boldsymbol{w}-{{\boldsymbol{\pi}}}_{F}^{k}\boldsymbol{w}\|_{L^{4}(\sigma)^{3}}$
$\displaystyle\leq\|\boldsymbol{w}-{{\boldsymbol{\pi}}}_{F}^{k}\boldsymbol{w}\|_{L^{4}(F)^{3}}$
$\displaystyle\leq\|{{\boldsymbol{\pi}}}_{T}^{k}\boldsymbol{w}-{{\boldsymbol{\pi}}}_{F}^{k}\boldsymbol{w}\|_{L^{4}(F)^{3}}+\|\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T}^{k}\boldsymbol{w}\|_{L^{4}(F)^{3}}$
$\displaystyle\leq\|{{\boldsymbol{\pi}}}_{F}^{k}({{\boldsymbol{\pi}}}_{T}^{k}\boldsymbol{w}-\boldsymbol{w})\|_{L^{4}(F)^{3}}+\|\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T}^{k}\boldsymbol{w}\|_{L^{4}(F)^{3}}\lesssim\|\boldsymbol{w}-{{\boldsymbol{\pi}}}_{T}^{k}\boldsymbol{w}\|_{L^{4}(F)^{3}},$

the approximation properties (5b) of the $L^{2}$-orthogonal projector with
$(l,m,r,s)=(k,0,4,k+1)$, and the bound (3), we obtain

$|\mathfrak{T}_{2,2,3}|\lesssim
h^{k+1}|\boldsymbol{w}|_{W^{k+1,4}({\mathcal{T}}_{h})^{3}}\|\boldsymbol{w}\|_{W^{1,4}(\Omega)^{3}}\|\uline{\boldsymbol{z}}_{h}\|_{1,h}.$
(71)

Plugging the estimates (67), (70), and (71) into (65), and, combining the
resulting estimate with (64), we finally obtain

$|\mathfrak{T}_{2}|\lesssim
h^{k+1}|\boldsymbol{w}|_{W^{k+1,4}({\mathcal{T}}_{h})^{3}}\|\boldsymbol{w}\|_{W^{1,4}(\Omega)^{3}}\|\uline{\boldsymbol{z}}_{h}\|_{1,h}.$
(72)

(iii.C) _Estimate of $\mathfrak{T}_{3}$ and $\mathfrak{T}_{4}$._ To bound
$\mathfrak{T}_{3}$, we follow the same steps as in [11, Lemma 7.iii.C] along
with the boundedness (47) of $\boldsymbol{G}_{\mathfrak{T}_{T}}^{2(k+1)}$ to
obtain

$|\mathfrak{T}_{3}|\lesssim|\boldsymbol{w}|_{H^{1}(\Omega)^{3}}\left(\sum_{T\in{\mathcal{T}}_{h}}\|\boldsymbol{w}-\boldsymbol{R}_{T}^{k}\uline{\hat{\boldsymbol{w}}}_{T}\|_{L^{4}(T)^{3}}^{4}\right)^{\frac{1}{4}}\|\uline{\boldsymbol{z}}_{h}\|_{1,h}.$

To estimate each addend in the second factor, we first insert
$\pm{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{w}$ and then use a triangle
inequality to write

$\displaystyle\|\boldsymbol{w}-\boldsymbol{R}_{T}^{k}\uline{\hat{\boldsymbol{w}}}_{T}\|_{L^{4}(T)^{3}}$
$\displaystyle\leq\|\boldsymbol{w}-{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{w}\|_{L^{4}(T)^{3}}+\|{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{w}-\boldsymbol{R}_{T}^{k}\uline{\hat{\boldsymbol{w}}}_{T}\|_{L^{4}(T)^{3}}$
(73) $\displaystyle\lesssim
h_{T}^{k+1}|\boldsymbol{w}|_{W^{k+1,4}(T)^{3}}+\|{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{w}-\boldsymbol{R}_{T}^{k}\uline{\hat{\boldsymbol{w}}}_{T}\|_{L^{4}(T)^{3}},$

where we have used the approximation properties (5a) of
${\boldsymbol{\pi}}_{T}^{k}$ with $(l,m,r,s)=(k,0,4,k+1)$ to conclude. To
estimate the second term in the right-hand side of (73), we proceed as
follows:

$\displaystyle\|{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{w}-\boldsymbol{R}_{T}^{k}\uline{\hat{\boldsymbol{w}}}_{T}\|_{L^{4}(T)^{3}}^{4}$
$\displaystyle=\sum_{\tau\in\mathfrak{T}_{T}}\|{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{w}-\boldsymbol{R}_{T}^{k}\uline{\hat{\boldsymbol{w}}}_{T}\|_{L^{4}(\tau)^{3}}^{4}$
$\displaystyle\lesssim\sum_{\tau\in\mathfrak{T}_{T}}\left(h_{\tau}^{-\frac{3}{4}}\|{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{w}-\boldsymbol{R}_{T}^{k}\uline{\hat{\boldsymbol{w}}}_{T}\|_{L^{2}(\tau)^{3}}\right)^{4}$
$\displaystyle\leq\sum_{\tau\in\mathfrak{T}_{T}}h_{\tau}^{-3}\left(\|{\boldsymbol{\pi}}_{T}^{k}\boldsymbol{w}-\boldsymbol{w}\|_{L^{2}(T)^{3}}+\|\boldsymbol{w}-\boldsymbol{R}_{T}^{k}\uline{\hat{\boldsymbol{w}}}_{T}\|_{L^{2}(T)^{3}}\right)^{4}$
$\displaystyle\lesssim\sum_{\tau\in\mathfrak{T}_{T}}h_{\tau}^{-3}h_{T}^{4(k+1)}|\boldsymbol{w}|_{H^{k+1}(T)^{3}}^{4}\lesssim
h_{T}^{4(k+1)}|\boldsymbol{w}|_{W^{k+1,4}(T)^{3}}^{4},$

where: to pass to the second line we have used the reverse Lebesgue embedding
(43) for $(X,\alpha,\beta)=(T,4,2)$; to pass to the third line, we have
inserted $\pm\boldsymbol{w}$ and used a triangle inequality along with
$\tau\subset T$; to pass to the fourth line, we have used the approximation
properties (5a) of the $L^{2}$-orthogonal projector with
$(l,m,r,s)=(k,0,2,k+1)$ for the first addend and the approximation property
(16) for the second addend; the conclusion follows from
$h_{\tau}^{-1}h_{T}\lesssim 1$ (consequence of mesh regularity), the bound (3)
on $\operatorname{card}(\mathfrak{T}_{T})$, and the Lebesgue embedding
$\|\zeta\|_{L^{2}(T)}\lesssim h_{T}^{\frac{3}{4}}\|\zeta\|_{L^{4}(T)}$ valid
for all $\zeta\in L^{4}(T)$. Plugging the above bound into (73), we get

$\|\boldsymbol{w}-\boldsymbol{R}_{T}^{k}\uline{\hat{\boldsymbol{w}}}_{T}\|_{L^{4}(T)^{3}}\lesssim
h_{T}^{k+1}|\boldsymbol{w}|_{W^{k+1,4}(T)^{3}}.$

In conclusion, we have that

$\displaystyle|\mathfrak{T}_{3}|$ $\displaystyle\lesssim
h^{k+1}|\boldsymbol{w}|_{H^{1}(\Omega)^{3}}|\boldsymbol{w}|_{W^{k+1,4}({\mathcal{T}}_{h})^{3}}\|\uline{\boldsymbol{z}}_{h}\|_{1,h}.$
(74)

Using similar arguments as for $\mathfrak{T}_{3}$, we have for the last term

$|\mathfrak{T}_{4}|\lesssim
h^{k+1}|\boldsymbol{w}|_{H^{1}(\Omega)^{3}}|\boldsymbol{w}|_{W^{k+1,4}({\mathcal{T}}_{h})^{3}}\|\uline{\boldsymbol{z}}_{h}\|_{1,h}.$
(75)

(iv.D) _Conclusion._ Taking absolute values in (61), recalling the definition
(9) of the dual norm, and using the estimates (62), (72), (74), and (75), and
additionally noticing that
$|\boldsymbol{w}|_{H^{1}(\Omega)^{3}}\lesssim|\boldsymbol{w}|_{W^{1,4}(\Omega)^{3}}$,
the conclusion follows. ∎

###  3.4 Discrete problem and main results

The HHO discretization of problem (1) reads: Find
$(\uline{\boldsymbol{u}}_{h},p_{h})\in\uline{\boldsymbol{U}}_{h,0}^{k}\times
P_{h}^{k}$ such that

$\displaystyle\nu
a_{h}(\uline{\boldsymbol{u}}_{h},\uline{\boldsymbol{v}}_{h})+t_{h}(\uline{\boldsymbol{u}}_{h},\uline{\boldsymbol{u}}_{h},\uline{\boldsymbol{v}}_{h})+b_{h}(\uline{\boldsymbol{v}}_{h},p_{h})$
$\displaystyle=\ell_{h}(\boldsymbol{f},\uline{\boldsymbol{v}}_{h})$
$\displaystyle\forall\uline{\boldsymbol{v}}_{h}\in\uline{\boldsymbol{U}}^{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}k}_{h,0},\marginpar{\normalcolor\raggedright\scriptsize\textsf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}
R\\#1 minor point: Typo}}\@add@raggedright$ (76a) $\displaystyle-
b_{h}(\uline{\boldsymbol{u}}_{h},q_{h})$ $\displaystyle=0$
$\displaystyle\forall q_{h}\in\mathbb{P}^{k}({\mathcal{T}}_{h}).$ (76b)

The existence of a solution to (76) for any $\boldsymbol{f}\in
L^{2}(\Omega)^{3}$ can be proved using a topological degree argument as in
[23, Theorem 1]. Similarly, uniqueness can be proved along the lines of
Theorem 2 therein under a smallness condition on $\boldsymbol{f}$.

Recalling the Hodge decomposition (2) and denoting by $C_{P}$ a Poincaré
constant in $\Omega$, Proposition 11 below is the discrete equivalent of the
following a priori continuous bound (see [11, Section 2.3])

$|\boldsymbol{u}|_{H^{1}(\Omega)^{3}}\leq\nu^{-1}C_{\rm
P}\|\boldsymbol{g}\|_{L^{2}(\Omega)^{3}}.$ (77)

######  Proposition 11 (Uniform a priori bound on the discrete velocity).

Let
$(\uline{\boldsymbol{u}}_{h},p_{h})\in\uline{\boldsymbol{U}}_{h,0}^{k}\times
P_{h}^{k}$ be a solution to (76). Then, given the Hodge decomposition (2) of
$\boldsymbol{f}$, we have the following uniform a priori bound for the
velocity:

$\|\uline{\boldsymbol{u}}_{h}\|_{1,h}\lesssim\nu^{-1}\|\boldsymbol{g}\|_{L^{2}(\Omega)^{3}}.$

###### Proof.

The proof follows the same reasoning as [11, Proposition 8] with Lemmas 8 and
10 replacing, respectively, [11, Eqs. (41)–(42) and Lemma 7]. ∎

######  Remark 12 (Efficient implementation).

When solving the algebraic problem corresponding to (76) by a first order
iterative algorithm, all element-based velocity unknowns and all but one
pressure unknowns per element can be statically condensed at each iteration in
the spirit of [20, Section 6.2] ; see [7] for a study of the effect of static
condensation strategies on the multigrid resolution of the global algebraic
systems arising from HHO discretizations of incompressible flow problems.

######  Remark 13 (The two-dimensional case).

The two-dimensional version of the method (76) will be considered numerically
in Section 4. Denoting by $u_{i}$, $i=1,\ldots,3$, the component of the
velocity field along the Cartesian axis $x_{i}$, the two-dimensional plane
velocity problem can be recovered from (1) setting $u_{3}=0$ and assuming that
$u_{1}$ and $u_{2}$ do not depend on $x_{3}$.

We next consider the discretization error defined as the difference between
the solution to the HHO scheme and the interpolate of the exact solution.

######  Theorem 14 (Error estimate for small data).

Recalling the Hodge decomposition (2) of the forcing term $\boldsymbol{f}$, we
assume that it holds, for some $\alpha\in(0,1),$

$\|\boldsymbol{g}\|_{{L}^{2}({\Omega})^{3}}\leq\alpha\frac{\nu^{2}C_{a}}{C_{t}C_{I}C_{P}},$

where $C_{a}$ and $C_{t}$ are defined in (53) and (59), while $C_{I}$ denotes
the continuity constant of the HHO interpolator in the discrete $H^{1}$-like
norm (see [16, Proposition 2.2]) and $C_{P}$ is the Poincaré constant in (77).
††margin:  Referees #1 and #2 remark: The error estimate is valid now for
arbitrary $k\geq 0$ Let $k\geq 0$ and let
$(\boldsymbol{u},p)\in\boldsymbol{U}\times P$ and
$(\underline{\boldsymbol{u}}_{h},p_{h})\in\uline{\boldsymbol{U}}_{h}^{k}\times
P_{h}^{k}$ solve (1) and (76), respectively. Assuming the additional
regularity $\boldsymbol{u}\in H^{k+2}({\mathcal{T}}_{h})^{3}$ and $p\in
H^{1}(\Omega)\cap H^{k+1}({\mathcal{T}}_{h})$, it holds:

$\|\uline{\boldsymbol{u}}_{h}-\uline{\boldsymbol{I}}_{h}^{k}\boldsymbol{u}\|_{1,h}+\nu^{-1}\|p_{h}-\pi_{h}^{k}p\|_{{L}^{2}({\Omega})}\\\
\lesssim
h^{k+1}(1-\alpha)^{-1}\left(|\boldsymbol{u}|_{H^{k+2}({\mathcal{T}}_{h})^{3}}+{\nu}^{-1}\|\boldsymbol{u}\|_{W^{1,4}(\Omega)^{3}}|\boldsymbol{u}|_{W^{k+1,4}({\mathcal{T}}_{h})^{3}}\right).$
(78)

where the hidden constant is independent of $\nu$, $\lambda$, $h$, as well as
$(\boldsymbol{u},p)$.

###### Proof.

Analogous to that of [11, Theorem 11]. ∎

######  Remark 15 (Pressure robustness).

The error estimate (78) is pressure-robust since the right-hand side does not
depend on $\lambda$ in (2) nor on the pressure.

##  4 Numerical tests

In this section we verify numerically the proposed method for general meshes
with convex elements for $\Omega\subset\mathbb{R}^{2}$. For each element
$T\in\mathcal{T}_{h}$, we construct its simplicial submesh $\mathfrak{T}_{T}$
††margin:  R#2 point: Explaining how $\mathfrak{T}_{T}$ is constructed using
an ear clipping algorithm, i.e., we construct $\mathfrak{T}_{T}$ in such a way
that no additional internal nodes are introduced and that
$\mathfrak{F}_{T}=\mathcal{F}_{T}$ (this construction fulfills the assumptions
made in Section 2.1). For the sake of completeness, we also include
comparisons with the original HHO method of [8]. Our implementation is based
on the HArDCore library333https://github.com/jdroniou/HArDCore and makes
extensive use of the linear algebra Eigen open-source library [32]. All the
steady-state computations presented hereafter are done by means of the pseudo-
transient-continuation algorithm analyzed by [34] employing the Selective
Evolution Relaxation (SER) strategy [44] for evolving the pseudo-time step
according to the Newton’s equations residual. Convergence to steady-state is
attained when the Euclidean norm of the residual for the momentum equation
drops below $10^{-11}$. At each pseudo-time step, the linearized equations are
exactly solved by means of the direct solver Pardiso [48]. Accordingly, the
Euclidean norm of the residual for the continuity equation is comparable to
the machine epsilon at all pseudo-time steps.

(a) Cartesian.

(b) Hexagonal.

(c) Kershaw.

Figure 3: Coarsest meshes used in Section 4.1.

###  4.1 Kovasznay flow

We start by assessing the convergence properties of the method using the well
known analytical solution of Kovasznay [36] with $\nu=0.025$; see, e.g., [16,
Section 6.1] for the expression of the velocity and pressure fields. We
consider computations over three $h$-refined mesh families (Cartesian,
hexagonal and Kershaw type). Figure 3 shows the coarsest mesh for each family.
We monitor the following quantities in Table 1: $N_{\rm dof}$ and $N_{\rm nz}$
denoting, respectively, the number of discrete unknowns and nonzero entries of
the statically condensed linearized problem;
$\|\uline{\boldsymbol{e}}_{h}\|_{\nu,h}\coloneq\left[\nu\mathrm{a}_{h}(\uline{\boldsymbol{e}}_{h},\uline{\boldsymbol{e}}_{h})\right]^{\nicefrac{{1}}{{2}}}$,
the energy norm of the error
$\uline{\boldsymbol{e}}_{h}\coloneq\uline{\boldsymbol{u}}_{h}-\uline{\boldsymbol{I}}_{h}^{k}\boldsymbol{u}$
on the velocity (using the global norm equivalence (53), an estimate in
$h^{k+1}$ for this quantity is readily inferred from (78));
$\|\boldsymbol{e}_{h}\|_{L^{2}(\Omega)^{2}}$ and
$\|\epsilon_{h}\|_{L^{2}(\Omega)}$, the $L^{2}$-errors on the velocity and the
pressure, respectively. Each error measure is accompanied by the corresponding
Estimated Order of Convergence (EOC) computed using successive refinement
steps. The results collected in Table 1 show that both the energy norm of the
error on the velocity and the $L^{2}$-norm of error on the pressure converge
as $h^{k+1}$ as expected. Additionally, the $L^{2}$-norm of the error of the
velocity converges with rates close to $h^{k+2}$.

$N_{\rm dof}$ | $\|\uline{\boldsymbol{e}}_{h}\|_{\nu,h}$ | EOC | $\|{\boldsymbol{e}}_{h}\|_{L^{2}(\Omega)^{2}}$ | EOC | $\|\epsilon_{h}\|_{L^{2}(\Omega)}$ | EOC
---|---|---|---|---|---|---
Cartesian, $k=0$
540 | 5.92E-01 | – | 1.11E-01 | – | 2.54E-01 | –
2080 | 3.35E-01 | 0.822 | 3.54E-02 | 1.644 | 8.58E-02 | 1.565
8160 | 1.78E-01 | 0.911 | 9.92E-03 | 1.835 | 2.46E-02 | 1.800
32320 | 9.10E-02 | 0.969 | 2.59E-03 | 1.938 | 6.50E-03 | 1.923
Cartesian, $k=1$
980 | 2.04E-01 | – | 1.97E-02 | – | 4.91E-02 | –
3760 | 5.73E-02 | 1.831 | 2.28E-03 | 3.114 | 5.86E-03 | 3.067
14720 | 1.51E-02 | 1.926 | 2.91E-04 | 2.970 | 7.75E-04 | 2.918
58240 | 3.85E-03 | 1.96 | 3.75E-05 | 2.957 | 1.08E-04 | 2.845
Hexagonal, $k=0$
3241 | 8.26E-01 | – | 5.46E-02 | – | 1.53E-01 | –
12081 | 4.42E-01 | 0.901 | 1.71E-02 | 1.674 | 4.21E-02 | 1.861
46561 | 2.27E-01 | 0.964 | 4.78E-03 | 1.838 | 1.10E-02 | 1.929
182721 | 1.14E-01 | 0.986 | 1.25E-03 | 1.931 | 2.85E-03 | 1.954
Hexagonal, $k=1$
6041 | 4.80E-01 | – | 3.04E-02 | – | 1.03E-01 | –
22481 | 8.40E-02 | 2.515 | 1.42E-03 | 4.421 | 4.53E-03 | 4.507
86561 | 1.78E-02 | 2.237 | 1.23E-04 | 3.529 | 4.20E-04 | 3.432
339521 | 4.20E-03 | 2.085 | 1.37E-05 | 3.164 | 6.84E-05 | 2.618
Kershaw, $k=0$
5577 | 5.76E-01 | – | 1.78E-01 | – | 2.10E-01 | –
22044 | 2.46E-01 | 1.231 | 6.36E-02 | 1.488 | 8.74E-02 | 1.269
49401 | 1.50E-01 | 1.230 | 2.93E-02 | 1.921 | 4.27E-02 | 1.774
87648 | 1.08E-01 | 1.142 | 1.66E-02 | 1.983 | 2.47E-02 | 1.908
136785 | 8.46E-02 | 1.092 | 1.06E-02 | 1.995 | 1.60E-02 | 1.951
Kershaw, $k=1$
10065 | 4.51E-01 | – | 5.88E-02 | – | 1.50E-01 | –
39732 | 7.00E-02 | 2.701 | 2.35E-03 | 4.672 | 5.12E-03 | 4.897
89001 | 3.02E-02 | 2.080 | 5.08E-04 | 3.785 | 1.09E-03 | 3.830
157872 | 1.64E-02 | 2.131 | 1.78E-04 | 3.653 | 4.23E-04 | 3.299
246345 | 1.04E-02 | 2.032 | 7.80E-05 | 3.703 | 2.11E-04 | 3.112
Table 1: Convergence rates for the numerical test of Section 4.1.

###  4.2 Robustness of the velocity error estimate

The second numerical example, inspired by [41, Benchmark 3.3], is meant to
demonstrate the robustness of the proposed method for large irrotational body
forces. Specifically, we verify numerically the fact that the approximation of
the velocity is independent of both $\lambda$ and $p$. Letting
$\Omega=(0,1)^{2}$ and $\lambda\geq 0$, we solve the Dirichlet problem
corresponding to the exact solution $(\boldsymbol{u},p)$ in (1) with velocity
components given by
$\boldsymbol{u}(\boldsymbol{x})\coloneq\begin{pmatrix}-x_{2}\\\
x_{1}\end{pmatrix}$ and pressure given by
$p(\boldsymbol{x})\coloneq{\lambda}x_{1}^{3}+\frac{x_{1}^{2}+x_{2}^{2}}{2}-\frac{1}{4}$.
We set $\nu=1$, then observe that the force in (1a) is purely irrotational,
i.e., $\boldsymbol{f}(\boldsymbol{x})=\begin{pmatrix}3\lambda x_{1}^{2}\\\
0\end{pmatrix}$. In the computations, we take $\lambda=10^{6}$ and consider a
sequence of uniformly $h$-refined meshes equivalent (by scaling and
translation) to the three mesh families used in the previous section, see
Figure 3. Table 2 collects the results for the Cartesian and hexagonal mesh
families, and Table 3 for the Kershaw mesh family. For the sake of comparison,
we also report in these tables the corresponding results obtained using the
original HHO method of [8]. It can be noticed that the solution is exactly
reproduced by the present method with $k=1$ on all the meshes, while a quick
convergence is observed for $k=0$ on the hexagonal and Kershaw meshes, most
likely due to the quadratic nature of the pressure. By contrast, the HHO
method of [8] shows large errors on the velocity due to the lack of pressure-
robustness.

$N_{\rm dof}$ | $\|\uline{\boldsymbol{e}}_{h}\|_{\nu,h}$ | EOC | $\|{\boldsymbol{e}}_{h}\|_{L^{2}(\Omega)^{2}}$ | EOC | $\|\epsilon_{h}\|_{L^{2}(\Omega)}$ | EOC | $\|\uline{\boldsymbol{e}}_{h}\|_{\nu,h}$ | EOC | $\|{\boldsymbol{e}}_{h}\|_{L^{2}(\Omega)^{2}}$ | EOC | $\|\epsilon_{h}\|_{L^{2}(\Omega)}$ | EOC
---|---|---|---|---|---|---|---|---|---|---|---|---
|  Cartesian, $k=0$, proposed method |  Cartesian, $k=0$, HHO method of [8]
540 | 1.94E-11 | – | 1.74E-12 | – | 8.28E-16 | – | 2.77E+05 | – | 1.79E+04 | – | 1.22E+02 | –
2080 | 3.06E-11 | – | 1.68E-12 | 0.052 | 1.09E-15 | – | 2.34E+04 | 3.562 | 9.93E+02 | 4.170 | 1.56E+00 | 6.290
8160 | 3.61E-11 | – | 6.22E-12 | – | 1.88E-15 | – | 1.18E+04 | 0.991 | 2.57E+02 | 1.951 | 1.13E-01 | 3.789
32320 | 3.30E-11 | 0.129 | 1.77E-12 | 1.816 | 1.01E-15 | 0.904 | 5.92E+03 | 0.994 | 6.54E+01 | 1.973 | 8.67E-03 | 3.705
|  Cartesian, $k=1$, proposed method |  Cartesian, $k=1$, HHO method of [8]
980 | 1.44E-10 | – | 5.16E-12 | – | 6.07E-05 | – | 2.01E+03 | – | 9.91E+01 | – | 3.09E-02 | –
3760 | 1.41E-10 | 0.034 | 3.22E-12 | 0.682 | 7.59E-06 | 3.000 | 5.06E+02 | 1.987 | 1.26E+01 | 2.977 | 5.06E-04 | 5.936
14720 | 1.81E-10 | – | 2.89E-12 | 0.155 | 9.49E-07 | 3.000 | 1.27E+02 | 1.994 | 1.58E+00 | 2.990 | 8.66E-06 | 5.867
58240 | 1.47E-10 | 0.306 | 2.42E-12 | 0.257 | 1.19E-07 | 3.000 | 3.18E+01 | 1.997 | 1.99E-01 | 2.995 | 5.06E-07 | 4.097
|  Hexagonal, $k=0$, proposed method |  Hexagonal, $k=0$, HHO method of [8]
3241 | 1.06E-01 | – | 4.18E-03 | – | 5.17E-05 | – | 1.97E+04 | – | 7.80E+02 | – | 1.04E+00 | –
12081 | 9.78E-03 | 3.434 | 2.51E-04 | 4.056 | 9.34E-06 | 2.468 | 1.48E+04 | 0.415 | 4.37E+02 | 0.838 | 3.33E-01 | 1.640
46561 | 8.85E-04 | 3.466 | 1.53E-05 | 4.038 | 1.67E-06 | 2.484 | 8.18E+03 | 0.854 | 1.36E+02 | 1.686 | 3.52E-02 | 3.244
182721 | 7.92E-05 | 3.482 | 9.42E-07 | 4.021 | 2.97E-07 | 2.492 | 4.15E+03 | 0.979 | 3.51E+01 | 1.950 | 2.92E-03 | 3.591
|  Hexagonal, $k=1$, proposed method |  Hexagonal, $k=1$, HHO method of [8]
6041 | 2.32E-10 | – | 2.13E-11 | – | 7.13E-06 | – | 6.72E+02 | – | 1.81E+01 | – | 1.01E-03 | –
22481 | 2.24E-10 | 0.050 | 1.08E-11 | 0.982 | 9.17E-07 | 2.959 | 1.76E+02 | 1.933 | 2.23E+00 | 3.019 | 2.74E-05 | 5.207
86561 | 2.60E-10 | – | 1.08E-11 | 0.007 | 1.16E-07 | 2.980 | 4.46E+01 | 1.981 | 2.81E-01 | 2.989 | 4.00E-06 | 2.774
339521 | 6.67E-10 | – | 1.36E-10 | – | 1.46E-08 | 3.0 | 1.12E+01 | 1.990 | 3.53E-02 | 2.991 | 6.96E-07 | 2.525

Table 2: Convergence rates for the numerical test of Section 4.2 for
$\lambda=10^{6}$ using the Cartesian and hexagonal mesh families.

.

$N_{\rm dof}$ | $\|\uline{\boldsymbol{e}}_{h}\|_{\nu,h}$ | EOC | $\|{\boldsymbol{e}}_{h}\|_{L^{2}(\Omega)^{2}}$ | EOC | $\|\epsilon_{h}\|_{L^{2}(\Omega)}$ | EOC | $\|\uline{\boldsymbol{e}}_{h}\|_{\nu,h}$ | EOC | $\|{\boldsymbol{e}}_{h}\|_{L^{2}(\Omega)^{2}}$ | EOC | $\|\epsilon_{h}\|_{L^{2}(\Omega)}$ | EOC
---|---|---|---|---|---|---|---|---|---|---|---|---
|  Kershaw, $k=0$, proposed method |  Kershaw, $k=0$, HHO method of [8]
5577 | 4.55E-03 | – | 7.46E-04 | – | 2.72E-07 | – | 1.42E+04 | – | 3.74E+02 | – | 2.43E-01 | –
22044 | 2.83E-04 | 4.027 | 4.73E-05 | 3.999 | 1.76E-08 | 3.970 | 7.17E+03 | 0.991 | 9.60E+01 | 1.970 | 1.80E-02 | 3.771
49401 | 5.59E-05 | 4.015 | 9.37E-06 | 4.006 | 3.50E-09 | 3.995 | 4.79E+03 | 1.000 | 4.30E+01 | 1.986 | 4.32E-03 | 3.532
87648 | 1.77E-05 | 4.010 | 2.97E-06 | 4.005 | 1.11E-09 | 3.992 | 3.59E+03 | 1.001 | 2.43E+01 | 1.985 | 1.73E-03 | 3.187
136785 | 7.23E-06 | 4.008 | 1.22E-06 | 4.005 | 4.59E-10 | 3.967 | 2.87E+03 | 1.001 | 1.56E+01 | 1.986 | 9.19E-04 | 2.845
|  Kershaw, $k=1$, proposed method |  Kershaw, $k=1$, HHO method of [8]
10065 | 5.35E-10 | – | 8.69E-11 | – | 1.69E-06 | – | 1.86E+02 | – | 2.88E+00 | – | 3.56E-05 | –
39732 | 5.45E-10 | – | 7.84E-11 | 0.150 | 2.11E-07 | 3.014 | 4.67E+01 | 2.006 | 3.56E-01 | 3.033 | 2.31E-06 | 3.965
89001 | 1.09E-09 | – | 1.76E-10 | – | 6.26E-08 | 3.008 | 2.08E+01 | 2.004 | 1.05E-01 | 3.013 | 5.99E-07 | 3.343
157872 | 1.60E-09 | – | 2.72E-10 | – | 2.64E-08 | 3.006 | 1.17E+01 | 2.003 | 4.43E-02 | 3.008 | 2.36E-07 | 3.247
246345 | 4.53E-10 | 5.665 | 4.84E-11 | 7.756 | 1.35E-08 | 3.005 | 7.49E+00 | 2.002 | 2.27E-02 | 3.005 | 1.16E-07 | 3.189
Table 3: Convergence rates for the numerical test of Section 4.2 for
$\lambda=10^{6}$ using the Kershaw mesh family.

###  4.3 Two-dimensional lid-driven cavity flow

The final numerical test is the classical two-dimensional lid-driven cavity
problem. The computational domain is the unit square $\Omega=(0,1)^{2}$ and we
initially set $\boldsymbol{f}=\boldsymbol{0}$. Homogeneous (wall) boundary
conditions are enforced at all but the top horizontal wall (at $x_{2}=1$),
where we enforce a unit tangential velocity $\boldsymbol{u}=(1,0)$ instead. In
Figure 4 we report the horizontal component $u_{1}$ of the velocity along the
vertical centerline $x_{1}=\frac{1}{2}$ and the vertical component $u_{2}$ of
the velocity along the horizontal centerline $x_{2}=\frac{1}{2}$ for a global
Reynolds number $\mathrm{Re}\coloneq\frac{1}{\nu}=1000$. The computation is
carried out setting $k=1$ for the finest meshes of the Cartesian, hexagonal,
and Kershaw sequences used in the previous section. Reference solutions from
the literature [30, 25] are also included for the sake of comparison. The
numerical solution obtained using the proposed method is in agreement with the
reference results for the selected value of the Reynolds number.

To check the robustness of the method with respect to irrotational body
forces, we then run the same test case but with
$\boldsymbol{f}=\lambda\nabla\psi$ where $\psi=\frac{1}{3}(x^{3}+y^{3})$. This
body force is completely irrotational, so the velocity approximation obtained
using the proposed method (76) should not be affected (and, therefore, should
not depend on $\lambda$). To verify this, we report in Figure 5 computations
for $\lambda=10^{6}$, using $k=1$ and the same meshes as before. As expected,
the velocity profiles are not affected by the value of $\lambda$. The same
plot also contains the results obtained with the original HHO formulation of
[8], but only for the Cartesian mesh and $\lambda=10^{3}$ (convergence was not
achieved for $\lambda=10^{6}$). It can be checked that the non-pressure-robust
version of the method converges to a complete different solution.

$-1$$-0.8$$-0.6$$-0.4$$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$u_{1}$$x_{2}$$k=1$,
cartesian mesh$k=1$, hexagonal mesh$k=1$, Kershaw meshGhia et al.Erturk et
al.$0$$0.2$$0.4$$0.6$$0.8$$1$$-1$$-0.8$$-0.6$$-0.4$$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$x_{1}$$u_{2}$$-1$$-0.8$$-0.6$$-0.4$$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$u_{1}$$x_{2}$$k=1$,
cartesian mesh$k=1$, hexagonal mesh$k=1$, Kershaw meshGhia et al.Erturk et
al.$0$$0.2$$0.4$$0.6$$0.8$$1$$-1$$-0.8$$-0.6$$-0.4$$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$x_{1}$$u_{2}$$-1$$-0.8$$-0.6$$-0.4$$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$u_{1}$$x_{2}$$k=1$,
cartesian mesh$k=1$, hexagonal mesh$k=1$, Kershaw meshGhia et al.Erturk et
al.$0$$0.2$$0.4$$0.6$$0.8$$1$$-1$$-0.8$$-0.6$$-0.4$$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$x_{1}$$u_{2}$$-1$$-0.8$$-0.6$$-0.4$$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$u_{1}$$x_{2}$$k=1$,
cartesian mesh$k=1$, hexagonal mesh$k=1$, Kershaw meshGhia et al.Erturk et
al.$0$$0.2$$0.4$$0.6$$0.8$$1$$-1$$-0.8$$-0.6$$-0.4$$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$x_{1}$$u_{2}$$-1$$-0.8$$-0.6$$-0.4$$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$u_{1}$$x_{2}$$k=1$,
cartesian mesh$k=1$, hexagonal mesh$k=1$, Kershaw meshGhia et al.Erturk et
al.$0$$0.2$$0.4$$0.6$$0.8$$1$$-1$$-0.8$$-0.6$$-0.4$$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$x_{1}$$u_{2}$
Figure 4: Two-dimensional lid-driven cavity flow, horizontal component $u_{1}$
of the velocity along the vertical centerline $x_{1}=\frac{1}{2}$ and the
vertical component $u_{2}$ of the velocity along the horizontal centerline
$x_{2}=\frac{1}{2}$ for $\mathrm{Re}=1{,}000$.
$-1$$-0.8$$-0.6$$-0.4$$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$u_{1}$$x_{2}$$\lambda=10^{6}$,
Cartesian, present m.$\lambda=10^{6}$, hexagonal, present m.$\lambda=10^{6}$,
Kershaw, present m.$\lambda=10^{3}$, Cartesian, [8]Ghia et al.Erturk et
al.$0$$0.2$$0.4$$0.6$$0.8$$1$$-1$$-0.8$$-0.6$$-0.4$$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$x_{1}$$u_{2}$
Figure 5: Two-dimensional lid-driven cavity flow with irrotational force
$\boldsymbol{f}=\lambda\nabla\psi$ with $\lambda=10^{6}$. Comparison between
the present method and the original HHO formulation of [8] both using $k=1$.
The plot represents the horizontal component $u_{1}$ of the velocity along the
vertical centerline $x_{1}=\frac{1}{2}$ and the vertical component $u_{2}$ of
the velocity along the horizontal centerline $x_{2}=\frac{1}{2}$ for
$\mathrm{Re}=1{,}000$.

## Acknowledgements

††margin:  Acknowledging Referee #2 remarks

The authors are grateful for the suggestions of an anonymous referee which
contributed to improving the quality of the manuscript. Daniele Di Pietro
acknowledges the partial support of _Agence Nationale de la Recherche_ grant
ANR-20-MRS2-0004 “NEMESIS” and I-Site MUSE grant ANR-16-IDEX-0006 “RHAMNUS”.

## References

  * [1] Naveed Ahmed, Alexander Linke and Christian Merdon “Towards pressure-robust mixed methods for the incompressible Navier-Stokes equations” In _Comput. Methods Appl. Math._ 18.3, 2018, pp. 353–372 DOI: 10.1515/cmam-2017-0047
  * [2] D. Arnold “Finite Element Exterior Calculus” SIAM, 2018
  * [3] M. Bebendorf “A Note on the Poincaré Inequality for Convex Domains” In _Z. Anal. Anwend._ 22.4, 2003, pp. 751–756 DOI: 10.4171/ZAA/1170
  * [4] L. Beirão da Veiga, F. Dassi and G. Vacca “The Stokes complex for Virtual Elements in three dimensions” In _Math. Models Methods Appl. Sci._ 30.03, 2020, pp. 477–512 DOI: 10.1142/S0218202520500128
  * [5] L. Beirão da Veiga, C. Lovadina and G. Vacca “Virtual Elements for the Navier–Stokes Problem on Polygonal Meshes” In _SIAM J. Numer. Anal._ 56.3, 2018, pp. 1210–1242 DOI: 10.1137/17M1132811
  * [6] D. Boffi, F. Brezzi and M. Fortin “Mixed finite element methods and applications” 44, Springer Series in Computational Mathematics Heidelberg: Springer, 2013, pp. xiv+685 DOI: 10.1007/978-3-642-36519-5
  * [7] L. Botti and D. A. Di Pietro “$p$-Multilevel preconditioners for HHO discretizations of the Stokes equations with static condensation” In _Commun. Appl. Math. Comput._ 4.3, 2022, pp. 783–822 DOI: 10.1007/s42967-021-00142-5
  * [8] L. Botti, Daniele A. Di Pietro and J. Droniou “A Hybrid High-Order method for the incompressible Navier-Stokes equations based on Temam’s device” In _J. Comput. Phys._ 376, 2019, pp. 786–816 DOI: 10.1016/j.jcp.2018.10.014
  * [9] L. Botti and F. C. Massa “HHO methods for the incompressible Navier-Stokes and the incompressible Euler equations”, 2021 arXiv:2112.09777 [math.NA]
  * [10] M. Botti, D. Castanon Quiroz, D. A. Di Pietro and A. Harnist “A hybrid high-order method for creeping flows of non-Newtonian fluids” In _ESAIM: Math. Model. Numer. Anal._ 55.5, 2021, pp. 2045–2073 DOI: 10.1051/m2an/2021051
  * [11] D. Castanon Quiroz and D. A. Di Pietro “A Hybrid High-Order method for the incompressible Navier–Stokes problem robust for large irrotational body forces” In _Comput. Math. Appl._ 79.9, 2020 DOI: 10.1016/j.camwa.2019.12.005
  * [12] Daniel Castanon Quiroz, Daniele A Di Pietro and André Harnist “A Hybrid High-Order method for incompressible flows of non-Newtonian fluids with power-like convective behaviour” In _IMA J. Numer. Anal._ , 2021 DOI: 10.1093/imanum/drab087
  * [13] W. Chen and Y. Wang “Minimal degree H(curl) and H(div) conforming finite elements on polytopal meshes” In _Math. Comput._ 86.307, 2017, pp. 2053–2087 DOI: 10.1090/mcom/3152
  * [14] P.G. Ciarlet and J.L. Lions “Handbook of Numerical Analysis: VOL II: Finite Element Methods. (Part 1).” North-Holland, 1991
  * [15] D. A. Di Pietro and J. Droniou “An arbitrary-order discrete de Rham complex on polyhedral meshes: Exactness, Poincaré inequalities, and consistency” Published online (open access) In _Found. Comput. Math._ , 2021 DOI: 10.1007/s10208-021-09542-8
  * [16] D. A. Di Pietro and J. Droniou “The Hybrid High-Order Method for Polytopal Meshes - Design, Analysis and Applications” 19, Springer Series in Modeling, Simulation and Applications, 2020 DOI: 10.1007/978-3-030-37203-3
  * [17] D. A. Di Pietro, J. Droniou and F. Rapetti “Fully discrete polynomial de Rham sequences of arbitrary degree on polygons and polyhedra” In _Math. Models Methods Appl. Sci._ 30.9, 2020, pp. 1809–1855 DOI: 10.1142/S0218202520500372
  * [18] D. A. Di Pietro and A. Ern “Discrete functional analysis tools for discontinuous Galerkin methods with application to the incompressible Navier–Stokes equations” In _Math. Comp._ 79, 2010, pp. 1303–1330 DOI: 10.1090/S0025-5718-10-02333-1
  * [19] D. A. Di Pietro and A. Ern “Mathematical aspects of discontinuous Galerkin methods” 69, Mathématiques & Applications Springer, Heidelberg, 2012 DOI: 10.1007/978-3-642-22980-0
  * [20] D. A. Di Pietro, A. Ern, A. Linke and F. Schieweck “A discontinuous skeletal method for the viscosity-dependent Stokes problem” In _Comput. Meth. Appl. Mech. Engrg._ 306, 2016, pp. 175–195 DOI: 10.1016/j.cma.2016.03.033
  * [21] Daniele A. Di Pietro and J. Droniou “A Hybrid High-Order method for Leray-Lions elliptic equations on general meshes” In _Math. Comp._ 86, 2017, pp. 2159–2191 DOI: 10.1090/mcom/3180
  * [22] Daniele A. Di Pietro and A. Ern “A hybrid high-order locking-free method for linear elasticity on general meshes” In _Meth. Appl. Mech. Engrg._ 283, 2015, pp. 1–21
  * [23] Daniele A. Di Pietro and Stella Krell “A Hybrid High-Order Method for the Steady Incompressible Navier–Stokes Problem” In _J. Sci. Comput._ 74.3, 2018, pp. 1677–1705 DOI: 10.1007/s10915-017-0512-x
  * [24] A. Ern and J.-L. Guermond “Finite Elements I, Approximation and Interpolation” Texts in Applied Mathematics 72, Springer-Verlag, New York, 2021
  * [25] E. Erturk, T. C. Corke and C. Gökçöl “Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds” In _Int. J. Numer. Meth. Fluids_ 48.7, 2005, pp. 747–774 DOI: 10.1002/fld.953
  * [26] Richard S. Falk and Michael Neilan “Stokes complexes and the construction of stable finite elements with pointwise mass conservation” In _SIAM J. Numer. Anal._ 51.2, 2013, pp. 1308–1326 DOI: 10.1137/120888132
  * [27] D. Frerichs and C. Merdon “Divergence-preserving reconstructions on polygons and a really pressure-robust virtual element method for the Stokes problem” In _IMA J. Numer. Anal._ , 2020 DOI: 10.1093/imanum/draa073
  * [28] G. P. Galdi “An Introduction to the Mathematical Theory of the Navier–Stokes Equations. Steady–State Problems. ” Second Ed., Springer Monographs in Mathematics. Springer, 2011
  * [29] Gabriel N. Gatica, Mauricio Munar and Filánder A. Sequeira “A mixed virtual element method for the Navier–Stokes equations” In _Math. Models Methods Appl. Sci._ 28.14, 2018, pp. 2719–2762 DOI: 10.1142/S0218202518500598
  * [30] U. Ghia, K.N. Ghia and C.T. Shin “High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method” In _J. Comput. Phys._ 48.3, 1982, pp. 387–411 DOI: 10.1016/0021-9991(82)90058-4
  * [31] P. Grisvard “Elliptic Problems in Nonsmooth Domains” Society for IndustrialApplied Mathematics, 2011 DOI: 10.1137/1.9781611972030
  * [32] B. Guennebaud and al. “Eigen v3”, 2010 URL: http://eigen.tuxfamily.org
  * [33] V. John, A. Linke, C. Merdon, M. Neilan and L. G. Rebholz “On the divergence constraint in mixed finite element methods for incompressible flows” In _SIAM Rev._ 59.3, 2017, pp. 492–544
  * [34] C. Kelley and D. Keyes “Convergence analysis of pseudo-transient continuation” In _SIAM J. Numer. Anal._ 35.2, 1998, pp. 508–523 DOI: 10.1137/S0036142996304796
  * [35] D. Kim, L. Zhao, E. Chung and E.-J. Park “Pressure-robust staggered DG methods for the Navier-Stokes equations on general meshes”, 2021 arXiv:2107.09226 [math.NA]
  * [36] L. I. G. Kovasznay “Laminar flow behind a two-dimensional grid” In _Proceedings of the Cambridge Philosophical Society_ 44.1, 1948, pp. 58–62 DOI: 10.1017/S0305004100023999
  * [37] Christian Kreuzer, Rüdiger Verfürth and Pietro Zanotti “Quasi-Optimal and Pressure Robust Discretizations of the Stokes Equations by Moment- and Divergence-Preserving Operators” In _Computational Methods in Applied Mathematics_ 21, 2021, pp. 423 –443
  * [38] Y. Kuznetsov and S. Repin “Mixed Finite Element Method on Polygonal and Polyhedral Meshes” In _Numerical Mathematics and Advanced Applications_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2004, pp. 615–622
  * [39] Y. A. Kuznetsov and S. I. Repin “Convergence analysis and error estimates for mixed finite element method on distorted meshes” In _J. Num. Math._ 13, 2005, pp. 33–51
  * [40] Philip L. Lederer, Alexander Linke, Christian Merdon and Joachim Schöberl “Divergence-free Reconstruction Operators for Pressure-Robust Stokes Discretizations with Continuous Pressure Finite Elements” In _SIAM J. Numer. Anal._ 55, 2017, pp. 1291–1314
  * [41] A. Linke and C. Merdon “On velocity errors due to irrotational forces in the Navier-Stokes momentum balance” In _J. Comput. Phys._ 313, 2016, pp. 654–661 DOI: 10.1016/j.jcp.2016.02.070