Virtual Element Methods on Meshes with Small Edges or Faces
We consider a model Poisson problem in $\R^d$ ($d=2,3$) and establish error estimates for virtual element methods on polygonal or polyhedral meshes that can contain small edges ($d=2$) or small faces ($d=3$).
\section{Introduction}\label{sec:Introduction}
 Let ${\Omega}\subset \mathbb{R}^d$ ($d=2,3$) be a  bounded polygonal/polyhedral domain
 and $f\in L_2({\Omega})$. The Poisson problem
 with the homogeneous Dirichlet boundary condition is to find
 $u\in H^1_0({\Omega})$ such that
\begin{equation}\label{eq:Poisson}
 a(u,v)=(f,v) \qquad\forall\,v\in H^1_0({\Omega}),
\end{equation}
 where
\begin{equation}\label{eq:aDef}
  a(u,v)=\int_{\Omega} \nabla u\cdot\nabla v\,dx
\end{equation}
 and $(\cdot,\cdot)$ is the  inner product for $L_2({\Omega})$.
 Here and throughout the paper we follow standard notation for differential
 operators, function
 spaces and norms that can be found for example in
 \cite{Ciarlet:1978:FEM,ADAMS:2003:Sobolev,BScott:2008:FEM}.
\par
 Problem \eqref{eq:Poisson} can be solved by virtual element methods
 \cite{BBCMMR:2013:VEM,AABMR:2013:Projector} on
 polygonal or polyhedral meshes.
 It has been observed in numerical experiments that the convergence rates
 for the virtual element methods do not deteriorate noticeably even
 in the presence of small edges or faces
 (cf. \cite{AABMR:2013:Projector,BLR:2016:VEM,BDR:2017:VEM}).
 Our goal is to establish error estimates that justify these numerical results
 for the virtual element methods introduced in \cite{AABMR:2013:Projector}.
\par
 We will develop error estimates that are based on general shape regularity
 assumptions on the subdomains in
 the polygonal or polyhedral meshes.  For the two dimensional problem, we assume
 that (i) each polygonal subdomain  is
 star-shaped with respect to a disc whose diameter is comparable to the
 diameter of the subdomain and (ii)
 the number of edges of the subdomains is uniformly bounded.  For the three
 dimensional problem, we assume
 that (i) each polyhedral subdomain is star-shaped with respect to a ball
  whose diameter is comparable to the diameter of the subdomain;
  (ii) the number of faces of the subdomains is uniformly bounded; and
 (iii) the faces of the subdomains satisfy the two dimensional shape
 regularity assumptions.
 Our error estimates  are optimal up to at most a logarithmic
 factor that involves the ratio of the lengths of the longest edge and
 the shortest edge of each subdomain
  in two dimensions and a similar ratio over the edges of the faces on the
  subdomains in three dimensions.
\par
 The rest of the paper is organized as follows.  We  begin with a
 star-shaped condition
 in Section~\ref{sec:SS}.  Then we treat the two dimensional case in
 Section~\ref{sec:LocalVEM2D} and Section~\ref{sec:Poisson2D},
 where the analysis benefits from the
 techniques developed in \cite{BLR:2016:VEM} and \cite{BGS:2017:VEM2}.
 The extension to the three dimensional Poisson problem is presented in
 Section~\ref{sec:Poisson3D}.
 We end with some concluding remarks in Section~\ref{sec:Conclusions}.
\par
 In order to avoid the proliferation of constants, we will often use the notation
 $A\lesssim B$ to represent the statement that $A\leq (\text{constant}) B$,
 where the positive constant is independent
 of mesh sizes.  The notation $A\approx B$ is equivalent to $A\lesssim B$
 and $B\lesssim A$.
 The precise dependence of the hidden constants will be declared in the text.
\par
 To minimize the technicalities, we also assume that ${\Omega}$ is convex so that the
 solution of \eqref{eq:Poisson} belongs to $H^2({\Omega})$ by
 elliptic regularity \cite{Grisvard:1985:EPN,Dauge:1988:EBV}.
\section{A Star-Shaped Condition}\label{sec:SS}
 Let $D$ be a bounded open polygon ($d=2$) or a bounded open
 polyhedron $(d=3)$, and
 $\hspace{1pt}h_{{\scriptscriptstyle D}}$ be the diameter of $D$.
\par
 We assume that
\begin{equation}\label{eq:SA}
  \text{$D$ is star-shaped with respect to a disc/ball $\fB_{\!\ssD}\subset D$ with radius $\rho_{\!{_\ssD}}\hspace{1pt}h_{{\scriptscriptstyle D}}$}.
\end{equation}
\par
 We will denote by $\tilde\fB_{\!\ssD}$ the disc/ball concentric with $\fB_{\!\ssD}$ whose radius is $\hspace{1pt}h_{{\scriptscriptstyle D}}$.
 It is clear that
\begin{equation}\label{eq:discs}
  \fB_{\!\ssD}\subset D\subset \tilde\fB_{\!\ssD}.
\end{equation}
 \par
 Below are some consequences of the star-shaped condition \eqref{eq:SA}.
 The hidden constants in Section~\ref{subsec:Sobolev}--Section~\ref{subsec:PF}
 only depend on $\rho_{\!{_\ssD}}$, while those in Section~\ref{subsec:DE} also depend on $k$.
\subsection{Sobolev Inequalities}\label{subsec:Sobolev}
 It follows from \eqref{eq:SA} that
\begin{alignat}{3}
    \|\zeta\|_{{L_\infty(D)}}&\lesssim  \hspace{1pt}h_{{\scriptscriptstyle D}}^{-(d/2)}\|\zeta\|_{L_2(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}^{1-(d/2)}|\zeta|_{{H^1(D)}}
      +\hspace{1pt}h_{{\scriptscriptstyle D}}^{2-(d/2)}|\zeta|_{H^2(D)}&\qquad&\forall\,
      \zeta\in H^2(D),\label{eq:Sobolev}\\
 \intertext{and in the case where $d=2$,}
      \|\zeta\|_{{L_\infty(D)}}&\lesssim  \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|\zeta\|_{L_2(D)}+|\zeta|_{{H^1(D)}}
      +\hspace{1pt}h_{{\scriptscriptstyle D}}^{1/2}|\zeta|_{H^{3/2}(D)}&\qquad&\forall\,
      \zeta\in H^{3/2}(D).\label{eq:Sobolev2}
\end{alignat}
 Details can be found in \cite[Lemma~4.3.4]{BScott:2008:FEM} and
 \cite[Section~6]{DS:1980:BH}.
\subsection{Bramble-Hilbert Estimates}\label{subsec:BH}
 Condition \eqref{eq:SA} also implies the following Bramble-Hilbert
 estimates \cite{BH:1970:Lemma}:
\begin{alignat}{3}
  \inf_{q\in\mathbb{P}_\ell}|\zeta-q|_{H^m(D)}&\lesssim
   h_D^{\ell+1-m}|\zeta|_{H^{\ell+1}(D)}
  &\qquad&\forall\,\zeta\in H^{\ell+1}(D),
  \,\ell=0,\ldots,k,\;\text{and}\; m\leq \ell,
  \label{eq:BHEstimates}\\
   \inf_{q\in\mathbb{P}_\ell}|\zeta-q|_{H^m(D)}&
   \lesssim h_D^{\ell+\frac12-m}|\zeta|_{H^{\ell+\frac12}(D)}
  &\qquad&\forall\,\zeta\in H^{\ell+\frac12}(D),
   \,\ell=0,\ldots,k,\;\text{and}\; m\leq \ell.
  \label{eq:BHEstimates2}
\end{alignat}
 Details can be found in \cite[Lemma~4.3.8]{BScott:2008:FEM}  and
 \cite[Section~6]{DS:1980:BH}.
\subsection{A Lipschitz Isomorphism between $D$ and $\fB_{\!\ssD}$}\label{subsec:Lipschitz}
\par
 In view of the star-shaped condition \eqref{eq:SA},
 there exists a Lipschitz isomorphism $\Phi:\fB_{\!\ssD}\longrightarrow D$ such that
 both $|\Phi|_{W^{1,\infty}(\fB_{\!\ssD})}$ and $|\Phi^{-1}|_{W^{1,\infty}(D)}$
 are bounded by a constant that
 only depends on $\rho_{\!{_\ssD}}$ (cf.  \cite[Section~1.1.8]{Mazya:2011:Sobolev}).
\par
 It follows that
 %
 \begin{align}
  |D|&\approx \hspace{1pt}h_{{\scriptscriptstyle D}}^d,\label{eq:G1}\\
  |\partial D|&\approx \hspace{1pt}h_{{\scriptscriptstyle D}}^{d-1},\label{eq:G2}
\end{align}
 where $|D|$ is the area of $D$ ($d=2$) or the volume of $D$ ($d=3$), and $|\partial D|$ is the
 arclength  of $\partial D$ ($d=2$) or the surface area of $D$ ($d=3$).
\par
 Moreover we have  (cf. \cite[Theorem~4.1]{Wloka:1987:PDE})
\begin{alignat}{3}
  \|\zeta\|_{L_2(\partial D)}&\approx \|\zeta\circ\Phi\|_{L_2(\partial\fB_{\!\ssD})}
  &\qquad&\forall\,\zeta\in L_2(\partial D),\label{eq:L2Iso1}\\
  \|\zeta\|_{L_2(D)}&\approx \|\zeta\circ\Phi\|_{L_2(\fB_{\!\ssD})}
  &\qquad&\forall\,\zeta\in L_2(D),\label{eq:L2Iso2}\\
  |\zeta|_{H^1(\partial D)}&\approx|\zeta\circ\Phi|_{H^1(\partial\fB_{\!\ssD})}
  &\qquad&\forall\,\zeta\in H^1(\partial D),\label{eq:H1Iso1}\\
  |\zeta|_{{H^1(D)}}&\approx|\zeta\circ\Phi|_{H^1(\fB_{\!\ssD})}
  &\qquad&\forall\,\zeta\in {H^1(D)},\label{eq:H1Iso2}\\
  |\zeta|_{H^{1/2}(\partial D)}&\approx |\zeta\circ\Phi|_{H^{1/2}(\partial\fB_{\!\ssD})}
  &\qquad&\forall\,\zeta\in H^{1/2}(\partial D).
  \label{eq:HalfIso}
 \end{alignat}

\subsection{Poincar\'e-Friedrichs Inequalities}\label{subsec:PF}
 The Bramble-Hilbert estimate \eqref{eq:BHEstimates} and the geometric estimates
 \eqref{eq:G1}--\eqref{eq:G2}  imply the following Poincar\'e-Friedrichs inequalities:
\begin{alignat}{3}
 h_D^{-(d/2)} \|\zeta\|_{L_2(D)}&\lesssim h_D^{-d}\Big|\int_{D} \zeta\,dx\Big|
 +\hspace{1pt}h_{{\scriptscriptstyle D}}^{1-(d/2)}|\zeta|_{{H^1(D)}}
    &\qquad&\forall\,\zeta\in {H^1(D)}, \label{eq:PF1}\\
    h_D^{-(d/2)} \|\zeta\|_{L_2(D)}&\lesssim h_D^{-(d-1)}\Big|\int_{\partial D} \zeta\,ds\Big|
    +\hspace{1pt}h_{{\scriptscriptstyle D}}^{1-(d/2)}|\zeta|_{{H^1(D)}}
    &\qquad&\forall\,\zeta\in {H^1(D)}.\label{eq:PF2}
\end{alignat}

\subsection{Estimates for $|\cdot|_{H^{1/2}(\partial D)}$}\label{subsec:Half}
 On the circle $\partial \fB_{\!\ssD}$, we have a standard estimate
\begin{equation*}
   |\zeta|_{H^{1/2}(\partial\fB_{\!\ssD})}\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1/2}\|\zeta\|_{L_2(\partial\fB_{\!\ssD})}
   +\hspace{1pt}h_{{\scriptscriptstyle D}}^{1/2}|\zeta|_{H^1(\partial \fB_{\!\ssD})}
   \qquad\forall\,\zeta\in H^1(\partial \fB_{\!\ssD}).
\end{equation*}
  It then follows from a Poincar\'e-Friedrichs inequality
  for a circle that
\begin{equation*}
  |\zeta|_{H^{1/2}(\partial\fB_{\!\ssD})}=|\zeta-\bar\zeta|_{H^{1/2}(\partial \fB_{\!\ssD})}
        \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1/2}\|\zeta-\bar\zeta\|_{L_2(\partial\fB_{\!\ssD})}
        +\hspace{1pt}h_{{\scriptscriptstyle D}}^{1/2}|\zeta|_{H^1(\partial \fB_{\!\ssD})}
        \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{1/2}|\zeta|_{H^1(\partial\fB_{\!\ssD})},
\end{equation*}
  where $\bar\zeta$ is the mean of $\zeta$ over $\partial \fB_{\!\ssD}$.
  Therefore, in view of \eqref{eq:H1Iso1} and \eqref{eq:HalfIso}, we have
 %
 \begin{equation}\label{eq:HalfAndOne}
    |\zeta|_{H^{1/2}(\partial D)}\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{1/2}|\zeta|_{H^1(\partial D)}
    \qquad\forall\,\zeta\in H^1(\partial D).
\end{equation}
\par
 Similarly, it follows from \eqref{eq:H1Iso2}, \eqref{eq:HalfIso}
 and the trace theorem for $H^1(\fB_{\!\ssD})$ that
\begin{alignat}{3}
  |\zeta|_{H^{1/2}(\partial D)}&\lesssim |\zeta|_{{H^1(D)}}
  &\qquad&\forall\,\zeta\in {H^1(D)}.\label{eq:HalfTrace}
\end{alignat}
\subsection{Trace Inequalities}\label{subsec:Trace}
  It follows from \eqref{eq:L2Iso1}, \eqref{eq:L2Iso2}, \eqref{eq:H1Iso2}
  and  standard (scaled) trace inequalities
  for $H^1(\fB_{\!\ssD})$    that
\begin{alignat}{3}
  \|\zeta\|_{L_2(\partial D)}^2&\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|\zeta\|_{L_2(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{{H^1(D)}}^2
  &\qquad&\forall\,\zeta\in {H^1(D)}.\label{eq:Trace}
\end{alignat}
\par
 We also have trace inequalities for the $H^1$ norm on $\partial D$ that require a different derivation.
\begin{lemma}\label{lem:CZ1}
  Let $e$ be an edge of $D\subset \mathbb{R}^2$.  We have
\begin{equation*
  \hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{H^1(e)}^2\lesssim |\zeta|_{H^1(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{H^{3/2}(D)}^2 \qquad
  \forall\,\zeta\in H^{3/2}(D).
\end{equation*}
\end{lemma}
\begin{proof} By scaling we can assume $\hspace{1pt}h_{{\scriptscriptstyle D}}=1$. Without
loss of generality we may also assume that
\begin{equation}\label{eq:CZ1Est0}
   \int_D \zeta\,dx=0.
\end{equation}
\par
  The existence of the Lipschitz isomorphism $\Phi:D\longrightarrow \fB_{\!\ssD}$ implies that the domain
 $D$ satisfies
 a uniform cone condition \cite[Section~4.8]{ADAMS:2003:Sobolev},
 with one reference cone and a finite cover of $\bar D$ that contains a
 fixed number of congruent open discs.  Furthermore, the angle and the height of the
 reference cone and the radius of the open discs only depend on $\rho_{\!{_\ssD}}$.  It follows that there
 exists a Calderon-Zygmund extension operator
 $E:H^1(D)\longrightarrow H^1(\mathbb{R}^2)$
 (cf. \cite[Theorem~5.28]{ADAMS:2003:Sobolev}) such that
 $E$ maps $H^2(D)$ into $H^2(\mathbb{R}^2)$ and the restriction of
 $E\zeta$ to $D$ equals $\zeta$.  Moreover we have
\begin{equation*}
  \|E\zeta\|_{H^1(\mathbb{R}^2)}\lesssim \|\zeta\|_{H^1(D)} \quad\forall\,\zeta\in{H^1(D)}
  \quad\text{and}\quad
  \|E\zeta\|_{H^2(\mathbb{R}^2)}\lesssim \|\zeta\|_{H^2(D)} \quad\forall\,\zeta\in{H^2(D)}.
\end{equation*}
 It follows from the interpolation of Sobolev spaces \cite[Chapter~7]{ADAMS:2003:Sobolev}
 that
\begin{align}\label{eq:CZ1Est1}
  \|E\zeta\|_{H^{3/2}(\mathbb{R}^2)}\lesssim
  \|\zeta\|_{H^{3/2}(D)}
  \lesssim  |\zeta|_{H^1(D)}
  +|\zeta|_{H^{3/2}(D)}
  \qquad\forall\,\zeta\in H^{3/2}(D),
\end{align}
 where we have used \eqref{eq:PF2} and \eqref{eq:CZ1Est0}.
\par
 Let $e$ be an edge of $D$, $\tilde e$ be the infinite line that contains $e$
  and $G$ be a half-plane that borders $\tilde e$.
 Then we have, by \eqref{eq:CZ1Est1} and
 the trace theorem \cite[Theorem~8.1]{Wloka:1987:PDE},
\begin{align*
 |\zeta|_{H^1(e)}&=|E\zeta|_{H^1(e)}
          \leq |E\zeta|_{H^1(\tilde e)}
          \lesssim \|E\zeta\|_{H^{3/2}(G)}
          \lesssim  |\zeta|_{H^1(D)}+|\zeta|_{H^{3/2}(D)}.
\end{align*}
\end{proof}
\par
 The proof of the following result is similar.
\begin{lemma}\label{lem:CZ2}
  Let $F$ be a face of $D\subset \mathbb{R}^3$.  We have
\begin{equation*
  \hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{H^1(F)}^2\lesssim |\zeta|_{H^1(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}^{2}|\zeta|_{H^2(D)}^2 \qquad
  \forall\,\zeta\in H^{2}(D).
\end{equation*}

\end{lemma}

\subsection{A Lifting Operator}\label{subsec:Lifting}
 It follows from \eqref{eq:L2Iso2}, \eqref{eq:H1Iso2}, \eqref{eq:HalfIso} and the
 inverse trace theorem for $H^1(\fB_{\!\ssD})$
 (cf. \cite[Theorem~8.8]{Wloka:1987:PDE}) that
 there exists a linear operator $\mathrm {Tr}^\dag:H^{1/2}(\partial D)\longrightarrow {H^1(D)}$
 such that
\begin{equation*
 \mathrm {Tr}^\dag \zeta=\zeta\;\text{on}\;\partial D

 \quad\text{and}\quad
 \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|\mathrm {Tr}^\dag \zeta\|_{L_2(D)}+|\mathrm {Tr}^\dag \zeta|_{{H^1(D)}}\leq C|\zeta|_{H^{1/2}(\partial D)},
\end{equation*}
 where the constant $C$ depends only on $\rho_{\!{_\ssD}}$.
 %
\subsection{Some Estimates for Polynomials}\label{subsec:DE}
 Let $\mathbb{P}_k$ be the space of polynomials of total degree $\leq k$ in $d$ variables.
 We obtain the following estimates by using the equivalence of norms
 on finite dimensional vector spaces and scaling.
\begin{lemma}\label{lem:DiscreteEstimates}
 We have
\begin{alignat}{3}
   \|p\|_{L_2(\partial D)}^2&\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|p\|_{L_2(D)}^2&\qquad&\forall\,p\in\mathbb{P}_k,
   \label{eq:DiscreteEstimate0}\\
    |p|_{{H^1(D)}}&\lesssim  h_D^{-1}\|p\|_{L_2(D)}
      &\qquad&\forall\,p\in \mathbb{P}_k,
      \label{eq:DiscreteEstimate1}\\
      \|p\|_{L_\infty(D)}&\lesssim |\bar{p}_{_{\partial D}}|+\hspace{1pt}h_{{\scriptscriptstyle D}}^{1-(d/2)}
      |p|_{{H^1(D)}}&\qquad&\forall\,p\in\mathbb{P}_k,
       \label{eq:DiscreteEstimate2}\\
          \|p\|_{L_\infty(D)}&\lesssim |\bar{p}_{_{D}}|+\hspace{1pt}h_{{\scriptscriptstyle D}}^{1-(d/2)}
          |p|_{{H^1(D)}}&\qquad&\forall\,p\in\mathbb{P}_k,
       \label{eq:DiscreteEstimate3}
\end{alignat}
 where $\bar p_{_{\partial D}}$ is the mean of $p$ over $\partial D$ and
 $\bar p_{_D}$ is the mean of $p$ over $D$.
\end{lemma}
\begin{proof}
 In view of \eqref{eq:discs}, we have
\begin{align}\label{eq:DS2}
 \|p\|_{L_\infty(D)}\leq \|p\|_{L_\infty(\tilde\mathfrak{B}_{{\scriptscriptstyle D}})}
 &\lesssim \|p\|_{L_\infty(\mathfrak{B}_{{\scriptscriptstyle D}})}\\
 &\lesssim (\text{diam}\,\mathfrak{B}_{{\scriptscriptstyle D}})^{-d/2}\|p\|_{L_2(\mathfrak{B}_{{\scriptscriptstyle D}})}
 \leq (\text{diam}\,\mathfrak{B}_{{\scriptscriptstyle D}})^{-d/2}\|p\|_{L_2(D)}.\notag
\end{align}
 The estimate \eqref{eq:DiscreteEstimate0}
 then follows from \eqref{eq:G2} and \eqref{eq:DS2}:
\begin{align*}
  \|p\|_{L_2(\partial D)}^2
                     \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{d-1}\|p\|_{L_\infty(D)}^2
                     \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|p\|_{L_2(D)}^2.
\end{align*}
\par
 Similarly, we have
\begin{equation*}
  |p|_{{H^1(D)}}\leq |p|_{H^1(\tilde\mathfrak{B}_{{\scriptscriptstyle D}})}\lesssim |p|_{H^1(\mathfrak{B}_{{\scriptscriptstyle D}})}
  \lesssim (\text{diam}\,\mathfrak{B}_{{\scriptscriptstyle D}})^{-1}\|p\|_{L_2(\mathfrak{B}_{{\scriptscriptstyle D}})}\leq
  (\text{diam}\,\mathfrak{B}_{{\scriptscriptstyle D}})^{-1}\|p\|_{L_2(D)},
\end{equation*}
 which together with \eqref{eq:SA} implies \eqref{eq:DiscreteEstimate1}.
\par
 The estimates \eqref{eq:DiscreteEstimate2} and \eqref{eq:DiscreteEstimate3}
 follow immediately from \eqref{eq:G1}--\eqref{eq:G2}, \eqref{eq:PF1}--\eqref{eq:PF2}  and
 \eqref{eq:DS2}.
\end{proof}
 %
\begin{lemma}\label{lem:Polynomial}
  Given any $p\in \mathbb{P}_{k-2}$ $(k\geq2)$, there exists $q\in\mathbb{P}_k$ such that
  $\Delta q=p$ and
\begin{equation*}
    \|\nabla q\|_{L_2(B)}\leq C (\mathrm{diam}\, B)\|p\|_{L_2(B)}
\end{equation*}
 where $B\subset \mathbb{R}^d$ is any ball and the positive constant $C$ depends only on $k$.
\end{lemma}
\begin{proof}  By scaling
  it suffices to treat the case where $B$ is a unit ball.
  Since $\Delta$ maps $\mathbb{P}_k$ onto $\mathbb{P}_{k-2}$, there exists an operator
 $\Delta^\dag:\mathbb{P}_{k-2}\longrightarrow \mathbb{P}_k$ such that $\Delta \Delta^\dag$ is the
 identity operator on $\mathbb{P}_{k-2}$, and we can take $q=\Delta^\dag p$.
   The lemma follows from the observation that
 both $\|p\|_{L_2(B)}$ and $\|\Delta^\dag p\|_{L_2(B)}$ are norms on
 $\mathbb{P}_{k-2}$ together with a standard inverse estimate
 \cite{Ciarlet:1978:FEM,BScott:2008:FEM}.
\end{proof}
\par

\section{Local Virtual Element Spaces in Two Dimensions}\label{sec:LocalVEM2D}
 In this section we obtain properties of the local virtual element spaces
 that will be used in the stability
 and error analyses in Section~\ref{sec:Poisson2D}.
\par
   Let the space $\mathbb{P}_k(D)$ be the restriction of $\mathbb{P}_k$
 to $D$ and the space of $\mathbb{P}_k(\p D)$ be defined by
\begin{equation*}
  \mathbb{P}_k(\p D)=\{v\in C(\partial D):\,v\big|_e\in \mathbb{P}_k(e) \;\text{for all}\;e\in\cE_{\ssD}\},
\end{equation*}
 where $C(\partial D)$ is the space of continuous functions on $\partial D$,
 $\mathbb{P}_k(e)$
 is the restriction of $\mathbb{P}_k$ to the edges $e$, and $\cE_{\ssD}$ is the set
 of the edges of $D$.  The length of an edge $e$ is denoted by $h_e$.
\subsection{The Projection $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}$ and the Space $\cQ^k(D)$}\label{subsec:PODCQD}
 The Sobolev space ${H^1(D)}$ is a Hilbert space under the inner product
\begin{equation}\label{eq:InnerProduct}
 (\!(\zeta,\eta)\!)=(\nabla \zeta,\nabla \eta)
     +\Big(\int_{\partial D}\zeta\,ds\Big)\Big(\int_{\partial D} \eta\,ds\Big).
\end{equation}
 The projection operator  from ${H^1(D)}$ onto $\mathbb{P}_k(D)$ with respect
 to $(\!(\cdot,\cdot)\!)$ is denoted by $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}$, i.e.,
 $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} \zeta\in \mathbb{P}_k(D)$ satisfies
\begin{equation}\label{eq:PODDef}
  (\!(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta,q)\!)=(\!(\zeta,q)\!) \qquad\forall\,q\in\mathbb{P}_k(D).
\end{equation}
 In particular we have
\begin{equation}\label{eq:Projection}
 \Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} p=p \qquad\forall\,p\in\mathbb{P}_k(D).
\end{equation}
\par
 It is straight-forward to check that \eqref{eq:PODDef} is equivalent to
\begin{alignat}{3}
   \int_D \nabla (\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta)\cdot\nabla q\,dx&=\int_D
  \nabla \zeta\cdot \nabla q\,dx\label{eq:POD1}\\
   &=\int_{\partial D} \zeta(n\cdot\nabla q)\,ds-\int_D \zeta(\Delta q)\,dx
  &\qquad&\forall\,q\in \mathbb{P}_k(D),\notag
\end{alignat}
 together with
\begin{alignat}{3}
   \int_{\partial D} \Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta\,ds&=\int_{\partial D}\zeta\,ds.
   \label{eq:POD2}
\end{alignat}
\par
 For $k\geq 1$,
 the virtual element space $\cQ^k(D)\subset {H^1(D)}$ is defined by the
 following conditions:  $v\in {H^1(D)}$ belongs to $\cQ^k(D)$ if and only if
 (i) the trace of $v$ on $\partial D$ belongs to $\mathbb{P}_k(\p D)$,
  (ii) the distribution $-\Delta v$ belongs to $\mathbb{P}_k(D)$,
 and (iii) we have
\begin{equation}\label{eq:Condition3}
 \Pi_{k,D}^0\hspace{1pt} v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v\in \mathbb{P}_{k-2}(D),
\end{equation}
 where $\Pi_{k,D}^0\hspace{1pt}$ is the projection from $L_2(D)$ onto $\mathbb{P}_k(D)$ and
 $\mathbb{P}_{-1}(D)=\{0\}$.
\begin{remark}\label{rem:Continuity}{\rm
 It follows from elliptic regularity for bounded Lipschitz domains
 \cite[Section~1.2]{Kenig:1994:CBMS}
 and conditions (i) and (ii) in the definition of $\cQ^k(D)$ that
 $\cQ^k(D)\subset C(\bar D)$.}
\end{remark}
\begin{remark}\label{rem:dof}{\rm The dimension of $\cQ^k(D)$ is
 the sum of the dimension of $\mathbb{P}_k(\p D)$ and
 the dimension of $\mathbb{P}_{k-2}(D)$ (cf. \cite{AABMR:2013:Projector}).
 The degrees of freedom consist of (i) the values of $v$ at
 the vertices of $D$ and
 at the points in the interior of the edges  that  together determine
 $\mathbb{P}_k(\p D)$,  and (ii) the moments of
 $\Pi_{k-2,D}^0 v$.  The set of the nodes in (i) will be denoted by $\mathcal{N}_{\p D}$.}
\end{remark}
\begin{remark}\label{rem:Computable}
{\rm
 It follows from \eqref{eq:POD1} and \eqref{eq:POD2}
  that the polynomial $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v$ can be computed in terms of the
  degrees of freedom of
  $v\in\cQ^k(D)$.  Moreover, the polynomial
  $\Pi_{k,D}^0\hspace{1pt} v$ can also be computed through \eqref{eq:Condition3}.}
\end{remark}
\subsection{A Minimum Energy Principle}\label{subsec:MEP}
 The following minimum energy principle is useful for bounding the $H^1$
 norm of a virtual element function.
\begin{lemma}\label{lem:MEP}
   The inequality
\begin{equation*
  |v|_{{H^1(D)}}\leq |\zeta|_{{H^1(D)}}
\end{equation*}
 holds for any $v\in\cQ^k(D)$ and $\zeta\in {H^1(D)}$ such that $\zeta-v=0$
 on $\partial D$ and $\Pi_{k,D}^0\hspace{1pt}(\zeta-v)=0$.
\end{lemma}
\begin{proof}  It follows from condition (ii) in the definition of $\cQ^k(D)$ that
\begin{equation*}
   \int_D \nabla v\cdot\nabla (\zeta-v)\,dx=\int_D (-\Delta v) (\zeta-v)dx=0
\end{equation*}
 and hence
 $
   |\zeta|_{H^1(D)}^2=|\zeta-v|_{H^1(D)}^2+|v|_{H^1(D)}^2.
 $
\end{proof}
\subsection{A Maximum Principle}\label{subsec:MP}
 We begin with a result from \cite[Lemma~3.3]{BLR:2016:VEM}.
\begin{lemma}\label{lem:BLR}
   There exists a positive constant $C$, depending only on
   $\rho_{\!{_\ssD}}$ and $k$, such that
\begin{equation}\label{eq:BLR}
 \|\Delta v\|_{L_2(D)}\leq C\hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|\nabla v\|_{L_2(D)} \qquad\forall\,v\in\cQ^k(D).
\end{equation}
\end{lemma}
\begin{proof} By scaling we may assume $\hspace{1pt}h_{{\scriptscriptstyle D}}=1$.
\par
  Let $\phi\geq0$ be a smooth (bump) function supported on the disc
  $\fB_{\!\ssD}$ with radius $\rho_{\!{_\ssD}}$
  such that
\begin{equation}\label{eq:BumpFunction}
  \int_D \phi\,dx=1, \quad |\phi|\lesssim 1 \quad \text{and}
  \quad |\nabla\phi(x)|\lesssim 1.
\end{equation}
 We have,
 by  the equivalence of norms on finite dimensional vector spaces, scaling
 and \eqref{eq:discs},
\begin{equation}\label{eq:Fundamental1}
  \|p\|_{L_2(D)}^2\leq  \|p\|_{L_2(\tilde\fB_{\!\ssD})}^2\lesssim \|p\|_{L_2(\fB_{\!\ssD})}^2
  \lesssim \int_{\fB_{\!\ssD}} p^2\phi\,dx
  \qquad\forall\,p\in\mathbb{P}_k.
\end{equation}
\par
 Since $\Delta v\in\mathbb{P}_k$,  it follows from \eqref{eq:DiscreteEstimate1}
 (with $\hspace{1pt}h_{{\scriptscriptstyle D}}=1$), \eqref{eq:Fundamental1}
 and integration by parts that
\begin{align*}
  \|\Delta v\|_{L_2(D)}^2&\lesssim \int_{\fB_{\!\ssD}} (\Delta v)^2\phi\,dx\\
    &= -\int_{\fB_{\!\ssD}} \nabla v\cdot(\phi\nabla(\Delta v)+(\Delta v)\nabla\phi\big)dx\\
          &\lesssim \|\nabla v\|_{L_2(D)} \big(\|\nabla(\Delta v)\|_{L_2(D)}+\|\Delta v\|_{L_2(D)}\big)
           \lesssim  \|\nabla v\|_{L_2(D)}\big(\|\Delta v\|_{L_2(D)}\big),
\end{align*}
 which implies \eqref{eq:BLR} (with $\hspace{1pt}h_{{\scriptscriptstyle D}}=1$).
\end{proof}
  The following maximum principle will be used in the analysis of the
  interpolation operator in Section~\ref{subsec:Interpolation}
  (cf. Lemma~\ref{lem:1DInterpolationError}),
  and in the stability and error analyses for
  virtual element methods in three dimensions (cf. \eqref{eq:3DIDtbar}
  and Lemma~\ref{lem:3DSDBdd}).
\begin{lemma}\label{lem:MaximumPrinciple}
  There exists a positive constant $C$, depending only on $\rho_{\!{_\ssD}}$ and $k$, such that
\begin{equation*}
  \|v\|_{{L_\infty(D)}}\leq C\big[ \|v\|_{L_\infty(\partial D)}+|v|_{H^1(D)}\big]
  \qquad\forall\,v\in \cQ^k(D).
\end{equation*}
\end{lemma}
\begin{proof}
 There exists $q\in \mathbb{P}_{k+2}$ such that
\begin{equation}\label{eq:MP2}
 \Delta q=\Delta v \quad\text{and}\quad \|\nabla q\|_{L_2(\fB_{\!\ssD})}
 \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}\|\Delta v\|_{L_2(\fB_{\!\ssD})}
\end{equation}
 by Lemma~\ref{lem:Polynomial} (with $p=\Delta v\in\mathbb{P}_k$).
\par
 Without loss of generality we may
 assume that the mean of $q$ over $D$ is zero.  Therefore we have
\begin{equation*
 \|q\|_{{L_\infty(D)}}\lesssim |q|_{{H^1(D)}}
       \lesssim |q|_{H^1(\tilde\fB_{\!\ssD})}
       \lesssim |q|_{H^1(\fB_{\!\ssD})}
       \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}\|\Delta v\|_{L_2(\fB_{\!\ssD})}
     
       \lesssim |v|_{{H^1(D)}}\notag
\end{equation*}
 by \eqref{eq:discs}, \eqref{eq:DiscreteEstimate3}, Lemma~\ref{lem:BLR},
 \eqref{eq:MP2} and scaling.
\par
 It then follows from
 the maximum principle for the harmonic function $v-q$
 (cf. \cite{Evans:2010:PDE}) that
\begin{align*}
 \|v\|_{{L_\infty(D)}}&\leq \|v-q\|_{{L_\infty(D)}}+\|q\|_{{L_\infty(D)}}\\
   &\leq \|v-q\|_{L_\infty(\partial D)}+\|q\|_{{L_\infty(D)}}
       \lesssim \|v\|_{L_\infty(\partial D)}+|v|_{H^1(D)}.
\end{align*}
\end{proof}
\subsection{The Semi-norm ${|\!|\!|\cdot|\!|\!|_{k,D}}$}\label{subsec:Norm}
 The semi-norm $|\!|\!|\cdot|\!|\!|_{k,D}$ on ${H^1(D)}$ is defined by
\begin{equation}\label{eq:tbarNormDef}
  |\!|\!| \zeta|\!|\!|_{k,D}^2=\|\Pi_{k-2,D}^0\zeta\|_{L_2(D)}^2+
  \hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{e\in\cE_{\ssD}}\|\Pi_{k-1,e}^0 \zeta\|_{L_2(e)}^2,
 \end{equation}
 where
 $\Pi_{k-1,e}^0$ is the orthogonal projection from $L_2(e)$ onto $\mathbb{P}_{k-1}(e)$.
\par
 It follows from \eqref{eq:Trace} and \eqref{eq:tbarNormDef} that
\begin{equation}\label{eq:tbarBdd}
  |\!|\!|\zeta|\!|\!|_{k,D}\leq C\big(\|\zeta\|_{L_2(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{{H^1(D)}}\big)
  \qquad\forall\,
  \zeta\in {H^1(D)},
\end{equation}
 where the positive constant $C$ depends only on $\rho_{\!{_\ssD}}$ and $k$.
\par
 We also have, by \eqref{eq:G2} and a standard estimate for polynomials
 in one variable,
\begin{align}\label{eq:tbarBdd2}
 |\!|\!| v|\!|\!|_{k,D}&\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}\|v\|_{L_\infty(\partial D)}
 +\|\Pi_{k-2,D}^0 v\|_{L_2(D)}\\
    &\lesssim  \hspace{1pt}h_{{\scriptscriptstyle D}}\Big(\sum_{p\in\mathcal{N}_{\p D}} v^2(p)\Big)^\frac12
    +\|\Pi_{k-2,D}^0 v\|_{L_2(D)}
   \qquad \forall\,v\in\cQ^k(D),\notag
\end{align}
 where the hidden constant depends only on $\rho_{\!{_\ssD}}$ and $k$.
\subsection{Estimates for $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}$}\label{subsec:PODEstimates}
\par
 All the hidden constants in this subsection depend only on $\rho_{\!{_\ssD}}$ and $k$.
 Besides the obvious stability estimate
\begin{equation}\label{eq:PODStability1}
  |\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta|_{{H^1(D)}}\leq |\zeta|_{{H^1(D)}}\qquad\forall\,\zeta\in {H^1(D)}
\end{equation}
 that follows from \eqref{eq:POD1}, we also have
 a stability estimate  for $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}$  in terms of $\|\cdot\|_{L_2(D)}$ and the semi-norm
 $|\!|\!|\cdot|\!|\!|_{k,D}$.
\begin{lemma}\label{lem:PODStability2}
 We have
\begin{equation*
 \|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta\|_{L_2(D)}\lesssim |\!|\!| \zeta|\!|\!|_{k,D} \qquad\forall\,\zeta\in {H^1(D)}.
\end{equation*}
\end{lemma}
\begin{proof}
   It follows from \eqref{eq:POD1} that
\begin{align*}
   &\int_D \nabla(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} \zeta)\cdot\nabla(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta)\,dx
   =\int_{\partial D}\zeta n\cdot\nabla(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta)ds
   -\int_D \zeta \Delta (\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta)dx\\
  &\hspace{40pt}\leq \Big(\sum_{e\in\cE_{\ssD}}\|\Pi_{k-1,e}^0\zeta\|_{L_2(e)}^2\Big)^\frac12
    \Big(\sum_{e\in\cE_{\ssD}}\|\nabla\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta\|_{L_2(e)}^2\Big)^\frac12\\
     &\hspace{80pt}+\|\Pi_{k-2,D}^0\zeta\|_{L_2(D)}\|\Delta(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta)\|_{L_2(D)},
\end{align*}
 and we have, by \eqref{eq:DiscreteEstimate0} and \eqref{eq:DiscreteEstimate1},
\begin{align*}
  \sum_{e\in\cE_{\ssD}}\|\nabla\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta\|_{L_2(e)}^2\lesssim
         \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|\nabla \Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta\|_{L_2(D)}^2 \quad\text{and}\quad
  \|\Delta(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta)\|_{L_2(D)}&\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|\nabla \Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta\|_{L_2(D)}.
\end{align*}
 It follows that
\begin{equation*}
  \|\nabla\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta\|_{L_2(D)}
  \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\Big(\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{e\in\cE_{\ssD}}\|\Pi_{k-1,e}^0\zeta\|_{L_2(e)}^2
  +\|\Pi_{k-2,D}^0\zeta\|_{L_2(D)}^2\Big)^\frac12
    =\hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|\!|\!|\zeta|\!|\!|_{k,D}.
\end{equation*}
\par
 Moreover \eqref{eq:G2}, \eqref{eq:POD2} and  \eqref{eq:tbarNormDef}, imply
\begin{align*}
    \Big|\int_{\partial D} \Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} \zeta\,ds\Big|
  =\Big|\sum_{e\in\cE_{\ssD}}\int_e \Pi_{0,e}^0\zeta\,ds\Big|
  &\leq
  \sum_{e\in\cE_{\ssD}}\sqrt{h_e}\|\Pi_{k-1,e}^0\zeta\|_{L_2(e)}\\
  &\lesssim \sqrt{\hspace{1pt}h_{{\scriptscriptstyle D}}}
    \Big(\sum_{e\in\cE_{\ssD}}\|\Pi_{k-1,e}^0\zeta\|_{L_2(e)}^2\Big)^\frac12\leq
    |\!|\!|\zeta|\!|\!|_{k,D}.
\end{align*}
\par
 Finally we have, by \eqref{eq:PF2},
\begin{align*}
  \|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta\|_{L_2(D)}&\leq \Big|\int_{\partial D} \Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} \zeta\,ds\Big| +
      \hspace{1pt}h_{{\scriptscriptstyle D}} \|\nabla\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta\|_{L_2(D)} \lesssim |\!|\!| \zeta|\!|\!|_{k,D}.
\end{align*}
\end{proof}
\par
 We can now establish error estimates for $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}$.
\par
\begin{lemma}\label{lem:PODErrors}
 We have
\begin{alignat}{3}
   \|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} \zeta\|_{L_2(D)}&\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell+1}|\zeta|_{H^{\ell+1}(D)}
 &\qquad&\forall\,\zeta\in H^{\ell+1}(D),\,0\leq\ell\leq k,
 \label{eq:LTwoPOD}\\
   |\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta|_{{H^1(D)}}&\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell}|\zeta|_{H^{\ell+1}(D)}
 &\qquad&\forall\,\zeta\in H^{\ell+1}(D),\,1\leq\ell\leq k,
 \label{eq:HOnePOD}\\
 |\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} \zeta|_{H^2(D)}&\lesssim h_D^{\ell-1}|\zeta|_{H^{\ell+1}(D)}&\qquad&
  \forall\,\zeta\in H^{\ell+1}(D),\,1\leq\ell\leq k.\label{eq:HTwoPOD}
\end{alignat}
\end{lemma}
\begin{proof}
  The estimate \eqref{eq:HOnePOD} follows immediately from \eqref{eq:BHEstimates},
  \eqref{eq:Projection} and
  \eqref{eq:PODStability1}.
\par
In view of \eqref{eq:DiscreteEstimate1} and \eqref{eq:PODStability1}, we have
\begin{equation*}
   |\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta|_{H^2(D)}\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta|_{{H^1(D)}}\leq\hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|\zeta|_{{H^1(D)}},
\end{equation*}
 which together with \eqref{eq:BHEstimates} and
  \eqref{eq:Projection} implies \eqref{eq:HTwoPOD}.
\par
 Similarly we have, by \eqref{eq:tbarBdd}
  and Lemma~\ref{lem:PODStability2}
\begin{equation*}
  \|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta\|_{L_2(D)}\lesssim
  |\!|\!| \zeta|\!|\!|_{k,D}
  \lesssim \|\zeta\|_{L_2(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{{H^1(D)}},
\end{equation*}
 which together with \eqref{eq:BHEstimates} and \eqref{eq:Projection} implies
 \eqref{eq:LTwoPOD}.
\end{proof}
\subsection{Estimates for $\Pi_{k,D}^0\hspace{1pt}$}\label{subsec:PDZEstimates}
 All the hidden constants in this subsection only depend on $\rho_{\!{_\ssD}}$ and $k$.
 We have an obvious stability estimate
\begin{equation}\label{eq:PDZLTwo}
   \|\Pi_{k,D}^0\hspace{1pt} \zeta\|_{L_2(D)}\leq \|\zeta\|_{L_2(D)} \qquad\forall\,\zeta\in L_2(D)
 \end{equation}
 %
  and an obvious relation
\begin{equation}\label{eq:PDZInvariance}
 \Pi_{k,D}^0\hspace{1pt} q=q\qquad\forall\,q\in \mathbb{P}_k(D).
\end{equation}
 %
 \par
  It follows from \eqref{eq:BHEstimates}, \eqref{eq:PDZLTwo} and
  \eqref{eq:PDZInvariance} that
\begin{equation} \label{eq:PDZLTwoError}
  \|\zeta-\Pi_{k,D}^0\hspace{1pt}\zeta\|_{L_2(D)}\lesssim
   \hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell+1}|\zeta|_{H^{\ell+1}(D)} \qquad\forall\,\zeta\in H^{\ell+1}(D),\,
      0\leq\ell\leq k.
\end{equation}
\par
 We also have a  stability estimate for $\Pi_{k,D}^0\hspace{1pt}$ in $|\cdot|_{{H^1(D)}}$.
\begin{lemma}\label{lem:PDZHOne}
   We have
 %
\begin{equation}\label{eq:PDZHOne}
     |\Pi_{k,D}^0\hspace{1pt} \zeta|_{{H^1(D)}}\lesssim |\zeta|_{{H^1(D)}} \qquad\forall\,\zeta\in {H^1(D)}.
\end{equation}
\end{lemma}
\begin{proof}  This is a consequence of  \eqref{eq:DiscreteEstimate1},
 \eqref{eq:PODStability1},
 \eqref{eq:LTwoPOD} and \eqref{eq:PDZLTwoError}:
\begin{align*}
|\Pi_{k,D}^0\hspace{1pt} \zeta|_{{H^1(D)}}&\leq |\Pi_{k,D}^0\hspace{1pt} \zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta|_{{H^1(D)}}+|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta|_{{H^1(D)}}\\
    &\lesssim  \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|\Pi_{k,D}^0\hspace{1pt} \zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta\|_{L_2(D)}+|\zeta|_{{H^1(D)}}\\
 &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\big(\|\Pi_{k,D}^0\hspace{1pt} \zeta-\zeta\|_{L_2(D)}+\|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta\|_{L_2(D)}\big)
  +|\zeta|_{{H^1(D)}}  \lesssim |\zeta|_{{H^1(D)}}.
\end{align*}
\end{proof}
\par
 We can then derive error estimates for $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}$ by combining the
 Bramble-Hilbert estimates
 and Lemma~\ref{lem:PDZHOne}.
\begin{lemma}\label{lem:PDZErrors}
 We have
\begin{alignat}{3}
       |\zeta-\Pi_{k,D}^0\hspace{1pt}\zeta|_{{H^1(D)}}&\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^\ell|\zeta|_{H^{\ell+1}(D)}
  &\qquad&\forall\,\zeta\in H^{\ell+1}(D), \, 1\leq\ell\leq k, \label{eq:PDZHOneError}\\
  |\zeta-\Pi_{k,D}^0\hspace{1pt}\zeta|_{H^2(D)}&\lesssim h^{\ell-1}|\zeta|_{H^{\ell+1}(D)}
     &\qquad&\forall\,\zeta\in H^{\ell+1}(D), \, 1\leq\ell\leq k.\label{eq:PDZHTwoError}
\end{alignat}
\end{lemma}
\begin{proof} In view of \eqref{eq:PDZInvariance},
 the estimate \eqref{eq:PDZHOneError} follows from
 \eqref{eq:BHEstimates} and \eqref{eq:PDZHOne}.
\par
Similarly the estimate  \eqref{eq:PDZHTwoError} follows from \eqref{eq:BHEstimates},
 \eqref{eq:PDZInvariance} and the inequality
\begin{equation*}
 |\Pi_{k,D}^0\hspace{1pt}\zeta|_{H^2(D)}\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|\Pi_{k,D}^0\hspace{1pt}\zeta|_{H^1(D)}
 \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|\zeta|_{H^1(D)}
\end{equation*}
 obtained from \eqref{eq:DiscreteEstimate1} and \eqref{eq:PDZHOne}.
\end{proof}
 The following is another useful estimate.
\begin{lemma}\label{lem:PDZcQD}
  We have
\begin{equation*}
  \|\Pi_{k,D}^0\hspace{1pt} v\|_{L_2(D)}\lesssim C |\!|\!| v|\!|\!|_{k,D} \qquad\forall\,v\in \cQ^k(D).
\end{equation*}
\end{lemma}
\begin{proof}    Let $v\in\cQ^k(D)$ be arbitrary.  It follows from
\eqref{eq:Condition3} that
\begin{align*}
   \|\Pi_{k,D}^0\hspace{1pt} v\|_{L_2(D)}^2&=\|\Pi_{k-2,D}^0 v\|_{L_2(D)}^2+\|(\Pi_{k,D}^0\hspace{1pt}-\Pi_{k-2,D}^0)v\|_{L_2(D)}^2\\
     &=\|\Pi_{k-2,D}^0 v\|_{L_2(D)}^2+\|(\Pi_{k,D}^0\hspace{1pt}-\Pi_{k-2,D}^0)\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v\|_{L_2(D)}^2\\
      &\leq \|\Pi_{k-2,D}^0 v\|_{L_2(D)}^2+\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v\|_{L_2(D)}^2,
\end{align*}
 which together with \eqref{eq:tbarNormDef} and Lemma~\ref{lem:PODStability2}
 completes the proof.
\end{proof}
\subsection{Inverse Estimates}\label{subsec:InverseEstimate}
 These are estimates that bound the norm $|v|_{H^1(D)}$ of a virtual element function
 $v\in\cQ^k(D)$ in terms of $\|\Pi_{k-2,D}^0 v\|_{L_2(D)}$ and norms that only involve
 the boundary data of $v$.
 They are crucial for the stability analysis of virtual element methods
in Section~\ref{subsec:DiscreteProblem}.
\par
 We begin with a key lemma.
\begin{lemma}\label{lem:Fundamental}
 There exists a positive constant $C$ depending only on $\rho_{\!{_\ssD}}$ and $k$, such that
\begin{equation}\label{eq:Fundamental}
  |v|_{{H^1(D)}}\leq C\big[ \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|\!|\!| v|\!|\!|_{k,D}+|v|_{H^{1/2}(\partial D)}\big].
\end{equation}
\end{lemma}
\begin{proof}  By scaling we may assume  $\hspace{1pt}h_{{\scriptscriptstyle D}}=1$.
\par
 Let $\mathrm {Tr}^\dag$ be the lifting operator from Section~\ref{subsec:Lifting}.
 The function $w=\mathrm {Tr}^\dag v\in H^1(D)$ satisfies $w=v$ on $\partial D$ and
\begin{equation}\label{eq:w}
   \|w\|_{H^1(D)}\lesssim |v|_{H^{1/2}(D)}.
\end{equation}
\par
 Let $\phi$ be the same (bump) function in the proof of Lemma~\ref{lem:MaximumPrinciple}
  and  $\zeta=w+p\phi$, where the polynomial $p\in \mathbb{P}_k(D)$ is determined by
\begin{equation*}
   \int_D (\zeta-v)q\,dx=0 \qquad\forall\,q\in \mathbb{P}_k(D),
\end{equation*}
 or equivalently
\begin{equation}\label{eq:Fundamental2}
  \int_D  pq\phi\,dx=\int_D (v-w)q\,dx=\int_D (\Pi_{k,D}^0\hspace{1pt} v-w)q\,dx \qquad\forall\,q\in \mathbb{P}_k(D).
\end{equation}
 Then we have
\begin{equation}\label{eq:Fundamental3}
 |v|_{{H^1(D)}}\leq |\zeta|_{{H^1(D)}}
\end{equation}
 by Lemma~\ref{lem:MEP} and, in view of \eqref{eq:DiscreteEstimate1},
 \eqref{eq:BumpFunction} and \eqref{eq:w},
\begin{equation}\label{eq:Fundamental4}
 |\zeta|_{{H^1(D)}}\leq |w|_{{H^1(D)}}+|p\phi|_{{H^1(D)}}
   \lesssim |w|_{{H^1(D)}}+\|p\|_{L_2(D)}\lesssim
   |v|_{H^{1/2}(\partial D)}+\|p\|_{L_2(D)}.
\end{equation}
\par
 Note that  \eqref{eq:Fundamental1} and \eqref{eq:Fundamental2} imply
\begin{equation*}
 \|p\|_{L_2(D)}\lesssim  \|\Pi_{k,D}^0\hspace{1pt} v-w\|_{L_2(D)}
\end{equation*}
 and hence
\begin{equation}\label{eq:Fundamental5}
  \|p\|_{L_2(D)} \lesssim \|\Pi_{k,D}^0\hspace{1pt} v\|_{L_2(D)}+\|w\|_{L_2(D)}\lesssim
   |\!|\!| v|\!|\!|_{k,D}+|v|_{H^{1/2}(\partial D)}
\end{equation}
 by Lemma~\ref{lem:PDZcQD} and \eqref{eq:w}.
\par
 The estimate \eqref{eq:Fundamental} (with $\hspace{1pt}h_{{\scriptscriptstyle D}}=1$) follows from
 \eqref{eq:Fundamental3}--\eqref{eq:Fundamental5}.
\end{proof}
\begin{lemma}\label{lem:InverseEstimate1}
 There exists a positive constant $C$, depending only $\rho_{\!{_\ssD}}$ and $k$, such that
\begin{equation}\label{eq:InverseEstimate1}
  |v|_{{H^1(D)}}\leq C \big[\hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|\!|\!| v|\!|\!|_{k,D}
  +\hspace{1pt}h_{{\scriptscriptstyle D}}^{1/2}\|\partial v/\partial s\|_{L_2(\partial D)}\big]
  \qquad\forall\,v\in\cQ^k(D),
\end{equation}
 where $\partial v/\partial s$ is a tangential derivative of $v$.
\end{lemma}
\begin{proof}
 The estimate \eqref{eq:InverseEstimate1}
 follows immediately from
 \eqref{eq:HalfAndOne} and  Lemma~\ref{lem:Fundamental} .
\end{proof}
\begin{lemma}\label{lem:InverseEstimate2}
  There exists a positive constant $C$, depending only on $\rho_{\!{_\ssD}}$, $|\cE_{\ssD}|$ and $k$, such
  that
\begin{equation}\label{eq:InverseEstimate2}
 |v|_{{H^1(D)}}\leq C \big[ \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|\!|\!| v|\!|\!|_{k,D}
  +\sqrt{\ln(1+\ \tau_{\!\ssD})}\|v\|_{L_\infty(\partial D)} \big]\qquad\forall\,v\in \cQ^k(D),
\end{equation}
 where
\begin{equation}\label{eq:TauD}
 \tau_{\!\ssD}= \frac{\max_{e\in\cE_{\ssD}}h_e}{\min_{e\in\cE_{\ssD}}h_e}.
\end{equation}
\end{lemma}
\begin{proof}
 According to \cite[Lemma~5.1]{BLR:2016:VEM}, we have
\begin{equation}\label{eq:HalfAndInfty}
 |v|_{H^{1/2}(\partial D)}\leq C \sqrt{\ln(1+\tau_{\!\ssD})}\|v\|_{L_\infty(\partial D)}
 \qquad\forall\,v\in\cQ^k(D),
\end{equation}
 where the positive constant $C$ only depends  on $\rho_{\!{_\ssD}}$, $|\cE_{\ssD}|$ and $k$.
 \par
 The estimate \eqref{eq:InverseEstimate2} follows from
  Lemma~\ref{lem:Fundamental}
 and \eqref{eq:HalfAndInfty}.
\end{proof}
\par
 Combining \eqref{eq:tbarBdd2} and \eqref{eq:InverseEstimate2},
 we have the following corollary.
\begin{corollary}\label{cor:InverseEstimate3}
   There exists a positive constant $C$, depending only on
   $\rho_{\!{_\ssD}}$, $|\cE_{\ssD}|$ and $k$, such
 that
\begin{equation*}
 |v|_{H^1(D)}\leq C\big[ \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|\Pi_{k-2,D}^0 v\|_{L_2(D)}+\sqrt{\ln(1+\tD)}
\|v\|_{L_\infty(\partial D)}\big]
\qquad\forall\,v\in\cQ^k(D).
\end{equation*}
\end{corollary}
\subsection{The Interpolation Operator}\label{subsec:Interpolation}
 For $s>1$
 the interpolation operator $I_{k,D}:H^s(D)\longrightarrow\cQ^k(D)$
 is defined by the
 condition that $\zeta$ and $I_{k,D}\zeta$ share the same degrees
 of freedom, i.e.,
 $I_{k,D}\zeta$ agrees with $\zeta$ at the nodes in $\mathcal{N}_{\p D}$ and
\begin{equation}\label{eq:IDDef}
 \Pi_{k-2,D}^0I_{k,D}\zeta=\Pi_{k-2,D}^0\zeta.
\end{equation}
\par
  It is clear that
\begin{equation}\label{eq:IDInvariance}
  I_{k,D} q=q\qquad\forall\,q\in\mathbb{P}_k(D),
\end{equation}
 and by a standard estimate for polynomials in one variable,
\begin{equation}\label{eq:TrivialBdd}
  \|I_{k,D} \zeta\|_{L_\infty(\partial D)}\leq C \max_{p\in\mathcal{N}_{\p D}}|\zeta(p)|
  \leq\|\zeta\|_{L_\infty(\partial D)}\qquad
  \forall\,\zeta\in H^s(D) \;\text{and}\; s>1,
 \end{equation}
 where the positive constant $C$ only depends on $k$.
\par
  For the three dimensional Poisson problem, if the solution belongs to
 $H^{\ell+1}({\Omega})$, then its restriction to a face $F$
 of a polyhedral subdomain belongs to $H^{\ell+\frac12}(F)$.  Therefore
 below we also consider the interpolants of functions in $H^{\ell+\frac12}(D)$.
\par
 We begin with several stability estimates for the interpolation operator.
\begin{lemma}\label{lem:IDtbar}
  There exists a positive constant $C$, depending only on $\rho_{\!{_\ssD}}$ and $k$, such that
\begin{alignat}{3}
   \tbarI_{k,D}\zeta|\!|\!|_{k,D}&\leq C\big[\|\zeta\|_{L_2(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{H^1(D)}
   +\hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta|_{H^2(D)}\big]
   &\qquad&\forall\,\zeta\in H^2(D),\label{eq:IDtbar1}\\
    \tbarI_{k,D}\zeta|\!|\!|_{k,D}&\leq C\big[\|\zeta\|_{L_2(D)}
    +\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{H^1(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}^{3/2}|\zeta|_{H^{3/2}(D)}\big]
   &\qquad&\forall\,\zeta\in H^{3/2}(D).\label{eq:IDtbar2}
\end{alignat}
\end{lemma}
\begin{proof} Let $\zeta\in {H^2(D)}$ (resp., $H^{3/2}(D)$)
 be arbitrary.  From \eqref{eq:G2},
  \eqref{eq:tbarBdd2}, \eqref{eq:IDDef} and \eqref{eq:TrivialBdd}, we have
\begin{align*
  \tbarI_{k,D}\zeta|\!|\!|_{k,D}&\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}\|I_{k,D}\zeta\|_{L_\infty(\partial D)}
  +\|\Pi_{k-2,D}^0I_{k,D}\zeta\|_{L_2(D)}\\
  &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}\|\zeta\|_{L_\infty(\partial D)}+\|\Pi_{k-2,D}^0 \zeta\|_{L_2(D)}\lesssim
  \hspace{1pt}h_{{\scriptscriptstyle D}}\|\zeta\|_{L_\infty(D)},\notag
\end{align*}
 which together with \eqref{eq:Sobolev} (resp., \eqref{eq:Sobolev2})
 implies \eqref{eq:IDtbar1} (resp., \eqref{eq:IDtbar2}).
\end{proof}
\begin{lemma}\label{lem:IDFundamental1}
 We have
\begin{align}
 |I_{k,D}\zeta|_{{H^1(D)}}&\lesssim |\zeta|_{{H^1(D)}}+\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{H^2(D)}
 \label{eq:IDHOne}\\
\intertext{for all $\zeta\in{H^2(D)}$, and}
  |I_{k,D}\zeta|_{{H^1(D)}}&\lesssim
   |\zeta|_{{H^1(D)}}+\hspace{1pt}h_{{\scriptscriptstyle D}}^{1/2}|\zeta|_{H^{3/2}(D)}\label{eq:IDHOneHalf}
\end{align}
 for all $\zeta\in H^{3/2}(D)$,
  where the hidden constants only depend on $\rho_{\!{_\ssD}}$, $|\cE_{\ssD}|$ and $k$.
\end{lemma}
\begin{proof}  Let $\zeta\in H^2(D)$ be arbitrary and $\bar\zeta_D$
be the mean of $\zeta$ over $D$.
 Since $I_{k,D}\bar\zeta_D=\bar\zeta_D$,
  it follows from \eqref{eq:InverseEstimate1} that
\begin{align*}
  |I_{k,D}\zeta|_{{H^1(D)}}^2&=|I_{k,D}(\zeta-\bar\zeta_D)|_{H^1(D)}^2\\
  &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-2}\tbarI_{k,D}(\zeta-\bar\zeta_D)|\!|\!|_{k,D}^2
  +\hspace{1pt}h_{{\scriptscriptstyle D}}\|\partial (I_{k,D}\zeta)/\partial s\|_{L_2(\partial D)}^2.
\end{align*}
 We have, by a standard interpolation  estimate
 in one dimension and \eqref{eq:Trace} (applied to the first order
 derivatives of $\zeta$),
\begin{align*}
 \hspace{1pt}h_{{\scriptscriptstyle D}}\|\partial(I_{k,D} \zeta)/\partial s\|_{L_2(\partial D)}^2&\lesssim\sum_{e\in\cE_{\ssD}}
 \hspace{1pt}h_{{\scriptscriptstyle D}}\|\partial\zeta/\partial s\|_{L_2(e)}^2\\
 &\lesssim \sum_{e\in\cE_{\ssD}}\big[|\zeta|_{H^1(D)}^2+
 \hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta|_{H^2(D)}^2\big]\lesssim
 |\zeta|_{H^1(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta|_{H^2(D)}^2.
\end{align*}
\par
 These two estimates together with \eqref{eq:PF2} and \eqref{eq:IDtbar1}
 imply \eqref{eq:IDHOne}.
\par
 Similarly we obtain \eqref{eq:IDHOneHalf} by replacing \eqref{eq:Trace} with
 the estimate in Lemma~\ref{lem:CZ1}.
\end{proof}
\begin{lemma}\label{lem:IDFundamental2}
 We have
\begin{align}
 \|I_{k,D}\zeta\|_{L_2(D)}&\lesssim  \|\zeta\|_{L_2(D)}
 +\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{{H^1(D)}}+\hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta|_{H^2(D)}
 \label{eq:IDLTwo}\\
\intertext{for all $\zeta\in {H^2(D)}$, and}
 \|I_{k,D}\zeta\|_{L_2(D)}&\lesssim
  \|\zeta\|_{L_2(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{{H^1(D)}}+\hspace{1pt}h_{{\scriptscriptstyle D}}^{3/2}|\zeta|_{H^{3/2}(D)}
  \label{eq:IDLTwoHalf}
\end{align}
 for all $\zeta\in H^{3/2}(D)$,  where the hidden constants only
 depend on $\rho_{\!{_\ssD}}$, $|\cE_{\ssD}|$ and $k$.
\end{lemma}
\begin{proof}
  From Lemma~\ref{lem:PODStability2} and \eqref{eq:LTwoPOD} we have
\begin{align}\label{eq:LTwoID}
 \|I_{k,D}\zeta\|_{L_2(D)}&\lesssim \|I_{k,D}\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta\|_{L_2(D)}
 +\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta\|_{L_2(D)}\\
 &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}|I_{k,D}\zeta|_{H^1(D)}+\tbarI_{k,D}\zeta|\!|\!|_{k,D}.\notag
\end{align}
 The estimates \eqref{eq:IDLTwo} and \eqref{eq:IDLTwoHalf} follow from
 \eqref{eq:LTwoID}, Lemma~\ref{lem:IDtbar}
 and Lemma~\ref{lem:IDFundamental1}.
\end{proof}
\par
 We can now derive error estimates for the interpolation operator.
\begin{lemma}\label{lem:InterpolationError}
 We have, for $1\leq \ell\leq k$,
\begin{alignat}{3}
 \|\zeta-I_{k,D}\zeta\|_{L_2(D)}+\|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta\|_{L_2(D)}
 &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell+1}|\zeta|_{H^{\ell+1}(D)}&\qquad&
 \forall\,\zeta\in H^{\ell+1}(D),
 \label{eq:LTwoPODID}\\
 |\zeta-I_{k,D}\zeta|_{{H^1(D)}}+|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta|_{{H^1(D)}}
 &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell}|\zeta|_{H^{\ell+1}(D)}&\qquad&
 \forall\,\zeta\in H^{\ell+1}(D),
 \label{eq:HOnePODID}\\
 |\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta|_{H^2(D)}&\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell-1}|\zeta|_{H^{\ell+1}(D)}
 &\qquad&\forall\,\zeta\in H^{\ell+1}(D),
  \label{eq:HTwoPODID}
\end{alignat}
 and
\begin{alignat}{3}
 \|\zeta-I_{k,D} \zeta\|_{L_2(D)}&\lesssim
  \hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell+\frac12}|\zeta|_{H^{\ell+\frac12}(D)}
  &\qquad&\forall\,\zeta\in H^{\ell+\frac12}(D),
  \label{eq:HalfIDLTwoError}\\
   |\zeta-I_{k,D} \zeta|_{H^1(D)}&\lesssim
  \hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell-\frac12}|\zeta|_{H^{\ell+\frac12}(D)}
  &\qquad&\forall\,\zeta\in H^{\ell+\frac12}(D),
  \label{eq:HalfIDHOneError}
\end{alignat}
 where the hidden constants  only depend on $\rho_{\!{_\ssD}}$, $|\cE_{\ssD}|$ and $k$.
\end{lemma}
\begin{proof}
   In view of  \eqref{eq:tbarBdd2},
   Lemma~\ref{lem:PODStability2}, \eqref{eq:IDtbar1} and \eqref{eq:IDLTwo},
   we have,
   for any $\zeta\in H^2(D)$,
\begin{align*}
   \|I_{k,D}\zeta\|_{L_2(D)}+\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} I_{k,D}\zeta\|_{L_2(D)}
   &\lesssim \|I_{k,D}\zeta\|_{L_2(D)}+\tbarI_{k,D}\zeta|\!|\!|_{k,D}\\
     &\lesssim \|\zeta\|_{L_2(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{{H^1(D)}}+\hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta|_{H^2(D)},
\end{align*}
 which together with \eqref{eq:BHEstimates},
 \eqref{eq:Projection} and \eqref{eq:IDInvariance} implies
 \eqref{eq:LTwoPODID}.
\par
 From \eqref{eq:PODStability1} and \eqref{eq:IDHOne} we also have,
 for any $\zeta\in H^2(D)$,
\begin{align*}
 & |I_{k,D}\zeta|_{H^1({\Omega})}+|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta|_{{H^1(D)}}\leq 2|I_{k,D}\zeta|_{{H^1(D)}}
  \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|\zeta\|_{L_2(D)}+|\zeta|_{{H^1(D)}}+\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{H^2(D)},
\end{align*}
 which together with \eqref{eq:BHEstimates},
 \eqref{eq:Projection} and \eqref{eq:IDInvariance} implies
 \eqref{eq:HOnePODID}.
\par
 Similarly, by using the relation
\begin{align*}
 |\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta|_{H^2(D)}&=|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta-\Pi_{1,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta|_{H^2(D)}\\
   &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta-\Pi_{1,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta|_{H^1(D)}\\
   &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|I_{k,D}\zeta|_{H^1(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|\zeta|_{H^1(D)}
   \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-2}\|\zeta\|_{L_2(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|\zeta|_{H^1(D)}+|\zeta|_{H^2(D)}
\end{align*}
 that follows from \eqref{eq:DiscreteEstimate1},
 \eqref{eq:PODStability1} and
 \eqref{eq:IDHOne}, we can establish \eqref{eq:HTwoPODID}
 through \eqref{eq:BHEstimates}, \eqref{eq:Projection}
 and \eqref{eq:IDInvariance}.
\par
 Finally we obtain \eqref{eq:HalfIDLTwoError} and
 \eqref{eq:HalfIDHOneError}
 by replacing \eqref{eq:BHEstimates}
 (resp., \eqref{eq:IDtbar1}, \eqref{eq:IDHOne} and  \eqref{eq:IDLTwo})
 with \eqref{eq:BHEstimates2}
  (resp., \eqref{eq:IDtbar2}, \eqref{eq:IDHOneHalf} and
  \eqref{eq:IDLTwoHalf})
  in the arguments for \eqref{eq:LTwoPODID} and \eqref{eq:HOnePODID}.
\end{proof}
\par
 The proof for the following result is similar.
\begin{lemma}\label{lem:LTwoIDSpecial}
  There exists a positive constant $C$, depending only on $\rho_{\!{_\ssD}}$, $|\cE_{\ssD}|$
  and $k$, such that
 %
 \begin{equation*
   \|I_{k,D} \zeta-\Pi_{1,D}^0I_{k,D} \zeta\|_{L_2(D)}\leq C \hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta|_{H^2(D)}
   \qquad\forall\,\zeta\in H^2(D).
 \end{equation*}
\end{lemma}
\par
 We also have interpolation error estimates in the $L_\infty$ norm.
\begin{lemma}\label{lem:1DInterpolationError}
   There exists a positive constant $C$, depending only on $\rho_{\!{_\ssD}}$, $N$ and $k$,
   such that
\begin{alignat}{3}
  \|\zeta-I_{k,D}\zeta\|_{L_\infty(D)}&\leq
  C \hspace{1pt}h_{{\scriptscriptstyle D}}^\ell|\zeta|_{H^{\ell+1}(D)}&\qquad&
  \forall\,\zeta\in H^{\ell+1}(D) \;\text{and}\;
  1\leq\ell\leq k,\label{eq:1DInterpolationError}\\
    \|\zeta-I_{k,D}\zeta\|_{L_\infty(D)}&\leq C\hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell-\frac12}|\zeta|_{H^{\ell+\frac12}(D)}
    &\qquad&
  \forall\,\zeta\in H^{\ell+\frac12}(D) \;\text{and}\;1\leq\ell\leq k.
  \label{eq:1DInterpolationErrorHalf}
\end{alignat}
\end{lemma}
\begin{proof}  It follows from \eqref{eq:Sobolev}, Lemma~\ref{lem:MaximumPrinciple},
 \eqref{eq:TrivialBdd} and \eqref{eq:IDHOne} that, for any $\zeta\in {H^2(D)}$,
\begin{equation*}
  \|I_{k,D}\zeta\|_{L_\infty(D)}\lesssim \|I_{k,D}\zeta\|_{L_\infty(\partial D)}+|I_{k,D}\zeta|_{H^1(D)}
          \lesssim\hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|\zeta\|_{L_2(D)}+|\zeta|_{H^1(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{H^2(D)},
\end{equation*}
 which together with \eqref{eq:BHEstimates} and \eqref{eq:IDInvariance} imply
 \eqref{eq:1DInterpolationError}.
\par
 The proof for \eqref{eq:1DInterpolationErrorHalf} is similar, but with \eqref{eq:Sobolev}
 (resp., \eqref{eq:BHEstimates} and \eqref{eq:IDHOne}) replaced by \eqref{eq:Sobolev2}
 (resp., \eqref{eq:BHEstimates2} and \eqref{eq:IDHOneHalf}).
\end{proof}
\subsection{The Null Space of $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}$}\label{subsec:Kernel}
  Let $\mathcal{N}(\POD)=\{v\in\cQ^k(D):\,\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v=0\}$ be the null space of the projection $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}$.
  The inverse estimates in Section~\ref{subsec:InverseEstimate} can be simplified
  for functions in $\mathcal{N}(\POD)$.
\begin{lemma}\label{lem:KerPOD}
  There exists a positive constant $C$, depending only on $\rho_{\!{_\ssD}}$
  and $k$, such that
\begin{equation}\label{eq:KerPOD}
   |\!|\!| v|\!|\!|_{k,D}^2\leq C
   \hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{e\in\cE_{\ssD}}\|\Pi_{k-1,e}^0 v\|_{L_2(e)}^2
   \qquad\forall\,v\in \mathcal{N}(\POD).
\end{equation}
\end{lemma}
\begin{proof}  We can assume $k\geq2$ since
 $\Pi_{k-2,D}\hspace{1pt}v=0$ for $k=1$.
  \par
 Let $v\in\mathcal{N}(\POD)$ be arbitrary.  It follows from
 \eqref{eq:POD1} that
\begin{equation}\label{eq:KerPOD1}
  \int_D v(\Delta q)\,dx=\int_{\partial D}v(n\cdot\nabla q)\,ds
  \qquad\forall \,q\in \mathbb{P}_k(D).
\end{equation}
\par
 According to \eqref{eq:discs} and Lemma~\ref{lem:Polynomial},
 given any $p\in\mathbb{P}_{k-2}$,
  there exists $q\in\mathbb{P}_k$ such that $\Delta q=p$ and
\begin{equation}\label{eq:KerPOD12}
  \|\nabla q\|_{L_2(D)}\leq\|\nabla q\|_{L_2(\tilde\fB_{\!\ssD})}\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}\|p\|_{L_2(\tilde\fB_{\!\ssD})}
  \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}\|p\|_{L_2(\fB_{\!\ssD})}\leq \hspace{1pt}h_{{\scriptscriptstyle D}}\|p\|_{L_2(D)}.
\end{equation}
\par
 It follows from
 \eqref{eq:DiscreteEstimate0}, \eqref{eq:KerPOD1} and \eqref{eq:KerPOD12} that
\begin{align*}
    \int_D vp\,dx
 &\leq \Big(\sum_{e\in\cE_{\ssD}}\|\Pi_{k-1,e}^0 v\|_{L_2(e)}^2\Big)^\frac12
  \Big(\sum_{e\in\cE_{\ssD}}\|\nabla q\|_{L_2(e)}^2\Big)^\frac12\\
 &\lesssim \Big(\sum_{e\in\cE_{\ssD}}\|\Pi_{k-1,e}^0 v\|_{L_2(e)}^2\Big)^\frac12 \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1/2}
    \|\nabla q\|_{L_2(D)}\\
 &\lesssim \Big(\sum_{e\in\cE_{\ssD}}\|\Pi_{k-1,e}^0 v\|_{L_2(e)}^2\Big)^\frac12 \hspace{1pt}h_{{\scriptscriptstyle D}}^{1/2}
      \|p\|_{L_2(D)}
\end{align*}
  and hence
\begin{equation}\label{eq:KerPod13}
  \|\Pi_{k-2,D}^0v\|_{L_2(D)}=\max_{p\in\mathbb{P}_{k-2}}\int_D vp\,dx /\|p\|_{L_2(D)}
   \lesssim \Big(\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{e\in\cE_{\ssD}}\|\Pi_{k-1,e}^0 v\|_{L_2(e)}^2\Big)^\frac12.
\end{equation}
\par
 The estimate \eqref{eq:KerPOD} follows from \eqref{eq:tbarNormDef}
 and  \eqref{eq:KerPod13}.
\end{proof}
\begin{lemma}\label{lem:BdryEst}   For any $v\in\mathbb{P}_k(\p D)$ that vanishes at
 some point on $\partial D$,
   we have
\begin{equation*}
     \|v\|_{L_2(\partial D)} \leq C  \hspace{1pt}h_{{\scriptscriptstyle D}}
            \|\partial v/\partial s\|_{L_2(\partial D)},
\end{equation*}
 where $\partial v /\partial s$ denotes a tangential derivative of $v$ along $\partial D$ and the
 positive constant $C$ only depends on $k$.
\end{lemma}
\begin{proof}  It follows from a Poincar\'e-Friedrichs inequality on
the circle $\partial\fB_{\!\ssD}$
 that
\begin{equation*}
  \|\zeta\|_{L_2(\partial\fB_{\!\ssD})}\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{H^1(\partial\fB_{\!\ssD})}
\end{equation*}
 for any $\zeta\in H^1(\partial\fB_{\!\ssD})$ that vanishes at some point on $\partial\fB_{\!\ssD}$.
 The lemma then follows from \eqref{eq:L2Iso1} and \eqref{eq:H1Iso1}.
\end{proof}
\par
 Note that every $v\in \mathcal{N}(\POD)$ must vanish at some point on $\partial{\Omega}$ because
 $\int_{\partial D}v\,ds=0$ by \eqref{eq:POD2}.
 Therefore, in view of Lemma~\ref{lem:KerPOD} and Lemma~\ref{lem:BdryEst},
 the inverse estimates
 \eqref{eq:InverseEstimate1} and \eqref{eq:InverseEstimate2} can be simplified to
\begin{alignat*}{3}
  |v|_{{H^1(D)}}&\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{1/2}\|\partial v/\partial s\|_{L_2(\partial D)}
  &\qquad&\forall\,v\in \mathcal{N}(\POD)\\
\intertext{with a hidden constant depending only on $\rho_{\!{_\ssD}}$ and $k$, and}
 |v|_{{H^1(D)}}&\lesssim \sqrt{\ln(1+\tau_{\!\ssD})}\|v\|_{L_\infty(\partial D)}
 &\qquad&\forall\,v\in \mathcal{N}(\POD)
\end{alignat*}
 with a hidden constant that also depends on $|\cE_{\ssD}|$.
\par
 Hence we have
\begin{alignat}{3}
 |v|_{H^1(D)}^2&=|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v|_{H^1(D)}^2+|v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v|_{H^1(D)}^2\label{eq:EnergyBdd1}\\
  &\lesssim |\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v|_{H^1(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}\|\partial (v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)/\partial s\|_{L_2(\partial D)}^2
 &\qquad&\forall\,v\in\cQ^k(D),\notag\\
 \intertext{with a hidden constant depending only on $\rho_{\!{_\ssD}}$ and $k$, and also}
 |v|_{H^1(D)}^2&\lesssim |\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v|_{H^1(D)}^2+\ln(1+\tau_{\!\ssD})\|v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v\|_{L_\infty(\partial D)}^2
 &\qquad&\forall\,v\in\cQ^k(D),\label{eq:EnergyBdd2}
\end{alignat}
 with a hidden constant that also depends on $|\cE_{\ssD}|$.
\section{The Poisson Problem in Two Dimensions}\label{sec:Poisson2D}
 In this section
 we consider virtual element methods for the Poisson problem
 \eqref{eq:Poisson} in two dimensions
 and establish error estimates under two global shape regularity assumptions.
\par
 Let $\mathcal{T}_h$ be a triangulation of the convex polygon ${\Omega}\subset\mathbb{R}^2$
 by polygonal subdomains,
 where $h=\displaystyle\max_{D\in\mathcal{T}_h}\hspace{1pt}h_{{\scriptscriptstyle D}}$ is the mesh parameter.
 The global virtual finite element space $\cQ^k_h$ is given by
 $$\cQ^k_h=\{v\in H^1_0({\Omega}):\, v\big|_D\in\cQ^k(D) \quad \forall\,D\in\mathcal{T}_h\}.$$
 The space of (discontinuous) piecewise
  polynomials of degree $\leq k$ with respect to $\mathcal{T}_h$ is denoted by $\cP^k_h$.
\par
 The operators $\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} :H^1({\Omega})\longrightarrow \cP^k_h$,
 $\Pi_{k,h}^0:L_2({\Omega})\longrightarrow \cP^k_h$
 and $I_{k,h}:H^2({\Omega})\cap H^1_0({\Omega})\longrightarrow \cQ^k_h$ are defined in terms of
 their local counterparts, i.e.,
 $$(\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\zeta)\big|_D=\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(\zeta\big|_D), \quad
 (\Pi_{k,h}^0\zeta)\big|_D=\Pi_{k,D}^0\hspace{1pt}(\zeta\big|_D) \quad
   \text{and} \quad (I_{k,h}\zeta)\big|_D=I_{k,D}(\zeta\big|_D).
 $$
\par
 The piecewise $H^1$ norm with respect to $\mathcal{T}_h$ is given by
\begin{equation}\label{eq:PiecewiseHOneNOrm}
 |v|_{h,1}=\Big(\sum_{D\in\cT_h} |v|_{{H^1(D)}}^2\Big)^\frac12.
\end{equation}
\subsection{Global Shape Regularity Assumptions}
\label{subsec:GlobalShapeRegularity}
 We assume that the local shape regularity assumption \eqref{eq:SA}
 is satisfied by all
  $D\in\mathcal{T}_h$ and impose the following global regularity assumptions.
\par\noindent
{\em Assumption 1}\quad There exists a positive number $\rho\in (0,1)$,
 independent of $h$, such that
\begin{equation}\label{eq:rho}
  \rho_{\!{_\ssD}}\geq\rho \qquad\forall\,D\in\mathcal{T}_h.
\end{equation}
\par\noindent
{\em Assumption 2}\quad There exists a positive integer $N$,
 independent of $h$, such that
\begin{equation}\label{eq:N}
 |\cE_{\ssD}|\leq N\qquad\forall\,D\in\mathcal{T}_h.
\end{equation}
\par
 The hidden constants in the rest of Section~\ref{sec:Poisson2D}
 will only depend on $\rho$, $N$ and $k$.
\subsection{The Discrete Problem}\label{subsec:DiscreteProblem}
 Let the local stabilizing bilinear form $S^D_1(\cdot,\cdot)$ and
 $S^D_2(\cdot,\cdot)$ be defined by
\begin{align}
 S^D_1(w,v)&=\sum_{\mathcal{N}_{\p D}} w(p)v(p),\label{eq:SD1}\\
  S^D_2(w,v)&=\hspace{1pt}h_{{\scriptscriptstyle D}}(\partial w/\partial s,\partial v/\partial s)_{L_2(\partial D)},\label{eq:SD2}
\end{align}
 where $\partial v/\partial s$ denotes a tangential derivative of $v$ along $\partial D$.
\begin{remark}\label{rem:SD}{\rm
 The local stability bilinear form $S^D_1(\cdot,\cdot)$ is the boundary part of
 the local stability bilinear form in \cite{BBCMMR:2013:VEM}.
 The bilinear form $S^D_2(\cdot,\cdot)$ was introduced in \cite{WRR:2016:VEM}.}
\end{remark}
\begin{remark}\label{rem:AlternativeSD2}{\rm
  We can also use the bilinear form $\tilde S^D_2(\cdot,\cdot)$ defined by
\begin{equation*}
   \tilde S^D_2(w,v)=\sum_{e\in\cE_{\ssD}}h_e(\partial w/\partial s,\partial v/\partial s)_{L_2(e)}
\end{equation*}
}
\end{remark}
\begin{remark}\label{rem:NormEquivalence}
{\rm
  By the equivalence of norms on finite dimensional vector spaces, we have
 $$\sum_{p\in\mathcal{N}_e}v^2(p)\approx \sum_{e\in\cE_{\ssD}}\|v\|_{L_\infty(e)}^2
   \approx \|v\|_{L_\infty(\partial D)}^2 \qquad \forall\, v\in \mathbb{P}_k(\p D).$$
}
\end{remark}
\par
 The discrete problem for \eqref{eq:Poisson} is to
 find $u_h\in\cQ^k_h$ such that
\begin{equation}\label{eq:DiscreteProblem}
   a_h(u_h,v)=(f,\Xi_h v) \qquad \forall\,v\in\cQ^k_h,
\end{equation}
 where
\begin{align}
  a_h(w,v)&=\sum_{D\in\cT_h} \big[a^D(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} w,\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)
  +S^D(w-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} w,v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)\big],\label{eq:ah}\\
  a^D(w,v)&=\int_D \nabla w\cdot\nabla v\,dx,\label{eq:aD}
\end{align}
 $S^D(\cdot,\cdot)$ is either $S^D_1(\cdot,\cdot)$ or
 $S^D_2(\cdot,\cdot)$,
 and $\Xi_h:\cQ^k_h\longrightarrow\cP^k_h$ is given by
\begin{equation}\label{eq:Xih}
    \Xi_h=\begin{cases} \Pi_{1,h}^0&\qquad\text{if $k=1,2,$}\\[4pt]
                                 \Pi_{k-2,h}^0 &\qquad\text{if $k\geq 3$}.\\
      \end{cases}
\end{equation}
\par
 It follows from \eqref{eq:EnergyBdd1} and \eqref{eq:EnergyBdd2} that
\begin{equation}\label{eq:Stability}
  |v|_{H^1({\Omega})}^2 \lesssim
    \alpha_h a_h(v,v) \qquad\forall\,v\in \cQ^k_h,
\end{equation}
 where
\begin{equation}\label{eq:kappah}
  \alpha_h=\begin{cases}\displaystyle
          \ln(1+\max_{D\in\mathcal{T}_h}\tau_{\!\ssD})&\qquad
          \text{if $S^D(\cdot,\cdot)=S^D_1(\cdot,\cdot)$}\\
     1&\qquad\text{if $S^D(\cdot,\cdot)=S^D_2(\cdot,\cdot)$}
\end{cases}\;.
\end{equation}
 The well-posedness of the discrete problem follows from the
 stability estimate \eqref{eq:Stability}.
\par
 We will use the following properties of $\Xi_h$ in the error analysis.
\begin{lemma}\label{lem:Xih}
 We have, for $1\leq\ell\leq k$,
\begin{alignat}{3}
  (f,w-\Xi_h w)&\lesssim h^\ell|f|_{H^{\ell-1}({\Omega})} |w|_{H^1({\Omega})}
  &\qquad&\forall\,f\in H^{\ell-1}({\Omega}),\;w\in\cQ^k_h, \label{eq:Xih1}\\
    (f,I_{k,h}\zeta-\Xi_hI_{k,h}\zeta)&\lesssim
    h^{\ell+1}|f|_{H^{\ell-1}({\Omega})}|\zeta|_{H^2({\Omega})}&\qquad&
   \forall f\in H^{\ell-1}({\Omega}),\;\zeta\in H^2({\Omega}).\label{eq:Xih2}
\end{alignat}
\end{lemma}
\begin{proof}
In view of the relation
\begin{align*}
  (f,w-\Xi_h w)&=(f-\Pi_{k-2,h}^0 f,w-\Xi_h w)\\
            &\leq \|f-\Pi_{k-2,h}^0 f\|_{L_2({\Omega})}\| w-\Xi_h w\|_{L_2({\Omega})}
           \leq \|f-\Pi_{\ell-2,h}^0 f\|_{L_2({\Omega})}
           \|w-\Pi_{0,h}^0 w\|_{L_2({\Omega})},
\end{align*}
 the estimate \eqref{eq:Xih1} follows from \eqref{eq:PDZLTwoError}.
\par
  Similarly we have
\begin{align*}
   (f,I_{k,h}\zeta-\Xi_hI_{k,h}\zeta)&=(f-\Pi_{k-2,h}^0 f,I_{k,h}\zeta-\Xi_hI_{k,h}\zeta)\\
    &\leq \|f-\Pi_{k-2,h}^0 f\|_{L_2({\Omega})}\|I_{k,h}\zeta-\Xi_hI_{k,h}\zeta\|_{L_2({\Omega})}\\
    &\leq \|f-\Pi_{\ell-2,h}^0f\|_{L_2({\Omega})}
    \|I_{k,h}\zeta-\Pi_{1,h}^0I_{k,h}\zeta\|_{L_2({\Omega})},
\end{align*}
 and the estimate \eqref{eq:Xih2} follows from \eqref{eq:PDZLTwoError} and
 Lemma~\ref{lem:LTwoIDSpecial}.
  Note that this is the reason why $\Xi_h$ is chosen to be
 $\Pi_{1,h}^0$ for $k=2$ instead of $\Pi_{k-2,h}^0=\Pi_{0,h}^0$.
\end{proof}
\subsection{An Abstract Error Estimate in the Energy Norm}
\label{subsec:AbstractEnergyError}
 Let $\|\cdot\|_h=\sqrt{a_h(\cdot,\cdot)}$ be the mesh-dependent energy norm.
 Note that \eqref{eq:Stability} implies
\begin{equation}\label{eq:HOneVEM}
 |v|_{H^1({\Omega})}\lesssim
    \sqrt{\alpha_h}\|v\|_h \qquad\forall\,v\in\cQ^k_h.
\end{equation}
\par
 The discrete problem \eqref{eq:DiscreteProblem} is defined in
 terms of a non-inherited
 symmetric positive definite bilinear form.
 We have a standard error estimate (cf. \cite[Lemma~10.1.7]{BScott:2008:FEM}
 and \cite{BSS:1972:NC})
\begin{equation}\label{eq:Standard}
 \|u-u_h\|_h\leq \inf_{v\in\cQ^k_h}\|u-v\|_h+\sup_{v\in \cQ^k_h}
 \frac{a_h(u,v)-(f,\Xi_h v)}{\|v\|_h}.
\end{equation}
 The key is to control the numerator on the
 right-hand side of \eqref{eq:Standard}.
\par
 In view of \eqref{eq:Poisson}, \eqref{eq:aDef} and \eqref{eq:POD1} we
 can write, for any $v\in\cQ^k_h$,
\begin{align*}
  a_h(u,v)&=\sum_{D\in\cT_h} \big[a^D(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u,\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)+S^D(u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u,v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)\big]\\
  &=\sum_{D\in\cT_h} \big[a^D(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u,v)+S^D(u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u,v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)\big]\\
  &=\sum_{D\in\cT_h} \big[a^D(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u-u,v)+S^D(u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u,v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)\big]+(f,v)\\
   &=\sum_{D\in\cT_h} \big[a^D(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u-u,v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)+S^D(u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u,v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)\big]+(f,v),
\end{align*}
 and hence, by \eqref{eq:PODStability1}, \eqref{eq:ah} and \eqref{eq:HOneVEM},
\begin{align}\label{eq:Numerator}
 a_h(u,v)-(f,\Xi_h v)&=\sum_{D\in\cT_h} \big[a^D(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u-u,v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)
 +S^D(u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u,v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)\big]\notag\\
 &\hspace{30pt}+(f,v-\Xi_h v),\notag\\
 &\lesssim \Big(\sum_{D\in\cT_h}|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u-u|_{{H^1(D)}}|v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v|_{{H^1(D)}}\Big)
 +\|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u\|_h\|v\|_h\\
  &\hspace{30pt}+\Big(\sup_{w\in\cQ^k_h}\frac{(f,w-\Xi_h w)}{|w|_{H^1({\Omega})}}\Big)
   \sqrt{\alpha_h}\|v\|_h\notag\\
   &\lesssim \|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u\|_h\|v\|_h+
  \Big(|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u|_{h,1}+\sup_{w\in\cQ^k_h}
  \frac{(f,w-\Xi_h w)}{|w|_{H^1({\Omega})}}\Big)\sqrt{\alpha_h}\|v\|_h.\notag
\end{align}
\par
 Putting \eqref{eq:Standard} and \eqref{eq:Numerator} together
 we arrive at the
 estimate
\begin{equation}\label{eq:EnergyError}
 \|u-u_h\|_h\lesssim \|u-I_{k,h} u\|_h+\|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u\|_h
 +\sqrt{\alpha_h}\Big(|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u|_{h,1}+\sup_{w\in\cQ^k_h}
 \frac{(f,w-\Xi_h w)}{|w|_{H^1({\Omega})}}\Big).
\end{equation}
\par
 Below we will derive concrete error estimates under the assumption
 that the solution $u$ of \eqref{eq:Poisson}  belongs to $H^{\ell+1}({\Omega})$
 for $1\leq\ell\leq k$.
\subsection{Concrete Error Estimates in the Energy Norm}
\label{subsec:ConcreteEnergyErrors}
 The terms on the right-hand side of \eqref{eq:EnergyError} are
 estimated as follows.
\par\medskip\noindent{\em Estimate for the Term Involving $f$}
\par\smallskip
 Since $u\in H^{\ell+1}({\Omega})$ and $f=-\Delta u$, we have, by \eqref{eq:Xih1},
\begin{equation}\label{eq:RHSError}
 \sup_{w\in\cQ^k_h}\frac{(f,w-\Xi_h w)}{|w|_{H^1({\Omega})}}
 \lesssim h^{\ell}|u|_{H^{\ell+1}({\Omega})}.
\end{equation}
\par\medskip\noindent{\em Estimate for $|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u|_{h,1}$}
\par\smallskip
 It follows directly from \eqref{eq:HOnePOD} that
\begin{equation}\label{eq:PODError}
 |u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u|_{h,1}=\Big(\sum_{D\in\cT_h} |u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u|_{H^1(D)}^2\Big)^\frac12
 \lesssim h^\ell|u|_{H^{\ell+1}({\Omega})}.
\end{equation}
\par\medskip\noindent{\em Estimate for $\|u-I_{k,h} u\|_h$}
\par\smallskip
 We will establish the estimate
\begin{equation}\label{eq:InterpolationError}
 \|u-I_{k,h} u\|_h\lesssim h^{\ell}|u|_{H^{\ell+1}(\O)}
\end{equation}
 for both choices of $S^D(\cdot,\cdot)$.
\par
 In the case where $S^D(\cdot,\cdot)=S^D_1(\cdot,\cdot)$, it follows from
 \eqref{eq:PODStability1},
 \eqref{eq:SD1} and \eqref{eq:ah} that
\begin{align}\label{eq:InterError1}
 \|u-I_{k,h} u\|_h^2&\lesssim \sum_{D\in\cT_h} |\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)|_{H^1(D)}^2+\sum_{D\in\cT_h}
  \|(u-I_{k,D} u)-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)\|_{L_\infty(\partial D)}^2\notag\\
 &\lesssim \sum_{D\in\cT_h} \big(|u-I_{k,D} u|_{H^1(D)}^2+\|u-I_{k,D} u\|_{L_\infty(\partial D)}^2\big)\\
   &\hspace{40pt}+\sum_{D\in\cT_h} |\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)|_{L_\infty(D)}^2,\notag
\end{align}
 and we have
\begin{align}\label{eq:InterError2}
 \sum_{D\in\cT_h} |\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)|_{L_\infty(D)}^2&\lesssim
 \sum_{D\in\cT_h}\big(\,(\,\overline{\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)}\,)_{\partial D}^2
 +|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)|_{H^1(D)}^2\big)\notag\\
 &\leq
 \sum_{D\in\cT_h}\big(\,(\,\overline{u-I_{k,D} u}\,)_{\partial D}^2+|u-I_{k,D} u|_{H^1(D)}^2\big)\\
 &\lesssim \sum_{D\in\cT_h} \big(\|u-I_{k,D} u\|_{L_\infty(\partial D)}^2
 +|u-I_{k,D} u|_{H^1(D)}^2\big),\notag
\end{align}
 by \eqref{eq:DiscreteEstimate2}, \eqref{eq:Condition3} and
 \eqref{eq:PODStability1}.
\par
 The estimate \eqref{eq:InterpolationError} now follows from
 \eqref{eq:HOnePODID}, \eqref{eq:1DInterpolationError},
  \eqref{eq:InterError1} and \eqref{eq:InterError2}.
\par
 In the case where $S^D(\cdot,\cdot)=S^D_2(\cdot,\cdot)$, it follows from
 \eqref{eq:SD2} and \eqref{eq:ah} that
\begin{align*}
  \|u-I_{k,h} u\|_h^2&\lesssim \sum_{D\in\cT_h}|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)|_{H^1(D)}^2
     + \sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}} \|\partial [\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)]/\partial s\|_{L_2(\partial D)}^2\notag\\
       &\hspace{30pt}+{\sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}} \|\partial(u-I_{k,D} u)/\partial s\|_{L_2(\partial D)}^2}.
\end{align*}
 We have
\begin{equation*}
 \sum_{D\in\cT_h}|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)|_{H^1(D)}^2\leq \sum_{D\in\cT_h}|u-I_{k,D} u|_{H^1(D)}^2\lesssim
 h^{2\ell}|u|_{H^{\ell+1}(\O)}^2
\end{equation*}
 by \eqref{eq:PODStability1} and \eqref{eq:HOnePODID},
 and hence, in view of
 \eqref{eq:DiscreteEstimate0}  (applied to the first order derivatives
 of the polynomial
 $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)$),
\begin{align*}
 \sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}} \|\partial [\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)]/\partial s\|_{L_2(\partial D)}^2&\lesssim
    \sum_{D\in\cT_h} \|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)|_{H^1(D)}^2\lesssim h^{2\ell}|u|_{H^{\ell+1}(\O)}^2.
\end{align*}
 Finally it follows from
 a standard interpolation error estimate in one variable
 and \eqref{eq:HalfTrace} (applied to the $\ell$-th order
 derivatives of $u$) that
\begin{align*}
  \sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}} \|\partial(u-I_{k,D} u)/\partial s\|_{L_2(\partial D)}^2&\lesssim
  \sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{e\in\cE_{\ssD}}h_e^{2\ell-1}
        \,|\partial^\ell u/\partial s^\ell|_{H^{1/2}(e)}^2
        \lesssim h^{2\ell}|u|_{H^{\ell+1}(\O)}^2.
\end{align*}
\par
 Together these estimates imply \eqref{eq:InterpolationError}.
\goodbreak
\par\medskip\noindent{\em Estimate for $\|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u\|_h$}
\par\smallskip
 We will show that the estimate
\begin{equation}\label{eq:ProjectionError}
 \|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u\|_h\lesssim h^\ell|u|_{H^{\ell+1}(\O)}
\end{equation}
 holds for both choices of $S^D(\cdot,\cdot)$.
  \par
 In the case where $S^D(\cdot,\cdot)=S^D_1(\cdot,\cdot)$, it follows from
 \eqref{eq:Sobolev}, Lemma~\ref{lem:PODErrors}, \eqref{eq:SD1} and \eqref{eq:ah}
 that
\begin{align*}\label{eq:ProjectionEst2}
 \|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u\|_h^2
  &\lesssim \sum_{D\in\cT_h} \|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u\|_{L_\infty(\partial D)}^2\\
     &\lesssim \sum_{D\in\cT_h}\big[\hspace{1pt}h_{{\scriptscriptstyle D}}^{-2}\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u\|_{L_2(D)}^2+|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u|_{{H^1(D)}}^2
         +\hspace{1pt}h_{{\scriptscriptstyle D}}^2|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u|_{H^2(D)}^2\big]\notag\\
             &\lesssim h^{2\ell}|u|_{H^{\ell+1}({\Omega})}^2.\notag
\end{align*}
\par
 In the case where $S^D(\cdot,\cdot)=S^D_2(\cdot,\cdot)$,  we have
\begin{align*}
   \|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u\|_h^2&=\sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}} \|\partial(u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u)/\partial s\|_{L_2(D)}^2\\
 &\lesssim \sum_{D\in\cT_h} \big[|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u|_{{H^1(D)}}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}^2|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u|_{H^2(D)}^2\big]
  \lesssim h^{2\ell}|u|_{H^{\ell+1}({\Omega})}^2\notag
\end{align*}
 by \eqref{eq:Trace} (applied to the first order derivatives of
  $u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u$) and Lemma~\ref{lem:PODErrors}.
\par
  Putting \eqref{eq:EnergyError}--\eqref{eq:InterpolationError}
  and \eqref{eq:ProjectionError} together, we
 arrive at the following result.
\begin{theorem}\label{thm:ConcreteEnergyError}
  Assuming the solution $u$ of \eqref{eq:Poisson} belongs to
 $H^{\ell+1}({\Omega})$ for $\ell$ between $1$ and $k$, we have
\begin{equation*}
  \|u-u_h\|_h\leq C\sqrt{\alpha_h} h^\ell|u|_{H^{\ell+1}({\Omega})},
\end{equation*}
 where $\alpha_h$ is defined in \eqref{eq:kappah} and the positive
 constant $C$ only depends
on $\rho$, $k$ and $N$.
\end{theorem}
\par We have similar estimates for the computable approximate
 solutions $\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h$ and $\Pi_{k,h}^0 u_h$.
\begin{theorem}\label{thm:ComputableEnergyError}
  Assuming the solution $u$ of \eqref{eq:Poisson} belongs
  to $H^{\ell+1}({\Omega})$
   for some $\ell$ between $1$ and $k$, we have
\begin{equation}\label{eq:ComputableEnergyError}
  |u-u_h|_{H^1({\Omega})}+\sqrt{\alpha_h}\big[|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h|_{h,1}
  +|u-\Pi_{k,h}^0 u|_{h,1}\big]
 \leq C \alpha_h h^\ell |u|_{H^{\ell+1}({\Omega})},
\end{equation}
 where $\alpha_h$ is defined in \eqref{eq:kappah} and the positive
 constant $C$ only depends on $\rho$, $N$ and $k$.
\end{theorem}
\begin{proof} In view of \eqref{eq:ah} and Theorem~\ref{thm:ConcreteEnergyError},
 we have
\begin{equation*}
 |\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} (u-u_h)|_{h,1}\leq \|u-u_h\|_h\lesssim
 \sqrt{\alpha_h}h^{\ell}|u|_{H^{\ell+1}({\Omega})},
\end{equation*}
 which together with  \eqref{eq:PODError}
 implies
\begin{equation*}
 |u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h|_{h,1}\leq |u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u|_{h,1}+|\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}(u-u_h)|_{h,1}
 \lesssim \sqrt{\alpha_h}h^{\ell}|u|_{H^{\ell+1}({\Omega})}.
\end{equation*}
 It follows from this estimate and
   \eqref{eq:PDZHOne} that
\begin{equation*}
  |u-\Pi_{k,h}^0 u|_{h,1}\leq |u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u|_{h,1}+|\Pi_{k,h}^0(\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u-u)|_{h,1}
  \lesssim |u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u|_{h,1}\lesssim
  \sqrt{\alpha_h}h^{\ell}|u|_{H^{\ell+1}(\O)}.
\end{equation*}
\par
 Finally we have, by \eqref{eq:HOnePODID}, \eqref{eq:HOneVEM},
 \eqref{eq:InterpolationError}
 and Theorem~\ref{thm:ConcreteEnergyError},
\begin{align*}
 |u-u_h|_{H^1({\Omega})}&\leq |u-I_{k,h} u|_{H^1({\Omega})}+|I_{k,h} u-u_h|_{H^1({\Omega})}\\
 &\lesssim |u-I_{k,h} u|_{H^1({\Omega})}+\sqrt{\alpha_h}\|I_{k,h} u-u_h\|_h\\
  &\leq |u-I_{k,h} u|_{H^1({\Omega})}+\sqrt{\alpha_h}
  \big(\|u-I_{k,h} u\|_h+\|u-u_h\|_h\big)
  \lesssim \alpha_h h^{\ell}|u|_{H^{\ell+1}(\O)}.
\end{align*}
\end{proof}
\goodbreak
\begin{remark}\label{rem:HOneErrors}{\rm
 In the case where $S^D(\cdot,\cdot)=S^D_1(\cdot,\cdot)$,
 the estimates for $|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h|_{h,1}$ and $|u-\Pi_{k,h}^0 u_h|_{h,1}$
 are better than the
 estimate for $|u-u_h|_{H^1({\Omega})}$.
}
\end{remark}
\subsection{Error Estimates in the $L_2$ Norm}\label{subsec:L2Error}
 We begin with two lemmas involving $S^D(\cdot,\cdot)$.
\begin{lemma}\label{lem:SDEstimate}
 We have
\begin{equation}\label{eq:SDEstimate}
 \sum_{D\in\cT_h}
 S^D(\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D} \zeta,\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D} \zeta)\lesssim
 h^2|\zeta|_{H^2({\Omega})}^2\qquad\forall\,
  \zeta\in H^2({\Omega})\cap H^1_0({\Omega}).
\end{equation}
\end{lemma}
\begin{proof}
  It follows from \eqref{eq:Trace} (applied to the first
  order partial derivatives of
  $\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta$), Lemma~\ref{lem:InterpolationError}
  and \eqref{eq:SD2} that
\begin{align*}
   &\sum_{D\in\cT_h}
   S^D_2(\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D} \zeta,\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D} \zeta)
   =\sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}}\|\partial(\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D} \zeta)/\partial s\|_{L_2(\partial D)}^2\\
   &\hspace{60pt}\lesssim \sum_{D\in\cT_h} \big[|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta|_{{H^1(D)}}^2+
         \hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta|_{H^2(D)}^2\big]
         \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta|_{H^2({\Omega})}^2,
\end{align*}
 and we have
\begin{align*}
  &\sum_{D\in\cT_h}
  S^D_1(\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta,\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D} \zeta)\lesssim
 \sum_{D\in\cT_h}\|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D} \zeta\|_{L_\infty(\partial D)}^2\\
  &\hspace{30pt}\lesssim \sum_{D\in\cT_h}\big[\hspace{1pt}h_{{\scriptscriptstyle D}}^{-2}\|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D} \zeta\|_{L_2(D)}^2
  +|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D} \zeta|_{{H^1(D)}}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D} \zeta|_{H^2(D)}^2\big]\\
  &\hspace{30pt}\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta|_{H^2({\Omega})}^2
\end{align*}
 by \eqref{eq:Sobolev}, Lemma~\ref{lem:InterpolationError} and
 \eqref{eq:SD1}.
\end{proof}
\begin{lemma}\label{lem:POSD}
    Assuming that $u\in H^{\ell+1}({\Omega})$ for some $\ell$
    between $1$ and $k$, we have
\begin{equation}\label{eq:POSD}
  \sum_{D\in\cT_h} S^D(u_h-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h,u_h-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h)\lesssim
  \alpha_h^2 h^{2\ell}|u|_{H^{\ell+1}(\O)}^2.
\end{equation}
\end{lemma}
\begin{proof}
  This is a consequence of \eqref{eq:ah},
\eqref{eq:ProjectionError}, Theorem~\ref{thm:ConcreteEnergyError}
  and Theorem~\ref{thm:ComputableEnergyError}:
\begin{align*}
  &\sum_{D\in\cT_h} S^D(u_h-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h,u_h-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h)= \|u_h-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h\|_h^2\\
  &\hspace{100pt}\lesssim \|u_h-u\|_h^2+\|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u\|_h^2
  +\|\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}(u-u_h)\|_h^2\\
  &\hspace{100pt}=\|u_h-u\|_h^2+\|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u\|_h^2+|u-u_h|_{H^1({\Omega})}^2
  \lesssim \alpha_h^2 h^{2\ell}|u|_{H^{\ell+1}(\O)}^2.
\end{align*}
\end{proof}
\par
 We can now prove a consistency estimate.
\begin{lemma}\label{lem:ConsistencyError}
  Assuming that $u\in H^{\ell+1}({\Omega})$ for some $\ell$
  between $1$ and $k$, we have
\begin{equation*}
  a(u-u_h,I_{k,h}\zeta)\leq C{\alpha_h}
   h^{\ell+1}|u|_{H^{\ell+1}({\Omega})}
     |\zeta|_{H^2({\Omega})} \qquad \forall\,
     \zeta\in H^2({\Omega})\cap H^1_0({\Omega}),
\end{equation*}
 where $\alpha_h$ is defined in \eqref{eq:kappah} and
 the positive constant $C$
 only depends on $\rho$, $k$ and $N$.
\end{lemma}
\begin{proof}
 We have, by \eqref{eq:Poisson}, \eqref{eq:aDef}, \eqref{eq:POD1} and
 \eqref{eq:DiscreteProblem}--\eqref{eq:aD},
\begin{align*
 &a(u-u_h,I_{k,h}\zeta)=a(u,I_{k,h}\zeta)-\sum_{D\in\cT_h} a^D(u_h,I_{k,D}\zeta)\notag\\
           &\hspace{40pt}=(f,I_{k,h}\zeta)-\sum_{D\in\cT_h} a^D(u_h,\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta)
           -\sum_{D\in\cT_h} a^D(u_h,I_{k,D}\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta)\\
           &\hspace{40pt}=(f,I_{k,h}\zeta)-a_h(u_h,I_{k,D}\zeta)
           +\sum_{D\in\cT_h} S^D(u_h-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h,I_{k,D}\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta)\\
           &\hspace{80pt}+\sum_{D\in\cT_h} a^D(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h-u_h,I_{k,D}\zeta
           -\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta)\\
 &\hspace{40pt}=(f,I_{k,h}\zeta-\Xi_hI_{k,h}\zeta)
 +\sum_{D\in\cT_h} S^D(u_h-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h,I_{k,D}\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta)\\
      &\hspace{80pt}+\sum_{D\in\cT_h} a^D(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h-u_h,I_{k,D}\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta),
\end{align*}
 and the three terms on the right-hand side  can be estimated as follows.
\par
 We have
\begin{equation*}
 |(f,I_{k,h}\zeta-\Xi_hI_{k,h}\zeta)|\lesssim h^{\ell+1}|u|_{H^{\ell+1}(\O)}|\zeta|_{H^2({\Omega})}
\end{equation*}
 by \eqref{eq:Xih2},
\begin{equation*}
\sum_{D\in\mathcal{T}_h} {S^D(u_h-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h,I_{k,D}\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta)}
 \lesssim {\alpha_h}h^{\ell+1}|u|_{H^{\ell+1}(\O)}|\zeta|_{H^2({\Omega})}
\end{equation*}
 by Lemma~\ref{lem:SDEstimate} and Lemma~\ref{lem:POSD},
 and
\begin{align*}
& \sum_{D\in\mathcal{T}_h} a^D(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u-u_h,I_{k,D}\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} I_{k,D}\zeta)\\
    &\hspace{40pt}\leq \big[|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u-u|_{h,1}+|u-u_h|_{H^1({\Omega})}\big]
   \big[|I_{k,D}\zeta-\zeta|_{H^1({\Omega})}+|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta|_{h,1}\big]\\
     &\hspace{40pt}\lesssim {\alpha_h}
      h^{\ell+1}|u|_{H^{\ell+1}(\O)}|\zeta|_{H^2({\Omega})}
\end{align*}
 by Lemma~\ref{lem:InterpolationError}, \eqref{eq:PODError}
 and Theorem~\ref{thm:ComputableEnergyError}.
\end{proof}
\begin{theorem}\label{thm:uhLTwoError}
 Assuming $u\in H^{\ell+1}({\Omega})$ for some $\ell$ between $1$ and $k$,
 there exists a positive constant $C$, depending only on
 $\rho$, $N$ and $k$, such that
\begin{equation}\label{eq:uhLTwoError}
 \|u-u_h\|_{L_2({\Omega})}\leq C{\alpha_h}h^{\ell+1}|u|_{H^{\ell+1}({\Omega})},
\end{equation}
 where $\alpha_h$ is defined in \eqref{eq:kappah}.
\end{theorem}
\begin{proof}
 Let $\zeta\in H^1_0({\Omega})$ be defined by
\begin{equation*}
   a(v,\zeta)=(v,u-u_h)  \qquad\forall\,v\in H^1_0({\Omega}).
\end{equation*}
 We have
 %
\begin{equation}\label{eq:Duality1}
   \|u-u_h\|_{L_2({\Omega})}^2=a(u-u_h,\zeta)
                          = a(u-u_h,\zeta-I_{k,h}\zeta)+a(u-u_h,I_{k,h}\zeta),
\end{equation}
 and since ${\Omega}$ is convex,
\begin{equation}\label{eq:Duality2}
  \|\zeta\|_{H^2({\Omega})}\leq C_{\Omega} \|u-u_h\|_{L_2({\Omega})}
\end{equation}
 by elliptic regularity \cite{Grisvard:1985:EPN,Dauge:1988:EBV}.
\par
 The first term on the right-hand side of  \eqref{eq:Duality1} satisfies
\begin{equation}\label{eq:Duality3}
   a(u-u_h,\zeta-I_{k,h}\zeta)\leq |u-u_h|_{H^1({\Omega})}
   |\zeta-I_{k,h}\zeta|_{H^1({\Omega})}\\
          \lesssim  h|u-u_h|_{H^1({\Omega})}|\zeta|_{H^2({\Omega})}
\end{equation}
 by \eqref{eq:HOnePODID}, and then
 \eqref{eq:uhLTwoError} follows from
 Theorem~\ref{thm:ComputableEnergyError},
 Lemma~\ref{lem:ConsistencyError}
 and \eqref{eq:Duality1}--\eqref{eq:Duality3}.
\end{proof}
\par
 We have similar $L_2$ error estimates for the computable
 approximations  $\Pi_{k,h}^0 u_h$ and $\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h$.
\begin{theorem}\label{thm:ComputableLTwoErrors}
 Assuming $u\in H^{\ell+1}({\Omega})$ for some $\ell$ between $1$ and $k$,
  there exists a positive constant $C$, depending only on $\rho$,
  $N$ and $k$,  such that
\begin{equation*}
 \|u-\Pi_{k,h}^0 u_h\|_{L_2({\Omega})}
 +\|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h\|_{L_2({\Omega})}\leq C{\alpha_h}h^{\ell+1}|u|_{H^{\ell+1}({\Omega})},
\end{equation*}
 where $\alpha_h$ is defined in \eqref{eq:kappah}.
\end{theorem}
\begin{proof} The estimate for $\Pi_{k,h}^0 u_h$ follows from
 \eqref{eq:PDZLTwoError}, Theorem~\ref{thm:uhLTwoError}
 and the relation
\begin{align*}
  \|u-\Pi_{k,h}^0 u_h\|_{L_2({\Omega})}&\leq \|u-\Pi_{k,h}^0 u\|_{L_2({\Omega})}
  +\|\Pi_{k,h}^0(u-u_h)\|_{L_2({\Omega})}\\
     &\leq \|u-\Pi_{k,h}^0 u\|_{L_2({\Omega})}+\|u-u_h\|_{L_2({\Omega})}.
\end{align*}
\par
 For the approximation $\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h$, we have
\begin{equation}\label{eq:LTwoPODError1}
 \|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h\|_{L_2({\Omega})}\leq\|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}I_{k,h} u\|_{L_2({\Omega})}
 +\|\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}(I_{k,h} u-u_h)\|_{L_2({\Omega})}
\end{equation}
 and, in view of \eqref{eq:Trace}, \eqref{eq:tbarNormDef},
  and Lemma~\ref{lem:PODStability2},
\begin{align}\label{eq:LTwoPODError2}
\|\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}(I_{k,h} u-u_h)\|_{L_2({\Omega})}^2&\lesssim
\sum_{D\in\mathcal{T}_h}|\!|\!| I_{k,D} u-u_h|\!|\!|_{k,D}^2\notag\\
    &\lesssim \sum_{D\in\mathcal{T}_h}\big(\hspace{1pt}h_{{\scriptscriptstyle D}}\|I_{k,D} u-u_h\|_{L_2(\partial D)}^2
    +\|\Pi_{k-2,D}^0 (I_{k,D} u-u_h)\|_{L_2(D)}^2\big)\\
    &\lesssim\sum_{D\in\mathcal{T}_h} \big(\|I_{k,D} u-u_h\|_{L_2(D)}^2
    +\hspace{1pt}h_{{\scriptscriptstyle D}}^2|I_{k,D} u-u_h|_{{H^1(D)}}^2\big)\notag\\
    &\lesssim \|u-u_h\|_{L_2({\Omega})}^2+ \|u-I_{k,h} u\|_{L_2({\Omega})}^2
    +h^2|u-I_{k,h} u|_{H^1({\Omega})}^2\notag\\
       &\hspace{50pt}+h^2|u-u_h|_{H^1({\Omega})}^2.\notag
\end{align}
\par
 The estimate for $\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h$ follows from
   Lemma~\ref{lem:InterpolationError},
 Theorem~\ref{thm:ComputableEnergyError},
 Theorem~\ref{thm:uhLTwoError}
 and \eqref{eq:LTwoPODError1}--\eqref{eq:LTwoPODError2}.
\end{proof}
\subsection{Error Estimates for $u_h$ in the $L_\infty$ Norm}\label{subsec:uhLInftyError}
 Here we consider a $L_\infty$ error estimate for $u_h$ over the edges of $\mathcal{T}_h$,
 where $u_h$ is computable.
 We will treat the two choices of $S^D(\cdot,\cdot)$ separately.  The set of all the
 edges in $\mathcal{T}_h$ will be denoted by $\mathcal{E}_h$.
\subsubsection{The Case where $S^D(\cdot,\cdot)=S^D_2(\cdot,\cdot)$}
\label{eq:subsubsec:uhLInftySD2}
 We have the following result for this choice of $S^D(\cdot,\cdot)$.
\begin{theorem}\label{thm:MaxNormBdryEst2}
 Assuming that
 the solution $u$ of \eqref{eq:Poisson}
  belongs to $\in H^{\ell+1}({\Omega})$ for some $\ell$ between $1$ and $k$,
 we have
\begin{equation*}
  \max_{e\in\mathcal{E}_h}\|u-u_h\|_{L_\infty(e)}\leq C
  h^{\ell}|u|_{H^{\ell+1}(\O)},
\end{equation*}
 where the positive constant $C$ only depends on $\rho$, $N$ and  $k$.
\end{theorem}
\begin{proof} First we observe that,
  by \eqref{eq:SD2}, \eqref{eq:ah} and Theorem~\ref{thm:ConcreteEnergyError},
\begin{equation}\label{eq:EdgeEstimate}
 \sum_{D\in\cT_h}\sum_{e\in\cE_{\ssD}}\hspace{1pt}h_{{\scriptscriptstyle D}}\|\partial[(u-u_h)-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} (u-u_h)]/\partial s\|_{L_2(e)}^2
 \lesssim \|u-u_h\|_h^2\lesssim h^{2\ell}|u|_{H^{\ell+1}(\O)}^2.
\end{equation}
\par
 We can connect any point in $e\in\mathcal{E}_h$
 to $\partial{\Omega}$, where $u-u_h=0$, by a path along the edges in $\mathcal{E}_h$.  Therefore it follows from
 a direct calculation (or a Sobolev inequality in one variable)
 that
\begin{align}\label{eq:PathEstimate}
  \|u-u_h\|_{L_\infty(e)}^2&\lesssim \sum_{D\in\cT_h}\sum_{e\in\cE_{\ssD}}
   h_e\|\partial(u-u_h)/\partial s\|_{L_2(e)}^2\notag\\
    &\lesssim \sum_{D\in\cT_h}\sum_{e\in\cE_{\ssD}}\hspace{1pt}h_{{\scriptscriptstyle D}}\|\partial[(u-u_h)-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} (u-u_h)]/\partial s\|_{L_2(e)}^2\\
    &\hspace{40pt}+\sum_{D\in\cT_h}\sum_{e\in\cE_{\ssD}}
   \hspace{1pt}h_{{\scriptscriptstyle D}}\|\partial[\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-u_h)]/\partial s\|_{L_2(e)}^2\qquad\qquad\forall\,e\in\mathcal{E}_h,\notag
\end{align}
 and we have,
 by \eqref{eq:Trace}, \eqref{eq:DiscreteEstimate1},
 \eqref{eq:PODStability1} and Theorem~\ref{thm:ComputableEnergyError},
\begin{align}\label{eq:TrickyEstimate}
  &\sum_{D\in\cT_h}\sum_{e\in\cE_{\ssD}} \hspace{1pt}h_{{\scriptscriptstyle D}}\|\partial[\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-u_h)]/\partial s\|_{L_2(e)}^2\notag\\
  &\hspace{80pt}\lesssim \sum_{D\in\cT_h} \big(|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-u_h)|_{H^1(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}^2|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-u_h)|_{H^2(D)}^2\big)\\
  &\hspace{80pt}\lesssim\sum_{D\in\cT_h} |\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-u_h)|_{H^1(D)}^2\lesssim  h^{2\ell}|u|_{H^{\ell+1}(\O)}^2.
  \notag
\end{align}
\par
 The estimates \eqref{eq:EdgeEstimate}--\eqref{eq:TrickyEstimate} together imply
\begin{equation*}
   \|u-u_h\|_{L_\infty(e)}\leq h^\ell|u|_{H^{\ell+1}(\O)} \qquad\forall\,e\in\mathcal{E}_h.
\end{equation*}
\end{proof}
\subsubsection{The Case where $S^D(\cdot,\cdot)=S^D_1(\cdot,\cdot)$}
\label{eq:subsubsec:uhLInftySD1}
 We will establish an analog of Theorem~\ref{thm:MaxNormBdryEst2} under the
 additional assumption that
 $\mathcal{T}_h$ is quasi-uniform, i.e., there exists a positive constant
 $\gamma$  independent of $h$ such that
\begin{equation}\label{eq:gamma}
  \hspace{1pt}h_{{\scriptscriptstyle D}}\geq \gamma h \qquad\forall\,\gamma\in\mathcal{T}_h.
\end{equation}
\begin{theorem}\label{thm:MaxNormBdryEst1}
 Assuming $\mathcal{T}_h$ is quasi-uniform and the solution $u$ of
 \eqref{eq:Poisson} belongs to $H^{\ell+1}({\Omega})$ for some $\ell$ between $1$ and $k$,
 we have
\begin{equation*}
  \max_{e\in\mathcal{E}_h}\|u-u_h\|_{L_\infty(e)}\leq
  C\ln(1+\max_{D\in\mathcal{T}_h}\tau_{\!\ssD})
  h^{\ell}|u|_{H^{\ell+1}(\O)},
\end{equation*}
 where the positive constant $C$ only depends on $\rho$,
 $N$, $\gamma$ and  $k$.
\end{theorem}
\begin{proof} Let $D\in\mathcal{T}_h$ be arbitrary.
 First we observe that,
 by \eqref{eq:SD1}, Remark~\ref{rem:NormEquivalence},
 \eqref{eq:ah}, \eqref{eq:InterpolationError},
 Theorem~\ref{thm:ConcreteEnergyError} and
 Theorem~\ref{thm:ComputableEnergyError},
\begin{align}\label{eq:SD1Edge1}
 &\|(I_{k,D} u-u_h)-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(I_{k,D} u-u_h)\|_{L_\infty(\partial D)}\notag\\
 &\hspace{80pt}\lesssim \|I_{k,h} u-u_h\|_h\\
 &\hspace{80pt}
 \lesssim \|I_{k,h} u-u\|_h+\|u-u_h\|_h
  \lesssim [\ln(1+\max_{D\in\mathcal{T}_h}\tau_{\!\ssD})]^\frac12h^\ell|u|_{H^{\ell+1}(\O)}.\notag
\end{align}
\par
 Furthermore, it follows from \eqref{eq:G2}, \eqref{eq:DiscreteEstimate3},
 \eqref{eq:PODStability1}, \eqref{eq:LTwoPODID},
  Theorem~\ref{thm:ConcreteEnergyError},
 Theorem~\ref{thm:ComputableLTwoErrors} and \eqref{eq:gamma} that
\begin{align}\label{eq:SD1Edge2}
  &\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(I_{k,D} u-u_h)\|_{L_\infty(\partial D)}\lesssim
  \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(I_{k,D} u-u_h)\|_{L_2(D)}+
   |\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(I_{k,D} u-u_h)|_{H^1(D)}\notag\\
   &\hspace{60pt}\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\big(\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D} u-u\|_{L_2(D)}
   +\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h\|_{L_2(D)}\big)+|I_{k,D} u-u_h|_{H^1(D)}\\
   &\hspace{60pt}\lesssim \ln(1+\max_{D\in\mathcal{T}_h}\tau_{\!\ssD}) h^\ell|u|_{H^{\ell+1}(\O)}.\notag
\end{align}
\par
 The theorem follows from \eqref{eq:1DInterpolationError},
 \eqref{eq:SD1Edge1}, \eqref{eq:SD1Edge2} and the triangle inequality.
\end{proof}
\subsection{Error Estimates for $\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h$ and $\Pi_{k,h}^0 u_h$ in the $L_\infty$ Norm}
\label{subsec:ComputableLInftyError}
\par
 Again we treat the two choices of $S^D(\cdot,\cdot)$ separately.
\subsubsection{The Case where $S^D(\cdot,\cdot)=S^D_2(\cdot,\cdot)$}
\label{eq:subsubsec:ComputableLInftySD2}
 For this choice of $S^D(\cdot,\cdot)$, we can establish the
 following result without
 assuming that $\mathcal{T}_h$ is quasi-uniform.
\begin{theorem}\label{thm:LInftyErrorsSD2}
  Assuming the solution $u$ of \eqref{eq:Poisson} belongs to
  $H^{\ell+1}({\Omega})$ for some $\ell$ between $1$ and $k$,
 there exists a positive constant $C$, depending only on $\rho$,
 $N$ and $k$, such that
\begin{equation*}
 \|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h\|_{L_\infty({\Omega})}+\|u-\Pi_{k,h}^0 u_h\|_{L_\infty({\Omega})}
  \leq C  h^{\ell}|u|_{H^{\ell+1}(\O)}.
\end{equation*}
\end{theorem}
\begin{proof}  For any $D\in\mathcal{T}_h$, we have,  by \eqref{eq:Sobolev},
 \eqref{eq:DiscreteEstimate1} and \eqref{eq:PODStability1},
\begin{align*}
 \|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h\|_{{L_\infty(D)}}
 &\lesssim
 \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h\|_{L_2(D)}+|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h|_{H^1(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h|_{H^2(D)}\\
 &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h\|_{L_2(D)}+|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h|_{H^1(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u|_{H^2(D)}\\
             &\hspace{50pt}+  \hspace{1pt}h_{{\scriptscriptstyle D}}|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} (u-u_h)|_{H^2(D)}\\
  &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h\|_{L_2(D)}+|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h|_{H^1(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u|_{H^2(D)}\\
  &\hspace{40pt}+|u-u_h|_{H^1(D)},
\end{align*}
 and
\begin{align*}
  \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h\|_{L_2(D)}&\leq \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|u-u_h\|_{L_2(D)}+
      \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|u_h-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h\|_{L_2(D)}\\
      &\lesssim \big(\|u-u_h\|_{L_\infty(\partial D)}+|u-u_h|_{H^1(D)}\big)
        +|u_h-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h|_{H^1(D)}
\end{align*}
 by \eqref{eq:G2}, \eqref{eq:PF2} and \eqref{eq:POD2}.  These two
 estimates together with
 Lemma~\ref{lem:PODErrors},
 Theorem~\ref{thm:ComputableEnergyError} and
 Theorem~\ref{thm:MaxNormBdryEst2} imply
 the estimate for $u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h$.
\par
 The estimate for $u-\Pi_{k,h}^0 u_h$ can be derived similarly.
 We have, by \eqref{eq:Sobolev},
 \eqref{eq:DiscreteEstimate1} and
 \eqref{eq:PDZHOne},
\begin{align*}
 \|u-\Pi_{k,D}^0\hspace{1pt} u_h\|_{{L_\infty(D)}}
 &\lesssim
 \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|u-\Pi_{k,D}^0\hspace{1pt} u_h\|_{L_2(D)}+|u-\Pi_{k,D}^0\hspace{1pt} u_h|_{H^1(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}|u-\Pi_{k,D}^0\hspace{1pt} u|_{H^2(D)}\\
 &\hspace{40pt}+|u-u_h|_{H^1(D)},
\end{align*}
 and
\begin{align*}
  \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|u-\Pi_{k,D}^0\hspace{1pt} u_h\|_{L_2(D)}&\leq \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|u-u_h\|_{L_2(D)}+
           \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|u_h-\Pi_{k,D}^0\hspace{1pt} u_h\|_{L_2(D)}\\
  &\lesssim \big(\|u-u_h\|_{L_\infty(\partial D)}+|u-u_h|_{H^1(D)}\big)+
  |u_h-\Pi_{k,D}^0\hspace{1pt} u_h|_{H^1(D)}
\end{align*}
 by \eqref{eq:G2}, \eqref{eq:PF1} and \eqref{eq:PF2}.
 The estimate for $u-\Pi_{k,h}^0 u_h$ now follows from
 Lemma~\ref{lem:PDZErrors},
  Theorem~\ref{thm:ComputableEnergyError} and
 Theorem~\ref{thm:MaxNormBdryEst2}.
\end{proof}
\subsubsection{The Case where $S^D(\cdot,\cdot)=S^D_1(\cdot,\cdot)$}
\label{eq:subsubsec:ComputableLInftySD1}
 The following analog of Theorem~\ref{thm:LInftyErrorsSD2} is
 proved by the same arguments but with
 Theorem~\ref{thm:MaxNormBdryEst2} replaced by
 Theorem~\ref{thm:MaxNormBdryEst1}.
\begin{theorem}\label{thm:LInftyErrorsSD1}
  Assuming $\mathcal{T}_h$ is quasi-uniform and the solution $u$ of \eqref{eq:Poisson}
  belongs to $H^{\ell+1}({\Omega})$ for some $\ell$ between $1$ and $k$,
 there exists a positive constant $C$, depending only on $\rho$,
  $N$, $\gamma$ and $k$, such that
\begin{equation*}
 \|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h\|_{L_\infty({\Omega})}+\|u-\Pi_{k,h}^0 u_h\|_{L_\infty({\Omega})}
  \leq C \ln(1+\max_{D\in\mathcal{T}_h}\tau_{\!\ssD})  h^{\ell}|u|_{H^{\ell+1}(\O)}.
\end{equation*}
\end{theorem}
\section{Virtual Element Methods for the Poisson Problem in
Three Dimensions}\label{sec:Poisson3D}
 The analysis of virtual element methods in three dimensions
 follows the same strategy as in two dimensions
 and many of the results in Section~\ref{sec:LocalVEM2D} and
 Section~\ref{sec:Poisson2D} carry over by identical arguments.
 We will only provide details for estimates that require different derivations.
\par
 Let $\mathcal{T}_h$ be a polyhedral mesh on ${\Omega}$. The set of the faces of a
 subdomain $D\in\mathcal{T}_h$ is denoted by
 $\mathcal{F}_\ssD$ and the set of the edges of
 $F$ is denoted by $\cE_{\ssF}$.    The set of all the faces of
 $\mathcal{T}_h$ is denoted by $\mathcal{F}_h$ and the set of all the edges of
  $\mathcal{T}_h$ is denoted by
 $\mathcal{E}_h$.
\subsection{Shape Regularity Assumptions in Three Dimensions}
\label{subsec:Shape3D}
 We impose the following shape regularity assumptions on $\mathcal{T}_h$,
 where $\hspace{1pt}h_{{\scriptscriptstyle D}}$ is the diameter of $D$.
\par\smallskip\noindent
{\em Assumption 1} \quad
 There exists $\rho\in (0,1)$, independent of $h$,
  such that every polyhedron $D\in\mathcal{T}_h$ is star-shaped with
 respect to a ball $\fB_{\!\ssD}$ with radius $\geq \rho \hspace{1pt}h_{{\scriptscriptstyle D}}$.
\par\smallskip\noindent
 {\em Assumption 2} \quad There exists a positive integer $N$,
 independent of $h$, such that
  $|\mathcal{F}_\ssD|\leq N$ for all $D\in\mathcal{T}_h$.
 \par\smallskip\noindent
 {\em Assumption 3}\quad The shape regularity assumptions in
 Section~\ref{subsec:GlobalShapeRegularity}
 are satisfied by all the faces in $\mathcal{F}_h$, with the same $\rho$ from
 Assumption~1 and the same
 $N$ from Assumption~2.
\par\smallskip
 All the hidden constants below will only depend on $\rho$, $N$ and $k$.
 \par
 Let $D$ be a polyhedron in $\mathcal{T}_h$.
 We can define the inner product $(\!(\cdot,\cdot)\!)$ by
 \eqref{eq:InnerProduct} where the infinitesimal arc-length $ds$
 is replaced by
 the infinitesimal
 surface area $dS$.  Then the projection operator
 $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}:H^1(D)\longrightarrow\mathbb{P}_k(D)$ is defined by
 \eqref{eq:PODDef} and
\begin{equation}\label{eq:3DPODStability1}
  |\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta|_{H^1(D)}\leq |\zeta|_{H^1(D)}\qquad\forall\,\zeta\in H^1(D).
\end{equation}
 The projection from ${L_2(D)}$ to $\mathbb{P}_k(D)$ is again denoted by $\Pi_{k,D}^0\hspace{1pt}$.
\par
 The results in Section~\ref{sec:SS}
  are valid for $D\in\mathcal{T}_h$ under Assumption 1.
  Consequently the results in Section~\ref{subsec:PODEstimates} and
  Section~\ref{subsec:PDZEstimates}  are also valid
  provided the semi-norm $|\!|\!|\cdot|\!|\!|_{k,D}$ is defined by the
  following analog of
  \eqref{eq:tbarNormDef}:
\begin{equation}\label{eq:tbarNorm3D}
  |\!|\!| \zeta|\!|\!|_{k,D}^2=\|\Pi_{k-2,D}^0 \zeta\|_{L_2(D)}^2
  +\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\mathcal{F}_\ssD}\|\Pi_{k-1,F}^0 \zeta\|_{L_2(F)}^2,
\end{equation}
 where $\Pi_{k-1,F}^0$ is the projection from $L_2(F)$ onto $\mathbb{P}_{k-1}(F)$.
\par
 We have
 the following estimates for $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}$ and  $\Pi_{k,D}^0\hspace{1pt}$:
\begin{equation}\label{eq:3DLTwoPODPDZ}
   \|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} \zeta\|_{L_2(D)}+ \|\zeta-\Pi_{k,D}^0\hspace{1pt}\zeta\|_{L_2(D)}
   \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell+1}|\zeta|_{H^{\ell+1}(D)}
 \quad\forall\,\zeta\in H^{\ell+1}(D),\,0\leq\ell\leq k,
\end{equation}
 and for $1\leq\ell \leq k$,
\begin{alignat}{3}
   |\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}\zeta|_{{H^1(D)}}+ |\zeta-\Pi_{k,D}^0\hspace{1pt}\zeta|_{{H^1(D)}}&\lesssim
   \hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell}|\zeta|_{H^{\ell+1}(D)}
 &\quad&\forall\,\zeta\in H^{\ell+1}(D), \label{eq:3DHOnePODPDZ}\\
  |\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} \zeta|_{H^2(D)}+|\zeta-\Pi_{k,D}^0\hspace{1pt}\zeta|_{H^2(D)}&
  \lesssim h_D^{\ell-1}|\zeta|_{H^{\ell+1}(D)}
 &\quad&
  \forall\,\zeta\in H^{\ell+1}(D).\label{eq:3DHTwoPODPDZ}
\end{alignat}
\par
 The analogs of Lemma~\ref{lem:PODStability2} and
 \eqref{eq:tbarNorm3D} lead to the
 estimate
\begin{equation}\label{eq:3DPODStability2}
  \|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} \zeta\|_{L_2(D)}^2\lesssim |\!|\!|\zeta|\!|\!|_{k,D}^2
  \lesssim \|\zeta\|_{L_2(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}\|\zeta\|_{L_2(\partial D)}^2\qquad
  \forall\,\zeta\in H^1(D),
\end{equation}
 and we also have the following analog of \eqref{eq:PDZHOne}:
\begin{equation}\label{eq:3DPDZHOne}
 |\Pi_{k,D}^0\hspace{1pt} \zeta|_{H^1(D)}\lesssim |\zeta|_{H^1(D)} \qquad\forall\,\zeta\in{H^1(D)}.
\end{equation}
\subsection{The Local Virtual Element Space $\bm\cQ^k(D)$}
\label{subsec:Local3D}
 The space $\cQ^k(\p D)$ of
 continuous piecewise (two dimensional) virtual element functions of order
 $\leq k$ on $\partial D$ is defined by
\begin{equation}\label{eq:cQbDDef}
  \cQ^k(\p D)=\{v\in C(\partial D):\,v\big|_F\in \mathcal{Q}^k(F) \quad\forall\,F\in \mathcal{F}_\ssD\}.
\end{equation}
\par
 For $k\geq 1$,
 the virtual element space $\cQ^k(D)\subset {H^1(D)}$ is defined by the
 following conditions:  $v\in {H^1(D)}$ belongs to $\cQ^k(D)$ if and only if
 (i) the trace of $v$ on $\partial D$ belongs to $\cQ^k(\p D)$,
  (ii) the distribution $-\Delta v$ belongs to $\mathbb{P}_k(D)$,
 and (iii) condition \eqref{eq:Condition3} is satisfied.
\begin{remark}\label{rem:3DContinuity}{\rm
  Since the restriction of $v\in\cQ^k(D)$ to $\partial D$ belongs to
  $C(\partial D)$ and $-\Delta v\in \mathbb{P}_k(D)$, the virtual
  element function $v$ is also continuous
 on $\bar D$ (cf. \cite[Section~1.2]{Kenig:1994:CBMS}).}
\end{remark}
\begin{remark}\label{rem:3Ddofs}{\rm
 The degrees of freedom of $\cQ^k(D)$ (cf. \cite{AABMR:2013:Projector})
 consist of (i) the values of $v$ at the
 vertices of $D$ and nodes on
 the interior of each edge of $D$ that determine a polynomial of degree $k$
 on each edge of $D$,
 (ii) the moments of $\Pi_{k-2,F}^0 v$ on each face $F$ of $D$, and
 (iii) the moments of
 $\Pi_{k-2,D}^0 v$ on $D$.}
\end{remark}
\begin{remark}\label{rem:3DComputable}{\rm
  For $v\in\cQ^k(D)$ and $F\in\mathcal{F}_\ssD$, the polynomial
  $\Pi_{k,F}^0 v$ can be computed in terms of the degrees of freedom of $v\big|_F$
   (cf. Remark~\ref{rem:Computable}).  Therefore the polynomial
  $\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v$ can be computed in terms of the degrees of freedom of
  $v\in\cQ^k(D)$ through \eqref{eq:POD1} and \eqref{eq:POD2}.
  The polynomial $\Pi_{k,D}^0\hspace{1pt} v$ can then be computed through
  \eqref{eq:Condition3}.
}
\end{remark}
\begin{remark}\label{rem:DAndF}{\rm
 Under Assumption~3 in Section~\ref{subsec:Shape3D}, the results in
 Section~\ref{sec:LocalVEM2D} (with $D$ replaced by $F$)
 are valid for the restriction of $v\in \cQ^k(D)$
 to any face $F$ of $D$.
 }
 \end{remark}
\par
 The three dimensional analogs of Lemma~\ref{lem:MEP},
 Lemma~\ref{lem:Fundamental}
 and Lemma~\ref{lem:InverseEstimate1}
 lead to the estimate
\begin{equation}\label{eq:3DInverse}
 |v|_{H^1(D)}^2\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-2}|\!|\!|\zeta|\!|\!|_{k,D}^2
      +\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\mathcal{F}_\ssD}\|\nabla_{\!\!F} v\|_{L_2(F)}^2 \qquad\forall\,v\in\cQ^k(D),
\end{equation}
 where $\nabla_{\!\!F}$ is the two dimensional gradient operator on the face $F$,
 and we also have an analog  of \eqref{eq:KerPod13}:
\begin{equation}\label{eq:3DKerPOD}
  \|\Pi_{k-2,D}^0 v\|_{L_2(D)}^2\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD} \|\Pi_{k-1,F}^0 v\|_{L_2(F)}^2
  \qquad\forall\,v\in \mathcal{N}(\POD).
\end{equation}
 Hence we have, by \eqref{eq:tbarNorm3D}, \eqref{eq:3DInverse} and
 \eqref{eq:3DKerPOD},
\begin{equation}\label{eq:3DInverse2}
  |v|_{H^1(D)}^2\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\sum_{F\in\mathcal{F}_\ssD}\|\Pi_{k-1,F}^0 v\|_{L_2(F)}^2\\
         +\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\mathcal{F}_\ssD}\|\nabla_{\!\!F} v\|_{L_2(F)}^2\qquad\forall\,v\in\mathcal{N}(\POD).
\end{equation}
\par
 The interpolation operator $I_{k,D}:H^2(D)\longrightarrow \cQ^k(D)$ is defined by the
 condition that $I_{k,D}\zeta$ and $\zeta$ share the same degrees of freedom.
 In particular we have
\begin{equation}\label{eq:3DIDInvariance}
 I_{k,D} q=q\qquad\forall\,q\in \mathbb{P}_k(D).
\end{equation}
\par
 Note that
\begin{equation}\label{eq:ConsistentID}
  I_{k,F} (\zeta\big|_F)=(I_{k,D}\zeta)\big|_F \qquad\forall\,F\in\mathcal{F}_\ssD,
\end{equation}
 and hence, in view of \eqref{eq:G2}, Lemma~\ref{lem:MaximumPrinciple}
 and \eqref{eq:tbarNorm3D},
\begin{align}\label{eq:3DIDtbar}
  |\!|\!| I_{k,D}\zeta|\!|\!|_{k,D}^2&=\|\Pi_{k-2,D}^0 (I_{k,D}\zeta)\|_{L_2(D)}^2
  +\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD}\|\Pi_{k-1,F}^0 (I_{k,D}\zeta)\|_{L_2(F)}^2\notag\\
  &=\|\Pi_{k-2,D}^0\zeta\|_{L_2(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD} \|\Pi_{k-1,F}^0(I_{k,F}\zeta)\|_{L_2(F)}^2\\
  &\lesssim \|\zeta\|_{L_2(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}\SumFh_F^2\|I_{k,F}\zeta\|_{L_\infty(F)}^2\notag\\
  &\lesssim \|\zeta\|_{L_2(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}\SumFh_F^2\big(\|I_{k,F}\zeta\|_{L_\infty(\partial F)}^2
  +\|\nabla_{\!\!F}(I_{k,F}\zeta)\|_{L_2(F)}^2\big).\notag
\end{align}
\par
 The error estimates for $I_{k,D}$ rely on the following analog of \eqref{eq:IDHOne},
 where $\tau_{\!\ssF}$ is defined by replacing $D$ by $F$ in \eqref{eq:TauD}.
\begin{lemma}\label{lem:3DIDHOne}
 We have
\begin{equation}\label{eq:3DIDHOne}
 |I_{k,D}\zeta|_{H^1(D)}\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}\|\zeta\|_{L_2(D)}
 +|\zeta|_{H^1(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{H^2(D)}
\end{equation}
 for all $\zeta\in{H^2(D)}$.
\end{lemma}
\begin{proof}
 Let $\zeta \in H^2(D)$ be arbitrary.  It follows from
 \eqref{eq:3DInverse}
 and \eqref{eq:3DIDtbar} that
\begin{align*}
  |I_{k,D}\zeta|_{H^1(D)}^2
  &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-2}\|\zeta\|_{L_2(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}} \sum_{F\in\FD} \|I_{k,F} \zeta\|_{L_\infty(\partial F)}^2
  +\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD} \|\nabla_{\!\!F}(I_{k,F}\zeta)\|_{L_2(F)}^2.
\end{align*}
 We have, by \eqref{eq:Sobolev} and \eqref{eq:TrivialBdd},
\begin{align}\label{eq:IFLinfty}
  \hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD}\|I_{k,F}\zeta\|_{L_\infty(\partial F)}^2&\lesssim
  \hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD} \|\zeta\|_{L_\infty(\partial F)}^2\\
     &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}\|\zeta\|_{L_\infty(D)}^2\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-2}\|\zeta\|_{L_2(D)}^2
     +|\zeta|_{H^1(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta|_{H^2(D)}^2,\notag
\end{align}
 and by \eqref{eq:HalfTrace}, Lemma~\ref{lem:CZ2} and \eqref{eq:IDHOneHalf},
\begin{align}\label{eq:IFHone}
  \hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD} \|\nabla_{\!\!F}(I_{k,F}\zeta)\|_{L_2(F)}^2&\lesssim
    \hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD}\big(|\zeta|_{H^1(F)}^2+h_F |\zeta|_{H^{3/2}(F)}^2\big)\\
   &\lesssim \sum_{D\in\cT_h}\big( |\zeta|_{H^1(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}^{2}|\zeta|_{H^2(D)}^2\big).
   \notag
\end{align}
\end{proof}
\par
 Note that \eqref{eq:3DIDtbar}, \eqref{eq:IFLinfty} and \eqref{eq:IFHone} imply
\begin{equation}\label{eq:3DIDtbar2}
  \tbarI_{k,D}\zeta|\!|\!|_{k,D}\lesssim
  \|\zeta\|_{L_2(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{H^1(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta|_{H^2(D)} \quad\forall\,\zeta\in{H^2(D)},
\end{equation}
 and hence we have, in view of \eqref{eq:3DLTwoPODPDZ},
 \eqref{eq:3DPODStability2}, \eqref{eq:3DIDtbar} and \eqref{eq:3DIDHOne},
\begin{align}\label{eq:3DIDLTwo}
  \|I_{k,D}\zeta\|_{L_2(D)}&\leq \|I_{k,D}\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta\|_{L_2(D)}+\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta\|_{L_2(D)}\notag\\
    &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}|I_{k,D}\zeta|_{H^1(D)}+\tbarI_{k,D}\zeta|\!|\!|_{k,D}\\
    &\lesssim \big(\|\zeta\|_{L_2(D)}
 +\hspace{1pt}h_{{\scriptscriptstyle D}}|\zeta|_{H^1(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta|_{H^2(D)}\big)\notag
\end{align}
 for all $\zeta\in H^2(D)$.
\par
 In view of \eqref{eq:3DIDInvariance},
 the following analogs of \eqref{eq:LTwoPODID}--\eqref{eq:HTwoPODID},
 where $\zeta\in H^{\ell+1}(D)$ and $1\leq\ell\leq k$,
  can be obtained by
 combining the Bramble-Hilbert estimates \eqref{eq:BHEstimates} with the stability
 estimates \eqref{eq:3DPODStability1},
 \eqref{eq:3DPODStability2}, \eqref{eq:3DIDHOne},
 \eqref{eq:3DIDtbar2} and \eqref{eq:3DIDLTwo}.
\begin{alignat}{3}
 \|\zeta-I_{k,D}\zeta\|_{L_2(D)}+\|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta\|_{L_2(D)}
 &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell+1}|\zeta|_{H^{\ell+1}(D)}
     \label{eq:3DLTwoPODID}\\
  |\zeta-I_{k,D}\zeta|_{{H^1(D)}}+|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta|_{{H^1(D)}}
 &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell}|\zeta|_{H^{\ell+1}(D)}
   \label{eq:3DHOnePODID}\\
    |\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta|_{H^2(D)}&\lesssim
    \hspace{1pt}h_{{\scriptscriptstyle D}}^{\ell-1}|\zeta|_{H^{\ell+1}(D)}
  \label{eq:3DHTwoPODID}
\end{alignat}
\begin{remark}\label{rem:3DLInftyInterpolation}{\rm
  We also have the following analog of \eqref{eq:1DInterpolationError}:
\begin{equation*}
  \|\zeta-I_{k,D}\zeta\|_{L_\infty(D)}\lesssim h^{\ell-\frac12}|u|_{H^{\ell+1}(D)}
\end{equation*}
 for all $\zeta\in H^{\ell+1}(D)$ and $1\leq \ell\leq k$.  The proof uses
 Lemma~\ref{lem:MaximumPrinciple} (which is valid in three dimensions) and the
 arguments for \eqref{eq:1DInterpolationError}.  But we do not need this estimate in
 the error analysis.
}
\end{remark}
\subsection{The Discrete Problem}\label{subsec:DP3D}
 Let the global virtual element space $\cQ^k_h$ be defined by
 $$\cQ^k_h=\{v\in H^1_0({\Omega}):\,v\big|_D\in\cQ^k(D)\quad\forall\,D\in\mathcal{T}_h\}.$$
 The discrete problem for \eqref{eq:Poisson} is to find $u_h\in\cQ^k_h$ such that
\begin{equation*}
  a_h(u_h,v)=(f,\Xi_h v) \qquad\forall\,v\in\cQ^k_h,
\end{equation*}
 where $\Xi_h$ is defined as in \eqref{eq:Xih},
\begin{equation*}
  a_h(w,v)=\sum_{D\in\cT_h}\big[\big(\nabla(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} w),\nabla(\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)\big)_{L_2(D)}+S^D(w-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} w,v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)\big],
\end{equation*}
 and the local stabilizing bilinear form $S^D(\cdot,\cdot)$ is given by
\begin{align}
    S^D(w,v)&=\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD}\Big(h_F^{-2}(\Pi_{k-2,F}^0 w,\Pi_{k-2,F}^0 v)_{L_2(F)}+
   \sum_{p\in\mathcal{N}_{\p F}}w(p)v(p)\Big).\label{eq:3DSD}
\end{align}
 Here $\mathcal{N}_{\p F}$ is the set of the nodes along $\partial F$ associated with the degrees of freedom
 of a virtual element function.
\begin{lemma}\label{lem:3DSDBdd}
 There exists a positive constant $C$, depending only on $\rho$, $N$ and $k$, such that
\begin{alignat}{3}
 |v|_{H^1(D)}^2&\leq C \big [\ln(1+\max_{F\in\mathcal{F}_\ssD}\tau_{\!\ssF})\big] S^D(v,v)
 &\qquad&\forall\,v\in\mathcal{N}(\POD). \label{eq:3DSDBdd}
\end{alignat}
\end{lemma}
\begin{proof}  Let $v\in\mathcal{N}(\POD)$ be arbitrary.  We have, by  \eqref{eq:G2},
 \eqref{eq:PF2}, Lemma~\ref{lem:MaximumPrinciple}, Corollary~\ref{cor:InverseEstimate3}
  and \eqref{eq:3DInverse2},
\begin{align*}
 |v|_{H^1(D)}^2&\lesssim  \hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD} \big(\hspace{1pt}h_{{\scriptscriptstyle D}}^{-2}\|v\|_{L_2(F)}^2
 +\|\nabla_{\!\!F} v\|_{L_2(F)}^2\big)\notag\\
      &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD}\Big(\hspace{1pt}h_{{\scriptscriptstyle D}}^{-2}h_F^2\big(\|v\|_{L_\infty(\partial F)}^2
      +\|\nabla_{\!\!F} v\|_{L_2(F)}^2\big)+\|\nabla_{\!\!F} v\|_{L_2(F)}^2\Big)\\
      &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD}\big(\|v\|_{L_\infty(\partial F)}^2+\|\nabla_{\!\!F} v\|_{L_2(F)}^2\big)\\
      &\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD}\big(h_F^{-2}\|\Pi_{k-2,F}^0 v\|_{L_2(F)}^2
      +\ln(1+\tau_{\!\ssF})\|v\|_{L_\infty(\partial F)}^2\big),
\end{align*}
  which together with Remark~\ref{rem:NormEquivalence} and \eqref{eq:3DSD}
  implies \eqref{eq:3DSDBdd}.
\end{proof}
\par
 It follows from \eqref{eq:3DSDBdd} that we have an analog of \eqref{eq:Stability}:
 \begin{equation}\label{eq:3DStability}
  |v|_{H^1({\Omega})}^2\leq 2\sum_{D\in\cT_h}\big[|\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} v|_{H^1({\Omega})}^2
  +|v-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v|_{{H^1(D)}}^2\big]\lesssim
    \beta_h a_h(v,v) \qquad\forall\,v\in \cQ^k_h,
\end{equation}
 where
\begin{equation}\label{eq:lambdah}
  \beta_h=     \ln(1+\max_{F\in\mathcal{F}_h}\tau_{\!\ssF}).
\end{equation}
 Hence the discrete problem is well-posed.
\begin{remark}\label{rem:3DInfo}{\rm
 The constants in the error estimates for the virtual element methods
  will only depend on
 $\rho$, $N$, $k$ and $\beta_h$. Therefore the existence of small
 faces in $\mathcal{T}_h$ does not
 affect the performance of the method.  It is only the relative sizes of
 the edges on each face that matter.
 }
\end{remark}
\par
 Note that the estimates in Lemma~\ref{lem:Xih} are also valid
 for ${\Omega}\subset \mathbb{R}^3$.
\subsection{Error Estimates in the Energy Norm}\label{subsec:3DEnergyError}
 The abstract error estimate
\begin{equation}\label{eq:3DAbstractEnergyError}
 \|u-u_h\|_h\lesssim \|u-I_{k,h} u\|_h+\|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u\|_h
  +\sqrt{\beta_h}\Big(|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u|_{h,1}+
  \sup_{w\in\cQ^k_h}\frac{(f,w-\Xi_hw)}{|w|_{H^1({\Omega})}}\Big)
\end{equation}
 is obtained by the same arguments as in Section~\ref{subsec:AbstractEnergyError},
 where $|\cdot|_{h,1}$ is defined  in \eqref{eq:PiecewiseHOneNOrm}.
\par
 We will derive concrete error estimates under the assumption that
 $u$ belongs to $H^{\ell+1}({\Omega})$ for $1\leq\ell\leq k$.
  Since the estimate
\begin{equation}\label{eq:3DEasy}
  |u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u|_{h,1}+ \sup_{w\in\cQ^k_h}\frac{(f,w-\Xi_hw)}{|w|_{H^1({\Omega})}}
  \lesssim h^\ell|u|_{H^{\ell+1}(\O)}
\end{equation}
  remains the same,
  we only need to estimate $\|u-I_{k,h} u\|_h$ and $\|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u\|_h$.
\par
 It follows from \eqref{eq:G1},  \eqref{eq:3DPODStability1},
 \eqref{eq:ConsistentID} and
 \eqref{eq:3DSD} that
\begin{align}\label{eq:3DInterpolationError1}
  \|u-I_{k,h} u\|_h^2&\lesssim \sum_{D\in\cT_h}|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)|_{H^1(D)}^2
  +\sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD} h_F^{-2}\|u-I_{k,D} u\|_{L_2(F)}^2\notag\\
  &\hspace{60pt}+\sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}}\SumFh_F^{-2}\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-I_{k,D} u)\|_{L_2(F)}^2\notag\\
    &\hspace{90pt}+\sum_{D\in\cT_h} \hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD} \|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} (u-I_{k,D} u)\|_{L_\infty(\partial F)}^2\\
        &\lesssim \sum_{D\in\cT_h} |u-I_{k,D} u|_{H^1(D)}^2
        +\sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD} h_F^{-2}\|u-I_{k,F} u\|_{L_2(F)}^2\notag\\
        &\hspace{60pt}+\sum_{D\in\cT_h} \hspace{1pt}h_{{\scriptscriptstyle D}}\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} (u-I_{k,D} u)\|_{L_\infty(D)}^2,\notag
\end{align}
 and we have, by \eqref{eq:HalfTrace} and \eqref{eq:HalfIDLTwoError},
\begin{equation}\label{eq:3DInterpolationError2}
  \sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD} h_F^{-2}\|u-I_{k,F} u\|_{L_2(F)}^2
  \lesssim \sum_{D\in\cT_h} \hspace{1pt}h_{{\scriptscriptstyle D}} \sum_{F\in\FD} h_F^{2\ell-1}|u|_{H^{\ell+(1/2)}(F)}^2
   \lesssim h^{2\ell}|u|_{H^{\ell+1}(\O)}^2.
\end{equation}
 Moreover the estimates \eqref{eq:G1}, \eqref{eq:DiscreteEstimate3} and
 \eqref{eq:3DPODStability1} imply
\begin{align}\label{eq:3DInterpolationError3}
 &\sum_{D\in\cT_h} \hspace{1pt}h_{{\scriptscriptstyle D}}\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} (u-I_{k,D} u)\|_{L_\infty(D)}^2 \notag\\
 &\hspace{60pt}
  \lesssim   \sum_{D\in\cT_h} \hspace{1pt}h_{{\scriptscriptstyle D}}^{-2}\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} (u-I_{k,D} u)\|_{L_2(D)}^2
  +\sum_{D\in\cT_h}|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} (u-I_{k,D} u)|_{H^1(D)}^2\\
  &\hspace{60pt}\lesssim
      \sum_{D\in\cT_h} \hspace{1pt}h_{{\scriptscriptstyle D}}^{-2}\big(\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u-u)\|_{L_2(D)}^2
      +\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D} u)\|_{L_2(D)}^2\big)\notag\\
      &\hspace{100pt}\notag+\sum_{D\in\cT_h}|u-I_{k,D} u|_{H^1(D)}^2.\notag
\end{align}
\par
 Combining \eqref{eq:3DLTwoPODPDZ}, \eqref{eq:3DLTwoPODID},
 \eqref{eq:3DHOnePODID} and
  \eqref{eq:3DInterpolationError1}--\eqref{eq:3DInterpolationError3},
 we obtain
\begin{equation}\label{eq:3DInterpolationError}
        \|u-I_{k,h} u\|_h^2\lesssim  h^{2\ell}|u|_{H^{\ell+1}(\O)}^2,
\end{equation}
 which is the analog of \eqref{eq:InterpolationError}.
\par
 From \eqref{eq:Sobolev}, \eqref{eq:G1},
 \eqref{eq:3DLTwoPODPDZ}--\eqref{eq:3DHTwoPODPDZ}
  and \eqref{eq:3DSD},  we find
\begin{align}\label{eq:3DProjectionError}
  \|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u\|_h^2&\lesssim
  \sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}}\SumFh_F^{-2}\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u\|_{L_2(F)}^2
   +\sum_{D\in\cT_h} \hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD} \|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u\|_{L_\infty(\partial F)}^2\notag\\
   &\lesssim \sum_{D\in\cT_h} \hspace{1pt}h_{{\scriptscriptstyle D}}\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u\|_{L_\infty(D)}^2\\
   &\lesssim    \sum_{D\in\cT_h}\big( \hspace{1pt}h_{{\scriptscriptstyle D}}^{-2}\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u\|_{L_2(D)}^2
      +|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u|_{H^1(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}^2|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u|_{H^2(D)}^2\big)\notag\\
      &\lesssim h^{2\ell}|u|_{H^{\ell+1}(\O)}^2,\notag
 \end{align}
 which is the analog of \eqref{eq:ProjectionError}.
\par
 The estimates \eqref{eq:3DAbstractEnergyError}, \eqref{eq:3DEasy},
 \eqref{eq:3DInterpolationError} and \eqref{eq:3DProjectionError}
 lead to the following analog of Theorem~\ref{thm:ConcreteEnergyError}.
\begin{theorem}\label{thm:3DConcreteEnergyError}
   Assuming the solution $u$ of \eqref{eq:Poisson} belongs to
 $H^{\ell+1}({\Omega})$ for $\ell$ between $1$ and $k$, we have
\begin{equation}\label{eq:3DConcreteEnergyError}
  \|u-u_h\|_h\leq C \sqrt{\beta_h}h^\ell|u|_{H^{\ell+1}(\O)}.
\end{equation}
 where $\beta_h$ is defined in \eqref{eq:lambdah} and
 the positive constant $C$ depends only on $\rho$, $N$ and $k$.
\end{theorem}
\par
 The following analog of
 Theorem~\ref{thm:ComputableEnergyError} on the computable approximate
 solutions $\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}  u_h$ and $\Pi_{k,h}^0 u_h$
 is obtained by the same arguments.
 \begin{theorem}\label{thm:3DComputableEnergyError}
    Assuming the solution $u$ of \eqref{eq:Poisson} belongs to
    $H^{\ell+1}({\Omega})$
   for $\ell$ between $1$ and $k$, there exists a positive constant $C$,
   depending only on
   $\rho$, $N$ and $k$, such that
 \begin{equation}\label{eq:3DComputableEnergyError}
  |u-u_h|_{H^1({\Omega})}+\sqrt{\beta_h}\big[|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h|_{h,1}
  +|u-\Pi_{k,h}^0 u|_{h,1}\big]\leq
  C \beta_h h^\ell |u|_{H^{\ell+1}({\Omega})},
\end{equation}
 where $\beta_h$ is defined in \eqref{eq:lambdah}.
\end{theorem}
\subsection{Error Estimates in the $L_2$ Norm}\label{subsec:3DL2Error}
 We begin with an analog of Lemma~\ref{lem:SDEstimate}.
\begin{lemma}\label{lem:3DSDEstimate}
 We have
\begin{equation*
\sum_{D\in\cT_h} S^D(\zeta-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}I_{k,h} \zeta,\zeta-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}I_{k,h} \zeta)\lesssim
  h^2|\zeta|_{H^2({\Omega})}^2\qquad\forall
  \zeta\in H^2({\Omega})\cap H^1_0({\Omega}).
\end{equation*}
\end{lemma}
\begin{proof} It follows from \eqref{eq:Sobolev},
\eqref{eq:3DLTwoPODID}--\eqref{eq:3DHTwoPODID} and
\eqref{eq:3DSD} that
\begin{align*}
 &\sum_{D\in\cT_h} S^D(\zeta-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} I_{k,h}\zeta,\zeta-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} I_{k,h}\zeta)\\
 &\hspace{30pt}\lesssim
 \sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD}\big(h_F^{-2}\|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta\|_{L_2(F)}^2
   +\|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta\|_{L_\infty(\partial F)}^2\big)\\
   &\hspace{30pt}\lesssim
   \sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}}\|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta\|_{L_\infty(\partial D)}^2\\
   &\hspace{30pt}\lesssim
   \sum_{D\in\cT_h} \big(\hspace{1pt}h_{{\scriptscriptstyle D}}^{-2}\|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta\|_{L_2(D)}^2+
   |\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta|_{H^1(D)}^2+\hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D}\zeta|_{H^2(D)}^2\big)\\
   &\hspace{30pt}\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^2|\zeta|_{H^2({\Omega})}^2.
\end{align*}
\end{proof}
\goodbreak
\par
  The same arguments as in the proof of Lemma~\ref{lem:POSD} lead
   to the following result.
\begin{lemma}\label{lem:3DPOSD}
  We have
\begin{equation*}
  \sum_{D\in\cT_h} S^D(u_h-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h,u_h-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h)\lesssim \beta_h^2 h^{2\ell}|u|_{H^{\ell+1}(\O)}^2.
\end{equation*}
\end{lemma}
\par
 With Lemma~\ref{lem:3DSDEstimate} and Lemma~\ref{lem:3DPOSD} in hand,
 we obtain
 the following analog of Theorem~\ref{thm:uhLTwoError}
  and Theorem~\ref{thm:ComputableLTwoErrors}
 by identical arguments.

 \begin{theorem}\label{thm:3DLTwoErrors}
 Assuming $u\in H^{\ell+1}({\Omega})$ for some $\ell$ between $1$ and $k$,
 there exists a positive constant $C$, depending only on $N$, $k$ and $\rho$,
  such that
\begin{align*}
 \|u-u_h\|_{L_2({\Omega})}+\|u-\Pi_{k,h}^0 u_h\|_{L_2({\Omega})}+ \|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h\|_{L_2({\Omega})}\leq C\beta_h
 h^{\ell+1}|u|_{H^{\ell+1}({\Omega})},
\end{align*}
 where $\beta_h$ is defined in \eqref{eq:lambdah}.
\end{theorem}
\subsection{Error Estimate in the $L_\infty$ Norm}\label{subsec:3DInftyError}
 We will derive $L_\infty$ error estimates under the additional
 assumption that $\mathcal{T}_h$ is quasi-uniform (cf. \eqref{eq:gamma}).
 We begin with an analog of Theorem~\ref{thm:MaxNormBdryEst1}.
\begin{theorem}\label{thm:uhInftyBdryEst}
      Assuming $\mathcal{T}_h$ is quasi-uniform and the solution $u$ of \eqref{eq:Poisson} belongs to
 $H^{\ell+1}({\Omega})$ for $\ell$ between $1$ and $k$, we have
\begin{equation}\label{eq:InftyBdryEst}
  \max_{e\in\mathcal{E}_h}\|u-u_h\|_{L_\infty(e)}\leq C\beta_h h^{\ell}|u|_{H^{\ell+1}(\O)},
\end{equation}
 where $\beta_h$ is defined in \eqref{eq:lambdah} and
  the positive constant $C$ only depends on $\rho$, $N$, $\gamma$ and $k$.
\end{theorem}
\begin{proof}  It follows from Remark~\ref{rem:NormEquivalence}
 that
\begin{equation*}
  \sum_{D\in\cT_h}\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD}\|(I_{k,D} u-u_h)-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} (I_{k,D} u-u_h)\|_{L_\infty(\partial F)}^2
   \lesssim \|I_{k,h} u-u_h\|_h^2
\end{equation*}
 and hence, for any $D\in\mathcal{T}_h$ and $F\in\mathcal{F}_\ssD$,
\begin{equation}\label{eq:3DBdryEst1}
  \|(I_{k,D} u-u_h)-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} (I_{k,D} u-u_h)\|_{L_\infty(\partial F)}^2\lesssim
    \beta_h\, h^{2\ell-1}|u|_{H^{\ell+1}(\O)}^2
\end{equation}
 by  \eqref{eq:3DInterpolationError}, Theorem~\ref{thm:3DConcreteEnergyError}
 and the quasi-uniformity of $\mathcal{T}_h$,
\par
 For any $D\in\mathcal{T}_h$ and $F\in\mathcal{F}_\ssD$,
   we have, by \eqref{eq:G1}, \eqref{eq:DiscreteEstimate3},
   \eqref{eq:3DPODStability1},
   \eqref{eq:3DLTwoPODID}, \eqref{eq:3DHOnePODID},
   Theorem~\ref{thm:3DLTwoErrors}
   and the quasi-uniformity of $\mathcal{T}_h$,
\begin{align}\label{eq:3DBdryEst2}
  &\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(I_{k,D} u-u_h)\|_{L_\infty(F)}^2\notag\\
    &\hspace{40pt}\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-3}\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(I_{k,D} u-u_h)\|_{L_2(D)}^2+
      \hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(I_{k,D} u-u_h)|_{H^1(D)}^2\notag\\
    &\hspace{40pt}\lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{-3}\big(\|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}I_{k,D} u-u\|_{L_2(D)}^2
    +\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h\|_{L_2(D)}^2\big)
   
    +\hspace{1pt}h_{{\scriptscriptstyle D}}^{-1}|I_{k,D} u-u_h|_{H^1(D)}^2\\
    &\hspace{40pt}\lesssim \beta_h^2 h^{2\ell-1}|u|_{H^{\ell+1}(\O)}^2.\notag
\end{align}
\par
 Finally we have, for any $D\in\mathcal{T}_h$ and $F\in\mathcal{F}_\ssD$,
\begin{equation}\label{eq:3DBdryEst3}
  \|u-I_{k,D} u\|_{L_\infty(F)}^2=\|u-I_{k,F} u\|_{L_\infty(F)}^2\lesssim
  h_F^{2\ell-1}|u|_{H^{\ell+\frac12}(F)}^2
   \lesssim \hspace{1pt}h_{{\scriptscriptstyle D}}^{2\ell-1}|u|_{H^{\ell+1}(D)}^2
\end{equation}
 by \eqref{eq:HalfTrace}, \eqref{eq:1DInterpolationErrorHalf} and
 \eqref{eq:ConsistentID}.
\par
 %
  The estimate \eqref{eq:InftyBdryEst} then follows from
  \eqref{eq:3DBdryEst1}--\eqref{eq:3DBdryEst3} and the triangle inequality.
\end{proof}
\par
 We also have estimates for
 the computable approximate solutions $\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h$ and
 $\Pi_{k,h}^0 u_h$.
\begin{theorem}\label{thm:3DMaxNormBdryEst}
     Assuming $\mathcal{T}_h$ is quasi-uniform and the solution $u$ of \eqref{eq:Poisson} belongs to
 $H^{\ell+1}({\Omega})$ for $\ell$ between $1$ and $k$, there exists a positive constant
 $C$, depending only on $\rho$, $N$, $\gamma$ and $k$, such that
\begin{equation}\label{eq:3DLInftyError}
   \|u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h\|_{L_\infty({\Omega})}+\|u-\Pi_{k,h}^0 u_h\|_{L_\infty({\Omega})}\leq C
   \beta_h h^{\ell-(1/2)}|u|_{H^{\ell+1}(\O)},
\end{equation}
 where $\beta_h$ is defined in \eqref{eq:lambdah}.
\end{theorem}
\begin{proof}
   For any $D\in\mathcal{T}_h$, we have,  by \eqref{eq:Sobolev},
 \eqref{eq:DiscreteEstimate1} and \eqref{eq:3DPODStability1},
\begin{align*}
 &\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h\|_{{L_\infty(D)}}\\
 &\hspace{40pt}\lesssim
 \hspace{1pt}h_{{\scriptscriptstyle D}}^{-\frac32}\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h\|_{L_2(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}^{-\frac12}|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h|_{H^1(D)}
 +\hspace{1pt}h_{{\scriptscriptstyle D}}^\frac12|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h|_{H^2(D)}\\
 &\hspace{40pt}\lesssim
 \hspace{1pt}h_{{\scriptscriptstyle D}}^{-\frac32}\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h\|_{L_2(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}^{-\frac12}|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h|_{H^1(D)}
 +\hspace{1pt}h_{{\scriptscriptstyle D}}^\frac12|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u|_{H^2(D)}\\
    &\hspace{70pt}+\hspace{1pt}h_{{\scriptscriptstyle D}}^\frac12|\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt}(u-u_h)|_{H^2(D)}\\
     &\hspace{40pt}\lesssim
 \hspace{1pt}h_{{\scriptscriptstyle D}}^{-\frac32}\|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h\|_{L_2(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}^{-\frac12}|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u_h|_{H^1(D)}
 +\hspace{1pt}h_{{\scriptscriptstyle D}}^\frac12|u-\Pi_{k,D}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} u|_{H^2(D)}\\
 &\hspace{60pt}+\hspace{1pt}h_{{\scriptscriptstyle D}}^{-\frac12}|u-u_h|_{H^1(D)},
\end{align*}
 which together with \eqref{eq:3DHTwoPODPDZ},
 Theorem~\ref{thm:3DComputableEnergyError}, Theorem~\ref{thm:3DLTwoErrors}
 and the quasi-uniformity of $\mathcal{T}_h$ implies the estimate for $u-\Pi_{k,h}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}} u_h$.
\par
 Similarly we have, by \eqref{eq:Sobolev},
 \eqref{eq:DiscreteEstimate1} and
 \eqref{eq:3DPDZHOne},
\begin{align*}
 &\|u-\Pi_{k,D}^0\hspace{1pt} u_h\|_{{L_\infty(D)}}\\
  &\hspace{40pt}\lesssim
 \hspace{1pt}h_{{\scriptscriptstyle D}}^{-\frac32}\|u-\Pi_{k,D}^0\hspace{1pt} u_h\|_{L_2(D)}+\hspace{1pt}h_{{\scriptscriptstyle D}}^{-\frac12}|u-\Pi_{k,D}^0\hspace{1pt} u_h|_{H^1(D)}
 +\hspace{1pt}h_{{\scriptscriptstyle D}}^\frac12|u-\Pi_{k,D}^0\hspace{1pt} u|_{H^2(D)}\\
 &\hspace{70pt}+\hspace{1pt}h_{{\scriptscriptstyle D}}^{-\frac12}|u-u_h|_{H^1(D)},
\end{align*}
 which together with \eqref{eq:3DHTwoPODPDZ},
 Theorem~\ref{thm:3DComputableEnergyError},
 Theorem~\ref{thm:3DLTwoErrors}
 and the quasi-uniformity of $\mathcal{T}_h$ implies the estimate for $u-\Pi_{k,h}^0 u_h$.
\end{proof}
\section{Concluding Remarks}\label{sec:Conclusions}
 We have developed error estimates for virtual element methods for the model
 Poisson problem in two and three dimensions that provide justifications for
  existing numerical results
 for polygonal (or polyhedral) meshes with small edges (or faces).
\par
 For the two dimensional problem, the convergence of the virtual element
 method based on the stabilizing bilinear form $S^D_2(\cdot,\cdot)$ is optimal under the
 shape regularity assumptions in Section~\ref{subsec:GlobalShapeRegularity}.
 Under the additional assumption that the edges of any subdomain
 in a polygonal mesh are comparable to one another,
 convergence of the virtual element method based on
 the stabilizing bilinear form $S^D_1(\cdot,\cdot)$ is also optimal.
\par
 For the three dimensional problem, the convergence of the virtual element method
  is optimal if, in addition to the assumptions in Section~\ref{subsec:Shape3D}, we
  also assume that the edges of any face in the polyhedral mesh are comparable to one another.
\par
  The results in this paper can be extended to virtual element methods with the
  stabilizing bilinear form
\begin{equation*}
  S^D(w,v)=(\Pi_{k-2,D}^0 w,\Pi_{k-2,D}^0 v)+\sum_{p\in\mathcal{N}_{\p D}}w(p)v(p)
\end{equation*}
 in two dimensions, and  the  stabilizing bilinear form
\begin{align*}
    S^D(w,v)&=(\Pi_{k-2,D}^0 w,\Pi_{k-2,D}^0 v)+\hspace{1pt}h_{{\scriptscriptstyle D}}\sum_{F\in\FD}\Big(h_F^{-2}(\Pi_{k-2,F}^0 w,\Pi_{k-2,F}^0 v)_{L_2(F)}+
   \sum_{p\in\mathcal{N}_{\p F}}w(p)v(p)\Big).\label{eq:3DSD}
\end{align*}
 in three dimensions.  The stability for these virtual element methods is automatic and the
 error analysis also does not pose any new difficulties.
\par
 The results in this paper can also be extended to virtual element methods ($k\geq2$)
 where the inner product \eqref{eq:InnerProduct} is replaced by the inner product
\begin{equation*}
 (\!(\zeta,\eta)\!)=(\nabla \zeta,\nabla \eta)
     +\Big(\int_{D}\zeta\,dx\Big)\Big(\int_{D} \eta\,dx\Big).
\end{equation*}
\par
 We note that error estimates for the Poisson problem on
 general polygonal or polyhedral domains can also
 be obtained by the techniques developed in this paper.
\par
 Finally it would be interesting to construct a three dimensional
 analog of the stabilizing bilinear form
 $S^D_2(\cdot,\cdot)$ defined in Section~\ref{subsec:DiscreteProblem}
  so that the convergence of the virtual element methods is optimal
 for polyhedral meshes with arbitrarily small faces and edges, and
 $L_\infty$ error estimates can be established without assuming the meshes are quasi-uniform.
  We conjecture that such a bilinear form can be defined by
\begin{align*}
  S^D(v,w)&=\SumFh_F(\nabla_{\!\!F}\Pi_{k,F}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v,\nabla_{\!\!F}\Pi_{k,F}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} w)_{L_2(F)}\\
   &\hspace{40pt}+\SumFh_F\sum_{e\in\cE_{\ssF}}
   h_e\big(\partial (v-\Pi_{k,F}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} v)/\partial s,\partial (w-\Pi_{k,F}^{\raise 1pt\hbox{$\scriptscriptstyle\nabla$}}\hspace{1pt} w)/\partial s\big)_{L_2(e)}.
\end{align*}

2017-10-03T02:12:36.000
1710
1710.00442
en
https://arxiv.org/abs/1710.00442
We consider a model Poisson problem in $\R^d$ ($d=2,3$) and establish error estimates for virtual element methods on polygonal or polyhedral meshes that can contain small edges ($d=2$) or small faces ($d=3$).
Numerical Analysis (math.NA)
Virtual Element Methods on Meshes with Small Edges or Faces
Qwen/Qwen-72B
1. YES
2. YES