A Hybrid High-Order method for Kirchhoff-Love plate bending problems

TL;DR

This paper introduces a new Hybrid High-Order method for accurately solving fourth-order Kirchhoff-Love plate bending problems on polygonal meshes, with proven convergence and robustness supported by numerical tests.

Contribution

The paper develops a novel HHO discretization supporting arbitrary orders and polygonal meshes, with proven convergence and error estimates for plate bending problems.

Findings

Convergence rate of h^{k+1} in energy norm

Error estimate of h^{k+3} in L^2 norm under regularity

Numerical experiments confirm robustness and accuracy

Abstract

We present a novel Hybrid High-Order (HHO) discretization of fourth-order elliptic problems arising from the mechanical modeling of the bending behavior of Kirchhoff-Love plates, including the biharmonic equation as a particular case. The proposed HHO method supports arbitrary approximation orders on general polygonal meshes, and reproduces the key mechanical equilibrium relations locally inside each element. When polynomials of degree $k \ge 1$ are used as unknowns, we prove convergence in $h^{k+1}$ (with $h$ denoting, as usual, the meshsize) in an energy-like norm. A key ingredient in the proof are novel approximation results for the energy projector on local polynomial spaces. Under biharmonic regularity assumptions, a sharp estimate in $h^{k+3}$ is also derived for the $L^2$-norm of the error on the deflection. The theoretical results are supported by numerical experiments, which additionally show the robustness of the method with respect to the choice of the stabilization.

Tables

Table 1. Table 1: Convergence of E ​ ( u ¯ h ) E subscript ¯ u h \rm E(\underline{\mathrm{u}}_{h}) with four uniform mesh refinements for each polynomial degree k ∈ { 1 , 2 , 3 , 4 } 𝑘 1 2 3 4 k\in\{1,2,3,4\} . The number of triangular elements is given by N 𝑁 N .

	$N=56$	$N=224$	$N=896$	$N=3584$
$k=1$	-1.662960e-03	-1.635846e-03	-1.632895e-03	-1.632670e-03
$k=2$	-1.637918e-03	-1.632707e-03	-1.632652e-03	-1.632653e-03
$k=3$	-1.632412e-03	-1.632634e-03	-1.632652e-03	-1.632653e-03
$k=4$	-1.632433e-03	-1.632638e-03	-1.632652e-03	-1.632653e-03

Table 2. Table 2: Convergence of E ​ ( u ¯ h ) E subscript ¯ u h \rm E(\underline{\mathrm{u}}_{h}) with uniform mesh refinements for each polynomial degree k ∈ { 1 , 2 , 3 , 4 } 𝑘 1 2 3 4 k\in\{1,2,3,4\} . The number of triangular elements is given by N 𝑁 N .

	$N=34$	$N=136$	$N=544$	$N=2176$	$N=8704$
$k=1$	-4.208744e-05	-3.071276e-05	-2.885137e-05	-2.833621e-05	-2.813136e-05
$k=2$	-3.167765e-05	-2.945556e-05	-2.858722e-05	-2.824400e-05	-2.809123e-05
$k=3$	-2.944060e-05	-2.868230e-05	-2.828896e-05	-2.811198e-05	-2.803012e-05
$k=4$	-2.899953e-05	-2.845505e-05	-2.818988e-05	-2.806669e-05	-2.800896e-05

Table 3. Table 3: Size of the matrix associated with a h subscript a ℎ \mathrm{a}_{h} with uniform mesh refinements for each polynomial degree k ∈ { 1 , 2 , 3 , 4 } 𝑘 1 2 3 4 k\in\{1,2,3,4\} . The number of triangular elements is given by N 𝑁 N .

	$N=34$	$N=136$	$N=544$	$N=2176$	$N=8704$
$k=1$	354	1320	5088	19968	79104
$k=2$	531	1980	7632	29952	118656
$k=3$	708	2640	10176	39936	158208
$k=4$	885	3300	12720	49920	197760

Table 4. Table 4: Number of nonzero elements of the matrix associated with a h subscript a ℎ \mathrm{a}_{h} with uniform mesh refinements for each polynomial degree k ∈ { 1 , 2 , 3 , 4 } 𝑘 1 2 3 4 k\in\{1,2,3,4\} . The number of triangular elements is given by N 𝑁 N .

	$N=34$	$N=136$	$N=544$	$N=2176$	$N=8704$
$k=1$	9468	37296	148032	589824	2354688
$k=2$	21303	83916	333072	1327104	5298048
$k=3$	37872	149184	592128	2359296	9418752
$k=4$	59175	233100	925200	3686400	14716800

Equations

$$-\text{mbox}\mathrm{div}\text{mbox}div\mathbf{M}$$

$$u$$

$${\partial }_{\mathbf{n}}u$$

$$\mathbf{M}=-\mathbb{A}{\nabla }^{2}u,$$

$$P_{\Omega }=\{{\Omega }_i:i\in I\}$$

$$(\mathbb{A}{\nabla }^{2}u,{\nabla }^{2}v)-:a(u,v)=(f,v)\forall v\in H_0^2(\Omega ).$$

$${\mathbb{A}}_T:-{\mathbb{A}}_{|T}\forall T\in {\mathcal{T}}_h.$$

$$H^l({\mathcal{T}}_h):-\left\{v\in L^2(\Omega ):v_{|T}\in H^l(T)\forall T\in {\mathcal{T}}_h\right\},$$

$$\forall v\in H^l({\mathcal{T}}_h),\Vert v{\Vert }_{H^l({\mathcal{T}}_h)}:-{\left(\sum _{T\in {\mathcal{T}}_h}\Vert v{\Vert }_{H^l(T)}^2\right)}^{\text{nicefrac}12}.$$

$${\varrho }^2{h}_T\le h_F\le h_T,$$

$$\underset{h\in \mathcal{H}}{\max }\underset{T\in {\mathcal{T}}_h}{\max }\mathrm{card}({\mathcal{F}}_T)\le N_{\partial }.$$

$$\Vert w{\Vert }_F$$

$$h_T^{\text{nicefrac}12}\Vert w{\Vert }_{\partial T}$$

$$\Vert \nabla w{\Vert }_T$$

$$\Vert w{\Vert }_T\le C_{\mathrm{p}}{h}_T\Vert \nabla w{\Vert }_T,$$

$$({\pi }_X^{l}v-v,w{)}_X=0\forall w\in {\mathbb{P}}^l(X).$$

$$|v-{\pi }_T^{l}v{|}_{H^m(T)}\le C{h}_T^{s-m}|v{|}_{H^s(T)}\forall m\in \{0,\dots ,s\},$$

$$|v-{\pi }_T^{l}v{|}_{H^m(\partial T)}\le C{h}_T^{s-m-\text{nicefrac}12}|v{|}_{H^s(T)}\forall m\in \{0,\dots ,s-1\}.$$

$$|\phi {|}_{H^m(T)}:-\sum _{\mathbf{\alpha }\in {\mathbb{N}}^2,\Vert \mathbf{\alpha }{\Vert }_1=m}\Vert {\partial }^{\mathbf{\alpha }}\phi {\Vert }_{L^2(T)},|\phi {|}_{H^m(\partial T)}:-\sum _{\mathbf{\alpha }\in {\mathbb{N}}^2,\Vert \mathbf{\alpha }{\Vert }_1=m}\Vert {\partial }^{\mathbf{\alpha }}\phi {\Vert }_{L^2(\partial T)},$$

$$a_{|T}({\varpi }_T^{l}v-v,w{)}_T=0$$

$$|v-{\varpi }_T^{l}v{|}_{H^m(T)}\le C{h}_T^{s-m}|v{|}_{H^s(T)}\forall m\in \{0,\dots ,s\},$$

$$|v-{\varpi }_T^{l}v{|}_{H^m(\partial T)}\le C{h}_T^{s-m-\text{nicefrac}12}|v{|}_{H^s(T)}\forall m\in \{0,\dots ,s-1\}.$$

$${\underline{\mathrm{U}}}_T^k:-{\mathbb{P}}^k(T)\times \left({\text{bigtimes}}_{F\in {\mathcal{F}}_T}{\mathbb{P}}^k(F{)}^2\right)\times \left({\text{bigtimes}}_{F\in {\mathcal{F}}_T}{\mathbb{P}}^k(F)\right).$$

$${\underline{\mathrm{v}}}_T=(v_T,{\left({\mathbf{v}}_{\nabla ,F}\right)}_{F\in {\mathcal{F}}_T},(v_F{)}_{F\in {\mathcal{F}}_T}),$$

$${\underline{\mathrm{I}}}_T^{k}v:-\left({\pi }_T^{k}v,({\mathbf{\pi }}_F^k(\nabla v{)}_{|F}{)}_{F\in {\mathcal{F}}_T},({\pi }_F^k(v_{|F}){)}_{F\in {\mathcal{F}}_T}\right).$$

$$a_{|T}(p^{k+2}_{T}\underline{\mathrm{v}}_{T},w)=\\ -(v_{T},\mbox{\rm{div}}\,\mbox{{div}}\bm{M}_{w,T})_{T}-\sum_{F\in\mathcal{F}_{T}}\left(\bm{v}_{\bm{\nabla},F},\bm{M}_{w,T}\bm{n}_{TF}\right)_{F}+\sum_{F\in\mathcal{F}_{T}}\big{(}v_{F},\mbox{{div}}\bm{M}_{w,T}\cdot\bm{n}_{TF}\big{)}_{F},$$

$${\pi }_T^1{p}_T^{k+2}{\underline{\mathrm{v}}}_T={\pi }_T^1{v}_T.$$

$$a_{|T}(p^{k+2}_{T}\underline{\mathrm{v}}_{T},w)=\\ a_{|T}(v_{T},w)-\sum_{F\in\mathcal{F}_{T}}\left(\bm{v}_{\bm{\nabla},F}-\bm{\nabla}v_{T},\bm{M}_{w,T}\bm{n}_{TF}\right)_{F}+\sum_{F\in\mathcal{F}_{T}}\big{(}v_{F}-v_{T},\mbox{{div}}\bm{M}_{w,T}\cdot\bm{n}_{TF}\big{)}_{F}.$$

$$p_T^{k+2}\circ {\underline{\mathrm{I}}}_T^k={\varpi }_T^{k+2}.$$

$$\text{mbox}\mathrm{div}\text{mbox}div{\mathbf{M}}_{w,T}\in {\mathbb{P}}^{k-2}(T)\subset {\mathbb{P}}^k(T)$$

$$({\mathbf{M}}_{w,T}{)}_{|F}{\mathbf{n}}_{TF}\in {\mathbb{P}}^k(F{)}^2,(\text{mbox}div{\mathbf{M}}_{w,T}{)}_{|F}\cdot {\mathbf{n}}_{TF}\in {\mathbb{P}}^{k-1}(F)\subset {\mathbb{P}}^k(F).$$

$$({\pi }_T^{k}v,\text{mbox}\mathrm{div}\text{mbox}div{\mathbf{M}}_{w,T}{)}_T$$

$$({\mathbf{\pi }}_F^k(\nabla v{)}_{|F},{\mathbf{M}}_{w,T}{\mathbf{n}}_{TF}{)}_F$$

$$({\pi }_F^k{v}_{|F},\text{mbox}div{\mathbf{M}}_{w,T}\cdot {\mathbf{n}}_{TF}{)}_F$$

Full text

MSC2010: 65N30, 65N12, 74K20

Keywords: Hybrid High-Order methods, Kirchhoff–Love plates, biharmonic problems, energy projector

1 Introduction

As remarked by O. C. Zienkiewicz [39], “one of the early requirements of the Finite Element (FE) approximation was the choice of shape functions which did not lead to infinite strains on element interfaces and which therefore preserved a necessary degree of continuity”. This requirement (also called of conformity) appeared easy to satisfy for simple self-adjoint problems governed by second-order equations, where $C^{0}$-continuity at interfaces is enough. The situation is different as far as it concerns the knowledge, essential in structural engineering, of the bending of plates, whose numerical treatment has always been a goal of FE computations. Since thin plate bending in the Kirchhoff–Love approximation is governed by a fourth-order equation, $C^{1}$-continuity has to be introduced (and the continuity of both the function and of its normal derivative assured at interfaces). This was difficult to achieve and computationally expensive in the classical FE framework, see e.g. Zienkiewicz [40] for a first engineering-oriented discussion and Ciarlet [20] for a mathematically-oriented one. In order to relax such $C^{1}$-continuity condition, many non-conforming, mixed, hybrid plates elements have been studied and tested all over the last fifty years, and the literature on this subject is very broad; a minimal and by far non-exhaustive sample includes the seminal paper by Lascaux–Lesaint [35], as well as the classical works of Amara–Capatina–Chatti [2] (based on a decomposition of the constraints imposed on the bending moments by applying twice the Tartar lemma and using the symmetry of the tensor), Bathe [6], Boffi–Brezzi–Fortin [11], Brenner [13], Brenner–Scott [14], Brezzi–Fortin [16], Ciarlet[20], Comodi [23], Hughes [33], Johnson [34]; see also references therein. More recent nonconforming methods which have similarities (and differences) with the one presented here include the Hybridizable Discontinuous Galerkin method [21] of Cockburn–Dong–Guzmán and the Weak Galerkin method [37] of Lin–Wang–Ye; see also [22] concerning the passage from Discontinuous Galerkin to hybrid methods. We also cite here the mixed method of Behrens–Guzmán [8] based on a system of first-order equations, and the HHO method of [18], where the fourth-order operator in the Cahn–Hilliard equations is treated as a system of second-order operators.

A recent approach to the construction of FE spaces with $C^{1}$-regularity, on the other hand, has been developed in the context of the Virtual Element Method (VEM) [3, 9, 15]. Here, global continuity requirements are enforced by renouncing an explicit expression of the basis functions at each point, and local contributions are built using computable projections thereof (a stabilization term therefore has to be added). We refer the reader to [17] [19] for an application of $C^{1}$-conforming virtual spaces to plate-bending problems similar to the ones considered here. Nonconforming versions of the VEM have also been developed for fourth order operators, see, e.g., the very recent contributions by Antonietti–Manzini–Verani [4] (including nodal unknowns) and Zhao–Chen–Zhang [41] (with $C^{0}$-continuous virtual functions).

The Kirchhoff–Love plate bending model problem considered in this work reads

[TABLE]

where $\Omega\subset\mathbb{R}^{2}$ denotes a two-dimensional bounded and connected polygonal domain, representing the middle surface of a plate in its reference configuration, and the divergence operator is denoted by div or $\mathbf{div}$, as to whether it acts on vector- or tensor-valued fields, respectively. In (1a), $f$ represents a surface load orthogonal to the plane of the plate, and $\bm{M}$ is the moment tensor, a second-order symmetric tensor field related to the scalar unknown $u$, the deflection of the plate, by the constitutive law

[TABLE]

where $\mathbb{A}$ is a fourth-order, symmetric and uniformly elliptic tensor field, and $-\bm{\nabla}^{2}u$ is referred to as the curvature tensor. For the sake of simplicity, we assume in what follows that $\mathbb{A}$ is piecewise constant on a finite polygonal partition

[TABLE]

of $\Omega$, and that $f\in L^{2}(\Omega)$. Variational formulations are classical for problem (1). For $X\subset\overline{\Omega}$, we denote by $(\cdot,\cdot)_{X}$ the scalar product in $L^{2}(X)$, $L^{2}(X)^{2}$ or $L^{2}(X)^{2\times 2}$, depending on the context, and by $\|{\cdot}\|_{X}$ the associated norm; we omit the subscript $X$ whenever $X=\Omega$. The primal variational formulation of problem (1) reads: Find $u\in H^{2}_{0}(\Omega)$ such that

[TABLE]

Owing to the Lax–Milgram Lemma, problem (3) is well-posed (see, e.g., [38, 16]).

In this work, we propose and analyze a novel Hybrid High-Order (HHO) method for the approximation of problem (3) which sits at the far end of the spectrum of nonconforming methods, since the underlying space does not even embed $C^{0}$-continuity. HHO methods, introduced in [28] in the context of quasi-incompressible linear elasticity, are a class of new-generation discretization methods for partial differential equations with several advantageous features. The most relevant in the context of plate bending problems are:

(i) the support of arbitrary approximation orders on general polygonal meshes;

(ii) the reproduction of key continuous properties (such as, e.g., local equilibrium relations) at the discrete level;

(iii) competitive computational cost thanks to static condensation and compact stencil.

We refer the reader to [29] for an introduction covering the salient aspects of HHO methods for linear and nonlinear problems.

The HHO method for problem (3) is formulated in terms of discrete unknowns defined on mesh faces and elements (whence the term hybrid), and such unknowns are polynomials of arbitrary degree $k\geq 1$ (whence the expression high-order). The construction is conceived so that only face-based unknowns are globally coupled, whereas element-based unknowns can be eliminated by static condensation; see Remark 11 below for further details. Element-based unknowns play the role of the deflection $u$ inside elements, whereas face unknowns play the role of the traces of $u$ and of its gradient on faces. From these unknowns, a reconstruction of the deflection of degree $(k+2)$ is obtained by solving a local problem inside each element. This reconstruction is conceived so that, composed with a local reduction map, it coincides with the local energy projector and, as such, has optimal approximation properties in the space of polynomials of total degree (up to) $(k+2)$; see Theorem 12 below, whose proof hinges on the recent results of [26]. The high-order deflection reconstruction is used to formulate a local contribution, which includes a carefully tailored stabilization term. The role of the latter is to ensure coercivity with respect to a $H^{2}$-like seminorm while, at the same time, preserving the approximation properties of the local deflection reconstruction.

An extensive convergence analysis of the method is carried out. Specifically, in Theorem 12 below we prove convergence in $h^{k+1}$ (with $h$ denoting, as usual, the meshsize) in an energy-like norm and, in Theorem 16 below, a sharp estimate in $h^{k+3}$ for the $L^{2}$-norm under biharmonic regularity assumptions. The latter result highlights a salient feature of HHO methods, namely the fact that, by construction, element-based unknowns superconverge to the $L^{2}$-orthogonal projection of the exact solution on general meshes. As this happens by design (i.e., this behavior is not serendipitous), this phenomenon is henceforth referred to as supercloseness rather than superconvergence. We also show that the method satisfies locally inside each element a discrete version of the principle of virtual work with moments and shear forces obeying a law of action and reaction. The performance of the method is showcased on numerical examples, including a study of the robustness with respect to the choice of the stabilization.

The rest of the paper is organized as follows. In Section 2 we introduce the discrete setting: regularity for polygonal meshes, basic results thereon, and local projectors. A novel general result contained in this section is Theorem 1, where optimal approximation properties for the local energy projector on local polynomial spaces are studied. The proof of this theorem is given in Section 6. In Section 3 we introduce the HHO method, state the main results corresponding to Theorems 12 and 16, and provide a few numerical examples. In Section 4 we prove the local equilibrium properties of the HHO method and identify discrete equilibrated counterparts of moments and shear forces at interfaces. Section 5 collects the technical proofs of the properties of the discrete bilinear form relevant to the analysis. Conclusions and perspectives are discussed in Section 7.

2 Discrete setting

In this section we introduce some assumptions on the mesh, recall a few known results, and define two projectors on local polynomial spaces that will play a key role in the analysis of the method.

2.1 Mesh

The HHO method is built upon a polygonal mesh of the domain $\Omega$ defined prescribing a set of elements $\mathcal{T}_{h}$ and a set of faces $\mathcal{F}_{h}$.

The set of elements $\mathcal{T}_{h}$ is a finite collection of open disjoint polygons $T$ with nonzero area such that $\overline{\Omega}=\bigcup_{T\in\mathcal{T}_{h}}\overline{T}$ and $h\coloneq\max_{T\in\mathcal{T}_{h}}h_{T}$, with $h_{T}$ denoting the diameter of $T$. The set of faces $\mathcal{F}_{h}$ is a finite collection of open disjoint line segments in $\overline{\Omega}$ with nonzero length such that, for all $F\in\mathcal{F}_{h}$, (i) either there exist two distinct mesh elements $T_{1},T_{2}\in\mathcal{T}_{h}$ such that $F\subset\partial T_{1}\cap\partial T_{2}$ (and $F$ is called an interface) or (ii) there exists a mesh element $T\in\mathcal{T}_{h}$ such that $F\subset\partial T\cap\partial\Omega$ (and $F$ is called a boundary face). We assume that $\mathcal{F}_{h}$ is a partition of the mesh skeleton in the sense that $\bigcup_{T\in\mathcal{T}_{h}}\partial T=\bigcup_{F\in\mathcal{F}_{h}}\overline{F}$.

We denote by $\mathcal{F}_{h}^{\rm i}$ the set of all interfaces and by $\mathcal{F}_{h}^{\rm b}$ the set of all boundary faces, so that $\mathcal{F}_{h}=\mathcal{F}_{h}^{\rm i}\cup\mathcal{F}_{h}^{\rm b}$. The length of a face $F\in\mathcal{F}_{h}$ is denoted by $h_{F}$. For any $T\in\mathcal{T}_{h}$, $\mathcal{F}_{T}$ is the set of faces that lie on $\partial T$ (the boundary of $T$) and, for any $F\in\mathcal{F}_{T}$, $\bm{n}_{TF}$ is the unit normal to $F$ pointing out of $T$. Symmetrically, for any $F\in\mathcal{F}_{h}$, $\mathcal{T}_{F}$ is the set containing the mesh elements sharing the face $F$ (two if $F$ is an interface, one if $F$ is a boundary face).

The notion of geometric regularity for polygonal meshes is more subtle than for standard meshes. To formulate it, we assume the existence of a matching simplicial submesh, meaning that there is a conforming triangulation $\mathfrak{T}_{h}$ of the domain such that each mesh element $T\in\mathcal{T}_{h}$ is decomposed into a finite number of triangles from $\mathfrak{T}_{h}$ and each mesh face $F\in\mathcal{F}_{h}$ is decomposed into a finite number of edges from the skeleton of $\mathfrak{T}_{h}$. We denote by $\varrho$ the regularity parameter such that (i) for any triangle $S\in\mathfrak{T}_{h}$ of diameter $h_{S}$ and inradius $r_{S}$, $\varrho h_{S}\leq r_{S}$ and (ii) for any mesh element $T\in\mathcal{T}_{h}$ and any triangle $S\in\mathfrak{T}_{h}$ such that $S\subset T$, $\varrho h_{T}\leq h_{S}$. When considering refined mesh sequences, the regularity parameter should remain bounded away from zero.

In what follows, we also assume that the mesh is compliant with the data, i.e., for each mesh element $T\in\mathcal{T}_{h}$ there exist a unique polygon $\Omega_{i}\in P_{\Omega}$ (see (2)) such that $T\subset\Omega_{i}$. As a result, the material tensor field $\mathbb{A}$ is element-wise constant, and we set for the sake of brevity

[TABLE]

We also denote by $\mathcal{A}_{T}^{-}$ and $\mathcal{A}_{T}^{+}$, respectively, the smallest and largest eigenvalues of $\mathbb{A}_{T}$, regarded as an endomorphism of $\mathbb{R}^{2\times 2}_{\mathrm{sym}}$. For $l\geq 0$ we also introduce, for later use, the broken Sobolev space

[TABLE]

equipped, unless noted otherwise, with the broken norm $\|{\cdot}\|_{H^{l}(\mathcal{T}_{h})}$ defined by

[TABLE]

2.2 Basic results

We next recall a few geometric and functional inequalities, whose proofs are straightforward adaptations of the results collected in [27, Chapter 1] (where a slightly different notion of mesh faces is considered). For any mesh element $T\in\mathcal{T}_{h}$ and any face $F\in\mathcal{F}_{T}$ it holds that

[TABLE]

which expresses the fact that we are working on isotropic meshes. Moreover, the maximum number of faces of a mesh element is uniformly bounded: There is an integer $N_{\partial}\geq 3$ only depending on $\varrho$ such that

[TABLE]

Let a polynomial degree $l\geq 0$ be fixed, let $X$ be a mesh element or face, and denote by $\mathbb{P}^{l}(X)$ the space spanned by the restrictions to $X$ of two-variate polynomials of total degree at most $l\geq 0$. There exist three real numbers $C_{\rm tr}>0$, $C_{\rm tr,c}>0$, and $C_{\rm inv}$ depending on $\varrho$ and possibly on $l$, but independent of $h$, such that for any $T\in\mathcal{T}_{h}$ and $F\in\mathcal{F}_{T}$, the following discrete trace, continuous trace, and inverse inequalities hold:

[TABLE]

We also recall the following Poincaré inequality, valid for all $T\in\mathcal{T}_{h}$ and all $w\in H^{1}(T)$ such that $(w,1)_{T}=0$:

[TABLE]

where the real number $C_{\rm p}$ is independent of both $h_{T}$ and $T$, but possibly depends on $\varrho$ (for instance, $C_{\rm p}=\pi^{-1}$ for convex elements [7]).

2.3 Projectors on local polynomial spaces

Projectors on local polynomial spaces are an essential ingredient in the construction and analysis of our method. Let a polynomial degree $l\geq 0$ be fixed, and let $X$ denote a mesh element or face. The $L^{2}$-orthogonal projector $\pi_{X}^{l}\colon L^{2}(X)\to\mathbb{P}^{l}(X)$ is such that, for all $v\in L^{2}(X)$, $\pi_{X}^{l}v$ is the unique polynomial satisfying the relation

[TABLE]

The corresponding vector-valued version, denoted by $\bm{\pi}_{X}^{l}$, acts component-wise. We recall the following approximation results that are a special case of the ones proved in [25, Lemmas 3.4 and 3.6]: There exists a real number $C>0$ independent of $h$, but possibly depending on $\varrho$ and $l$, such that, for all $T\in\mathcal{T}_{h}$, all $s\in\{0,\dots,l+1\}$, and all $v\in H^{s}(T)$,

[TABLE]

Here we have set, for any $\varphi\in H^{s}(T)$,

[TABLE]

with $m$ respectively as in (11a) and (11b), $\|\bm{\alpha}\|_{1}\coloneq\alpha_{1}+\alpha_{2}$ and $\partial^{\bm{\alpha}}\coloneq\partial_{1}^{\alpha_{1}}\partial_{2}^{\alpha_{2}}$. Notice that, in the second definition, $\varphi$ and $\partial^{\bm{\alpha}}\varphi$ stand for the boundary traces of the function and of its derivatives up to order $m$, respectively.

Let a mesh element $T\in\mathcal{T}_{h}$ be fixed. For $u,v\in H^{2}(T)$, we let $a_{|T}(u,v)\coloneq\left(\mathbb{A}_{T}\bm{\nabla}^{2}u,\bm{\nabla}^{2}v\right)_{T}$ and introduce the local energy projector $\varpi_{T}^{l}:H^{2}(T)\to\mathbb{P}^{l}(T)$ such that, for any integer $l\geq 2$ and any function $v\in H^{2}(T)$,

[TABLE]

Optimal approximation properties for the local energy projector are stated in the following theorem, whose proof is given in Section 6.

Theorem 1** (Optimal approximation properties of the local energy projector).**

There is a real number $C>0$ independent of $h$, but possibly depending on $\mathbb{A}$, $\varrho$ and $l$, such that, for all $T\in\mathcal{T}_{h}$, all $s\in\{2,\ldots,l+1\}$, and all $v\in H^{s}(T)$, it holds

[TABLE]

Remark 2** (Dependence on the material tensor).**

It can be checked that the constant $C$ in the right-hand side of (13) actually depends on $\mathbb{A}_{T}$ only through the square root of the ratio between $\mathcal{A}_{T}^{+}$ and $\mathcal{A}_{T}^{-}$.

3 The Hybrid High-Order method

In this section, we present the construction underlying the HHO method, state the discrete problem, and discuss the main results. Henceforth, we fix once and for all a polynomial degree $k\geq 1$.

3.1 Local discrete unknowns and interpolation

Let a mesh element $T\in\mathcal{T}_{h}$ be fixed. The local space of discrete unknowns is defined as the set

[TABLE]

For a general collection of discrete unknowns $\underline{\mathrm{v}}_{T}\in\underline{\mathrm{U}}_{T}^{k}$, we use the standard underlined HHO notation

[TABLE]

where $v_{T}$ contains the element-based discrete unknowns, $\bm{v}_{\bm{\nabla},F}$ the discrete unknowns related to the trace of the gradient on the face $F$, and $v_{F}$ the discrete unknowns related to the trace on $F$.

The local interpolation operator $\underline{\mathrm{I}}_{T}^{k}\colon H^{2}(T)\to\underline{\mathrm{U}}_{T}^{k}$ is such that, for all $v\in H^{2}(T)$,

[TABLE]

Since the boundary of $T$ is piecewise smooth (see, e.g., [38]), the trace theorem ensures that the restrictions $v_{|F}$ and $(\bm{\nabla}v)_{|F}$ of $v$ appearing in (15) are both well-defined.

3.2 Local deflection reconstruction

Let again a mesh element $T\in\mathcal{T}_{h}$ and a polynomial degree $k\geq 1$ be fixed. We introduce the local deflection reconstruction operator $p^{k+2}_{T}\colon\underline{\mathrm{U}}_{T}^{k}\to\mathbb{P}^{k+2}(T)$ such that, for all $\underline{\mathrm{v}}_{T}\in\underline{\mathrm{U}}_{T}^{k}$, $p^{k+2}_{T}\underline{\mathrm{v}}_{T}\in\mathbb{P}^{k+2}(T)$ satisfies for all $w\in\mathbb{P}^{k+2}(T)$

[TABLE]

where $\bm{M}_{w,T}\coloneq-\mathbb{A}_{T}\bm{\nabla}^{2}w$. Here, the notation $\bm{M}_{w,T}$ is used to emphasize the fact that $\bm{M}_{w,T}$ is a moment tensor of virtual nature (with space of virtual deflections equal to $\mathbb{P}^{k+2}(T)$) unlike tensor $\bm{M}$ appearing in bilinear form $a$ introduced in (3). The right-hand side of (16) is conceived so as to resemble an integration by parts formula where the roles of the function represented by $\underline{\mathrm{v}}_{T}$ and of its gradient are played by element discrete unknowns inside volumetric integrals and by face-based discrete unknowns on boundary integrals.

Since ${\mbox{\rm ker}\,}\bm{\nabla}^{2}=\mathbb{P}^{1}(T)$, the compatibility condition for problem (16) requires that the linear form on the right-hand side vanish on the elements of $\mathbb{P}^{1}(T)$; since $\bm{M}_{w,T}=\bm{0}$ for all $w\in\mathbb{P}^{1}(T)$, this condition is satisfied. The solution of (16) is not unique: if $p^{k+2}_{T}\underline{\mathrm{v}}_{T}\in\mathbb{P}^{k+2}(T)$ is a solution, $p^{k+2}_{T}\underline{\mathrm{v}}_{T}+z_{T}$ for any $z_{T}\in\mathbb{P}^{1}(T)$ also is. To ensure uniqueness, we add the closure condition

[TABLE]

Notice, in passing, that element discrete unknowns do not contribute to the right-hand side of (16) for $k=1$, and they only appear in the closure condition (17).

For further use, we also observe that, since $v_{T}$ is smooth, performing an integration by parts on the first term in the right-hand side of (16) and using the symmetry of $\mathbb{A}_{T}$ leads to the following reformulation, which points out the non-conformity of the method:

[TABLE]

The definition of $p^{k+2}_{T}$ is justified by the following proposition, which establishes a link with the local energy projector defined by (12).

Proposition 3** (Link with the local energy projector).**

It holds

[TABLE]

Proof.

We write (16) for $\underline{\mathrm{v}}_{T}=\underline{\mathrm{I}}_{T}^{k}v$ (cf. (15) for the definition of the local interpolator). Since $w\in\mathbb{P}^{k+2}(T)$ and $\mathbb{A}_{T}$ is a constant tensor, we infer that

[TABLE]

and, for all $F\in\mathcal{F}_{T}$,

[TABLE]

Consequently, recalling the definition (10) of $\pi^{k}_{T}$, $\bm{\pi}_{F}^{k}$, and $\pi_{F}^{k}$, we have

[TABLE]

Plugging the above identities into the right-hand side of (16), performing an integration by parts, and using the symmetry of $\mathbb{A}_{T}$, we arrive at the following orthogonality condition:

[TABLE]

Comparing (20) and (17) with the definition (12) of $\varpi_{T}^{k+2}$ concludes the proof. ∎

Remark 4** (Approximation properties for $p^{k+2}_{T}\circ\underline{\mathrm{I}}_{T}^{k}$).**

The above result implies that $p^{k+2}_{T}\circ\underline{\mathrm{I}}_{T}^{k}=\varpi_{T}^{k+2}$ has optimal approximation properties in $\mathbb{P}^{k+2}(T)$, in the sense made precise by Theorem 1.

3.3 Local contribution

We introduce the local discrete bilinear form $\mathrm{a}_{T}(\cdot,\cdot)$ on $\underline{\mathrm{U}}_{T}^{k}\times\underline{\mathrm{U}}_{T}^{k}$ given by

[TABLE]

Here, the first contribution is the usual Galerkin term responsible for consistency. The second contribution, in charge of stability, penalizes high-order differences between the reconstruction and the unknowns and is such that, for all $(\underline{\mathrm{u}}_{T},\underline{\mathrm{v}}_{T})\in\underline{\mathrm{U}}_{T}^{k}\times\underline{\mathrm{U}}_{T}^{k}$,

[TABLE]

Remark 5** (Stabilization).**

Other expressions are possible for the stabilization term, and the specific choice can affect the accuracy of the results. In particular, the discussion below remains true if we replace (21) by

[TABLE]

with $\eta>0$ denoting a user-dependent parameter independent of $h$. In practice, it is important that the numerical results be only marginally affected by the specific choice of the stabilization. We refer the reader to Section 3.7 below for a numerical study of the robustness of the method with respect to $\eta$.

The following proposition states a consistency result for the stabilization bilinear form (22).

Proposition 6** (Consistency of $\mathrm{s}_{T}$).**

There is a real number $C>0$ independent of $h$, but possibly depending on $\mathbb{A}$, $\varrho$ and $k$, such that, for all $v\in H^{k+3}(T)$,

[TABLE]

Proof.

We have

[TABLE]

where, recalling Proposition 3 and using the linearity and idempotency of projectors,

[TABLE]

By the boundedness of $L^{2}$-projectors, along with the approximation properties (13a)–(13b) of $\varpi^{k+2}_{T}$ with $s=k+3$ and, respectively, $m=0$ for $\mathfrak{T}_{1}$, $m=1$ for $\mathfrak{T}_{2}$, and again $m=0$ for $\mathfrak{T}_{3}$, the conclusion follows. ∎

We equip the space $\underline{\mathrm{U}}^{k}_{T}$ with the following local discrete seminorm:

[TABLE]

The following result shows that the bilinear form $\mathrm{a}_{T}$ induces on $\underline{\mathrm{U}}_{T}^{k}$ a seminorm $\|{\cdot}\|_{\mathrm{a},T}$ uniformly equivalent to $\|{\cdot}\|_{\mathbb{A},T}$.

Lemma 7** (Local coercivity and boundedness).**

There is a real number $C>0$ independent of $h$, but possibly depending on $\mathbb{A}$, $\varrho$ and $k$, such that, for all $T\in\mathcal{T}_{h}$, the following inequalities hold (expressing, respectively, the coercivity and boundedness of $\mathrm{a}_{T}$):

[TABLE]

Proof.

See Section 5.1. ∎

Remark 8** (Polynomial degree).**

The assumption $k\geq 1$ is essential in the proof of the above result. For this reason, the steps in which this hypothesis is used are pointed out accordingly.

3.4 Global space, interpolation, and norm

We define the following global space of discrete unknowns:

[TABLE]

Note that interface unknowns in $\underline{\mathrm{U}}_{h}^{k}$ are single-valued, i.e., their values match from one element to the adjacent one. For a collection of discrete unknowns in $\underline{\mathrm{U}}_{h}^{k}$, we use the notation

[TABLE]

and we denote by $\underline{\mathrm{v}}_{T}=\left(v_{T},(\bm{v}_{\bm{\nabla},F})_{F\in\mathcal{F}_{T}},(v_{F})_{F\in\mathcal{F}_{T}}\right)\in\underline{\mathrm{U}}^{k}_{T}$ its restriction to a mesh element $T\in\mathcal{T}_{h}$. We also denote by $v_{h}$ (no underline) the broken polynomial function on $\mathcal{T}_{h}$ such that

[TABLE]

We define the global interpolator $\underline{\mathrm{I}}_{h}^{k}\colon H^{2}(\Omega)\to\underline{\mathrm{U}}_{h}^{k}$ such that, for all $v\in H^{2}(\Omega)$,

[TABLE]

The space $\underline{\mathrm{U}}_{h}^{k}$ is equipped with the following seminorm (cf. (25) for the definition of $\|{\cdot}\|_{\mathbb{A},T}$):

[TABLE]

We notice that the couple of boundary conditions (1b)–(1c) is equivalent to the couple $u=0$ on $\partial\Omega$ and $\bm{\nabla}u=\bm{0}$ on $\partial\Omega$. Indeed, the fact that $u$ vanishes on $\partial\Omega$ implies its tangential derivative to vanish on $\partial\Omega$ as well. Accounting for this remark, we introduce the following subspace that incorporates the latter couple of boundary conditions in a strong manner:

[TABLE]

It is a simple matter to check that the image of the restriction of $\underline{\mathrm{I}}_{h}^{k}$ to $H^{2}_{0}(\Omega)$ is contained in $\underline{\mathrm{U}}_{h,0}^{k}$.

Proposition 9** (Norm $\|\underline{\mathrm{v}}_{h}\|_{\mathbb{A},h}$).**

The mapping $\underline{\mathrm{U}}_{h,0}^{k}\ni\underline{\mathrm{v}}_{h}\mapsto\|\underline{\mathrm{v}}_{h}\|_{\mathbb{A},h}\in\mathbb{R}$ defines a norm on $\underline{\mathrm{U}}_{h,0}^{k}$.

Proof.

The seminorm property is trivial. It then suffices to show that $\|\underline{\mathrm{v}}_{h}\|_{\mathbb{A},h}=0\Longrightarrow\underline{\mathrm{v}}_{h}=\underline{0}\in\underline{\mathrm{U}}^{k}_{h,0}$. Clearly, $\|\underline{\mathrm{v}}_{h}\|_{\mathbb{A},h}=0$ implies $\bm{\nabla}^{2}v_{T}\equiv\bm{0}$ for all $T\in\mathcal{T}_{h}$ and $\bm{v}_{\bm{\nabla},F}-\bm{\nabla}v_{T}\equiv\bm{0}$ and $v_{F}-v_{T}\equiv 0$ for all $F\in\mathcal{F}_{h}$. By definition (30), we have $\bm{v}_{\bm{\nabla},F}\equiv\bm{0}$ and $v_{F}=0$ for all $F\in\mathcal{F}^{b}_{h}$; thus, for any $T\in\mathcal{T}_{h}$, if $\mathcal{F}_{T}\cap\mathcal{F}^{\mathrm{b}}_{h}\neq\emptyset$ then there exists $F\in\mathcal{F}^{\mathrm{b}}_{h}$ such that $\bm{\nabla}v_{T}\equiv\bm{0}$ and $v_{T}\equiv 0$ on $F$. Since $\bm{\nabla}^{2}v_{T}\equiv\bm{0}$ in $T$, these facts imply that $v_{T}\equiv 0$ in $T$, which in turn implies that $v_{F}\equiv 0$ and $\bm{v}_{\bm{\nabla},F}\equiv\bm{0}$ for all $F\in\mathcal{F}_{T}$. Repeating this argument for inner layers of elements yields the assertion. ∎

3.5 Discrete problem

The discrete problem is formulated as follows: Find $\underline{\mathrm{u}}_{h}\in\underline{\mathrm{U}}^{k}_{h,0}$ such that

[TABLE]

with global bilinear form $\mathrm{a}_{h}$ on $\underline{\mathrm{U}}_{h}^{k}\times\underline{\mathrm{U}}_{h}^{k}$ obtained by element-by-element assembly setting

[TABLE]

The following lemma summarizes the properties of the global bilinear form $\mathrm{a}_{h}$.

Lemma 10** (Properties of $\mathrm{a}_{h}$).**

The bilinear form $\mathrm{a}_{h}$ defined by (32) has the following properties:

1. (i)

Coercivity and boundedness.* There is a real number $C>0$ independent of $h$, but possibly depending on $\mathbb{A}$, $\varrho$ and $k$, such that*

[TABLE] 2. (ii)

Consistency.* There is a real number $C>0$ independent of $h$, but possibly depending on $\mathbb{A}$, $\varrho$ and $k$, such that, for all $v\in H_{0}^{2}(\Omega)\cap H^{4}(\Omega)\cap H^{k+3}(\mathcal{T}_{h})$, it holds that*

[TABLE]

Proof.

See Section 5.1. ∎

As a consequence of the first inequality in (33), the discrete problem (31) admits a unique solution. This solution minimizes the following discrete energy:

[TABLE]

The discrete energy will play an important role in numerical experiments (cf. Section 3.7 below).

Remark 11** (Implementation).**

Proceeding as in standard FE methods, to write an algebraic version of problem (31) we associate to each mesh element or face a set of degrees of freedom (DOFs) that form a basis for the dual space of the local polynomial space supported by it. Let a basis $\mathcal{B}_{h}$ for the space $\underline{\mathrm{U}}_{h,0}^{k}$ be fixed such that every basis function is supported by only one mesh element or face. To fix the ideas, we take as DOFs the coefficients of the expansion of a HHO function $\underline{\mathrm{v}}_{h}\in\underline{\mathrm{U}}_{h,0}^{k}$ in $\mathcal{B}_{h}$, and we collect them in the vector $\mathsf{V}_{h}$ partitioned as

[TABLE]

where the subvector $\mathsf{V}_{\mathcal{T}_{h}}$ collects the coefficients associated with element-based DOFs, while the remaining coefficients (associated to face-based DOFs) are collected in $\mathsf{V}_{\mathcal{F}_{h}}$. Denote by $\mathsf{A}_{h}$ the matrix representation of the bilinear form $\mathrm{a}_{h}$ and by $\mathsf{B}_{h}$ the vector representation of the linear form $\underline{\mathrm{v}}_{h}\mapsto(f,v_{h})$, both partitioned in a similar way. The algebraic problem corresponding to (31) reads

[TABLE]

The submatrix $\mathsf{A}_{\mathcal{T}_{h}\mathcal{T}_{h}}$ is block-diagonal and symmetric positive definite, and is therefore inexpensive to invert. In the practical implementation, this remark can be exploited by solving the linear system (36) in two steps:

[TABLE]

1. (i)

First, element-based coefficients in $\mathsf{U}_{\mathcal{T}_{h}}$ are expressed in terms of $\mathsf{B}_{\mathcal{T}_{h}}$ and $\mathsf{U}_{\mathcal{F}_{h}}$ by the inexpensive solution of the first block equation:

[TABLE]

This step is referred to as static condensation in the FE literature; 2. (ii)

Second, face-based coefficients in $\mathsf{U}_{\mathcal{F}_{h}}$ are obtained solving the following global problem involving quantities attached to the mesh skeleton:

[TABLE]

This computationally more intensive step requires to invert the symmetric matrix $\tilde{\mathsf{A}}_{\mathcal{F}_{h}\mathcal{F}_{h}}$, whose stencil involves neighbours through faces, and which has size $N_{\rm dof}\times N_{\rm dof}$ with $N_{\rm dof}=2\operatorname{card}(\mathcal{F}_{h}^{\rm i}){k+d-1\choose k}$. Observing that $\tilde{\mathsf{A}}_{\mathcal{F}_{h}\mathcal{F}_{h}}$ is in fact the Schur complement of $\mathsf{A}_{\mathcal{T}_{h}\mathcal{T}_{h}}$ in $\mathsf{A}_{h}$, and since $\mathsf{A}_{h}$ is symmetric and both $\mathsf{A}_{h}$ and $\mathsf{A}_{\mathcal{T}_{h}\mathcal{T}_{h}}$ are positive definite, a classical result in linear algebra yields that also $\tilde{\mathsf{A}}_{\mathcal{F}_{h}\mathcal{F}_{h}}$ is positive definite (see, e.g., **[32]**).

3.6 Main results

We next present the main results of the analysis, namely error estimates in an energy-like norm, in a jump-seminorm, and in the $L^{2}$-norm. Inside the proofs of this section, we often abridge as $a\lesssim b$ the inequality $a\leq Cb$ with $C>0$ independent of $h$, but possibly depending on $\mathbb{A}$, $\varrho$, and $k$.

3.6.1 Energy error estimate

We introduce the global deflection reconstruction operator $p^{k+2}_{h}:\underline{\mathrm{U}}_{h}^{k}\to L^{2}(\Omega)$ such that, for all $\underline{\mathrm{v}}_{h}\in\underline{\mathrm{U}}_{h}^{k}$,

[TABLE]

We also define the stabilization seminorm $|{\cdot}|_{\mathrm{s},h}$ on $\underline{\mathrm{U}}_{h}^{k}$ setting, for all $\underline{\mathrm{v}}_{h}\in\underline{\mathrm{U}}_{h}^{k}$,

[TABLE]

Theorem 12** (Energy error estimate).**

Let $u\in H^{2}_{0}(\Omega)$ and $\underline{\mathrm{u}}_{h}\in\underline{\mathrm{U}}^{k}_{h,0}$ denote the unique solutions to the continuous (3) and discrete (31) problems, respectively. Assume the additional regularity $u\in H^{4}(\Omega)\cap H^{k+3}(\mathcal{T}_{h})$. Then, it holds that

[TABLE]

where $\bm{\nabla}_{h}$ denotes the usual broken gradient operator on $\mathcal{T}_{h}$ and the real number $C>0$ is independent of $h$ (but possibly depends on $\mathbb{A}$, $\varrho$, and $k$).

Remark 13** (Regularity of the solution).**

Concerning the regularity assumptions on $u$, we mention as an example that, for $k=1$, the regularity $u\in H^{4}(\Omega)$ is satisfied by the solution of the biharmonic problem with Dirichlet boundary conditions (obtained taking $\mathbb{A}=\mathbb{I}$ in (1)) posed on a three-dimensional cubic domain, provided the load $f$ is square-integrable (see, e.g., Maz’ya [36, Chapter 4]). In two dimensions, under the weaker assumption that $f\in H^{-1}(\Omega)$, it holds that $u\in H^{3}(\Omega)$ provided $\Omega$ is convex (see, e.g., Grisvard [30, Chapter 3]). In general, a regularity assumption on the exact solution is actually the consequence of a compatibility condition between the datum regularity and the domain geometry. When $f\in L^{2}(\Omega)$ in two dimensions, one can have $u\in H^{4}(\Omega)$ under the condition of Kondratiev on the opening of each corner (see, e.g., Blum–Rannacher [10], Grisvard [31]). As a further reference on the regularity for the solution of fourth-order elliptic problems, we also refer the reader to Dauge [24, Chapter 4]. To close this remark, we emphasize that, since $u$ needs only be locally regular inside each element, the presence of corner singularities and layers can be accounted for by a judicious choice of $h$ and, possibly, of $k$.

Proof of Theorem 12.

Let, for the sake of brevity, $\widehat{\underline{\mathrm{u}}}_{h}\coloneq\underline{\mathrm{I}}^{k}_{h}u$. We start by proving that

[TABLE]

with norm $\|{\cdot}\|_{\mathrm{a},h}$ defined by (33). Using the linearity of $\mathrm{a}_{h}$ in its first argument together with the discrete problem (31), and recalling that $\mbox{\rm{div}}\,\mbox{{div}}\mathbb{A}\bm{\nabla}^{2}u=f$ a.e. in $\Omega$, we have, for all $\underline{\mathrm{v}}_{h}\in\underline{\mathrm{U}}^{k}_{h,0}$,

[TABLE]

Thus, choosing $\underline{\mathrm{v}}_{h}=\underline{\mathrm{u}}_{h}-\widehat{\underline{\mathrm{u}}}_{h}$ and using the consistency (34) of $\mathrm{a}_{h}$ to bound the supremum in the right-hand side, the basic estimate (39) follows.

Let us now prove (38). Using the triangle inequality, we infer that

[TABLE]

where we have used the discrete Cauchy–Schwarz inequality together with the definition (33) of the $\|{\cdot}\|_{\mathrm{a},h}$-norm in the last line. The conclusion follows using (39) to estimate the first term in the right-hand side, the optimal approximation properties (13a) of $p^{k+2}_{T}\widehat{\underline{\mathrm{u}}}_{T}=\varpi_{T}^{k+2}u$ with $s=k+3$ and $m=2$ for all $T\in\mathcal{T}_{h}$ to estimate the second term, and the consistency (24) of $\mathrm{s}_{T}$ for all $T\in\mathcal{T}_{h}$ to estimate the third term. ∎

Remark 14** (Convergence of face unknowns).**

Combining the norm equivalence (33) with (39), we readily infer that

[TABLE]

which shows, in particular, that the face variables converge in an energy-like norm to the corresponding projections of the exact solution and its normal derivative. This is in itself a supercloseness result for the face variables, since, replacing $\bm{\pi}_{F}^{k}(\bm{\nabla}u)$ by $\bm{\nabla}u$ and $\pi_{F}^{k}u$ by $u$ in the left-hand side of the above inequality, one would only obtain a suboptimal estimate in $h^{k-1}$ (which would only converge if $k\geq 2$). An optimal error estimate in $h^{k+1}$ for the trace of $u$ and its gradient can be recovered using the deflection reconstruction instead of the face variables:

[TABLE]

Remark 15** (Convergence of the jumps).**

From the estimate of Theorem 12, one can prove that the jumps of $p^{k+2}_{h}\underline{\mathrm{u}}_{h}$ and of its gradient converge to zero with optimal rate. To this end, define on $H^{2}(\mathcal{T}_{h})$ (cf. definition (4)) the following jump seminorm:

[TABLE]

where $[{\cdot}]_{F}$ is the usual jump operator if $F$ is an interface (the sign is irrelevant), whereas $[\varphi]_{F}\coloneq\varphi_{|F}$ if $F$ is a boundary face, and $\mathcal{A}_{F}\coloneq\min_{T\in\mathcal{T}_{F}}\mathcal{A}_{T}^{+}$. Then, observing that $|p^{k+2}_{h}\underline{\mathrm{u}}_{h}|_{\mathrm{J},h}^{2}\leq 2\varrho^{-6}|\underline{\mathrm{u}}_{h}|^{2}_{\mathrm{s},h}$ as a consequence of the triangle inequality together with (6), it is inferred from (38) that

[TABLE]

with real number $C>0$ independent of $h$, but possibly depending on $\mathbb{A}$, $\varrho$, and $k$.

3.6.2 $\bm{L}^{2}$-error estimate

A sharp $L^{2}$-norm error estimate can also be inferred assuming biharmonic regularity, in the following form: For all $q\in L^{2}(\Omega)$, the unique solution $z\in H^{2}_{0}(\Omega)$ to

[TABLE]

satisfies the a priori estimate (see, e.g., [10])

[TABLE]

with $C_{\mathrm{bihar}}>0$ only depending on $\Omega$ and on $\mathbb{A}$.

Theorem 16** ($L^{2}$-error estimate).**

Let $u\in H^{2}_{0}(\Omega)$ and $\underline{\mathrm{u}}_{h}\in\underline{\mathrm{U}}^{k}_{h,0}$ denote the unique solutions to the continuous (3) and discrete (31) problems, respectively. Assume $u\in H^{4}(\Omega)\cap{H^{k+3}(\mathcal{T}_{h})}$, biharmonic regularity, and $f\in H^{k+1}(\mathcal{T}_{h})$. Then, there exists a real number $C>0$ depending on $\mathbb{A}$, $\varrho$, and $k$, but independent of $h$, such that

[TABLE]

Proof.

Let, for the sake of brevity, $\widehat{\underline{\mathrm{u}}}_{h}\coloneq\underline{\mathrm{I}}^{k}_{h}u$. By the triangle inequality, we have that

[TABLE]

By the approximation properties (13a) of $p^{k+2}_{T}\circ\underline{\mathrm{I}}^{k}_{T}=\varpi^{k+2}_{T}$ (cf. Remark 4) with $s=k+3$ and $m=0$, we immediately have that

[TABLE]

For the second term, on the other hand, we observe that

[TABLE]

where we have used the triangle inequality and the approximation properties of $\pi_{T}^{1}$ for $s=2$ and $m=0$, as well as the closure condition (17) to pass to the second line, and the definition of the $\|{\cdot}\|_{\mathrm{a},h}$-norm to conclude. Using (39) and Lemma 17 below to bound the first and second addend in the right-hand side, respectively, the conclusion follows. ∎

The following lemma, used in the proof of Theorem 16 above, shows that element-based discrete unknowns behave “almost” like the $L^{2}$-orthogonal projection of the exact solution on the space of broken polynomials of total degree at most $k$ on $\mathcal{T}_{h}$.

Lemma 17** (Supercloseness of element discrete unknowns).**

Under the assumptions and notations of Theorem 16, it holds that

[TABLE]

where $\widehat{u}_{h}$ and $u_{h}$ are the broken polynomial functions of total degree at most $k$ such that $\widehat{u}_{h|T}\coloneq\widehat{u}_{T}=\pi^{k}_{T}u$ and $u_{h|T}\coloneq u_{T}$ for any mesh element $T\in\mathcal{T}_{h}$.

Proof.

Set, for the sake of brevity, $\underline{\mathrm{e}}_{h}\coloneq\widehat{\underline{\mathrm{u}}}_{h}-\underline{\mathrm{u}}_{h}$ and $e_{h}\coloneq\widehat{u}_{h}-u_{h}$. Let $z$ solve (41) with $q=e_{h}$ and set $\widehat{\underline{\mathrm{z}}}_{h}\coloneq\underline{\mathrm{I}}^{k}_{h}z$. Integrating by parts, using the linearity of $\mathrm{a}_{h}$ in its first argument, as well as the continuity of moments and shear forces at interfaces, and letting $\widecheck{z}_{T}\coloneq\varpi^{k+2}_{T}(z_{|T})$ for all $T\in\mathcal{T}_{h}$, we have that $\|e_{h}\|^{2}=\mathfrak{T}_{1}+\mathfrak{T}_{2}$, with

[TABLE]

where $\pi^{k}_{h}$ is such that $(\pi^{k}_{h}v)_{|T}=\pi^{k}_{T}(v_{|T})$ for all $T\in\mathcal{T}_{h}$ and all $v\in H^{2}(\Omega)$. The Cauchy–Schwarz inequality then yields

[TABLE]

The approximation properties (13) of $\varpi^{k+2}_{T}$ with $s=4$, the consistency property (24) of the stabilization bilinear form, and the stability of $\mathrm{a}_{h}$ together with the energy error estimate (39) allow to conclude that

[TABLE]

where in the last estimate we have used the biharmonic regularity hypothesis. Turning to $\mathfrak{T}_{2}$, using the fact that $(f,\pi^{k}_{T}z)_{T}=(\pi^{k}_{T}f,z)_{T}$ and exploiting the orthogonality property (20), we have

[TABLE]

We have that $|\mathfrak{T}_{2,1}|\lesssim h^{k+3}\|u\|_{{H^{k+3}(\mathcal{T}_{h})}}\|e_{h}\|$ by the Cauchy–Schwarz inequality, the approximation properties of $\varpi^{k+2}_{T}$, and biharmonic regularity. An analogous bound can be obtained for $|\mathfrak{T}_{2,2}|$. Finally, we observe that $\mathfrak{T}_{2,3}=(f-\pi^{k}_{h}f,z-\pi^{k}_{h}z)$ by the definition (10) of the $L^{2}$-orthogonal projector. Using the approximation properties (11a) of $\pi^{k}_{T}$ with $l=k$, $m=0$, and $s=k+1$ for the first factor, $s=2$ for the second, we obtain

[TABLE]

This concludes the proof. ∎

3.7 Numerical examples

In this section we solve problem (1) for $\mathbb{A}=\mathbb{I}$ (i.e., the biharmonic equation) in the unit square and, with a view towards testing the convergence of the method in the case of more complex geometries, in a L-shaped domain as well.

3.7.1 Unit square

In this first case, the domain under consideration is $\Omega=(0,1)\times(0,1)$. The right-hand side $f$ is set in agreement with the exact solution

[TABLE]

on three different meshes: triangular, cartesian and hexagonal (cf. Fig. 1). Figures 2 and 3 show convergence results in the energy norm and in the $L^{2}$-norm, respectively, for different meshes and polynomial degrees, up to three. We consider $\|\widehat{\underline{\mathrm{u}}}_{h}-\underline{\mathrm{u}}_{h}\|_{\mathrm{a},h}$ and $\|\pi^{k}_{h}u-u_{h}\|$ as measures of the error in the energy norm and in the $L^{2}$-norm, respectively. Since biharmonic regularity (42) is satisfied (the domain is convex and the exact solution is of class $C^{\infty}$), the numerical results show asymptotic convergence rates that match those predicted by the theory, i.e. estimates (39) and (44), in all of the three cases.

Also, we check the numerical convergence of the discrete energy (35) with four uniformly refined triangular meshes, and a polynomial degree $k$ ranging from 1 to 4. As Table 1 shows, only three refinements are necessary when $k\in\{2,3,4\}$ to achieve a five-significant-digit precision for the limit of the discrete energy.

We finally test the robustness of the variant of the HHO method based on the local bilinear form (23) with respect to the user-dependent parameter $\eta$. In Figures 4 and 5 we plot, respectively, the energy- and $L^{2}$-norms of the error when $\eta$ varies from $10^{-3}$ to $10^{3}$ on fixed meshes corresponding to the third refinement level of the ones in Figure 1. From these plots, the robustness of the method can be appreciated, as the energy error spans only two orders of magnitude and the $L^{2}$-error spans four orders of magnitude, while the user-dependent parameter $\eta$ spans six orders of magnutide.

3.7.2 L-shaped domain

We now consider the domain $\Omega=\left((0,1)\times(0,1)\right)\setminus\left(\left({1}/{2},1\right)\times\left({1}/{2},1\right)\right)$, and a uniform load $f\equiv 1$. Figure 6a shows the numerical solution obtained for $k=3$ on five nested, uniformly refined triangular meshes. Since a closed-form solution is not available in this case, we check the numerical convergence of the discrete energy on the above-mentioned meshes, again for a polynomial degree $k$ ranging from $1$ to $4$. As Table 2 shows, this energy converges numerically towards a value given by -2.80e-05 to two significant digits. This allows to conclude that the method converges even in situations where such singular geometries are considered. As expected, since biharmonic regularity is not satisfied in this case (because of the domain geometry), convergence is slower than in Table 1, and five mesh refinements are required to achieve a two-significant-digit precision for the limit. For further details, we refer the reader to Section 7.1.

4 Local principle of virtual work and laws of action-reaction

Let a mesh element $T\in\mathcal{T}_{h}$ be fixed. At the continuous level, the deflection field $u$ satisfies, for all $v\in\mathbb{P}^{k}(T)$,

[TABLE]

For an interface $F\in\mathcal{F}_{T_{1}}\cap\mathcal{F}_{T_{2}}$, with $T_{1}$, $T_{2}$ distinct elements of $\mathcal{T}_{h}$, since $\bm{n}_{T_{2}F}=-\bm{n}_{T_{1}F}$, both moments and shear forces obey the following laws of action-reaction:

[TABLE]

The denomination for equations (46b) emphasizes the fact that the moment (resp., shear force) exerted on element $T_{1}$ by element $T_{2}$ through the common interface $F$ is the opposite of the moment (resp., shear force) exerted on $T_{2}$ by $T_{1}$ through $F$.

We next show that the solution to discrete problem (31) satisfies discrete counterparts of (46a) and (46b). This requires a reformulation of the stabilization contribution in terms of the differences between face-based and element-based discrete unknowns. Define the space

[TABLE]

and the boundary difference operator $\underline{\mathrm{\delta}}_{\partial T}^{k}\colon\underline{\mathrm{U}}_{T}^{k}\to\underline{\mathrm{D}}_{\partial T}^{k}$ such that, for all $\underline{\mathrm{v}}_{T}\in\underline{\mathrm{U}}_{T}^{k}$,

[TABLE]

Proposition 18** (Boundary difference reformulation of $\mathrm{s}_{T}$).**

The local stabilization bilinear form $\mathrm{s}_{T}$ defined by (22) can be rewritten, for all $\underline{\mathrm{u}}_{T},\underline{\mathrm{v}}_{T}\in\underline{\mathrm{U}}_{T}^{k}$,

[TABLE]

Proof.

As a consequence of (19), for all $v_{T}\in\mathbb{P}^{k}(T)\subset\mathbb{P}^{k+2}(T)$ it holds

[TABLE]

where we have used the fact that, as a projector, $\varpi^{k+2}_{T}$ preserves polynomials up to degree $(k+2)$. Now, using (48) and the linearity of $p^{k+2}_{T}$, we have

[TABLE]

Also, for all $F\in\mathcal{F}_{T}$, it holds

[TABLE]

and, analogously,

[TABLE]

Using (49), (50), and (51) respectively in the first, second, and third term in the right-hand side of (22), the conclusion follows. ∎

Define now the residual operator

[TABLE]

such that, for all $\underline{\mathrm{v}}_{T}\in\underline{\mathrm{U}}^{k}_{T}$ and all $\underline{\mathrm{\alpha}}_{\partial T}\equiv((\bm{\alpha}_{\bm{\nabla},F})_{F\in\mathcal{F}_{T}},(\alpha_{F})_{F\in\mathcal{F}_{T}})\in\underline{\mathrm{D}}^{k}_{\partial T}$,

[TABLE]

Problem (52) is well-posed as a consequence of the Riesz representation theorem for the $L^{2}$-like product in the left-hand side.

Lemma 19** (Local principle of virtual work and laws of action-reaction).**

Denote by $\underline{\mathrm{u}}_{h}\in\underline{\mathrm{U}}^{k}_{h,0}$ the unique solution to (31) and, for all $T\in\mathcal{T}_{h}$ and all $F\in\mathcal{F}_{T}$, define the discrete moment and shear force

[TABLE]

Then, the following discrete counterparts of (46a) and (46b) hold, respectively: For any mesh element $T\in\mathcal{T}_{h}$,

[TABLE]

Proof.

Recalling the definition (21) of $\mathrm{a}_{T}$, and using the reformulation (47) of ${\mathrm{s}}_{T}$ together with the definition (52) of the residual operator, it is inferred from the discrete problem (31) that, for all $\underline{\mathrm{v}}_{h}\in\underline{\mathrm{U}}_{h,0}^{k}$, it holds

[TABLE]

Using the definition (18) of $p^{k+2}_{T}\underline{\mathrm{v}}_{T}$ with $w=p^{k+2}_{T}\underline{\mathrm{u}}_{T}$ for the first term, and recalling (52) and (53), we can rewrite (55) as

[TABLE]

Thus, for a given mesh element $T\in\mathcal{T}_{h}$, choosing in (56) $\underline{\mathrm{v}}_{h}$ such that $v_{T}$ spans $\mathbb{P}^{k}(T)$, $v_{T^{\prime}}\equiv 0$ for all $T^{\prime}\in\mathcal{T}_{h}\setminus\{T\}$, $\bm{v}_{\bm{\nabla},F}\equiv\bm{0}$ and $v_{F}\equiv 0$ for all $F\in\mathcal{F}_{h}$ immediately yields (54a). Next, for a given interface $F\in\mathcal{F}_{T_{1}}\cap\mathcal{F}_{T_{2}}$, choosing in (56) $\underline{\mathrm{v}}_{h}$ such that $v_{T}\equiv 0$ for all $T\in\mathcal{T}_{h}$, $\bm{v}_{\bm{\nabla},F^{\prime}}\equiv\bm{0}$ for all $F^{\prime}\in\mathcal{F}_{h}\setminus\{F\}$, $v_{F}\equiv 0$ for all $F\in\mathcal{F}_{h}$, and letting $\bm{v}_{\bm{\nabla},F}$ span $\mathbb{P}^{k}(F)^{2}$ yields the first equation in (54b). Similarly, choosing in (56) $\underline{\mathrm{v}}_{h}$ such that $v_{T}\equiv 0$ for all $T\in\mathcal{T}_{h}$, $\bm{v}_{\bm{\nabla},F}\equiv\bm{0}$ for all $F\in\mathcal{F}_{h}\setminus\{F\}$, $v_{F^{\prime}}\equiv 0$ for all $F\in\mathcal{F}_{h}\setminus\{F\}$, and letting $v_{F}$ span $\mathbb{P}^{k}(F)$ yields the second equation in (54b). ∎

5 Properties of the discrete bilinear form

This section contains the proofs of the technical Lemmas 7 and 10.

5.1 Local coercivity and boundedness

Proof of Lemma 7.

Let a mesh element $T\in\mathcal{T}_{h}$ be fixed, and let $\underline{\mathrm{v}}_{T}\in\underline{\mathrm{U}}_{T}^{k}$.

(i) Coercivity. Taking $w=v_{T}\in\mathbb{P}^{k}(T)\subset\mathbb{P}^{k+2}(T)$ in (18) gives

[TABLE]

Using the Cauchy–Schwarz inequality to bound the first term in the right-hand side, the Cauchy–Schwarz and discrete trace (8a) inequalities to bound the second, and the Cauchy–Schwarz, discrete trace (8a) and inverse (8c) inequalities to bound the third, and simplifying we obtain:

[TABLE]

It remains to estimate the boundary terms inside the parentheses using the $\|{\cdot}\|_{\mathrm{a},T}$-seminorm.

(i.a) Bound on $\frac{\mathcal{A}_{T}^{+}}{h_{T}}\sum_{F\in\mathcal{F}_{T}}\|\bm{v}_{\bm{\nabla},F}-\bm{\nabla}v_{T}\|_{F}^{2}$. For all $F\in\mathcal{F}_{T}$, inserting $\pm\bm{\pi}_{F}^{k}\bm{\nabla}\left(p^{k+2}_{T}\underline{\mathrm{v}}_{T}-\pi^{k}_{T}p^{k+2}_{T}\underline{\mathrm{v}}_{T}\right)$ into the norm and using the linearity of $\bm{\pi}_{F}^{k}$ and the fact that it preserves polynomials in $\mathbb{P}^{k}(F)^{2}$ as a projector, we obtain

[TABLE]

where we have used the triangle inequality to pass to the second line. By the definition (22) of $\mathrm{s}_{T}$, we readily infer that

[TABLE]

Using the $L^{2}(F)^{2}$-boundedness of $\bm{\pi}_{F}^{k}$ followed by the discrete trace inequality (8a), we can write $|\mathfrak{T}_{2}|\lesssim h_{T}^{-\nicefrac{{1}}{{2}}}\|\bm{\nabla}(p^{k+2}_{T}\underline{\mathrm{v}}_{T}-\pi^{k}_{T}p^{k+2}_{T}\underline{\mathrm{v}}_{T})\|_{T}$. Then, by the approximation properties (11a) of $\pi^{k}_{T}$ with $l=k$, $m=1$, and $s=2$, we infer that

[TABLE]

so that

[TABLE]

Notice, in passing, that in the second bound in (59) we have used the fact that $k\geq 1$. Finally, the third term in the right-hand side of (58) can be estimated using the discrete trace (8a) and inverse (8c) inequalities together with the definition (22) of $\rm s_{T}$ as follows:

[TABLE]

Hence, multiplying (58) by $h_{T}^{-\nicefrac{{1}}{{2}}}\sqrt{\mathcal{A}^{+}_{T}}$, squaring, summing over $F\in\mathcal{F}_{T}$, using the above estimates for $\mathfrak{T}_{1}$, $\mathfrak{T}_{2}$, $\mathfrak{T}_{3}$, and recalling the uniform bound (7) on $\operatorname{card}(\mathcal{F}_{T})$, we have

[TABLE]

(i.b) Bound on $\frac{\mathcal{A}_{T}^{+}}{h_{T}^{3}}\sum_{F\in\mathcal{F}_{T}}\|v_{F}-v_{T}\|_{F}^{2}$. For all $F\in\mathcal{F}_{T}$, inserting $\pm\pi_{F}^{k}\left(p^{k+2}_{T}\underline{\mathrm{v}}_{T}-\pi_{T}^{k}p^{k+2}_{T}\underline{\mathrm{v}}_{T}\right)$ into the norm, and using the linearity of $\pi_{F}^{k}$ and $\pi_{T}^{k}$ together with the fact that they preserve polynomials up to degree $k$ as projectors, we have that

[TABLE]

By the definition (22) of $\mathrm{s}_{T}$, it is readily inferred that

[TABLE]

The second term can be estimated as follows:

[TABLE]

where we have used the $L^{2}(F)$-boundedness of $\pi^{k}_{F}$, the discrete trace inequality (8a), the uniform equivalence of face and element diameters (6) to replace $h_{F}$ with $h_{T}$, and the approximation property (11a) with $l=k$, $s=2$, and $m=0$. Again, here the hypothesis $k\geq 1$ is necessary to infer the second bound. Hence,

[TABLE]

Finally, using the discrete trace inequality (8a) followed by the definition (22) of $\mathrm{s}_{T}$, we have

[TABLE]

Multiplying (61) by $h_{T}^{-\nicefrac{{3}}{{2}}}\sqrt{\mathcal{A}^{+}_{T}}$, squaring, summing over $F\in\mathcal{F}_{T}$, using the above estimates for $\mathfrak{T}_{1}$, $\mathfrak{T}_{2}$, $\mathfrak{T}_{3}$, and recalling the uniform bound (7) on $\operatorname{card}(\mathcal{F}_{T})$, we arrive at

[TABLE]

(i.c) Conclusion. Combining (57), (62), and (60), the first inequality in (26) follows.

(ii) Boundedness. Taking $w=p^{k+2}_{T}\underline{\mathrm{v}}_{T}$ in (18), using the Cauchy–Schwarz, discrete trace (8a) and inverse inequalities (8c), and simplifying, we get

[TABLE]

which bounds the portion of $\|\underline{\mathrm{v}}_{T}\|_{\mathrm{a},T}$ stemming from the consistency term in (21). It remains to bound on the local stabilization terms in $\mathrm{s}_{T}(\underline{\mathrm{v}}_{T},\underline{\mathrm{v}}_{T})$.

(ii.a) Bound on $\frac{\mathcal{A}_{T}^{+}}{h_{T}^{4}}\|\pi^{k}_{T}(p^{k+2}_{T}\underline{\mathrm{v}}_{T}-v_{T})\|_{T}^{2}$. Inserting $\pm p^{k+2}_{T}\underline{\mathrm{v}}_{T}$ into the norm and using the triangle inequality, we have that

[TABLE]

For the first term, using the approximation property (11a) with $l=k$, $m=0$, and $s=2$, and (63), we get

[TABLE]

Once more, we use here the fact that $k\geq 1$. For the second term, inserting $0=-\pi_{T}^{1}p^{k+2}_{T}\underline{\mathrm{v}}_{T}+\pi_{T}^{1}v_{T}$ into the norm (see (17)) and using the triangle inequality, we obtain

[TABLE]

The approximation property (11a) with $l=1$, $m=0$, and $s=2$ gives $\|p^{k+2}_{T}\underline{\mathrm{v}}_{T}-\pi^{1}_{T}p^{k+2}_{T}\underline{\mathrm{v}}_{T}\|_{T}\lesssim h_{T}^{2}\|\bm{\nabla}^{2}p^{k+2}_{T}\underline{\mathrm{v}}_{T}\|_{T}$ and $\|v_{T}-\pi^{1}_{T}v_{T}\|_{T}\lesssim h_{T}^{2}\|\bm{\nabla}^{2}v_{T}\|_{T}$ so that, accounting for (63),

[TABLE]

Squaring (64), multiplying the resulting inequality by $\mathcal{A}_{T}^{+}/h_{T}^{4}$, and using the above estimates for $\mathfrak{T}_{1}$ and $\mathfrak{T}_{2}$ together with the uniform bound (7) on $\operatorname{card}(\mathcal{F}_{T})$, we conclude that

[TABLE]

(ii.b) Bound on $\frac{\mathcal{A}_{T}^{+}}{h_{T}}\sum_{F\in\mathcal{F}_{T}}\|\bm{\pi}^{k}_{F}(\bm{\nabla}p^{k+2}_{T}\underline{\mathrm{v}}_{T}-\bm{v}_{\bm{\nabla},F})\|_{F}^{2}$. For any $F\in\mathcal{F}_{T}$, inserting $\pm\bm{\nabla}v_{T}$ into the norm, invoking the linearity of $\bm{\pi}_{F}^{k}$ together with the fact that it preserves polynomials in $\mathbb{P}^{k}(F)^{2}$ as a projector, and using the triangle inequality, we have that

[TABLE]

where to pass to the second line we have used the $L^{2}(F)^{2}$-boundedness of $\bm{\pi}_{F}^{k}$, the discrete trace inequality (8a), and the inverse inequality (8c), while to pass to the third line we have estimated the first addend as the term $\mathfrak{T}_{2}$ in (64) (which requires again $k\geq 1$). Thus, squaring the above inequality, summing over $F\in\mathcal{F}_{T}$, multiplying it by $\mathcal{A}_{T}^{+}/h_{T}$, and using the uniform bound (7) on $\operatorname{card}(\mathcal{F}_{T})$, we finally infer

[TABLE]

(ii.c) Bound on $\frac{\mathcal{A}_{T}^{+}}{h_{T}^{3}}\sum_{F\in\mathcal{F}_{T}}\|\pi^{k}_{F}(p^{k+2}_{T}\underline{\mathrm{v}}_{T}-v_{F})\|_{F}^{2}$. For any $F\in\mathcal{F}_{T}$, inserting $\pm v_{T}$ into the norm, invoking the linearity of $\pi_{F}^{k}$ together with the fact that it preserves polynomials in $\mathbb{P}^{k}(F)$ as a projector, and using the triangle inequality, we infer that

[TABLE]

where to pass to the second line we have used the $L^{2}(F)$-boundedness of $\pi^{k}_{F}$ followed by the discrete trace inequality (8a) and the uniform equivalence of the element and face diameters expressed by (6), while to pass to the third line we have estimated the first addend as the term $\mathfrak{T}_{2}$ in (64) and, once more, we used the fact that $k\geq 1$. Hence, multiplying (67) by $h_{T}^{-\nicefrac{{3}}{{2}}}\sqrt{\mathcal{A}_{T}^{+}}$, squaring, summing over $F\in\mathcal{F}_{T}$, recalling (64), and using the uniform bound (7) on $\operatorname{card}(\mathcal{F}_{T})$, we conclude that

[TABLE]

(ii.d) Conclusion. The second inequality in (26) then follows combining (63), (68), and (66) and recalling the definition (25) of $\|{\cdot}\|_{\mathbb{A},T}$. ∎

5.2 Global coercivity, boundedness, and consistency

Proof of Lemma 10.

(i) Coercivity and boundedness. The norm equivalence (33) is an immediate consequence of Lemma 7 together with the definition (29) of the $\|{\cdot}\|_{\mathbb{A},h}$-norm.

(ii) Consistency. Let us prove (34). An element-wise integration by parts yields

[TABLE]

where we have used the fact that moments and Kirchhoff shear forces are continuous at interfaces owing to the regularity of $v$ (see (46b) for the expression of these continuity properties for the exact solution $u$) and that homogeneous boundary conditions are embedded in $\underline{\mathrm{U}}^{k}_{h,0}$. Now, let

[TABLE]

we have

[TABLE]

Thus, letting $\mathcal{E}_{h}(\underline{\mathrm{w}}_{h})\coloneq(\mbox{\rm{div}}\,\mbox{{div}}\mathbb{A}\bm{\nabla}^{2}v,w_{h})-\mathrm{a}_{h}(\widehat{\underline{\mathrm{v}}}_{h},\underline{\mathrm{w}}_{h})$, (69) and (71) yield

[TABLE]

By the definition (12) of the local energy projector, we have that

[TABLE]

Using the approximation properties (13) with $l=k+2$, $s=k+3$, and $m=2,3$, we infer that

[TABLE]

Moreover, for all $T\in\mathcal{T}_{h}$, we have $\mathrm{s}_{T}(\widehat{\underline{\mathrm{v}}}_{T},\underline{\mathrm{w}}_{T})\leq\mathrm{s}_{T}(\widehat{\underline{\mathrm{v}}}_{T},\widehat{\underline{\mathrm{v}}}_{T})^{\nicefrac{{1}}{{2}}}\mathrm{s}_{T}(\underline{\mathrm{w}}_{T}\underline{\mathrm{w}}_{T})^{\nicefrac{{1}}{{2}}}$; as for the first factor, by (24) we have $\mathrm{s}_{T}(\widehat{\underline{\mathrm{v}}}_{T},\widehat{\underline{\mathrm{v}}}_{T})^{\nicefrac{{1}}{{2}}}\lesssim\sqrt{\mathcal{A}_{T}^{+}}h_{T}^{k+1}|v|_{H^{k+3}(T)},$ whereas the second inequality in (26) gives $\mathrm{s}_{T}(\underline{\mathrm{w}}_{T},\underline{\mathrm{w}}_{T})^{\nicefrac{{1}}{{2}}}\lesssim\|\underline{\mathrm{w}}_{T}\|_{\mathbb{A},T}$, so that

[TABLE]

Using (72), (73), and (74) to estimate $\mathcal{E}_{h}(\underline{\mathrm{w}}_{h})$, and using the resulting bound in the supremum in (34) concludes the proof. ∎

6 Proof of Theorem 1

(i) Proof of (13a). We apply [26, Lemma 3]. Therefore, proving (13a) amounts to proving the following estimates:

[TABLE]

where $a\lesssim b$ means $a\leq Cb$ with $C>0$ as in (13).

We start by proving (75a). The definition (12) of $\varpi^{l}_{T}$ implies that

[TABLE]

where we have used the definition of $\mathcal{A}_{T}^{-}$ (see Section 2.1) in the first line, the characterization of $\varpi^{l}_{T}$ as $\arg\min_{z\in\mathbb{P}^{l}(T)}\|\mathbb{A}_{T}^{\nicefrac{{1}}{{2}}}\bm{\nabla}^{2}(v-z)\|_{T}$ in the second line, along with the definition of $\mathcal{A}_{T}^{+}$ and the $H^{2}$-stability of the $L^{2}$-orthogonal projector (resulting from (11a) with $s=m=2$) to conclude. Thus, using again the triangle inequality, we have that

[TABLE]

and (75a) is proved.

To prove (75b), we introduce the quantities $0=-\bm{\nabla}\pi^{1}_{T}\varpi^{l}_{T}v+\bm{\nabla}\pi^{1}_{T}v$ (recall the second condition in (12)) and $\pm\bm{\nabla}v$ inside the $L^{2}(T)$-norm of $\bm{\nabla}\varpi^{l}_{T}v$ to infer that

[TABLE]

where we have used the approximation estimate (11a) for $\pi^{1}_{T}$ with $m=1$ and $s=2$ together with the fact that, for any $w\in H^{2}(T)$, $|w|_{H^{2}(T)}\lesssim\|\bm{\nabla}^{2}w\|_{T}$ to estimate the first two terms, and (75a) to conclude.

The proof of (75c) is completely analogous. We obtain

[TABLE]

where we have used (11a) to estimate the first two addends in the first line, with $m=0$ and $s=2$ for the first one and with $m=0$ and $s=1$ for the second one. This concludes the proof of (13a).

(ii) Proof of (13b). For $m\leq s-1$, by applying the continuous trace inequality (8b) to $w=\partial^{\bm{\alpha}}(v-\varpi^{l}_{T}v)\in H^{1}(T)$ for all $\bm{\alpha}\in\mathbb{N}^{2}$ such that $\alpha_{1}+\alpha_{2}=m$, we have

[TABLE]

The conclusion follows using (13a) for $m$ and $m+1$ to bound the terms in the right-hand side.

7 Concluding remarks

Some concluding remarks are in order.

7.1 Computational cost of the method

It is worth to draw some conclusions from the numerical tests set forth in Section 3.7, with particular reference to the L-shaped domain case. Indeed, in many applications, as well as from a theoretical viewpoint, it is interesting to estimate the computational cost of a given numerical method. Here, we can evaluate the computational cost of our method by comparing the sizes of the matrices associated with the bilinear form $\mathrm{a}_{h}$, as well as the number of nonzero elements of such matrices666 The latter, in particular, gives an insight into the stencil of the method., upon varying the polynomial degree $k$ and the number of elements $N$. These two quantities are represented in Tables 3 and 4, respectively. As Table 3 shows, in certain cases (compare, for instance, the results given by the two choices $k=2$, $N=136$ and $k=3$, $N=34$ in Table 2) using polynomials of high order on coarse triangulations is more convenient than using polynomials of lower order on finer triangulations to obtain a given numerical value of the discrete energy to two significant digits.

7.2 Mixed formulations

The results of this paper concern the primal formulation (1) of the Kirchhoff–Love plate bending model problem. As it is well known, this problem admits dual and mixed formulations that have been the basis for the development of mixed and hybrid nonconforming finite elements (see, e.g., [16]). A HHO discretization based on a mixed formulation will make the object of a future work, as well as the study of its relation with the method presented here and its variations. We notice, in passing, that a similar study for a second-order elliptic problem has been carried out in [1] and, in a more general setting, in [12]. The latter works can be regarded as a generalization to new generation polytopal methods of the classical hybridization techniques of Arnold–Brezzi [5].

Acknowledgements. The authors are grateful to Franco Brezzi (IMATI Pavia) for the fruitful discussions that have helped shape up this work.