# Polynomial robust stability analysis for $$H$$ (div)-conforming finite elements for the Stokes equations

Polynomial robust stability analysis for $$H$$ (div)-conforming finite elements for the Stokes... Abstract In this article, we consider a discontinuous Galerkin method for the discretization of the Stokes problem. We use $$H(\textrm{div})$$-conforming finite elements, as they provide major benefits such as exact mass conservation and pressure-independent error estimates. The main aspect of this work lies in the analysis of high-order approximations. We show that the considered method is uniformly stable with respect to the polynomial order $$k$$ and provides optimal error estimates $$\| {\textbf{u} - \textbf{u}_h} \|_{1_h} + \| {{\it {\Pi}}^{Q_h} p-p_h} \|_{0} \le c \left(h/k \right)^s \| \textbf{u} \|_{s+1}$$. To derive these estimates, we prove a $$k$$-robust Ladyženskaja-Babuška-Brezzi (LBB) condition. This proof is based on a polynomial $$H^2$$-stable extension operator. This extension operator itself is of interest for the numerical analysis of $$C^0$$-continuous discontinuous Galerkin methods for fourth-order problems. 1. Introduction In this article, we consider the numerical solution of the Stokes equations on a bounded domain $${\it {\Omega}} \subset \mathbb{R}^2$$,   \begin{align} \label{int::stokesequations} \begin{array}{rcll} -\nu {\it {\Delta}} \textbf{u} + \nabla p &=& \textbf{f} &\textrm{in } {\it {\Omega}}, \\ \text{div}~ \textbf{u} & = & 0 & \textrm{in } {\it {\Omega}}, \end{array} \end{align} (1.1) with boundary conditions $$\textbf{u} = 0$$ on $$\partial {\it {\Omega}}$$, where $$\nu =\mathrm{const}$$ is the kinematic viscosity, $$\textbf{u}$$ is the velocity field, $$p$$ is the pressure and $$\textbf{f}$$ are external forces. The approximation of the Navier–Stokes problem is well analysed, and many different finite element methods have been introduced (see, for example, Girault & Raviart, 1986; Donea & Huerta, 2003; Glowinski, 2003; Elman et al., 2005). Furthermore, discontinuous Galerkin (DG) finite element methods for elliptic problems became popular (see, for example, Arnold et al., 2002; Houston et al., 2002; Rivière, 2008) and also for flow problems (Baker et al., 1990; Karakashian & Katsaounis, 2000; Cockburn et al., 2002; Schötzau et al., 2002; Toselli, 2002; Cockburn et al., 2004, 2005, 2007; Girault et al., 2005). In the last two mentioned works, an elementwise divergence-free discretization is constructed. In contrast, we consider a normal-continuous $$H(\operatorname{div})$$-conforming method introduced by Cockburn et al. (2007). This has different advantages such as local conservation, the possibility to use an upwinding scheme for convection-dominated flows and pressure robust (independent) error estimates due to exactly divergence-free velocity test functions (see Linke, 2014; Brennecke et al., 2015; Linke et al., 2016; Linke & Merdon, 2016). Note that this discretization was also used for the Brinkman problem in Könnö & Stenberg (2011). To reduce the computational costs of DG methods, we also want to mention hybrid DG (HDG) methods, where new variables are introduced on the skeleton and a static condensation technique is used for the element unknowns (see Cockburn et al., 2010, 2011; Egger & Schöberl, 2010; Nguyen et al., 2010, 2011; Egger & Waluga, 2013) and for $$H(\operatorname{div})$$-conforming methods (Lehrenfeld, 2010; Fu et al., 2016; Lehrenfeld & Schöberl, 2016). The method we use is well analysed with respect to the mesh size $$h$$ and provides optimal error estimates. The main contribution of this article is to show that the method is also uniformly stable with respect to the polynomial order $$k$$. For this purpose, we prove that the constant $$\beta$$ of the LBB condition   \begin{align*} \sup_{\textbf{0} \neq \textbf{v}_h \in \textbf{V}_h} \frac{b(\textbf{v}_h,q_h)}{\| \textbf{v}_h \|_{\textbf{V}_h}} \geq \beta \| q_h \|_{Q_h} \quad \forall\, q_h \in Q_h, \end{align*} is independent of the order $$k$$. Together with standard continuity and ellipticity estimates this leads to a stable high-order method. Note that small adaptations of our results also lead to polynomial robustness for HDG methods. High-order methods for incompressible flow problems are of theoretical and practical importance. In Bernardi & Maday (1997) and Karniadakis & Sherwin (2005), the authors consider a spectral method on the unit cube using polynomials of order $$k$$ and $$k-2$$ for the velocity and the pressure, respectively. The resulting method leads to $$\beta(k) = \mathscr{O}(k^{-\frac{d-1}{2}})$$, where $$d$$ is the space dimension. The same method on triangles is discussed in Schwab (1998) with similar results. Furthermore, the bad influence of a dependency on $$k$$ of the LBB constant for an iterative method for solving the Navier–Stokes equation was analysed. Understanding the problem, an improvement was achieved in Bernardi & Maday (1999). They used polynomials of partial order $$k$$ for the velocity and polynomials of total order at most $$k-1$$ for the pressure, resulting in a uniformly stable method. Similar achievements for $$hp$$ mixed finite element methods are accomplished in Stenberg & Suri (1996). Therein, different combinations of elements on quadrilaterals like continuous polynomials of order $$k$$ for the velocity and discontinuous polynomials of order $$k-2$$ for the pressure are discussed. An exact analysis is presented but again revealed a dependency on $$k$$. They also considered different tensor product spaces for each component of the velocity. A similar approach leading to an optimal exactly divergence-free method was presented in Zhang (2009) using polynomials of order $$k+1$$ in the $$x$$-direction and polynomials of order $$k$$ in the $$y$$-direction for the velocity in $$x$$-direction, and vice versa for the velocity in $$y$$-direction. Using proper degrees of freedom, this leads to a similar method on quadrilaterals as we use on triangles. Another approach, combining the tensor product structure on quadrilaterals and the advantage of approximating more complex geometries using triangles, is analysed in Su et al. (2016). The key of this method is to use the Duffy transformation and a proper pair of approximation spaces, which leads to $$\beta(k) = \mathscr{O}(k^{-\frac{1}{2}})$$ with the drawback of using rational functions for the approximation. We also want to mention the method considered in Ainsworth & Coggins (2002), where a uniformly stable approximation using a continuous ansatz for the velocity and pressure is presented. They adopted the ideas of Bernardi & Maday (1999) but enriched the pressure space to overcome the lack of convergence order that would appear using just a continuous version of this method. Considering continuous approximations, also the famous Taylor–Hood elements on triangles and quadrilaterals (see Brezzi & Falk, 1991; Boffi et al., 2013) have to be discussed. Although these methods were shown to be stable with respect to the mesh size $$h$$, numerical evidence predicts that they are not uniformly stable with respect to $$k$$. Of course, high-order methods were also used for discontinuous finite element methods. We want to mention the work of Toselli (2002) and Schötzau et al. (2002), where an analysis for $$hp$$-DG methods on quadrilaterals is presented but revealed a dependency on the order $$k$$, and also the work of Egger & Waluga (2013), where an HDG method on triangles and quadrilaterals with similar results is introduced. The rest of the article is structured in the following way. In Section 2 we present the Stokes equation and the considered discretization method. Furthermore, we present a short proof of the continuous divergence stability to motivate the existence of an $$H^2$$-stable polynomial extension operator, which is used to prove the main theorem in Section 3. In Section 4 we take a look at some numerical examples and finally present the construction of an $$H^2$$-stable polynomial extension operator in Section 5. 1.1 Preliminaries We assume an open bounded domain $${\it {\Omega}} \subset \mathbb{R}^2$$ with a Lipschitz boundary $${\it {\Gamma}}$$; thus for every point on the boundary there exists a Lipschitz-continuous mapping $${\it {\Phi}}^x$$. If this mapping is furthermore differentiable up to order $$m$$, we say $${\it {\Gamma}} \in \mathscr{C}^{m,1}$$. On $${\it {\Omega}}$$, we define a shape-regular triangulation $$\mathscr{T}$$. Furthermore, we assume $$\mathscr{T}$$ to be quasi-uniform; thus, there exists one global mesh size $$h$$ such that $$h \approx \textrm{diam}(T)\ \textrm{for all}\ T \in \mathscr{T}$$. The set of edges, with respect to the triangulation $$\mathscr{T}$$, will be defined as $$\mathscr{F}$$. We call $$\widehat{T} := \{(x,y): 0\le x \le 1, 0 \le y \le 1, x+y \le 1 \}$$ the reference element with the edges $$E_1:= \{(x,0): 0\le x \le 1 \}$$, $$E_2:= \{(x,1-x): 0\le x \le 1 \}$$ and $$E_3:= \{(0,y): 0\le y \le 1 \}$$, and define the interval $$E:=\{(x,0): -1\le x \le 1 \}$$. On all triangles we use $$\textbf{n}$$ and $$\tau$$ as symbols for the normal and tangential vectors. For all subsets $$\omega \subseteq{\it {\Omega}}$$ with $$\gamma := \partial \omega$$ we have the space $$L^2(\omega)$$ with the norm $$\| \cdot \|_{0,\omega}$$ and the Sobolev spaces $$H^1(\omega)$$, $$H^2(\omega), H^{1/2}(\gamma)$$ with the corresponding Sobolev norms $$\| \cdot \|_{s,\omega}$$. For better readability, we skip the index $$\omega$$ if it is clear which domain is considered. On the edge $$E_1$$ we define the weighted $$L^2$$- and $$H^{1/2}$$-norm of a function $$u$$ as   \begin{align*} \| u \|_{0^{*},E_1}^2 := \int_0^1 \left(\frac{1}{x} + \frac{1}{1-x}\right) u(x)^2 ~\mathrm{d} s \quad \textrm{and} \quad \| u \|_{1/2^*,E _1}^2 := | u |^2_{1/2,E _1} + \| u \|_{0^{*},E_1}^2. \end{align*} Furthermore, we use the closed subspaces   \begin{align*} L^2_0({\it {\Omega}}) &:= \{ q \in L^2({\it {\Omega}}): \textstyle\int_{\it {\Omega}} q ~\mathrm{d} x = 0\} \quad \textrm{and} \quad H^1_0({\it {\Omega}}) := \{ u \in H^1({\it {\Omega}}): \mathrm{tr} ~ u = 0 ~\textrm{on}~ \partial {\it {\Omega}}\}, \end{align*} and the polynomial spaces   \begin{align*} \mathscr{P}^m(\mathscr{T}) &:= \{ v : v|_T \in \mathscr{P}^m(T)~ \forall\, T \in \mathscr{T} \} = \prod_{ T \in \mathscr{T}} \mathscr{P}^m(T) \quad \textrm{and} \quad \pmb{\mathscr{P}}^m(\mathscr{T}) := [\mathscr{P}^m(\mathscr{T})]^2, \\ \mathscr{P}_{00}^m(E_1) &:= \{v \in \mathscr{P}^m(E_1): v(0) = v'(0) = v(1) = v'(1) = 0\}. \end{align*} Also we define the following subspaces of vectorial polynomials on the reference triangle:   \begin{align*} \pmb{\mathscr{P}}^m_{\tau}(\widehat{T}):= \{\textbf{v} \in \pmb{\mathscr{P}}^m(\widehat{T}): \textstyle\int_{\partial \widehat{T}} \textbf{v} \cdot \tau = 0\} \quad \textrm{and} \quad \pmb{\mathscr{P}}^m_{\textbf{n}}(\widehat{T}):= \{\textbf{v} \in \pmb{\mathscr{P}}^m(\widehat{T}): \textstyle\int_{\partial \widehat{T}} \textbf{v} \cdot \textbf{n} = 0\}. \end{align*} In this article, we use an index notation for partial derivatives; thus for an arbitrary function $$u$$ we write $$u_{,x} := \frac{\partial u}{\partial x}$$ and $$u_{,y} := \frac{\partial u}{\partial y}$$ and in a similar manner for second-order derivatives. We write $$(x,y)^{\mathrm{t}}$$ as the transposed vector of $$(x,y)$$ and use $$\perp$$ as symbol for a counterclockwise rotation by $$\pi/2$$, thus $$(x,y)^\perp := (-y,x)$$. By this we also define the two-dimensional operator $$\text{curl}~ = \nabla^\perp$$. Finally, we introduce the expression $$a \preccurlyeq b$$ when there exists a constant $$c$$ independent of $$a, b$$, the polynomial order and the mesh size such that $$a \le c b$$. 2. Discretization of the Stokes problem In this section, we present the discretization of the stationary incompressible Stokes equations (1.1) from Cockburn et al. (2007). They use mixed-order finite element spaces with the polynomial orders $$k$$ and $$k-1$$ for the velocity and the pressure, respectively. To ensure a local conservative and energy-stable method, we provide exactly divergence-free velocity fields using an $$H(\operatorname{div})$$-conforming method; thus every discrete velocity field $$\textbf{u}_h$$ is in   \begin{align*} H(\operatorname{div})({\it {\Omega}}) := \{ \textbf{u} \in [L^2({\it {\Omega}})]^2: \text{div}~ \textbf{u} \in L^2({\it {\Omega}}) \}. \end{align*} To ensure $$\textbf{u}_h \in H(\operatorname{div})({\it {\Omega}})$$, we demand normal continuity across each edge resulting in the approximation space for the velocity   \begin{align*} \textbf{V}_h := \{ \textbf{u}_h \in \pmb{\mathscr{P}}^k(T): [\![{\textbf{u}_h \cdot n}]\!] = 0 ~ \forall\, E \in \mathscr{F} \} \subset H(\operatorname{div},{\it {\Omega}}), \end{align*} where $$[\![{\cdot}]\!]$$ is the usual jump operator on interior edges and just the identity on boundary edges. For the pressure space, we assume no continuity across edges:   \begin{align*} Q_h :=\prod\limits_{T \in \mathscr{T}}\mathscr{P}^{k-1}(T) \cap L^2_0({\it {\Omega}}). \end{align*} Note, that this pair of finite element spaces fulfils $$\text{div}~ \textbf{V}_h = Q_h$$ and thus, a weakly incompressible velocity field $$\textbf{u}_h \in \textbf{V}_h$$ is also exactly divergence-free:   \begin{align*} \int_{\it {\Omega}} \text{div}~{\textbf{u}_h} ~ q ~\mathrm{d} x= 0 \quad \forall\, q \in Q_h \quad \Rightarrow \quad \text{div}~{\textbf{u}_h} = 0 \quad \textrm{in } {\it {\Omega}}. \end{align*} Furthermore, using $$\pmb{\{} \cdot \pmb{\}}$$ as the symbol for the mean value on an edge $$E \in \mathscr{F}$$ we define the bilinear form are   \begin{align}\label{bilinearform::a} a(\textbf{u}_h,\textbf{v}_h) &:= \sum\limits_{T \in \mathscr{T}} \int_T \nu \nabla \textbf{u}_h : \nabla \textbf{v}_h ~\mathrm{d} x - \sum\limits_{E \in \mathscr{F}} \int_{E} \nu \pmb{\{} {\nabla \textbf{u}_h \cdot n} \pmb{\}} [\![{\textbf{v}_h \cdot \tau}]\!] ~\mathrm{d} s\nonumber\\ &\quad{} -\sum\limits_{E \in \mathscr{F}} \int_{E} \nu \pmb{\{} {\nabla \textbf{v}_h \cdot n} \pmb{\}} [\![{\textbf{u}_h \cdot \tau}]\!] ~\mathrm{d} s + \sum\limits_{E \in \mathscr{F}} \int_{E} \nu \frac{\alpha k^2}{h} [\![{\textbf{u}_h \cdot \tau}]\!] [\![{\textbf{v}_h \cdot \tau}]\!] ~\mathrm{d} s, \end{align} (2.1) where $$\alpha >0$$ and $$k^2/h$$ is a standard stability coefficient (see Houston et al., 2002) and the bilinear form and linear form   \begin{align} \label{bilinearform::b} b(\textbf{u}_h, q_h) := - \sum\limits_{T \in \mathscr{T}} \int_T \text{div}~{\textbf{u}_h} ~ q_h ~\mathrm{d} x \quad \textrm{and} \quad l(\textbf{v}_h):= \sum\limits_{T \in \mathscr{T}} \int_T \textbf{f} \cdot \textbf{v}_h ~\mathrm{d} x. \end{align} (2.2) The discrete Stokes problem now reads, find $$(\textbf{u}_h, p_h)$$ in $$\textbf{V}_h\times Q_h$$ such that   \begin{align} \label{dis::stokesproblem} \begin{array}{rcll} a(\textbf{u}_h,\textbf{v}_h) + b(\textbf{v}_h, p_h) &=& l(\textbf{v}_h)\quad &\forall\, \textbf{v}_h \in \textbf{V}_h, \\ b(\textbf{u}_h, q_h) & = & 0 \quad & \forall\, q_h \in Q_h. \end{array} \end{align} (2.3) On the spaces $$Q_h$$ and $$\textbf{V}_h$$, we use the $$L^2$$-norm $$\| {\cdot} \|_{0}$$ and $$\| {\textbf{v}_h} \|_{1_h}^2 := \sum\limits_{T \in \mathscr{T}} \| {\nabla \textbf{v}_h} \|_{0,T}^2 + \sum\limits_{E \in \mathscr{F}} \frac{k^2}{h} \| {[\![{\textbf{v}_h \cdot \tau}]\!]} \|_{0,E}^2$$, respectively. We have chosen homogeneous Dirichlet boundary conditions, but other choices are possible as well. In the $$H(\operatorname{div})$$-conforming DG method, we distinguish between essential boundary conditions for the normal component and for the tangential components. The essential normal boundary condition is incorporated into the finite element space, whereas an essential tangential boundary condition is treated in the DG sense (see Georgoulis et al., 2010). Natural boundary conditions are treated in the common way as well. This distinction is particularly useful for the so-called slip boundary conditions, where only the normal velocity is constrained to $$\textbf{u} \cdot \textbf{n} = 0$$. Lemma 2.1 For a proper choice of the stabilization parameter $$\alpha >0$$ in (2.1), there exist constants $$\alpha_1>0, \alpha_2>0$$ and $$\alpha_3>0$$ independent of the mesh size $$h$$ and the polynomial order $$k$$ such that $$a(\cdot,\cdot)$$ is coercive,   \begin{align*} a(\textbf{v}_h,\textbf{v}_h) \ge \nu \alpha_3 \| {\textbf{v}_h} \|_{1_h}^2 \quad \forall\, \textbf{v}_h \in \textbf{V}_h, \end{align*} and $$a(\cdot,\cdot)$$ and $$b(\cdot,\cdot)$$ are continuous,   \begin{align*} |a(\textbf{u}_h,\textbf{v}_h)| \le \nu \alpha_1 \| {\textbf{u}_h} \|_{1_h}\| {\textbf{v}_h} \|_{1_h} \quad \forall\, \textbf{v}_h, \textbf{u}_h \in \textbf{V}_h, \quad |b(\textbf{u}_h,q_h)| \le \alpha_2 \| {\textbf{u}_h} \|_{1_h} \| {q_h} \|_{0} \quad \forall\, \textbf{u}_h \in \textbf{V}_h, \forall\, q_h \in Q_h. \end{align*} Proof. The continuity properties of $$a(\cdot,\cdot)$$ and $$b(\cdot,\cdot)$$ follow by the definition of the norm $$\| {\cdot} \|_{1_h}$$. The coercivity follows with similar arguments to Stamm & Wihler (2010). □ Theorem 2.2 There exists a constant $$\beta > 0$$ independent of the polynomial order $$k$$ and the mesh size $$h$$ such that   \begin{align*} \sup\limits_{\textbf{0} \neq \textbf{v}_h \in \textbf{V}_h} \frac{b(\textbf{v}_h, q_h)}{\| {\textbf{v}_h} \|_{1_h}} \ge \beta \| {q_h} \|_{0} \quad \forall\, q_h \in Q_h. \end{align*} Proof. The proof is presented in Section 3. □ Theorem 2.3 Assume $$(\textbf{u}_h, p_h)$$ in $$\textbf{V}_h\times Q_h$$ is the solution of the discrete problem (2.3) and $$(u,p)$$ is the exact solution of the Stokes problem (1.1). Furthermore, let us assume regularity $$\textbf{u} \in [H^{s+1}({\it {\Omega}})]^2$$ and $$p \in H^s({\it {\Omega}})$$. Then, there exists a constant $$c_{\rm err}$$ independent of the mesh size $$h$$ and the polynomial order $$k$$ such that for $$s \ge 1$$ and $$s \le k$$ it holds that   \begin{align}\label{optimalerrorest} \| {\textbf{u} - \textbf{u}_h} \|_{1_h} + \| {{\it {\Pi}}^{Q_h} p-p_h} \|_{0} \le c_{\rm err} \left(\frac{h}{k} \right)^{s} \| \textbf{u} \|_{s+1}, \end{align} (2.4) where $${\it {\Pi}}^{Q_h}$$ is the $$L^2$$ projector onto $$Q_h$$. Proof. In a first step, we discretize the Poisson problem $$-\nu {\it {\Delta}} \textbf{u} = f - \nabla p$$ using a DG approximation, i.e., let $$\textbf{w}_h$$ be the solution of   \begin{align*} a(\textbf{w}_h, \textbf{v}_h) = l(v_h) - b(\textbf{v}_h, p) \quad \forall\, \textbf{v}_h \in \textbf{V}_h, \end{align*} where we used integration by parts for $$\int_{\it {\Omega}} \nabla p \cdot \textbf{v}_h ~\mathrm{d} x$$ and $$\textbf{v}_h \in H(\operatorname{div})({\it {\Omega}})$$. Using the estimate from Stamm & Wihler (2010, Section 3.2), including the properties of the $$L^2$$ projector on triangles in Chernov (2012, Theorem 1.1), we get   \begin{align} \label{errorestlaplace} \| {\textbf{u} - \textbf{w}_h} \|_{1_h} \le \tilde{c}_{\mathrm{err}} \left(\frac{h}{k} \right)^{s} \| \textbf{u} \|_{s+1}. \end{align} (2.5) Since $$(\textbf{u}_h, p_h)$$ is the solution of the discrete problem (2.3) we have   \begin{align*} \begin{array}{rcll} a(\textbf{u}_h,\textbf{v}_h) + b(\textbf{v}_h, p_h) &=& b(\textbf{v}_h, p) + a(\textbf{w}_h, \textbf{v}_h) \quad &\forall\, \textbf{v}_h \in \textbf{V}_h, \\ b(\textbf{u}_h, q_h) & = & 0 \quad & \forall\, q_h \in Q_h, \end{array} \end{align*} and thus   \begin{align*} \begin{array}{rcll} a(\textbf{u}_h - \textbf{w}_h,\textbf{v}_h) + b(\textbf{v}_h, p_h-p) &=& 0 \quad &\forall\, \textbf{v}_h \in \textbf{V}_h, \\ b(\textbf{u}_h- \textbf{w}_h, q_h) & = & -b(\textbf{w}_h, q_h) \quad & \forall\, q_h \in Q_h. \end{array} \end{align*} Due to the property $$\text{div}~ \textbf{V}_h = Q_h$$, we replace the term $$b(\textbf{v}_h, p_h-p)$$ in the first row by $$b(\textbf{v}_h, p_h - {\it {\Pi}}^{Q_h} p)$$. Using $$p_h - {\it {\Pi}}^{Q_h} p \in Q_h$$ and the standard stability estimate of saddle point problems (see, for example, Boffi et al., 2013, Theorem 4.2.3), we get   \begin{align*} \| {\textbf{u}_h - \textbf{w}_h} \|_{1_h} + \frac{\beta}{\alpha_1} \| p_h - {\it {\Pi}}^{Q_h} p \|_0 &\le \left(1+ \frac{2\alpha_1}{\alpha_3\beta}\right) \| \text{div}~{\textbf{w}_h} \|_0 \\ &= \left(1+ \frac{2\alpha_1}{\alpha_3\beta}\right) \| \text{div}~{\textbf{w}_h} - \underbrace{\text{div}~{\textbf{u}}}_{=0} \|_0 \le \left(1+ \frac{2\alpha_1}{\alpha_3\beta}\right) \| {\textbf{u} - \textbf{w}_h} \|_{1_h}. \end{align*} The estimation (2.4) follows with (2.5), the triangle inequality and the robustness of the constants $$\alpha_1, \alpha_3$$ and $$\beta$$ with respect to the mesh size $$h$$ and $$k$$. □ Remark 2.4 The introduced method can also be used for a non-quasi-uniform triangulation, thus using triangles with different sizes $$h_T$$ and furthermore individual polynomial degrees $$k_T$$. By that we get a similar local error estimation to Stamm & Wihler (2010),   \begin{align*} \| {\textbf{u} - \textbf{u}_h} \|_{1_h} + \| {{\it {\Pi}}^{Q_h} p-p_h} \|_{0} \le c_{\rm err} \sqrt{ \sum\limits_{T \in \mathscr{T}}\left(\frac{h_T}{k_T} \right)^{2s} \| \textbf{u} \|_{s+1,T}^2}. \end{align*} 2.1 Continuous LBB condition as motivation for an $$H^2$$-extension In this section, we present a proof for the infinite-dimensional version of the LBB condition of the Stokes problem as can be found in Bernardi & Maday (1997). For this, we define the velocity space $$\textbf{V} := [H^1_0({\it {\Omega}})]^2$$, the pressure space $$Q := L^2_0({\it {\Omega}})$$ and show that   \begin{align*} \sup\limits_{\textbf{0} \neq \textbf{v} \in \textbf{V}} \frac{b(\textbf{v}, q)}{\|\textbf{v}\|_1} \ge \beta_\infty \| {q} \|_{0} \quad \forall\, q \in Q. \end{align*} Note that the LBB condition is equivalent to the existence of an $$H^1$$-stable right inverse of the divergence operator. Theorem 2.5 Let $${\it {\Omega}} \subset \mathbb{R}^2$$ be a bounded domain with a smooth Lipschitz boundary $$\partial {\it {\Omega}} \in \mathscr{C}^{1,1}$$. The divergence operator from $$\textbf{V}$$ to $$Q$$ is onto, so for every $$q \in Q$$ there exists a $$\textbf{v} \in \textbf{V}$$ such that   \begin{align*} \text{div}~{\textbf{v}} = q \quad \textrm{and} \quad \|\textbf{v}\|_1 \preccurlyeq \| {q} \|_{0}. \end{align*} Proof. Let $$q$$ be an arbitrary function in $$Q$$. In the first step, we consider the Poisson problem $${\it {\Delta}} \varphi = q$$ in $${\it {\Omega}}$$ with Neumann boundary conditions $$\nabla \varphi \cdot \textbf{n} = 0$$ on $$\partial {\it {\Omega}}$$. Due to the zero mean value of $$q$$ this problem has a unique solution in $$H^1({\it {\Omega}}) / \mathbb{R}$$. Now set $$\textbf{v} := \nabla \varphi$$. It follows that $$\text{div}~{\textbf{v}} = {\it {\Delta}} \varphi = q$$ and, using a regularity result for the Poisson problem, also $$\| v\|_1 = \| \varphi \|_2 \preccurlyeq \| {q} \|_{0}$$. Furthermore, note that we already have $$\textbf{v} \cdot n = \nabla \varphi \cdot n = 0$$ on the boundary $$\partial {\it {\Omega}}$$, so the idea is to construct a correction for the tangential component to satisfy the zero boundary values of $$\textbf{V}$$. For this, we seek for a function $$\psi \in H^2({\it {\Omega}})$$, such that   \begin{align*} \psi = 0 \quad \textrm{on } \partial {\it {\Omega}} \quad \textrm{and} \quad \frac{\partial \psi}{\partial n} = - \textbf{v} \cdot \tau \quad \textrm{on } \partial {\it {\Omega}} \quad \textrm{and} \quad \| \psi \|_2 \preccurlyeq \| \textbf{v} \|_1. \end{align*} The existence of $$\psi$$ holds true as we assumed a smooth boundary (see Bernardi & Maday, 1997, Theorem 1.12). Now set $$\tilde{\textbf{v}}:= \textbf{v} + \text{curl}~{\psi}$$ to get $$\text{div}~{\tilde{\textbf{v}}} = \text{div}~{\textbf{v}} + \text{div}~{\text{curl}~{\psi}} = q$$ in $${\it {\Omega}}$$ and on the boundary $$\partial {\it {\Omega}}$$:   \begin{align*} \tilde{\textbf{v}} \cdot \textbf{n} = \textbf{v} \cdot \textbf{n} + \text{curl}~{\psi} \cdot \textbf{n} = \nabla \psi \cdot \tau = 0 \quad \textrm{and} \quad \tilde{\textbf{v}} \cdot \tau = \textbf{v} \cdot \tau + \text{curl}~{\psi} \cdot \tau = \textbf{v} \cdot \tau + \nabla \psi \cdot \textbf{n} = 0. \end{align*} Finally, due to the $$H^2$$-continuity of $$\psi$$, we get $$\|\tilde{\textbf{v}}\|_1 = \|\textbf{v}\|_1 + \|\text{curl}~{\psi}\|_1 \preccurlyeq \|\textbf{v}\|_1 \preccurlyeq \| {q} \|_{0}$$. □ It follows immediately (with a scaling of $$\tilde{\textbf{v}}$$ by $$-1$$) that   \begin{align*} \sup\limits_{\textbf{0} \neq \textbf{v} \in \textbf{V}} \frac{b(\textbf{v}, q)}{\|\textbf{v}\|_1} \succcurlyeq \frac{-\int_{\it {\Omega}} \text{div}~{\tilde{\textbf{v}}}{q} ~\mathrm{d} x }{\|\tilde{\textbf{v}}\|_1} \succcurlyeq \frac{\| {q} \|_{0}^2}{\| {q} \|_{0}} \succcurlyeq \| {q} \|_{0}. \end{align*} The crucial part of this proof was the existence of the correction $$\psi$$ that is stable in the $$H^2$$-norm. This holds true in the case of a smooth boundary of $${\it {\Omega}}$$, but as the boundary of an element $$T \in \mathscr{T}$$ is just in $$\mathscr{C}^{0,1}$$, we get a problem when we want to adapt this proof to show the main Theorem 2.2. Such $$H^2$$-stable extensions of boundary values for nonregular boundaries $$\partial {\it {\Omega}}$$, which are defined as a union of a finite number of regular parts, are considered in Grisvard (1985) and Bernardi & Maday (1992). For this, they assume that compatibility conditions are satisfied at the points where two parts of the boundary gather. Those conditions are quite restrictive, thus do not hold true for all traces of polynomials. For this reason, such a proof cannot be used, for example, in the case of continuous velocity elements. Considering only normal continuous approximations with a tangential continuity only in a DG sense creates enough freedom to construct such an $$H^2$$-extension. 3. Robust high-order LBB estimation In this section, we present the proof of Theorem 2.2 in similar steps as in the proof of Theorem 2.5. For this, we assume the existence of a stable $$H^2$$-extension which is then presented in Section 5. Theorem 3.1 ($$H^2$$-extension) For every $$k$$ there exists an operator $$\mathscr{E} : \pmb{\mathscr{P}}_{\textbf{n}}^k(\widehat{T}) \rightarrow \mathscr{P}^{k+1}(\widehat{T}),$$ such that for $$\textbf{u}_h \in \pmb{\mathscr{P}}_{\textbf{n}}^k(\widehat{T})$$ it holds that   \begin{align} &\text{curl}~{\mathscr{E}(\textbf{u}_h)} \cdot \textbf{n} = \textbf{u}_h \cdot \textbf{n} \quad \textrm{on } \partial \widehat{T}, \label{lbb:h2ext:propA} \\ \end{align} (3.1)  \begin{align} & \|(\textbf{u}_h - \text{curl}~{\mathscr{E}(\textbf{u}_h)}) \cdot \tau \|_{0,\partial \widehat{T}} \preccurlyeq \frac{1}{k} \| \textbf{u}_h \|_{1,\widehat{T}}, \label{lbb:h2ext:propB} \\ \end{align} (3.2)  \begin{align} & \| \mathscr{E}(\textbf{u}_h) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u}_h \|_{1,\widehat{T}}.\label{lbb:h2ext:propC} \end{align} (3.3) Proof. See Section 5. □ We first show the LBB condition on the reference triangle $$\widehat{T}$$. For this, we define the spaces $$\widehat{\textbf{V}}_h := \pmb{\mathscr{P}}^k(\widehat{T})$$ with $$\widehat{\textbf{V}}_{h,0} := \{ \hat{\textbf{v}}_h \in \widehat{\textbf{V}}_h: \hat{\textbf{v}}_h \cdot n = 0 ~ \textrm{on } \partial \widehat{T}\}$$ and $$\widehat{Q}_h := \mathscr{P}^{k-1}(\widehat{T}) \cap L^2_0(\widehat{T})$$. The norm $$\| {\cdot} \|_{1_h}$$ on $$\widehat{\textbf{V}}_h$$ now reads   \begin{align*} \|\hat{\textbf{v}}_h\|^2_{1_h,\widehat{T}} = \|{\nabla \hat{\textbf{v}}_h}\|^2_{0,\widehat{T}} + \sum\limits_{E \in \partial \widehat{T}} k^2 \| \hat{\textbf{v}}_h \cdot \tau \|^2_{0,E} \quad \forall\, \hat{\textbf{v}}_h \in \widehat{\textbf{V}}_h. \end{align*} Theorem 3.2 (Local LBB condition) There exists a constant $$\beta > 0$$ independent of the polynomial order $$k$$ such that   \begin{align*} \sup\limits_{\textbf{0} \neq \hat{\textbf{v}}_h \in \widehat{\textbf{V}}_{h,0}} \frac{b(\hat{\textbf{v}}_h, \hat{q}_h)}{\|\hat{\textbf{v}}_h\|_{1_h,\widehat{T}}} \ge \beta \|\hat{q}_h\|_{0,\widehat{T}} \quad \forall\, \hat{q}_h \in \widehat{Q}_h. \end{align*} Proof. Let $$\hat{q}_h \in \widehat{Q}_h$$ be an arbitrary function. For a fixed point $$\hat{x} \in \widehat{T}$$, we define a local Poincaré operator $$\mathscr{Z}_{\hat{x}}:\mathscr{P}^{k-1}(\widehat{T}) \rightarrow [\mathscr{P}^{k}(\widehat{T})]^2$$ as introduced in Costabel & McIntosh (2010) by   \begin{align*} \hat{q}_h(x) \mapsto \mathscr{Z}_{\hat{x}}(\hat{q}_h)(x):= -(x-\hat{x}) \int_0^1 t \hat{q}_h(\gamma(t)) ~\mathrm{d} t, \end{align*} with $$\gamma(t) := \hat{x} + t(x-\hat{x})$$, and by integrating over all points in $$\widehat{T}$$ we furthermore define   \begin{align*} \hat{\textbf{u}}_{h}^1(x) := -\int_{\widehat{T}} \theta(\hat{x}) \mathscr{Z}_{\hat{x}}(\hat{q}_h)(x) ~\mathrm{d} \hat{x}, \end{align*} where $$\theta \in C^\infty_0(\widehat{T})$$ is a smooth function. We observe that $$-\text{div}~{\hat{\textbf{u}}_{h}^1} = \hat{q}_h$$ and $$\| \hat{\textbf{u}}_{h}^1\|_{1,\widehat{T}} \preccurlyeq \|\hat{q}_h\|_{0,\widehat{T}}$$ (see Costabel & McIntosh, 2010, Corollary 3.4). As $$\hat{q}_h \in L^2_0(\widehat{T})$$, we see that $$\hat{\textbf{u}}_{h}^1 \in \pmb{\mathscr{P}}_{\textbf{n}}^k(\widehat{T}),$$ so we apply the extension operator of Theorem 3.1 to define $$\hat{\textbf{u}}_{h}^2:= \text{curl}~{\mathscr{E}(\hat{\textbf{u}}_{h}^1)}$$. By that we get for $$\hat{\textbf{u}}_h := \hat{\textbf{u}}_{h}^1 - \hat{\textbf{u}}_{h}^2$$ and property (3.1),   \begin{align*} -\text{div}~{\hat{\textbf{u}}_h} = -\text{div}~{\hat{\textbf{u}}_{h}^1} + \text{div}~{\hat{\textbf{u}}_{h}^2} = \hat{q}_h \textrm{ in } \widehat{T} \quad \textrm{and} \quad \hat{\textbf{u}}_h \cdot n = \hat{\textbf{u}}_{h}^1 \cdot n - \hat{\textbf{u}}_{h}^2 \cdot n = 0 \textrm{ on } \partial \widehat{T}. \end{align*} Together with (3.3) and (3.2) we get   \begin{align*} \|\hat{\textbf{u}}_h\|^2_{1_h,\widehat{T}} &=\|{\nabla \hat{\textbf{u}}_{h}^1 - \nabla \hat{\textbf{u}}_{h}^2}\|^2_{0,\widehat{T}} + \sum\limits_{E \in \partial \widehat{T}} k^2 \| (\hat{\textbf{u}}_{h}^1 - \hat{\textbf{u}}_{h}^2) \cdot \tau \|^2_{0,E} \preccurlyeq \|\nabla \hat{\textbf{u}}_{h}^1\|^2_{0,\widehat{T}} + \|{\nabla \hat{\textbf{u}}_{h}^2}\|^2_{0,\widehat{T}} +\| \hat{\textbf{u}}_h \|_{1,\widehat{T}}^2 \\ &\preccurlyeq \| \hat{\textbf{u}}_h \|_{1,\widehat{T}}^2 + \|{\nabla \hat{\textbf{u}}_{h}^2}\|^2_{0,\widehat{T}} = \| \hat{\textbf{u}}_h \|_{1,\widehat{T}}^2 + | \mathscr{E}(\hat{\textbf{u}}_{h}^1) |^2_{2,\widehat{T}} \preccurlyeq \| \hat{\textbf{u}}_h \|_{1,\widehat{T}}^2 \preccurlyeq \|\hat{q}_h\|_{0,\widehat{T}}^2. \end{align*} As $$\hat{\textbf{u}}_h \in \widehat{\textbf{V}}_{h,0}$$ we can bound the supremum from below   \begin{align*} \sup\limits_{\textbf{0} \neq \hat{\textbf{v}}_h \in \widehat{\textbf{V}}_{h,0}} \frac{b(\hat{\textbf{v}}_h, \hat{q}_h)}{\|\hat{\textbf{v}}_h\|_{1_h,\widehat{T}}} \succcurlyeq \frac{-\int_{\widehat{T}} \text{div}~{\hat{\textbf{u}}_h} ~{\hat{q}_h} ~\mathrm{d} x }{\|\hat{\textbf{u}}_h\|_{1_h,\widehat{T}}} \succcurlyeq \frac{ \|\hat{q}_h\|_{0,\widehat{T}}^2}{ \|\hat{q}_h\|_{0,\widehat{T}}} \succcurlyeq \|\hat{q}_h\|_{0,\widehat{T}}.\quad\square \end{align*} In the following, we derive a similar LBB estimate for each physical triangle. For this purpose, we define for all $$T \in \mathscr{T}$$ the spaces $$\textbf{V}_h(T) := \pmb{\mathscr{P}}^k(T)$$ with $$\textbf{V}_{h,0}(T) := \{ \textbf{v}_h \in \textbf{V}_h(T): \textbf{v}_h \cdot n = 0 ~ \textrm{on } \partial T\}$$ and $$Q_h(T) := \mathscr{P}^{k-1}(T) \cap L^2_0(T)$$. Corollary 3.3 For each triangle $$T \in \mathscr{T}$$ there exists a constant $$\beta > 0$$ independent of the polynomial order $$k$$ and the mesh size $$h$$ such that   \begin{align*} \sup\limits_{\textbf{0} \neq \textbf{v}_h \in \textbf{V}_{h,0}(T)} \frac{b(\textbf{v}_h, q_h)}{\|\textbf{v}_h\|_{1_h}} \ge \beta \|q_h\|_{0,T} \quad \forall\, q_h \in Q_h(T). \end{align*} Proof. Let $$F$$ be the mapping from the reference triangle $$\widehat{T}$$ to $$T$$. For an arbitrary function $$q_h \in Q_h(T)$$ we define $$\hat{q}_h := \mathscr{J} q_h \in \widehat{Q}_h$$ where $$\mathscr{J} := \det F'$$. On the reference triangle we use Theorem 3.2 to find a function $$\hat{\textbf{u}}_h \in \widehat{\textbf{V}}_{h,0}$$ with $$-\text{div}~{\hat{\textbf{u}}_h} = \hat{q}_h$$ and $$\|\hat{\textbf{u}}_h\|_{1_h,\widehat{T}} \preccurlyeq \| \hat{q}_h\|_{0,\widehat{T}}$$ (compare Theorem 2.5). Using the Piola transformation $$\mathbb{P}$$ we define $$\textbf{u}_h := \mathbb{P}(\hat{\textbf{u}}_h) = \frac{1}{\mathscr{J}} F' \hat{\textbf{u}}_h$$. As $$\mathbb{P}$$ preserves the normal component we have $$\textbf{u}_h \in \textbf{V}_{h,0}(T)$$ and due to the proper mapping of the divergence also $$-\text{div}~{\textbf{u}_h} = q_h$$. Next note that due to the shape regularity of $$\mathscr{T}$$ we have $$|| F'|| \approx h$$ and $$| \mathscr{J}| \approx h^2$$. By standard scaling arguments it follows that $$\frac{1}{h} \|\hat{\textbf{u}}_h\|_{1_h,\widehat{T}} \approx \|\textbf{u}_h\|_{1_h}$$ and $$\| \hat{q}_h\|_{0,\widehat{T}} \approx h \|q_h\|_{0,T}$$; thus $$\|\textbf{u}_h\|_{1_h} \preccurlyeq \| q_h\|_{0,T}$$. The proof is concluded with similar steps as at the end of the proof of Theorem 3.2. □ Proof of Theorem 2.2. For the proof we assume that $$k \ge 2$$. For the analysis of the lowest-order case, we refer to Lehrenfeld (2010). We construct a Fortin operator $${\it {\Pi}}_{F}$$ with   \begin{align}\label{fortin} b({\it {\Pi}}_{F} \textbf{u} -\textbf{u}, q_h) = 0 \quad \forall\, q_h \in Q_h \quad \textrm{and} \quad \| {{\it {\Pi}}_{F} \textbf{u} } \|_{1_h} \le c_{F} \|\textbf{u} \|_1, \end{align} (3.4) where $$c_{\mathrm{F}}$$ is a robust constant with respect to $$h$$ and $$k$$. Then, Theorem 2.2 follows from Boffi et al. (2013, Proposition 5.4.2). We define the operator $${\it {\Pi}}_{F}$$ as the sum of a low-order operator $${\it {\Pi}}_{F}^1$$ and a correction operator $${\it {\Pi}}_{F}^2$$. The first operator is the Fortin operator for the pair $${\left(\pmb{\mathscr{P}}^2(\mathscr{T}) \cap [C^0({\it {\Omega}})]^2 \right) \times \mathscr{P}^0(\mathscr{T})}$$, (see Boffi et al., 2013, Section 8.4), which is uniformly continuous in $$h$$. For the second operator, let $$(\textbf{w}_h^{\mathrm{T}},r_h^{\mathrm{T}}) \in \textbf{V}_{h,0}(T) \times Q_h(T)$$ be the solution of the local correction problem   \begin{align*} \begin{array}{rcll} a^{\mathrm{T}}(\textbf{w}_h^{\mathrm{T}},\textbf{v}_h) + b^{\mathrm{T}}(\textbf{v}_h, r_h^{\mathrm{T}}) &=& 0 \quad &\forall\, \textbf{v}_h \in \textbf{V}_{h,0}(T), \\ b^{\mathrm{T}}(\textbf{w}_h^{\mathrm{T}}, q_h) & = & \int_T \text{div}~{(\textbf{u} - {\it {\Pi}}_{F}^1 \textbf{u})} q_h ~\mathrm{d} x \quad & \forall\, q_h \in Q_h(T), \end{array} \end{align*} where $$a^{\mathrm{T}}(\cdot,\cdot)$$ and $$b^{\mathrm{T}}(\cdot,\cdot)$$ are the restrictions of the bilinear forms (2.1) and (2.2) on each triangle $$T \in \mathscr{T}$$. Note that this problem is solvable as $$\text{div}~{(\textbf{u} - {\it {\Pi}}_{F}^1 \textbf{u})} \perp \mathscr{P}^0(\mathscr{T})$$. Using Corollary 3.3 and the stability result of saddle point problems, we get   \begin{align*} \| {\textbf{w}^{\mathrm{T}}_h} \|_{1_h} \preccurlyeq \| \text{div}~{(\textbf{u} - {\it {\Pi}}_{F}^1 \textbf{u})} \|_0 \le c_{2} \| \textbf{u} \|_1, \end{align*} where the constant $$c_{2}$$ is independent of $$h$$ and $$k$$. Now we set $${\it {\Pi}}_{F}^2 \textbf{u} := \sum\limits_{T \in \mathscr{T}} \textbf{w}^{\mathrm{T}}_h$$ and see that $${\it {\Pi}}_{F} = {\it {\Pi}}_{F}^1 + {\it {\Pi}}_{F}^2$$ fulfils (3.4) as for all $$T \in \mathscr{T}$$ it holds that   \begin{align*} \text{div}~{{\it {\Pi}}_{F} \textbf{u}}|_T &= \text{div}~{{\it {\Pi}}_{F}^1 \textbf{u}}|_T + \text{div}~{{\it {\Pi}}_{F}^2 \textbf{u}}|_T = \text{div}~{{\it {\Pi}}_{F}^1 \textbf{u}}|_T + \sum\limits_{T \in \mathscr{T}} {\it {\Pi}}^{Q_h} (\text{div}~{\textbf{u}}|_T - \text{div}~{{\it {\Pi}}_{F}^1 \textbf{u}}|_T) = {\it {\Pi}}^{Q_h} \text{div}~{\textbf{u}}|_T, \\ \| {{\it {\Pi}}_{F} \textbf{u} } \|_{1_h} &\le \| {{\it {\Pi}}_{F}^1 \textbf{u} } \|_{1_h} + \| {{\it {\Pi}}_{F}^2 \textbf{u} } \|_{1_h} \preccurlyeq \| \textbf{u} \|_1 + \sum\limits_{T \in \mathscr{T}}\| {\textbf{w}^{\mathrm{T}}_h} \|_{1_h} \preccurlyeq \| \textbf{u} \|_1.\quad\square \end{align*} Remark 3.4 In the sense of Remark 2.4 one can choose the polynomial order in Theorem 3.1 as $$k_T$$; thus to show Theorem 2.2 one uses Theorem 3.2 and Corollary 3.3 with the proper order. Remark 3.5 As mentioned in Boffi et al. (2013, Remark 5.2.7), the existence of a Fortin operator $${\it {\Pi}}_{F}$$ can be used to construct an error estimation independent of the LBB constant $$\beta_\infty$$ of the infinite-dimensional Stokes problem (1.1). This may be an advantage for problems where $$\beta_\infty$$ is large. 4. Numerical examples Theorem 2.2 shows that the LBB constant is independent of the polynomial order $$k$$. Besides the analysis in Section 3, numerical tests also show this independence. In Table 1, one sees the different values of $$\beta$$ on the reference triangle $$\widehat{T}$$ and the reference quadrilateral $$\widehat{Q} = (0,1) \times (0,1)$$ for different polynomial orders $$k$$. Where the first line supports our results, the second line predicts that the polynomial robustness seems also to hold true in the case of quadrilaterals. The LBB constant was computed as the square root of the minimal positive eigenvalue of $$B G^{-1} B^{\mathrm{T}} M^{-1}$$, where $$B$$ is the discrete matrix of $$b(\cdot, \cdot)$$ and $$G$$ and $$M$$ are the Gramian matrices generated by $$(\textbf{V}_h, \| {\cdot} \|_{1_h})$$ and $$(Q_h, \| \cdot \|_0),$$ respectively. Table 1 LBB constant dependence of $$k$$ on the reference triangle $$\widehat{T}$$ and the reference quadrilateral $$\widehat{Q}$$ k  4  8  16  32  Triangle  0.167  0.190  0.201  0.205  Quadrilateral  0.305  0.313  0.315  0.315  k  4  8  16  32  Triangle  0.167  0.190  0.201  0.205  Quadrilateral  0.305  0.313  0.315  0.315  Table 1 LBB constant dependence of $$k$$ on the reference triangle $$\widehat{T}$$ and the reference quadrilateral $$\widehat{Q}$$ k  4  8  16  32  Triangle  0.167  0.190  0.201  0.205  Quadrilateral  0.305  0.313  0.315  0.315  k  4  8  16  32  Triangle  0.167  0.190  0.201  0.205  Quadrilateral  0.305  0.313  0.315  0.315  In the second example, we take a closer look at the stability of the method introduced in Section 2. For this, we solve problem (2.3) on the unit square $${\it {\Omega}}=(0,1)^2$$, where the exact solution is given by $$\textbf{u}:=\text{curl}~(\zeta) = (\zeta_{,y}, -\zeta_{,x})$$ with $$\zeta = \sin(x)^2\cos(x)^2$$, $$p=0$$, and set $$\nu=1$$. For this purpose, we use a triangulation with $$52$$ elements with $$h \approx 0.2$$ and a stabilization parameter $$\alpha = 4$$. In Table 2, we see the behaviour of the error $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$, the error of the best approximation $$\| {\textbf{u}-\textbf{u}_{\mathrm{best}}} \|_{1_h}$$ with $$\textbf{u}_{\mathrm{best}} := \arg \min\limits_{\textbf{v}_h \in \textbf{V}_h} \| {\textbf{u}-\textbf{v}_h} \|_{1_h}$$ and the ratio $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h} / \| {\textbf{u}-\textbf{u}_{\mathrm{best}}} \|_{1_h}$$. The values support that the discrete Stokes solution is close to the best approximation as the ratio is practically close to 1. Table 2 The error $$\| {\textbf{u}-\textbf{u}_{\mathrm{best}}} \|_{1_h}$$, $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$ and the ratio for different $$k$$ $$k$$  $$2$$  $$4$$  $$6$$  $$8$$  $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$  $$1.587$$  $$3.343e-02$$  $$3.430e-04$$  $$2.258e-06$$  $$\| {\textbf{u}-\textbf{u}_{\rm{best}}} \|_{1_h}$$  $$1.180$$  $$2.302e-02$$  $$2.309e-04$$  $$1.597e-06$$  Ratio  $$1.344$$  $$1.452$$  $$1.486$$  $$1.414$$  $$k$$  $$2$$  $$4$$  $$6$$  $$8$$  $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$  $$1.587$$  $$3.343e-02$$  $$3.430e-04$$  $$2.258e-06$$  $$\| {\textbf{u}-\textbf{u}_{\rm{best}}} \|_{1_h}$$  $$1.180$$  $$2.302e-02$$  $$2.309e-04$$  $$1.597e-06$$  Ratio  $$1.344$$  $$1.452$$  $$1.486$$  $$1.414$$  Table 2 The error $$\| {\textbf{u}-\textbf{u}_{\mathrm{best}}} \|_{1_h}$$, $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$ and the ratio for different $$k$$ $$k$$  $$2$$  $$4$$  $$6$$  $$8$$  $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$  $$1.587$$  $$3.343e-02$$  $$3.430e-04$$  $$2.258e-06$$  $$\| {\textbf{u}-\textbf{u}_{\rm{best}}} \|_{1_h}$$  $$1.180$$  $$2.302e-02$$  $$2.309e-04$$  $$1.597e-06$$  Ratio  $$1.344$$  $$1.452$$  $$1.486$$  $$1.414$$  $$k$$  $$2$$  $$4$$  $$6$$  $$8$$  $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$  $$1.587$$  $$3.343e-02$$  $$3.430e-04$$  $$2.258e-06$$  $$\| {\textbf{u}-\textbf{u}_{\rm{best}}} \|_{1_h}$$  $$1.180$$  $$2.302e-02$$  $$2.309e-04$$  $$1.597e-06$$  Ratio  $$1.344$$  $$1.452$$  $$1.486$$  $$1.414$$  All implementations of the numerical examples were performed with the finite element library NGSolve/Netgen (see Schöberl, 1997; 2014). 5. $$\boldsymbol{H^2}$$-stable polynomial extension In this section, we prove the existence of a polynomial-preserving $$H^2$$-stable extension operator. Note that in the two-dimensional case, the $$\text{curl}~$$ operator can be represented as a rotated gradient; thus on the boundary we have   \begin{align*} \text{curl}~{\textbf{u}_h} \cdot \textbf{n} = \nabla^\perp \textbf{u}_h \cdot \textbf{n} = \nabla \textbf{u}_h \cdot \tau \quad \textrm{on } \partial T. \end{align*} For the ease of notation and readability, we switch to the gradient and change the tangential and normal vector in Theorem 3.1. Furthermore, we skip the subscript $$h$$ of $$\textbf{u}_h$$; thus we write $$\textbf{u}$$ for a vectorial polynomial and show instead the following theorem. Theorem 5.1 ($$H^2$$-extension) For every $$k$$ there exists an operator $$\mathscr{E} : \pmb{\mathscr{P}}_{\tau}^k(\widehat{T}) \rightarrow \mathscr{P}^{k+1}(\widehat{T})$$, such that for $$\textbf{u} \in \pmb{\mathscr{P}}_{\tau}^k(\widehat{T})$$ it holds that   \begin{align} &\nabla {\mathscr{E}(\textbf{u})} \cdot \tau = \textbf{u} \cdot \tau \quad \textrm{on } \partial \widehat{T} \label{ext:h2ext:propA}, \\ \end{align} (5.1)  \begin{align} & \|(\textbf{u} - \nabla \mathscr{E}(\textbf{u})) \cdot \textbf{n} \|_{0,\partial \widehat{T}} \preccurlyeq \frac{1}{k} \| u \|_{1, \widehat{T}} \label{ext:h2ext:propB}, \\ \end{align} (5.2)  \begin{align} & \| \mathscr{E}(\textbf{u}) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}. \label{ext:h2ext:propC} \end{align} (5.3) Proof. The proof is provided in Section 5.5. □ 5.1 Literature and structure of this section Extension or lifting operators are a well-discussed topic, as they present the inverse map of trace operators that are known to be continuous, for example, in the scalar case from $$H^1(T)$$ onto $$H^{1/2}(\partial T)$$. The challenge then is to construct an operator that maps functions from $$H^{1/2}(\partial T)$$ into $$H^1(T)$$. Furthermore, a polynomial extension has the additional property that the operator maps a given polynomial on the boundary onto a polynomial on the element. The importance of such extensions has its origin mostly in the analysis of $$p$$- and $$hp$$- finite element methods (see, for example, Demkowicz & Babuška, 2003) and spectral methods and preconditioning (see, for example, Schöberl et al., 2008). Although the existence of polynomial $$H^1$$-, $$H(\operatorname{div})$$- and $$H(\text{curl}~op)$$-extensions are already well analysed, to the best of our knowledge, a stable polynomial $$H^2$$-extension, which is presented in this section, is the first result of this kind. Besides the application of the $$H^2$$-stable extension operator in the proof of the local LBB condition Theorem 3.2, this operator may be of interest, for example, to construct proper interpolation operators for fourth-order $$C^0$$-continuous DG methods (see, for example, Brenner et al., 2010, 2012). Before we present our results, we want to mention some literature, as many techniques we use are motivated by their accomplishments. The first work we consider is the pioneering article of Babuška & Suri (1987), which contains major ideas that are developed and adapted in later contributions including the technique of splitting the operator into primary and proper correction liftings. For polynomial-preserving extensions, we want to mention the work of Maday (1989) and Bernardi & Maday (1990) for two dimensions and the work of Belgacem (1994), Muñoz-Sola (1997) and, more recently, Ainsworth & Demkowicz (2009) for three dimensions. Finally, we want to mention the works of Demkowicz et al. (2008, 2009, 2012), where the polynomial extensions also provide a commuting diagram property for the corresponding spaces. For the proof of Theorem 5.1, we proceed in several steps. We start with an extension operator for the tangential values in Section 5.2 to provide the properties (5.1) and (5.2). After that we show in Section 5.3 how to proceed with the normal component. For this purpose, we assume that the polynomial on the edge has a zero of order 2 in the vertices to show an estimate in the $$H^{1/2}_{00}$$-norm. To lift the normal trace also for arbitrary polynomials, we show in Section 5.4 that the error, hence the part of the polynomial that does not satisfy the assumptions needed for the extension before, is bounded with a proper order of $$k$$ to show (5.2). Finally, we combine all estimates in Section 5.5 to prove Theorem 5.1. For the continuity estimates, we refer several times to Lemma 5.3 which is proved below. Remark 5.2 In this section, we only present the major techniques to show the $$H^2$$-continuity for the most difficult terms due to the similar structure of the rest. For a more detailed analysis, we refer to Lederer (2016, Chapter 6). Lemma 5.3 Let $$a \in C^1(E_1)$$ and $$u \in H^{1/2}(E_1)$$. For $$(x,y) \in \widehat{T}$$ we define $$g(x,y):=\int_0^1 a(s) u(x+sy) ~\mathrm{d} s$$. There holds   \begin{align*} \| g\|_{1, \widehat{T}} \preccurlyeq \| u \|_{1/2,E_1}. \end{align*} Proof. To derive the estimate of the $$L^2$$-norm we define for a fixed $$y \in [0,1]$$ the line $${l_y := \{(x,y): x \in [0,1-y] \}}$$ and use the Cauchy–Schwarz inequality to get   \begin{align*} \| g \|_{0,l_y}^2 &= \int_0^{1-y}\left(\int_0^1 a(s) u(x+sy) ~\mathrm{d} s \right)^2~\mathrm{d} x \preccurlyeq \int_0^{1-y} \int_0^1 u(x + sy)^2 ~\mathrm{d} s ~\mathrm{d} x. \end{align*} With the substitution $$t = x+sy$$ and Fubini’s theorem this leads to   \begin{align*} \| g \|_{0,\widehat{T}}^2 &=\int_0^1 \| g \|_{0,l_y}^2 ~\mathrm{d} y \preccurlyeq \int_0^1 \frac{1}{y} \int_0^{1-y} \int_x^{x+y} u(t)^2 ~\mathrm{d} t ~\mathrm{d} x ~\mathrm{d} y = \int_0^1 \frac{1}{y}\iint\limits_{\substack{0 \le x \le 1-y \\ x \le t \le x+y }} u(t)^2 ~\mathrm{d}(x,t) ~\mathrm{d} y \\ &\preccurlyeq \int_0^1 \frac{1}{y}\iint\limits_{\substack{0 \le t \le 1 \\ t-y \le x \le t }} u(t)^2 ~\mathrm{d}(x,t) ~\mathrm{d} y = \int_0^1 \frac{1}{y} \int_0^{1} u(t)^2 \underbrace{\int_{t-y}^{\mathrm{t}} ~\mathrm{d} x}_{=y} ~\mathrm{d} t ~\mathrm{d} y \preccurlyeq \int_0^1 \| u \|_{0,E_1}^2 ~\mathrm{d} y \preccurlyeq \| u \|_{1/2,E_1}^2. \end{align*} To bound the first-order derivatives, we use the real method of interpolation of spaces (see Bergh & Löfström, 1976) and Peetre’s K-functional technique (see Peetre, 1963). It is well known that we have an equivalent norm on the space $$H^{1/2}(E_1)$$ given by   \begin{align*} \| u \|_{1/2,E_1}^2 = \int_0^\infty y^{-2} |K(y,u)|^2 ~\mathrm{d} y \quad \textrm{with} \quad K(y,u) := \inf\limits_{\substack{u_0, u_1 \\u = u_0 + u_1 } } \sqrt{\|u_0 \|_{0,E_1}^2 + y^2 \|u_1'\|_{0,E_1}^2}. \end{align*} The derivative with respect to $$x$$ reads, $$g_{,x} = \int_0^1 u'(x+sy) a(s) ~\mathrm{d} s$$; thus with a similar estimate as for the $$L^2$$-norm we get $$\| g_{,x} \|_{0,l_y}^2 \preccurlyeq \| u'\|_{0,E_1}^2$$. Using integration by parts, we can also write   \begin{align*} g_{,x} = \underbrace{\frac{1}{y} \int_0^1 a'(s) u(x+sy) ~\mathrm{d} s}_{=:A_1} + \underbrace{\frac{1}{y}(a(1) u(x+y) - a(0)u(x))}_{=:B_2}. \end{align*} For $$A_1$$ we proceed similarly as for the $$L^2$$ estimate to get $$\| A_1 \|_{0,l_y}^2 \preccurlyeq \frac{1}{y^2} \| u \|^2_{0,E_1}$$. For $$B_1$$ we get with the substitution $$t= x+y$$,   \begin{align*} \| B_1 \|_{0,l_y}^2 &= \frac{1}{y^2} \int_0^{1-y} (a(1)u(x+y) - a(0)u(x))^2 ~\mathrm{d} x \preccurlyeq \frac{1}{y^2} \int_0^{1-y} u(x+y)^2 + u(x)^2 ~\mathrm{d} x \\ &\preccurlyeq \frac{1}{y^2} \int_y^{1} u(t)^2 ~\mathrm{d} t + \frac{1}{y^2}\int_0^{1-y} u(x)^2 ~\mathrm{d} x \preccurlyeq \frac{1}{y^2} \| u \|^2_{0,E_1}. \end{align*} Combining both estimates, so $$\| g_{,x} \|_{0,l_y} \preccurlyeq 1/y \| u\|_{0,E_1}$$ and $$\| g_{,x} \|_{0,l_y} \preccurlyeq \| u'\|_{0,E_1}$$, leads to   \begin{align*} \| g_{,x} \|^2_{0,l_y} &\le \inf\limits_{\substack{u_0, u_1 \\u = u_0 + u_1 } } \| u_0 \|^2_{0,l_y} + \| u_1 \|^2_{0,l_y} \preccurlyeq \inf\limits_{\substack{u_0, u_1 \\u = u_0 + u_1 } } \frac{1}{y^2} \| u_0 \|^2_{0,E_1} + \| u_1' \|^2_{0,E_1} \preccurlyeq \frac{1}{y^2} K(y,u)^2, \end{align*} and thus   \begin{align*} \| g_{,x} \|_{0,\widehat{T}}^2 &= \int_0^1 \| g \|_{0,l_y}^2 ~\mathrm{d} y \preccurlyeq \int_0^1 \frac{1}{y^2} K(y,u)^2 ~\mathrm{d} y \preccurlyeq \int_0^\infty \frac{1}{y^2} K(y,u)^2 ~\mathrm{d} y = \|u \|_{1/2,E_1}^2. \end{align*} A similar estimation for the derivative with respect to $$y$$ concludes the proof. □ 5.2 Tangential extension Theorem 5.4 For every $$k$$ there exists an operator $$\mathscr{E}^{\tau} : \pmb{\mathscr{P}}_{\tau}^k(\widehat{T}) \rightarrow \mathscr{P}^{k+1}(\widehat{T})$$ such that for $$\textbf{u} \in \pmb{\mathscr{P}}_{\tau}^k(\widehat{T})$$ it holds that   \begin{align} &\nabla {\mathscr{E}^{\tau}(\textbf{u})} \cdot \tau = \textbf{u} \cdot \tau \quad \textrm{on } \partial \widehat{T}, \label{ext:h2exttang:propA} \\ \end{align} (5.4)  \begin{align} & \| \mathscr{E}^{\tau}(\textbf{u}) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}. \label{ext:h2exttang:propC} \end{align} (5.5) Proof. For the proof we proceed in several steps. First, we construct an extension from the lower edge $$E_1$$ onto the triangle $$\widehat{T}$$ without any restrictions on the values on the two other edges $$E_2$$ and $$E_3$$. After that, we construct two more extensions with proper values on the other edges and combine them afterwards to define $$\mathscr{E}^{\tau}$$. Step 1. For $$u_\tau(x) := \textbf{u}(x,0) \cdot \tau$$, where $$\tau := (1,0)^{\mathrm{t}}$$ is the tangential vector on the lower edge $$E_1$$, we define   \begin{align} \label{ext::tangextdefstepone} \psi(x) := \int_0^x u_\tau(s) ~\mathrm{d} s \qquad \textrm{and} \qquad \mathscr{E}^{\tau}_1(\textbf{u})(x,y) := \int_0^1 \psi(x + sy) ~\mathrm{d} s. \end{align} (5.6) Note that the derivatives read   \begin{align*} \mathscr{E}^{\tau}_{1,x}(\textbf{u}) = \int_0^1 \textbf{u}t(x+sy) ~\mathrm{d} s \quad \textrm{and} \quad \mathscr{E}^{\tau}_{1,y}(\textbf{u}) = \int_0^1 \textbf{u}t(x+sy)s ~\mathrm{d} s, \end{align*} thus   \begin{align*} \nabla \mathscr{E}^{\tau}_1(\textbf{u}) \cdot \tau = \int_0^1 \textbf{u}t(x) ~\mathrm{d} s = \textbf{u} \cdot \tau \quad \textrm{on } E_1. \end{align*} Next we show the $$H^2$$-continuity, so $$\| \mathscr{E}^{\tau}_1(\textbf{u}) \|_{0,\widehat{T}}^2 + \| \nabla \mathscr{E}^{\tau}_1(\textbf{u}) \|_{1,\widehat{T}}^2 \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}^2$$. For the estimate of the $$L^2$$-norm, we proceed similarly to the first steps of the proof of Lemma 5.3, thus we estimate the $$L^2$$-norm on a line $$l_y$$ using the Cauchy–Schwarz inequality and Fubini’s theorem and then integrate from $$0$$ to $$1$$ with respect to $$y$$. For the first-order derivatives, we can use Lemma 5.3 to get   \begin{align*} \| \mathscr{E}^{\tau}_{1,x}(\textbf{u}) \|_{1,\widehat{T}}^2 \preccurlyeq \| \textbf{u}t \|_{1/2,E_1}^2 \quad \text{and} \quad \| \mathscr{E}^{\tau}_{1,y}(\textbf{u}) \|_{1,\widehat{T}}^2 \preccurlyeq \| \textbf{u}t \|_{1/2,E_1}^2. \end{align*} Altogether we have   \begin{align} \label{ext::tangextsteponeest} \| \mathscr{E}^{\tau}_{1}(\textbf{u}) \|_{2,\widehat{T}}^2 \preccurlyeq \| \textbf{u}t \|_{1/2,E_1}^2 \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}^2. \end{align} (5.7) Step 2. In the next step, we want to construct an extension from the lower edge $$E_1$$ with the restriction that the gradient of the extension has zero tangential values on the edge $$E_2$$. Similarly to before, we define for $$u_\tau(x) := \textbf{u}(x,0) \cdot \tau$$ on $$E_1$$,   \begin{align*} \psi(x) := \int_0^x u_\tau(s) ~\mathrm{d} s - \overline{\psi(x)} \qquad \textrm{with} \qquad \overline\psi(x) := \int_0^1 u_\tau(s) ~\mathrm{d} s, \end{align*} and   \begin{align} \label{ext::tangextdefsteptwo} \mathscr{E}^{\tau}_2(\textbf{u})(x,y) := \int_0^1 \psi(x + sy) ~\mathrm{d} s - \frac{y}{1-x} \int_0^1 \psi(x+s(1-x)) ~\mathrm{d} s. \end{align} (5.8) Using integration by parts for the second integral we furthermore use the representation   \begin{align} \label{ext::tangextdefsteptwoalt} \mathscr{E}^{\tau}_2(\textbf{u})(x,y) &= \int_0^1 \psi(x + sy) ~\mathrm{d} s + y \int_0^1 \psi'(x+s(1-x)) s ~\mathrm{d} s - \overbrace{\psi(1)}^{=0} \nonumber\\ &= \int_0^1 \psi(x + sy) ~\mathrm{d} s + y \int_0^1 u_\tau(x+s(1-x)) s ~\mathrm{d} s. \end{align} (5.9) On the edges $$E_1, E_2$$ we have   \begin{align*} \mathscr{E}^{\tau}_2(\textbf{u})|_{E_2} = 0 \Rightarrow \nabla \mathscr{E}^{\tau}_2(\textbf{u}) \cdot \tau|_{E_2} = 0 \quad \textrm{and} \quad \mathscr{E}^{\tau}_2(\textbf{u}) \cdot \tau|_{E_1} = \frac{\partial \psi}{\partial x} = u_\tau = \textbf{u} \cdot \tau|_{E_1}. \end{align*} For the $$H^2$$-continuity we proceed as before. Note that the first part of $$\mathscr{E}^{\tau}_2(\textbf{u})$$ is the same as in Step 1, so we consider only the second integral from (5.9); thus the correction term is   \begin{align*} \mathscr{E}^{\tau,c}_2(\textbf{u}) := y \int_0^1 u_\tau(x+s(1-x)) s ~\mathrm{d} s. \end{align*} We present only the estimate for the second-order derivative with respect to $$x$$ as some new techniques have to be used there. The rest follows with similar steps to before. Using integration by parts for $$\textbf{u}t'$$ we observe   \begin{align*} \mathscr{E}^{\tau,c}_{2,xx}(\textbf{u}) &= \frac{y}{(1-x)^2} \int_0^1 \textbf{u}t(x+s(1-x))(2s-1) ~\mathrm{d} s \\ &\quad{} + \frac{y}{1-x} \int_0^1 \textbf{u}t'(x+s(1-x))(2s-1) ~\mathrm{d} s \\ &= \frac{y}{(1-x)^2} \int_0^1 \textbf{u}t(x+s(1-x))(2s-1) ~\mathrm{d} s \\ & \quad{}+ \frac{y}{(1-x)^2} \int_0^1 \textbf{u}t(x+s(1-x))(4s-3) ~\mathrm{d} s + \frac{y}{(1-x)^2} \textbf{u}t(x), \end{align*} and by splitting the second integral into two terms, finally   \begin{align*} \mathscr{E}^{\tau,c}_{2,xx}(\textbf{u}) = \underbrace{\frac{y}{(1-x)^2} \int_0^1 \textbf{u}t(x+s(1-x))(6s-3) ~\mathrm{d} s}_{=:A_2} + \underbrace{ \frac{y}{(1-x)^2} \int_0^1 \textbf{u}t(x) - \textbf{u}t(x+s(1-x)) ~\mathrm{d} s}_{=: B_2}. \end{align*} We start with the estimate of $$A_2$$. For this, note that the polynomial $$6s-3$$ has a zero integral value on $$E_1$$, so we subtract the mean value $$\overline{\textbf{u}t(x)} := \frac{1}{1-x} \int_x^1 \textbf{u}t(s) ~\mathrm{d} s$$, and get   \begin{align*} A_2 &= \frac{y}{(1-x)^2} \int_0^1 \textbf{u}t(x+s(1-x))(6s-3) ~\mathrm{d} s = \frac{y}{(1-x)^2} \int_0^1 (\textbf{u}t(x+s(1-x)) - \overline{\textbf{u}t(x)}) (6s-3) ~\mathrm{d} s. \end{align*} Using the Cauchy–Schwarz inequality and $$\frac{y}{1-x} \le 1$$ on $$\widehat{T}$$ and the substitution $$t = x + s(1-x)$$ leads to   \begin{align*} \| A_2 \|^2_{0,\widehat{T}} &= \int_0^1 \int_0^{1-x} A_2^2 ~\mathrm{d} y ~\mathrm{d} x \\ & \preccurlyeq \int_0^1 \frac{1}{(1-x)^2} \left(\int_0^1 (\textbf{u}t(x+s(1-x)) - \overline{\textbf{u}t(x)}) (6s-3) ~\mathrm{d} s \right)^2 \int_0^{1-x} ~\mathrm{d} y ~\mathrm{d} x\\ & \preccurlyeq \int_0^1 \frac{1}{1-x} \int_0^1 (\textbf{u}t(x+s(1-x)) - \overline{\textbf{u}t(x)})^2 ~\mathrm{d} s ~\mathrm{d} x \\ & = \int_0^1 \frac{1}{(1-x)^2} \int_x^1 (\textbf{u}t(t) - \overline{\textbf{u}t(x)})^2 ~\mathrm{d} t ~\mathrm{d} x. \end{align*} For the next step, we use the following identity for the inner integral:   \begin{align*} \int_x^1 (\textbf{u}t(t) - \overline{\textbf{u}t(x)})^2 ~\mathrm{d} t = \frac{1}{2(1-x)} \int_x^1 \int_x^1 (\textbf{u}t(t) - \textbf{u}t(s))^2 ~\mathrm{d} t ~\mathrm{d} s. \end{align*} Similarly to Step 1, we use Fubini’s theorem to handle the $$(1-x)^3$$ in the denominator, so   \begin{align*} \| A_2 \|^2_{0,\widehat{T}} &\preccurlyeq \int_0^1 \frac{1}{(1-x)^3} \int_x^1 \int_x^1 (\textbf{u}t(t) - \textbf{u}t(s))^2 ~\mathrm{d} t ~\mathrm{d} s ~\mathrm{d} x \\ & \preccurlyeq \iiint\limits_{\substack{x \le s\\ x \le t \\ 0 \le s,t \le 1 }} \frac{(\textbf{u}t(t) - \textbf{u}t(s))^2 }{(1-x)^3} ~\mathrm{d} (s,t,x) \preccurlyeq \int_0^1 \int_0^1 \int_0^{\min(s,t)} \frac{1}{(1-x)^3} ~\mathrm{d} x (\textbf{u}t(t) - \textbf{u}t(s))^2 ~\mathrm{d} t ~\mathrm{d} s. \end{align*} Without loss of generality assuming $$s < t$$ and using $$1-s \ge t-s$$ for $$s < t < 1$$, it follows for the inner integral that   \begin{align*} \int_0^{\min(s,t)} \frac{1}{(1-x)^3}~\mathrm{d} x = \left(\frac{1}{1-\min(s,t) } \right) ^2 - 1 \le \left(\frac{1}{1-s} \right)^2 \le \left (\frac{1}{t-s} \right)^2. \end{align*} Together with the definition of the $$H^{1/2}$$-seminorm (see, Sobolev–Slobodeckij norm, for example, in Schwab, 1998, Theorem A.7) we get   \begin{align*} \| A_2 \|^2_{0,\widehat{T}} & \preccurlyeq \int_0^1 \int_0^1 \frac{(\textbf{u}t(t) - \textbf{u}t(s))^2}{(t-s)^2} ~\mathrm{d} t ~\mathrm{d} s = | \textbf{u}t|^2_{1/2,E_1} \preccurlyeq \| \textbf{u} \|^2_{1,\widehat{T}}. \end{align*} For $$B_2$$, we proceed similarly using the Cauchy–Schwarz inequality, Fubini’s theorem, the substitution $$t = x + s(1-x)$$ and $$1-x \ge t-x$$ for $$x < t < 1$$; thus   \begin{align*} \| B_2 \|^2_{0,\widehat{T}} &= \int_0^1 \int_0^{1-x} B_2^2 ~\mathrm{d} y ~\mathrm{d} x \preccurlyeq \int_0^1 \frac{1}{1-x} \int_0^1 \left(\textbf{u}t(x) - \textbf{u}t(x+s(1-x)) \right)^2 ~\mathrm{d} s ~\mathrm{d} x \\ & = \int_0^1 \frac{1}{(1-x)^2} \int_x^1 \left(\textbf{u}t(x) - \textbf{u}t(t) \right)^2 ~\mathrm{d} t ~\mathrm{d} x = \iint\limits_{\substack{x \le t \le 1\\0\le x \le 1 }} \frac{\left(\textbf{u}t(x) - \textbf{u}t(t) \right)^2}{(1-x)^2} ~\mathrm{d}(x,t)\\ &\preccurlyeq \iint\limits_{\substack{0 \le x \le t \\0\le t \le 1 }} \frac{\left(\textbf{u}t(x) - \textbf{u}t(t) \right)^2}{(x-t)^2} ~\mathrm{d}(x,t) \preccurlyeq \int_0^1 \int_0^1 \frac{ \left(\textbf{u}t(x) - \textbf{u}t(t) \right)^2}{(x-t)^2} ~\mathrm{d} x ~\mathrm{d} t \preccurlyeq | \textbf{u}t|^2_{1/2,E_1} \preccurlyeq \| \textbf{u} \|^2_{1,\widehat{T}}. \end{align*} Altogether we have $$\| \mathscr{E}^{\tau,c}_{2,xx}(\textbf{u})\|_{0,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}$$ and assuming similar estimates for the other derivatives, finally   \begin{align} \label{ext::exttangtwoest} \| \mathscr{E}^{\tau}_{2}(\textbf{u})\|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}. \end{align} (5.10) Step 3. For this step, we assume that the tangential values of the input function $$\textbf{u}$$ are zero on the edges $$E_2$$ and $$E_3$$ and that it has a zero tangential integral value; thus   \begin{align} \label{ext::assumptionsstepthree} \int_{\partial \widehat{T}} \textbf{u} \cdot \tau = 0 \quad \textrm{and} \quad \textbf{u} \cdot \tau|_{E_2} = \textbf{u} \cdot \tau|_{E_3} = 0. \end{align} (5.11) We set $$u_\tau(x) := \textbf{u}(x,0) \cdot \tau$$ on $$E_1$$ and $$\psi(x) := \int_0^x u_\tau(s) ~\mathrm{d} s$$, to define   \begin{align} \label{ext::tangextdefstepthree} \mathscr{E}^{\tau}_{3}(\textbf{u}) &:= \int_0^1 \psi(x+sy) ~\mathrm{d} s - \frac{y}{1-x} \int_0^1 \psi(x+s(1-x)) ~\mathrm{d} s \nonumber\\ &\quad\ \underbrace{- \frac{y}{x+y} \int_0^1 \psi(s(x+y)) ~\mathrm{d} s}_{=: A_3} + \underbrace{ y \int_0^1 \psi(s) ~\mathrm{d} s}_{=: B_3}. \end{align} (5.12) As in Step 2, we observe $$\nabla \mathscr{E}^{\tau}_{3}(\textbf{u}) \cdot \tau |_{E_2} = \nabla \mathscr{E}^{\tau}_{3}(\textbf{u}) \cdot \tau |_{E_3} = 0$$ and $$\nabla \mathscr{E}^{\tau}_{3}(\textbf{u}) \cdot \tau |_{E_1} = \textbf{u}t$$. For the $$H^2$$-continuity we have only to estimate the terms $$A_3$$ and $$B_3$$ as the other terms are the same as in Steps 1 and 2. For this purpose, note that due to the assumptions on $$\textbf{u}$$, we have $$\psi(1) = \psi(0) = 0$$. Using the identity $$\frac{~\mathrm{d}}{~\mathrm{d} s} \psi(s(x+y)) \frac{1}{x+y} = \psi'(s(x+y))$$ and integration by parts we write   \begin{align} \label{ext::tangextdefstepthreealtone} A_3 &= - \frac{y}{x+y} \int_0^1 \psi(s(x+y)) ~\mathrm{d} s = y \int_0^1 \psi'(s(x+y))(s-1) ~\mathrm{d} s = y \int_0^1 \textbf{u}t(s(x+y))(s-1) ~\mathrm{d} s, \end{align} (5.13) and   \begin{align} \label{ext::tangextdefstepthreealttwo} B_3 = y \int_0^1 \psi(s) ~\mathrm{d} s = -y \int_0^1 \psi'(s)s ~\mathrm{d} s = -y \int_0^1 \textbf{u}t(s)s ~\mathrm{d} s. \end{align} (5.14) Using these representations and the same techniques as in Steps 2 and 3, we estimate the $$H^2$$-norm of $$A_3$$ and $$B_3$$ to show   \begin{align} \label{ext::exttangthreeest} \| \mathscr{E}^{\tau}_{3}(\textbf{u}) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}. \end{align} (5.15) Step 4. We finally combine the three extensions to show Theorem 5.4. For this purpose, let us assume we have a given function $$\textbf{u} \in \pmb{\mathscr{P}}_{\tau}^k(\widehat{T})$$ with $$\int_{\partial \widehat{T}} \textbf{u} \cdot \tau = 0$$. We first introduce two mappings from the reference triangle $$\widehat{T}$$ to itself by   \begin{align*} F_2: (x,y) \mapsto (x,1-x-y) \quad \textrm{and} \quad F_3: (x,y) \mapsto (y,x), \end{align*} where $$F_2$$ maps the values from $$E_2$$ to $$E_1$$ and vice versa, and the mapping $$F_3$$ from $$E_3$$ to $$E_1$$ and vice versa. Furthermore, we define for $$F_2$$ and $$F_3$$ the corresponding covariant mappings $$\mathscr{C}_2$$ and $$\mathscr{C}_3$$. Using those transformations we now introduce the extensions from Steps 2 and 3 also from the other edges; thus we define   \begin{align*} \tilde{\mathscr{E}^{\tau}_2}(\textbf{u})(x,y) = \mathscr{E}^{\tau}_2(\mathscr{C}_2\textbf{u})(F_2(x,y)) \quad \textrm{and} \quad \tilde{\mathscr{E}^{\tau}_3}(\textbf{u})(x,y) = \mathscr{E}^{\tau}_3(\mathscr{C}_3\textbf{u})(F_3(x,y)), \end{align*} with the properties   \begin{align} &(\nabla \tilde{\mathscr{E}^{\tau}_2}(\textbf{u})\cdot \tau)|_{E_2} = (\textbf{u} \cdot \tau)|_{E_2}, \quad (\nabla \tilde{\mathscr{E}^{\tau}_2}(\textbf{u})\cdot \tau)|_{E_1} = 0, \quad \| \tilde{\mathscr{E}^{\tau}_2}(\textbf{u}) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}, \label{ext::exttangtwotildeest}\\ \end{align} (5.16)  \begin{align} &(\nabla \tilde{\mathscr{E}^{\tau}_3}(\textbf{u})\cdot \tau)|_{E_3} = (\textbf{u} \cdot \tau)|_{E_3}, \quad (\nabla \tilde{\mathscr{E}^{\tau}_3}(\textbf{u})\cdot \tau)|_{E_1} = (\nabla \tilde{\mathscr{E}^{\tau}_3}(\textbf{u})\cdot \tau)|_{E_2}= 0, \quad \| \tilde{\mathscr{E}^{\tau}_3}(\textbf{u}) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}},\label{ext::exttangthreetildeest} \end{align} (5.17) which follow from the proper transformation of tangential values due to the use of the covariant transformations $$\mathscr{C}_2$$ and $$\mathscr{C}_3$$ and estimates (5.10) and (5.15). We define the final extension by setting $$e_1 := \mathscr{E}^{\tau}_1(\textbf{u})$$, $$e_2 := e_1 + \tilde{\mathscr{E}^{\tau}_2}(\textbf{u}- \nabla e_1)$$ and $$\mathscr{E}^{\tau}(\textbf{u}) := e_2 + \tilde{\mathscr{E}^{\tau}_3}(\textbf{u} - \nabla e_2)$$. Note that due to Green’s theorem, the surface integral over the boundary of the reference element $$\partial \widehat{T}$$ of $$(\textbf{u} - \nabla e_2) \cdot \tau$$ is equal to zero, and due to the properties of $$\tilde{\mathscr{E}^{\tau}_2}$$ and $$\mathscr{E}^{\tau}_1$$, the tangential values on $$E_1$$ and $$E_2$$ also vanish; thus assumptions (5.11) are fulfilled. Using (5.16) and (5.17) we observe   \begin{align*} \nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \tau|_{E_1} &= \nabla e_2 \cdot \tau|_{E_1} + \underbrace{\nabla \tilde{\mathscr{E}^{\tau}}_3(\textbf{u} - \nabla e_2) \cdot \tau|_{E_1}}_{=0} = \nabla e_1 \cdot \tau|_{E_1} + \underbrace{ \nabla \tilde{\mathscr{E}^{\tau}}_2(\textbf{u} - \nabla e_1) \cdot \tau|_{E_1}}_{=0} = (\textbf{u} \cdot \tau)|_{E_1}, \end{align*} and similarly $$\nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \tau|_{E_2} = \textbf{u} \cdot \tau|_{E_2}$$ and $$\nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \tau|_{E_3} = \textbf{u} \cdot \tau|_{E_3}$$; thus property (5.4) is fulfilled. With (5.16), (5.17) and (5.7) and the linearity of the operators, we furthermore have $$\| \mathscr{E}^{\tau}(\textbf{u}) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}$$; thus also the $$H^2$$-continuity (5.5) holds true. It remains to show that $$\mathscr{E}^{\tau}(\textbf{u}) \in \mathscr{P}^{k+1}(\widehat{T})$$. First, note that from $$\textbf{u} \in [\mathscr{P}^k(\widehat{T})]^2$$, it follows that $$\textbf{u} \cdot \tau \in \mathscr{P}^k(\partial \widehat{T})$$. Looking at the definition of the first extension (5.6) we see that we integrate $$u$$ from $$0$$ to $$x$$ to define $$\psi$$; thus, here, we increase the order by 1 resulting in $$\mathscr{E}^{\tau}_1(\textbf{u}) \in \mathscr{P}^{k+1}(\widehat{T})$$. For the other two extensions (5.8) and (5.12), this may not hold due to the fractional factors, but using the alternative representations (5.9), (5.13) and (5.14) we see that also the corrections are polynomial liftings and it follows that $$\mathscr{E}^{\tau}(\textbf{u}) \in \mathscr{P}^{k+1}(\widehat{T})$$. □ 5.3 Normal extension Theorem 5.5 For every $$k$$ there exists an operator $$\mathscr{E}^{\textbf{n}} : \mathscr{P}_{00}^k(E_1) \rightarrow \mathscr{P}^{k+1}(\widehat{T})$$ such that for $$u \in \mathscr{P}_{00}^k(E_1)$$ it holds that   \begin{align} &\mathscr{E}^{\textbf{n}}(u) = 0 \quad\hspace{20mm} \textrm{on } \partial \widehat{T}, \label{ext:h2extnorm:propA} \\ \end{align} (5.18)  \begin{align} &\nabla \mathscr{E}^{\textbf{n}}(u) \cdot \textbf{n} = u \quad \textrm{on } E_1, \label{ext:h2extnorm:propB} \\ \end{align} (5.19)  \begin{align} &\nabla \mathscr{E}^{\textbf{n}}(u) \cdot \textbf{n} = 0 \quad \textrm{on } \partial \widehat{T} \setminus E_1, \label{ext:h2extnorm:propD} \\ \end{align} (5.20)  \begin{align} & \| \mathscr{E}^{\textbf{n}}(u) \|_{2,\widehat{T}} \preccurlyeq \| u \|_{1/2^*,E _1}. \label{ext:h2extnorm:propC} \end{align} (5.21) Proof. Similarly to the proof of Theorem 5.4 we proceed in several steps. We first construct an extension from the lower edge and correct the value and the normal derivative on the other two edges afterwards. Step 1. We start with the first extension, so we define   \begin{align} \label{ext::extnormstepone} \mathscr{E}^{\textbf{n}}_1(u) := -y \int_0^1a(s)u(x+sy) ~\mathrm{d} s \quad \textrm{with} \quad a(s):=6s(1-s). \end{align} (5.22) It immediately follows that $$\mathscr{E}^{\textbf{n}}_1(u)|_{E_1} = 0$$. For the derivatives, we observe due to $$a(0) = a(1) = 0$$, $$\int_0^1a(s) ~\mathrm{d} s = 1$$, $$\int_0^1a'(s)s ~\mathrm{d} s = -1$$ and using integration by parts with $$u'(x+sy) = \frac{1}{y}\frac{d}{ds} u(x+sy)$$, that   \begin{align*} \mathscr{E}^{\textbf{n}}_{1,x}(u) &= -y \int_0^1 a(s) u'(x+sy) ~\mathrm{d} s = \int_0^1 a'(s) u(x+sy) ~\mathrm{d} s, \\ \mathscr{E}^{\textbf{n}}_{1,y}(u) &= -\int_0^1 a(s) u(x+sy) ~\mathrm{d} s - y\int_0^1 a(s) u'(x+sy)s ~\mathrm{d} s = \int_0^1 a'(s) u(x+sy)s ~\mathrm{d} s; \end{align*} thus   \begin{align*} \nabla \mathscr{E}^{\textbf{n}}_1(u) \cdot n|_{E_1} =- \mathscr{E}^{\textbf{n}}_{1,y}(u)|_{E_1} = u. \end{align*} The $$H^2$$-continuity estimate follows by applying Lemma 5.3 for the derivatives and the Cauchy–Schwarz inequality for the other terms; thus we have   \begin{align*} \| \mathscr{E}^{\textbf{n}}_1(u) \|_{2,\widehat{T}} \preccurlyeq \| u \|_{1/2,E_1} \preccurlyeq \| u \|_{1/2^*,E_1}. \end{align*} Step 2. In this step, we want to correct the values and the normal derivative on the second edge $$E_2$$ without changing the values and the normal derivative on the bottom edge $$E_1$$. For this purpose, we introduce the polynomials $$b(s) = 3s^2-2s^3$$ and $$c(s)=s^3-s^2$$, with the properties   \begin{align*} b(0) = b'(0) = b'(1) = 0,\quad b(1) = 1\quad \textrm{and} \quad c(0) = c'(0) = c(1) = 0,\quad c'(1) = 1, \end{align*} and use them as blending coefficients to define   \begin{align*} \mathscr{E}^{\textbf{n}}_2(u)(x,y):= \mathscr{E}^{\textbf{n}}_1(u)(x,y) - b\left(\frac{y}{1-x}\right)\mathscr{E}^{\textbf{n}}_1(u) (x,1-x) - c\left(\frac{y}{1-x}\right)(1-x) \mathscr{E}^{\textbf{n}}_{1,y}(u)(x,1-x). \end{align*} The idea is that the second term corrects the values and the last term corrects the normal derivative on the edge $$E_2$$. Indeed, we observe on $$E_1$$ as $$y=0$$ and on $$E_2$$ as $$y = 1-x$$ that   \begin{align} \mathscr{E}^{\textbf{n}}_2(u)|_{E_1}&= \mathscr{E}^{\textbf{n}}_1(u)|_{E_1}- b(0) \mathscr{E}^{\textbf{n}}_1(u)(x,1-x) - c(0) \mathscr{E}^{\textbf{n}}_{1,y}(u)(x,1-x) = 0,\nonumber \\ \mathscr{E}^{\textbf{n}}_2(u)|_{E_2}&= \mathscr{E}^{\textbf{n}}_1(u)|_{E_1}- b(1) \mathscr{E}^{\textbf{n}}_1(u)(x,1-x) - c(1) \mathscr{E}^{\textbf{n}}_{1,y}(u)(x,1-x) = 0.\label{ext::extnormstep2onedgetwoB} \end{align} (5.23) The derivative with respect to $$y$$ reads   \begin{align*} \mathscr{E}^{\textbf{n}}_{2,y}(u) &= \mathscr{E}^{\textbf{n}}_{1,y}(u) - \frac{1}{1-x}b'\left(\frac{y}{1-x}\right) \mathscr{E}^{\textbf{n}}_1(u) (x,1-x) - \frac{1}{1-x}c'\left(\frac{y}{1-x}\right)(1-x) \mathscr{E}^{\textbf{n}}_{1,y}(u) (x,1-x). \end{align*} On $$E_1$$ it follows, due to $$b'(0) = c'(0) = 0$$, that $$\mathscr{E}^{\textbf{n}}_{2,y}(u)|_{E_1}= \mathscr{E}^{\textbf{n}}_{1,y}(u)|_{E_1} = u$$, so the normal derivative has not changed in the second step. In a similar way, we also observe that $$\mathscr{E}^{\textbf{n}}_{2,y}(u)|_{E_2}=0$$. Now note that due to the constant zero value on the edge, see equation (5.23), we derive that the tangential derivative on the edge $$\nabla \mathscr{E}^{\textbf{n}}_2(u) \cdot \tau$$ on $$E_2$$ has to be zero. As $$\nabla \mathscr{E}^{\textbf{n}}_2(u) \cdot \tau|_{E_2} = 0 \Leftrightarrow - \mathscr{E}^{\textbf{n}}_{2,x}(u)|_{E_2} = \mathscr{E}^{\textbf{n}}_{2,y}(u)|_{E_2}$$ and $$\mathscr{E}^{\textbf{n}}_{2,y}(u)|_{E_2} = 0$$, it follows that $$\nabla \mathscr{E}^{\textbf{n}}_{2}(u) \cdot n|_{E_2} = 0$$, so the correction term induced a zero normal derivative on the second edge $$E_2$$. The $$H^2$$-estimate remains. As the first term of $$\mathscr{E}^{\textbf{n}}_2(u)$$ has already been analysed in the first step, we focus on just the correction terms. We start with first term $$A_4:= b(\frac{y}{1-x})\mathscr{E}^{\textbf{n}}_1(u) (x,1-x)$$ and the estimate for the $$y$$ derivative   \begin{align*} A_{4,y} = 6\frac{y}{1-x}\left(1-\frac{y}{1-x}\right) \int_0^1 u(x+s(1-x))a(s) ~\mathrm{d} s. \end{align*} Using the Cauchy–Schwarz inequality, $$\frac{y}{1-x} \le 1$$ on $$\widehat{T}$$ and the substitution $$t := x+s(1-x)$$, we get   \begin{align*} \| A_{4,y} \|_{0,\widehat{T}}^2 &\preccurlyeq \int_0^1\int_0^{1-x}\int_0^1 u(x+s(1-x))^2 ~\mathrm{d} s ~\mathrm{d} y ~\mathrm{d} x = \int_0^1\int_0^{1-x} \frac{1}{1-x}\int_x^1 u(t)^2 ~\mathrm{d} t ~\mathrm{d} y ~\mathrm{d} x \\ &= \int_0^1 \int_x^1 u(t)^2 ~\mathrm{d} t ~\mathrm{d} x \preccurlyeq \| u \|_{0,E_1}^2 \preccurlyeq \| u \|_{1/2^*,E_1}^2, \end{align*} and in a similar way we also bound $$\| A_{4,x} \|_{0,\widehat{T}}$$ and $$\| A_4 \|_{0,\widehat{T}}$$. The crucial point in this estimate was that we were able to use property $$\frac{y}{1-x} \le 1$$ twice, thus there were no bad coefficients any more. The estimates of the second-order derivatives are a little bit more tricky, as there remain some fractions with singularities. We start with the second-order derivative with respect to $$y$$ given by   \begin{align*} A_{4,yy} = 6\frac{1}{1-x}\left(1-\frac{2y}{1-x}\right) \int_0^1 u(x+s(1-x))a(s) ~\mathrm{d} s. \end{align*} The idea is to use Fubini’s theorem,   \begin{align*} \| A_{4,yy} \|_{0,\widehat{T}}^2 &\preccurlyeq \int_0^1\int_0^{1-x} \frac{1}{(1-x)^2}\int_0^1 u(x+s(1-x))^2 ~\mathrm{d} s ~\mathrm{d} y ~\mathrm{d} x \\ &= \int_0^1 \frac{1}{(1-x)^2}\int_x^1 u(x+s(1-x))^2 ~\mathrm{d} s ~\mathrm{d} y ~\mathrm{d} x = \iint\limits_{\substack{0 \le x \le 1\\ x \le t}} \frac{u(t)^2}{(1-x)^2}~\mathrm{d}(x,t) \\ &= \iint\limits_{\substack{0 \le t \le 1\\ t \le x}} \frac{u(t)^2}{(1-x)^2}~\mathrm{d}(x,t) = \int_0^1\int_0^{\mathrm{t}} \frac{1}{(1-x)^2} ~\mathrm{d} x~ u(t)^2 ~\mathrm{d} t \\ &= \int_0^1 \frac{1}{1-t}u(t)^2 ~\mathrm{d} t \preccurlyeq \| u \|_{0^*,E_1}^2 \preccurlyeq \| u \|_{1/2^*,E_1}^2. \end{align*} With the techniques just presented and similar estimates as in the proof in Theorem 5.4 and Lemma 5.3, all other derivatives of the second correction are bounded and we get   \begin{align*} \| \mathscr{E}^{\textbf{n}}_2(u) \|_{2,\widehat{T}} \preccurlyeq \| u \|_{1/2^*,E_1}. \end{align*} Step 3. Similarly to Step 2 we correct now the values on the last edge $$E_3$$ with two more corrections using the same blending coefficients; thus we define   \begin{align*} \mathscr{E}^{\textbf{n}}(u)(x,y)&:= \mathscr{E}^{\textbf{n}}_2(u)(x,y) - b\left(\frac{y}{x+y}\right)\mathscr{E}^{\textbf{n}}_2(u) (0,x+y) \\ &\quad\ - c\left(\frac{y}{x+y}\right)(x+y) \mathscr{E}^{\textbf{n}}_{2,x}(u)(0,x+y) +c\left(\frac{y}{x+y}\right)(x+y) \mathscr{E}^{\textbf{n}}_{2,y}(u)(0,x+y). \end{align*} As in Step 2 we use the properties of the blending coefficients and similar techniques to bound the integrals to derive (5.18), (5.19), (5.20) and (5.21):   \begin{gather*} \mathscr{E}^{\textbf{n}}(u)|_{E_1} = \mathscr{E}^{\textbf{n}}(u)|_{E_2} = \mathscr{E}^{\textbf{n}}(u)|_{E_3} = 0 \\ (\nabla \mathscr{E}^{\textbf{n}}(u)\cdot n)|_{E_2} = (\nabla \mathscr{E}^{\textbf{n}}(u)\cdot n)|_{E_3} = 0 \quad \textrm{and} \quad (\nabla \mathscr{E}^{\textbf{n}}(u)\cdot n)|_{E_1} = u, \\ \| \mathscr{E}^{\textbf{n}}(u) \|_{2,\widehat{T}} \preccurlyeq \| u \|_{1/2^*,E_1}. \end{gather*} It remains to show that $$\mathscr{E}^{\textbf{n}}(u)$$ belongs to $$\mathscr{P}^{k+1}(\widehat{T})$$. The idea is similar to the tangential extension. Looking at the definition of the first step (5.22), we increase the order by multiplying with $$y$$. The crucial parts are the correction terms as the blending polynomials $$b(y/(1-x))$$ and $$c(y/(1-x))$$ for the second step, and $$b(y/(x+y))$$ and $$c(y/(x+y))$$ for the third step produce singularities of order 3 in the points $$(0,0)$$ and $$(1,0)$$. To overcome this problem, note that the given polynomial has a zero of order 2 in the vertices; thus there exists a polynomial $$v \in \mathscr{P}^{k-2}(E_1)$$, such that $$u(x) = (1-x)^2v(x)$$. Using the definitions of the polynomials $$b$$ and $$c$$, the extension of the second step $$\mathscr{E}^{\textbf{n}}_2(u)$$ reads   \begin{align*} \mathscr{E}^{\textbf{n}}_2(u)(x,y) & = y \int_0^1 a(s)u(x+sy) ~\mathrm{d} s - \frac{3 y^2(1-x)-2y^3}{(1-x)^2}\int_0^1a(s) u(x+s(1-x)) ~\mathrm{d} s \\ &\quad{} - \frac{y^3-y^2(1-x)}{(1-x)^2}\int_0^1a'(s) u(x+s(1-x))s ~\mathrm{d} s. \end{align*} As $$u(x+s(1-x)) = (1-x)^2(1-s)^2v(x+s(1-x))$$ this leads to   \begin{align*} \mathscr{E}^{\textbf{n}}_2(u)(x,y) & = y \int_0^1 a(s)u(x+sy) ~\mathrm{d} s - (3 y^2(1-x)-2y^3)\int_0^1a(s) (1-s)^2v(x+s(1-x)) ~\mathrm{d} s \\ &\quad{} - (y^3-y^2(1-x))\int_0^1a'(s) (1-s)^2v(x+s(1-x))s ~\mathrm{d} s, \end{align*} thus $$\mathscr{E}^{\textbf{n}}_2(u) \in \mathscr{P}^{k+1}(\widehat{T})$$. For the third step, we do the same by writing $$u(x) = x^2w(x)$$ with $$w \in \mathscr{P}^{k-2}(E_1)$$, finally leading to $$\mathscr{E}^{\textbf{n}}(u) \in \mathscr{P}^{k+1}(\widehat{T})$$. □ Remark 5.6 In a similar way to the last step of the proof of Theorem 5.4, it is possible to define the normal extension $$\mathscr{E}^{\textbf{n}}$$ for the other two edges $$E_2$$ and $$E_3$$ using proper transformations. We then use the subscript $$\mathscr{E}^{\textbf{n}}_{E_i}(\cdot)$$ with $$i \in \{ 1,2,3 \}$$ to symbolize which extension is used. 5.4 Splitting into compatible and incompatible polynomials Using Theorem 5.5, it would now be possible to correct the normal derivative after a first extension using Theorem 5.4. The crucial point is that the polynomial would need a zero of order 2 in the vertices. The following theorem helps us later to provide a stable splitting of the correction into two parts. Theorem 5.7 Assume a given function $$u \in \mathscr{P}^k(E_1)$$, with $$u =0$$ on $$\partial E_1$$. There holds   \begin{align} \label{ext::divtheorem:A} |u'(1)| \preccurlyeq k^2 \| u \|_{1/2^*,E_1}. \end{align} (5.24) Furthermore, there exists a function $$e \in \mathscr{P}^k(E_1)$$ with $$e'(1) = 1$$ and $$e'(0) = e(0) = e(1) = 0$$ such that   \begin{align} \label{ext::divtheorem:B} \| e\|_{1/2^*,E_1} \preccurlyeq \frac{1}{k^2} \quad \textrm{and} \quad \| e\|_{0,E_1} \preccurlyeq \frac{1}{k^3}. \end{align} (5.25) Proof. We start with the second statement. For this purpose, we present the proof on the interval $$E$$ from which (5.25) follows with a linear transformation. We want to remind the reader of the definition of Jacobi polynomials with respect to the weight function $$(1-x)^\alpha$$ (see, for example, Abramowitz, 1974; Andrews et al., 1999),   \begin{align*} p^\alpha_n(x) := \frac{1}{2^nn!(1-x)^\alpha} \frac{{\rm{d}}}{{\rm{d}} x^n}\left((1-x)^\alpha(x^2-1)^n \right)\!, \quad n \in \mathbb{N}_0, \alpha > -1, \end{align*} where in the special case of $$\alpha = 0$$ the polynomials are called Legendre polynomials. For our proof, we use integrated Jacobi polynomials with $$\alpha=1$$,   \begin{align*} \widehat{p}_n(x) &:= -\int_x^1 p_{n-1}^1(s) ~\mathrm{d} s, \quad n\ge 1 \quad \textrm{and} \quad \widehat{p}_0 (x) := 1, \end{align*} and integrated Legendre polynomials,   \begin{align*} l_{n+1}(x) := -\int_x^1 p_n^0(s) ~\mathrm{d} s, \quad n \ge 0 \quad \textrm{and} \quad l_0 := -x + 1. \end{align*} These polynomials fulfil the following properties (see Andrews et al., 1999; Beuchler & Schöberl, 2006):   \begin{gather} \widehat{p}_n (1) = 0, \quad 1 \le n \le k \quad \textrm{and} \quad \widehat{p}'_n (1) = n, \quad 0 \le n \le k, \end{gather} (5.26)  \begin{gather} (2n+1) p_n^0 = (n+1) p^1_n - n p^1_{n-1}, \quad n\ge 0 \quad \textrm{and} \quad p_m^1 = \frac{1}{m+1} \sum\limits_{n=0}^m (2n+1) p_n^0, \quad m\ge 0, \end{gather} (5.27)  \begin{gather} (2n+1) l_{n+1} = p^0_{n+1} - p^0_{n-1}, \quad n>0, \end{gather} (5.28) where we have used $$p^0_{-1} := -1$$. Furthermore, we have a weighted $$L^2$$ orthogonality for the Jacobi polynomials, and due to the definition, a weighted orthogonality in the $$H^1$$-seminorm for the integrated Jacobi polynomials,   \begin{align} \label{ext::dividingtheorem::orthogonality} \int_{-1}^1 (1-x) \widehat{p}'_n(x) \widehat{p}'_m(x) ~\mathrm{d} x = \int_{-1}^1 (1-x) p^1_{n-1}(x) p^1_{m-1}(x) ~\mathrm{d} x = \delta_{n,m} \frac{2}{n+1}, \end{align} (5.29) where $$\delta$$ is the Kronecker delta. To find a proper candidate that fulfils the bounds (5.25), we first seek for the minimum of the weighted $$H^1$$-seminorm with proper restrictions; thus   \begin{align*} \tilde{e} := \arg \min\limits_{\substack{v \in \mathscr{P}^k \\ v(1) = 0 \\v'(1) = 1}} \int_{-1}^1(1-x)v'(x)^2 ~\mathrm{d} s = \arg \min\limits_{\substack{v \in \mathscr{P}^k \\ v(1) = 0 \\v'(1) = 1}} | v|^2_{1^*,E}. \end{align*} Using integrated Jacobi polynomials as the basis for $$\mathscr{P}^k(E),$$ we use the representation of $$\tilde{e}$$ with coefficients $$c_j$$ as $$\tilde{e}(x) = \sum\limits_{j=0}^k c_j \widehat{p}_j(x)$$. To determine the coefficients, so to explicitly solve the minimization problem, we first note that due to the boundary restrictions $$\tilde{e}(1) = 0$$, it is clear that $$c_0 =0$$ and with (5.26) we get $$\tilde{e}'(1) = \sum\limits_{j=1}^k c_j j = 1$$. Using (5.29) we furthermore have $$| \tilde{e}|^2_{1^*,E} = \sum\limits_{j=1}^k c_j^2\frac{2}{j+1}$$. Now we use the technique of Lagrangian multipliers; thus we define the function   \begin{align*} L(c_1,\dots,c_k,\lambda) = \sum\limits_{j=1}^k c_j^2\frac{2}{j+1} + \lambda \left(\sum\limits_{j=1}^k c_j j -1\right) \quad \textrm{with} \quad \frac{\partial L}{ \partial c_j} \stackrel{!}{=} 0 \quad \forall\, j \in \{1,\dots,k \} \quad \textrm{and} \quad \frac{\partial L}{\partial \lambda} \stackrel{!}{=} 0. \end{align*} Solving this leads to   \begin{align} \lambda = \frac{-48}{k(k+1)(k+2)(3k+1)} \quad \textrm{and} \quad c_j= \frac{-12 j(1+j)}{k(k+1)(k+2)(3k+1)}\label{ext::divi::coeffs} \end{align} (5.30) and   \begin{align} \label{ext:divtheorem:approx} | \tilde{e} |^2_{1*,E} = \sum\limits_{j=1}^k c_j^2 \frac{2}{1+j} = \frac{24}{3k^4+10k^3+9k^2+2k} \approx \frac{1}{k^4}. \end{align} (5.31) For the $$L^2$$-norm we observe using (5.27) and (5.28) that   \begin{align*} \tilde{e}(x) &=\sum\limits_{j=1}^k -c_j\int_x^1 p_{j-1}^1(s) ~\mathrm{d} s = \sum\limits_{j=1}^k -c_j\int_x^1 \frac{1}{j} \sum_{i=0}^{j-1}(2i+1)p_{i}^0(s) ~\mathrm{d} s = \sum\limits_{j=1}^k \frac{c_j}{j} \sum_{i=0}^{j-1}\underbrace{-(2i+1)\int_x^1p_{i}^0(s)}_{= (2i+1)l_{i+1}} ~\mathrm{d} s\\ & = \sum\limits_{j=1}^k \frac{c_j}{j} \sum_{i=0}^{j-1}\left(p_{i+1}^0(x) - p_{i-1}^0(x)\right) = \sum\limits_{j=1}^k \frac{c_j}{j} \left(p_{j}^0(x) + p_{j-1}^0(x)\right), \end{align*} and so with the definition of the coefficients (5.30) and using an inverse inequality (for example, Bernardi & Maday, 1997, page 253), also   \begin{align*} \| \tilde{e} \|^2_{0,E} = \sum_{j=1}^k \frac{c_j^2}{j^2} \| p_{j}^0 + p_{j-1}^0 \|_{0,E}^2 \preccurlyeq \sum_{j=1}^k \frac{j^2}{k^8} \underbrace{ \| p_{j}^0 \|_{0,E}^2}_{\preccurlyeq \frac{2}{2j+1}} \preccurlyeq \frac{j}{k^8} \sum_{j=1}^k1 \preccurlyeq \frac{1}{k^6} \quad \text{and} \quad \| \tilde{e}\|^2_{1,E} \preccurlyeq \frac{1}{k^2}. \end{align*} Using a linear transformation $$F$$ from $$E$$ to $$E_1$$, we set $$e(x) := \tilde{e}(F^{-1}(x))\frac{x^2}{2}$$ to see $$e(0) = e(1) = e'(0) = 0$$ and $$e'(1)=1$$, and   \begin{align} \label{ext:divtheorem:twoestimates} \| e \|_{0,E_1} \preccurlyeq \frac{1}{k^3} \quad \text{and} \quad \| e \|_{1,E_1} \preccurlyeq \frac{1}{k}. \end{align} (5.32) Similarly to the proof of Lemma 5.3, we now use the real method of interpolation of spaces. As $$u=0$$ on $$\partial E_1$$ we have $$u \in H^1_0(E_1)$$; thus together with $$H^{1/2}_{00}(E_1) = [L^2(E_1),H_0^1(E_1)]$$ and the definition of the norm on an interpolated space we have with (5.32),   \begin{align*} \| e\|_{1/2^*, E_1} \preccurlyeq \sqrt{ \| e\|_{0, E_1}\| e\|_{1, E_1} } \preccurlyeq \frac{1}{k^2}, \end{align*} so (5.25) is proved. For the proof of (5.24), we first define an extension from the edge to the triangle by   \begin{align*} \psi(u)(x,y):= \int_0^1 a(s) u(x+sy) ~\mathrm{d} s \quad \textrm{with} \quad a(s) = 4-6s, \end{align*} and so $$u'(1) = \frac{\partial \psi}{\partial x}(1,0)$$. Using Lemma 5.3 we get $$\| \psi \|_{1,\widehat{T}} \preccurlyeq \| u \|_{1/2,E_1}$$. Next we define the mean value along the line $$l_x := \{(x,y): 0 \le x \le 1, y \in [0,1-x] \}$$,   \begin{align*} \overline{u}(x,y) := \frac{1}{1-x} \int_x^1 \psi(x,s) ~\mathrm{d} s = \int_0^1 \psi(x,(1-x)s) ~\mathrm{d} s. \end{align*} Due to $$\frac{\partial \overline{u}}{\partial y} = 0$$, it follows with $$\frac{\partial \overline{u}}{\partial x} := \overline{u}'$$ that   \begin{align*} | \overline{u} |_{1,\widehat{T}}^2 = \int_0^1 \int_0^{1-x} \overline{u}'(x)^2 ~\mathrm{d} y ~\mathrm{d} x = \int_0^1 (1-x) \overline{u}'(x)^2 ~\mathrm{d} x = | \overline{u} |_{1*,E_1}^2 \end{align*} and so   \begin{align*} | \overline{u} |_{1*,E_1} = | \overline{u} |_{1,\widehat{T}} \preccurlyeq \| \psi \|_{1,\widehat{T}} \preccurlyeq \| u \|_{1/2, E_1}. \end{align*} Using (5.31) and $$\tilde{e}'(1)=1$$, we furthermore show that $$| \overline{u}'(1) | \preccurlyeq k^2 | \overline{u} |_{1/2^*,E_1} \preccurlyeq k^2 \| u \|_{1/2, E_1}$$, and as   \begin{align*} \overline{u}'(1) = u'(1) \overbrace{ \int_0^1a(s) ~\mathrm{d} s}^{=1} - \frac{1}{2} u'(1) \overbrace{\int_0^1a(s)s ~\mathrm{d} s}^{=0} = u'(1), \end{align*} we finally get   \begin{align*} | u'(1) | \preccurlyeq k^2 \| u \|_{1/2, E_1} \preccurlyeq k^2 \| u \|_{1/2^*, E_1}.\quad\square \end{align*} 5.5 Proof of Theorem 5.1 Proof. In the first step, we use Theorem 5.4 to find a function $$\mathscr{E}^{\tau}(\textbf{u})$$ with a proper tangential derivative. For the difference $$\textbf{u}_c:= \textbf{u} - \nabla \mathscr{E}^{\tau}(\textbf{u})$$ it follows that on the boundary $$\partial \widehat{T}$$ we have $$\textbf{u}_c \cdot \tau = 0$$. Now let $$\textbf{n}_i$$ be the normal vector on the edge $$E_i$$ and $$u_{\textbf{n}_i}:= \textbf{u}_c \cdot \textbf{n}_i$$, so the remaining difference of the normal component. The idea is now to split this error in two parts to use Theorems 5.5 and 5.7. We start with the lower edge $$E_1$$ and define $$u_{1}:= \textbf{u}_{c} \cdot ((x,y) - V_2) \in \mathscr{P}^{k+1}(\widehat{T})$$, where $$V_2 = (0,1)$$ is the vertex opposite to $$E_1$$. As $$((x,y) - V_2) \approx \tau$$ on the edges $$E_2$$ and $$E_3$$ we have $$u_{1}|_{E_2} = u_{1}|_{E_3} = 0$$. On the lower edge we have $$((x,y) - V_2) = (x,-1)$$ and as $$\textbf{u}_c \cdot \tau = 0$$, thus the first component of $$\textbf{u}_c$$ is equal to zero, we get $$u_{1}|_{E_1} = u_{\textbf{n}_1} \in \mathscr{P}^k(E_1)$$. Using Theorem 5.7, we find two functions $$e_0,e_1 \in \mathscr{P}^k(E_1)$$ with   \begin{align*} e'_1(1) =1,\quad e'_1(0) = e_1(0) = e_1(1) = 0 \quad \textrm{and} \quad e'_0(0) =1,\quad e'_0(1) = e_0(0) = e_0(1) = 0, \end{align*} where we mirrored the edge $$E$$ in Theorem 5.7 to find $$e_0$$. We are now able to split the error to define a good and a bad part on the edge $$E_1$$ by   \begin{align*} u_{\textbf{n}_1}^b := (u_{1}|_{E_1})'(1) e_1 + (u_{1}|_{E_1})'(0) e_0 \quad \textrm{and} \quad u_{\textbf{n}_1}^g := u_{\textbf{n}_1} - u_{\textbf{n}_1}^b. \end{align*} The second function $$u_{\textbf{n}_1}^g$$ is good in the sense of having a zero of order 2 in the vertices, so $$u_{\textbf{n}_1}^g \in \mathscr{P}_{00}^k(E_1)$$; thus we apply Theorem 5.5. For the other two edges, we proceed similarly (see Remark 5.6) to finally define   \begin{align*} \mathscr{E}(\textbf{u}):= \mathscr{E}^{\tau}(\textbf{u}) + \mathscr{E}^{\textbf{n}}_{E_1}(u_{\textbf{n}_1}^g) + \mathscr{E}^{\textbf{n}}_{E_2}(u_{\textbf{n}_2}^g) + \mathscr{E}^{\textbf{n}}_{E_3}(u_{\textbf{n}_3}^g). \end{align*} Note that due to (5.19) and (5.20) the normal derivatives of the different corrections do not interfere. As $$\mathscr{E}^{\textbf{n}}(u_{\textbf{n}_i}^g) = 0$$ (see (5.18)) on the boundary $$\partial \widehat{T}$$ for $$i=1,2,3$$ the corresponding tangential derivative is also zero; thus we have   \begin{align*} \nabla \mathscr{E}(\textbf{u}) \cdot \tau = \nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \tau + \nabla \mathscr{E}^{\textbf{n}}_{E_1}(u_{\textbf{n}_1}^g)\cdot \tau + \nabla \mathscr{E}^{\textbf{n}}_{E_2}(u_{\textbf{n}_2}^g)\cdot \tau + \nabla \mathscr{E}^{\textbf{n}}_{E_3}(u_{\textbf{n}_3}^g)\cdot \tau = \nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \tau = \textbf{u} \cdot \tau, \end{align*} so property (5.1) is proved. For $$\mathscr{E}^{\textbf{n}}_{E_1}(u_{\textbf{n}_1}^g),$$ we get using (5.21), (5.25) and (5.24) as $$u_{1}|_{E_1}=0$$ on $$\partial E_1$$,   \begin{align*} \| \mathscr{E}^{\textbf{n}}_{E_1}(u_{\textbf{n}_1}^g) \|_{2,\widehat{T}} &\preccurlyeq \| u_{\textbf{n}_1}^g \|_{1/2^*,E_1} = \| u_{\textbf{n}_1} - u_{\textbf{n}_1}^b \|_{1/2^*,E_1} \\ &\preccurlyeq \| u_{\textbf{n}_1} \|_{1/2^*,E_1} + |(u_{1}|_{E_1})'(1)| \| e_1\|_{1/2^*,E_1} + |(u_{1}|_{E_1})'(0)| \| e_0\|_{1/2^*,E_1} \\ &\preccurlyeq \| u_{\textbf{n}_1} \|_{1/2^*,E_1} + \|u_{1}\|_{1/2^*,E_1} \frac{k^2}{k^2} + \|u_{1}\|_{1/2^*,E_1} \frac{k^2}{k^2} \preccurlyeq \|u_{1}\|_{1/2^*,E_1}. \end{align*} As $$u_{1}|_{E_2} = u_{1}|_{E_3} = 0$$ we bound $$\|u_{1}\|_{1/2^*,E_1}$$ by the $$H^1$$-norm on the triangle; thus we get the estimate $$\| \mathscr{E}^{\textbf{n}}_{E_1}(u_{\textbf{n}_1}^g) \|_{2,\widehat{T}} \preccurlyeq \| u_{1} \|_{1,\widehat{T}} \preccurlyeq \| \textbf{u}_{c} \|_{1,\widehat{T}}$$. Using the same arguments for the other two normal extensions and inequality (5.5) implies   \begin{align*} \| \mathscr{E}(\textbf{u}) \|_{2,\widehat{T}} \le \| \mathscr{E}^{\tau}(\textbf{u}) \|_{2,\widehat{T}} + 3 \| \textbf{u}_{c} \|_{1,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}} + \| \textbf{u} - \nabla \mathscr{E}^{\tau}(\textbf{u}) \|_{1,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}; \end{align*} thus property (5.3) is proved. To show (5.2) first note that on the boundary $$\partial \widehat{T}$$ we have   \begin{align*} \nabla \mathscr{E}(\textbf{u}) \cdot \textbf{n} &= \nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \textbf{n} + \sum\limits_{i=1}^3 \nabla \mathscr{E}^{\textbf{n}}_{E_i}(u_{\textbf{n}_i}^g) \cdot n =\nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \textbf{n} + \overbrace{\sum\limits_{i=1}^3 u_{\textbf{n}_i}}^{=\textbf{u}_c \cdot n} - \sum\limits_{i=1}^3 u_{\textbf{n}_i}^b\\ & = \nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \textbf{n} + \textbf{u} \cdot \textbf{n} - \nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \textbf{n} - \sum\limits_{i=1}^3 u_{\textbf{n}_i}^b= \textbf{u} \cdot \textbf{n} - \sum\limits_{i=1}^3 u_{\textbf{n}_i}^b, \end{align*} and as $$u_{\textbf{n}_i}^b|_{E_j} = 0$$ for $$j \neq i$$ it follows that $$\| \left(\textbf{u} - \nabla \mathscr{E}(\textbf{u}) \right) \cdot \textbf{n} \|_{0,\partial \widehat{T}} \le \sum\limits_{i=1}^3 \| u_{\textbf{n}_i}^b \|_{0,E_i}$$. As before, we use (5.25) and (5.24) to get   \begin{align*} \| u_{\textbf{n}_1}^b \|_{0,E_1} &\preccurlyeq |(u_{1}|_{E_1})'(1)| \| e_1 \|_{0,E_1} + |(u_{1}|_{E_1})'(0)| \| e_0 \|_{0,E_1} \preccurlyeq \frac{1}{k} \|u_1\|_{1/2^*,E_1}\\ &\preccurlyeq \frac{1}{k} \|u_1\|_{1,\widehat{T}} \preccurlyeq \frac{1}{k} \|\textbf{u}_c\|_{1,\widehat{T}} = \frac{1}{k} \|\textbf{u} - \nabla \mathscr{E}^{\tau}(\textbf{u}) \|_{1,\widehat{T}} \preccurlyeq \frac{1}{k}\|\textbf{u}\|_{1,\widehat{T}}, \end{align*} and with a similar estimate for $$u_{\textbf{n}_2}^b$$ and $$u_{\textbf{n}_3}^b$$ we finally get (5.2):   \begin{align*} \| \left(\textbf{u} - \nabla \mathscr{E}(\textbf{u}) \right) \cdot \textbf{n} \|_{0,\partial \widehat{T}} \preccurlyeq \frac{1}{k}\|\textbf{u}\|_{1,\widehat{T}}.\quad\square \end{align*} References Abramowitz, M. & Stegun, I. A. ( 1974) Handbook of Mathematical Functions, With Formulas, Graphs, and Mathematical Tables . Dover Publications. Ainsworth, M. & Coggins, P. ( 2002) A uniformly stable family of mixed $$hp$$-finite elements with continuous pressures for incompressible flow. IMA J. Numer. Anal. , 22, 307– 327. Google Scholar CrossRef Search ADS   Ainsworth, M. & Demkowicz, L. ( 2009) Explicit polynomial preserving trace liftings on a triangle. Math. Nachr. , 282, 640– 658. Google Scholar CrossRef Search ADS   Andrews, G., Askey, R. & Roy, R. ( 1999) Special Functions . Encyclopedia of Mathematics and its Applications. Cambridge: Cambridge University Press. Google Scholar CrossRef Search ADS   Arnold, D. N., Brezzi, F., Cockburn, B. & Marini, L. D. ( 2001/02) Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J. Numer. Anal. , 39, 1749– 1779. Google Scholar CrossRef Search ADS   Babuška, I. & Suri, M. ( 1987) The $$h-p$$ version of the finite element method with quasiuniform meshes. ESAIM Math. Model. Numer. Anal. , 21, 199– 238. Google Scholar CrossRef Search ADS   Baker, G. A., Jureidini, W. N. & Karakashian, O. A. ( 1990) Piecewise solenoidal vector fields and the Stokes problem. SIAM J. Numer. Anal. , 27, 1466– 1485. Google Scholar CrossRef Search ADS   Belgacem, F. B. ( 1994) Polynomial extensions of compatible polynomial traces in three dimensions. Comput. Methods Appl. Mech. Engrg. , 116, 235– 241. Google Scholar CrossRef Search ADS   Bergh, J. & Löfström, J. ( 1976) Interpolation Spaces: An Introduction . Grundlehren der mathematischen Wissenschaften. Berlin-New York: Springer. Google Scholar CrossRef Search ADS   Bernardi, C. & Maday, Y. ( 1990) Relèvement polynomial de traces et applications. RAIRO Modél. Math. Anal. Numér. , 24, 557– 611. Google Scholar CrossRef Search ADS   Bernardi, C. & Maday, Y. ( 1992) Approximations spectrales de problèmes aux limites elliptiques.  Mathématiques et Applications. Berlin: Springer. Bernardi, C. & Maday, Y. ( 1997) Spectral methods. Handbook of Numerical Analysis  ( Ciarlet P. G. & Lions J. L. eds), vol. 5. Amsterdam: North-Holland, pp. 209– 485. Google Scholar CrossRef Search ADS   Bernardi, C. & Maday, Y. ( 1999) Uniform inf-sup conditions for the spectral discretization of the Stokes problem. Math. Models Methods Appl. Sci. , 9, 395– 414. Google Scholar CrossRef Search ADS   Beuchler, S. & Schöberl, J. ( 2006) New shape functions for triangular $$p$$-FEM using integrated Jacobi polynomials. Numer. Math. , 103, 339– 366. Google Scholar CrossRef Search ADS   Boffi, D., Fortin, M. & Brezzi, F. ( 2013) Mixed Finite Element Methods and Applications . Computational Mathematics. Berlin: Springer. Google Scholar CrossRef Search ADS   Brennecke, C., Linke, A., Merdon, C. & Schöberl, J. ( 2015) Optimal and pressure-independent $$L^2$$ velocity error estimates for a modified Crouzeix–Raviart Stokes element with BDM reconstructions. J. Comput. Math. , 33, 191– 208. Google Scholar CrossRef Search ADS   Brenner, S. C., Gu, S., Gudi, T. & Sung, L.-y. ( 2012) A quadratic $$C^{\circ}$$ interior penalty method for linear fourth order boundary value problems with boundary conditions of the Cahn–Hilliard type. SIAM J. Numer. Anal. , 50, 2088– 2110. Google Scholar CrossRef Search ADS   Brenner, S. C., Gudi, T. & Sung, L.-y. ( 2010) An a posteriori error estimator for a quadratic $$C^0$$-interior penalty method for the biharmonic problem. IMA J. Numer. Anal. , 30, 777– 798. Google Scholar CrossRef Search ADS   Brezzi, F. & Falk, R. S. ( 1991) Stability of higher-order Hood–Taylor method. SIAM J. Numer. Anal. , 28, 581– 590. Google Scholar CrossRef Search ADS   Chernov, A. ( 2012) Optimal convergence estimates for the trace of the polynomial $$L^2$$-projection operator on a simplex. Math. Comp. , 81, 765– 787. Google Scholar CrossRef Search ADS   Cockburn, B., Gopalakrishnan, J., Nguyen, N. C., Peraire, J. & Sayas, F.-J. ( 2011) Analysis of HDG methods for Stokes flow. Math. Comp. , 80, 723– 760. Google Scholar CrossRef Search ADS   Cockburn, B., Kanschat, G. & Schötzau, D. ( 2004) The local discontinuous Galerkin method for the Oseen equations. Math. Comp. , 73, 569– 593. Google Scholar CrossRef Search ADS   Cockburn, B., Kanschat, G. & Schotzau, D. ( 2005) A locally conservative LDG method for the incompressible Navier–Stokes equations. Math. Comp. , 74, 1067– 1095. Google Scholar CrossRef Search ADS   Cockburn, B., Kanschat, G. & Schötzau, D. ( 2007) A note on discontinuous Galerkin divergence-free solutions of the Navier–Stokes equations. J. Sci. Comput. , 31, 61– 73. Google Scholar CrossRef Search ADS   Cockburn, B., Kanschat, G., Schötzau, D. & Schwab, C. ( 2002) Local discontinuous Galerkin methods for the Stokes system. SIAM J. Numer. Anal. , 40, 319– 343. Google Scholar CrossRef Search ADS   Cockburn, B., Nguyen, N. C. & Peraire, J. ( 2010) A comparison of HDG methods for Stokes flow. J. Sci. Comput. , 45, 215– 237. Google Scholar CrossRef Search ADS   Costabel, M. & McIntosh, A. ( 2010) On Bogovskiĭ and regularized Poincaré integral operators for de Rham complexes on Lipschitz domains. Math. Z. , 265, 297– 320. Google Scholar CrossRef Search ADS   Demkowicz, L. & Babuška, I. ( 2003) $$p$$ interpolation error estimates for edge finite elements of variable order in two dimensions. SIAM J. Numer. Anal. , 41, 1195– 1208. Google Scholar CrossRef Search ADS   Demkowicz, L., Gopalakrishnan, J. & Schöberl, J. ( 2008) Polynomial extension operators. I. SIAM J. Numer. Anal. , 46, 3006– 3031. Google Scholar CrossRef Search ADS   Demkowicz, L., Gopalakrishnan, J. & Schöberl, J. ( 2009) Polynomial extension operators. II. SIAM J. Numer. Anal. , 47, 3293– 3324. Google Scholar CrossRef Search ADS   Demkowicz, L., Gopalakrishnan, J. & Schöberl, J. ( 2012) Polynomial extension operators. Part III. Math. Comp. , 81, 1289– 1326. Google Scholar CrossRef Search ADS   Donea, J. & Huerta, D. ( 2003) Finite Element Methods for Flow Problems . Hoboken, NJ: John Wiley. Google Scholar CrossRef Search ADS   Egger, H. & Schöberl, J. ( 2010) A hybrid mixed discontinuous Galerkin finite-element method for convection-diffusion problems. IMA J. Numer. Anal. , 30, 1206– 1234. Google Scholar CrossRef Search ADS   Egger, H. & Waluga, C. ( 2013) $$hp$$ analysis of a hybrid DG method for Stokes flow. IMA J. Numer. Anal. , 33, 687– 721. Google Scholar CrossRef Search ADS   Elman, H. C., Silvester, D. J. & Wathen, A. J. ( 2005) Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics . Numerical Mathematics and Scientific Computation. Oxford, NY: Oxford University Press. Fu, G., Jin, Y. & Qiu, W. ( 2016) Parameter-free superconvergent $$H(\mathrm{div})$$-conforming HDG methods for the Brinkman equations. ArXiv e-prints . Georgoulis, E. H., Hall, E. & Melenk, J. M. ( 2010) On the suboptimality of the $$p$$-version interior penalty discontinuous Galerkin method. J. Sci. Comput. , 42, 54– 67. Google Scholar CrossRef Search ADS   Girault, V. & Raviart, P.-A. ( 1986) Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms . Computational Mathematics. Berlin: Springer. Extended version. Google Scholar CrossRef Search ADS   Girault, V., Rivière, B. & Wheeler, M. F. ( 2005) A discontinuous Galerkin method with nonoverlapping domain decomposition for the Stokes and Navier–Stokes problems. Math. Comp. , 74, 53– 84. Google Scholar CrossRef Search ADS   Glowinski, R. ( 2003) Finite Element Methods for Incompressible Viscous Flow . Handbook of Numerical Analysis. Amsterdam: Elsevier, pp. 3– 1176. Grisvard, P. ( 1985) Elliptic Problems in Nonsmooth Domains . Classics in Applied Mathematics. Boston, MA: Society for Industrial and Applied Mathematics. Houston, P., Schwab, C. & Süli, E. ( 2002) Discontinuous $$hp$$-finite element methods for advection-diffusion-reaction problems. SIAM J. Numer. Anal. , 39, 2133– 2163. Google Scholar CrossRef Search ADS   Karakashian, O. & Katsaounis, T. ( 2000) A discontinuous Galerkin method for the incompressible Navier–Stokes equations. Discontinuous Galerkin Methods (Newport, RI, 1999) . Lecture Notes in Computer Science and Engineering, vol. 11. Berlin: Springer, pp. 157– 166. Google Scholar CrossRef Search ADS   Karniadakis, G. E. & Sherwin, S. J. ( 2005) Spectral/hp Element Methods for Computational Fluid Dynamics . Numerical Mathematics and Scientific Computation. Oxford: Oxford University Press. Google Scholar CrossRef Search ADS   Könnö, J. & Stenberg, R. ( 2011) H (div)-conforming finite elements for the Brinkman problem. Math. Models Methods Appl. Sci. , 21, 2227– 2248. Google Scholar CrossRef Search ADS   Lederer, P. ( 2016) Pressure-robust discretizations for Navier–Stokes equations: divergence-free reconstruction for Taylor–Hood elements and high order hybrid discontinuous Galerkin methods. Master’s Thesis , Vienna Technical University. Lehrenfeld, C. ( 2010) Hybrid discontinuous Galerkin methods for solving incompressible flow problems. Master’s Thesis, Rheinisch Westfalischen Technischen Hochschule Aachen . Lehrenfeld, C. & Schöberl, J. ( 2016) High order exactly divergence-free hybrid discontinuous Galerkin methods for unsteady incompressible flows. Comput. Methods Appl. Mech. Engrg. , 307, 339– 361. Google Scholar CrossRef Search ADS   Linke, A. ( 2014) On the role of the Helmholtz decomposition in mixed methods for incompressible flows and a new variational crime. Comput. Methods Appl. Mech. Engrg. , 268, 782– 800. Google Scholar CrossRef Search ADS   Linke, A., Matthies, G. & Tobiska, L. ( 2016) Robust arbitrary order mixed finite element methods for the incompressible Stokes equations with pressure independent velocity errors. ESAIM Math. Model. Numer. Anal. , 50, 289– 309. Google Scholar CrossRef Search ADS   Linke, A. & Merdon, C. ( 2016) On velocity errors due to irrotational forces in the Navier–Stokes momentum balance. J. Comput. Phys. , 313, 654– 661. Google Scholar CrossRef Search ADS   Maday, Y. ( 1989) Relèvements de traces polynomiales et interpolations hilbertiennes entre espaces de polynômes. C. R. Acad. Sci. Paris Sér. I Math. , 309, 463– 468. Muñoz-Sola, R. ( 1997) Polynomial liftings on a tetrahedron and applications to the $$h$$-$$p$$ version of the finite element method in three dimensions. SIAM J. Numer. Anal. , 34, 282– 314. Google Scholar CrossRef Search ADS   Nguyen, N. C., Peraire, J. & Cockburn, B. ( 2010) A hybridizable discontinuous Galerkin method for Stokes flow. Comput. Methods Appl. Mech. Engrg. , 199, 582– 597. Google Scholar CrossRef Search ADS   Nguyen, N. C., Peraire, J. & Cockburn, B. ( 2011) An implicit high-order hybridizable discontinuous Galerkin method for the incompressible Navier–Stokes equations. J. Comput. Phys. , 230, 1147– 1170. Google Scholar CrossRef Search ADS   Peetre, J. ( 1963) Nouvelles propriétés d’espaces d’interpolation. C. R. Acad. Sci. Paris , 256, 1424– 1426. Rivière, B. ( 2008) Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation.  Frontiers in Applied Mathematics, vol. 35. Philadelphia: Society for Industrial and Applied Mathematics (SIAM), pp. xxii+190. Google Scholar CrossRef Search ADS   Schöberl, J. ( 1997) NETGEN An advancing front 2D/3D-mesh generator based on abstract rules. Comput. Vis. Sci. , 1, 41– 52. Google Scholar CrossRef Search ADS   Schöberl, J. ( 2014) C++11 implementation of finite elements in NGSolve. Technical Report.  Institute for Analysis and Scientific Computing, Vienna University of Technology. Schöberl, J., Melenk, J. M., Pechstein, C. & Zaglmayr, S. ( 2008) Additive Schwarz preconditioning for $$p$$-version triangular and tetrahedral finite elements. IMA J. Numer. Anal. , 28, 1– 24. Google Scholar CrossRef Search ADS   Schötzau, D., Schwab, C. & Toselli, A. ( 2002) Mixed hp-DGFEM for incompressible flows. SIAM J. Numer. Anal. , 40, 2171– 2194. Google Scholar CrossRef Search ADS   Schwab, C. C. ( 1998) p- and hp- Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics . Numerical Mathematics and Scientific Computation. Oxford, NY: Clarendon Press. Stamm, B. & Wihler, T. P. ( 2010) $$hp$$-optimal discontinuous Galerkin methods for linear elliptic problems. Math. Comp. , 79, 2117– 2133. Google Scholar CrossRef Search ADS   Stenberg, R. & Suri, M. ( 1996) Mixed hp finite element methods for problems in elasticity and Stokes flow. Num. Math. , 72, 367– 389. Google Scholar CrossRef Search ADS   Su, Y., Chen, L., Li, X. & Xu, C. ( 2016) On the inf-sup constant of a triangular spectral method for the Stokes equations. Comput. Methods Appl. Math. , 16, 507– 522. Google Scholar CrossRef Search ADS   Toselli, A. ( 2002) $$hp$$ discontinuous Galerkin approximations for the Stokes problem. Math. Models Methods Appl. Sci. , 12, 1565– 1597. Google Scholar CrossRef Search ADS   Zhang, S. ( 2009) A family of $$Q_{k+1,k}\times Q_{k,k+1}$$ divergence-free finite elements on rectangular grids. SIAM J. Numer. Anal. , 47, 2090– 2107. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# Polynomial robust stability analysis for $$H$$ (div)-conforming finite elements for the Stokes equations

, Volume Advance Article – Aug 29, 2017
29 pages

/lp/ou_press/polynomial-robust-stability-analysis-for-h-div-conforming-finite-uEfdk9xuLX
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx051
Publisher site
See Article on Publisher Site

### Abstract

Abstract In this article, we consider a discontinuous Galerkin method for the discretization of the Stokes problem. We use $$H(\textrm{div})$$-conforming finite elements, as they provide major benefits such as exact mass conservation and pressure-independent error estimates. The main aspect of this work lies in the analysis of high-order approximations. We show that the considered method is uniformly stable with respect to the polynomial order $$k$$ and provides optimal error estimates $$\| {\textbf{u} - \textbf{u}_h} \|_{1_h} + \| {{\it {\Pi}}^{Q_h} p-p_h} \|_{0} \le c \left(h/k \right)^s \| \textbf{u} \|_{s+1}$$. To derive these estimates, we prove a $$k$$-robust Ladyženskaja-Babuška-Brezzi (LBB) condition. This proof is based on a polynomial $$H^2$$-stable extension operator. This extension operator itself is of interest for the numerical analysis of $$C^0$$-continuous discontinuous Galerkin methods for fourth-order problems. 1. Introduction In this article, we consider the numerical solution of the Stokes equations on a bounded domain $${\it {\Omega}} \subset \mathbb{R}^2$$,   \begin{align} \label{int::stokesequations} \begin{array}{rcll} -\nu {\it {\Delta}} \textbf{u} + \nabla p &=& \textbf{f} &\textrm{in } {\it {\Omega}}, \\ \text{div}~ \textbf{u} & = & 0 & \textrm{in } {\it {\Omega}}, \end{array} \end{align} (1.1) with boundary conditions $$\textbf{u} = 0$$ on $$\partial {\it {\Omega}}$$, where $$\nu =\mathrm{const}$$ is the kinematic viscosity, $$\textbf{u}$$ is the velocity field, $$p$$ is the pressure and $$\textbf{f}$$ are external forces. The approximation of the Navier–Stokes problem is well analysed, and many different finite element methods have been introduced (see, for example, Girault & Raviart, 1986; Donea & Huerta, 2003; Glowinski, 2003; Elman et al., 2005). Furthermore, discontinuous Galerkin (DG) finite element methods for elliptic problems became popular (see, for example, Arnold et al., 2002; Houston et al., 2002; Rivière, 2008) and also for flow problems (Baker et al., 1990; Karakashian & Katsaounis, 2000; Cockburn et al., 2002; Schötzau et al., 2002; Toselli, 2002; Cockburn et al., 2004, 2005, 2007; Girault et al., 2005). In the last two mentioned works, an elementwise divergence-free discretization is constructed. In contrast, we consider a normal-continuous $$H(\operatorname{div})$$-conforming method introduced by Cockburn et al. (2007). This has different advantages such as local conservation, the possibility to use an upwinding scheme for convection-dominated flows and pressure robust (independent) error estimates due to exactly divergence-free velocity test functions (see Linke, 2014; Brennecke et al., 2015; Linke et al., 2016; Linke & Merdon, 2016). Note that this discretization was also used for the Brinkman problem in Könnö & Stenberg (2011). To reduce the computational costs of DG methods, we also want to mention hybrid DG (HDG) methods, where new variables are introduced on the skeleton and a static condensation technique is used for the element unknowns (see Cockburn et al., 2010, 2011; Egger & Schöberl, 2010; Nguyen et al., 2010, 2011; Egger & Waluga, 2013) and for $$H(\operatorname{div})$$-conforming methods (Lehrenfeld, 2010; Fu et al., 2016; Lehrenfeld & Schöberl, 2016). The method we use is well analysed with respect to the mesh size $$h$$ and provides optimal error estimates. The main contribution of this article is to show that the method is also uniformly stable with respect to the polynomial order $$k$$. For this purpose, we prove that the constant $$\beta$$ of the LBB condition   \begin{align*} \sup_{\textbf{0} \neq \textbf{v}_h \in \textbf{V}_h} \frac{b(\textbf{v}_h,q_h)}{\| \textbf{v}_h \|_{\textbf{V}_h}} \geq \beta \| q_h \|_{Q_h} \quad \forall\, q_h \in Q_h, \end{align*} is independent of the order $$k$$. Together with standard continuity and ellipticity estimates this leads to a stable high-order method. Note that small adaptations of our results also lead to polynomial robustness for HDG methods. High-order methods for incompressible flow problems are of theoretical and practical importance. In Bernardi & Maday (1997) and Karniadakis & Sherwin (2005), the authors consider a spectral method on the unit cube using polynomials of order $$k$$ and $$k-2$$ for the velocity and the pressure, respectively. The resulting method leads to $$\beta(k) = \mathscr{O}(k^{-\frac{d-1}{2}})$$, where $$d$$ is the space dimension. The same method on triangles is discussed in Schwab (1998) with similar results. Furthermore, the bad influence of a dependency on $$k$$ of the LBB constant for an iterative method for solving the Navier–Stokes equation was analysed. Understanding the problem, an improvement was achieved in Bernardi & Maday (1999). They used polynomials of partial order $$k$$ for the velocity and polynomials of total order at most $$k-1$$ for the pressure, resulting in a uniformly stable method. Similar achievements for $$hp$$ mixed finite element methods are accomplished in Stenberg & Suri (1996). Therein, different combinations of elements on quadrilaterals like continuous polynomials of order $$k$$ for the velocity and discontinuous polynomials of order $$k-2$$ for the pressure are discussed. An exact analysis is presented but again revealed a dependency on $$k$$. They also considered different tensor product spaces for each component of the velocity. A similar approach leading to an optimal exactly divergence-free method was presented in Zhang (2009) using polynomials of order $$k+1$$ in the $$x$$-direction and polynomials of order $$k$$ in the $$y$$-direction for the velocity in $$x$$-direction, and vice versa for the velocity in $$y$$-direction. Using proper degrees of freedom, this leads to a similar method on quadrilaterals as we use on triangles. Another approach, combining the tensor product structure on quadrilaterals and the advantage of approximating more complex geometries using triangles, is analysed in Su et al. (2016). The key of this method is to use the Duffy transformation and a proper pair of approximation spaces, which leads to $$\beta(k) = \mathscr{O}(k^{-\frac{1}{2}})$$ with the drawback of using rational functions for the approximation. We also want to mention the method considered in Ainsworth & Coggins (2002), where a uniformly stable approximation using a continuous ansatz for the velocity and pressure is presented. They adopted the ideas of Bernardi & Maday (1999) but enriched the pressure space to overcome the lack of convergence order that would appear using just a continuous version of this method. Considering continuous approximations, also the famous Taylor–Hood elements on triangles and quadrilaterals (see Brezzi & Falk, 1991; Boffi et al., 2013) have to be discussed. Although these methods were shown to be stable with respect to the mesh size $$h$$, numerical evidence predicts that they are not uniformly stable with respect to $$k$$. Of course, high-order methods were also used for discontinuous finite element methods. We want to mention the work of Toselli (2002) and Schötzau et al. (2002), where an analysis for $$hp$$-DG methods on quadrilaterals is presented but revealed a dependency on the order $$k$$, and also the work of Egger & Waluga (2013), where an HDG method on triangles and quadrilaterals with similar results is introduced. The rest of the article is structured in the following way. In Section 2 we present the Stokes equation and the considered discretization method. Furthermore, we present a short proof of the continuous divergence stability to motivate the existence of an $$H^2$$-stable polynomial extension operator, which is used to prove the main theorem in Section 3. In Section 4 we take a look at some numerical examples and finally present the construction of an $$H^2$$-stable polynomial extension operator in Section 5. 1.1 Preliminaries We assume an open bounded domain $${\it {\Omega}} \subset \mathbb{R}^2$$ with a Lipschitz boundary $${\it {\Gamma}}$$; thus for every point on the boundary there exists a Lipschitz-continuous mapping $${\it {\Phi}}^x$$. If this mapping is furthermore differentiable up to order $$m$$, we say $${\it {\Gamma}} \in \mathscr{C}^{m,1}$$. On $${\it {\Omega}}$$, we define a shape-regular triangulation $$\mathscr{T}$$. Furthermore, we assume $$\mathscr{T}$$ to be quasi-uniform; thus, there exists one global mesh size $$h$$ such that $$h \approx \textrm{diam}(T)\ \textrm{for all}\ T \in \mathscr{T}$$. The set of edges, with respect to the triangulation $$\mathscr{T}$$, will be defined as $$\mathscr{F}$$. We call $$\widehat{T} := \{(x,y): 0\le x \le 1, 0 \le y \le 1, x+y \le 1 \}$$ the reference element with the edges $$E_1:= \{(x,0): 0\le x \le 1 \}$$, $$E_2:= \{(x,1-x): 0\le x \le 1 \}$$ and $$E_3:= \{(0,y): 0\le y \le 1 \}$$, and define the interval $$E:=\{(x,0): -1\le x \le 1 \}$$. On all triangles we use $$\textbf{n}$$ and $$\tau$$ as symbols for the normal and tangential vectors. For all subsets $$\omega \subseteq{\it {\Omega}}$$ with $$\gamma := \partial \omega$$ we have the space $$L^2(\omega)$$ with the norm $$\| \cdot \|_{0,\omega}$$ and the Sobolev spaces $$H^1(\omega)$$, $$H^2(\omega), H^{1/2}(\gamma)$$ with the corresponding Sobolev norms $$\| \cdot \|_{s,\omega}$$. For better readability, we skip the index $$\omega$$ if it is clear which domain is considered. On the edge $$E_1$$ we define the weighted $$L^2$$- and $$H^{1/2}$$-norm of a function $$u$$ as   \begin{align*} \| u \|_{0^{*},E_1}^2 := \int_0^1 \left(\frac{1}{x} + \frac{1}{1-x}\right) u(x)^2 ~\mathrm{d} s \quad \textrm{and} \quad \| u \|_{1/2^*,E _1}^2 := | u |^2_{1/2,E _1} + \| u \|_{0^{*},E_1}^2. \end{align*} Furthermore, we use the closed subspaces   \begin{align*} L^2_0({\it {\Omega}}) &:= \{ q \in L^2({\it {\Omega}}): \textstyle\int_{\it {\Omega}} q ~\mathrm{d} x = 0\} \quad \textrm{and} \quad H^1_0({\it {\Omega}}) := \{ u \in H^1({\it {\Omega}}): \mathrm{tr} ~ u = 0 ~\textrm{on}~ \partial {\it {\Omega}}\}, \end{align*} and the polynomial spaces   \begin{align*} \mathscr{P}^m(\mathscr{T}) &:= \{ v : v|_T \in \mathscr{P}^m(T)~ \forall\, T \in \mathscr{T} \} = \prod_{ T \in \mathscr{T}} \mathscr{P}^m(T) \quad \textrm{and} \quad \pmb{\mathscr{P}}^m(\mathscr{T}) := [\mathscr{P}^m(\mathscr{T})]^2, \\ \mathscr{P}_{00}^m(E_1) &:= \{v \in \mathscr{P}^m(E_1): v(0) = v'(0) = v(1) = v'(1) = 0\}. \end{align*} Also we define the following subspaces of vectorial polynomials on the reference triangle:   \begin{align*} \pmb{\mathscr{P}}^m_{\tau}(\widehat{T}):= \{\textbf{v} \in \pmb{\mathscr{P}}^m(\widehat{T}): \textstyle\int_{\partial \widehat{T}} \textbf{v} \cdot \tau = 0\} \quad \textrm{and} \quad \pmb{\mathscr{P}}^m_{\textbf{n}}(\widehat{T}):= \{\textbf{v} \in \pmb{\mathscr{P}}^m(\widehat{T}): \textstyle\int_{\partial \widehat{T}} \textbf{v} \cdot \textbf{n} = 0\}. \end{align*} In this article, we use an index notation for partial derivatives; thus for an arbitrary function $$u$$ we write $$u_{,x} := \frac{\partial u}{\partial x}$$ and $$u_{,y} := \frac{\partial u}{\partial y}$$ and in a similar manner for second-order derivatives. We write $$(x,y)^{\mathrm{t}}$$ as the transposed vector of $$(x,y)$$ and use $$\perp$$ as symbol for a counterclockwise rotation by $$\pi/2$$, thus $$(x,y)^\perp := (-y,x)$$. By this we also define the two-dimensional operator $$\text{curl}~ = \nabla^\perp$$. Finally, we introduce the expression $$a \preccurlyeq b$$ when there exists a constant $$c$$ independent of $$a, b$$, the polynomial order and the mesh size such that $$a \le c b$$. 2. Discretization of the Stokes problem In this section, we present the discretization of the stationary incompressible Stokes equations (1.1) from Cockburn et al. (2007). They use mixed-order finite element spaces with the polynomial orders $$k$$ and $$k-1$$ for the velocity and the pressure, respectively. To ensure a local conservative and energy-stable method, we provide exactly divergence-free velocity fields using an $$H(\operatorname{div})$$-conforming method; thus every discrete velocity field $$\textbf{u}_h$$ is in   \begin{align*} H(\operatorname{div})({\it {\Omega}}) := \{ \textbf{u} \in [L^2({\it {\Omega}})]^2: \text{div}~ \textbf{u} \in L^2({\it {\Omega}}) \}. \end{align*} To ensure $$\textbf{u}_h \in H(\operatorname{div})({\it {\Omega}})$$, we demand normal continuity across each edge resulting in the approximation space for the velocity   \begin{align*} \textbf{V}_h := \{ \textbf{u}_h \in \pmb{\mathscr{P}}^k(T): [\![{\textbf{u}_h \cdot n}]\!] = 0 ~ \forall\, E \in \mathscr{F} \} \subset H(\operatorname{div},{\it {\Omega}}), \end{align*} where $$[\![{\cdot}]\!]$$ is the usual jump operator on interior edges and just the identity on boundary edges. For the pressure space, we assume no continuity across edges:   \begin{align*} Q_h :=\prod\limits_{T \in \mathscr{T}}\mathscr{P}^{k-1}(T) \cap L^2_0({\it {\Omega}}). \end{align*} Note, that this pair of finite element spaces fulfils $$\text{div}~ \textbf{V}_h = Q_h$$ and thus, a weakly incompressible velocity field $$\textbf{u}_h \in \textbf{V}_h$$ is also exactly divergence-free:   \begin{align*} \int_{\it {\Omega}} \text{div}~{\textbf{u}_h} ~ q ~\mathrm{d} x= 0 \quad \forall\, q \in Q_h \quad \Rightarrow \quad \text{div}~{\textbf{u}_h} = 0 \quad \textrm{in } {\it {\Omega}}. \end{align*} Furthermore, using $$\pmb{\{} \cdot \pmb{\}}$$ as the symbol for the mean value on an edge $$E \in \mathscr{F}$$ we define the bilinear form are   \begin{align}\label{bilinearform::a} a(\textbf{u}_h,\textbf{v}_h) &:= \sum\limits_{T \in \mathscr{T}} \int_T \nu \nabla \textbf{u}_h : \nabla \textbf{v}_h ~\mathrm{d} x - \sum\limits_{E \in \mathscr{F}} \int_{E} \nu \pmb{\{} {\nabla \textbf{u}_h \cdot n} \pmb{\}} [\![{\textbf{v}_h \cdot \tau}]\!] ~\mathrm{d} s\nonumber\\ &\quad{} -\sum\limits_{E \in \mathscr{F}} \int_{E} \nu \pmb{\{} {\nabla \textbf{v}_h \cdot n} \pmb{\}} [\![{\textbf{u}_h \cdot \tau}]\!] ~\mathrm{d} s + \sum\limits_{E \in \mathscr{F}} \int_{E} \nu \frac{\alpha k^2}{h} [\![{\textbf{u}_h \cdot \tau}]\!] [\![{\textbf{v}_h \cdot \tau}]\!] ~\mathrm{d} s, \end{align} (2.1) where $$\alpha >0$$ and $$k^2/h$$ is a standard stability coefficient (see Houston et al., 2002) and the bilinear form and linear form   \begin{align} \label{bilinearform::b} b(\textbf{u}_h, q_h) := - \sum\limits_{T \in \mathscr{T}} \int_T \text{div}~{\textbf{u}_h} ~ q_h ~\mathrm{d} x \quad \textrm{and} \quad l(\textbf{v}_h):= \sum\limits_{T \in \mathscr{T}} \int_T \textbf{f} \cdot \textbf{v}_h ~\mathrm{d} x. \end{align} (2.2) The discrete Stokes problem now reads, find $$(\textbf{u}_h, p_h)$$ in $$\textbf{V}_h\times Q_h$$ such that   \begin{align} \label{dis::stokesproblem} \begin{array}{rcll} a(\textbf{u}_h,\textbf{v}_h) + b(\textbf{v}_h, p_h) &=& l(\textbf{v}_h)\quad &\forall\, \textbf{v}_h \in \textbf{V}_h, \\ b(\textbf{u}_h, q_h) & = & 0 \quad & \forall\, q_h \in Q_h. \end{array} \end{align} (2.3) On the spaces $$Q_h$$ and $$\textbf{V}_h$$, we use the $$L^2$$-norm $$\| {\cdot} \|_{0}$$ and $$\| {\textbf{v}_h} \|_{1_h}^2 := \sum\limits_{T \in \mathscr{T}} \| {\nabla \textbf{v}_h} \|_{0,T}^2 + \sum\limits_{E \in \mathscr{F}} \frac{k^2}{h} \| {[\![{\textbf{v}_h \cdot \tau}]\!]} \|_{0,E}^2$$, respectively. We have chosen homogeneous Dirichlet boundary conditions, but other choices are possible as well. In the $$H(\operatorname{div})$$-conforming DG method, we distinguish between essential boundary conditions for the normal component and for the tangential components. The essential normal boundary condition is incorporated into the finite element space, whereas an essential tangential boundary condition is treated in the DG sense (see Georgoulis et al., 2010). Natural boundary conditions are treated in the common way as well. This distinction is particularly useful for the so-called slip boundary conditions, where only the normal velocity is constrained to $$\textbf{u} \cdot \textbf{n} = 0$$. Lemma 2.1 For a proper choice of the stabilization parameter $$\alpha >0$$ in (2.1), there exist constants $$\alpha_1>0, \alpha_2>0$$ and $$\alpha_3>0$$ independent of the mesh size $$h$$ and the polynomial order $$k$$ such that $$a(\cdot,\cdot)$$ is coercive,   \begin{align*} a(\textbf{v}_h,\textbf{v}_h) \ge \nu \alpha_3 \| {\textbf{v}_h} \|_{1_h}^2 \quad \forall\, \textbf{v}_h \in \textbf{V}_h, \end{align*} and $$a(\cdot,\cdot)$$ and $$b(\cdot,\cdot)$$ are continuous,   \begin{align*} |a(\textbf{u}_h,\textbf{v}_h)| \le \nu \alpha_1 \| {\textbf{u}_h} \|_{1_h}\| {\textbf{v}_h} \|_{1_h} \quad \forall\, \textbf{v}_h, \textbf{u}_h \in \textbf{V}_h, \quad |b(\textbf{u}_h,q_h)| \le \alpha_2 \| {\textbf{u}_h} \|_{1_h} \| {q_h} \|_{0} \quad \forall\, \textbf{u}_h \in \textbf{V}_h, \forall\, q_h \in Q_h. \end{align*} Proof. The continuity properties of $$a(\cdot,\cdot)$$ and $$b(\cdot,\cdot)$$ follow by the definition of the norm $$\| {\cdot} \|_{1_h}$$. The coercivity follows with similar arguments to Stamm & Wihler (2010). □ Theorem 2.2 There exists a constant $$\beta > 0$$ independent of the polynomial order $$k$$ and the mesh size $$h$$ such that   \begin{align*} \sup\limits_{\textbf{0} \neq \textbf{v}_h \in \textbf{V}_h} \frac{b(\textbf{v}_h, q_h)}{\| {\textbf{v}_h} \|_{1_h}} \ge \beta \| {q_h} \|_{0} \quad \forall\, q_h \in Q_h. \end{align*} Proof. The proof is presented in Section 3. □ Theorem 2.3 Assume $$(\textbf{u}_h, p_h)$$ in $$\textbf{V}_h\times Q_h$$ is the solution of the discrete problem (2.3) and $$(u,p)$$ is the exact solution of the Stokes problem (1.1). Furthermore, let us assume regularity $$\textbf{u} \in [H^{s+1}({\it {\Omega}})]^2$$ and $$p \in H^s({\it {\Omega}})$$. Then, there exists a constant $$c_{\rm err}$$ independent of the mesh size $$h$$ and the polynomial order $$k$$ such that for $$s \ge 1$$ and $$s \le k$$ it holds that   \begin{align}\label{optimalerrorest} \| {\textbf{u} - \textbf{u}_h} \|_{1_h} + \| {{\it {\Pi}}^{Q_h} p-p_h} \|_{0} \le c_{\rm err} \left(\frac{h}{k} \right)^{s} \| \textbf{u} \|_{s+1}, \end{align} (2.4) where $${\it {\Pi}}^{Q_h}$$ is the $$L^2$$ projector onto $$Q_h$$. Proof. In a first step, we discretize the Poisson problem $$-\nu {\it {\Delta}} \textbf{u} = f - \nabla p$$ using a DG approximation, i.e., let $$\textbf{w}_h$$ be the solution of   \begin{align*} a(\textbf{w}_h, \textbf{v}_h) = l(v_h) - b(\textbf{v}_h, p) \quad \forall\, \textbf{v}_h \in \textbf{V}_h, \end{align*} where we used integration by parts for $$\int_{\it {\Omega}} \nabla p \cdot \textbf{v}_h ~\mathrm{d} x$$ and $$\textbf{v}_h \in H(\operatorname{div})({\it {\Omega}})$$. Using the estimate from Stamm & Wihler (2010, Section 3.2), including the properties of the $$L^2$$ projector on triangles in Chernov (2012, Theorem 1.1), we get   \begin{align} \label{errorestlaplace} \| {\textbf{u} - \textbf{w}_h} \|_{1_h} \le \tilde{c}_{\mathrm{err}} \left(\frac{h}{k} \right)^{s} \| \textbf{u} \|_{s+1}. \end{align} (2.5) Since $$(\textbf{u}_h, p_h)$$ is the solution of the discrete problem (2.3) we have   \begin{align*} \begin{array}{rcll} a(\textbf{u}_h,\textbf{v}_h) + b(\textbf{v}_h, p_h) &=& b(\textbf{v}_h, p) + a(\textbf{w}_h, \textbf{v}_h) \quad &\forall\, \textbf{v}_h \in \textbf{V}_h, \\ b(\textbf{u}_h, q_h) & = & 0 \quad & \forall\, q_h \in Q_h, \end{array} \end{align*} and thus   \begin{align*} \begin{array}{rcll} a(\textbf{u}_h - \textbf{w}_h,\textbf{v}_h) + b(\textbf{v}_h, p_h-p) &=& 0 \quad &\forall\, \textbf{v}_h \in \textbf{V}_h, \\ b(\textbf{u}_h- \textbf{w}_h, q_h) & = & -b(\textbf{w}_h, q_h) \quad & \forall\, q_h \in Q_h. \end{array} \end{align*} Due to the property $$\text{div}~ \textbf{V}_h = Q_h$$, we replace the term $$b(\textbf{v}_h, p_h-p)$$ in the first row by $$b(\textbf{v}_h, p_h - {\it {\Pi}}^{Q_h} p)$$. Using $$p_h - {\it {\Pi}}^{Q_h} p \in Q_h$$ and the standard stability estimate of saddle point problems (see, for example, Boffi et al., 2013, Theorem 4.2.3), we get   \begin{align*} \| {\textbf{u}_h - \textbf{w}_h} \|_{1_h} + \frac{\beta}{\alpha_1} \| p_h - {\it {\Pi}}^{Q_h} p \|_0 &\le \left(1+ \frac{2\alpha_1}{\alpha_3\beta}\right) \| \text{div}~{\textbf{w}_h} \|_0 \\ &= \left(1+ \frac{2\alpha_1}{\alpha_3\beta}\right) \| \text{div}~{\textbf{w}_h} - \underbrace{\text{div}~{\textbf{u}}}_{=0} \|_0 \le \left(1+ \frac{2\alpha_1}{\alpha_3\beta}\right) \| {\textbf{u} - \textbf{w}_h} \|_{1_h}. \end{align*} The estimation (2.4) follows with (2.5), the triangle inequality and the robustness of the constants $$\alpha_1, \alpha_3$$ and $$\beta$$ with respect to the mesh size $$h$$ and $$k$$. □ Remark 2.4 The introduced method can also be used for a non-quasi-uniform triangulation, thus using triangles with different sizes $$h_T$$ and furthermore individual polynomial degrees $$k_T$$. By that we get a similar local error estimation to Stamm & Wihler (2010),   \begin{align*} \| {\textbf{u} - \textbf{u}_h} \|_{1_h} + \| {{\it {\Pi}}^{Q_h} p-p_h} \|_{0} \le c_{\rm err} \sqrt{ \sum\limits_{T \in \mathscr{T}}\left(\frac{h_T}{k_T} \right)^{2s} \| \textbf{u} \|_{s+1,T}^2}. \end{align*} 2.1 Continuous LBB condition as motivation for an $$H^2$$-extension In this section, we present a proof for the infinite-dimensional version of the LBB condition of the Stokes problem as can be found in Bernardi & Maday (1997). For this, we define the velocity space $$\textbf{V} := [H^1_0({\it {\Omega}})]^2$$, the pressure space $$Q := L^2_0({\it {\Omega}})$$ and show that   \begin{align*} \sup\limits_{\textbf{0} \neq \textbf{v} \in \textbf{V}} \frac{b(\textbf{v}, q)}{\|\textbf{v}\|_1} \ge \beta_\infty \| {q} \|_{0} \quad \forall\, q \in Q. \end{align*} Note that the LBB condition is equivalent to the existence of an $$H^1$$-stable right inverse of the divergence operator. Theorem 2.5 Let $${\it {\Omega}} \subset \mathbb{R}^2$$ be a bounded domain with a smooth Lipschitz boundary $$\partial {\it {\Omega}} \in \mathscr{C}^{1,1}$$. The divergence operator from $$\textbf{V}$$ to $$Q$$ is onto, so for every $$q \in Q$$ there exists a $$\textbf{v} \in \textbf{V}$$ such that   \begin{align*} \text{div}~{\textbf{v}} = q \quad \textrm{and} \quad \|\textbf{v}\|_1 \preccurlyeq \| {q} \|_{0}. \end{align*} Proof. Let $$q$$ be an arbitrary function in $$Q$$. In the first step, we consider the Poisson problem $${\it {\Delta}} \varphi = q$$ in $${\it {\Omega}}$$ with Neumann boundary conditions $$\nabla \varphi \cdot \textbf{n} = 0$$ on $$\partial {\it {\Omega}}$$. Due to the zero mean value of $$q$$ this problem has a unique solution in $$H^1({\it {\Omega}}) / \mathbb{R}$$. Now set $$\textbf{v} := \nabla \varphi$$. It follows that $$\text{div}~{\textbf{v}} = {\it {\Delta}} \varphi = q$$ and, using a regularity result for the Poisson problem, also $$\| v\|_1 = \| \varphi \|_2 \preccurlyeq \| {q} \|_{0}$$. Furthermore, note that we already have $$\textbf{v} \cdot n = \nabla \varphi \cdot n = 0$$ on the boundary $$\partial {\it {\Omega}}$$, so the idea is to construct a correction for the tangential component to satisfy the zero boundary values of $$\textbf{V}$$. For this, we seek for a function $$\psi \in H^2({\it {\Omega}})$$, such that   \begin{align*} \psi = 0 \quad \textrm{on } \partial {\it {\Omega}} \quad \textrm{and} \quad \frac{\partial \psi}{\partial n} = - \textbf{v} \cdot \tau \quad \textrm{on } \partial {\it {\Omega}} \quad \textrm{and} \quad \| \psi \|_2 \preccurlyeq \| \textbf{v} \|_1. \end{align*} The existence of $$\psi$$ holds true as we assumed a smooth boundary (see Bernardi & Maday, 1997, Theorem 1.12). Now set $$\tilde{\textbf{v}}:= \textbf{v} + \text{curl}~{\psi}$$ to get $$\text{div}~{\tilde{\textbf{v}}} = \text{div}~{\textbf{v}} + \text{div}~{\text{curl}~{\psi}} = q$$ in $${\it {\Omega}}$$ and on the boundary $$\partial {\it {\Omega}}$$:   \begin{align*} \tilde{\textbf{v}} \cdot \textbf{n} = \textbf{v} \cdot \textbf{n} + \text{curl}~{\psi} \cdot \textbf{n} = \nabla \psi \cdot \tau = 0 \quad \textrm{and} \quad \tilde{\textbf{v}} \cdot \tau = \textbf{v} \cdot \tau + \text{curl}~{\psi} \cdot \tau = \textbf{v} \cdot \tau + \nabla \psi \cdot \textbf{n} = 0. \end{align*} Finally, due to the $$H^2$$-continuity of $$\psi$$, we get $$\|\tilde{\textbf{v}}\|_1 = \|\textbf{v}\|_1 + \|\text{curl}~{\psi}\|_1 \preccurlyeq \|\textbf{v}\|_1 \preccurlyeq \| {q} \|_{0}$$. □ It follows immediately (with a scaling of $$\tilde{\textbf{v}}$$ by $$-1$$) that   \begin{align*} \sup\limits_{\textbf{0} \neq \textbf{v} \in \textbf{V}} \frac{b(\textbf{v}, q)}{\|\textbf{v}\|_1} \succcurlyeq \frac{-\int_{\it {\Omega}} \text{div}~{\tilde{\textbf{v}}}{q} ~\mathrm{d} x }{\|\tilde{\textbf{v}}\|_1} \succcurlyeq \frac{\| {q} \|_{0}^2}{\| {q} \|_{0}} \succcurlyeq \| {q} \|_{0}. \end{align*} The crucial part of this proof was the existence of the correction $$\psi$$ that is stable in the $$H^2$$-norm. This holds true in the case of a smooth boundary of $${\it {\Omega}}$$, but as the boundary of an element $$T \in \mathscr{T}$$ is just in $$\mathscr{C}^{0,1}$$, we get a problem when we want to adapt this proof to show the main Theorem 2.2. Such $$H^2$$-stable extensions of boundary values for nonregular boundaries $$\partial {\it {\Omega}}$$, which are defined as a union of a finite number of regular parts, are considered in Grisvard (1985) and Bernardi & Maday (1992). For this, they assume that compatibility conditions are satisfied at the points where two parts of the boundary gather. Those conditions are quite restrictive, thus do not hold true for all traces of polynomials. For this reason, such a proof cannot be used, for example, in the case of continuous velocity elements. Considering only normal continuous approximations with a tangential continuity only in a DG sense creates enough freedom to construct such an $$H^2$$-extension. 3. Robust high-order LBB estimation In this section, we present the proof of Theorem 2.2 in similar steps as in the proof of Theorem 2.5. For this, we assume the existence of a stable $$H^2$$-extension which is then presented in Section 5. Theorem 3.1 ($$H^2$$-extension) For every $$k$$ there exists an operator $$\mathscr{E} : \pmb{\mathscr{P}}_{\textbf{n}}^k(\widehat{T}) \rightarrow \mathscr{P}^{k+1}(\widehat{T}),$$ such that for $$\textbf{u}_h \in \pmb{\mathscr{P}}_{\textbf{n}}^k(\widehat{T})$$ it holds that   \begin{align} &\text{curl}~{\mathscr{E}(\textbf{u}_h)} \cdot \textbf{n} = \textbf{u}_h \cdot \textbf{n} \quad \textrm{on } \partial \widehat{T}, \label{lbb:h2ext:propA} \\ \end{align} (3.1)  \begin{align} & \|(\textbf{u}_h - \text{curl}~{\mathscr{E}(\textbf{u}_h)}) \cdot \tau \|_{0,\partial \widehat{T}} \preccurlyeq \frac{1}{k} \| \textbf{u}_h \|_{1,\widehat{T}}, \label{lbb:h2ext:propB} \\ \end{align} (3.2)  \begin{align} & \| \mathscr{E}(\textbf{u}_h) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u}_h \|_{1,\widehat{T}}.\label{lbb:h2ext:propC} \end{align} (3.3) Proof. See Section 5. □ We first show the LBB condition on the reference triangle $$\widehat{T}$$. For this, we define the spaces $$\widehat{\textbf{V}}_h := \pmb{\mathscr{P}}^k(\widehat{T})$$ with $$\widehat{\textbf{V}}_{h,0} := \{ \hat{\textbf{v}}_h \in \widehat{\textbf{V}}_h: \hat{\textbf{v}}_h \cdot n = 0 ~ \textrm{on } \partial \widehat{T}\}$$ and $$\widehat{Q}_h := \mathscr{P}^{k-1}(\widehat{T}) \cap L^2_0(\widehat{T})$$. The norm $$\| {\cdot} \|_{1_h}$$ on $$\widehat{\textbf{V}}_h$$ now reads   \begin{align*} \|\hat{\textbf{v}}_h\|^2_{1_h,\widehat{T}} = \|{\nabla \hat{\textbf{v}}_h}\|^2_{0,\widehat{T}} + \sum\limits_{E \in \partial \widehat{T}} k^2 \| \hat{\textbf{v}}_h \cdot \tau \|^2_{0,E} \quad \forall\, \hat{\textbf{v}}_h \in \widehat{\textbf{V}}_h. \end{align*} Theorem 3.2 (Local LBB condition) There exists a constant $$\beta > 0$$ independent of the polynomial order $$k$$ such that   \begin{align*} \sup\limits_{\textbf{0} \neq \hat{\textbf{v}}_h \in \widehat{\textbf{V}}_{h,0}} \frac{b(\hat{\textbf{v}}_h, \hat{q}_h)}{\|\hat{\textbf{v}}_h\|_{1_h,\widehat{T}}} \ge \beta \|\hat{q}_h\|_{0,\widehat{T}} \quad \forall\, \hat{q}_h \in \widehat{Q}_h. \end{align*} Proof. Let $$\hat{q}_h \in \widehat{Q}_h$$ be an arbitrary function. For a fixed point $$\hat{x} \in \widehat{T}$$, we define a local Poincaré operator $$\mathscr{Z}_{\hat{x}}:\mathscr{P}^{k-1}(\widehat{T}) \rightarrow [\mathscr{P}^{k}(\widehat{T})]^2$$ as introduced in Costabel & McIntosh (2010) by   \begin{align*} \hat{q}_h(x) \mapsto \mathscr{Z}_{\hat{x}}(\hat{q}_h)(x):= -(x-\hat{x}) \int_0^1 t \hat{q}_h(\gamma(t)) ~\mathrm{d} t, \end{align*} with $$\gamma(t) := \hat{x} + t(x-\hat{x})$$, and by integrating over all points in $$\widehat{T}$$ we furthermore define   \begin{align*} \hat{\textbf{u}}_{h}^1(x) := -\int_{\widehat{T}} \theta(\hat{x}) \mathscr{Z}_{\hat{x}}(\hat{q}_h)(x) ~\mathrm{d} \hat{x}, \end{align*} where $$\theta \in C^\infty_0(\widehat{T})$$ is a smooth function. We observe that $$-\text{div}~{\hat{\textbf{u}}_{h}^1} = \hat{q}_h$$ and $$\| \hat{\textbf{u}}_{h}^1\|_{1,\widehat{T}} \preccurlyeq \|\hat{q}_h\|_{0,\widehat{T}}$$ (see Costabel & McIntosh, 2010, Corollary 3.4). As $$\hat{q}_h \in L^2_0(\widehat{T})$$, we see that $$\hat{\textbf{u}}_{h}^1 \in \pmb{\mathscr{P}}_{\textbf{n}}^k(\widehat{T}),$$ so we apply the extension operator of Theorem 3.1 to define $$\hat{\textbf{u}}_{h}^2:= \text{curl}~{\mathscr{E}(\hat{\textbf{u}}_{h}^1)}$$. By that we get for $$\hat{\textbf{u}}_h := \hat{\textbf{u}}_{h}^1 - \hat{\textbf{u}}_{h}^2$$ and property (3.1),   \begin{align*} -\text{div}~{\hat{\textbf{u}}_h} = -\text{div}~{\hat{\textbf{u}}_{h}^1} + \text{div}~{\hat{\textbf{u}}_{h}^2} = \hat{q}_h \textrm{ in } \widehat{T} \quad \textrm{and} \quad \hat{\textbf{u}}_h \cdot n = \hat{\textbf{u}}_{h}^1 \cdot n - \hat{\textbf{u}}_{h}^2 \cdot n = 0 \textrm{ on } \partial \widehat{T}. \end{align*} Together with (3.3) and (3.2) we get   \begin{align*} \|\hat{\textbf{u}}_h\|^2_{1_h,\widehat{T}} &=\|{\nabla \hat{\textbf{u}}_{h}^1 - \nabla \hat{\textbf{u}}_{h}^2}\|^2_{0,\widehat{T}} + \sum\limits_{E \in \partial \widehat{T}} k^2 \| (\hat{\textbf{u}}_{h}^1 - \hat{\textbf{u}}_{h}^2) \cdot \tau \|^2_{0,E} \preccurlyeq \|\nabla \hat{\textbf{u}}_{h}^1\|^2_{0,\widehat{T}} + \|{\nabla \hat{\textbf{u}}_{h}^2}\|^2_{0,\widehat{T}} +\| \hat{\textbf{u}}_h \|_{1,\widehat{T}}^2 \\ &\preccurlyeq \| \hat{\textbf{u}}_h \|_{1,\widehat{T}}^2 + \|{\nabla \hat{\textbf{u}}_{h}^2}\|^2_{0,\widehat{T}} = \| \hat{\textbf{u}}_h \|_{1,\widehat{T}}^2 + | \mathscr{E}(\hat{\textbf{u}}_{h}^1) |^2_{2,\widehat{T}} \preccurlyeq \| \hat{\textbf{u}}_h \|_{1,\widehat{T}}^2 \preccurlyeq \|\hat{q}_h\|_{0,\widehat{T}}^2. \end{align*} As $$\hat{\textbf{u}}_h \in \widehat{\textbf{V}}_{h,0}$$ we can bound the supremum from below   \begin{align*} \sup\limits_{\textbf{0} \neq \hat{\textbf{v}}_h \in \widehat{\textbf{V}}_{h,0}} \frac{b(\hat{\textbf{v}}_h, \hat{q}_h)}{\|\hat{\textbf{v}}_h\|_{1_h,\widehat{T}}} \succcurlyeq \frac{-\int_{\widehat{T}} \text{div}~{\hat{\textbf{u}}_h} ~{\hat{q}_h} ~\mathrm{d} x }{\|\hat{\textbf{u}}_h\|_{1_h,\widehat{T}}} \succcurlyeq \frac{ \|\hat{q}_h\|_{0,\widehat{T}}^2}{ \|\hat{q}_h\|_{0,\widehat{T}}} \succcurlyeq \|\hat{q}_h\|_{0,\widehat{T}}.\quad\square \end{align*} In the following, we derive a similar LBB estimate for each physical triangle. For this purpose, we define for all $$T \in \mathscr{T}$$ the spaces $$\textbf{V}_h(T) := \pmb{\mathscr{P}}^k(T)$$ with $$\textbf{V}_{h,0}(T) := \{ \textbf{v}_h \in \textbf{V}_h(T): \textbf{v}_h \cdot n = 0 ~ \textrm{on } \partial T\}$$ and $$Q_h(T) := \mathscr{P}^{k-1}(T) \cap L^2_0(T)$$. Corollary 3.3 For each triangle $$T \in \mathscr{T}$$ there exists a constant $$\beta > 0$$ independent of the polynomial order $$k$$ and the mesh size $$h$$ such that   \begin{align*} \sup\limits_{\textbf{0} \neq \textbf{v}_h \in \textbf{V}_{h,0}(T)} \frac{b(\textbf{v}_h, q_h)}{\|\textbf{v}_h\|_{1_h}} \ge \beta \|q_h\|_{0,T} \quad \forall\, q_h \in Q_h(T). \end{align*} Proof. Let $$F$$ be the mapping from the reference triangle $$\widehat{T}$$ to $$T$$. For an arbitrary function $$q_h \in Q_h(T)$$ we define $$\hat{q}_h := \mathscr{J} q_h \in \widehat{Q}_h$$ where $$\mathscr{J} := \det F'$$. On the reference triangle we use Theorem 3.2 to find a function $$\hat{\textbf{u}}_h \in \widehat{\textbf{V}}_{h,0}$$ with $$-\text{div}~{\hat{\textbf{u}}_h} = \hat{q}_h$$ and $$\|\hat{\textbf{u}}_h\|_{1_h,\widehat{T}} \preccurlyeq \| \hat{q}_h\|_{0,\widehat{T}}$$ (compare Theorem 2.5). Using the Piola transformation $$\mathbb{P}$$ we define $$\textbf{u}_h := \mathbb{P}(\hat{\textbf{u}}_h) = \frac{1}{\mathscr{J}} F' \hat{\textbf{u}}_h$$. As $$\mathbb{P}$$ preserves the normal component we have $$\textbf{u}_h \in \textbf{V}_{h,0}(T)$$ and due to the proper mapping of the divergence also $$-\text{div}~{\textbf{u}_h} = q_h$$. Next note that due to the shape regularity of $$\mathscr{T}$$ we have $$|| F'|| \approx h$$ and $$| \mathscr{J}| \approx h^2$$. By standard scaling arguments it follows that $$\frac{1}{h} \|\hat{\textbf{u}}_h\|_{1_h,\widehat{T}} \approx \|\textbf{u}_h\|_{1_h}$$ and $$\| \hat{q}_h\|_{0,\widehat{T}} \approx h \|q_h\|_{0,T}$$; thus $$\|\textbf{u}_h\|_{1_h} \preccurlyeq \| q_h\|_{0,T}$$. The proof is concluded with similar steps as at the end of the proof of Theorem 3.2. □ Proof of Theorem 2.2. For the proof we assume that $$k \ge 2$$. For the analysis of the lowest-order case, we refer to Lehrenfeld (2010). We construct a Fortin operator $${\it {\Pi}}_{F}$$ with   \begin{align}\label{fortin} b({\it {\Pi}}_{F} \textbf{u} -\textbf{u}, q_h) = 0 \quad \forall\, q_h \in Q_h \quad \textrm{and} \quad \| {{\it {\Pi}}_{F} \textbf{u} } \|_{1_h} \le c_{F} \|\textbf{u} \|_1, \end{align} (3.4) where $$c_{\mathrm{F}}$$ is a robust constant with respect to $$h$$ and $$k$$. Then, Theorem 2.2 follows from Boffi et al. (2013, Proposition 5.4.2). We define the operator $${\it {\Pi}}_{F}$$ as the sum of a low-order operator $${\it {\Pi}}_{F}^1$$ and a correction operator $${\it {\Pi}}_{F}^2$$. The first operator is the Fortin operator for the pair $${\left(\pmb{\mathscr{P}}^2(\mathscr{T}) \cap [C^0({\it {\Omega}})]^2 \right) \times \mathscr{P}^0(\mathscr{T})}$$, (see Boffi et al., 2013, Section 8.4), which is uniformly continuous in $$h$$. For the second operator, let $$(\textbf{w}_h^{\mathrm{T}},r_h^{\mathrm{T}}) \in \textbf{V}_{h,0}(T) \times Q_h(T)$$ be the solution of the local correction problem   \begin{align*} \begin{array}{rcll} a^{\mathrm{T}}(\textbf{w}_h^{\mathrm{T}},\textbf{v}_h) + b^{\mathrm{T}}(\textbf{v}_h, r_h^{\mathrm{T}}) &=& 0 \quad &\forall\, \textbf{v}_h \in \textbf{V}_{h,0}(T), \\ b^{\mathrm{T}}(\textbf{w}_h^{\mathrm{T}}, q_h) & = & \int_T \text{div}~{(\textbf{u} - {\it {\Pi}}_{F}^1 \textbf{u})} q_h ~\mathrm{d} x \quad & \forall\, q_h \in Q_h(T), \end{array} \end{align*} where $$a^{\mathrm{T}}(\cdot,\cdot)$$ and $$b^{\mathrm{T}}(\cdot,\cdot)$$ are the restrictions of the bilinear forms (2.1) and (2.2) on each triangle $$T \in \mathscr{T}$$. Note that this problem is solvable as $$\text{div}~{(\textbf{u} - {\it {\Pi}}_{F}^1 \textbf{u})} \perp \mathscr{P}^0(\mathscr{T})$$. Using Corollary 3.3 and the stability result of saddle point problems, we get   \begin{align*} \| {\textbf{w}^{\mathrm{T}}_h} \|_{1_h} \preccurlyeq \| \text{div}~{(\textbf{u} - {\it {\Pi}}_{F}^1 \textbf{u})} \|_0 \le c_{2} \| \textbf{u} \|_1, \end{align*} where the constant $$c_{2}$$ is independent of $$h$$ and $$k$$. Now we set $${\it {\Pi}}_{F}^2 \textbf{u} := \sum\limits_{T \in \mathscr{T}} \textbf{w}^{\mathrm{T}}_h$$ and see that $${\it {\Pi}}_{F} = {\it {\Pi}}_{F}^1 + {\it {\Pi}}_{F}^2$$ fulfils (3.4) as for all $$T \in \mathscr{T}$$ it holds that   \begin{align*} \text{div}~{{\it {\Pi}}_{F} \textbf{u}}|_T &= \text{div}~{{\it {\Pi}}_{F}^1 \textbf{u}}|_T + \text{div}~{{\it {\Pi}}_{F}^2 \textbf{u}}|_T = \text{div}~{{\it {\Pi}}_{F}^1 \textbf{u}}|_T + \sum\limits_{T \in \mathscr{T}} {\it {\Pi}}^{Q_h} (\text{div}~{\textbf{u}}|_T - \text{div}~{{\it {\Pi}}_{F}^1 \textbf{u}}|_T) = {\it {\Pi}}^{Q_h} \text{div}~{\textbf{u}}|_T, \\ \| {{\it {\Pi}}_{F} \textbf{u} } \|_{1_h} &\le \| {{\it {\Pi}}_{F}^1 \textbf{u} } \|_{1_h} + \| {{\it {\Pi}}_{F}^2 \textbf{u} } \|_{1_h} \preccurlyeq \| \textbf{u} \|_1 + \sum\limits_{T \in \mathscr{T}}\| {\textbf{w}^{\mathrm{T}}_h} \|_{1_h} \preccurlyeq \| \textbf{u} \|_1.\quad\square \end{align*} Remark 3.4 In the sense of Remark 2.4 one can choose the polynomial order in Theorem 3.1 as $$k_T$$; thus to show Theorem 2.2 one uses Theorem 3.2 and Corollary 3.3 with the proper order. Remark 3.5 As mentioned in Boffi et al. (2013, Remark 5.2.7), the existence of a Fortin operator $${\it {\Pi}}_{F}$$ can be used to construct an error estimation independent of the LBB constant $$\beta_\infty$$ of the infinite-dimensional Stokes problem (1.1). This may be an advantage for problems where $$\beta_\infty$$ is large. 4. Numerical examples Theorem 2.2 shows that the LBB constant is independent of the polynomial order $$k$$. Besides the analysis in Section 3, numerical tests also show this independence. In Table 1, one sees the different values of $$\beta$$ on the reference triangle $$\widehat{T}$$ and the reference quadrilateral $$\widehat{Q} = (0,1) \times (0,1)$$ for different polynomial orders $$k$$. Where the first line supports our results, the second line predicts that the polynomial robustness seems also to hold true in the case of quadrilaterals. The LBB constant was computed as the square root of the minimal positive eigenvalue of $$B G^{-1} B^{\mathrm{T}} M^{-1}$$, where $$B$$ is the discrete matrix of $$b(\cdot, \cdot)$$ and $$G$$ and $$M$$ are the Gramian matrices generated by $$(\textbf{V}_h, \| {\cdot} \|_{1_h})$$ and $$(Q_h, \| \cdot \|_0),$$ respectively. Table 1 LBB constant dependence of $$k$$ on the reference triangle $$\widehat{T}$$ and the reference quadrilateral $$\widehat{Q}$$ k  4  8  16  32  Triangle  0.167  0.190  0.201  0.205  Quadrilateral  0.305  0.313  0.315  0.315  k  4  8  16  32  Triangle  0.167  0.190  0.201  0.205  Quadrilateral  0.305  0.313  0.315  0.315  Table 1 LBB constant dependence of $$k$$ on the reference triangle $$\widehat{T}$$ and the reference quadrilateral $$\widehat{Q}$$ k  4  8  16  32  Triangle  0.167  0.190  0.201  0.205  Quadrilateral  0.305  0.313  0.315  0.315  k  4  8  16  32  Triangle  0.167  0.190  0.201  0.205  Quadrilateral  0.305  0.313  0.315  0.315  In the second example, we take a closer look at the stability of the method introduced in Section 2. For this, we solve problem (2.3) on the unit square $${\it {\Omega}}=(0,1)^2$$, where the exact solution is given by $$\textbf{u}:=\text{curl}~(\zeta) = (\zeta_{,y}, -\zeta_{,x})$$ with $$\zeta = \sin(x)^2\cos(x)^2$$, $$p=0$$, and set $$\nu=1$$. For this purpose, we use a triangulation with $$52$$ elements with $$h \approx 0.2$$ and a stabilization parameter $$\alpha = 4$$. In Table 2, we see the behaviour of the error $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$, the error of the best approximation $$\| {\textbf{u}-\textbf{u}_{\mathrm{best}}} \|_{1_h}$$ with $$\textbf{u}_{\mathrm{best}} := \arg \min\limits_{\textbf{v}_h \in \textbf{V}_h} \| {\textbf{u}-\textbf{v}_h} \|_{1_h}$$ and the ratio $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h} / \| {\textbf{u}-\textbf{u}_{\mathrm{best}}} \|_{1_h}$$. The values support that the discrete Stokes solution is close to the best approximation as the ratio is practically close to 1. Table 2 The error $$\| {\textbf{u}-\textbf{u}_{\mathrm{best}}} \|_{1_h}$$, $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$ and the ratio for different $$k$$ $$k$$  $$2$$  $$4$$  $$6$$  $$8$$  $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$  $$1.587$$  $$3.343e-02$$  $$3.430e-04$$  $$2.258e-06$$  $$\| {\textbf{u}-\textbf{u}_{\rm{best}}} \|_{1_h}$$  $$1.180$$  $$2.302e-02$$  $$2.309e-04$$  $$1.597e-06$$  Ratio  $$1.344$$  $$1.452$$  $$1.486$$  $$1.414$$  $$k$$  $$2$$  $$4$$  $$6$$  $$8$$  $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$  $$1.587$$  $$3.343e-02$$  $$3.430e-04$$  $$2.258e-06$$  $$\| {\textbf{u}-\textbf{u}_{\rm{best}}} \|_{1_h}$$  $$1.180$$  $$2.302e-02$$  $$2.309e-04$$  $$1.597e-06$$  Ratio  $$1.344$$  $$1.452$$  $$1.486$$  $$1.414$$  Table 2 The error $$\| {\textbf{u}-\textbf{u}_{\mathrm{best}}} \|_{1_h}$$, $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$ and the ratio for different $$k$$ $$k$$  $$2$$  $$4$$  $$6$$  $$8$$  $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$  $$1.587$$  $$3.343e-02$$  $$3.430e-04$$  $$2.258e-06$$  $$\| {\textbf{u}-\textbf{u}_{\rm{best}}} \|_{1_h}$$  $$1.180$$  $$2.302e-02$$  $$2.309e-04$$  $$1.597e-06$$  Ratio  $$1.344$$  $$1.452$$  $$1.486$$  $$1.414$$  $$k$$  $$2$$  $$4$$  $$6$$  $$8$$  $$\| {\textbf{u}-\textbf{u}_h} \|_{1_h}$$  $$1.587$$  $$3.343e-02$$  $$3.430e-04$$  $$2.258e-06$$  $$\| {\textbf{u}-\textbf{u}_{\rm{best}}} \|_{1_h}$$  $$1.180$$  $$2.302e-02$$  $$2.309e-04$$  $$1.597e-06$$  Ratio  $$1.344$$  $$1.452$$  $$1.486$$  $$1.414$$  All implementations of the numerical examples were performed with the finite element library NGSolve/Netgen (see Schöberl, 1997; 2014). 5. $$\boldsymbol{H^2}$$-stable polynomial extension In this section, we prove the existence of a polynomial-preserving $$H^2$$-stable extension operator. Note that in the two-dimensional case, the $$\text{curl}~$$ operator can be represented as a rotated gradient; thus on the boundary we have   \begin{align*} \text{curl}~{\textbf{u}_h} \cdot \textbf{n} = \nabla^\perp \textbf{u}_h \cdot \textbf{n} = \nabla \textbf{u}_h \cdot \tau \quad \textrm{on } \partial T. \end{align*} For the ease of notation and readability, we switch to the gradient and change the tangential and normal vector in Theorem 3.1. Furthermore, we skip the subscript $$h$$ of $$\textbf{u}_h$$; thus we write $$\textbf{u}$$ for a vectorial polynomial and show instead the following theorem. Theorem 5.1 ($$H^2$$-extension) For every $$k$$ there exists an operator $$\mathscr{E} : \pmb{\mathscr{P}}_{\tau}^k(\widehat{T}) \rightarrow \mathscr{P}^{k+1}(\widehat{T})$$, such that for $$\textbf{u} \in \pmb{\mathscr{P}}_{\tau}^k(\widehat{T})$$ it holds that   \begin{align} &\nabla {\mathscr{E}(\textbf{u})} \cdot \tau = \textbf{u} \cdot \tau \quad \textrm{on } \partial \widehat{T} \label{ext:h2ext:propA}, \\ \end{align} (5.1)  \begin{align} & \|(\textbf{u} - \nabla \mathscr{E}(\textbf{u})) \cdot \textbf{n} \|_{0,\partial \widehat{T}} \preccurlyeq \frac{1}{k} \| u \|_{1, \widehat{T}} \label{ext:h2ext:propB}, \\ \end{align} (5.2)  \begin{align} & \| \mathscr{E}(\textbf{u}) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}. \label{ext:h2ext:propC} \end{align} (5.3) Proof. The proof is provided in Section 5.5. □ 5.1 Literature and structure of this section Extension or lifting operators are a well-discussed topic, as they present the inverse map of trace operators that are known to be continuous, for example, in the scalar case from $$H^1(T)$$ onto $$H^{1/2}(\partial T)$$. The challenge then is to construct an operator that maps functions from $$H^{1/2}(\partial T)$$ into $$H^1(T)$$. Furthermore, a polynomial extension has the additional property that the operator maps a given polynomial on the boundary onto a polynomial on the element. The importance of such extensions has its origin mostly in the analysis of $$p$$- and $$hp$$- finite element methods (see, for example, Demkowicz & Babuška, 2003) and spectral methods and preconditioning (see, for example, Schöberl et al., 2008). Although the existence of polynomial $$H^1$$-, $$H(\operatorname{div})$$- and $$H(\text{curl}~op)$$-extensions are already well analysed, to the best of our knowledge, a stable polynomial $$H^2$$-extension, which is presented in this section, is the first result of this kind. Besides the application of the $$H^2$$-stable extension operator in the proof of the local LBB condition Theorem 3.2, this operator may be of interest, for example, to construct proper interpolation operators for fourth-order $$C^0$$-continuous DG methods (see, for example, Brenner et al., 2010, 2012). Before we present our results, we want to mention some literature, as many techniques we use are motivated by their accomplishments. The first work we consider is the pioneering article of Babuška & Suri (1987), which contains major ideas that are developed and adapted in later contributions including the technique of splitting the operator into primary and proper correction liftings. For polynomial-preserving extensions, we want to mention the work of Maday (1989) and Bernardi & Maday (1990) for two dimensions and the work of Belgacem (1994), Muñoz-Sola (1997) and, more recently, Ainsworth & Demkowicz (2009) for three dimensions. Finally, we want to mention the works of Demkowicz et al. (2008, 2009, 2012), where the polynomial extensions also provide a commuting diagram property for the corresponding spaces. For the proof of Theorem 5.1, we proceed in several steps. We start with an extension operator for the tangential values in Section 5.2 to provide the properties (5.1) and (5.2). After that we show in Section 5.3 how to proceed with the normal component. For this purpose, we assume that the polynomial on the edge has a zero of order 2 in the vertices to show an estimate in the $$H^{1/2}_{00}$$-norm. To lift the normal trace also for arbitrary polynomials, we show in Section 5.4 that the error, hence the part of the polynomial that does not satisfy the assumptions needed for the extension before, is bounded with a proper order of $$k$$ to show (5.2). Finally, we combine all estimates in Section 5.5 to prove Theorem 5.1. For the continuity estimates, we refer several times to Lemma 5.3 which is proved below. Remark 5.2 In this section, we only present the major techniques to show the $$H^2$$-continuity for the most difficult terms due to the similar structure of the rest. For a more detailed analysis, we refer to Lederer (2016, Chapter 6). Lemma 5.3 Let $$a \in C^1(E_1)$$ and $$u \in H^{1/2}(E_1)$$. For $$(x,y) \in \widehat{T}$$ we define $$g(x,y):=\int_0^1 a(s) u(x+sy) ~\mathrm{d} s$$. There holds   \begin{align*} \| g\|_{1, \widehat{T}} \preccurlyeq \| u \|_{1/2,E_1}. \end{align*} Proof. To derive the estimate of the $$L^2$$-norm we define for a fixed $$y \in [0,1]$$ the line $${l_y := \{(x,y): x \in [0,1-y] \}}$$ and use the Cauchy–Schwarz inequality to get   \begin{align*} \| g \|_{0,l_y}^2 &= \int_0^{1-y}\left(\int_0^1 a(s) u(x+sy) ~\mathrm{d} s \right)^2~\mathrm{d} x \preccurlyeq \int_0^{1-y} \int_0^1 u(x + sy)^2 ~\mathrm{d} s ~\mathrm{d} x. \end{align*} With the substitution $$t = x+sy$$ and Fubini’s theorem this leads to   \begin{align*} \| g \|_{0,\widehat{T}}^2 &=\int_0^1 \| g \|_{0,l_y}^2 ~\mathrm{d} y \preccurlyeq \int_0^1 \frac{1}{y} \int_0^{1-y} \int_x^{x+y} u(t)^2 ~\mathrm{d} t ~\mathrm{d} x ~\mathrm{d} y = \int_0^1 \frac{1}{y}\iint\limits_{\substack{0 \le x \le 1-y \\ x \le t \le x+y }} u(t)^2 ~\mathrm{d}(x,t) ~\mathrm{d} y \\ &\preccurlyeq \int_0^1 \frac{1}{y}\iint\limits_{\substack{0 \le t \le 1 \\ t-y \le x \le t }} u(t)^2 ~\mathrm{d}(x,t) ~\mathrm{d} y = \int_0^1 \frac{1}{y} \int_0^{1} u(t)^2 \underbrace{\int_{t-y}^{\mathrm{t}} ~\mathrm{d} x}_{=y} ~\mathrm{d} t ~\mathrm{d} y \preccurlyeq \int_0^1 \| u \|_{0,E_1}^2 ~\mathrm{d} y \preccurlyeq \| u \|_{1/2,E_1}^2. \end{align*} To bound the first-order derivatives, we use the real method of interpolation of spaces (see Bergh & Löfström, 1976) and Peetre’s K-functional technique (see Peetre, 1963). It is well known that we have an equivalent norm on the space $$H^{1/2}(E_1)$$ given by   \begin{align*} \| u \|_{1/2,E_1}^2 = \int_0^\infty y^{-2} |K(y,u)|^2 ~\mathrm{d} y \quad \textrm{with} \quad K(y,u) := \inf\limits_{\substack{u_0, u_1 \\u = u_0 + u_1 } } \sqrt{\|u_0 \|_{0,E_1}^2 + y^2 \|u_1'\|_{0,E_1}^2}. \end{align*} The derivative with respect to $$x$$ reads, $$g_{,x} = \int_0^1 u'(x+sy) a(s) ~\mathrm{d} s$$; thus with a similar estimate as for the $$L^2$$-norm we get $$\| g_{,x} \|_{0,l_y}^2 \preccurlyeq \| u'\|_{0,E_1}^2$$. Using integration by parts, we can also write   \begin{align*} g_{,x} = \underbrace{\frac{1}{y} \int_0^1 a'(s) u(x+sy) ~\mathrm{d} s}_{=:A_1} + \underbrace{\frac{1}{y}(a(1) u(x+y) - a(0)u(x))}_{=:B_2}. \end{align*} For $$A_1$$ we proceed similarly as for the $$L^2$$ estimate to get $$\| A_1 \|_{0,l_y}^2 \preccurlyeq \frac{1}{y^2} \| u \|^2_{0,E_1}$$. For $$B_1$$ we get with the substitution $$t= x+y$$,   \begin{align*} \| B_1 \|_{0,l_y}^2 &= \frac{1}{y^2} \int_0^{1-y} (a(1)u(x+y) - a(0)u(x))^2 ~\mathrm{d} x \preccurlyeq \frac{1}{y^2} \int_0^{1-y} u(x+y)^2 + u(x)^2 ~\mathrm{d} x \\ &\preccurlyeq \frac{1}{y^2} \int_y^{1} u(t)^2 ~\mathrm{d} t + \frac{1}{y^2}\int_0^{1-y} u(x)^2 ~\mathrm{d} x \preccurlyeq \frac{1}{y^2} \| u \|^2_{0,E_1}. \end{align*} Combining both estimates, so $$\| g_{,x} \|_{0,l_y} \preccurlyeq 1/y \| u\|_{0,E_1}$$ and $$\| g_{,x} \|_{0,l_y} \preccurlyeq \| u'\|_{0,E_1}$$, leads to   \begin{align*} \| g_{,x} \|^2_{0,l_y} &\le \inf\limits_{\substack{u_0, u_1 \\u = u_0 + u_1 } } \| u_0 \|^2_{0,l_y} + \| u_1 \|^2_{0,l_y} \preccurlyeq \inf\limits_{\substack{u_0, u_1 \\u = u_0 + u_1 } } \frac{1}{y^2} \| u_0 \|^2_{0,E_1} + \| u_1' \|^2_{0,E_1} \preccurlyeq \frac{1}{y^2} K(y,u)^2, \end{align*} and thus   \begin{align*} \| g_{,x} \|_{0,\widehat{T}}^2 &= \int_0^1 \| g \|_{0,l_y}^2 ~\mathrm{d} y \preccurlyeq \int_0^1 \frac{1}{y^2} K(y,u)^2 ~\mathrm{d} y \preccurlyeq \int_0^\infty \frac{1}{y^2} K(y,u)^2 ~\mathrm{d} y = \|u \|_{1/2,E_1}^2. \end{align*} A similar estimation for the derivative with respect to $$y$$ concludes the proof. □ 5.2 Tangential extension Theorem 5.4 For every $$k$$ there exists an operator $$\mathscr{E}^{\tau} : \pmb{\mathscr{P}}_{\tau}^k(\widehat{T}) \rightarrow \mathscr{P}^{k+1}(\widehat{T})$$ such that for $$\textbf{u} \in \pmb{\mathscr{P}}_{\tau}^k(\widehat{T})$$ it holds that   \begin{align} &\nabla {\mathscr{E}^{\tau}(\textbf{u})} \cdot \tau = \textbf{u} \cdot \tau \quad \textrm{on } \partial \widehat{T}, \label{ext:h2exttang:propA} \\ \end{align} (5.4)  \begin{align} & \| \mathscr{E}^{\tau}(\textbf{u}) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}. \label{ext:h2exttang:propC} \end{align} (5.5) Proof. For the proof we proceed in several steps. First, we construct an extension from the lower edge $$E_1$$ onto the triangle $$\widehat{T}$$ without any restrictions on the values on the two other edges $$E_2$$ and $$E_3$$. After that, we construct two more extensions with proper values on the other edges and combine them afterwards to define $$\mathscr{E}^{\tau}$$. Step 1. For $$u_\tau(x) := \textbf{u}(x,0) \cdot \tau$$, where $$\tau := (1,0)^{\mathrm{t}}$$ is the tangential vector on the lower edge $$E_1$$, we define   \begin{align} \label{ext::tangextdefstepone} \psi(x) := \int_0^x u_\tau(s) ~\mathrm{d} s \qquad \textrm{and} \qquad \mathscr{E}^{\tau}_1(\textbf{u})(x,y) := \int_0^1 \psi(x + sy) ~\mathrm{d} s. \end{align} (5.6) Note that the derivatives read   \begin{align*} \mathscr{E}^{\tau}_{1,x}(\textbf{u}) = \int_0^1 \textbf{u}t(x+sy) ~\mathrm{d} s \quad \textrm{and} \quad \mathscr{E}^{\tau}_{1,y}(\textbf{u}) = \int_0^1 \textbf{u}t(x+sy)s ~\mathrm{d} s, \end{align*} thus   \begin{align*} \nabla \mathscr{E}^{\tau}_1(\textbf{u}) \cdot \tau = \int_0^1 \textbf{u}t(x) ~\mathrm{d} s = \textbf{u} \cdot \tau \quad \textrm{on } E_1. \end{align*} Next we show the $$H^2$$-continuity, so $$\| \mathscr{E}^{\tau}_1(\textbf{u}) \|_{0,\widehat{T}}^2 + \| \nabla \mathscr{E}^{\tau}_1(\textbf{u}) \|_{1,\widehat{T}}^2 \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}^2$$. For the estimate of the $$L^2$$-norm, we proceed similarly to the first steps of the proof of Lemma 5.3, thus we estimate the $$L^2$$-norm on a line $$l_y$$ using the Cauchy–Schwarz inequality and Fubini’s theorem and then integrate from $$0$$ to $$1$$ with respect to $$y$$. For the first-order derivatives, we can use Lemma 5.3 to get   \begin{align*} \| \mathscr{E}^{\tau}_{1,x}(\textbf{u}) \|_{1,\widehat{T}}^2 \preccurlyeq \| \textbf{u}t \|_{1/2,E_1}^2 \quad \text{and} \quad \| \mathscr{E}^{\tau}_{1,y}(\textbf{u}) \|_{1,\widehat{T}}^2 \preccurlyeq \| \textbf{u}t \|_{1/2,E_1}^2. \end{align*} Altogether we have   \begin{align} \label{ext::tangextsteponeest} \| \mathscr{E}^{\tau}_{1}(\textbf{u}) \|_{2,\widehat{T}}^2 \preccurlyeq \| \textbf{u}t \|_{1/2,E_1}^2 \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}^2. \end{align} (5.7) Step 2. In the next step, we want to construct an extension from the lower edge $$E_1$$ with the restriction that the gradient of the extension has zero tangential values on the edge $$E_2$$. Similarly to before, we define for $$u_\tau(x) := \textbf{u}(x,0) \cdot \tau$$ on $$E_1$$,   \begin{align*} \psi(x) := \int_0^x u_\tau(s) ~\mathrm{d} s - \overline{\psi(x)} \qquad \textrm{with} \qquad \overline\psi(x) := \int_0^1 u_\tau(s) ~\mathrm{d} s, \end{align*} and   \begin{align} \label{ext::tangextdefsteptwo} \mathscr{E}^{\tau}_2(\textbf{u})(x,y) := \int_0^1 \psi(x + sy) ~\mathrm{d} s - \frac{y}{1-x} \int_0^1 \psi(x+s(1-x)) ~\mathrm{d} s. \end{align} (5.8) Using integration by parts for the second integral we furthermore use the representation   \begin{align} \label{ext::tangextdefsteptwoalt} \mathscr{E}^{\tau}_2(\textbf{u})(x,y) &= \int_0^1 \psi(x + sy) ~\mathrm{d} s + y \int_0^1 \psi'(x+s(1-x)) s ~\mathrm{d} s - \overbrace{\psi(1)}^{=0} \nonumber\\ &= \int_0^1 \psi(x + sy) ~\mathrm{d} s + y \int_0^1 u_\tau(x+s(1-x)) s ~\mathrm{d} s. \end{align} (5.9) On the edges $$E_1, E_2$$ we have   \begin{align*} \mathscr{E}^{\tau}_2(\textbf{u})|_{E_2} = 0 \Rightarrow \nabla \mathscr{E}^{\tau}_2(\textbf{u}) \cdot \tau|_{E_2} = 0 \quad \textrm{and} \quad \mathscr{E}^{\tau}_2(\textbf{u}) \cdot \tau|_{E_1} = \frac{\partial \psi}{\partial x} = u_\tau = \textbf{u} \cdot \tau|_{E_1}. \end{align*} For the $$H^2$$-continuity we proceed as before. Note that the first part of $$\mathscr{E}^{\tau}_2(\textbf{u})$$ is the same as in Step 1, so we consider only the second integral from (5.9); thus the correction term is   \begin{align*} \mathscr{E}^{\tau,c}_2(\textbf{u}) := y \int_0^1 u_\tau(x+s(1-x)) s ~\mathrm{d} s. \end{align*} We present only the estimate for the second-order derivative with respect to $$x$$ as some new techniques have to be used there. The rest follows with similar steps to before. Using integration by parts for $$\textbf{u}t'$$ we observe   \begin{align*} \mathscr{E}^{\tau,c}_{2,xx}(\textbf{u}) &= \frac{y}{(1-x)^2} \int_0^1 \textbf{u}t(x+s(1-x))(2s-1) ~\mathrm{d} s \\ &\quad{} + \frac{y}{1-x} \int_0^1 \textbf{u}t'(x+s(1-x))(2s-1) ~\mathrm{d} s \\ &= \frac{y}{(1-x)^2} \int_0^1 \textbf{u}t(x+s(1-x))(2s-1) ~\mathrm{d} s \\ & \quad{}+ \frac{y}{(1-x)^2} \int_0^1 \textbf{u}t(x+s(1-x))(4s-3) ~\mathrm{d} s + \frac{y}{(1-x)^2} \textbf{u}t(x), \end{align*} and by splitting the second integral into two terms, finally   \begin{align*} \mathscr{E}^{\tau,c}_{2,xx}(\textbf{u}) = \underbrace{\frac{y}{(1-x)^2} \int_0^1 \textbf{u}t(x+s(1-x))(6s-3) ~\mathrm{d} s}_{=:A_2} + \underbrace{ \frac{y}{(1-x)^2} \int_0^1 \textbf{u}t(x) - \textbf{u}t(x+s(1-x)) ~\mathrm{d} s}_{=: B_2}. \end{align*} We start with the estimate of $$A_2$$. For this, note that the polynomial $$6s-3$$ has a zero integral value on $$E_1$$, so we subtract the mean value $$\overline{\textbf{u}t(x)} := \frac{1}{1-x} \int_x^1 \textbf{u}t(s) ~\mathrm{d} s$$, and get   \begin{align*} A_2 &= \frac{y}{(1-x)^2} \int_0^1 \textbf{u}t(x+s(1-x))(6s-3) ~\mathrm{d} s = \frac{y}{(1-x)^2} \int_0^1 (\textbf{u}t(x+s(1-x)) - \overline{\textbf{u}t(x)}) (6s-3) ~\mathrm{d} s. \end{align*} Using the Cauchy–Schwarz inequality and $$\frac{y}{1-x} \le 1$$ on $$\widehat{T}$$ and the substitution $$t = x + s(1-x)$$ leads to   \begin{align*} \| A_2 \|^2_{0,\widehat{T}} &= \int_0^1 \int_0^{1-x} A_2^2 ~\mathrm{d} y ~\mathrm{d} x \\ & \preccurlyeq \int_0^1 \frac{1}{(1-x)^2} \left(\int_0^1 (\textbf{u}t(x+s(1-x)) - \overline{\textbf{u}t(x)}) (6s-3) ~\mathrm{d} s \right)^2 \int_0^{1-x} ~\mathrm{d} y ~\mathrm{d} x\\ & \preccurlyeq \int_0^1 \frac{1}{1-x} \int_0^1 (\textbf{u}t(x+s(1-x)) - \overline{\textbf{u}t(x)})^2 ~\mathrm{d} s ~\mathrm{d} x \\ & = \int_0^1 \frac{1}{(1-x)^2} \int_x^1 (\textbf{u}t(t) - \overline{\textbf{u}t(x)})^2 ~\mathrm{d} t ~\mathrm{d} x. \end{align*} For the next step, we use the following identity for the inner integral:   \begin{align*} \int_x^1 (\textbf{u}t(t) - \overline{\textbf{u}t(x)})^2 ~\mathrm{d} t = \frac{1}{2(1-x)} \int_x^1 \int_x^1 (\textbf{u}t(t) - \textbf{u}t(s))^2 ~\mathrm{d} t ~\mathrm{d} s. \end{align*} Similarly to Step 1, we use Fubini’s theorem to handle the $$(1-x)^3$$ in the denominator, so   \begin{align*} \| A_2 \|^2_{0,\widehat{T}} &\preccurlyeq \int_0^1 \frac{1}{(1-x)^3} \int_x^1 \int_x^1 (\textbf{u}t(t) - \textbf{u}t(s))^2 ~\mathrm{d} t ~\mathrm{d} s ~\mathrm{d} x \\ & \preccurlyeq \iiint\limits_{\substack{x \le s\\ x \le t \\ 0 \le s,t \le 1 }} \frac{(\textbf{u}t(t) - \textbf{u}t(s))^2 }{(1-x)^3} ~\mathrm{d} (s,t,x) \preccurlyeq \int_0^1 \int_0^1 \int_0^{\min(s,t)} \frac{1}{(1-x)^3} ~\mathrm{d} x (\textbf{u}t(t) - \textbf{u}t(s))^2 ~\mathrm{d} t ~\mathrm{d} s. \end{align*} Without loss of generality assuming $$s < t$$ and using $$1-s \ge t-s$$ for $$s < t < 1$$, it follows for the inner integral that   \begin{align*} \int_0^{\min(s,t)} \frac{1}{(1-x)^3}~\mathrm{d} x = \left(\frac{1}{1-\min(s,t) } \right) ^2 - 1 \le \left(\frac{1}{1-s} \right)^2 \le \left (\frac{1}{t-s} \right)^2. \end{align*} Together with the definition of the $$H^{1/2}$$-seminorm (see, Sobolev–Slobodeckij norm, for example, in Schwab, 1998, Theorem A.7) we get   \begin{align*} \| A_2 \|^2_{0,\widehat{T}} & \preccurlyeq \int_0^1 \int_0^1 \frac{(\textbf{u}t(t) - \textbf{u}t(s))^2}{(t-s)^2} ~\mathrm{d} t ~\mathrm{d} s = | \textbf{u}t|^2_{1/2,E_1} \preccurlyeq \| \textbf{u} \|^2_{1,\widehat{T}}. \end{align*} For $$B_2$$, we proceed similarly using the Cauchy–Schwarz inequality, Fubini’s theorem, the substitution $$t = x + s(1-x)$$ and $$1-x \ge t-x$$ for $$x < t < 1$$; thus   \begin{align*} \| B_2 \|^2_{0,\widehat{T}} &= \int_0^1 \int_0^{1-x} B_2^2 ~\mathrm{d} y ~\mathrm{d} x \preccurlyeq \int_0^1 \frac{1}{1-x} \int_0^1 \left(\textbf{u}t(x) - \textbf{u}t(x+s(1-x)) \right)^2 ~\mathrm{d} s ~\mathrm{d} x \\ & = \int_0^1 \frac{1}{(1-x)^2} \int_x^1 \left(\textbf{u}t(x) - \textbf{u}t(t) \right)^2 ~\mathrm{d} t ~\mathrm{d} x = \iint\limits_{\substack{x \le t \le 1\\0\le x \le 1 }} \frac{\left(\textbf{u}t(x) - \textbf{u}t(t) \right)^2}{(1-x)^2} ~\mathrm{d}(x,t)\\ &\preccurlyeq \iint\limits_{\substack{0 \le x \le t \\0\le t \le 1 }} \frac{\left(\textbf{u}t(x) - \textbf{u}t(t) \right)^2}{(x-t)^2} ~\mathrm{d}(x,t) \preccurlyeq \int_0^1 \int_0^1 \frac{ \left(\textbf{u}t(x) - \textbf{u}t(t) \right)^2}{(x-t)^2} ~\mathrm{d} x ~\mathrm{d} t \preccurlyeq | \textbf{u}t|^2_{1/2,E_1} \preccurlyeq \| \textbf{u} \|^2_{1,\widehat{T}}. \end{align*} Altogether we have $$\| \mathscr{E}^{\tau,c}_{2,xx}(\textbf{u})\|_{0,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}$$ and assuming similar estimates for the other derivatives, finally   \begin{align} \label{ext::exttangtwoest} \| \mathscr{E}^{\tau}_{2}(\textbf{u})\|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}. \end{align} (5.10) Step 3. For this step, we assume that the tangential values of the input function $$\textbf{u}$$ are zero on the edges $$E_2$$ and $$E_3$$ and that it has a zero tangential integral value; thus   \begin{align} \label{ext::assumptionsstepthree} \int_{\partial \widehat{T}} \textbf{u} \cdot \tau = 0 \quad \textrm{and} \quad \textbf{u} \cdot \tau|_{E_2} = \textbf{u} \cdot \tau|_{E_3} = 0. \end{align} (5.11) We set $$u_\tau(x) := \textbf{u}(x,0) \cdot \tau$$ on $$E_1$$ and $$\psi(x) := \int_0^x u_\tau(s) ~\mathrm{d} s$$, to define   \begin{align} \label{ext::tangextdefstepthree} \mathscr{E}^{\tau}_{3}(\textbf{u}) &:= \int_0^1 \psi(x+sy) ~\mathrm{d} s - \frac{y}{1-x} \int_0^1 \psi(x+s(1-x)) ~\mathrm{d} s \nonumber\\ &\quad\ \underbrace{- \frac{y}{x+y} \int_0^1 \psi(s(x+y)) ~\mathrm{d} s}_{=: A_3} + \underbrace{ y \int_0^1 \psi(s) ~\mathrm{d} s}_{=: B_3}. \end{align} (5.12) As in Step 2, we observe $$\nabla \mathscr{E}^{\tau}_{3}(\textbf{u}) \cdot \tau |_{E_2} = \nabla \mathscr{E}^{\tau}_{3}(\textbf{u}) \cdot \tau |_{E_3} = 0$$ and $$\nabla \mathscr{E}^{\tau}_{3}(\textbf{u}) \cdot \tau |_{E_1} = \textbf{u}t$$. For the $$H^2$$-continuity we have only to estimate the terms $$A_3$$ and $$B_3$$ as the other terms are the same as in Steps 1 and 2. For this purpose, note that due to the assumptions on $$\textbf{u}$$, we have $$\psi(1) = \psi(0) = 0$$. Using the identity $$\frac{~\mathrm{d}}{~\mathrm{d} s} \psi(s(x+y)) \frac{1}{x+y} = \psi'(s(x+y))$$ and integration by parts we write   \begin{align} \label{ext::tangextdefstepthreealtone} A_3 &= - \frac{y}{x+y} \int_0^1 \psi(s(x+y)) ~\mathrm{d} s = y \int_0^1 \psi'(s(x+y))(s-1) ~\mathrm{d} s = y \int_0^1 \textbf{u}t(s(x+y))(s-1) ~\mathrm{d} s, \end{align} (5.13) and   \begin{align} \label{ext::tangextdefstepthreealttwo} B_3 = y \int_0^1 \psi(s) ~\mathrm{d} s = -y \int_0^1 \psi'(s)s ~\mathrm{d} s = -y \int_0^1 \textbf{u}t(s)s ~\mathrm{d} s. \end{align} (5.14) Using these representations and the same techniques as in Steps 2 and 3, we estimate the $$H^2$$-norm of $$A_3$$ and $$B_3$$ to show   \begin{align} \label{ext::exttangthreeest} \| \mathscr{E}^{\tau}_{3}(\textbf{u}) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}. \end{align} (5.15) Step 4. We finally combine the three extensions to show Theorem 5.4. For this purpose, let us assume we have a given function $$\textbf{u} \in \pmb{\mathscr{P}}_{\tau}^k(\widehat{T})$$ with $$\int_{\partial \widehat{T}} \textbf{u} \cdot \tau = 0$$. We first introduce two mappings from the reference triangle $$\widehat{T}$$ to itself by   \begin{align*} F_2: (x,y) \mapsto (x,1-x-y) \quad \textrm{and} \quad F_3: (x,y) \mapsto (y,x), \end{align*} where $$F_2$$ maps the values from $$E_2$$ to $$E_1$$ and vice versa, and the mapping $$F_3$$ from $$E_3$$ to $$E_1$$ and vice versa. Furthermore, we define for $$F_2$$ and $$F_3$$ the corresponding covariant mappings $$\mathscr{C}_2$$ and $$\mathscr{C}_3$$. Using those transformations we now introduce the extensions from Steps 2 and 3 also from the other edges; thus we define   \begin{align*} \tilde{\mathscr{E}^{\tau}_2}(\textbf{u})(x,y) = \mathscr{E}^{\tau}_2(\mathscr{C}_2\textbf{u})(F_2(x,y)) \quad \textrm{and} \quad \tilde{\mathscr{E}^{\tau}_3}(\textbf{u})(x,y) = \mathscr{E}^{\tau}_3(\mathscr{C}_3\textbf{u})(F_3(x,y)), \end{align*} with the properties   \begin{align} &(\nabla \tilde{\mathscr{E}^{\tau}_2}(\textbf{u})\cdot \tau)|_{E_2} = (\textbf{u} \cdot \tau)|_{E_2}, \quad (\nabla \tilde{\mathscr{E}^{\tau}_2}(\textbf{u})\cdot \tau)|_{E_1} = 0, \quad \| \tilde{\mathscr{E}^{\tau}_2}(\textbf{u}) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}, \label{ext::exttangtwotildeest}\\ \end{align} (5.16)  \begin{align} &(\nabla \tilde{\mathscr{E}^{\tau}_3}(\textbf{u})\cdot \tau)|_{E_3} = (\textbf{u} \cdot \tau)|_{E_3}, \quad (\nabla \tilde{\mathscr{E}^{\tau}_3}(\textbf{u})\cdot \tau)|_{E_1} = (\nabla \tilde{\mathscr{E}^{\tau}_3}(\textbf{u})\cdot \tau)|_{E_2}= 0, \quad \| \tilde{\mathscr{E}^{\tau}_3}(\textbf{u}) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}},\label{ext::exttangthreetildeest} \end{align} (5.17) which follow from the proper transformation of tangential values due to the use of the covariant transformations $$\mathscr{C}_2$$ and $$\mathscr{C}_3$$ and estimates (5.10) and (5.15). We define the final extension by setting $$e_1 := \mathscr{E}^{\tau}_1(\textbf{u})$$, $$e_2 := e_1 + \tilde{\mathscr{E}^{\tau}_2}(\textbf{u}- \nabla e_1)$$ and $$\mathscr{E}^{\tau}(\textbf{u}) := e_2 + \tilde{\mathscr{E}^{\tau}_3}(\textbf{u} - \nabla e_2)$$. Note that due to Green’s theorem, the surface integral over the boundary of the reference element $$\partial \widehat{T}$$ of $$(\textbf{u} - \nabla e_2) \cdot \tau$$ is equal to zero, and due to the properties of $$\tilde{\mathscr{E}^{\tau}_2}$$ and $$\mathscr{E}^{\tau}_1$$, the tangential values on $$E_1$$ and $$E_2$$ also vanish; thus assumptions (5.11) are fulfilled. Using (5.16) and (5.17) we observe   \begin{align*} \nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \tau|_{E_1} &= \nabla e_2 \cdot \tau|_{E_1} + \underbrace{\nabla \tilde{\mathscr{E}^{\tau}}_3(\textbf{u} - \nabla e_2) \cdot \tau|_{E_1}}_{=0} = \nabla e_1 \cdot \tau|_{E_1} + \underbrace{ \nabla \tilde{\mathscr{E}^{\tau}}_2(\textbf{u} - \nabla e_1) \cdot \tau|_{E_1}}_{=0} = (\textbf{u} \cdot \tau)|_{E_1}, \end{align*} and similarly $$\nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \tau|_{E_2} = \textbf{u} \cdot \tau|_{E_2}$$ and $$\nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \tau|_{E_3} = \textbf{u} \cdot \tau|_{E_3}$$; thus property (5.4) is fulfilled. With (5.16), (5.17) and (5.7) and the linearity of the operators, we furthermore have $$\| \mathscr{E}^{\tau}(\textbf{u}) \|_{2,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}$$; thus also the $$H^2$$-continuity (5.5) holds true. It remains to show that $$\mathscr{E}^{\tau}(\textbf{u}) \in \mathscr{P}^{k+1}(\widehat{T})$$. First, note that from $$\textbf{u} \in [\mathscr{P}^k(\widehat{T})]^2$$, it follows that $$\textbf{u} \cdot \tau \in \mathscr{P}^k(\partial \widehat{T})$$. Looking at the definition of the first extension (5.6) we see that we integrate $$u$$ from $$0$$ to $$x$$ to define $$\psi$$; thus, here, we increase the order by 1 resulting in $$\mathscr{E}^{\tau}_1(\textbf{u}) \in \mathscr{P}^{k+1}(\widehat{T})$$. For the other two extensions (5.8) and (5.12), this may not hold due to the fractional factors, but using the alternative representations (5.9), (5.13) and (5.14) we see that also the corrections are polynomial liftings and it follows that $$\mathscr{E}^{\tau}(\textbf{u}) \in \mathscr{P}^{k+1}(\widehat{T})$$. □ 5.3 Normal extension Theorem 5.5 For every $$k$$ there exists an operator $$\mathscr{E}^{\textbf{n}} : \mathscr{P}_{00}^k(E_1) \rightarrow \mathscr{P}^{k+1}(\widehat{T})$$ such that for $$u \in \mathscr{P}_{00}^k(E_1)$$ it holds that   \begin{align} &\mathscr{E}^{\textbf{n}}(u) = 0 \quad\hspace{20mm} \textrm{on } \partial \widehat{T}, \label{ext:h2extnorm:propA} \\ \end{align} (5.18)  \begin{align} &\nabla \mathscr{E}^{\textbf{n}}(u) \cdot \textbf{n} = u \quad \textrm{on } E_1, \label{ext:h2extnorm:propB} \\ \end{align} (5.19)  \begin{align} &\nabla \mathscr{E}^{\textbf{n}}(u) \cdot \textbf{n} = 0 \quad \textrm{on } \partial \widehat{T} \setminus E_1, \label{ext:h2extnorm:propD} \\ \end{align} (5.20)  \begin{align} & \| \mathscr{E}^{\textbf{n}}(u) \|_{2,\widehat{T}} \preccurlyeq \| u \|_{1/2^*,E _1}. \label{ext:h2extnorm:propC} \end{align} (5.21) Proof. Similarly to the proof of Theorem 5.4 we proceed in several steps. We first construct an extension from the lower edge and correct the value and the normal derivative on the other two edges afterwards. Step 1. We start with the first extension, so we define   \begin{align} \label{ext::extnormstepone} \mathscr{E}^{\textbf{n}}_1(u) := -y \int_0^1a(s)u(x+sy) ~\mathrm{d} s \quad \textrm{with} \quad a(s):=6s(1-s). \end{align} (5.22) It immediately follows that $$\mathscr{E}^{\textbf{n}}_1(u)|_{E_1} = 0$$. For the derivatives, we observe due to $$a(0) = a(1) = 0$$, $$\int_0^1a(s) ~\mathrm{d} s = 1$$, $$\int_0^1a'(s)s ~\mathrm{d} s = -1$$ and using integration by parts with $$u'(x+sy) = \frac{1}{y}\frac{d}{ds} u(x+sy)$$, that   \begin{align*} \mathscr{E}^{\textbf{n}}_{1,x}(u) &= -y \int_0^1 a(s) u'(x+sy) ~\mathrm{d} s = \int_0^1 a'(s) u(x+sy) ~\mathrm{d} s, \\ \mathscr{E}^{\textbf{n}}_{1,y}(u) &= -\int_0^1 a(s) u(x+sy) ~\mathrm{d} s - y\int_0^1 a(s) u'(x+sy)s ~\mathrm{d} s = \int_0^1 a'(s) u(x+sy)s ~\mathrm{d} s; \end{align*} thus   \begin{align*} \nabla \mathscr{E}^{\textbf{n}}_1(u) \cdot n|_{E_1} =- \mathscr{E}^{\textbf{n}}_{1,y}(u)|_{E_1} = u. \end{align*} The $$H^2$$-continuity estimate follows by applying Lemma 5.3 for the derivatives and the Cauchy–Schwarz inequality for the other terms; thus we have   \begin{align*} \| \mathscr{E}^{\textbf{n}}_1(u) \|_{2,\widehat{T}} \preccurlyeq \| u \|_{1/2,E_1} \preccurlyeq \| u \|_{1/2^*,E_1}. \end{align*} Step 2. In this step, we want to correct the values and the normal derivative on the second edge $$E_2$$ without changing the values and the normal derivative on the bottom edge $$E_1$$. For this purpose, we introduce the polynomials $$b(s) = 3s^2-2s^3$$ and $$c(s)=s^3-s^2$$, with the properties   \begin{align*} b(0) = b'(0) = b'(1) = 0,\quad b(1) = 1\quad \textrm{and} \quad c(0) = c'(0) = c(1) = 0,\quad c'(1) = 1, \end{align*} and use them as blending coefficients to define   \begin{align*} \mathscr{E}^{\textbf{n}}_2(u)(x,y):= \mathscr{E}^{\textbf{n}}_1(u)(x,y) - b\left(\frac{y}{1-x}\right)\mathscr{E}^{\textbf{n}}_1(u) (x,1-x) - c\left(\frac{y}{1-x}\right)(1-x) \mathscr{E}^{\textbf{n}}_{1,y}(u)(x,1-x). \end{align*} The idea is that the second term corrects the values and the last term corrects the normal derivative on the edge $$E_2$$. Indeed, we observe on $$E_1$$ as $$y=0$$ and on $$E_2$$ as $$y = 1-x$$ that   \begin{align} \mathscr{E}^{\textbf{n}}_2(u)|_{E_1}&= \mathscr{E}^{\textbf{n}}_1(u)|_{E_1}- b(0) \mathscr{E}^{\textbf{n}}_1(u)(x,1-x) - c(0) \mathscr{E}^{\textbf{n}}_{1,y}(u)(x,1-x) = 0,\nonumber \\ \mathscr{E}^{\textbf{n}}_2(u)|_{E_2}&= \mathscr{E}^{\textbf{n}}_1(u)|_{E_1}- b(1) \mathscr{E}^{\textbf{n}}_1(u)(x,1-x) - c(1) \mathscr{E}^{\textbf{n}}_{1,y}(u)(x,1-x) = 0.\label{ext::extnormstep2onedgetwoB} \end{align} (5.23) The derivative with respect to $$y$$ reads   \begin{align*} \mathscr{E}^{\textbf{n}}_{2,y}(u) &= \mathscr{E}^{\textbf{n}}_{1,y}(u) - \frac{1}{1-x}b'\left(\frac{y}{1-x}\right) \mathscr{E}^{\textbf{n}}_1(u) (x,1-x) - \frac{1}{1-x}c'\left(\frac{y}{1-x}\right)(1-x) \mathscr{E}^{\textbf{n}}_{1,y}(u) (x,1-x). \end{align*} On $$E_1$$ it follows, due to $$b'(0) = c'(0) = 0$$, that $$\mathscr{E}^{\textbf{n}}_{2,y}(u)|_{E_1}= \mathscr{E}^{\textbf{n}}_{1,y}(u)|_{E_1} = u$$, so the normal derivative has not changed in the second step. In a similar way, we also observe that $$\mathscr{E}^{\textbf{n}}_{2,y}(u)|_{E_2}=0$$. Now note that due to the constant zero value on the edge, see equation (5.23), we derive that the tangential derivative on the edge $$\nabla \mathscr{E}^{\textbf{n}}_2(u) \cdot \tau$$ on $$E_2$$ has to be zero. As $$\nabla \mathscr{E}^{\textbf{n}}_2(u) \cdot \tau|_{E_2} = 0 \Leftrightarrow - \mathscr{E}^{\textbf{n}}_{2,x}(u)|_{E_2} = \mathscr{E}^{\textbf{n}}_{2,y}(u)|_{E_2}$$ and $$\mathscr{E}^{\textbf{n}}_{2,y}(u)|_{E_2} = 0$$, it follows that $$\nabla \mathscr{E}^{\textbf{n}}_{2}(u) \cdot n|_{E_2} = 0$$, so the correction term induced a zero normal derivative on the second edge $$E_2$$. The $$H^2$$-estimate remains. As the first term of $$\mathscr{E}^{\textbf{n}}_2(u)$$ has already been analysed in the first step, we focus on just the correction terms. We start with first term $$A_4:= b(\frac{y}{1-x})\mathscr{E}^{\textbf{n}}_1(u) (x,1-x)$$ and the estimate for the $$y$$ derivative   \begin{align*} A_{4,y} = 6\frac{y}{1-x}\left(1-\frac{y}{1-x}\right) \int_0^1 u(x+s(1-x))a(s) ~\mathrm{d} s. \end{align*} Using the Cauchy–Schwarz inequality, $$\frac{y}{1-x} \le 1$$ on $$\widehat{T}$$ and the substitution $$t := x+s(1-x)$$, we get   \begin{align*} \| A_{4,y} \|_{0,\widehat{T}}^2 &\preccurlyeq \int_0^1\int_0^{1-x}\int_0^1 u(x+s(1-x))^2 ~\mathrm{d} s ~\mathrm{d} y ~\mathrm{d} x = \int_0^1\int_0^{1-x} \frac{1}{1-x}\int_x^1 u(t)^2 ~\mathrm{d} t ~\mathrm{d} y ~\mathrm{d} x \\ &= \int_0^1 \int_x^1 u(t)^2 ~\mathrm{d} t ~\mathrm{d} x \preccurlyeq \| u \|_{0,E_1}^2 \preccurlyeq \| u \|_{1/2^*,E_1}^2, \end{align*} and in a similar way we also bound $$\| A_{4,x} \|_{0,\widehat{T}}$$ and $$\| A_4 \|_{0,\widehat{T}}$$. The crucial point in this estimate was that we were able to use property $$\frac{y}{1-x} \le 1$$ twice, thus there were no bad coefficients any more. The estimates of the second-order derivatives are a little bit more tricky, as there remain some fractions with singularities. We start with the second-order derivative with respect to $$y$$ given by   \begin{align*} A_{4,yy} = 6\frac{1}{1-x}\left(1-\frac{2y}{1-x}\right) \int_0^1 u(x+s(1-x))a(s) ~\mathrm{d} s. \end{align*} The idea is to use Fubini’s theorem,   \begin{align*} \| A_{4,yy} \|_{0,\widehat{T}}^2 &\preccurlyeq \int_0^1\int_0^{1-x} \frac{1}{(1-x)^2}\int_0^1 u(x+s(1-x))^2 ~\mathrm{d} s ~\mathrm{d} y ~\mathrm{d} x \\ &= \int_0^1 \frac{1}{(1-x)^2}\int_x^1 u(x+s(1-x))^2 ~\mathrm{d} s ~\mathrm{d} y ~\mathrm{d} x = \iint\limits_{\substack{0 \le x \le 1\\ x \le t}} \frac{u(t)^2}{(1-x)^2}~\mathrm{d}(x,t) \\ &= \iint\limits_{\substack{0 \le t \le 1\\ t \le x}} \frac{u(t)^2}{(1-x)^2}~\mathrm{d}(x,t) = \int_0^1\int_0^{\mathrm{t}} \frac{1}{(1-x)^2} ~\mathrm{d} x~ u(t)^2 ~\mathrm{d} t \\ &= \int_0^1 \frac{1}{1-t}u(t)^2 ~\mathrm{d} t \preccurlyeq \| u \|_{0^*,E_1}^2 \preccurlyeq \| u \|_{1/2^*,E_1}^2. \end{align*} With the techniques just presented and similar estimates as in the proof in Theorem 5.4 and Lemma 5.3, all other derivatives of the second correction are bounded and we get   \begin{align*} \| \mathscr{E}^{\textbf{n}}_2(u) \|_{2,\widehat{T}} \preccurlyeq \| u \|_{1/2^*,E_1}. \end{align*} Step 3. Similarly to Step 2 we correct now the values on the last edge $$E_3$$ with two more corrections using the same blending coefficients; thus we define   \begin{align*} \mathscr{E}^{\textbf{n}}(u)(x,y)&:= \mathscr{E}^{\textbf{n}}_2(u)(x,y) - b\left(\frac{y}{x+y}\right)\mathscr{E}^{\textbf{n}}_2(u) (0,x+y) \\ &\quad\ - c\left(\frac{y}{x+y}\right)(x+y) \mathscr{E}^{\textbf{n}}_{2,x}(u)(0,x+y) +c\left(\frac{y}{x+y}\right)(x+y) \mathscr{E}^{\textbf{n}}_{2,y}(u)(0,x+y). \end{align*} As in Step 2 we use the properties of the blending coefficients and similar techniques to bound the integrals to derive (5.18), (5.19), (5.20) and (5.21):   \begin{gather*} \mathscr{E}^{\textbf{n}}(u)|_{E_1} = \mathscr{E}^{\textbf{n}}(u)|_{E_2} = \mathscr{E}^{\textbf{n}}(u)|_{E_3} = 0 \\ (\nabla \mathscr{E}^{\textbf{n}}(u)\cdot n)|_{E_2} = (\nabla \mathscr{E}^{\textbf{n}}(u)\cdot n)|_{E_3} = 0 \quad \textrm{and} \quad (\nabla \mathscr{E}^{\textbf{n}}(u)\cdot n)|_{E_1} = u, \\ \| \mathscr{E}^{\textbf{n}}(u) \|_{2,\widehat{T}} \preccurlyeq \| u \|_{1/2^*,E_1}. \end{gather*} It remains to show that $$\mathscr{E}^{\textbf{n}}(u)$$ belongs to $$\mathscr{P}^{k+1}(\widehat{T})$$. The idea is similar to the tangential extension. Looking at the definition of the first step (5.22), we increase the order by multiplying with $$y$$. The crucial parts are the correction terms as the blending polynomials $$b(y/(1-x))$$ and $$c(y/(1-x))$$ for the second step, and $$b(y/(x+y))$$ and $$c(y/(x+y))$$ for the third step produce singularities of order 3 in the points $$(0,0)$$ and $$(1,0)$$. To overcome this problem, note that the given polynomial has a zero of order 2 in the vertices; thus there exists a polynomial $$v \in \mathscr{P}^{k-2}(E_1)$$, such that $$u(x) = (1-x)^2v(x)$$. Using the definitions of the polynomials $$b$$ and $$c$$, the extension of the second step $$\mathscr{E}^{\textbf{n}}_2(u)$$ reads   \begin{align*} \mathscr{E}^{\textbf{n}}_2(u)(x,y) & = y \int_0^1 a(s)u(x+sy) ~\mathrm{d} s - \frac{3 y^2(1-x)-2y^3}{(1-x)^2}\int_0^1a(s) u(x+s(1-x)) ~\mathrm{d} s \\ &\quad{} - \frac{y^3-y^2(1-x)}{(1-x)^2}\int_0^1a'(s) u(x+s(1-x))s ~\mathrm{d} s. \end{align*} As $$u(x+s(1-x)) = (1-x)^2(1-s)^2v(x+s(1-x))$$ this leads to   \begin{align*} \mathscr{E}^{\textbf{n}}_2(u)(x,y) & = y \int_0^1 a(s)u(x+sy) ~\mathrm{d} s - (3 y^2(1-x)-2y^3)\int_0^1a(s) (1-s)^2v(x+s(1-x)) ~\mathrm{d} s \\ &\quad{} - (y^3-y^2(1-x))\int_0^1a'(s) (1-s)^2v(x+s(1-x))s ~\mathrm{d} s, \end{align*} thus $$\mathscr{E}^{\textbf{n}}_2(u) \in \mathscr{P}^{k+1}(\widehat{T})$$. For the third step, we do the same by writing $$u(x) = x^2w(x)$$ with $$w \in \mathscr{P}^{k-2}(E_1)$$, finally leading to $$\mathscr{E}^{\textbf{n}}(u) \in \mathscr{P}^{k+1}(\widehat{T})$$. □ Remark 5.6 In a similar way to the last step of the proof of Theorem 5.4, it is possible to define the normal extension $$\mathscr{E}^{\textbf{n}}$$ for the other two edges $$E_2$$ and $$E_3$$ using proper transformations. We then use the subscript $$\mathscr{E}^{\textbf{n}}_{E_i}(\cdot)$$ with $$i \in \{ 1,2,3 \}$$ to symbolize which extension is used. 5.4 Splitting into compatible and incompatible polynomials Using Theorem 5.5, it would now be possible to correct the normal derivative after a first extension using Theorem 5.4. The crucial point is that the polynomial would need a zero of order 2 in the vertices. The following theorem helps us later to provide a stable splitting of the correction into two parts. Theorem 5.7 Assume a given function $$u \in \mathscr{P}^k(E_1)$$, with $$u =0$$ on $$\partial E_1$$. There holds   \begin{align} \label{ext::divtheorem:A} |u'(1)| \preccurlyeq k^2 \| u \|_{1/2^*,E_1}. \end{align} (5.24) Furthermore, there exists a function $$e \in \mathscr{P}^k(E_1)$$ with $$e'(1) = 1$$ and $$e'(0) = e(0) = e(1) = 0$$ such that   \begin{align} \label{ext::divtheorem:B} \| e\|_{1/2^*,E_1} \preccurlyeq \frac{1}{k^2} \quad \textrm{and} \quad \| e\|_{0,E_1} \preccurlyeq \frac{1}{k^3}. \end{align} (5.25) Proof. We start with the second statement. For this purpose, we present the proof on the interval $$E$$ from which (5.25) follows with a linear transformation. We want to remind the reader of the definition of Jacobi polynomials with respect to the weight function $$(1-x)^\alpha$$ (see, for example, Abramowitz, 1974; Andrews et al., 1999),   \begin{align*} p^\alpha_n(x) := \frac{1}{2^nn!(1-x)^\alpha} \frac{{\rm{d}}}{{\rm{d}} x^n}\left((1-x)^\alpha(x^2-1)^n \right)\!, \quad n \in \mathbb{N}_0, \alpha > -1, \end{align*} where in the special case of $$\alpha = 0$$ the polynomials are called Legendre polynomials. For our proof, we use integrated Jacobi polynomials with $$\alpha=1$$,   \begin{align*} \widehat{p}_n(x) &:= -\int_x^1 p_{n-1}^1(s) ~\mathrm{d} s, \quad n\ge 1 \quad \textrm{and} \quad \widehat{p}_0 (x) := 1, \end{align*} and integrated Legendre polynomials,   \begin{align*} l_{n+1}(x) := -\int_x^1 p_n^0(s) ~\mathrm{d} s, \quad n \ge 0 \quad \textrm{and} \quad l_0 := -x + 1. \end{align*} These polynomials fulfil the following properties (see Andrews et al., 1999; Beuchler & Schöberl, 2006):   \begin{gather} \widehat{p}_n (1) = 0, \quad 1 \le n \le k \quad \textrm{and} \quad \widehat{p}'_n (1) = n, \quad 0 \le n \le k, \end{gather} (5.26)  \begin{gather} (2n+1) p_n^0 = (n+1) p^1_n - n p^1_{n-1}, \quad n\ge 0 \quad \textrm{and} \quad p_m^1 = \frac{1}{m+1} \sum\limits_{n=0}^m (2n+1) p_n^0, \quad m\ge 0, \end{gather} (5.27)  \begin{gather} (2n+1) l_{n+1} = p^0_{n+1} - p^0_{n-1}, \quad n>0, \end{gather} (5.28) where we have used $$p^0_{-1} := -1$$. Furthermore, we have a weighted $$L^2$$ orthogonality for the Jacobi polynomials, and due to the definition, a weighted orthogonality in the $$H^1$$-seminorm for the integrated Jacobi polynomials,   \begin{align} \label{ext::dividingtheorem::orthogonality} \int_{-1}^1 (1-x) \widehat{p}'_n(x) \widehat{p}'_m(x) ~\mathrm{d} x = \int_{-1}^1 (1-x) p^1_{n-1}(x) p^1_{m-1}(x) ~\mathrm{d} x = \delta_{n,m} \frac{2}{n+1}, \end{align} (5.29) where $$\delta$$ is the Kronecker delta. To find a proper candidate that fulfils the bounds (5.25), we first seek for the minimum of the weighted $$H^1$$-seminorm with proper restrictions; thus   \begin{align*} \tilde{e} := \arg \min\limits_{\substack{v \in \mathscr{P}^k \\ v(1) = 0 \\v'(1) = 1}} \int_{-1}^1(1-x)v'(x)^2 ~\mathrm{d} s = \arg \min\limits_{\substack{v \in \mathscr{P}^k \\ v(1) = 0 \\v'(1) = 1}} | v|^2_{1^*,E}. \end{align*} Using integrated Jacobi polynomials as the basis for $$\mathscr{P}^k(E),$$ we use the representation of $$\tilde{e}$$ with coefficients $$c_j$$ as $$\tilde{e}(x) = \sum\limits_{j=0}^k c_j \widehat{p}_j(x)$$. To determine the coefficients, so to explicitly solve the minimization problem, we first note that due to the boundary restrictions $$\tilde{e}(1) = 0$$, it is clear that $$c_0 =0$$ and with (5.26) we get $$\tilde{e}'(1) = \sum\limits_{j=1}^k c_j j = 1$$. Using (5.29) we furthermore have $$| \tilde{e}|^2_{1^*,E} = \sum\limits_{j=1}^k c_j^2\frac{2}{j+1}$$. Now we use the technique of Lagrangian multipliers; thus we define the function   \begin{align*} L(c_1,\dots,c_k,\lambda) = \sum\limits_{j=1}^k c_j^2\frac{2}{j+1} + \lambda \left(\sum\limits_{j=1}^k c_j j -1\right) \quad \textrm{with} \quad \frac{\partial L}{ \partial c_j} \stackrel{!}{=} 0 \quad \forall\, j \in \{1,\dots,k \} \quad \textrm{and} \quad \frac{\partial L}{\partial \lambda} \stackrel{!}{=} 0. \end{align*} Solving this leads to   \begin{align} \lambda = \frac{-48}{k(k+1)(k+2)(3k+1)} \quad \textrm{and} \quad c_j= \frac{-12 j(1+j)}{k(k+1)(k+2)(3k+1)}\label{ext::divi::coeffs} \end{align} (5.30) and   \begin{align} \label{ext:divtheorem:approx} | \tilde{e} |^2_{1*,E} = \sum\limits_{j=1}^k c_j^2 \frac{2}{1+j} = \frac{24}{3k^4+10k^3+9k^2+2k} \approx \frac{1}{k^4}. \end{align} (5.31) For the $$L^2$$-norm we observe using (5.27) and (5.28) that   \begin{align*} \tilde{e}(x) &=\sum\limits_{j=1}^k -c_j\int_x^1 p_{j-1}^1(s) ~\mathrm{d} s = \sum\limits_{j=1}^k -c_j\int_x^1 \frac{1}{j} \sum_{i=0}^{j-1}(2i+1)p_{i}^0(s) ~\mathrm{d} s = \sum\limits_{j=1}^k \frac{c_j}{j} \sum_{i=0}^{j-1}\underbrace{-(2i+1)\int_x^1p_{i}^0(s)}_{= (2i+1)l_{i+1}} ~\mathrm{d} s\\ & = \sum\limits_{j=1}^k \frac{c_j}{j} \sum_{i=0}^{j-1}\left(p_{i+1}^0(x) - p_{i-1}^0(x)\right) = \sum\limits_{j=1}^k \frac{c_j}{j} \left(p_{j}^0(x) + p_{j-1}^0(x)\right), \end{align*} and so with the definition of the coefficients (5.30) and using an inverse inequality (for example, Bernardi & Maday, 1997, page 253), also   \begin{align*} \| \tilde{e} \|^2_{0,E} = \sum_{j=1}^k \frac{c_j^2}{j^2} \| p_{j}^0 + p_{j-1}^0 \|_{0,E}^2 \preccurlyeq \sum_{j=1}^k \frac{j^2}{k^8} \underbrace{ \| p_{j}^0 \|_{0,E}^2}_{\preccurlyeq \frac{2}{2j+1}} \preccurlyeq \frac{j}{k^8} \sum_{j=1}^k1 \preccurlyeq \frac{1}{k^6} \quad \text{and} \quad \| \tilde{e}\|^2_{1,E} \preccurlyeq \frac{1}{k^2}. \end{align*} Using a linear transformation $$F$$ from $$E$$ to $$E_1$$, we set $$e(x) := \tilde{e}(F^{-1}(x))\frac{x^2}{2}$$ to see $$e(0) = e(1) = e'(0) = 0$$ and $$e'(1)=1$$, and   \begin{align} \label{ext:divtheorem:twoestimates} \| e \|_{0,E_1} \preccurlyeq \frac{1}{k^3} \quad \text{and} \quad \| e \|_{1,E_1} \preccurlyeq \frac{1}{k}. \end{align} (5.32) Similarly to the proof of Lemma 5.3, we now use the real method of interpolation of spaces. As $$u=0$$ on $$\partial E_1$$ we have $$u \in H^1_0(E_1)$$; thus together with $$H^{1/2}_{00}(E_1) = [L^2(E_1),H_0^1(E_1)]$$ and the definition of the norm on an interpolated space we have with (5.32),   \begin{align*} \| e\|_{1/2^*, E_1} \preccurlyeq \sqrt{ \| e\|_{0, E_1}\| e\|_{1, E_1} } \preccurlyeq \frac{1}{k^2}, \end{align*} so (5.25) is proved. For the proof of (5.24), we first define an extension from the edge to the triangle by   \begin{align*} \psi(u)(x,y):= \int_0^1 a(s) u(x+sy) ~\mathrm{d} s \quad \textrm{with} \quad a(s) = 4-6s, \end{align*} and so $$u'(1) = \frac{\partial \psi}{\partial x}(1,0)$$. Using Lemma 5.3 we get $$\| \psi \|_{1,\widehat{T}} \preccurlyeq \| u \|_{1/2,E_1}$$. Next we define the mean value along the line $$l_x := \{(x,y): 0 \le x \le 1, y \in [0,1-x] \}$$,   \begin{align*} \overline{u}(x,y) := \frac{1}{1-x} \int_x^1 \psi(x,s) ~\mathrm{d} s = \int_0^1 \psi(x,(1-x)s) ~\mathrm{d} s. \end{align*} Due to $$\frac{\partial \overline{u}}{\partial y} = 0$$, it follows with $$\frac{\partial \overline{u}}{\partial x} := \overline{u}'$$ that   \begin{align*} | \overline{u} |_{1,\widehat{T}}^2 = \int_0^1 \int_0^{1-x} \overline{u}'(x)^2 ~\mathrm{d} y ~\mathrm{d} x = \int_0^1 (1-x) \overline{u}'(x)^2 ~\mathrm{d} x = | \overline{u} |_{1*,E_1}^2 \end{align*} and so   \begin{align*} | \overline{u} |_{1*,E_1} = | \overline{u} |_{1,\widehat{T}} \preccurlyeq \| \psi \|_{1,\widehat{T}} \preccurlyeq \| u \|_{1/2, E_1}. \end{align*} Using (5.31) and $$\tilde{e}'(1)=1$$, we furthermore show that $$| \overline{u}'(1) | \preccurlyeq k^2 | \overline{u} |_{1/2^*,E_1} \preccurlyeq k^2 \| u \|_{1/2, E_1}$$, and as   \begin{align*} \overline{u}'(1) = u'(1) \overbrace{ \int_0^1a(s) ~\mathrm{d} s}^{=1} - \frac{1}{2} u'(1) \overbrace{\int_0^1a(s)s ~\mathrm{d} s}^{=0} = u'(1), \end{align*} we finally get   \begin{align*} | u'(1) | \preccurlyeq k^2 \| u \|_{1/2, E_1} \preccurlyeq k^2 \| u \|_{1/2^*, E_1}.\quad\square \end{align*} 5.5 Proof of Theorem 5.1 Proof. In the first step, we use Theorem 5.4 to find a function $$\mathscr{E}^{\tau}(\textbf{u})$$ with a proper tangential derivative. For the difference $$\textbf{u}_c:= \textbf{u} - \nabla \mathscr{E}^{\tau}(\textbf{u})$$ it follows that on the boundary $$\partial \widehat{T}$$ we have $$\textbf{u}_c \cdot \tau = 0$$. Now let $$\textbf{n}_i$$ be the normal vector on the edge $$E_i$$ and $$u_{\textbf{n}_i}:= \textbf{u}_c \cdot \textbf{n}_i$$, so the remaining difference of the normal component. The idea is now to split this error in two parts to use Theorems 5.5 and 5.7. We start with the lower edge $$E_1$$ and define $$u_{1}:= \textbf{u}_{c} \cdot ((x,y) - V_2) \in \mathscr{P}^{k+1}(\widehat{T})$$, where $$V_2 = (0,1)$$ is the vertex opposite to $$E_1$$. As $$((x,y) - V_2) \approx \tau$$ on the edges $$E_2$$ and $$E_3$$ we have $$u_{1}|_{E_2} = u_{1}|_{E_3} = 0$$. On the lower edge we have $$((x,y) - V_2) = (x,-1)$$ and as $$\textbf{u}_c \cdot \tau = 0$$, thus the first component of $$\textbf{u}_c$$ is equal to zero, we get $$u_{1}|_{E_1} = u_{\textbf{n}_1} \in \mathscr{P}^k(E_1)$$. Using Theorem 5.7, we find two functions $$e_0,e_1 \in \mathscr{P}^k(E_1)$$ with   \begin{align*} e'_1(1) =1,\quad e'_1(0) = e_1(0) = e_1(1) = 0 \quad \textrm{and} \quad e'_0(0) =1,\quad e'_0(1) = e_0(0) = e_0(1) = 0, \end{align*} where we mirrored the edge $$E$$ in Theorem 5.7 to find $$e_0$$. We are now able to split the error to define a good and a bad part on the edge $$E_1$$ by   \begin{align*} u_{\textbf{n}_1}^b := (u_{1}|_{E_1})'(1) e_1 + (u_{1}|_{E_1})'(0) e_0 \quad \textrm{and} \quad u_{\textbf{n}_1}^g := u_{\textbf{n}_1} - u_{\textbf{n}_1}^b. \end{align*} The second function $$u_{\textbf{n}_1}^g$$ is good in the sense of having a zero of order 2 in the vertices, so $$u_{\textbf{n}_1}^g \in \mathscr{P}_{00}^k(E_1)$$; thus we apply Theorem 5.5. For the other two edges, we proceed similarly (see Remark 5.6) to finally define   \begin{align*} \mathscr{E}(\textbf{u}):= \mathscr{E}^{\tau}(\textbf{u}) + \mathscr{E}^{\textbf{n}}_{E_1}(u_{\textbf{n}_1}^g) + \mathscr{E}^{\textbf{n}}_{E_2}(u_{\textbf{n}_2}^g) + \mathscr{E}^{\textbf{n}}_{E_3}(u_{\textbf{n}_3}^g). \end{align*} Note that due to (5.19) and (5.20) the normal derivatives of the different corrections do not interfere. As $$\mathscr{E}^{\textbf{n}}(u_{\textbf{n}_i}^g) = 0$$ (see (5.18)) on the boundary $$\partial \widehat{T}$$ for $$i=1,2,3$$ the corresponding tangential derivative is also zero; thus we have   \begin{align*} \nabla \mathscr{E}(\textbf{u}) \cdot \tau = \nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \tau + \nabla \mathscr{E}^{\textbf{n}}_{E_1}(u_{\textbf{n}_1}^g)\cdot \tau + \nabla \mathscr{E}^{\textbf{n}}_{E_2}(u_{\textbf{n}_2}^g)\cdot \tau + \nabla \mathscr{E}^{\textbf{n}}_{E_3}(u_{\textbf{n}_3}^g)\cdot \tau = \nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \tau = \textbf{u} \cdot \tau, \end{align*} so property (5.1) is proved. For $$\mathscr{E}^{\textbf{n}}_{E_1}(u_{\textbf{n}_1}^g),$$ we get using (5.21), (5.25) and (5.24) as $$u_{1}|_{E_1}=0$$ on $$\partial E_1$$,   \begin{align*} \| \mathscr{E}^{\textbf{n}}_{E_1}(u_{\textbf{n}_1}^g) \|_{2,\widehat{T}} &\preccurlyeq \| u_{\textbf{n}_1}^g \|_{1/2^*,E_1} = \| u_{\textbf{n}_1} - u_{\textbf{n}_1}^b \|_{1/2^*,E_1} \\ &\preccurlyeq \| u_{\textbf{n}_1} \|_{1/2^*,E_1} + |(u_{1}|_{E_1})'(1)| \| e_1\|_{1/2^*,E_1} + |(u_{1}|_{E_1})'(0)| \| e_0\|_{1/2^*,E_1} \\ &\preccurlyeq \| u_{\textbf{n}_1} \|_{1/2^*,E_1} + \|u_{1}\|_{1/2^*,E_1} \frac{k^2}{k^2} + \|u_{1}\|_{1/2^*,E_1} \frac{k^2}{k^2} \preccurlyeq \|u_{1}\|_{1/2^*,E_1}. \end{align*} As $$u_{1}|_{E_2} = u_{1}|_{E_3} = 0$$ we bound $$\|u_{1}\|_{1/2^*,E_1}$$ by the $$H^1$$-norm on the triangle; thus we get the estimate $$\| \mathscr{E}^{\textbf{n}}_{E_1}(u_{\textbf{n}_1}^g) \|_{2,\widehat{T}} \preccurlyeq \| u_{1} \|_{1,\widehat{T}} \preccurlyeq \| \textbf{u}_{c} \|_{1,\widehat{T}}$$. Using the same arguments for the other two normal extensions and inequality (5.5) implies   \begin{align*} \| \mathscr{E}(\textbf{u}) \|_{2,\widehat{T}} \le \| \mathscr{E}^{\tau}(\textbf{u}) \|_{2,\widehat{T}} + 3 \| \textbf{u}_{c} \|_{1,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}} + \| \textbf{u} - \nabla \mathscr{E}^{\tau}(\textbf{u}) \|_{1,\widehat{T}} \preccurlyeq \| \textbf{u} \|_{1,\widehat{T}}; \end{align*} thus property (5.3) is proved. To show (5.2) first note that on the boundary $$\partial \widehat{T}$$ we have   \begin{align*} \nabla \mathscr{E}(\textbf{u}) \cdot \textbf{n} &= \nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \textbf{n} + \sum\limits_{i=1}^3 \nabla \mathscr{E}^{\textbf{n}}_{E_i}(u_{\textbf{n}_i}^g) \cdot n =\nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \textbf{n} + \overbrace{\sum\limits_{i=1}^3 u_{\textbf{n}_i}}^{=\textbf{u}_c \cdot n} - \sum\limits_{i=1}^3 u_{\textbf{n}_i}^b\\ & = \nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \textbf{n} + \textbf{u} \cdot \textbf{n} - \nabla \mathscr{E}^{\tau}(\textbf{u}) \cdot \textbf{n} - \sum\limits_{i=1}^3 u_{\textbf{n}_i}^b= \textbf{u} \cdot \textbf{n} - \sum\limits_{i=1}^3 u_{\textbf{n}_i}^b, \end{align*} and as $$u_{\textbf{n}_i}^b|_{E_j} = 0$$ for $$j \neq i$$ it follows that $$\| \left(\textbf{u} - \nabla \mathscr{E}(\textbf{u}) \right) \cdot \textbf{n} \|_{0,\partial \widehat{T}} \le \sum\limits_{i=1}^3 \| u_{\textbf{n}_i}^b \|_{0,E_i}$$. As before, we use (5.25) and (5.24) to get   \begin{align*} \| u_{\textbf{n}_1}^b \|_{0,E_1} &\preccurlyeq |(u_{1}|_{E_1})'(1)| \| e_1 \|_{0,E_1} + |(u_{1}|_{E_1})'(0)| \| e_0 \|_{0,E_1} \preccurlyeq \frac{1}{k} \|u_1\|_{1/2^*,E_1}\\ &\preccurlyeq \frac{1}{k} \|u_1\|_{1,\widehat{T}} \preccurlyeq \frac{1}{k} \|\textbf{u}_c\|_{1,\widehat{T}} = \frac{1}{k} \|\textbf{u} - \nabla \mathscr{E}^{\tau}(\textbf{u}) \|_{1,\widehat{T}} \preccurlyeq \frac{1}{k}\|\textbf{u}\|_{1,\widehat{T}}, \end{align*} and with a similar estimate for $$u_{\textbf{n}_2}^b$$ and $$u_{\textbf{n}_3}^b$$ we finally get (5.2):   \begin{align*} \| \left(\textbf{u} - \nabla \mathscr{E}(\textbf{u}) \right) \cdot \textbf{n} \|_{0,\partial \widehat{T}} \preccurlyeq \frac{1}{k}\|\textbf{u}\|_{1,\widehat{T}}.\quad\square \end{align*} References Abramowitz, M. & Stegun, I. A. ( 1974) Handbook of Mathematical Functions, With Formulas, Graphs, and Mathematical Tables . Dover Publications. Ainsworth, M. & Coggins, P. ( 2002) A uniformly stable family of mixed $$hp$$-finite elements with continuous pressures for incompressible flow. IMA J. Numer. Anal. , 22, 307– 327. Google Scholar CrossRef Search ADS   Ainsworth, M. & Demkowicz, L. ( 2009) Explicit polynomial preserving trace liftings on a triangle. Math. Nachr. , 282, 640– 658. Google Scholar CrossRef Search ADS   Andrews, G., Askey, R. & Roy, R. ( 1999) Special Functions . Encyclopedia of Mathematics and its Applications. Cambridge: Cambridge University Press. Google Scholar CrossRef Search ADS   Arnold, D. N., Brezzi, F., Cockburn, B. & Marini, L. D. ( 2001/02) Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J. Numer. Anal. , 39, 1749– 1779. Google Scholar CrossRef Search ADS   Babuška, I. & Suri, M. ( 1987) The $$h-p$$ version of the finite element method with quasiuniform meshes. ESAIM Math. Model. Numer. Anal. , 21, 199– 238. Google Scholar CrossRef Search ADS   Baker, G. A., Jureidini, W. N. & Karakashian, O. A. ( 1990) Piecewise solenoidal vector fields and the Stokes problem. SIAM J. Numer. Anal. , 27, 1466– 1485. Google Scholar CrossRef Search ADS   Belgacem, F. B. ( 1994) Polynomial extensions of compatible polynomial traces in three dimensions. Comput. Methods Appl. Mech. Engrg. , 116, 235– 241. Google Scholar CrossRef Search ADS   Bergh, J. & Löfström, J. ( 1976) Interpolation Spaces: An Introduction . Grundlehren der mathematischen Wissenschaften. Berlin-New York: Springer. Google Scholar CrossRef Search ADS   Bernardi, C. & Maday, Y. ( 1990) Relèvement polynomial de traces et applications. RAIRO Modél. Math. Anal. Numér. , 24, 557– 611. Google Scholar CrossRef Search ADS   Bernardi, C. & Maday, Y. ( 1992) Approximations spectrales de problèmes aux limites elliptiques.  Mathématiques et Applications. Berlin: Springer. Bernardi, C. & Maday, Y. ( 1997) Spectral methods. Handbook of Numerical Analysis  ( Ciarlet P. G. & Lions J. L. eds), vol. 5. Amsterdam: North-Holland, pp. 209– 485. Google Scholar CrossRef Search ADS   Bernardi, C. & Maday, Y. ( 1999) Uniform inf-sup conditions for the spectral discretization of the Stokes problem. Math. Models Methods Appl. Sci. , 9, 395– 414. Google Scholar CrossRef Search ADS   Beuchler, S. & Schöberl, J. ( 2006) New shape functions for triangular $$p$$-FEM using integrated Jacobi polynomials. Numer. Math. , 103, 339– 366. Google Scholar CrossRef Search ADS   Boffi, D., Fortin, M. & Brezzi, F. ( 2013) Mixed Finite Element Methods and Applications . Computational Mathematics. Berlin: Springer. Google Scholar CrossRef Search ADS   Brennecke, C., Linke, A., Merdon, C. & Schöberl, J. ( 2015) Optimal and pressure-independent $$L^2$$ velocity error estimates for a modified Crouzeix–Raviart Stokes element with BDM reconstructions. J. Comput. Math. , 33, 191– 208. Google Scholar CrossRef Search ADS   Brenner, S. C., Gu, S., Gudi, T. & Sung, L.-y. ( 2012) A quadratic $$C^{\circ}$$ interior penalty method for linear fourth order boundary value problems with boundary conditions of the Cahn–Hilliard type. SIAM J. Numer. Anal. , 50, 2088– 2110. Google Scholar CrossRef Search ADS   Brenner, S. C., Gudi, T. & Sung, L.-y. ( 2010) An a posteriori error estimator for a quadratic $$C^0$$-interior penalty method for the biharmonic problem. IMA J. Numer. Anal. , 30, 777– 798. Google Scholar CrossRef Search ADS   Brezzi, F. & Falk, R. S. ( 1991) Stability of higher-order Hood–Taylor method. SIAM J. Numer. Anal. , 28, 581– 590. Google Scholar CrossRef Search ADS   Chernov, A. ( 2012) Optimal convergence estimates for the trace of the polynomial $$L^2$$-projection operator on a simplex. Math. Comp. , 81, 765– 787. Google Scholar CrossRef Search ADS   Cockburn, B., Gopalakrishnan, J., Nguyen, N. C., Peraire, J. & Sayas, F.-J. ( 2011) Analysis of HDG methods for Stokes flow. Math. Comp. , 80, 723– 760. Google Scholar CrossRef Search ADS   Cockburn, B., Kanschat, G. & Schötzau, D. ( 2004) The local discontinuous Galerkin method for the Oseen equations. Math. Comp. , 73, 569– 593. Google Scholar CrossRef Search ADS   Cockburn, B., Kanschat, G. & Schotzau, D. ( 2005) A locally conservative LDG method for the incompressible Navier–Stokes equations. Math. Comp. , 74, 1067– 1095. Google Scholar CrossRef Search ADS   Cockburn, B., Kanschat, G. & Schötzau, D. ( 2007) A note on discontinuous Galerkin divergence-free solutions of the Navier–Stokes equations. J. Sci. Comput. , 31, 61– 73. Google Scholar CrossRef Search ADS   Cockburn, B., Kanschat, G., Schötzau, D. & Schwab, C. ( 2002) Local discontinuous Galerkin methods for the Stokes system. SIAM J. Numer. Anal. , 40, 319– 343. Google Scholar CrossRef Search ADS   Cockburn, B., Nguyen, N. C. & Peraire, J. ( 2010) A comparison of HDG methods for Stokes flow. J. Sci. Comput. , 45, 215– 237. Google Scholar CrossRef Search ADS   Costabel, M. & McIntosh, A. ( 2010) On Bogovskiĭ and regularized Poincaré integral operators for de Rham complexes on Lipschitz domains. Math. Z. , 265, 297– 320. Google Scholar CrossRef Search ADS   Demkowicz, L. & Babuška, I. ( 2003) $$p$$ interpolation error estimates for edge finite elements of variable order in two dimensions. SIAM J. Numer. Anal. , 41, 1195– 1208. Google Scholar CrossRef Search ADS   Demkowicz, L., Gopalakrishnan, J. & Schöberl, J. ( 2008) Polynomial extension operators. I. SIAM J. Numer. Anal. , 46, 3006– 3031. Google Scholar CrossRef Search ADS   Demkowicz, L., Gopalakrishnan, J. & Schöberl, J. ( 2009) Polynomial extension operators. II. SIAM J. Numer. Anal. , 47, 3293– 3324. Google Scholar CrossRef Search ADS   Demkowicz, L., Gopalakrishnan, J. & Schöberl, J. ( 2012) Polynomial extension operators. Part III. Math. Comp. , 81, 1289– 1326. Google Scholar CrossRef Search ADS   Donea, J. & Huerta, D. ( 2003) Finite Element Methods for Flow Problems . Hoboken, NJ: John Wiley. Google Scholar CrossRef Search ADS   Egger, H. & Schöberl, J. ( 2010) A hybrid mixed discontinuous Galerkin finite-element method for convection-diffusion problems. IMA J. Numer. Anal. , 30, 1206– 1234. Google Scholar CrossRef Search ADS   Egger, H. & Waluga, C. ( 2013) $$hp$$ analysis of a hybrid DG method for Stokes flow. IMA J. Numer. Anal. , 33, 687– 721. Google Scholar CrossRef Search ADS   Elman, H. C., Silvester, D. J. & Wathen, A. J. ( 2005) Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics . Numerical Mathematics and Scientific Computation. Oxford, NY: Oxford University Press. Fu, G., Jin, Y. & Qiu, W. ( 2016) Parameter-free superconvergent $$H(\mathrm{div})$$-conforming HDG methods for the Brinkman equations. ArXiv e-prints . Georgoulis, E. H., Hall, E. & Melenk, J. M. ( 2010) On the suboptimality of the $$p$$-version interior penalty discontinuous Galerkin method. J. Sci. Comput. , 42, 54– 67. Google Scholar CrossRef Search ADS   Girault, V. & Raviart, P.-A. ( 1986) Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms . Computational Mathematics. Berlin: Springer. Extended version. Google Scholar CrossRef Search ADS   Girault, V., Rivière, B. & Wheeler, M. F. ( 2005) A discontinuous Galerkin method with nonoverlapping domain decomposition for the Stokes and Navier–Stokes problems. Math. Comp. , 74, 53– 84. Google Scholar CrossRef Search ADS   Glowinski, R. ( 2003) Finite Element Methods for Incompressible Viscous Flow . Handbook of Numerical Analysis. Amsterdam: Elsevier, pp. 3– 1176. Grisvard, P. ( 1985) Elliptic Problems in Nonsmooth Domains . Classics in Applied Mathematics. Boston, MA: Society for Industrial and Applied Mathematics. Houston, P., Schwab, C. & Süli, E. ( 2002) Discontinuous $$hp$$-finite element methods for advection-diffusion-reaction problems. SIAM J. Numer. Anal. , 39, 2133– 2163. Google Scholar CrossRef Search ADS   Karakashian, O. & Katsaounis, T. ( 2000) A discontinuous Galerkin method for the incompressible Navier–Stokes equations. Discontinuous Galerkin Methods (Newport, RI, 1999) . Lecture Notes in Computer Science and Engineering, vol. 11. Berlin: Springer, pp. 157– 166. Google Scholar CrossRef Search ADS   Karniadakis, G. E. & Sherwin, S. J. ( 2005) Spectral/hp Element Methods for Computational Fluid Dynamics . Numerical Mathematics and Scientific Computation. Oxford: Oxford University Press. Google Scholar CrossRef Search ADS   Könnö, J. & Stenberg, R. ( 2011) H (div)-conforming finite elements for the Brinkman problem. Math. Models Methods Appl. Sci. , 21, 2227– 2248. Google Scholar CrossRef Search ADS   Lederer, P. ( 2016) Pressure-robust discretizations for Navier–Stokes equations: divergence-free reconstruction for Taylor–Hood elements and high order hybrid discontinuous Galerkin methods. Master’s Thesis , Vienna Technical University. Lehrenfeld, C. ( 2010) Hybrid discontinuous Galerkin methods for solving incompressible flow problems. Master’s Thesis, Rheinisch Westfalischen Technischen Hochschule Aachen . Lehrenfeld, C. & Schöberl, J. ( 2016) High order exactly divergence-free hybrid discontinuous Galerkin methods for unsteady incompressible flows. Comput. Methods Appl. Mech. Engrg. , 307, 339– 361. Google Scholar CrossRef Search ADS   Linke, A. ( 2014) On the role of the Helmholtz decomposition in mixed methods for incompressible flows and a new variational crime. Comput. Methods Appl. Mech. Engrg. , 268, 782– 800. Google Scholar CrossRef Search ADS   Linke, A., Matthies, G. & Tobiska, L. ( 2016) Robust arbitrary order mixed finite element methods for the incompressible Stokes equations with pressure independent velocity errors. ESAIM Math. Model. Numer. Anal. , 50, 289– 309. Google Scholar CrossRef Search ADS   Linke, A. & Merdon, C. ( 2016) On velocity errors due to irrotational forces in the Navier–Stokes momentum balance. J. Comput. Phys. , 313, 654– 661. Google Scholar CrossRef Search ADS   Maday, Y. ( 1989) Relèvements de traces polynomiales et interpolations hilbertiennes entre espaces de polynômes. C. R. Acad. Sci. Paris Sér. I Math. , 309, 463– 468. Muñoz-Sola, R. ( 1997) Polynomial liftings on a tetrahedron and applications to the $$h$$-$$p$$ version of the finite element method in three dimensions. SIAM J. Numer. Anal. , 34, 282– 314. Google Scholar CrossRef Search ADS   Nguyen, N. C., Peraire, J. & Cockburn, B. ( 2010) A hybridizable discontinuous Galerkin method for Stokes flow. Comput. Methods Appl. Mech. Engrg. , 199, 582– 597. Google Scholar CrossRef Search ADS   Nguyen, N. C., Peraire, J. & Cockburn, B. ( 2011) An implicit high-order hybridizable discontinuous Galerkin method for the incompressible Navier–Stokes equations. J. Comput. Phys. , 230, 1147– 1170. Google Scholar CrossRef Search ADS   Peetre, J. ( 1963) Nouvelles propriétés d’espaces d’interpolation. C. R. Acad. Sci. Paris , 256, 1424– 1426. Rivière, B. ( 2008) Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation.  Frontiers in Applied Mathematics, vol. 35. Philadelphia: Society for Industrial and Applied Mathematics (SIAM), pp. xxii+190. Google Scholar CrossRef Search ADS   Schöberl, J. ( 1997) NETGEN An advancing front 2D/3D-mesh generator based on abstract rules. Comput. Vis. Sci. , 1, 41– 52. Google Scholar CrossRef Search ADS   Schöberl, J. ( 2014) C++11 implementation of finite elements in NGSolve. Technical Report.  Institute for Analysis and Scientific Computing, Vienna University of Technology. Schöberl, J., Melenk, J. M., Pechstein, C. & Zaglmayr, S. ( 2008) Additive Schwarz preconditioning for $$p$$-version triangular and tetrahedral finite elements. IMA J. Numer. Anal. , 28, 1– 24. Google Scholar CrossRef Search ADS   Schötzau, D., Schwab, C. & Toselli, A. ( 2002) Mixed hp-DGFEM for incompressible flows. SIAM J. Numer. Anal. , 40, 2171– 2194. Google Scholar CrossRef Search ADS   Schwab, C. C. ( 1998) p- and hp- Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics . Numerical Mathematics and Scientific Computation. Oxford, NY: Clarendon Press. Stamm, B. & Wihler, T. P. ( 2010) $$hp$$-optimal discontinuous Galerkin methods for linear elliptic problems. Math. Comp. , 79, 2117– 2133. Google Scholar CrossRef Search ADS   Stenberg, R. & Suri, M. ( 1996) Mixed hp finite element methods for problems in elasticity and Stokes flow. Num. Math. , 72, 367– 389. Google Scholar CrossRef Search ADS   Su, Y., Chen, L., Li, X. & Xu, C. ( 2016) On the inf-sup constant of a triangular spectral method for the Stokes equations. Comput. Methods Appl. Math. , 16, 507– 522. Google Scholar CrossRef Search ADS   Toselli, A. ( 2002) $$hp$$ discontinuous Galerkin approximations for the Stokes problem. Math. Models Methods Appl. Sci. , 12, 1565– 1597. Google Scholar CrossRef Search ADS   Zhang, S. ( 2009) A family of $$Q_{k+1,k}\times Q_{k,k+1}$$ divergence-free finite elements on rectangular grids. SIAM J. Numer. Anal. , 47, 2090– 2107. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Aug 29, 2017

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations