High-order polygonal discontinuous Petrov-Galerkin (PolyDPG) methods using ultraweak formulations

TL;DR

This paper introduces a novel high-order polygonal finite element method called PolyDPG, utilizing ultraweak formulations within the discontinuous Petrov-Galerkin framework, enabling stable, conforming polygonal discretizations with adaptive capabilities.

Contribution

It is the first to apply ultraweak formulations to high-order polygonal DPG methods, eliminating the need for stabilization and providing convergence proofs and adaptive strategies.

Findings

PolyDPG achieves optimal convergence rates.

Method handles distorted and concave polygonal meshes.

Includes an open-source software implementation.

Abstract

This work represents the first endeavor in using ultraweak formulations to implement high-order polygonal finite element methods via the discontinuous Petrov-Galerkin (DPG) methodology. Ultraweak variational formulations are nonstandard in that all the weight of the derivatives lies in the test space, while most of the trial space can be chosen as copies of $L^2$-discretizations that have no need to be continuous across adjacent elements. Additionally, the test spaces are broken along the mesh interfaces. This allows one to construct conforming polygonal finite element methods, termed here as PolyDPG methods, by defining most spaces by restriction of a bounding triangle or box to the polygonal element. The only variables that require nontrivial compatibility across elements are the so-called interface or skeleton variables, which can be defined directly on the element boundaries. Unlike other high-order polygonal methods, PolyDPG methods do not require ad hoc stabilization terms thanks to the crafted stability of the DPG methodology. A proof of convergence of the form $h^p$ is provided and corroborated through several illustrative numerical examples. These include polygonal meshes with $n$-sided convex elements and with highly distorted concave elements, as well as the modeling of discontinuous material properties along an arbitrary interface that cuts a uniform grid. Since PolyDPG methods have a natural a posteriori error estimator a polygonal adaptive strategy is developed and compared to standard adaptivity schemes based on constrained hanging nodes. This work is also accompanied by an open-source $\texttt{PolyDPG}$ software supporting polygonal and conventional elements.

Equations

$$-\mathrm{div}(k\nabla u)=r,\iff \{\begin{array}{rl}\mathrm{div}\mathbf{q}& =r,\\ \frac{1}{k}\mathbf{q}+\nabla u& =0.\end{array}$$

$$\begin{array}{c}{b}^{\mathcal{P}}(u,v)={\ell }^{\mathcal{P}}(v)\forall v\in {\mathcal{V}}^{\mathcal{P}}={\mathcal{U}}^{\mathcal{P}}=H_0^1(\Omega ),\\ b^{\mathcal{P}}(u,v)=(k\nabla u,\nabla v{)}_{\Omega },{\ell }^{\mathcal{P}}(v)=(r,v{)}_{\Omega },\end{array}$$

$$\begin{gathered}b_{0}(\mathfrak{u}_{0},\mathfrak{v}_{0})=\ell(\mathfrak{v}_{0})\,\qquad\forall(v,\boldsymbol{\tau})=\mathfrak{v}_{0}\in\mathscr{V}_{0}=H_{0}^{1}(\Omega)\times\boldsymbol{H}(\operatorname{div},\Omega)\,,\\ b_{0}\big{(}(u,\boldsymbol{q}),(v,\boldsymbol{\tau})\big{)}=-(\boldsymbol{q},\nabla v)_{\Omega}+(\textstyle{\frac{1}{k}}\boldsymbol{q},\boldsymbol{\tau})_{\Omega}-(u,\operatorname{div}\boldsymbol{\tau})_{\Omega},\qquad\ell\big{(}(v,\boldsymbol{\tau})\big{)}=(r,v)_{\Omega}\,,\end{gathered}$$

$$H^1(\mathcal{T})$$

$$\mathbf{H}(\mathrm{div},\mathcal{T})$$

$$(u,v{)}_{\mathcal{T}}=\sum _{K\in \mathcal{T}}(u{|}_K,v{|}_K{)}_K.$$

$$H_0^{{}^1}(\partial \mathcal{T})$$

$$H^{{}^1}(\partial \mathcal{T})$$

$$⟨\hat{u},\hat{v}{⟩}_{\partial \mathcal{T}}=\sum _{K\in \mathcal{T}}⟨(\hat{u}{)}_K,(\hat{v}{)}_K{⟩}_{\partial K},$$

$$\begin{array}{c}({\mathfrak{u}}_0,\hat{\mathfrak{u}})=\mathfrak{u}\in \mathcal{U}={\mathcal{U}}_0\times \hat{\mathcal{U}},\\ (u,\mathbf{q})={\mathfrak{u}}_0\in {\mathcal{U}}_0=L^2(\Omega )\times {\mathbf{L}}^2(\Omega ),(\hat{u},{\hat{q}}_{\hat{\mathbf{n}}})=\hat{\mathfrak{u}}\in \hat{\mathcal{U}}=H_0^{{}^1}(\partial \mathcal{T})\times H^{{}^1}(\partial \mathcal{T}),\end{array}$$

$$b(\mathfrak{u},\mathfrak{v})=\ell (\mathfrak{v})\forall (v,\mathbf{\tau })=\mathfrak{v}\in \mathcal{V}=H^1(\mathcal{T})\times \mathbf{H}(\mathrm{div},\mathcal{T}),$$

$$\displaystyle b\big{(}(\mathfrak{u}_{0},\hat{\mathfrak{u}}),\mathfrak{v}\big{)}=b_{0}(\mathfrak{u}_{0},\mathfrak{v})+\hat{b}(\hat{\mathfrak{u}},\mathfrak{v})\,,\qquad\ell\big{(}(v,\boldsymbol{\tau})\big{)}=(r,v)_{\mathcal{T}}\,,$$

$$\displaystyle b_{0}\big{(}(u,\boldsymbol{q}),(v,\boldsymbol{\tau})\big{)}=-(\boldsymbol{q},\nabla v)_{\mathcal{T}}+(\textstyle{\frac{1}{k}}\boldsymbol{q},\boldsymbol{\tau})_{\mathcal{T}}-(u,\operatorname{div}\boldsymbol{\tau})_{\mathcal{T}}\,,$$

$$\displaystyle\hat{b}\big{(}(\hat{u},\hat{q}_{\hat{\mathbf{n}}}),(v,\boldsymbol{\tau})\big{)}=\langle\hat{q}_{\hat{\mathbf{n}}},v_{\partial\mathcal{T}}\rangle_{\partial\mathcal{T}}+\langle\hat{u},\boldsymbol{\tau}_{\partial\mathcal{T}}\rangle_{\partial\mathcal{T}}\,,$$

$${\mathsf{B}}^{\mathcal{P}}{\mathsf{u}}_h=\text{l},$$

$$\underset{\delta {\mathfrak{u}}_h\in {\mathcal{U}}_h\setminus \{0\}}{\inf }\underset{{\mathfrak{v}}_h\in {\mathcal{V}}_h\setminus \{0\}}{\sup }\frac{b(\delta {\mathfrak{u}}_h,{\mathfrak{v}}_h)}{\Vert \delta {\mathfrak{u}}_h{\Vert }_{\mathcal{U}}\Vert {\mathfrak{v}}_h{\Vert }_{\mathcal{V}}}={\gamma }_h>0,$$

$${\mathfrak{u}}_h^{\mathrm{opt}}=\underset{\delta {\mathfrak{u}}_h\in {\mathcal{U}}_h}{\mathrm{argmin}}\Vert \mathcal{B}\delta {\mathfrak{u}}_h-\ell {\Vert }_{{\mathcal{V}}^{\prime }},\iff b({\mathfrak{u}}_h^{\mathrm{opt}},{\mathfrak{v}}^{\mathrm{opt}})=\ell ({\mathfrak{v}}^{\mathrm{opt}})\forall {\mathfrak{v}}^{\mathrm{opt}}\in {\mathcal{V}}^{\mathrm{opt}}={\mathcal{R}}_{\mathcal{V}}^{-1}\mathcal{B}{\mathcal{U}}_h,$$

$${\mathsf{B}}^{\mathrm{n}\text{-}\mathrm{opt}}{\mathsf{u}}_h={\mathsf{B}}^{\mathsf{T}}{\mathsf{G}}^{-1}\mathsf{B}{\mathsf{u}}_h={\mathsf{B}}^{\mathsf{T}}{\mathsf{G}}^{-1}\text{l}=\text{l},$$

$$\begin{gathered}{}\\[-4.0pt] \text{\footnotesize$\{\mathfrak{v}_{i}\}_{i=1}^{M}$}\left\{\vphantom{\begin{pmatrix}|\\ |\\ |\end{pmatrix}}\right.\begin{bmatrix}\makebox[0.0pt][l]{$\smash{\overbrace{\phantom{\begin{matrix}\qquad\quad|\qquad\quad\end{matrix}}}^{\{(\mathfrak{u}_{0})_{j}\}_{j=1}^{N_{0}}}}$}\qquad\quad|\qquad\quad&\makebox[0.0pt][l]{$\smash{\overbrace{\phantom{\begin{matrix}\quad\qquad|\qquad\quad\end{matrix}}}^{\{\hat{\mathfrak{u}}_{j}\}_{j=1}^{\hat{N}}}}$}\quad\qquad|\qquad\quad\\ \quad\mathsf{B}_{0}&\quad\hat{\mathsf{B}}\\ \qquad\quad|\qquad\quad&\qquad\quad|\qquad\quad\end{bmatrix}=\mathsf{B}=\begin{bmatrix}\mathsf{B}_{uv}&\mathsf{B}_{\boldsymbol{q}v}&\mathsf{B}_{\hat{u}v}&\mathsf{B}_{\hat{q}_{\hat{\mathbf{n}}}v}\\ \mathsf{B}_{u\boldsymbol{\tau}}&\mathsf{B}_{\boldsymbol{q}\boldsymbol{\tau}}&\mathsf{B}_{\hat{u}\boldsymbol{\tau}}&\mathsf{B}_{\hat{q}_{\hat{\mathbf{n}}}\boldsymbol{\tau}}\end{bmatrix}\end{gathered}$$

$$\Vert \mathcal{B}{\mathfrak{u}}_h-\ell {\Vert }_{{\mathcal{V}}^{\prime }}^2\approx \Vert {\mathcal{R}}_{{\mathcal{V}}_r}^{-1}(\mathcal{B}{\mathfrak{u}}_h-\ell ){\Vert }_{\mathcal{V}}^2=(\mathsf{B}{\mathsf{u}}_h-\text{l}{)}^{\mathsf{T}}{\mathsf{G}}^{-1}(\mathsf{B}{\mathsf{u}}_h-\text{l}),$$

$$\Vert (v,\mathbf{\tau }){\Vert }_{\mathcal{V}}^2=\Vert v{\Vert }_{H^1(\mathcal{T})}^2+\Vert \mathbf{\tau }{\Vert }_{H(\mathrm{div},\mathcal{T})}^2=(v,v{)}_{\mathcal{T}}+(\nabla v,\nabla v{)}_{\mathcal{T}}+(\mathbf{\tau },\mathbf{\tau }{)}_{\mathcal{T}}+(\mathrm{div}\mathbf{\tau },\mathrm{div}\mathbf{\tau }{)}_{\mathcal{T}}.$$

$$\|(v,\boldsymbol{\tau})\|_{\mathscr{V}}^{2}=\|\textstyle{\frac{1}{k}}\boldsymbol{\tau}-\nabla v\|_{\boldsymbol{L}^{2}(\mathcal{T})}^{2}+\|-\operatorname{div}\boldsymbol{\tau}\|_{L^{2}(\mathcal{T})}^{2}+\varepsilon^{2}\big{(}\|v\|_{L^{2}(\mathcal{T})}^{2}+\|\boldsymbol{\tau}\|_{\boldsymbol{L}^{2}(\mathcal{T})}^{2}\big{)}\,,$$

$$\begin{gathered}\lx@xy@svg{\hbox{\raise 0.0pt\hbox{\kern 19.04604pt\hbox{\ignorespaces\ignorespaces\ignorespaces\hbox{\vtop{\kern 0.0pt\offinterlineskip\halign{\entry@#!@&&\entry@@#!@\cr&&\\&&\crcr}}}\ignorespaces{\hbox{\kern-19.04604pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{H^{1}(T_{K})\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 22.51686pt\raise 5.43056pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern 0.0pt\raise-2.43056pt\hbox{$\scriptstyle{\operatorname{curl}\quad}$}}}\kern 3.0pt}}}}}}\ignorespaces{\hbox{\kern 43.04604pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 43.04604pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\boldsymbol{H}(\operatorname{div},T_{K})\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 99.3932pt\raise 5.39166pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern 0.0pt\raise-2.39166pt\hbox{$\scriptstyle{\nabla\cdot}$}}}\kern 3.0pt}}}}}}\ignorespaces{\hbox{\kern 121.62473pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 121.62473pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\,\,L^{2}(T_{K})\,\,}$}}}}}}}{\hbox{\kern-17.8951pt\raise-23.13776pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\mathcal{P}^{p}(T_{K})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces{\hbox{\kern 0.0pt\raise-11.56888pt\hbox{{}{}\hbox{\kern-3.8889pt\raise-2.5pt\hbox{\xyRotate@@{-3072}\kern 0.0pt\kern 0.0pt}}}}}\ignorespaces{}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 22.51686pt\raise-17.7072pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern 0.0pt\raise-2.43056pt\hbox{$\scriptstyle{\operatorname{curl}\quad}$}}}\kern 3.0pt}}}}}}\ignorespaces{\hbox{\kern 47.93588pt\raise-23.13776pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 47.93588pt\raise-23.13776pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\mathcal{R}\mathcal{T}^{p}(T_{K})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces{\hbox{\kern 69.71986pt\raise-11.56888pt\hbox{{}{}\hbox{\kern-3.8889pt\raise-2.5pt\hbox{\xyRotate@@{-3072}\kern 0.0pt\kern 0.0pt}}}}}\ignorespaces{}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 99.3932pt\raise-17.7461pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern 0.0pt\raise-2.39166pt\hbox{$\scriptstyle{\nabla\cdot}$}}}\kern 3.0pt}}}}}}\ignorespaces{\hbox{\kern 120.39369pt\raise-23.13776pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 120.39369pt\raise-23.13776pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\mathcal{P}^{p-1}(T_{K})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\,,}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces{\hbox{\kern 142.84431pt\raise-11.56888pt\hbox{{}{}\hbox{\kern-3.8889pt\raise-2.5pt\hbox{\xyRotate@@{-3072}\kern 0.0pt\kern 0.0pt}}}}}\ignorespaces{}\ignorespaces}}}}\ignorespaces\end{gathered}$$

$$\begin{gathered}\lx@xy@svg{\hbox{\raise 0.0pt\hbox{\kern 20.90419pt\hbox{\ignorespaces\ignorespaces\ignorespaces\hbox{\vtop{\kern 0.0pt\offinterlineskip\halign{\entry@#!@&&\entry@@#!@\cr&&\\&&\crcr}}}\ignorespaces{\hbox{\kern-19.38249pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{H^{1}(Q_{K})\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 23.56874pt\raise 5.43056pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern 0.0pt\raise-2.43056pt\hbox{$\scriptstyle{\operatorname{curl}\qquad\qquad}$}}}\kern 3.0pt}}}}}}\ignorespaces{\hbox{\kern 45.81337pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 45.81337pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\qquad\boldsymbol{H}(\operatorname{div},Q_{K})\qquad\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 127.95772pt\raise 5.39166pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern 0.0pt\raise-2.39166pt\hbox{$\scriptstyle{\qquad\nabla\cdot}$}}}\kern 3.0pt}}}}}}\ignorespaces{\hbox{\kern 170.98016pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 170.98016pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\,\,L^{2}(Q_{K})\,\,}$}}}}}}}{\hbox{\kern-20.90419pt\raise-23.13776pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\mathcal{Q}^{p,p}(Q_{K})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces{\hbox{\kern 0.0pt\raise-11.56888pt\hbox{{}{}\hbox{\kern-3.8889pt\raise-2.5pt\hbox{\xyRotate@@{-3072}\kern 0.0pt\kern 0.0pt}}}}}\ignorespaces{}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 23.56874pt\raise-17.7072pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern 0.0pt\raise-2.43056pt\hbox{$\scriptstyle{\operatorname{curl}\qquad\qquad}$}}}\kern 3.0pt}}}}}}\ignorespaces{\hbox{\kern 44.90419pt\raise-23.13776pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 44.90419pt\raise-23.13776pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\mathcal{Q}^{p,p-1}(Q_{K})\!\times\!\mathcal{Q}^{p-1,p}(Q_{K})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces{\hbox{\kern 92.82367pt\raise-11.56888pt\hbox{{}{}\hbox{\kern-3.8889pt\raise-2.5pt\hbox{\xyRotate@@{-3072}\kern 0.0pt\kern 0.0pt}}}}}\ignorespaces{}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 127.95772pt\raise-17.7461pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern 0.0pt\raise-2.39166pt\hbox{$\scriptstyle{\qquad\nabla\cdot}$}}}\kern 3.0pt}}}}}}\ignorespaces{\hbox{\kern 164.74315pt\raise-23.13776pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 164.74315pt\raise-23.13776pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\mathcal{Q}^{p-1,p-1}(Q_{K})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\,,}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces{\hbox{\kern 192.53621pt\raise-11.56888pt\hbox{{}{}\hbox{\kern-3.8889pt\raise-2.5pt\hbox{\xyRotate@@{-3072}\kern 0.0pt\kern 0.0pt}}}}}\ignorespaces{}\ignorespaces}}}}\ignorespaces\end{gathered}$$

$${\mathcal{P}}^{p-1}(\partial K)$$

$${\mathcal{P}}_C^p(\partial K)$$

$$\displaystyle\mathscr{U}_{h}=\big{\{}(u,\boldsymbol{q},\hat{u},\hat{q}_{\hat{\mathbf{n}}})\in\mathscr{U}\mid u|_{K}$$

$$\displaystyle\qquad\hat{u}_{K}\in\mathcal{P}^{p}_{C}(\partial K),\,(\hat{q}_{\hat{\mathbf{n}}})_{K}\in\mathcal{P}^{p-1}(\partial K),\,\forall K\in\mathcal{T}\big{\}}\,.$$

$$\mathscr{V}_{r}=\big{\{}(v,\boldsymbol{\tau})\mid v|_{K}\in\mathcal{P}^{p+\Delta p_{K}}(K),\,\boldsymbol{\tau}|_{K}\in\mathcal{R}\mathcal{T}^{p+\Delta p_{K}}(K),\,\forall K\in\mathcal{T}\big{\}}\,.$$

$$\dim ({\mathcal{U}}_h(K))$$

$$\dim ({\mathcal{V}}_r(K))$$

$$\mathrm{ov}({\mathcal{T}}_{\mathcal{K}})=\underset{x\in {\mathbb{R}}^2}{\sup }\mathrm{ov}(x)<\infty ,\mathrm{ov}(x)=|\{\mathcal{K}\in {\mathcal{T}}_{\mathcal{K}}\mid x\in \mathcal{K}\}|.$$

$$b({\mathfrak{u}}_h,{\mathfrak{v}}_h)=\ell ({\mathfrak{v}}_h),\forall {\mathfrak{v}}_h\in {\mathcal{V}}_h={\mathcal{R}}_{{\mathcal{V}}_r}^{-1}\mathcal{B}{\mathcal{U}}_h,$$

$$\Vert \mathfrak{u}-{\mathfrak{u}}_h{\Vert }_{\mathcal{U}}\le C{h}^{\min \{s,p\}}\Vert \mathfrak{u}{\Vert }_{{\mathcal{U}}^s},$$

$$\text{Relative error}=\frac{\Vert {\mathfrak{u}}_0-({\mathfrak{u}}_0{)}_h{\Vert }_{{\mathcal{U}}_0}}{\Vert {\mathfrak{u}}_0{\Vert }_{{\mathcal{U}}_0}},\text{with}\Vert (u,\mathbf{q}){\Vert }_{{\mathcal{U}}_0}^2=\Vert u{\Vert }_{L^2(\Omega )}^2+\Vert \mathbf{q}{\Vert }_{{\mathbf{L}}^2(\Omega )}^2=(u,u{)}_{\Omega }+(\mathbf{q},\mathbf{q}{)}_{\Omega },$$

Full text

Keywords: discontinuous Petrov-Galerkin (DPG) methodology, ultraweak formulations, polygonal finite element methods, adaptivity, distortion tolerance, high-order discretization

1 Introduction

Numerical solutions of boundary value problems with meshes of general polytopes were first proposed by Wachspress [88], who introduced rational barycentric coordinates that formed a finite element basis over convex polygons, leading to a conforming finite element method (FEM) with new types of elements. Over the last two decades, there has been a growing collection of numerical methods using general polytopes which extend well beyond the original ideas of Wachspress. Among the reasons for this group of methods to thrive is a handful of advantages that polytopes offer over traditionally shaped elements (simplices, hexahedra, etc.). These include: matching complex interfaces (see e.g. [73, 25]); greater flexibility to mesh complex geometries and their role as transition elements [82]; avoiding the limitations of parametric elements for highly distorted or ill-shaped elements (see e.g. [27, 67]); handling multiple hanging nodes in local $h$-refinements [83]; and allowing for greater deformations and less tendency to mesh-locking in incompressible media [28].

The features just mentioned give polytopal FEMs a wide range of applicability, especially where conventional methods do not fare well. In fact, they are useful for resolving problems involving the deformation of materials with heterogeneous microstructure [54], modeling complex materials like elastomers and biomaterials [28, 35], creating meshes where interface fitting is required [25], and modeling fractured media [12]. Promising results have also been obtained in crack propagation modeling [80, 68, 13, 15] and in topology optimization [84, 53, 3, 86], since polygonal meshes combine the ability to mesh complex geometries with a reasonable number of elements while reducing mesh-induced bias in particular directions (which occurs in structured meshes of triangles or quadrilaterals) [84, 68, 3].

Many methods still utilize different types of generalized barycentric coordinates (including some valid in nonconvex polytopes), which have proliferated since Wachspress originally introduced them, as well as other choices of shape functions (see e.g. [14]). These methods are usually $H^{1}$-conforming Galerkin FEMs [82], but there are some extensions to mixed methods (see e.g. [28]). They mostly allow very flexible refinement schemes while avoiding constrained approximations [83], but they are typically limited by first order $h$-convergence. Some families of high-order shape functions have been proposed, but only for convex polytopes (see e.g. [78, 56]). As the barycentric coordinates are in general rational polynomials, another challenge is the choice of the quadrature scheme used for integration [72, 29].

Mimetic finite difference (MFD) methods are based on another discretization technique which also supports polygonal elements. The technique consists of designing discrete differential operators such that fundamental vector calculus identities and physical laws can be reproduced in a discrete context [66, 18, 17]. Later, the ideas of MFDs led to the development of virtual element methods (VEMs) [8]. In VEMs, appropriate spaces are tailored for each polytopal element, such that their functions have continuous and piecewise polynomial traces over the boundaries. The integrals over the cells can be computed exactly (i.e. up to machine precision) with quadrature points only on the boundary [69]. The power of VEMs lies partly in eliminating the need of explicitly constructing the shape functions in the element, and yet resulting in a FEM-like variational setting [11]. They are also high-order methods [7], and recent work has resulted in the construction of $\boldsymbol{H}(\operatorname{div})$- and $\boldsymbol{H}(\operatorname{curl})$-conforming spaces [10]. VEMs have been used for different problems like linear elasticity, plate bending, and second-order elliptic problems [9, 19, 11]. But it must be noted that VEMs need a problem-dependent stability operator to guarantee their convergence [69], and the solution at interior points of the elements is not accessible directly, so it has to be approximated [11].

Another method is the polytopal interior penalty $hp$ discontinuous Galerkin (IPDG) method [20]. It is a nonconforming high-order method, which uses restrictions of standard FE spaces associated to a bounding box of each element. Due to its nonconformity, the method has a thorough but nonstandard equation-dependent error analysis, and like VEMs, it needs adding extra terms to ensure stability. Lastly, other recent methods include hybrid mimetic mixed methods [47, 46], PFEM-VEM [69], the weak Galerkin (WG) method [73, 74, 89], hybrid high-order (HHO) methods [45], and hybridizable discontinuous Galerkin (HDG) methods [31, 33]. More details on the historical development can be found in the thorough review [69].

The objective of this article is to present a completely new family of high-order methods termed polygonal discontinuous Petrov-Galerkin (PolyDPG) methods. They are based on so-called “broken” ultraweak variational formulations discretized using the discontinuous Petrov-Galerkin (DPG) methodology [39]. These formulations, despite being well-defined at the infinite-dimensional level, admit a very large degree of discontinuities in both the trial and test spaces, since their test spaces are broken (i.e. they may be discontinuous across element interfaces) and part of their trial spaces is in $L^{2}$. In fact, the only communication between elements happens through the so-called skeleton (or interface) variables that live on the element boundaries. These nonstandard formulations can be systematically discretized in a conforming fashion (i.e., with discrete trial and test spaces that are subspaces of the infinite-dimensional ones) and solved using the variationally versatile DPG methodology, which always produces a positive-definite finite element stiffness matrix. The DPG methodology is essentially crafted to produce stability by using optimal test functions and without resorting to additional stabilization terms. DPG methods have been successfully used for equations involving numerical stability issues [34, 43, 23, 75, 64], and applied to various physical problems such as wave propagation [90, 57, 40, 77], transmission problems [59, 52], electromagnetism [21], elasticity [62, 16, 50, 49], fluid flow [79, 22, 48, 63] and optical fibers via Schrödinger’s equation [41].

In this paper we consider 2D problems, where the element boundaries are merely line segments, so high-order discretization of the skeleton variables is straightforward. As we will show, this makes the broken ultraweak formulations an ideal framework for defining polygonal elements, and it results in the conforming FEMs we refer to as PolyDPG methods. PolyDPG methods are competitive with other existing polygonal methods, since they arise from very different ideas and they inherit many advantages from the DPG methodology. For example, they can be easily generalized to different linear equations; they have a solid mathematical background in terms of proving stability and high-order convergence; they allow for discontinuous material properties while retaining stability; they result in positive-definite stiffness matrices; and they carry a completely natural arbitrary-order a posteriori error estimator, which facilitates implementation of adaptive refinement strategies. The last feature is particularly desirable when combined with polygonal elements, because there is no need for the constrained approximation technology to treat hanging nodes, paving the way for use in applications like dynamic fracture [80, 68, 13, 15] and topology optimization [84, 53, 3, 86]. We complement this article by providing an open-source software in MATLAB®, also named PolyDPG [87].

The outline of the article is as follows. In Section 2 we describe a PolyDPG method for a model problem (Poisson’s equation), along with the DPG solution scheme and the convergence theory (with the proof relegated to Appendix A). In Section 3 several illustrative examples are presented. High-order convergence for different $p$ is verified for both convex and highly distorted concave elements. Then, a physically relevant problem involving discontinuous material properties along an arbitrary interface is solved. Finally, an adaptive refinement strategy is described, successfully implemented, and compared to traditional adaptive schemes. Our concluding remarks are presented in Section 4.

2 PolyDPG methods

Typical FEMs map elements from the actual physical space to a known fixed master element space corresponding to the same element type. For example, in 2D a general quadrilateral in $\mathbb{R}^{2}$ is mapped to a master quadrilateral (typically $(0,1)^{2}$ or $(-1,1)^{2}$). This requires defining a master element for each element type, which is possible for limited types of elements (e.g. quadrilaterals and triangles in 2D, or hexahedra, tetrahedra, triangular prisms and pyramids in 3D), but is usually nonviable when dealing with general polytopes. Thus, as with any polytopal FEM, the idea is to circumvent any master elements by shifting the focus directly to the physical space.

The main issue in doing so is satisfying inter-element continuity of the basis functions, which is required for discretizing Sobolev spaces such as $H^{1}$. This is partly resolved by using generalized barycentric coordinates, but these techniques are usually limited to first order methods (in terms of convergence), and it becomes difficult to discretize other Sobolev spaces such as $\boldsymbol{H}(\operatorname{curl})$ and $\boldsymbol{H}(\operatorname{div})$ even for the lowest order cases [26]. Indeed, even with the “traditional” pyramid element, having high-order discretizations for different spaces is challenging to achieve [76, 51, 55, 1], and so is the case for 2D non-affine quadrilaterals [4]. To overcome this, VEMs concentrate on the boundaries while nonconforming polytopal discontinuous methods, like IPDG, HHO, WG, and HDG (which are closely related [32, 31]), remove the continuity requirements altogether. However, all of these methods need to carefully add (equation-dependent) stabilization or penalty terms [8, 20, 45, 89, 33], and they must account for these in the error analysis, leading to a nonstandard theory of convergence [30].

As will be seen, the discontinuous Petrov-Galerkin (DPG) methodology is very general from a variational standpoint, so it is not limited to the traditional primal and mixed formulations. Thus, without sacrificing any desirable stability properties, it is able to discretize “broken” ultraweak variational formulations, which avoid most inter-element continuity requirements. The only continuity requirements are met by skeleton variables which live on the element boundaries. Technically speaking, the resulting method is still a conforming FEM, and the “standard” error analysis can be applied. This is very useful, because it allows to generalize the method to any well-posed linear equation formulated with traditional functional spaces ($H^{1}$, $\boldsymbol{H}(\operatorname{curl})$, $\boldsymbol{H}(\operatorname{div})$ and $L^{2}$).

In 2D, the polygonal element boundaries are simply line segments, so it is easy to define high-order discretizations along the mesh skeleton. Given that this is less trivial for polyhedra in 3D, we only analyze 2D problems in this introductory paper. We now proceed by introducing the model problem and its corresponding ultraweak formulations in the next section.

2.1 Model problem and ultraweak variational formulations

As a model problem, consider Poisson’s equation coming from the steady-state heat equation in a (heterogeneous) domain $\Omega\subseteq\mathbb{R}^{2}$, where $u$ is the temperature, $\boldsymbol{q}$ is the heat flux, $k>0$ is the variable thermal conductivity, and $r$ is the internal heat source,

[TABLE]

Note that the equation can be written directly as a second order system (left) or as a first order system (right). For simplicity, we assume temperature boundary conditions along all of $\partial\Omega$, so that $u=g$ at $\partial\Omega$, where $g$ is a known function.

To solve the equation using FEMs, a variational form is required, and in this respect, there are many possibilities. For now assume vanishing temperature boundary conditions so that $g=0$. The classical approach stems directly from the second order equation by multiplying by a test function and integrating by parts once, leading to the primal formulation where the solution $u$ is sought in the trial space $\mathscr{U}^{\mathscr{P}}$ and must satisfy

[TABLE]

with $(u,v)_{K}=\int_{K}u\cdot v\operatorname{d}K$ for $K\subseteq\Omega$. Notice in this case $\mathscr{V}^{\mathscr{P}}=\mathscr{U}^{\mathscr{P}}$, so both spaces can be discretized in the same way, leading to the Galerkin method. The same property holds for standard mixed formulations which stem from the first order system. The ultraweak formulation is also derived from the first order system, but here all equations are integrated by parts to pass the derivatives to the test functions. The resulting ultraweak formulation seeks $(u,\boldsymbol{q})=\mathfrak{u}_{0}\in\mathscr{U}_{0}=L^{2}(\Omega)\times\boldsymbol{L}^{2}(\Omega)$ satisfying

[TABLE]

where $\boldsymbol{L}^{2}(\Omega)=(L^{2}(\Omega))^{2}$. Clearly the trial and test spaces in this case are completely different, $\mathscr{U}_{0}\neq\mathscr{V}_{0}$. Thus, to solve this system it is necessary to drift away from the traditional Galerkin method. As we will see, a discretization via minimum residual FEMs is a viable option. It is worth remarking that the primal and ultraweak formulations are mutually well-posed in the infinite-dimensional setting [62, 38, 21]. Since the primal formulation is known to be well-posed in view of the Lax-Milgram theorem and Poincaré’s inequality, so is the ultraweak formulation. This guarantees the existence of a unique solution in the trial space satisfying a stability estimate.

The ultraweak formulation has copies of $L^{2}(\Omega)$ as a trial space, thus its discretization does not require satisfying any inter-element continuity, which is very desirable for polygons. However, all the difficulties are passed to the test space for which inter-element continuity requirements are essential. Fortunately, it is possible to remove these requirements in the test space as well, but at the cost of introducing skeleton variables, as we will see shortly. In fact, the practicality of DPG methods relies on using broken (or discontinuous) test spaces, and this results in a slightly modified formulation called the broken ultraweak formulation, which will be derived in what follows. Consider a mesh (i.e. an open partition), $\mathcal{T}$, of $\Omega$ comprised of (disjoint) elements $K\in\mathcal{T}$, and define the broken spaces and piecewise integration,

[TABLE]

Then, element-wise, multiply by broken test functions $(v,\boldsymbol{\tau})=\mathfrak{v}\in\mathscr{V}=H^{1}(\mathcal{T})\times\boldsymbol{H}(\operatorname{div},\mathcal{T})$, integrate by parts, and sum across all elements. The result is very similar to the ultraweak formulation, but has new terms on the boundaries of the elements involving $u|_{\partial K}$ and $\boldsymbol{q}|_{\partial K}\!\cdot\!\hat{\mathbf{n}}_{K}$, where $\hat{\mathbf{n}}_{K}$ is the outward normal to the element $K$. These terms vanish if the test space is not broken (i.e. $\mathscr{V}_{0}$). Unfortunately, if we want $u\in L^{2}(\Omega)$ and $\boldsymbol{q}\in\boldsymbol{L}^{2}(\Omega)$, then the traces $u|_{\partial K}$ and $\boldsymbol{q}|_{\partial K}\!\cdot\!\hat{\mathbf{n}}_{K}$ technically do not exist [70] and to incorporate them it is necessary to add new skeleton (or interface) variables in the spaces

[TABLE]

where the duality $\langle\cdot,\cdot\rangle_{\partial K}$ can be thought of as a boundary integral (for smooth enough inputs it is actually a boundary integral). Therefore, the resulting broken ultraweak variational formulation seeks

[TABLE]

such that

[TABLE]

where v_{\partial\mathcal{T}}=\prod_{K\in\mathcal{T}}(v|_{K})\big{|}_{\partial K} and \boldsymbol{\tau}_{\partial\mathcal{T}}=\prod_{K\in\mathcal{T}}(\boldsymbol{\tau}|_{K})\big{|}_{\partial K}\!\cdot\!\hat{\mathbf{n}}_{K}. This formulation can also be proved to be well-posed, with stability properties independent of the choice of the mesh [21, 62]. With nontrivial boundary conditions, $g\neq 0$, simply consider $\ell(\mathfrak{v})=(r,v)_{\mathcal{T}}-\langle\widetilde{g}_{\partial\mathcal{T}},\boldsymbol{\tau}_{\partial\mathcal{T}}\rangle_{\partial\mathcal{T}}$ instead, where $\widetilde{g}\in H^{1}(\Omega)$ is an extension of g\in H^{\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}}(\partial\Omega)=\{f=\widetilde{f}|_{\partial\Omega}\mid\widetilde{f}\in H^{1}(\Omega)\}, and add $\widetilde{g}$ to the solution $u$ of (2.7) to obtain the final temperature.

Despite looking intricate, the broken ultraweak variational formulation has the advantage of removing much of the inter-element compatibility conditions, since some of its trial variables are in $L^{2}$ and its test variables are discontinuous along the elements. The only inter-element compatibility is due to the skeleton variables, which reside solely on the element boundaries. In 2D, as we mentioned before, this is extremely convenient since the element boundaries are simply 1D line segments.

2.2 Discretization and the DPG methodology

In this section we present the procedure of discretizing the ultraweak formulations. The Galerkin method is the widely used approach for conventional formulations. It employs the same test and trial spaces, leading to a square linear system of equations. Indeed, consider the primal formulation in (2.2), with $\{\mathfrak{u}^{\mathscr{P}}_{j}\}_{j=1}^{N}$ being a basis for the discrete subspaces $\mathscr{U}_{h}^{\mathscr{P}}=\mathscr{V}_{h}^{\mathscr{P}}\subseteq\mathscr{U}^{\mathscr{P}}=\mathscr{V}^{\mathscr{P}}$. Then, the discrete solution $u_{h}=\sum_{j=1}^{N}(\mathsf{u}_{h})_{j}\mathfrak{u}^{\mathscr{P}}_{j}\in\mathscr{U}_{h}^{\mathscr{P}}$ for $\mathsf{u}_{h}\in\mathbb{R}^{N}$, satisfies

[TABLE]

where $\mathsf{B}_{ij}^{\mathscr{P}}=b^{\mathscr{P}}(\mathfrak{u}^{\mathscr{P}}_{j},\mathfrak{v}^{\mathscr{P}}_{i})$ and \mathchoice{\raisebox{0.0pt}{\displaystyle{}\text{{l}}}}{\raisebox{0.0pt}{\textstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\text{{l}}}}_{i}^{\mathscr{P}}=\ell^{\mathscr{P}}(\mathfrak{v}^{\mathscr{P}}_{i}) with $\mathfrak{v}^{\mathscr{P}}_{i}=\mathfrak{u}^{\mathscr{P}}_{i}$, so that $\mathsf{B}^{\mathscr{P}}\in\mathbb{R}^{N\!\times\!N}$ and \mathchoice{\raisebox{0.0pt}{\displaystyle{}\text{{l}}}}{\raisebox{0.0pt}{\textstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\text{{l}}}}^{\mathscr{P}}\in\mathbb{R}^{N}. The basis functions, $\mathfrak{u}^{\mathscr{P}}_{j}$, are chosen with a very small support not exceeding a few neighboring elements, resulting in a computationally practical method due to the sparse structure of $\mathsf{B}^{\mathscr{P}}$.

In general, when the trial and test spaces are different, $\mathscr{U}\neq\mathscr{V}$, this approach is still possible but requires finding bases $\{\mathfrak{u}_{j}\}_{j=1}^{N}$ and $\{\mathfrak{v}_{i}\}_{i=1}^{N}$ for $\mathscr{U}_{h}\subseteq\mathscr{U}$ and $\mathscr{V}_{h}\subseteq\mathscr{V}$ respectively. However, two issues immediately arise. First, the canonical polynomial-based discrete basis of $\mathscr{V}_{h}\subseteq\mathscr{V}$ typically is not of size $N$ (the same size of the basis for $\mathscr{U}_{h}$). Second, even if a nonstandard basis for $\mathscr{V}_{h}$ of the right size is found, the resulting numerical method could very well be unstable, meaning that the inf-sup inequality,

[TABLE]

might not hold. In fact, depending on the equation and mesh size, even the Galerkin method can be unstable. Minimum residual finite element methods overcome these two difficulties by design.

Let $\mathscr{U}^{\prime}$ and $\mathscr{V}^{\prime}$ be the dual spaces to $\mathscr{U}$ and $\mathscr{V}$ respectively, and define $\mathscr{B}:\mathscr{U}\to\mathscr{V}^{\prime}$ and its adjoint $\mathscr{B}^{\prime}:\mathscr{V}\to\mathscr{U}^{\prime}$ through duality pairings as $\langle\mathscr{B}\mathfrak{u},\mathfrak{v}\rangle=b(\mathfrak{u},\mathfrak{v})=\langle\mathfrak{u},\mathscr{B}^{\prime}\mathfrak{v}\rangle$. Then, for a discrete trial space $\mathscr{U}_{h}\subseteq\mathscr{U}$, minimum residual methods seek the minimizer of the residual [39, 62],

[TABLE]

where $\mathscr{R}_{\mathscr{V}}:\mathscr{V}\to\mathscr{V}^{\prime}$ is the Riesz map, which is defined by duality as $\langle\mathscr{R}_{\mathscr{V}}\mathfrak{v},\delta\mathfrak{v}\rangle=(\mathfrak{v},\delta\mathfrak{v})_{\mathscr{V}}$, with $(\cdot,\cdot)_{\mathscr{V}}$ being the inner product of the Hilbert space $\mathscr{V}$. Here, $\mathscr{V}^{\mathrm{opt}}=\mathscr{R}_{\mathscr{V}}^{-1}\mathscr{B}\mathscr{U}_{h}$ is called the optimal test space, because this exact choice of discrete test space automatically results in the best inf-sup stable discrete method satisfying (2.9) [39]. Given an element of the basis for $\mathscr{U}_{h}$, $\mathfrak{u}_{i}\in\{\mathfrak{u}_{j}\}_{j=1}^{N}$, the corresponding optimal test function is $\mathfrak{v}_{i}^{\mathrm{opt}}=\mathscr{R}_{\mathscr{V}}^{-1}\mathscr{B}\mathfrak{u}_{i}$. With these choices the resulting matrix $\mathsf{B}^{\mathrm{opt}}_{ij}=b(\mathfrak{u}_{j},\mathfrak{v}_{i}^{\mathrm{opt}})$, called the optimal stiffness matrix, is always symmetric positive-definite.

Unfortunately, computing $\mathscr{R}_{\mathscr{V}}^{-1}$ is impossible since $\mathscr{V}$ is infinite-dimensional. Thus, minimum residual methods simply make a choice of an enriched test space $\mathscr{V}_{r}\subseteq\mathscr{V}$ (with $M=\dim(\mathscr{V}_{r})\geq\dim(\mathscr{U}_{h})=N$) over which the operator is inverted. The advantage is that this enriched space may be discretized with a standard canonical polynomial-based basis, $\{\mathfrak{v}_{i}\}_{i=1}^{M}$, and ultimately the resulting near-optimal space is $\mathscr{V}_{h}=\mathscr{V}^{\mathrm{n\text{-}opt}}=\mathscr{R}_{\mathscr{V}_{r}}^{-1}\mathscr{B}\mathscr{U}_{h}$ and its corresponding near-optimal basis is $\mathfrak{v}_{i}^{\mathrm{n\text{-}opt}}=\mathscr{R}_{\mathscr{V}_{r}}^{-1}\mathscr{B}\mathfrak{u}_{i}$ for every $\mathfrak{u}_{i}\in\{\mathfrak{u}_{j}\}_{j=1}^{N}$. The resulting discrete method can be shown to be equivalent to the linear system,

[TABLE]

where $\mathfrak{u}_{h}=\sum_{j=1}^{N}(\mathsf{u}_{h})_{j}\mathfrak{u}_{j}\in\mathscr{U}_{h}$ is the discrete solution; the Gram matrix $\mathsf{G}_{ij}=(\mathfrak{v}_{i},\mathfrak{v}_{j})_{\mathscr{V}}$ is a discretization of $\mathscr{R}_{\mathscr{V}_{r}}$; $\mathsf{B}_{ij}=b(\mathfrak{u}_{j},\mathfrak{v}_{i})$ and \mathchoice{\raisebox{0.0pt}{\displaystyle{}\text{{l}}}}{\raisebox{0.0pt}{\textstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\text{{l}}}}_{i}=\ell(\mathfrak{v}_{i}) are called the enriched stiffness matrix and load; and $\mathsf{B}^{\mathrm{n\text{-}opt}}_{ij}=b(\mathfrak{u}_{j},\mathfrak{v}_{i}^{\mathrm{n\text{-}opt}})$ and \mathchoice{\raisebox{0.0pt}{\displaystyle{}\text{{l}}}}{\raisebox{0.0pt}{\textstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\text{{l}}}}_{i}^{\mathrm{n\text{-}opt}}=\ell(\mathfrak{v}_{i}^{\mathrm{n\text{-}opt}}) are the near-optimal stiffness matrix and load. Clearly the enriched stiffness matrix is rectangular and tall, $\mathsf{B}\in\mathbb{R}^{M\!\times\!N}$ with $M\geq N$, while the near-optimal stiffness matrix is square and symmetric positive-definite, $\mathsf{B}^{\mathrm{n\text{-}opt}}\in\mathbb{R}^{N\!\times\!N}$. To implement, one has to form the Gram matrix ($\mathsf{G}\in\mathbb{R}^{M\!\times\!M}$), enriched stiffness matrix ($\mathsf{B}\in\mathbb{R}^{M\!\times\!N}$) and enriched load vector (\mathchoice{\raisebox{0.0pt}{\displaystyle{}\text{{l}}}}{\raisebox{0.0pt}{\textstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\text{{l}}}}\in\mathbb{R}^{M}) first, then calculate the near-optimal stiffness matrix ($\mathsf{B}^{\mathrm{n\text{-}opt}}=\mathsf{B}^{\mathsf{T}}\mathsf{G}^{-1}\mathsf{B}\in\mathbb{R}^{N\!\times\!N}$) and near-optimal load vector (\mathchoice{\raisebox{0.0pt}{\displaystyle{}\text{{l}}}}{\raisebox{0.0pt}{\textstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\text{{l}}}}^{\mathrm{n\text{-}opt}}=\mathsf{B}^{\mathsf{T}}\mathsf{G}^{-1}\mathchoice{\raisebox{0.0pt}{\displaystyle{}\text{{l}}}}{\raisebox{0.0pt}{\textstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptstyle{}\text{{l}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\text{{l}}}}\in\mathbb{R}^{N}), and finally solve for the basis coefficients of the discrete solution ($\mathsf{u}_{h}\in\mathbb{R}^{N}$).

All this derivation holds for any arbitrary linear variational formulation including the ultraweak formulations in (2.3) and (2.7). The method is near-optimal in that it is designed to approximate the optimal method (with $\mathsf{B}^{\mathrm{opt}}$), so in principle it is not known to be stable, but in practice it typically is or can be made stable (if it is not stable simply enrich $\mathscr{V}_{r}$ even more so that $M\gg N$). In fact, the stability of the near-optimal method can rigorously be proved by constructing a Fortin operator, $\Pi_{F}:\mathscr{V}\to\mathscr{V}_{r}$ [58, 21].

However, there are major differences between applying this method to the ultraweak formulation in (2.3) and the broken ultraweak formulation in (2.7). Namely, for the standard ultraweak formulation the enriched (sparse) stiffness matrix, $\mathsf{B}$, and the Gram matrix, $\mathsf{G}$, are assembled globally first and then the near-optimal stiffness matrix, $\mathsf{B}^{\mathrm{n\text{-}opt}}$ is computed using (2.11). This is very expensive, especially due to the inversion of $\mathsf{G}$. Thus, despite many advantages, the method is not very practical. On the other hand, when using broken test spaces, as in the broken ultraweak formulation, the matrix $\mathsf{G}$ has a disjoint diagonal block structure, where each block corresponds to one element. Hence, the Gram matrix can be inverted locally, allowing the local near-optimal stiffness matrices $(\mathsf{B}^{\mathrm{n\text{-}opt}})_{K}$ to be computed directly for each element $K\in\mathcal{T}$. This is turn allows $\mathsf{B}^{\mathrm{n\text{-}opt}}$ to be assembled as in any other FEM. Thus, using formulations with broken test spaces localizes the computations and parallelizes the assembly, making it a practical FEM. However, when compared to traditional FEMs, the local computations are more expensive due to the additional skeleton variables. Note that the broken ultraweak formulation in (2.7) has an enriched stiffness matrix with the structure,

[TABLE]

where $(\mathsf{B}_{0})_{ij}=b_{0}((\mathfrak{u}_{0})_{j},\mathfrak{v}_{i})$ and $\hat{\mathsf{B}}_{ij}=\hat{b}(\hat{\mathfrak{u}}_{j},\mathfrak{v}_{i})$, with the $\mathscr{U}_{h}$-basis $\{\mathfrak{u}_{j}\}_{j=1}^{N}=\{((\mathfrak{u}_{0})_{j},0)\}_{j=1}^{N_{0}}\cup\{(0,\hat{\mathfrak{u}}_{j})\}_{j=1}^{\hat{N}}$ so that $N=N_{0}+\hat{N}$, and similarly with the other sub-blocks.

In the literature, the application of minimum residual methods to variational formulations with broken test spaces is referred to as the DPG methodology. The methodology is quite general as it can be applied to variational formulations other than the broken ultraweak such as broken primal or broken mixed formulations [62, 21]. Each application case results in a different DPG method similar to how the Galerkin methodology can be applied to primal and mixed formulations (where $\mathscr{U}_{h}=\mathscr{V}_{h}$). Nonetheless, the lack of inter-element compatibility restrictions on the $\mathscr{U}_{0}$-part of the trial space (which lies in copies of $L^{2}$) makes the ultraweak formulation a natural candidate to develop a DPG method for polygonal elements.

It is worth mentioning that the DPG methodology carries a natural arbitrary-order residual-based a posteriori error estimator. The expression for the residual is,

[TABLE]

where $\mathfrak{u}_{h}$ (and $\mathsf{u}_{h}$) is the solution. Note that the test spaces are broken, so the computations can be performed locally. Therefore, (2.13) can serve as an a posteriori error estimator for driving different adaptive strategies [62, 65]. Adaptivity in its own right is a very interesting subject of study for polygonal elements, as they provide great flexibility for the implementation of such strategies without resorting to constrained approximations to deal with hanging nodes. More details on this will be given in Section 3.4.

A final comment on minimum residual FEMs, including all DPG methods, is that the choice of test norm (or inner product) for $\mathscr{V}$, which appears in the computation of $\mathsf{G}$, has a significant influence. Generally speaking, the standard norms are usually chosen as test norms. For example, the standard norm for the broken ultraweak formulation in (2.7) is,

[TABLE]

However, there are other norms that still make $\mathscr{V}$ a Hilbert space but lead to different results. Specifically for the broken ultraweak formulations, the adjoint graph norm has interesting properties [39]. Using the ultraweak formulation in (2.3), the first two terms in this norm can be derived as,

[TABLE]

where $\|\cdot\|_{L^{2}(\mathcal{T})}^{2}=(\cdot,\cdot)_{\mathcal{T}}$ and the same with $\|\cdot\|_{\boldsymbol{L}^{2}(\mathcal{T})}^{2}$. The third term, which has the $\varepsilon^{2}$ factor, makes the norm localizable, because otherwise (2.15) would not be a norm for arbitrary broken functions $v\in H^{1}(\mathcal{T})$ (although it would be a norm for $v\in H_{0}^{1}(\Omega)$). One can choose an arbitrary value for $\varepsilon>0$, but using small values of $\varepsilon$ (with the caveat of ill-conditioned local problems) is of particular interest for certain equations, such as Helmholtz [57]. Note that the corresponding inner products for the (real-valued) Hilbert space $\mathscr{V}$ can be derived from the polarization identity, (\mathfrak{v}_{1},\mathfrak{v}_{2})_{\mathscr{V}}=\frac{1}{4}\big{(}\|\mathfrak{v}_{1}+\mathfrak{v}_{2}\|_{\mathscr{V}}^{2}-\|\mathfrak{v}_{1}-\mathfrak{v}_{2}\|_{\mathscr{V}}^{2}\big{)}.

2.3 Choice of trial and test spaces

The choice of trial and test spaces is important to establish the method’s convergence. As mentioned before, strict inter-element compatibility requirements leaves very limited options. Particularly, the problem seems to be extremely complicated for general polygons with high-order discretizations. Fortunately, the $\mathscr{U}_{0}$ trial space component of the broken ultraweak formulation in (2.6) consists of copies of $L^{2}$, so its discretization can be discontinuous across the elements. Moreover, the test spaces are broken, so their discretization should be discontinuous across elements too. This freedom allows one to create bases locally, disregarding the neighboring elements. In particular, bases may be defined by restriction (to the polygonal element of interest), as we will see next.

Our procedure is similar to that in [20] where a bounding box was utilized, but we use a bounding triangle instead. First, the centroid of the polygon and the furthest vertex from the centroid are determined. Next, a bounding circle centered at the centroid and passing through the furthest vertex is defined. Then, the bounding equilateral triangle inscribing the circle is computed such that one of its edge-midpoints is the polygon’s furthest vertex. This is shown in Figure 1. Lastly, the “usual” high-order polynomial shape functions for the triangle are used and then restricted to the polygon. We use the term “usual” liberally, but to clarify, we include further details below.

There are several spaces at the infinite-dimensional level which we want to discretize using this technique. Namely, the test space components, $H^{1}(\mathcal{T})$ and $\boldsymbol{H}(\operatorname{div},\mathcal{T})$, and the $\mathscr{U}_{0}$ trial space component, which may be represented by $L^{2}(\Omega)$. Following our technique, the procedure reduces to finding the local discretizations of $H^{1}(T_{K})$, $\boldsymbol{H}(\operatorname{div},T_{K})$ and $L^{2}(T_{K})$, where $T_{K}$ is the bounding triangle of the polygonal element $K\in\mathcal{T}$. These three spaces actually form a differential de Rahm exact sequence, and it is convenient that their respective discretizations do too. For triangles, this is satisfied by the classical Nédélec sequence of the first type [44, 51],

[TABLE]

where $\mathcal{P}^{p}(T_{K})$ are the polynomials in $x=(x_{1},x_{2})$ of total order less than or equal $p\in\mathbb{N}$, the 2D Raviart-Thomas space is $\mathcal{R}\mathcal{T}^{p}(T_{K})=(\mathcal{P}^{p-1}(T_{K}))^{2}+x\mathcal{P}^{p-1}(T_{K})$ (a rotation of the 2D Nédélec space), and the 2D scalar-to-vector curl operator is defined as \operatorname{curl}(u)=\big{(}\begin{smallmatrix}0&1\\ -1&0\end{smallmatrix}\big{)}\nabla u for any $u\in H^{1}(T_{K})$. Notice that the parameter $p$ represents the order of the discrete sequence and does not necessarily coincide with the order of the polynomials of a particular discretization. For example if $p=3$, the discretization of $L^{2}(T_{K})$ are the polynomials of at most total order $p-1=2$.

This sequence has many desirable properties, and precisely because of these, we prefer to use a bounding triangle instead of a bounding box. In particular, the spaces are invariant under affine transformations (the spaces remain the same even if the bounding triangle is arbitrarily rotated about the polygon centroid); the overall drop of polynomial order across the sequence is one (from $\mathcal{P}^{p}(T_{K})$ to $\mathcal{P}^{p-1}(T_{K})$); the approximation properties are suitable (see Appendix A); and they are the smallest possible spaces with all these properties (see [5, §3.4]).

Having said that, a similar procedure can be carried out for a bounding box, $Q_{K}$ of $K\in\mathcal{T}$, where the spaces become

[TABLE]

with $\mathcal{Q}^{p,q}(Q_{K})=\mathcal{P}^{p}(x_{1})\otimes\mathcal{P}^{q}(x_{2})$.

In either case, the final spaces for the polygon $K\subseteq T_{K}$ (or $K\subseteq Q_{K}$) are defined by restricting the domain to $K\in\mathcal{T}$, so we denote them by $\mathcal{P}^{p}(K)$ and $\mathcal{R}\mathcal{T}^{p}(K)$ (or $\mathcal{Q}^{p,p}(K)$) instead.

The only remaining spaces to specify are those of the skeleton variables lying in the $\hat{\mathscr{U}}$ trial space component (see (2.6)). These can also be deduced using the same philosophy of exact sequences, but utilizing the traces instead. Indeed, the spaces H_{0}^{\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}}(\partial\mathcal{T}) and H^{-\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}}(\partial\mathcal{T}) are merely $\mathcal{T}$-tuples of compatible traces of $H^{1}(K)$ and normal-traces of $\boldsymbol{H}(\operatorname{div},K)$ respectively. If two elements of different type (a triangle and a quadrilateral) share an edge, the discrete spaces should be compatible across that edge. This is the case when considering the $H^{1}$-discretizations of triangles and quadrilaterals: even though the discretizations themselves are different ($\mathcal{P}^{p}$ and $\mathcal{Q}^{p,p}$), their restrictions to edges are exactly the same, $\mathcal{P}^{p}(e)$, where $e$ represents an edge parametrized linearly by $t_{e}$. The same occurs with the $\boldsymbol{H}(\operatorname{div})$-discretizations, which have $\mathcal{P}^{p-1}(e)$ as normal-trace along the edges. Additionally, the $H^{1}$-discretizations should be compatible at vertices. This is consistent with 1D discretizations of $H^{1}$ and $L^{2}$, which also form an exact sequence, but instead occurring along the boundary of each element and being edge-parametrized along all edges (see [51, §1.6]). This pattern should hold for arbitrary polygons as well. For this, let $\mathcal{E}(K)$ be the set of edges of a polygon $K\in\mathcal{T}$, and define the local discretizations,

[TABLE]

where $C^{0}(\partial K)$ are the continuous functions in $\partial K$ (the intersection ensures that values of neighboring edges coincide at a common vertex), and the local trace spaces are H^{\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}}(\partial K)=\{\hat{u}_{K}=u|_{\partial K}\mid u\in H^{1}(K)\} and H^{-\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}}(\partial K)=\{(\hat{q}_{\hat{\mathbf{n}}})_{K}=\boldsymbol{q}|_{\partial K}\!\cdot\!\hat{\mathbf{n}}_{K}\mid\boldsymbol{q}\in\boldsymbol{H}(\operatorname{div},K)\}.

Now we have enough information to actually globally define the discrete trial space. For a value of $p\in\mathbb{N}$, it is

[TABLE]

Notice that the condition $(u,\boldsymbol{q},\hat{u},\hat{q}_{\hat{\mathbf{n}}})\in\mathscr{U}$ (so $(\hat{u},\hat{q}_{\hat{\mathbf{n}}})\in\hat{\mathscr{U}}$) implies that $\hat{u}$ vanishes at the boundaries, that $\hat{u}_{K_{1}}|_{e}=\hat{u}_{K_{2}}|_{e}$, and that $(\hat{q}_{\hat{\mathbf{n}}})_{K_{1}}|_{e}=-(\hat{q}_{\hat{\mathbf{n}}})_{K_{2}}|_{e}$, where $e$ is a common edge between the elements $K_{1}$ and $K_{2}$. No such compatibility implications exist for $(u,\boldsymbol{q})\in\mathscr{U}_{0}$.

For the enriched test space, the discretizations are chosen from a sequence of order $p+\Delta p$, and we say the space is $p$-enriched, so that

[TABLE]

The notation $\Delta p_{K}$ indicates that this value is element-dependent. In fact, recall that for minimum residual methods to work, $M=\dim(\mathscr{V}_{r})\geq\dim(\mathscr{U}_{h})=N$, and this restriction on the dimensionality should hold locally as well. Thus, $\Delta p_{K}$ has to be chosen such that this condition holds. This is important for the polygonal element methods, because when a polygon has many sides, the size of the local trial space may be quite large and a large value of $\Delta p_{K}$ must be chosen for that particular element.

To elaborate, consider an interior $n$-sided polygonal element $K$ (so that $\partial K\cap\partial\Omega=\varnothing$). Its local trial and test space dimensions would be

[TABLE]

Thus, for $p=2$ and $n=3$ (a triangle), $\dim(\mathscr{U}_{h}(K))=21$, so that a value of $\Delta p_{K}=1$ is sufficient ($\dim(\mathscr{V}_{r}(K))=25$); but if $p=2$ and $n=8$ (an octagon), $\dim(\mathscr{U}_{h}(K))=41$, a value of at least $\Delta p_{K}=3$ (so that $\dim(\mathscr{V}_{r}(K))=56$) is required. Having said that, sometimes for simplicity a valid value of $\Delta p$ is chosen uniformly throughout the mesh (this is the case for all of our examples in Section 3).

To illustrate, some representative shape functions of the components of $\mathscr{U}_{h}(K)$ and $\mathscr{V}_{r}(K)$ are shown in Figure 2 for the different energy spaces and multiple values of $p$.

We refer to the high-order polygonal DPG method resulting from this choice of trial and enriched test spaces as a PolyDPG method for Poisson’s equation. However, it can easily be generalized to ultraweak formulations coming from other linear equations (see Remark 2.2 later), so it is more appropriate to allude to a family of PolyDPG methods. Note that the methods seem to be very expensive due to the large number of variables in the trial space $\mathscr{U}_{h}$, but this is deceiving. In fact, all of the $\mathscr{U}_{0}$ trial space components can be statically condensed locally for ultraweak formulations, meaning that this part of the near-optimal stiffness matrix, $\mathsf{B}^{\mathrm{n\text{-}opt}}$, can be effectively removed by taking Schur complements. Thus, the only remaining connectivity is that coming from the skeleton variables in $\hat{\mathscr{U}}$. So computationally speaking, solving with these variational formulations is not as costly as one might initially imagine.

2.4 Convergence

Since the subspaces used to discretize the ultraweak variational formulation are, rigorously speaking, subsets of the infinite dimensional trial and test spaces, PolyDPG methods are conforming FEMs. Thus, the “standard” convergence theory can be applied. However, this is an understatement because the skeleton variables are not standard, so they require a careful treatment. The details are left to Appendix A, but the main result is stated here along with the key assumptions.

Definition**.**

A collection of subsets of $\mathbb{R}^{2}$, $\mathcal{T}_{\mathcal{K}}$, is said to have the finite overlap condition if

[TABLE]

For a family of such collections given by a parameter $\mathfrak{h}\in\mathfrak{H}$, $\{\mathcal{T}_{\mathcal{K},\mathfrak{h}}\}_{\mathfrak{h}\in\mathfrak{H}}$, the finite overlap condition is said to be robust in $\mathfrak{h}$ if there exists an integer $M_{\mathrm{ov}}>0$, independent of $\mathfrak{h}$, such that $\mathrm{ov}(\mathcal{T}_{\mathcal{K},\mathfrak{h}})\leq M_{\mathrm{ov}}$ for any $\mathfrak{h}\in\mathfrak{H}$.

Definition**.**

A triangulation $\mathscr{T}(K)=\{\mathscr{T}_{i}(K)\}_{i\in I_{K}}$ (with $I_{K}$ finite) of a (simple) polygonal element $K$ is said to be edge-compatible if for each edge of $K$, only one $\mathscr{T}_{i}(K)$ shares that edge. For any polygon such a triangulation is known to exist [71, 24, 2]. The triangulation is additionally said to be shape-regular if all $\mathscr{T}_{i}(K)$ satisfy a kind of uniform shape-regularity condition (e.g. they satisfy a minimum angle condition or the ratio of their diameters to their incircle radii remains bounded).

Theorem 2.1**.**

Let $p\in\mathbb{N}$ be a polynomial order and $\{\mathcal{T}_{\mathfrak{h}}\}_{\mathfrak{h}\in\mathfrak{H}}$ be a family of polygonal meshes discretizing the domain $\Omega$, such that there exist shape-regular edge-compatible triangulations for all $K\in\mathcal{T}_{\mathfrak{h}}$ with a robust shape-regularity condition independent of $K\in\mathcal{T}_{\mathfrak{h}}$ across all $\mathfrak{h}\in\mathfrak{H}$. Assume that the associated collections of bounding triangles (see Section 2.3), $\{\mathcal{T}_{T,\mathfrak{h}}\}_{\mathfrak{h}\in\mathfrak{H}}=\{\{T_{K}\}_{K\in\mathcal{T}_{\mathfrak{h}}}\}_{\mathfrak{h}\in\mathfrak{H}}$, where $T_{K}$ is the bounding triangle of a polygonal element $K$, satisfy a robust finite overlap condition. Also, assume the existence of a linear and continuous Fortin operator, $\Pi_{F}:\mathscr{V}\to\mathscr{V}_{r}$, satisfying the orthogonality condition, $b(\mathfrak{u}_{h},\mathfrak{v}-\Pi_{F}\mathfrak{v})=0$, for all $\mathfrak{u}_{h}\in\mathscr{U}_{h}$ and $\mathfrak{v}\in\mathscr{V}$, and with a continuity bound, $M_{F}>0$ (so $\|\Pi_{F}\mathfrak{v}\|_{\mathscr{V}}\leq M_{F}\|\mathfrak{v}\|_{\mathscr{V}}$ for all $\mathfrak{v}\in\mathscr{V}$), where $b:\mathscr{U}\times\mathscr{V}\to\mathbb{R}$, $\ell:\mathscr{V}\to\mathbb{R}$, $\mathscr{U}_{h}$ and $\mathscr{V}_{r}$ are given in (2.6), (2.7), (2.19) and (2.20). Then, the problem of finding $\mathfrak{u}_{h}\in\mathscr{U}_{h}$ such that

[TABLE]

has a unique solution. When compared to the unique solution of the infinite dimensional problem, $\mathfrak{u}\in\mathscr{U}$ (so $b(\mathfrak{u},\mathfrak{v})=\ell(\mathfrak{v})$ for all $\mathfrak{v}\in\mathscr{V}$), and assuming it is regular enough, $\mathfrak{u}\in\mathscr{U}^{s}\subseteq\mathscr{U}$, for an $s>\frac{1}{2}$, the following $h$-convergence estimate holds provided $M_{F}$ is independent of $\mathfrak{h}$,

[TABLE]

where $h=\sup_{K\in\mathcal{T}_{\mathfrak{h}}}\operatorname{diam}(K)$ and $C=C(s,p,\Omega)>0$ is a constant independent of $\mathfrak{h}$ (and so of $h$ as well). For more details about $s>\frac{1}{2}$ and $\mathscr{U}^{s}$ see Appendix A. Moreover, if $M_{F}$ is $p$-independent as well, then in the $p$-asymptotic limit $C=\widetilde{C}(\ln p)^{2}p^{-s}$ where $\widetilde{C}=\widetilde{C}(s,\Omega)$ is independent of $p$.

Remark 2.1**.**

The robust finite overlap condition is also assumed in [20], and is not a very restrictive assumption. It is used in the proof to establish a robust finite constant for the global $L^{2}(\Omega)$ convergence estimates (details are in Appendix A). On the other hand, the robust shape-regular edge-compatible triangulation of all elements is a more restrictive assumption, but it is necessary to prove the convergence estimates of the skeleton variables.

Remark 2.2**.**

As shown in Appendix A, the theorem actually holds for any well-posed broken ultraweak variational formulation with trial variables in $L^{2}(\Omega)$ and skeleton (also trial) variables in subsets of H^{\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}}(\partial\mathcal{T}) and H^{-\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}}(\partial\mathcal{T}). Thus, this result also holds for other equations such as linear elasticity, acoustics, and convection-dominated diffusion.

Remark 2.3**.**

The arguments can be easily extended to a 3D mesh with polyhedral elements provided all the faces of the polyhedra are triangular. Then, the proof would even hold for equations involving skeleton variables representing the traces of $\boldsymbol{H}(\operatorname{curl},\Omega)$ spaces, like an ultraweak formulation of Maxwell’s equations (see [21]). However, the problem (and the corresponding numerical implementation) is more challenging for general polyhedra in 3D.

3 Numerical examples

In this section we consider several examples to examine the performance of the PolyDPG method. In all cases, Poisson’s equation representing the nondimensionalized steady-state heat equation was solved in the domain $\Omega=(0,1)^{2}$. Unless otherwise stated, bounding triangles were utilized (as opposed to bounding boxes) and the (nondimensional) conductivity was taken as $k=1$. Also, a default uniform value of $\Delta p=1$ was used, but was increased (uniformly across the mesh, for the sake of simplicity) if deemed necessary (see (2.21) in Section 2.3). For all computations, the adjoint graph norm written in (2.15) with $\varepsilon=1$ was used as the test space norm.

In the first example, we studied nontrivial meshes with $n$-sided convex polygons. In the second example, we considered highly distorted concave elements in the mesh. The third example was inspired by problems in geoscience, where arbitrary faults separating different material properties occur. To model this, we cut a uniform grid at an angle, so that the resulting mesh had different polygons (pentagons, quadrilaterals and triangles) with discontinuous material properties at each side of the cut. In these three examples, “uniform” refinements were analyzed for different values of $p\in\mathbb{N}$, in the sense that the largest element diameter was roughly cut in half with each refinement. In the final example, we described a polygonal adaptivity scheme by using the PolyDPG arbitrary-order a posteriori error estimator, and compared it with conventional adaptive methods (using standard element shapes). This is particularly important since adaptive refinement algorithms applied to polygonal elements have applications in topology optimization [84, 53, 3, 86] and crack propagation [80, 68].

Note that in all examples we only report the relative error in the $\mathscr{U}_{0}$ trial space component. This is because a rigorous computation of the norms in the $\hat{\mathscr{U}}$ trial space component is simply not viable. The $\mathscr{U}_{0}$ relative error is defined as

[TABLE]

where $\mathfrak{u}_{0}$ is the exact solution and $(\mathfrak{u}_{0})_{h}$ is the computed solution from the PolyDPG method.

Remark 3.1** (PolyDPG software).**

Implementation of PolyDPG methods may deceptively appear difficult when compared to typical FEM algorithms, so we developed an open-source code written in MATLAB® also called PolyDPG [87]. It can be run sequentially or in parallel, and it supports both conventional and polygonal elements. We hope this removes some qualms related to the implementation and makes DPG methods more accessible to other researchers. The shape functions used in the code were originally described in [51] (see Figure 2). The numerical integration was carried out by splitting the polygons into triangles (through Delaunay triangulation), followed by using Gaussian quadrature for each triangle (the Gaussian quadrature points and weights were carefully mapped back from a square), so that polynomial integrands of a certain order were computed up to machine precision.

3.1 Mesh with convex polygons

In this example, we investigated meshes with $n$-sided convex polygonal elements. The software PolyMesher [85] was used to generate the polygonal meshes. In Figure 3 an initial mesh and three subsequent refinements are displayed. The elements are colored according to their number of sides, ranging from $4$ (quadrilaterals) to $7$ (heptagons). We used the manufactured solution,

[TABLE]

for $(x,y)\in\Omega=(0,1)^{2}$ to determine the forcing, i.e. the internal heat source $r$ in (2.1), and the boundary conditions of $u$ at $\partial\Omega$.

As mentioned before, given a trial space associated to a parameter $p$, the corresponding (uniform) value of $\Delta p$ was calculated from (2.21) (using the polygon with the greatest number of sides). Given the presence of hexagons and heptagons, this meant that $\Delta p=2$ was required when $p=1,2$, while $\Delta p=3$ was needed when $p=3,4$. The numerical results are plotted and presented in Figure 4 for $p=4$, including the skeleton temperature, temperature, and heat flux. Additionally, the relative error, calculated using (3.1), is shown in Figure 5, where the expected $h$-convergence rates can be observed for all values of $p$ (the behavior is of the form $h^{p}$ as established by Theorem 2.1). Note that the number of degrees of freedom, $N_{\mathrm{dof}}$, is proportional to $h^{2}$. Thus, the log-log slope indicators in Figure 5 display a $2$ in the $N_{\mathrm{dof}}$-direction, while the other label corresponds to the $h$-convergence rate, $\widetilde{p}$ (so that $\frac{\widetilde{p}}{2}$ is the $N_{\mathrm{dof}}$-convergence rate).

3.2 Mesh with distorted elements

To demonstrate the distortion tolerance of PolyDPG methods, we considered a mesh with highly distorted quadrilaterals, including concave elements. The pattern was then scaled and tessellated to produce the refinements shown in Figure 6. This example is challenging in the sense that other numerical methods likely fail due to the degeneration of either the parametric mapping or the barycentric coordinates associated with the highly distorted elements [67, 61]. The same problem as in Section 3.1 was solved (see (3.2) for manufactured solution). The solution values and $h$-convergence rates for $1\leq p\leq 4$ are shown in Figures 7 and 8 respectively. The expected convergence behavior was observed, showing the flexibility of PolyDPG methods to deal with irregular elements.

3.3 Interface problem

The inspiration behind this example came from geoscience applications where faults abruptly separate the material properties within a domain. Here we considered a domain composed of two materials with different heat conductivities, which share an interface (for simplicity a straight line at an arbitrary angle dividing the square). The heat conductivities are assumed to be uniform on each side of the interface, taking values $k_{1}$ and $k_{2}$, as depicted in Figure 9.

To model certain interfaces one would need unstructured grids. However, by using PolyDPG methods we are able to consider a uniform background grid and simply cut the elements through the interface, leading to the creation of triangles, right trapezoids and pentagons near the interface. In fact, to refine the mesh, first the background mesh was uniformly refined, and then the elements were cut by the interface line. There is one caveat which is only evident for high values of $p$ or small values of $h$: when extremely small triangles (compared to their neighbors) are formed, the assembled stiffness matrix becomes ill conditioned (so the infinite-precision result in Theorem 2.1 seizes to hold). Thus, it is necessary to either relocate the nodes along the interface or to collapse the nodes of the small triangle into a single node on the interface. We chose to implement the latter approach whenever the area of the small triangle was less than 1% of the area of the background grid elements. The meshes obtained are shown in Figure 10.

For this problem we designed a manufactured solution that guarantees continuity of the temperature and the heat flux across the interface, taking into account the finite jump in the conductivity coefficient. By means of a translated and rotated system of coordinates, and following the notation in Figure 9, the exact solution is given by,

[TABLE]

where the coordinates $x^{\prime}$ and $y^{\prime}$ come from a translation and rotation of the reference system defined by the following transformation,

[TABLE]

The values of conductivity and the geometric data used for the numerical computation are $k_{1}=1$, $k_{2}=5$, $x_{0}=0.12$ and $\theta=\tan^{-1}(1/0.65)$. The nonzero boundary conditions were imposed using projection-based interpolation of the manufactured solution on the boundary edges [37, 44].

Figure 11 shows the appearance of the computed ultraweak solution. As it can be observed in Figure 12, the expected convergence rates were verified once again. It is remarkable that without collapsing any nodes in these meshes, the same data points were observed for $1\leq p\leq 3$, but the last data point for $p=4$ did behave unexpectedly, so collapsing the nodes is still recommended in general.

3.4 Adaptivity

In the last example, we aimed to present a polygonal adaptive strategy. This is of interest as it has direct applications in fracture dynamics [80, 68] and topology optimization [84, 53]. Implementing such a strategy was possible, because the DPG methodology carries a natural arbitrary-order a posteriori error estimator (see (2.13) and Section 2.2). Indeed, assuming that $\eta_{K}$ is the a posteriori error estimator (representing the square of the residual as in (2.13)) for $K\in\mathcal{T}$, and $\eta_{\max}=\max_{K\in\mathcal{T}}\eta_{K}$, then the criterion used to mark an element for refinement was if $\eta_{K}\geq 0.25\eta_{\max}$ [42].

In order to refine traditional quadrilateral elements, typically hanging nodes arise in the mesh. But in practice, only one “level” of refinement is possible per element (often edges cannot have more than one hanging node), resulting in so-called quadtree meshes [83]. To implement this strategy a constrained approximation technology is necessary to handle the hanging nodes. Additionally, under anisotropic refinements, sometimes dead-lock scenarios arise (where it is logically impossible to continue refining) and these must be avoided [36]. In short, it may be challenging to implement conventional refinement strategies used for adaptivity.

An important advantage of the polygonal elements is that they naturally embrace hanging nodes, because they merely represent that a polygon has an extra edge collinear with another edge. Thus, the polygonal methods do not require an extra level of difficulty in terms of implementing the adaptive refinements. We devised a practical convex polygonal refinement strategy as illustrated in Figure 13: (a) shows the initial mesh in which an element of interest is picked and split into quadrilaterals by using the centroid and edge midpoints as depicted in (b); next, any of the resulting elements can be subsequently refined into finer quadrilaterals as shown in (c); and lastly, as shown in (d), if a neighbor element needs to be refined too, it is split into quadrilaterals assuming all adjacent collinear edges constitute a single edge (i.e. the vertices of this combined edge are used in the calculation of the centroid and its midpoint used to place the new quadrilateral node).

The manufactured solution for this problem is the sum of two Gaussian surfaces, given by the function,

[TABLE]

where the standard deviation is $\sigma=\sqrt{10^{-3}}$ and the two means are $\mu_{1}=0.25$ and $\mu_{2}=0.75$. Again, projection-based interpolation [37, 44] was used to approximate the nearly vanishing temperature boundary conditions.

In order to compare with other adaptive schemes, a traditional adaptive strategy using quadtree meshes and constrained hanging nodes via quadrilateral elements was considered here [36]. Starting with the same initial mesh, the traditional refinement strategy and the polygonal refinement strategy were allowed to refine accordingly. When using the polygonal strategy on these quadrilateral meshes, we used the more natural choice of bounding boxes instead of the bounding triangles. Additionally, the same polygonal refinement strategy was applied to an initial polygonal mesh (using bounding triangles as usual). Figure 14 shows the results of the three different scenarios after several refinements. Clearly, the traditional adaptive strategy produces quadtree meshes (see Figure 14(a)), so it is forced to refine and create new elements in areas of the domain where the solution is nearly constant. However, the polygonal adaptive strategy applied to the same initial mesh produces a more localized refinement pattern which is not a quadtree mesh (see Figure 14(b)). Lastly, the polygonal adaptive strategy applied to a polygonal mesh produces a completely nonstandard, yet localized mesh (see Figure 14(c)).

The numerical solution for $p=6$ and $\Delta p=2$ using the mesh in Figure 14(c) is presented in Figure 15. The error convergence curves corresponding to the three refinement schemes in Figure 14 are also displayed in Figure 16. The proposed polygonal refinement technique generates more edges (each new sub-segment becomes an edge) resulting in more degrees of freedom. However, in the end the additional cost is compensated by producing less elements than traditional quadtree refinement schemes (compare (b) and (c) with (a) in Figure 14). It can be seen from Figure 16 that the convergence behavior in terms of degrees of freedom is very similar using both approaches. Therefore, the polygonal adaptive strategy proposed here is competitive with the existing strategies for traditional elements, whilst being more general in its applicability as it also works for polygonal elements.

4 Conclusions

A PolyDPG method discretized with high-order polygonal elements was successfully implemented using ultraweak formulations and the DPG methodology. Here, the PolyDPG method solves Poisson’s equation. However, like with the DPG methodology, the discretization and theory is quite general. Thus, it can be applied to a large family of equations including acoustics, convection-dominated diffusion and linear elasticity. PolyDPG methods are conforming FEMs, and as with many other polytopal methods, the spaces and integration schemes are defined directly in the physical space. Indeed, given that the ultraweak formulations avoid inter-element compatibility conditions, it is relatively straightforward to obtain many of the shape functions by restricting them from a bounding (triangular or quadrilateral) element to the polygonal element. Despite the greater computational cost compared to conventional methods, the resulting PolyDPG methods are naturally high-order, carry their own residual-based a posteriori error estimator, have no need of ad hoc stabilization terms, and always produce positive-definite stiffness matrices. Moreover, under reasonable assumptions, a rigorous proof demonstrating the convergence of PolyDPG methods was included. To complement this work, the PolyDPG software [87] written in MATLAB® is provided. We hope this will prove to be a practical tool for other researchers interested in polygonal FEMs and in DPG methods.

Different illustrative examples corroborated the expected results. In the first example, $n$-sided convex polygons were investigated, while in the second example, highly distorted concave elements were examined. In both cases, as predicted by the theory, convergence rates of the form $h^{p}$ were observed for different values of $p$, confirming that PolyDPG methods are distortion-tolerant. The third example was relevant to the field of geosciences, where faults cause heterogeneity in the domain. This was simulated by irregularly cutting a uniform grid with an interface and assigning different material properties on each side. Once again, the method converged as expected, displaying its robustness in resolving heterogeneous material properties. The final example explored a polygonal adaptivity scheme driven by the arbitrary-order a posteriori error estimator of PolyDPG methods. Even though polygonal and standard refinement strategies led to practically identical convergence curves, polygonal techniques are more general since they apply to polygonal elements and avoid the typical approaches of constrained approximations via hanging nodes. These techniques may be useful in applications such as crack propagation and topology optimization.

Extension of the presented technique to arbitrary 3D polyhedral elements is in progress. In principle, the current numerical method can be extended naturally to polyhedral elements, as long as all the faces are triangular, but the case of arbitrary faces is much more challenging and might lead to analyzing nonconforming numerical methods.

Acknowledgments.

This work was partially supported with grants by ONR (N00014-15-1-2496), NSF (DMS-1418822), and AFOSR (FA9550-12-1-0484). Co-author Jaime Mora is also sponsored by a 2015 Colciencias-Fulbright scholarship (Colombia).

Appendix A Convergence

A.1 Stability and Fortin operators

Since the numerical method is technically a conforming FEM, the “standard” theory of convergence can be applied. However, the issue of numerical stability, in the sense of (2.9), must be addressed first. The DPG methodology is basically crafted to “almost” satisfy this condition, and intuitively, the larger the enriched test space $\mathscr{V}_{r}\subseteq\mathscr{V}$, the more certainty there is that the condition is satisfied. This translates to increasing $\Delta p_{K}$ for all $K\in\mathcal{T}$ in (2.20), so that $\mathscr{V}_{r}$ becomes larger. Note that this increases the local (element-wise) computational burden. In practice, the numerical stability is observed even with very modest values of $\Delta p$.

However, to have a rigorous result, it is necessary to establish (2.9) theoretically. To do so, it is helpful to consider a linear and continuous Fortin operator, $\Pi_{F}:\mathscr{V}\to\mathscr{V}_{r}$, satisfying the orthogonality condition, $b(\mathfrak{u}_{h},\mathfrak{v}-\Pi_{F}\mathfrak{v})=0$, for all $\mathfrak{u}_{h}\in\mathscr{U}_{h}$ and $\mathfrak{v}\in\mathscr{V}$. If it exists, it follows that [58],

[TABLE]

where $M_{F}\geq\|\Pi_{F}\|=\sup_{\mathfrak{v}\in\mathscr{V}}\frac{\|\Pi_{F}\mathfrak{v}\|_{\mathscr{V}}}{\|\mathfrak{v}\|_{\mathscr{V}}}$, $\|b\|=\sup_{(\mathfrak{u},\mathfrak{v})\in\mathscr{U}\!\times\!\mathscr{V}}\frac{|b(\mathfrak{u},\mathfrak{v})|}{\|\mathfrak{u}\|_{\mathscr{U}}\|\mathfrak{v}\|_{\mathscr{V}}}$ and $\gamma=\inf_{\mathfrak{u}\in\mathscr{U}}\sup_{\mathfrak{v}\in\mathscr{V}}\frac{|b(\mathfrak{u},\mathfrak{v})|}{\|\mathfrak{u}\|_{\mathscr{U}}\|\mathfrak{v}\|_{\mathscr{V}}}$, where the infima and suprema are tacitly assumed to be taken over nonzero elements. Note that when $\mathscr{V}$ is a broken test space, as in this case, the Fortin operator can be separately constructed locally at each element $K\in\mathcal{T}$. Constructions of such Fortin operators do exist for triangles [21], but have not been constructed for other shapes yet. Nevertheless, numerical results show it is reasonable to expect them to exist, and this will be assumed in what follows. In any case, note that Fortin operators merely yield a conservative estimate, but in practice the results are better (i.e. instead of $M_{F}$, there is a moderate constant, $\mathcal{O}(1)$-$\mathcal{O}(10)$, multiplying $\frac{\|b\|}{\gamma}$ in (A.1)).

A.2 Fractional spaces

For a given $s\geq 0$ and any Lipschitz domain $A\subseteq\mathbb{R}^{2}$, the fractional Sobolev spaces, $H^{1+s}(A)$, $\boldsymbol{H}^{s}(\operatorname{div},A)$ and $H^{s}(A)$, are slightly smoother subspaces of $H^{1}(A)$, $\boldsymbol{H}(\operatorname{div},A)=\boldsymbol{H}^{0}(\operatorname{div},A)$ and $L^{2}(A)=H^{0}(A)$, respectively. As an obvious placeholder, let $\mathscr{W}(A)=\mathscr{W}^{0}(A)$ be one of these spaces, and for $s\geq 0$, let $\mathscr{W}^{s}(A)$ be its slightly smoother fractional counterpart. As one might expect, $\|w\|_{\mathscr{W}(A)}\leq\|w\|_{\mathscr{W}^{s_{1}}(A)}\leq\|w\|_{\mathscr{W}^{s_{2}}(A)}$ for $0\leq s_{1}\leq s_{2}$ and all $w\in\mathscr{W}^{s_{2}}(A)\subseteq\mathscr{W}^{s_{1}}(A)\subseteq\mathscr{W}(A)$. For more details on these spaces and their norms, see [70].

Using interpolation theory (see [70, Appendix B]) applied to the universal extension operators of Sobolev spaces of differential forms defined in [60] (which is even more general than the universal extension operator defined by Stein in [81]), it is possible to establish the existence of a continuous extension operator,

[TABLE]

where $s\geq 0$, $\Omega$ is the domain where the equations are being solved, and $C_{E}=C_{E}(s,\Omega)>0$.

The fractional skeleton spaces are better defined directly through fractional trace operators of Lipschitz elements $K\in\mathcal{T}$ (see [70]) as H^{\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}+s}(\partial K)=\mathrm{tr}_{H^{1+s}(K)}(H^{1+s}(K)) and H^{-\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}+s}(\partial K)=\mathrm{tr}_{\boldsymbol{H}^{s}(\operatorname{div},K)}(\boldsymbol{H}^{s}(\operatorname{div},K)) (see (2.18) for the explicit trace operators for the $s=0$ case). Again, using placeholders these are written as $\mathscr{W}^{s}(\partial K)=\{\hat{w}_{K}=\mathrm{tr}_{\mathscr{W}^{s}(K)}w\mid w\in\mathscr{W}^{s}(K)\}$, so that their minimum energy extension norm is

[TABLE]

At the global level, define the global trace operators as

[TABLE]

Note that $H_{0}^{1+s}(\Omega)=\mkern 1.5mu\overline{\mkern-1.5muC_{0}^{\infty}(\Omega)\mkern-1.5mu}\mkern 1.5mu{}^{\|\cdot\|_{H^{1+s}(\Omega)}}$, so that the global fractional skeleton spaces are (see (2.5) for the $s=0$ case),

[TABLE]

Analogous to (2.6), the fractional trial subspace for $s\geq 0$ is

[TABLE]

and it is easy to see $\|\mathfrak{u}\|_{\mathscr{U}}\leq\|\mathfrak{u}\|_{\mathscr{U}^{s_{1}}}\leq\|\mathfrak{u}\|_{\mathscr{U}^{s_{2}}}$ for $0\leq s_{1}\leq s_{2}$ and all $\mathfrak{u}\in\mathscr{U}^{s_{2}}\subseteq\mathscr{U}^{s_{1}}\subseteq\mathscr{U}$.

A.3 Approximation properties

Next, for a bounded Lipschitz domain $A\subseteq\mathbb{R}^{2}$ and polynomial order $p\in\mathbb{N}$, consider commuting exact sequence discretizations for $H^{1}(A)$, $\boldsymbol{H}(\operatorname{div},A)$ and $L^{2}(A)$, such that

[TABLE]

More abstractly, the discretizations are written as $\mathscr{W}_{hp}(A)\subseteq\mathscr{W}(A)$. Then, given the polynomials contained in the discretizations $\mathscr{W}_{hp}(A)$, it is well known that for $p\geq s>0$, there exists a constant $C_{h}=C_{h}(A,s)>0$, such that for all $w\in\mathscr{W}^{s}(A)$,

[TABLE]

For each $K\in\mathcal{T}$, the local trace spaces are supposed to be $\mathscr{W}_{hp}(\partial K)=\mathrm{tr}_{\mathscr{W}^{0}(K)}(\mathscr{W}_{hp}(K))$ for some $\mathscr{W}_{hp}(K)$.

We would like these approximation properties to hold for our choices of discrete trial spaces, $\mathscr{U}_{h}$ in (2.19), and this is indeed the case. The first two components of $\mathscr{U}_{h}$ when restricted to $K\in\mathcal{T}$ (representing $L_{hp}^{2}(K)$ and $(L_{hp}^{2}(K))^{2}$) are restrictions to $K$ of $\mathcal{P}^{p-1}(T_{K})$ and $(\mathcal{P}^{p-1}(T_{K}))^{2}$, so they do trivially contain $\mathcal{P}^{p-1}(K)$ and $(\mathcal{P}^{p-1}(K))^{2}$ respectively. This means that (A.8) holds for those two spaces, but as we will see soon, it suffices (and is preferable) to have this result for the bounding triangle $T_{K}$ (which is obviously true). For the third and fourth components of $\mathscr{U}_{h}$, representing the skeleton variables, locally at each $K\in\mathcal{T}$ it suffices to show that $\mathcal{P}^{p}_{C}(\partial K)=\mathrm{tr}_{H^{1}(K)}(H_{hp}^{1}(K))$ and $\mathcal{P}^{p-1}(\partial K)=\mathrm{tr}_{\boldsymbol{H}(\operatorname{div},K)}(\boldsymbol{H}_{hp}(\operatorname{div},K))$ for some $H_{hp}^{1}(K)$ and $\boldsymbol{H}_{hp}(\operatorname{div},K)$ satisfying the properties in (A.7), where $\mathcal{P}^{p}_{C}(\partial K)$ and $\mathcal{P}^{p-1}(\partial K)$ are defined in (2.18). For this, consider the shape-regular edge-compatible triangulations of each $K\in\mathcal{T}$, denoted by $\mathscr{T}(K)=\{\mathscr{T}_{i}(K)\}_{i\in I_{K}}$ (with $I_{K}$ finite), and define the spaces,

[TABLE]

It can easily be checked that $\hat{u}_{K}=u|_{\partial K}\in\mathcal{P}^{p}_{C}(\partial K)$ for all $u\in H_{hp}^{1}(K)$ and $(\hat{q}_{\hat{\mathbf{n}}})_{K}=\boldsymbol{q}|_{\partial K}\!\cdot\!\hat{\mathbf{n}}_{K}\in\mathcal{P}^{p-1}(\partial K)$ for all $\boldsymbol{q}\in\boldsymbol{H}_{hp}(\operatorname{div},K)$, and that these inclusions are surjective. Thus, $\mathcal{P}^{p}_{C}(\partial K)=\mathrm{tr}_{H^{1}(K)}(H_{hp}^{1}(K))$ and $\mathcal{P}^{p-1}(\partial K)=\mathrm{tr}_{\boldsymbol{H}(\operatorname{div},K)}(\boldsymbol{H}_{hp}(\operatorname{div},K))$ as desired. This implies that (A.8) also holds for $H_{hp}^{1}(K)$ and $\boldsymbol{H}_{hp}(\operatorname{div},K)$, which are closely related to the skeleton discretizations of $\mathscr{U}_{h}$.

A.4 Interpolation estimates

The idea is to define a bounded linear interpolation operator $\Pi_{\mathscr{U}^{s}}:\mathscr{U}^{s}\to\mathscr{U}_{h}$ such that $\Pi_{\mathscr{U}^{s}}\mathfrak{u}_{h}=\mathfrak{u}_{h}$ for every $\mathfrak{u}_{h}\in\mathscr{U}_{h}$ and $s>\frac{1}{2}$. Typically this implies constructing interpolation operators for every component of $\mathscr{U}$. Moreover, for each component this construction is done locally at every $K\in\mathcal{T}$ in such a way that the inter-element compatibility properties are satisfied.

The first two components of $\mathscr{U}$ are $L^{2}(\Omega)$ and $\boldsymbol{L}^{2}(\Omega)$, which are effectively three $L^{2}(\Omega)$ components. The last two skeleton components are H_{0}^{\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}+s}(\partial\mathcal{T}) and H^{-\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}+s}(\partial\mathcal{T}). The discretizations of these three spaces are (see (2.19)),

[TABLE]

where the definitions of $H_{hp}^{1}(K)$ and $\boldsymbol{H}_{hp}(\operatorname{div},K)$ are in (A.9). Thus, it suffices to construct,

[TABLE]

meaning that we must define $\Pi_{H^{s}(K)}$, \Pi_{H^{\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}+s}(\partial K)} and \Pi_{H^{-\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}+s}(\partial K)}.

The operator $\Pi_{H^{s}(K)}$ can be chosen as the $L^{2}(K)$-projection to $\mathcal{P}^{p-1}(K)$ directly on $K$ (so $\Pi_{H^{s}(K)}\delta u=\delta u$ for all $\delta u\in\mathcal{P}^{p-1}(K)$). Consider now a simple scaling by $h_{K}=\operatorname{diam}(K)$, so that $\hat{K}$ has $\operatorname{diam}(\hat{K})=1$. Using (A.8) for $p\geq s>\frac{1}{2}$ results in the abstract expression,

[TABLE]

for any $w\in\mathscr{W}^{s}(\hat{K})$, where $C_{\hat{K}}=C_{\hat{K}}(\hat{K},p,s)>0$. Scaling appropriately then yields for any $w\in\mathscr{W}^{s}(K)$,

[TABLE]

The issue with this estimate is that it depends on the element shape $K$ (via $\hat{K}$), so it is inconvenient as it may become much larger with mesh refinements. The solution is to use the bounding triangle and the extension operator defined in (A.2), so that (as in [20]) the interpolation operator is defined for any $w\in\mathscr{W}^{s}(\Omega)$ as,

[TABLE]

where $\Pi_{\mathscr{W}^{s}(T_{K})}$ is the $L^{2}(T_{K})$-projection. Scaling and rotating transforms the bounding triangle $T_{K}$ to a unique triangle $\hat{T}_{0}$ (independent of the element $K$) with $\operatorname{diam}(\hat{T}_{0})=1$. This means $T_{K}$ is scaled by $h_{T_{K}}=\operatorname{diam}(T_{K})=\frac{6}{\sqrt{3}}r_{\max}\leq\sqrt{12}h_{K}$, where $r_{\max}$ is the distance of the centroid to the furthest vertex and $h_{K}=\operatorname{diam}(K)$. Using the same reasoning gives,

[TABLE]

for every $w\in\mathscr{W}^{s}(K)$, where $C_{\hat{T}_{0}}=C_{\hat{T}_{0}}(p,s)>0$ is now independent of $K$.

Next, consider the skeleton variables for an element $K\in\mathcal{T}$ and its respective shape-regular and edge-compatible triangulation denoted by $\mathscr{T}(K)=\{\mathscr{T}_{i}(K)\}_{i\in I_{K}}$. The theory of projection-based interpolation [37] implies that for any polygonal domain $A\subseteq\mathbb{R}^{2}$ and for any $s>\frac{1}{2}$ there exist commuting operators,

[TABLE]

Thus, for any $p\geq s>\frac{1}{2}$ and triangle $\mathscr{T}_{i}(K)$ (so $\operatorname{diam}(\mathscr{T}_{i}(K))\leq\operatorname{diam}(K)=h_{K}$), the result in (A.13) applies and yields

[TABLE]

where the $K$-independent $C_{\hat{\mathscr{T}}_{0}}=C_{\hat{\mathscr{T}}_{0}}(p,s)>0$ exists due to the assumed uniform shape-regularity of the $\mathscr{T}_{i}(K)$ (across all $K\in\mathcal{T}$ and all meshes being considered). Adding among $\mathscr{T}(K)$ is valid due to the compatibility of the projection-based interpolation in the triangulation, so that

[TABLE]

Lastly, consider the well-defined trace interpolation,

[TABLE]

so that (see (A.3)),

[TABLE]

This is true for every $w\in\mathrm{tr}_{\mathscr{W}^{s}(K)}^{-1}\{\hat{w}_{K}\}$, so take the infimum to yield

[TABLE]

Putting everything together and generalizing for any $p\in\mathbb{N}$ and $s>\frac{1}{2}$, gives

[TABLE]

where the constants $C_{\hat{T}_{0}}$, $C_{H^{1+s}(\hat{\mathscr{T}}_{0})}$ and $C_{\boldsymbol{H}^{s}(\operatorname{div},\hat{\mathscr{T}}_{0})}$ only depend on $p$ and $s$, but not on $K$ (the last two constants depend on the uniform shape-regularity of the edge-compatible triangulations of all elements). Finally, since these constants come from triangles, the theory of projection-based interpolation [37] implies that in the $p$-asymptotic limit,

[TABLE]

where $\widetilde{C}_{\hat{T}_{0}}$, $\widetilde{C}_{H^{1+s}(\hat{\mathscr{T}}_{0})}$ and $\widetilde{C}_{\boldsymbol{H}^{s}(\operatorname{div},\hat{\mathscr{T}}_{0})}$ are constants independent of $p$ and of any $K\in\mathcal{T}$ across all possible meshes being considered.

A.5 Final convergence estimates

Use the global interpolation operators in (A.11) to construct the bounded linear global interpolation operator $\Pi_{\mathscr{U}^{s}}:\mathscr{U}^{s}\to\mathscr{U}_{h}$. Note that adding (A.22) associated with $u\in H^{s}(\Omega)$ among $K\in\mathcal{T}$, using the robust finite overlap condition, and the extension operator in (A.2), gives:

[TABLE]

where $h=\sup_{K\in\mathcal{T}}h_{K}$ and $C_{E}=C_{E}(s,\Omega)$ is not dependent on $p$. The same estimate holds for the variable $\boldsymbol{q}\in(H^{s}(\Omega))^{2}$, and similar bounds (even without using extension operators and the finite overlap condition) hold for \hat{u}\in H_{0}^{\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}+s}(\partial\mathcal{T}) and \hat{q}_{\hat{\mathbf{n}}}\in H^{-\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}+s}(\partial\mathcal{T}). Then, assume $M_{F}$ is independent of the family of meshes being considered, and choose the interpolant in (A.1) along with the estimates of the type in (A.24), so that

[TABLE]

where $C=C(p,s,\Omega)>0$, but is independent of the meshes being considered. Moreover, if $M_{F}$ is $p$-independent, then in the $p$-asymptotic limit, the following $hp$-convergence estimate holds (see (A.23)),

[TABLE]

where $\widetilde{C}=\widetilde{C}(s,\Omega)$ is independent of $p$. This concludes the results summarized in Theorem 2.1.

Remark A.1**.**

Starting directly from the quasi-optimal error estimate in (A.1), and avoiding interpolation, it is possible to get a better estimate for the variable \hat{u}\in H_{0}^{\mathchoice{\raisebox{0.0pt}{\displaystyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\textstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}{\raisebox{0.0pt}{\scriptscriptstyle{}\mathchoice{\raisebox{-0.2pt}{$\displaystyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\textstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptstyle{}^{1}\!$}}{\raisebox{-0.2pt}{$\scriptscriptstyle{}^{1}\!$}}/\mathchoice{\raisebox{-0.1pt}{$\displaystyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\textstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptstyle{}_{\!2}$}}{\raisebox{-0.1pt}{$\scriptscriptstyle{}_{\!2}$}}}}+s}(\partial\mathcal{T}) by using the results in [6], provided the triangulations $\mathscr{T}(K)$ are quasi-uniform across all $K\in\mathcal{T}$ and all meshes being considered. In that case, the $hp$-convergence estimate in (A.26) will have a $\ln p$ instead of $(\ln p)^{2}$.