<figure><img src="image_1.png" /></figure>

HAL
open science

# Hybrid high-order methods for elliptic PDEs on curved and complicated domains

Zhaonan Dong, Zuodong Wang

► To cite this version:

Zhaonan Dong, Zuodong Wang. Hybrid high-order methods for elliptic PDEs on curved and complicated domains. 13th ICOSAHOM 2020/2021 Conference Proceeding, Jul 2021, Vienna, Austria. hal-03476090v2

HAL Id: hal-03476090

https://inria.hal.science/hal-03476090v2

Submitted on 9 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Hybrid high-order methods for elliptic PDEs on curved and complicated domains

Zhaonan Dong and Zuodong Wang

**Abstract** We introduce a variant of the hybrid high-order method (HHO) employing Nitsche's boundary penalty techniques for the Poisson problem on the curved and complicated Lipschitz domain. The proposed method has two advantages: Firstly, there are no face unknowns introduced on the boundary of the domain, which avoids the computation of the parameterized mapping for the face unknowns on the curved domain boundary. Secondly, using Nitsche's boundary penalty techniques for weakly imposing Dirichlet boundary conditions one can obtain the stability and optimal error estimate independent of the number and measure of faces on the domain boundary. Finally, a numerical experiment is presented in this chapter to confirm the theoretical results.

## 1 Introduction

Hybrid high-order (HHO) methods were introduced for solving the linear diffusion problems in [10] and for locking-free linear elasticity problems in [9]. In the HHO method, the discrete unknowns are attached to the mesh cells, and the mesh faces. The two key ingredients for HHO methods are local reconstruction operators and local stabilization operators in each mesh cell. HHO methods offer various attractive features, such as the support of polytopal meshes, optimal error estimates, local conservation properties, and computational efficiency due to compact stencils and local elimination of the cell unknowns by static condensation. We refer to two monographs [6, 7] for a comprehensive review.

Zhaonan Dong

Inria, 2 rue Simone Iff, 75589 Paris, France, e-mail: zhaonan.dong@inria.fr

CERMICS, École des Ponts, 77455 Marne-la-Vallée, France,

Zuodong Wang

École des Ponts, 77455 Marne-la-Vallée, France. e-mail: zuodong.wang@eleves.enpc.fr
The present work aims to analyze and test a variant of the HHO method employing Nitsche's techniques for solving the Poisson problem with mixed boundary conditions on curved and complicated domains. The complicated domain in this work is denoted as the Lipschitz domain with a possibly arbitrary number of (d-1)-dimensional tiny faces. Using Nitsche's techniques for weakly enforcing the boundary conditions or the interface conditions for HHO methods has been applied for second-order, and fourth-order elliptic PDEs on polygonal and polyhedral meshes cells [2, 4, 10, 11]. The HHO methods mentioned above can deal with the curved interface or domain boundary without computing the parameterized mappings. However, to the best of the authors' knowledge, there is no HHO methods in the literature can solve the Poisson problem on a complicated Lipschitz domain containing a lot of small faces on the Neumann boundary without being influenced by the number of faces. In this work, we focus on the complicated domain with arbitrary number of planar faces only.

There are two main novelties of the HHO methods in this work: Firstly, the proposed HHO methods do not contain any face unknowns on the Dirichlet and Neumann boundary faces. There is no need to compute the parameterized mapping to satisfy the boundary condition on the curved domain. Secondly, the size of the linear system is independent of the number of faces on the domain boundary, which is essential for solving PDEs on a complicated domain with a lot of tiny faces. Although, the computational cost for computing the face integrals still depends on the number of faces on the domain boundary. Finally, we present one numerical example to test the proposed method.

## 2 Model and discrete setting

We follow the standard notations in Sobolev space, $||\cdot||$ means the standard $L^2$ norm and $|\cdot|_{H^{t+1}}$ means the standard $H^{t+1}$ seminorm, and we simplify the notation of $L^2$ inner product as $(\cdot, \cdot)$, if there is no special statement.

### 2.1 Model problem

Let $\Omega$ be a bounded Lipschitz domain in $\mathbb{R}^d$, $d \in \{2, 3\}$, $\Gamma = \Gamma_D \cup \Gamma_N$, where $\Gamma_D$ is the closed Dirichlet boundary satisfying $|\Gamma_D| > 0$ and $\Gamma_N$ is the Neumann boundary. We consider the solution $u \in H^1(\Omega)$ of the following elliptic boundary problem:

$$ -\Delta u = f \quad \text{in} \quad \Omega, \qquad (1) $$

$$ u = g_D \quad \text{on} \quad \Gamma_D, \quad \nabla u \cdot n_{\Omega} = g_N \quad \text{on} \quad \Gamma_N $$

with $f \in L^2(\Omega)$, $g_D \in H^{\frac{1}{2}}(\Gamma_D)$, $g_N \in L^2(\Gamma_N)$ and $n_{\Omega}$ to be the unit outward normal vector of domain $\Omega$.
The weak form of the above problem is defined as: Find $u \in H^1(\Omega)$, satisfying $u|_{\Gamma_D} = g_D$,

$$a(u, v) = \ell(v), \quad \forall v \in H_{\Gamma_D}^1(\Omega),$$

with $a(u, v) = (\nabla u, \nabla v)_\Omega$ and $\ell(v) = (f, v)_\Omega + (g_N, v)_{\Gamma_N}$, where $H_{\Gamma_D}^1(\Omega) = \{v \in H^1(\Omega)|v=0 \text{ on } \Gamma_D\}$.

## 2.2 Polytopal and curved meshes

Let $\{\mathcal{T}_h\}_{h>0}$ be a mesh family such that each mesh $\mathcal{T}_h$ covers $\Omega$ exactly. A generic mesh cell is denoted by $T \in \mathcal{T}_h$, its diameter by $h_T$, and its unit outward normal vector by $n_T$. We partition the boundary $\partial T$ of any mesh cell $T \in \mathcal{T}_h$ by means of the two subsets $\partial T^i = \partial T \cap \bar{\Omega}$, $\partial T^b = \partial T \cap \partial \Omega$. Similarly, we partition the mesh as $\mathcal{T}_h = \mathcal{T}_h^i \cup \mathcal{T}_h^b$, where $\mathcal{T}_h^b$ is the collection of all the mesh cells $T$ such that $\partial T^b$ has positive measure. Moreover, we further split $\partial T^b$ into two nonoverlapping subsets $\partial T^D = \partial T^b \cap \Gamma_D$ and $\partial T^N = \partial T^b \cap \Gamma_N$. The mesh faces are collected in the set $\mathcal{F}_h$, which is split as $\mathcal{F}_h = \mathcal{F}_h^i \cup \mathcal{F}_h^b$, where $\mathcal{F}_h^i$ is the collection of the interior faces (shared by two distinct mesh cells) and $\mathcal{F}_h^b$ the collection of the boundary faces. For all $F \in \mathcal{F}_h$, we orient $F$ by means of the fixed unit normal vector $n_F$ whose direction is arbitrary for all $F \in \mathcal{F}_h^i$ and $n_F = n_\Omega$ for all $F \in \mathcal{F}_h^b$. For any mesh cell $T \in \mathcal{T}_h$, the mesh faces composing its boundary $\partial T$ are collected in the set $\mathcal{F}_{\partial T}$, which is partitioned as $\mathcal{F}_{\partial T} = \mathcal{F}_{\partial T^i} \cup \mathcal{F}_{\partial T^D} \cup \mathcal{F}_{\partial T^N}$ with obvious notation.

**Assumption (1).** Any interior element $T \in \mathcal{T}_h^i$ is a polytope with planar faces, and the sequence of interior meshes $\mathcal{T}_h^i$ is shape-regular in the sense of [9].

(2). For any boundary mesh cell $T \in \mathcal{T}_h^b$, all the faces in $\mathcal{F}_{\partial T^i}$ are planar with diameter uniformly equivalent to $h_T$. Moreover, for each $\partial T^b$, $T$ can be decomposed into a finite union of nonoverlapping subsets, $\{T_{\partial T^b,m}\}_{m \in \{1, \dots, n_{T, \partial T^b}\}}$, such that $T_{\partial T^b,m}$ is star-shaped with respect to an interior ball of radius $r_{T,m}$ satisfying $\rho h_T \le r_{T,m}$, with the mesh-regularity parameter $\rho > 0$. We refer to Figure 1 for an illustration with $n_{T, \partial T^i} = 1$ (left) and $n_{T, \partial T^i} = 3$ (right). $\square$

**Remark 1 (Mesh assumptions)** The above mesh assumptions on the mesh sequence are pretty general. In this work, we assume the interior faces of the mesh are planar and shape-regular. The reason is that those face unknowns on interior faces will be coupled globally after the static condensation, and we try to avoid using the meshes with an unbounded number of tiny faces in the interior of the domain. However, HHO methods have been proven to be stable and converge optimally on the meshes with an arbitrary number of tiny faces under certain star-shapedness assumptions. We refer to the recent work [13] for a detailed discussion.

On the other hand, using a Nitsche-like penalty method avoids introducing discrete face unknowns on the boundary faces. Therefore, the face on the domain boundary can be curved or complicated. The star-shapedness assumption (2) on the mesh
<figure><img src="image_3.png" /><figcaption>Fig. 1 Examples for the mesh assumptions</figcaption></figure>

boundary cells is introduced to invoke a Poincaré-type inequality, discrete trace inequality and multiplicative trace inequality to hold in such cells.

**Remark 2 (Limitation of the domain)** In this work, the complicated domain is denoted as the Lipschitz domain with a possibly arbitrary number of $(d-1)$-dimensional tiny planar faces. Following the technical result in [1, Theorem 2.1], a Lipschitz domain can be covered by a finite number of subdomains that are star-shaped with respect to an open ball. This result implies that the above mesh assumption is not very restrictive for the complicated Lipschitz domains. However, if the domain is non-Lipschitz, the above meshes assumption may fail to hold on to such domains.

## 2.3 Analysis tools

Let us briefly review the main analysis tools used in this section. We only state the results. In what follows, we always consider a shape-regular mesh sequence satisfying the above assumptions. Moreover, in various bounds, we use $a \le b$ to denote $a \le Cb$ with $C$ to be any generic constant (its value can change at each occurrence) that is independent of mesh size $h > 0$, but may depend on the shape-regularity of the mesh sequence and the polynomial degree for all $T \in \mathcal{T}_h$. Let $k \in \mathbb{N}$ be the polynomial degree and $\mathbb{P}_{k,d}(T)$ to be the polynomial space of order at most $k$ in any cell $T \in \mathcal{T}_h$.

**Lemma 1 (Discrete trace inequalities)**

For any $T \in \mathcal{T}_h$ satisfying the mesh assumptions, the following relation holds for all $v \in \mathbb{P}_{k,d}(T)$:

$$ \|v\|_{\partial T} \le h_T^{-\frac{1}{2}} \|v\|_T. \qquad (2) $$

**Lemma 2 (Multiplicative trace inequality)**

Let element $T \in \mathcal{T}_h$. Then, the following statement holds true for all $v \in H^1(\Omega)$:
$$ \|v\|_{\partial T} \lesssim (h_T^{-\frac{1}{2}} \|v\|_T + \|v\|_T^{\frac{1}{2}} \|\nabla v\|_T^{\frac{1}{2}}). \quad (3) $$

**Lemma 3 (Polynomial approximation)**

Let $k \ge 0$ be the polynomial degree, and there is a real number $t \in [0, k+1]$, all $m \in \{0, \dots, \lfloor t \rfloor\}$, and $v \in H^t(T)$,

$$ |v - \Pi_T^k(v)|_{H^m(T)} \lesssim h_T^{t-m} |v|_{H^t(T)}, \quad (4) $$

where $\Pi_T^k$ denoted the $L^2$-orthogonal projection onto $\mathbb{P}_{k,d}(T)$.

The following approximation result holds by using (4) combined with the multiplicative trace inequality (3): Let $k \ge 0$ be the polynomial degree, the following estimate holds for all $v \in H^t(T)$,

$$ \|v - \Pi_T^k(v)\|_T + h_T^{\frac{1}{2}} \|v - \Pi_T^k(v)\|_{\partial T} + h_T \|\nabla(v - \Pi_T^k(v))\|_T \\ + h_T^{\frac{3}{2}} \|\nabla(v - \Pi_T^k(v))\|_{\partial T} \lesssim h_T^t |v|_{H^t(T)} \quad (5) $$

**Remark 3** We briefly comment on the proofs of the above lemmas. For the Lemma 1 and Lemma 2, the proof on mesh cells having flat faces can be found in [8, Sec. 1.4.3]. On the mesh cells having a curved face, these results can be found in [3, 17] assuming that the curved face is a $C^2$ manifold. More recently, these results were extended in [4] with fully explicit constants to Lipschitz manifolds satisfying the mild additional geometric assumptions.

Concerning Lemma 3, the key step is to establish the Poincaré inequality since (4) can then be derived by using recursively the Poincaré inequality. On the interior mesh cells, which can be decomposed as a finite union of subsimplices, this latter inequality is established by proceeding as in [16, 14]. On the boundary mesh cells, which can have a curved face, one invokes the star-shapedness assumption with respect to a ball. We refer the reader to [18] for the derivation of this inequality with an explicitly determined constant under such an assumption.

## 3 HHO discretization

In this section, we will introduce the mixed-order HHO method employing Nitsche's boundary penalty techniques. We first start to introduce the HHO space. Let $k \ge 0$ be the polynomial degree. We consider a pair $\hat{v}_T^k = (v_T^{k+1}, v_{\partial T^i}^k)$, where $v_T^{k+1}$ is defined on $T$ and $v_{\partial T^i}$ is defined on the inner faces (facewise) $F \in \mathcal{F}_{\partial T^i}$ composing the boundary $\partial T^i$ of $T$. We define the local HHO space

$$ \hat{v}_T^k = \mathbb{P}_{k+1,d}(T) \times \mathbb{P}_{k,d-1}(\mathcal{F}_{\partial T^i}), $$

noting the above local HHO space does not contain face unknowns on the domain boundary.
Next, we denote by $\partial_n$ the (scalar-valued) outward normal derivative on $\partial T$ for all $T \in \mathcal{T}_h$. The construction of the HHO methods are as follows: We start with the reconstruction of the gradient operator and then add the stabilization operator. The reconstruction operator is defined as: $R_T^i: \hat{V}_T^k \to \mathbb{P}_{k+1,d}(T)$ s.t. for every pair $\hat{v}_T^k = (v_T^{k+1}, v_{\partial T^i}^k) \in \hat{V}_T^k$, the function $R_T^i(\hat{v}_T^k)$ is s.t. for all $q \in \mathbb{P}_{k+1,d}(T)$,

$$ (\nabla R_T^i(\hat{v}_T^k), \nabla q)_T := (\nabla v_T^{k+1}, \nabla q)_T - (v_T^{k+1} - v_{\partial T^i}^k, \partial_n q)_{\partial T^i} - (v_T^{k+1}, \partial_n q)_{\partial T^D} \\ = -(v_T^{k+1}, \Delta q)_T + (v_{\partial T^i}^k, \partial_n q)_{\partial T^i} + (v_T^{k+1}, \partial_n q)_{\partial T^N} \quad (6) $$

together with $(R_T^i(\hat{v}_T^k) - v_T^{k+1}, 1)_T = 0$. In addition, we define the lifting operator $L_T^{k+1}: L^2(\partial T^D) \to \mathbb{P}_{k+1,d}(T)$ for all $T \in \mathcal{T}_h$ such that, for all $g_D \in L^2(\partial T^D)$ and all $q \in \mathbb{P}_{k+1,d}(T)$,

$$ (\nabla L_T^{k+1}(g_D), \nabla q)_T = (g_D, \partial_n q)_{\partial T^D}, \quad (7) $$

together with the condition $(L_T^{k+1}(g_D), 1)_T = 0$. Noticing that $L_T^{k+1}(g_D) = 0$ for all interior cell $T \in \mathcal{T}_h^i$.

Next, we consider the stabilization operators. For the interior faces $F \in \mathcal{F}_{\partial T^i}$, the LS stabilization operators are defined as

$$ S_T^i(\hat{v}_T^k) = \Pi_{\partial T^i}^k (v_T^{k+1} - v_{\partial T^i}^k). $$

where $\Pi_{\partial T^i}^k$ denotes the $L^2$-orthogonal projection onto the broken polynomial space $\mathbb{P}_{k,d-1}(\mathcal{F}_{\partial T^i})$.

For the Dirichlet face $F \in \mathcal{F}_{\partial T^D}$, we employ Nitsche's boundary penalty techniques, i.e., $S_T^D(\hat{v}_T^k) = v_T^{k+1}|_{\partial T^D}$.

Then, the local bilinear form $\hat{a}_T$ and local linear form $\ell_T$ in each cell $T$ is defined as follows:

$$ \hat{a}_T(\hat{u}_T^k, \hat{v}_T^k) = (\nabla R_T^i(\hat{u}_T^k), \nabla R_T^i(\hat{v}_T^k))_T \\ + h_T^{-1}(S_T^i(\hat{u}_T^k), S_T^i(\hat{v}_T^k))_{\partial T^i} + h_T^{-1}(u_T^{k+1}, v_T^{k+1})_{\partial T^D}, \\ \ell_T(\hat{v}_h^k) = (f, v_T^{k+1})_T + (g_N, v_T^{k+1})_{\partial T^N} - (g_D, \partial_n R_T^i(\hat{v}_T^k) - h_T^{-1}v_T^{k+1})_{\partial T^D}, $$

for all $\hat{u}_T^k, \hat{v}_T^k$ in $\hat{V}_T^k$.

Next, let us define the global HHO space as

$$ \hat{V}_h^k = \mathbb{P}_{k+1,d}(\mathcal{T}_h) \times \mathbb{P}_{k,d-1}(\mathcal{F}_h^i), $$

then the global bilinear forms $\hat{a}_h$, and linear form $\ell_h$ as follows:

$$ \hat{a}_h(\hat{u}_h^k, \hat{v}_h^k) = \sum_{T \in \mathcal{T}_h} \hat{a}_T(\hat{u}_T^k, \hat{v}_T^k), \\ \ell_h(\hat{v}_h^k) = \sum_{T \in \mathcal{T}_h} \ell_T(\hat{v}_h^k), \quad (8) $$
for all $\hat{u}_h^k, \hat{v}_h^k$ in $\hat{V}_h^k$

Finally, the discrete problem can be defined as follows:

$$ \left\{ \begin{array}{l} \text{Find } \hat{u}_h^k \in \hat{V}_h^k \text{ such that} \\ \hat{a}_h(\hat{u}_h^k, \hat{v}_h^k) = \ell_h(\hat{v}_h^k), \forall \hat{v}_h^k \in \hat{V}_h^k, \end{array} \right. \qquad (9) $$

*Remark 4* We mention that there is another version for HHO method with Nitsche's boundary penalty techniques which does not contain the lifting operator, see section 2 in [5]. Without the lifting operator, the stability of the method in [5] is achieved by choosing the penalty parameter large enough. On the contrary, the proposed HHO methods do not suffer from it.

# 4 Main results

## 4.1 Stability and well-posedness

The stability and approximation results presented in this section follows the same technique as in [15, Chapter 39]. For simplification of the presentation, the details of some proofs are not presented in this work.

For all $T \in \mathcal{T}_h$ and all $\hat{v}_h^k \in \hat{V}_h^k$, we define the local seminorm as:

$$ |\hat{v}_T^k|^2_{\hat{V}_T^k} = \|\nabla v_T^{k+1}\|_T^2 + h_T^{-1} \|v_{\partial T^i}^k - v_T^{k+1}\|_{\partial T^i}^2 + h_T^{-1} \|v_T^{k+1}\|_{\partial T^D}^2. $$

In addition, we define the norm for the space $\hat{V}_h^k$ as

$$ \|\hat{v}_h^k\|_{\hat{V}_h^k}^2 = \sum_{T \in \mathcal{T}_h} |\hat{v}_T^k|^2_{\hat{V}_T^k}. $$

**Lemma 4 (Local stability and boundedness)**

Let all $T \in \mathcal{T}_h$ satisfy the mesh assumptions. Then, the local stability holds

$$ |\hat{v}_T^k|^2_{\hat{V}_T^k} \lesssim \hat{a}_T(\hat{v}_T^k, \hat{v}_T^k) \lesssim |\hat{v}_T^k|^2_{\hat{V}_T^k}. \qquad (10) $$

*Proof* The proof follows the same techniques as in [6, Sec. 2.2]. □

An immediate consequence of Lemma 4 is the following bound establishing that the discrete bilinear form $\hat{a}_h$ is coercive on $\hat{V}_h^k$:

$$ \|\hat{v}_h^k\|_{\hat{V}_h^k}^2 \lesssim \hat{a}_h(\hat{v}_h^k, \hat{v}_h^k), \quad \forall \hat{v}_h^k \in \hat{V}_h^k $$

Invoking the Lax–Milgram lemma the (9) is well-posed.
## 4.2 Approximation and $H^1$ error estimate

For any $T \in \mathcal{T}_h$, we define the local reduction operator $\hat{I}_T^k : H^1(T) \to \mathbb{P}_{k+1,d}(T)$, for all $v \in H^1(T)$,

$$ \hat{I}_T^k(v) = (\Pi_T^{k+1}(v), \Pi_{\partial T^i}^k(v)) \in \hat{V}_T^k. $$

We then define the energy projection operator $\mathcal{E}_T : H^1(T) \to \mathbb{P}_{k+1,d}(T)$ such that

$$ \mathcal{E}_T(v) = R_T^i \circ \hat{I}_T^k(v) + L_T^{k+1}(v), \quad \forall v \in H^1(T). $$

Similarly, the global energy projection is defined by $\mathcal{E}(v)|_T = \mathcal{E}_T(v)$, for all $T \in \mathcal{T}_h$.

**Lemma 5 (Energy projection)**

Let $\mathcal{E}_T$ be the energy projection. Then, the following relation holds,

$$ (\nabla(\mathcal{E}_T(v) - v), \nabla q)_T = (\Pi_T^{k+1}(v) - v, \partial_n q)_{\partial T^N}, \quad (11) $$

and the condition $(\mathcal{E}_T(v), 1)_T = (v, 1)_T$ for all $q \in \mathbb{P}_{k+1,d}(T)$ and $v \in H^1(T)$.

**Proof** Let $v \in H^1(T)$ and $\phi = R_T^i \circ \hat{I}_T^k(v)$. Using the definition of the reconstruction operator (6), we infer that for all $q \in V_T^{k+1}$,

$$ \begin{aligned} (\nabla \phi, \nabla q)_T &= -(\Pi_T^{k+1}(v), \Delta q)_T + (\Pi_{\partial T^i}^k(v), \partial_n q)_{\partial T^i} + (\Pi_T^{k+1}(v), \partial_n q)_{\partial T^N} \\ &= -(v, \Delta q)_T + (v, \partial_n q)_{\partial T^i} + (v, \partial_n q)_{\partial T^N} + (\Pi_T^{k+1}(v) - v, \partial_n q)_{\partial T^N} \\ &= (\nabla v, \nabla q)_T - (v, \partial_n q)_{\partial T^D} + (\Pi_T^{k+1}(v) - v, \partial_n q)_{\partial T^N}, \end{aligned} $$

since $\partial_n q|_F$ is a polynomial of order $k$ on $F \in \partial T^i$. Then, using the definition of lifting operator, we have $(v, \partial_n q)_{\partial T^D} = (\nabla L_T^{k+1}(v), \nabla q)_T$, this gives $(\nabla(\mathcal{E}_T(v) - v), \nabla q)_T = (\Pi_T^{k+1}(v) - v, \partial_n q)_{\partial T^N}$. Moreover, $(R_T^i \circ \hat{I}_T^i(v) + L_T^{k+1}(v), 1)_T = (\Pi_T^{k+1}(v), 1)_T = (v, 1)_T$ by the definition of the reconstruction operator and the lifting operator. $\square$

**Remark 5** By noticing the term $(\Pi_T^{k+1}(v) - v, \partial_n q)_{\partial T^N}$ depends on the Neumann boundary, $\mathcal{E}_T(v)$ is the elliptic projection for the function $v$ on the cells which does not contain any face on the Neumann boundary.

We start with the approximation properties of operators.

**Lemma 6 (Approximation results)**

For all $v \in H^{t+1}(\Omega)$, $t > \frac{1}{2}$ and $v|_T \in H^{k+2}(T)$ for all $T \in \mathcal{T}_h$, we have

$$ \|\nabla(\mathcal{E}_T(v) - v)\|_T + h_T^{\frac{1}{2}} \|\nabla(\mathcal{E}_T(v) - v)\|_{\partial T^i \cup \partial T^D} \lesssim h_T^{k+1} |v|_{H^{k+2}(T)}, \quad (12) $$

$$ h_T^{-\frac{1}{2}} \|S_T^i(\hat{I}_T^k(v))\|_{\partial T^i} + h_T^{-\frac{1}{2}} \|S_T^D(\hat{I}_T^k(v)) - v\|_{\partial T^D} \lesssim h_T^{k+1} |v|_{H^{k+2}(T)}. \quad (13) $$

**Proof** Invoking the discrete trace inequality (2), multiplicative trace inequality (3) and approximation result (4) The above results follow the same steps as in [6, Sec. 2.3]. $\square$
We define the consistency error $\delta_h \in (\hat{V}_h^k)'$ which is the dual space of $\hat{V}_h^k$ such that, for all $\hat{v}_h^k \in \hat{V}_h^k$,

$$ \langle \delta_h, \hat{v}_h^k \rangle = \sum_{T \in \mathcal{T}_h} (\hat{a}_T^k (\hat{I}_T^k(u), \hat{v}_T^k) - \ell_T(\hat{v}_T^k)), \quad (14) $$

where the brackets refer to the duality pairing between $(\hat{V}_h^k)'$ and $\hat{V}_h^k$.

**Lemma 7 (Consistency and boundedness)**

Assume that $u \in H^{l+1}(\Omega)$ with $t > \frac{1}{2}$ and $\hat{v}_h^k \in \hat{V}_h^k$. The following holds true

$$ |\langle \delta_h, \hat{v}_h^k \rangle| \le \left( \sum_{T \in \mathcal{T}_h} \|g_T\|_{*T}^2 + \|\xi_T\|_{\#T}^2 \right)^{\frac{1}{2}} \times \|\hat{v}_h^k\|_{\hat{V}_h^k}, $$

where

$$ \begin{aligned} \|g_T\|_{*T}^2 &= \|\nabla(u - \mathcal{E}(u))\|_T^2 + h_T \|\nabla(u - \mathcal{E}(u))\|_{\partial T^i \cup \partial T^D}^2, \\ \|\xi_T\|_{\#T}^2 &= h_T^{-1} \|u - \Pi_T^{k+1}(u)\|_{\partial T^i}^2 + h_T^{-1} \|u - \Pi_T^{k+1}(u)\|_{\partial T^D}^2. \end{aligned} $$

**Proof** Let $\hat{v}_h^k \in \hat{V}_h^k$. Recalling that $-\Delta u = f$, we obtain $\langle \delta_h, \hat{v}_h^k \rangle = \Psi_1 + \Psi_2$, where

$$ \Psi_1 = \sum_{T \in \mathcal{T}_h} \left( (\nabla R_T^i \circ \hat{I}_T^i(u), \nabla R_T^i(\hat{v}_T^k))_T + (\Delta u, v_T^{k+1})_T \right. \\ \left. \qquad + (\nabla L_T^{k+1}(g_D), \nabla R_T^i(\hat{v}_T^k))_T - (g_N, v_T^{k+1})_{\partial T^N} \right), \quad (15) $$

$$ \Psi_2 = \sum_{T \in \mathcal{T}_h} (h_T^{-1}(S_T^i \circ \hat{I}_T^k(u), S_T^i(\hat{v}^k q w_T))_{\partial T^i} + h_T^{-1}(\Pi_T^{k+1}(u) - g_D, v_T^{k+1})_{\partial T^D}). \quad (16) $$

Integrating by parts in every mesh cell $T \in \mathcal{T}_h$, using the exact solution of PDE, $\sum_{T \in \mathcal{T}_h} (\nabla u \cdot n_T, v_{\partial T^i}^k)_{\partial T^i} = 0$, since $\partial_n u$ is meaningful and single valued on all interior faces $F \in \partial T^i$, and the definition of the reconstruction operator, we have

$$ \Psi_1 = \sum_{T \in \mathcal{T}_h} \left( -(\nabla(\mathcal{E}(u) - u) \cdot n_T, v_T^{k+1})_{\partial T^D} \right. \\ \left. \qquad + (\nabla(\mathcal{E}(u) - u), \nabla v_T^{k+1})_T + (\nabla(\mathcal{E}(u) - u) \cdot n_T, v_{\partial T^i}^k - v_T^{k+1})_{\partial T^i} \right). $$

Noticing $u = g_D$ on $\Gamma_D$, we have

$$ \Psi_2 = \sum_{T \in \mathcal{T}_h} (h_T^{-1}(\Pi_{\partial T^i}^k(u - \Pi_T^{k+1}(u)), v_{\partial T^i}^k - v_T^{k+1})_{\partial T^i} + h_T^{-1}(\Pi_T^{k+1}(u) - u, v_T^{k+1})_{\partial T^D}). $$

Then, the Cauchy-Schwarz inequality and discrete trace inequality (2) finish this proof. □

**Theorem 1 (Broken $H^1$ error estimate)**
In the framework of Lemma 7, with the $\hat{u}_T^k$ to be the solution of the discrete problem (9) on each cell $T \in \mathcal{T}_h$, we have

$$ \sum_{T \in \mathcal{T}_h} \|\nabla(u|_T - R_T^i(\hat{u}_T^k) - L_T^{k+1}(g_D))\|_T^2 \le \sum_{T \in \mathcal{T}_h} (\|g_T\|_{*T}^2 + \|\xi_T\|_{\#T}^2). \quad (17) $$

Moreover, if $u|_T \in H^{k+2}(T)$ for all $T \in \mathcal{T}_h$, with $k \ge 0$, we have

$$ \sum_{T \in \mathcal{T}_h} \|\nabla(u|_T - R_T^i(\hat{u}_T^k) - L_T^{k+1}(g_D))\|_T^2 \le \sum_{T \in \mathcal{T}_h} (h_T^{2(k+1)}|u|_{H^{k+2}(T)}^2). \quad (18) $$

**Proof** First, let us set $\hat{e}_h^k = \hat{I}_h^k(u) - \hat{u}_h^k \in \hat{V}_h^k$, so that $\hat{e}_T^k = \hat{I}_T^k(u|_T) - \hat{u}_T^k \in \hat{V}_T^k$ for all $T \in \mathcal{T}_h$. The property (10) and the identity $a_h(\hat{e}_h^k, \hat{e}_h^k) = -\langle \delta_h, \hat{e}_h^k \rangle$ imply that

$$ \sum_{T \in \mathcal{T}_h} \|\nabla(R_T^i(\hat{e}_T^k))\|_T^2 \le \|\hat{e}_h^k\|_{\hat{V}_h^k}^2 \le a_h(\hat{e}_h^k, \hat{e}_h^k) = -\langle \delta_h, \hat{e}_h^k \rangle. $$

Owing to the Lemma 7, we infer that

$$ \sum_{T \in \mathcal{T}_h} \|\nabla(R_T^i(\hat{e}_T^k))\|_T^2 \le \sum_{T \in \mathcal{T}_h} (\|g_T\|_{*T}^2 + \|\xi_T\|_{\#T}^2). $$

By adding and subtracting $R_T^i(\hat{I}_T^k(u|_T))$, we have

$$ u|_T - R_T^i(\hat{u}_T^k) - L_T^{k+1}(u|_T) = (u|_T - \mathcal{E}_T(u|_T)) + R_T^i(\hat{e}_T^k). $$

Using the triangle inequality, the error estimate (17) is derived. Furthermore, the error bound (18) results from (17) and the approximation results in Lemma 6. □

## 5 Numerical experiment

We will test the proposed HHO-N methods on the domain containing a lot of tiny faces. The computational domain is constructed as follows: 1). Construct a unit disc centred at origin, with a circular hole centred at origin and radius 0.3 removed. 2). Remove the domain on the right-hand side of the curve $x = 0.05 \sin(4\pi y - 4\pi) + 0.7$. Finally, randomly perturb the 4010 boundary vertices with the amplitude less or equal to 0.001. This results a polygon centred at (0,0) with saw-like boundary consisting of 4010 small edges which does not embedded on any smooth manifold, see Figure 2. The HHO-N method is tested with the exact solution $u = \sin(\pi x) \sin(\pi y)$, with Dirichlet boundary condition imposed on 753 boundary edges allocated inside the domain $-1 \le x \le -0.8$, and Neumann boundary condition on the remaining boundary edges.

We consider a quasi-uniform sequence of meshes composed of 217, 515, 2058, 5455, and 8992 triangular cells. All meshes fit the domain $\Omega$ exactly, and for every
<figure><img src="image_4.png" /><figcaption>Fig. 2 Complicated domain</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3 Convergence rates in the $L^2$-norm and energy seminorm for $k = 0, 1, 2, 3$</figcaption></figure>

mesh in the sequence, each interior cell has only straight edges, whereas each boundary cell has a lot of small edges that exactly fit the boundary of $\Omega$. The convergence result is given in Figure 3. The key observation is that the convergence rates are $O(h^{k+2})$ to the $O(h^{k+1})$ for the $L^2$ error and the energy error, respectively. We emphasize that the total number of global coupled unknowns for the proposed HHO methods is independent of the number of edges on the domain boundary.
Acknowledgments

The authors would like to thank Professor Alexandre Ern (École des Ponts & Inria, France) for his valuable advice and helpful discussion.

References

1. C. Amrouche, P. G. Ciarlet, and C. Mardare. On a lemma of Jacques-Louis Lions and its relation to other fundamental results. *J. Math. Pures Appl.* (9), 104(2):207–226, 2015.

2. E. Burman, M. Cicuttin, G. Delay, and A. Ern. An unfitted hybrid high-order method with cell agglomeration for elliptic interface problems. *SIAM J. Sci. Comput.*, 43(2):A859–A882, 2021.

3. E. Burman and A. Ern. An unfitted hybrid high-order method for elliptic interface problems. *SIAM J. Numer. Anal.*, 56(3):1525–1546, 2018.

4. A. Cangiani, Z. Dong, and E. H. Georgoulis. hp-version discontinuous Galerkin methods on essentially arbitrarily-shaped elements. *Math. Comp.*, 91(333):1–35, 2022.

5. K. L. Cascavita, F. Chouly, and A. Ern. Hybrid high-order discretizations combined with nitsche's method for dirichlet and signorini boundary conditions. *IMA J. Numer. Anal.*, 40(4):2189–2226, 2020.

6. M. Cicuttin, A. Ern, and N. Pignet. Hybrid high-order methods. a primer with application to solid mechanics. *SpringerBriefs in Mathematics*, 2021.

7. D. A. Di Pietro and J. Droniou. *The hybrid high-order method for polytopal meshes*, volume 19 of MS&A. *Modeling, Simulation and Applications*. Springer, Cham, [2020] ©2020. Design, analysis, and applications.

8. D. A. Di Pietro and A. Ern. *Mathematical Aspects of Discontinuous Galerkin Methods*, volume 69. Springer Berlin Heidelberg, Berlin, Heidelberg, 2012.

9. D. A. Di Pietro and A. Ern. A Hybrid High-Order locking-free method for linear elasticity on general meshes. *Comput. Meth. Appl. Mech. Engrg.*, 283:1–21, 2015.

10. D. A. Di Pietro, A. Ern, and S. Lemaire. An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators. *Comput. Meth. Appl. Math.*, 14(4):461–472, 2014.

11. Z. Dong and A. Ern. Hybrid high-order and weak Galerkin methods for the biharmonic problem. *arXiv preprint arXiv:2103.16404*. Accepted in *SIAM J. Numer. Anal*.

12. Z. Dong and A. Ern. Hybrid high-order method for singularly perturbed fourth-order problems on curved domains. *ESAIM Math. Model. Numer. Anal.*, 55(6):3091–3114, 2021.

13. L. Droniou, J. and Yemm. Robust hybrid high-order method on polytopal meshes with small faces. *Comput. Methods Appl. Math.*, 22(1):47–71, 2022.

14. A. Ern and J.-L. Guermond. Finite element quasi-interpolation and best approximation. *ESAIM Math. Model. Numer. Anal. (M2AN)*, 51(4):1367–1385, 2017.

15. A. Ern and J.-L. Guermond. *Finite Elements II: Galerkin Approximation, Elliptic and Mixed PDEs*, volume 73 of *Texts in Applied Mathematics*. Springer Nature, Cham, Switzerland, 2021.

16. A. Veeser and R. Verfürth. Poincaré constants for finite element stars. *IMA J. Numer. Anal.*, 32(1):30–47, 2012.

17. H. Wu and Y. Xiao. An unfitted hp-interface penalty finite element method for elliptic interface problems. *J. Comput. Math.*, 37(3):316–339, 2019.

18. W. Zheng and H. Qi. On Friedrichs-Poincaré-type inequalities. *J. Math. Anal. Appl.*, 304(2):542–551, 2005.