The harmonic virtual element method: stabilization and exponential convergence for the Laplace problem on polygonal domains

TL;DR

The paper introduces the harmonic virtual element method (harmonic VEM) for 2D Laplace problems on polygonal domains, achieving exponential convergence rates superior to traditional hp finite element and virtual element methods.

Contribution

It develops a boundary-only degrees of freedom harmonic VEM, stabilizes it, and proves exponential convergence, outperforming existing hp methods.

Findings

Exponential convergence rate of order exp(-b√N) achieved.

Harmonic VEM uses only boundary degrees of freedom.

Method outperforms hp FEM and VEM in convergence speed.

Abstract

We introduce the harmonic virtual element method (harmonic VEM), a modification of the virtual element method (VEM) for the approximation of the 2D Laplace equation using polygonal meshes. The main difference between the harmonic VEM and the VEM is that in the former method only boundary degrees of freedom are employed. Such degrees of freedom suffice for the construction of a proper energy projector on (piecewise harmonic) polynomial spaces. The harmonic VEM can also be regarded as an "$H^1$-conformisation" of the Trefftz discontinuous Galerkin-finite element method (TDG-FEM). We address the stabilization of the proposed method and develop an $hp$ version of harmonic VEM for the Laplace equation on polygonal domains. As in Trefftz DG-FEM, the asymptotic convergence rate of harmonic VEM is exponential and reaches order $\mathcal O ( \exp(-b\sqrt[2]{N}))$, where $N$ is the number of degrees of freedom. This result overperformes its counterparts in the framework of $hp$ FEM and $hp$ VEM, where the asymptotic rate of convergence is of order $\mathcal O ( \exp(-b\sqrt[3]{N}) )$.

Equations

$$D^{\mathbf{\alpha }}u={\partial }^{{\alpha }_1,{\alpha }_2}u,\text{where }\mathbf{\alpha }=({\alpha }_1,{\alpha }_2)\in {\mathbb{N}}_0^2.$$

$$|D^{k}u{|}^2=\sum _{\mathbf{\alpha }\in {\mathbb{N}}_0^2,|\mathbf{\alpha }|=k}|D^{\mathbf{\alpha }}u{|}^2.$$

$${\Phi }_{\mathbf{\beta }}(\text{x}):=\prod _{i=1}^{{\mathcal{V}}_{\Omega }}\min (1,|\text{x}-{\text{A}}_i|{)}^{{\beta }_i},$$

$$\Vert u{\Vert }_{H_{\mathbf{\beta }}^{m,\ell }(\Omega )}^2:=\Vert u{\Vert }_{\ell -1,\Omega }^2+|u{|}_{H_{\beta }^{m,\ell }(\Omega )}^2:=\Vert u{\Vert }_{\ell -1,\Omega }^2+\sum _{k=\ell }^m\Vert {\Phi }_{\mathbf{\beta }+k-\ell }|D^{k}u|{\Vert }_{0,\Omega }^2.$$

$$\mathbf{\beta }+k-\ell \in {\mathbb{R}}^{{\mathcal{V}}_{\Omega }},(\mathbf{\beta }+k-\ell {)}_i={\mathbf{\beta }}_i+k-\ell ,i=1,\dots ,{\mathcal{V}}_{\Omega }.$$

$$\begin{array}{rl}& {\mathcal{B}}_{\mathbf{\beta }}^{\ell }(\Omega ):=\left\{u\in H_{\mathbf{\beta }}^{m,\ell }(\Omega )\forall m\ge \ell \ge 0\text{ with }\Vert {\Phi }_{\mathbf{\beta }+k-\ell }|D^{k}u|{\Vert }_{0,\Omega }\le c_u{d}_u^{k-\ell }(k-\ell )!\forall k\in {\mathbb{N}}_0,k\ge \ell \right\},\\ & {\text{cacca}}_{\mathbf{\beta }}^2(\Omega ):=\left\{u\in H_{\mathbf{\beta }}^{m,2}(\Omega )\forall m\ge 2\text{ with }|D^{k}u(\text{x})|\le c_u{d}_u^{k}k!{\Phi }_{\mathbf{\beta }+k-1}^{-1}(\text{x})\forall k\in {\mathbb{N}}_0,\forall \text{x}\in \overline{\Omega }\right\},\end{array}$$

$$|D^{\mathbf{\alpha }}u({\text{x}}_0)|\le \left|D^{k}u({\text{x}}_0)\right|\le c_u\frac{d_u^{|\mathbf{\alpha }|}}{{\Phi }_k({\text{x}}_0)}|\mathbf{\alpha }|!\forall {\text{x}}_0\in \overline{\Omega },$$

$$\mathcal{N}(u):=\bigcup_{\textbf{x}_{0}\in\overline{\Omega};\,\textbf{x}_{0}\neq\textbf{A}_{i},\,i=1,\dots,\mathcal{V}_{\Omega}}\left\{\textbf{x}\in\mathbb{R}^{2}\,\Big{|}\,|\textbf{x}-\textbf{x}_{0}|<c\frac{\Phi(\textbf{x}_{0})}{d_{u}},\;\forall\,c\in\left(0,\frac{1}{2}\right)\right\}.$$

$$\sum _{\mathbf{\alpha }\in {\mathbb{N}}_0^2}\frac{D^{\mathbf{\alpha }}u({\text{x}}_0)}{\mathbf{\alpha }!}(\text{x}-{\text{x}}_0{)}^{\mathbf{\alpha }},{\text{x}}_0\in \overline{\Omega };{\text{x}}_0\ne {\text{A}}_i,i=1,\dots ,{\mathcal{V}}_{\Omega },$$

$$\sum _{k\in {\mathbb{N}}_0}\sum _{|\mathbf{\alpha }|=k}\frac{|D^{\mathbf{\alpha }}u({\text{x}}_0)|}{\mathbf{\alpha }!}|\text{x}-{\text{x}}_0{|}^{|\mathbf{\alpha }|}\le \overline{c}<\infty ,\forall \text{x}\in B\left({\text{x}}_0,c\frac{\Phi ({\text{x}}_0)}{d_u}\right),c\in \left(0,\frac{1}{2}\right),$$

$$\begin{array}{rl}& \sum _{k\in {\mathbb{N}}_0}\sum _{|\mathbf{\alpha }|=k}\frac{|D^{\mathbf{\alpha }}u({\text{x}}_0)}{\mathbf{\alpha }!}|\text{x}-{\text{x}}_0{|}^k\le \sum _{k\in {\mathbb{N}}_0}\sum _{|\mathbf{\alpha }|=k}\frac{1}{\mathbf{\alpha }!}{c}_u\frac{d_u^k}{{\Phi }_k({\text{x}}_0)}|\mathbf{\alpha }|!c^k\frac{{\Phi }_k({\text{x}}_0)}{d_u^k({\text{x}}_0)}\\ & =c_u\sum _{k\in {\mathbb{N}}_0}\sum _{|\mathbf{\alpha }|=k}\frac{|\mathbf{\alpha }|!}{\mathbf{\alpha }!}{c}^k=c_u\sum _{k\in {\mathbb{N}}_0}\sum _{\ell =0}^k\left(\genfrac{}{}{0ex}{}{k}{\ell }\right)c^k=c_u\sum _{k\in {\mathbb{N}}_0}{\left(2c\right)}^k\le \overline{c}<+\infty ,\end{array}$$

$$\{\begin{array}{ll}\Delta u=0& \text{in }\\ u=g& \text{on }\partial \Omega .\end{array}$$

$$\{\begin{array}{l}\text{find }u\in V_g\text{ such that}\\ a(u,v)=0\forall v\in V_0,\end{array}$$

$$\begin{array}{rl}& a(u,v)=(\nabla u,\nabla v{)}_{0,\Omega },V_{\tilde{g}}:=H_{\tilde{g}}^1(\Omega ):=\left\{v\in H^1(\Omega )\mid v_{{|}_{\partial K}}=\tilde{g}\right\}\text{for some }\tilde{g}\in {\mathcal{B}}_{\beta }^{\frac{3}{2}}(\partial \Omega ).\end{array}$$

$$\{\begin{array}{l}0\text{ is a vertex of }\Omega ,\\ u\text{, solution to eqrefweak continuous problem, has only a singularity, precisely at }0.\end{array}$$

$$\mathbb{B}(\partial K):=\left\{v\in {\mathcal{C}}^0(\partial K)\mid v{}_{{|}_s}\in {\mathbb{P}}_p(s),\forall s\in {\mathcal{E}}^K\right\},$$

$$V^{\Delta }(K):=\left\{v\in H^1(K)\mid \Delta v=0,v{}_{{|}_{\partial K}}\in \mathbb{B}(\partial K)\right\}.$$

$${\text{dof}}_i({\phi }_j)={\delta }_{i,j},i,j=1,\dots ,\dim (V^{\Delta }(K)),$$

$$V_n:=\left\{v_n\in {\mathcal{C}}^0(\overline{\Omega })\mid v_n{}_{{|}_K}\in V^{\Delta }(K),\forall K\in {\mathcal{T}}_n\right\},$$

$$V_{n,0}:=\left\{v_n\in V_n\mid v_n{}_{{|}_{\partial \Omega }}=0\right\}$$

$$V_{n,g}:=\left\{v_n\in V_n\mid v_n{}_{{|}_s}=g_{GL}^s\forall s\in {\mathcal{E}}_n^b\right\}.$$

$$\{\begin{array}{l}\text{find }{u}_n\in V_{n,g}\text{ such that}\\ a_n(u_n,v_n)=0\forall v_n\in V_{n,0},\end{array}$$

$$a(u,v)=\sum _{K\in {\mathcal{T}}_n}{a}^K(u_{{|}_K},v_{{|}_K}),a^K(u_{{|}_K},v_{{|}_K}):={\int }_K\nabla u\cdot \nabla vdx.$$

$$a_n(u_n,v_n)=\sum _{K\in {\mathcal{T}}_n}{a}_n^K(u_n{}_{{|}_K},v_n{}_{{|}_K})\forall u_n,v_n\in V_n$$

$$a^K(q,v)=a_n^K(q,v)\forall q\in {\mathbb{H}}_p(K),\forall v\in V^{\Delta }(K),$$

$${\alpha }_{\ast }(p)|v{|}_{1,K}^2\le a_n^K(v,v)\le {\alpha }^{\ast }(p)|v{|}_{1,K}^2\forall v\in V^{\Delta }(K),$$

$$\alpha (p):=\frac{1+{\alpha }^{\ast }(p)}{{\alpha }_{\ast }(p)},$$

$$|v{|}_{1,{\mathcal{T}}_n}^2:=\sum _{K\in {\mathcal{T}}_n}|v{|}_{1,K}^2\forall v\in L^2(\Omega )\text{such that }v{}_{{|}_K}\in H^1(K)\forall K\in {\mathcal{T}}_n.$$

$$|u-u_{n}|_{1,\Omega}\leq\alpha(p)\Big{\{}|u-u_{\pi}|_{1,\mathcal{T}_{n}}+|u-u_{I}|_{1,\Omega}\Big{\}}\qquad\forall u_{\pi}\in S^{p,\Delta}(\Omega,\mathcal{T}_{n}),\quad\forall u_{I}\in V_{n,g},$$

$$|u-u_n{|}_{1,\Omega }\le |u-u_I{|}_{1,\Omega }+|u_I-u_n{|}_{1,\Omega }\forall u_{\pi }\in S^{p,\Delta }(\Omega ,{\mathcal{T}}_n),\forall u_I\in V_{n,g}.$$

Full text

1 Introduction

In this work, we deal with the approximation of the Laplace equation on polygonal domains based on a novel method, whose main advantage is the fact of having a very small number of degrees of freedom if compared to standard finite element methods (FEM). This is not of course the first attempt to approximate a Laplace problem with methods based on approximation spaces having small dimension. Among the other methods available, we limit ourselves to recall only two of them.

The first one is the boundary element method (BEM) [32]. BE spaces consist of functions defined only on the boundary of the computational domain. Clearly, the BE space on a boundary mesh of characteristic mesh size $h$ contains many less degrees of freedom than the corresponding FE space on a volume mesh having the same characteristic mesh size $h$. This comes at a price of a fully populated matrix in the resulting system of linear equations, expensive quadrature rules needed for evaluation of matrix entries and expensive numerical reconstruction of the solution in the interior of the computational domain. These difficulties can be partially alleviated by using advanced fast boundary element methods (see e.g. [32] and references therein), that usually results in nontrivial algorithms that are not easy to implement.

A second (and more recent) approach is given by the so-called Trefftz discontinuous Galerkin-FEM (TDG-FEM), which was introduced in [22, 23] and was generalized to its $hp$ version in [21] (we recall that an $hp$ Galerkin method is a method where the convergence of the error is achieved by a proper combination of mesh refinements and an increase of the local degree of accuracy and thereby of the dimension of local spaces). TDG-FE spaces consist of piecewise harmonic polynomials over a decomposition of the computational domain into triangles and quadrilaterals. As a consequence, the resulting method has a DG structure, since the dimension of harmonic polynomial spaces is not large enough for enforcing global continuity of the discretization space. We also point out that in [21] it was provided a result concerning $hp$ approximation of harmonic functions by means of harmonic polynomials, following the ideas of the pioneering works [28, 26].

The advantage of TDG-FEM with respect to standard FEM is that the dimension of local spaces considerably reduces still keeping the optimal rate of convergence of the error. More precisely, for a fixed local polynomial degree $p$, the dimension of the local TDG-FE space is equal to $2p+1\approx 2p$, whereas the dimension of local FE spaces is $\frac{(p+1)(p+2)}{2}\approx\frac{p^{2}}{2}$. This advantage is possible since the degrees of freedom that are removed in TDG-FEM are superfluous for the approximation of a Laplace equation. We emphasize that employing piecewise harmonic polynomials leads inevitably to a discontinuous method, which is therefore not anymore $H^{1}$-conforming.

The approach in [21] can be generalized easily to polygonal TDG-FEM, following e.g. [3]. Polygonal methods received an outstanding interest in the last decade by the scientific community due to the high flexibility in dealing with nonstandard geometries. Among the others, we mention the following methods: hybrid high–order methods [18], mimetic finite differences [10], hybrid DG-methods [16], polygonal FEM [34, 29, 20], polygonal DG-FEM [14], BEM-based FEM [31].

The virtual element method (VEM) is an alternative approach enabling computation of polygonal (polyhedral in 3D) meshes [6, 7]. It is based on globally continuous discretization spaces that generally consist locally of Trefftz-like functions. More precisely, the key idea of VEM is that trial and test spaces consists of functions that are solutions to local PDE problems in each element. Since these local problems do not admit closed form solutions, the bilinear form, and thereby the entries of the stiffness matrix, are not computable in general. The computable version involves an approximate discrete bilinear form consisting of two additive parts: one that involves local projections on polynomial spaces and a computable stabilizing bilinear form. We emphasize that the approximate discrete bilinear form can be evaluated without explicit knowledge of local basis functions in the interior of the polygonal element: an indirect description via the associated set of internal degrees of freedom suffices.

In [8], the $hp$ version of VEM for the Poisson problem with quasi-uniform meshes and constant polynomial degree was studied, whereas, in [9], the $hp$ version of VEM for the approximation of corner singularities was discussed. Besides, a multigrid algorithm for the pure $p$ version of VEM was investigated in [4]. Also, a study regarding the condition number of the stiffness matrix for the $hp$ version of 2D and 3D VEM is the topic of [24] and [17], respectively.

The aim of the present work consists in modifying the $hp$ virtual element space of [9], trying at the same time to mimick the “harmonic” approach of TDG-FEM. The arising method, which goes under the name of harmonic VEM, makes use only of boundary degrees of freedom (the internal degrees of freedom of the standard VEM can be omitted). More precisely, functions in the harmonic virtual element space are harmonic reconstructions of piecewise continuous polynomial traces over the boundary of the polygons in the polygonal decomposition of the computational domain. It is immediate to check that the associated space contains (globally discontinuous) piecewise harmonic polynomials.

The stiffness matrix is not computed exactly on the harmonic virtual element space. Its construction is based on two ingredients: a local energy projector on the space of harmonic polynomials and a stabilizing bilinear form, which only approximates the continuous one and which is computable on the complete space. As in standard VEM, the projectors and stabilizing bilinear forms are computed only by means of the degrees of freedom, without the need of knowing trial and test functions in the interior of individual elements explicitly. Importantly, the implementation of the harmonic VEM does not require two-dimensional quadrature formulas.

The main result of the paper states that, similarly to the $hp$ version of TDG-FEM, the asymptotic convergence rate for the energy error is proportional to $\exp(-b\sqrt[2]{N})$, where $N$ is the dimension of the global discretization space. This result is an improvement of the analogous statement in the framework of the $hp$ FEM [33] and $hp$ VEM [9], where the rate of decay of the error is proportional to $\exp(-b\sqrt[3]{N})$. As a byproduct of the main result we prove in Sections 3.1 and 3.2 novel stabilization estimates that are much sharper than in the general $hp$-VEM [9] and that are interesting on their own.

We state the difference between the two approaches, namely the TDG-FEM [21] and the harmonic VEM. The TDG-FEM is a discontinuous method, but local spaces are made of explicitly known functions, i.e. (harmonic) polynomials; besides, only internal degrees of freedom on each element are considered. The harmonic VEM is a $H^{1}$-conforming method which only employs boundary degrees of freedom; the basis functions are not known explicitly, but the stiffness matrix can be efficiently built employing only the degrees of freedom. Both methods are characterized by the fact that the space of piecewise harmonic polynomials is contained in the local approximation spaces; in fact, the TDG-FE space is the space of (globally discontinuous) piecewise harmonic polynomials, while the harmonic virtual element space is richer, in general.

We emphasize that the formulation of the $hp$ harmonic VEM presents some improvements with respect to the standard $hp$ VEM [9]. Less assumptions on the geometry of the polygonal decomposition are required, better bounds for the stabilization are presented and a very tidy result, concerning approximation by functions in the harmonic virtual element space, is proven.

In this paper, we only investigate in details the $hp$ version of harmonic VEM, that is the method of choice for an efficient approximation of corner singularities. A modification of Section 4.3, along with the arguments in [28], leads to $h$ approximation results. For the $p$ version of harmonic VEM, instead, one has to deal with two issues that lie outside the scope of the present paper.

The first one is the pollution effect due to the stabilization of the method, which is typical also of the $p$ version of VEM, see Lemmata 3.1 and 3.3, which in fact can be overcome at the price of having a stabilization challenging to be computed, as explained in Section 3.2.

The second one is that the $p$ approximation estimates by harmonic functions depend on the shape of the domain of approximation via the so-called “exterior cone condition”, see [27, Theorem 2]. These matters introduce additional technicalities which will not be addressed in this paper.

It is worth mentioning that the theoretical framework of [28, 30] for the analysis of methods based on harmonic polynomials can be seen as an intermediate step towards more gruelling challenges. More precisely, one can use the so-called Vekua theory [36] to shift the results related to the approximation of harmonic functions (i.e. for functions belonging to the kernel of the Laplace operator) by means of harmonic polynomials, to the results related to the approximation of functions in the kernel of more general differential operators by means of generalized harmonic polynomials. A very important example is provided by the approximation of functions in the kernel of the Helmholtz operator; this was investigated in deep in the framework of partition of unity methods (PUM) [28] and TDG-FEM [30]. An extension of the harmonic VEM to Trefftz-like VEM for the Helmholtz equation has been recently investigated in [25].

The outline of the paper is the following. In Section 2, we present the model problem and we recall some regularity properties of its solution. In Section 3, we introduce the harmonic VEM; in particular, we discuss the contruction of the stiffness matrix along with the construction of a proper stabilization of the method and of an energy projector from local harmonic virtual element spaces into spaces of (piecewise) harmonic polynomials. After having defined the concept of “$hp$ graded polygonal meshes”, we prove approximation estimates by harmonic polynomials and functions in the harmonic virtual element space in Section 4. This approximation scheme results in exponential convergence of the energy error with respect to the total number of degrees of freedom. Numerical tests validating the theoretical results, together with a numerical comparison between the performances of the $hp$ harmonic VEM and the $hp$ VEM, are shown in Section 5.

Throughout the paper, we write $f\lesssim g$ for two positive quantities $f$ and $g$ depending on a discretization parameter (typically $h$ or $p$) if there exists a parameter-independent positive constant $c$ such that $f\leq cg$ holds for all values of the parameter. We write $f\approx g$ if $f\lesssim g$ and $g\lesssim f$ holds.

We denote by $\mathbb{P}_{p}(D)$ the spaces of polynomials of degree $p\in\mathbb{N}_{0}$ on a domain $D$ in one or two variables (depending on the Hausdorff dimension of $D$). Finally, we denote by $\mathbb{H}_{p}(D)$ the space of harmonic polynomials of degree $p\in\mathbb{N}_{0}$ on $D\subseteq\mathbb{R}^{2}$.

2 The model problem and the functional setting

Throughout the paper, we will employ the standard notation for Lebesgue and Sobolev spaces on a domain $D$, see [1]. In particular, we denote by $L^{2}(D)$ the Lebesgue space of square integrable functions and by $H^{s}(D)$, $s\in\mathbb{R}_{+}$, the Sobolev space $W^{2,s}(D)$. We set $\|\cdot\|_{0,D}$ the standard Lebesgue norm and $\|\cdot\|_{s,D}$ and $|\cdot|_{s,D}$ the Sobolev norms and seminorms, respectively.

We will use the following notation for partial derivatives:

[TABLE]

We will also write

[TABLE]

Moreover, we will employ the Sobolev weighted spaces and countably normed spaces defined e.g. in [5]. For the sake of completeness, we recall their definition. Given $\Omega\subseteq\mathbb{R}^{2}$ a bounded and simply connected polygonal domain, let $\mathcal{V}_{\Omega}$ be the number of vertices of the closure of $\Omega$ and let $\{\textbf{A}_{i}\}_{i=1}^{\mathcal{V}_{\Omega}}$ be the set of such vertices. We introduce the weight function

[TABLE]

where $|\cdot|$ denotes the Euclidean norm and $\boldsymbol{\beta}\in[0,1)^{\mathcal{V}_{\Omega}}$ is the weight vector. We will write $\Phi_{n}$, $n\in\mathbb{N}_{0}$, meaning that we will consider a weight vector $\boldsymbol{\beta}$ with constant entries $\boldsymbol{\beta}_{i}=n$, $\forall i=1,\dots,\mathcal{V}_{\Omega}$. Furthermore, we will denote the particular weight function $\Phi_{1}$ by $\Phi$.

The weighted Sobolev space $H_{\boldsymbol{\beta}}^{m,\ell}(\Omega)$, $\boldsymbol{\beta}\in[0,1)^{\mathcal{V}_{\Omega}}$, $m,\ell\in\mathbb{N}_{0}$, $m\geq\ell$, are defined as the completion of $\mathcal{C}^{\infty}(\overline{\Omega})$ with respect to the norms

[TABLE]

With an abuse of notation, the sum between the vector $\boldsymbol{\beta}$ and the scalar $k-\ell$ is meant to be

[TABLE]

Given $\ell\in\mathbb{N}_{0}$ and $\boldsymbol{\beta}\in[0,1)^{\mathcal{V}_{\Omega}}$, we define the countably normed spaces (or Babuška spaces) as

[TABLE]

where $c_{u}$ and $d_{u}$ are two constants greater than or equal to 1, depending only on the function $u$.

We define $\mathcal{B}_{\boldsymbol{\beta}}^{\frac{3}{2}}(\partial\Omega)$ and $\mathcal{\cacca}_{\boldsymbol{\beta}}^{\frac{3}{2}}(\partial K)$ as the set of the traces of functions belonging to $\mathcal{B}_{\boldsymbol{\beta}}^{2}(\Omega)$ and $\mathcal{\cacca}_{\boldsymbol{\beta}}^{2}(\Omega)$, respectively.

From (2) and (5), given $u\in\mathcal{\cacca}^{2}_{\boldsymbol{\beta}}(\Omega)$, for every $\boldsymbol{\alpha}\in\mathbb{N}_{0}^{2}$, $|\boldsymbol{\alpha}|=k\geq 1$, $k\in\mathbb{N}$, it holds that

[TABLE]

since $\beta-1\in[-1,0)$.

As a consequence, any function in $\mathcal{\cacca}^{2}_{\boldsymbol{\beta}}(\Omega)$ admits an analytic continuation on

[TABLE]

In order to see this, it suffices to show that the Taylor series

[TABLE]

converges uniformly in $\mathcal{N}(u)$.

In particular, we prove that it converges uniformly in the ball $B\left(\textbf{x}_{0},c\frac{\Phi(\textbf{x}_{0})}{d_{u}}\right)$ for all $c\in\left(0,\frac{1}{2}\right)$, where $\textbf{x}_{0}\in\overline{\Omega}$, $\textbf{x}_{0}\neq\textbf{A}_{i},$ $i=1,\dots,\mathcal{V}_{\Omega}$. In other words, we have to prove

[TABLE]

where $\overline{c}$ is a positive constant depending only on function $u$.

Using (6) and the fact that x belongs to $B\left(\textbf{x}_{0},c\frac{\Phi(\textbf{x}_{0})}{d_{u}}\right)$, we obtain

[TABLE]

since we are assuming that $c\in\left(0,\frac{1}{2}\right)$.

In this paper, we concentrate on the model problem given by the Laplace equation in $\Omega$ endowed with nonhomogeneous Dirichlet boundary conditions. Given $g:\partial\Omega\to\mathbb{R}$, find $u:\Omega\to\mathbb{R}$ satisfying

[TABLE]

The weak formulation reads:

[TABLE]

where

[TABLE]

It is well known that problem (11) is well-posed.

Assuming that the Dirichlet datum $g\in\mathcal{B}^{\frac{3}{2}}_{\boldsymbol{\beta}}(\partial\Omega)$, then the solution to problem (11) belongs to $\mathcal{B}^{2}_{\boldsymbol{\beta}}(\Omega)$ and $\mathcal{\cacca}^{2}_{\boldsymbol{\beta}}(\Omega)$ defined in (5), see [5, Theorem 2.2] and [33, Theorem 4.44]. Owing to (6) and the subsequent argument, $u$ is analytic on $\mathcal{N}(u)$ defined in (7).

Before concluding this section, we make the following simplifying assumption:

[TABLE]

As a consequence of (12), the solution $u$ of (11) is assumed to be analytic far from the singularity at $\mathbf{0}$. The general case of multiple corner singularities can be treated analogously. The main result of the paper Theorem 4.6, namely the exponential convergence of the energy error in terms of the number of degrees of accuracy, remains valid also if $u$ is singular at all the other vertices.

Henceforth, we also assume that the Babuška and weighted Sobolev spaces introduced above are defined taking into account in their definition only the singularity at $\mathbf{0}$. More precisely, we define such spaces by modifying the weight function in (3) into $\Phi_{\beta}(\textbf{x})=\min(1,|\textbf{x}|)^{\beta}$, for some $\beta\in[0,1)$, that is, the weight function associated only with the vertex $\mathbf{0}$.

3 Harmonic virtual element method with nonuniform degrees of accuracy

In this section, we introduce a method for the approximation of problem (11) employing polygonal meshes. This method takes the name of harmonic virtual element method (henceforth harmonic VEM) and is a modification of the standard virtual element method (henceforth VEM) tailored for the approximation of solutions to harmonic problem.

Let $\{\mathcal{T}_{n}\}$ be a sequence of polygonal decompositions of $\Omega$. Let $\mathcal{V}_{n}$ ($\mathcal{V}_{n}^{b}$) and $\mathcal{E}_{n}$ ($\mathcal{E}_{n}^{b}$) be the set of (boundary) vertices and edges of decomposition $\mathcal{T}_{n}$, respectively. We assume that $\mathcal{T}_{n}$ is a conforming decomposition for all $n\in\mathbb{N}_{0}$, that is to say that each boundary edge is an edge of only one element of $\mathcal{T}_{n}$, whereas each internal edge is an edge of precisely two elements of $\mathcal{T}_{n}$. Given $K\in\mathcal{T}_{n}$, we denote by $\mathcal{V}^{K}$ and $\mathcal{E}^{K}$ the set of vertices and edges of the polygon $K$. Moreover, we set $h_{K}:=\text{diam}(K)$ the diameter of polygon $K$, for all $K\in\mathcal{T}_{n}$, and $h_{s}:=|s|$ the length of edge $s$, for all $s\in\mathcal{E}^{K}$. Note that hanging nodes, i.e. multiple edges on a straight line, are allowed.

We require the following two assumptions on the polygonal decomposition $\mathcal{T}_{n}$.

• (D1)

Every $K\in\mathcal{T}_{n}$ is star-shaped with respect to a ball of radius greater than or equal to $\rho_{0}h_{K}$, $\rho_{0}$ being a universal positive constant belonging to $(0,\frac{1}{2})$. Since there are many possible balls satisfying the star-shapedeness condition we fix for each $K\in\mathcal{T}_{n}$ one ball $B=B(K)$. Furthermore, for all $K\in\mathcal{T}_{n}$ abutting $\mathbf{0}$, the subtriangulation $\widetilde{\mathcal{T}}=\widetilde{\mathcal{T}}(K)$ obtained by joining the vertices of $K$ to $\mathbf{0}$ is made of triangles that are star-shaped with respect to a ball of radius greater than or equal to $\rho_{0}h_{T}$, $h_{T}$ being diam($T$) for all $T\in\widetilde{\mathcal{T}}$. For all $T\in\widetilde{\mathcal{T}}(K)$, it holds $h_{K}\approx h_{T}$.

• (D2)

For all edges $s\in\mathcal{E}^{K}$, $K\in\mathcal{T}_{n}$, it holds $h_{s}\geq\rho_{0}h_{K}$, $\rho_{0}$ being the same constant in the assumption (D1). Besides, the number of edges in K, is uniformly bounded independently on the geometry of the domain.

We define the local harmonic virtual element spaces. Given $p\in\mathbb{N}$ and given the following space, defined on the boundary of a polygon $K\in\mathcal{T}_{n}$ as

[TABLE]

we set

[TABLE]

The functions in $V^{\Delta}(K)$ are then the solutions to local Laplace problems with piecewise polynomial Dirichlet data; therefore, they are not known explicitly in closed form.

Let us consider the following set of linear functionals on $V^{\Delta}(K)$. Given $v\in V^{\Delta}(K)$:

• •

the values of $v$ at the vertices of $K$;

• •

the values of $v$ at the $p-1$ internal Gauß-Lobatto nodes of $s$, for all $s$ edges of $K$.

This is a set of degrees of freedom, since (i) the dimension of $V^{\Delta}(K)$ is equal to the number of functionals defined above and (ii) such functionals are uninsolvent, owing to the fact that weak harmonic functions that vanish on $\partial K$, vanish also in the interior of $K$. Thus, the dimension of space $V^{\Delta}(K)$ is finite and is equal to $\sum_{s\in\mathcal{E}^{K}}p=p\cdot\#(\text{edges of }K)$.

We note that the definition of the edge degrees of freedom as the values of Gauß-Lobatto nodes is not the only possible; for instance, modal degrees of freedom of integrated Legendre polynomials is suitable as well.

By $\text{dof}_{i}$ we denote the $i$-th degree of freedom of $V^{\Delta}(K)$, whereas by $\{\varphi_{i}^{K}\}_{i=1}^{\dim(V^{\Delta}(K))}$ we denote the canonical basis of $V^{\Delta}(K)$, i.e. the set of basis functions in $V^{\Delta}(K)$ given by

[TABLE]

where $\delta_{i,j}$ is the Kronecker delta. We define the global harmonic virtual element space

[TABLE]

its subspace having vanishing boundary trace

[TABLE]

and its affine subspace containing interpolated essential boundary conditions

[TABLE]

Here $g_{GL}^{s}$ is the Gauß-Lobatto interpolant of degree $p$ of $g$ on the edge $s$ and where we recall $\mathcal{E}_{n}^{b}$ is the set of boundary edges of $\mathcal{T}_{n}$. We remark that $g_{GL}^{s}$ is well defined, since $g\in\mathcal{B}_{\beta}^{\frac{3}{2}}(\Omega)$, which implies $g\in\mathcal{C}^{0}(\overline{\Omega})$, see [33, Proposition 4.3].

The global degrees of freedom in the spaces (16), (17) and (18) are obtained by a standard continuous matching between the degrees of freedom of local spaces and, in the latter case, by imposing proper polynomial Dirichlet boundary conditions.

The space $V_{n,0}$ (17) and the affine space $V_{n,g}$ (18) consist then of piecewise harmonic functions on each element, piecewise continuous polynomials on the skeleton and piecewise Gauß-Lobatto interpolant of the Dirichlet datum $g$ on the boundary. The name component “virtual” emphasizes that such functions are not known explicitly at the interior of each $K\in\mathcal{T}_{n}$, since they are weak solutions to local Laplace problems with polynomial Dirichlet boundary conditions. On the other hand, the name component “harmonic” emphasizes that functions in $V_{n}$ are piecewise harmonic.

We point out that the choice of Gauß-Lobatto interpolation of the Dirichlet datum (18) will play a role in the $p$ approximation estimates. However, other choices in order to have a proper $p$ approximation of the boundary datum could be performed; for instance, one could use integrated Legendre polynomials interpolation of the Dirichlet datum as well. We stick here to the choice of Gauß-Lobatto interpolation for the sake of clarity.

Having defined the approximation spaces, we introduce the harmonic VEM associated with (11):

[TABLE]

where $a_{n}(\cdot,\cdot)$ is an approximate symmetric bilinear form defined on the unrestricted space $V_{n}\times V_{n}$, see (16). We require that the bilinear form $a_{n}(\cdot,\cdot)$ is explicitly computable by means of the degrees of freedom of the space and it must mimic the properties of its continuous counterpart $a(\cdot,\cdot)$; in particular, appropriate continuity and coercivity properties on $a_{n}$ are required. We argue and derive a suitable representation of $a_{n}(\cdot,\cdot)$ step-by-step.

First of all, we recall the representation

[TABLE]

Thus it is natural to seek for $a_{n}(\cdot,\cdot)$ as a sum of its local contributions

[TABLE]

Here, the $a_{n}^{K}(\cdot,\cdot)$ are local discrete bilinear forms defined on $V^{\Delta}(K)\times V^{\Delta}(K)$.

Next, we impose the validity of the two following assumptions on $a_{n}^{K}(u_{n},v_{n})$:

• (A1)

local harmonic polynomial consistency: for all $K\in\mathcal{T}_{n}$, it must hold

[TABLE]

where we recall that $\mathbb{H}_{p}(K)$ is the space of harmonic polynomials of degree $p$ over $K$;

• (A2)

local stability: for all $K\in\mathcal{T}_{n}$, it must hold

[TABLE]

where $0<\alpha_{*}(p)\leq\alpha^{*}(p)<+\infty$ are two constants which may depend on the local space $V^{\Delta}(K)$. In particular, $\alpha_{*}$ and $\alpha^{*}$ must be independent of $h_{K}$.

The assumption (A2) is required to guarantee that the discrete bilinear form scales like its continuous counterpart. In particular, it implies the coercivity and the continuity of the discrete bilinear form $a_{n}$. This, along with Lax Milgram lemma, implies the well-posedness of the problem (19).

On the other hand, the assumption (A1) implies that the problem (19) passes the patch test, meaning that, if the solution to the continuous problem (11) is a piecewise discontinuous harmonic polynomial, then the method described in (19) returns exactly, up to machine precision, the exact solution. For this reason, $p$ can be regarded as the degree of accuracy of the method.

We now investigate the behaviour of the error in the energy norm. The following variation of the quasioptimality result for the discrete solution is an adaptation of [9, Lemma 1]. We define

[TABLE]

where $\alpha_{*}(p)$ and $\alpha^{*}(p)$ are introduced in (22), and the $H^{1}$-broken Sobolev seminorm associated with the polygonal decomposition $\mathcal{T}_{n}$

[TABLE]

Lemma 3.1**.**

We assume that the assumptions (A1) and (A2) are satisfied. Let $u$ and $u_{n}$ be the solutions to problems (11) and (19), respectively. Then, the following holds true:

[TABLE]

where $\alpha(p)$ and $V_{n,g}$ are defined in (23) and (18), respectively, and where $S^{p,\Delta}(\Omega,\mathcal{T}_{n})$ is the space of (globally discontinuous) piecewise harmonic polynomials of degree $p$ on each $K\in\mathcal{T}_{n}$.

Proof.

A triangle inequality yields

[TABLE]

Owing to the assumptions (A1) and (A2), and the problems (11) and (19), one gets

[TABLE]

The claim follows plugging (27) in (26) and from simple algebra. ∎

Lemma 3.1 states that the energy error arising from the method can be bounded by a sum of local contributions of best local error terms with respect to the space of harmonic polynomials and to the space of functions in the harmonic virtual element space (14). We note that such best errors are weighted by the factor $\alpha(p)$ defined in (23).

We exhibit now an explicit choice for $a_{n}^{K}(\cdot,\cdot)$. To this end, we need to define a local energy projection from the local harmonic virtual element space $V^{\Delta}(K)$ defined (14) into $\mathbb{H}_{p}(K)$, which we recall is the space of harmonic polynomials of degree $p$ over $K$. We then introduce the projector $\Pi^{\nabla,K}_{p}$ defined as

[TABLE]

The second equation in (28) only fixes constants and can be substituted by other computable choices, see [7, 2]. Henceforth, when no confusion occurs, we will write $\Pi^{\nabla}$ in lieu of $\Pi^{\nabla,K}_{p}$.

We note that the projector $\Pi^{\nabla}$ can be computed by means of the dofs of space $V^{\Delta}(K)$. In fact, it suffices to apply an integration by parts to get

[TABLE]

where $\mathbf{n}$ denotes the exterior normal versor on the boundary of $K$, $\partial_{\mathbf{n}}q$ denotes the associated normal derivative and where we used that $q$ is harmonic, i.e. $\Delta q=0$. In order to conclude, it suffices to note that both $v$ and $\partial_{\mathbf{n}}q$ are explicity known on $\partial K$.

Let now $S^{K}:\ker(\Pi^{\nabla})\times\ker(\Pi^{\nabla})\rightarrow\mathbb{R}$ be any computable bilinear form satisfying the following stability assumption:

[TABLE]

where $0<c_{*}(p)\leq c^{*}(p)<+\infty$ are two constants which may depend on the local space $\ker(\Pi^{\nabla})$. An explicit selection for $S^{K}$ and a derivation of explicit bounds on $c_{*}(p)$ and $c^{*}(p)$ in terms of $p$ and $h_{K}$ are the topic of Section 3.1.

At this point, we are ready to define the local discrete bilinear form. We set

[TABLE]

We observe that the local stability property (29) implies the validity of the assumptions (A1) and (A2). In particular, the assumption (A2) holds with

[TABLE]

In Sections 3.1 and 3.2, we investigate the behaviour of $\alpha(p)$ in terms of $p$ for particular choices of the stabilization $S^{K}$ satisfying (29).

Remark 1*.*

So far, we have assumed that the Laplace problem (10) is endowed with Dirichlet boundary conditions. In the case of the Laplace problem $\Delta u=0$ with mixed boundary conditions

[TABLE]

over two parts of the boundary $\partial\Omega=\overline{\Gamma_{1}}\cap\overline{\Gamma_{2}}$ having nonzero measure, the right-hand side of the weak formulation (19) is augmented by the term $(g_{2},v)_{0,\Gamma_{2}}$.

3.1 A stabilization with the $L^{2}$-norm on the skeleton

In this section we introduce a computable local stabilizing bilinear form $S^{K}$ satisfying (29) and obtain explicit bounds in terms of the local degree of accuracy $p$ for the corresponding stabilization constants $c_{*}(p)$ and $c^{*}(p)$. Our first candidate is

[TABLE]

Since functions in $V^{\Delta}(K)$, defined in (14), are piecewise polynomials on the boundary of the element, then it is clear that the local stabilization introduced in (32) is explicitly computable.

For computational purposes, we substitute the edge integrals on the right-hand side of (32) with Gauß-Lobatto quadratures. This new choice is spectrally equivalent to the one in (32). Indeed, recalling [12, (2.14)] and setting $\widehat{I}=[-1,1]$, $\{\widehat{\eta}_{j}^{p}\}_{i=0}^{p}$ and $\{\widehat{\xi}_{j}^{p}\}_{i=0}^{p}$ the Gauß-Lobatto weights and nodes on $\widehat{I}$, then there exists a positive universal constant $c$ such that

[TABLE]

A scaling argument in addition to the assumption (D2) guarantees that the terms of the sum on the right-hand side of (32) can be replaced with Gauß-Lobatto quadrature formulas. This last choice is, from the computational point of view, more convenient than (32), since it results in diagonal matrix blocks. Thus, we emphasize our choice of $S^{K}$ by writing explicitly its definition. To each $s\in\mathcal{E}^{K}$ we associate the set of Gauß-Lobatto weights and nodes $\{\eta_{j}^{p,s}\}_{j=0}^{p}$ and $\{\xi_{j}^{p,s}\}_{j=0}^{p}$, respectively. The local stabilizing bilinear form associated with method (19) reads

[TABLE]

Next, we discuss the issue of showing explicit stability bounds (29) in terms of the local degree of accuracy.

Let us denote by

[TABLE]

the mean value of $v$ over $K\in\mathcal{T}_{n}$. Then the Poincaré inequality, see e.g. [13], implies

[TABLE]

Moreover, when $v\in\ker(\Pi^{\nabla})$, the following improved estimate is valid.

Lemma 3.2**.**

Let $K\in\mathcal{T}_{n}$ and let $\Pi^{\nabla}$ be defined in (28). For any $v\in\ker(\Pi^{\nabla})$, the following holds true:

[TABLE]

where $\lambda_{K}$ and $\omega_{K}$ denote the smallest exterior and the largest interior angle of $K$, respectively.

Proof.

We prove the assertion only for $K$ convex, i.e. $0<\omega_{K}<\pi$, since the nonconvex case can be treated analogously. Moreover, we assume without loss of generality that $h_{K}=1$. The general form of the assertion (37) follows then by a scaling argument.

The proof is based on an Aubin-Nitsche-type argument. For a fixed $v\in\ker(\Pi^{\nabla})$, we consider an auxiliary problem of finding $\eta$ such that

[TABLE]

where we recall that $\overline{v}$ is defined in (35).

We observe that by construction the right-hand side in (38) has vanishing mean and thus by the Lax-Milgram lemma the solution $\eta\in H^{1}(K)$ is well defined. The additional regularity of $\eta$ depends on the size of interior angles of $K$. In particular, if $K$ is convex, there holds $\eta\in H^{2}(K)$. More precisely,

[TABLE]

see e.g. [33, Section 4.2].

In the following, we utilize the additive splitting $\eta=\eta_{1}+\eta_{2}$, where the summands satisfy

[TABLE]

Again, standard a priori regularity results entail for a convex $K$

[TABLE]

Therefore, a combination of (39) and (40) with a triangle inequality, yields

[TABLE]

Besides, given any $w\in H^{1}(K)$ which is also harmonic, one has

[TABLE]

Recalling that $v\in\ker(\Pi^{\nabla})$ and applying sequentially (38), an integration by parts, (42), orthogonality of $\Pi^{\nabla}$, the Cauchy-Schwarz inequality and [27, Theorem 2], we deduce

[TABLE]

where $\lambda_{K}$ denotes the smallest exterior angle of $K$.

Plugging (43) in (41), we get the assertion. ∎

Now, we are ready to prove stability estimates for the spectrally equivalent $L^{2}$-norm stabilizations introduced in (32) and (34).

Lemma 3.3**.**

The bilinear forms $S^{K}$ defined in (32) and (34) fulfill the two-sided estimate (29) with constants satisfying

[TABLE]

where $\lambda_{K}$ and $\omega_{K}$ denote the smallest exterior and the largest interior angles of $K$, respectively.

Proof.

We prove the assertion only for $K$ convex, i.e. $0<\omega_{K}<\pi$, since the nonconvex case can be treated analogously. Moreover, in view of (33), it suffices to consider the bilinear form $S^{K}$ from (32). We also assume $h_{K}$=1 since the assertion will follow by a scaling argument.

We start by proving the lower bound for $c_{*}(p)$. Given $v\in\ker(\Pi^{\nabla})$, we write

[TABLE]

where we used an integration by parts and the fact that $v$ is harmonic in $K$. We apply now a Neumann trace inequality [33, Theorem A33] with $\Delta v=0$ in $K$, in order to show that

[TABLE]

Plugging (46) in (45) and using the polynomial $hp$ inverse inequality on an interval [33, Theorem 3.91] and interpolation theory [35], we obtain

[TABLE]

which is the asserted bound on $c_{*}(p)$.

Next, we investigate the behaviour of $c^{*}(p)$. Given $v\in\ker(\Pi^{\nabla})$ and $\overline{v}$ defined as in (35), one has

[TABLE]

We observe that, by (28), $v$ has zero boundary mean and therefore, by the Cauchy-Schwarz inequality,

[TABLE]

Hence by (47), (48), the multiplicative trace inequality and (37), we deduce

[TABLE]

where $\lambda_{K}$ denotes the smallest exterior angle of $K$. ∎

Lemma 3.3 and (31) imply that $\alpha(p)$ introduced in (23) admits the upper bound

[TABLE]

where $\lambda_{K}$ and $\omega_{K}$ denote the smallest exterior and largest interior angles of $K$, for all $K\in\mathcal{T}_{n}$, respectively.

We emphasize that the corresponding stability constant obtained for the standard (i.e. nonharmonic) $hp$ virtual element method, see [9, Theorem 2], grows much faster in $p$ than $\alpha(p)$. More precisely, it was proven that

[TABLE]

where, for all $K\in\mathcal{T}_{n}$, $\omega_{K}$ denotes the largest interior angle of $K$.

We conclude this section by noting that the stabilization introduced in (32) is basically, up to a $p$ scaling, the weighted (with Gauß-Lobatto weights) boundary contribution of the standard VEM stabilization introduced in [6, 7].

3.2 An optimal stabilization with the $H^{1/2}$-norm on the skeleton

In view of Theorem 4.6, which guarantees exponential convergence of the method in terms of the number of degrees of freedom, the mild behaviour of the stability constants $c_{*}(p)$ and $c^{*}(p)$ described in Lemma 3.3 in terms of $p$ has no effect on the asymptotic convergence rate of the method this remains exponential.

However, it is worth mentioning that there exists an optimal stabilization bilinear form $S^{K}$ with uniformly bounded stability constants $c_{*}$ and $c^{*}$; such stabilization reads

[TABLE]

where $(\cdot,\cdot)_{\frac{1}{2},\partial K}$ in the inner product on the Hilbert space $H^{\frac{1}{2}}(\partial K)$.

Lemma 3.4**.**

Let $S^{K}$ be defined as in (51). Then, for all $v\in\ker(\Pi^{\nabla})$, $\Pi^{\nabla}$ being defined in (28), the following holds true:

[TABLE]

Proof.

The statement follows from the proof of Lemma 3.3 and a scaling argument. ∎

It can be expected that the evaluation of (51) is more involved than the evaluation of the other variants of stabilization presented in Section 3.1, namely those in (32) and (34). In the following, we briefly discuss evaluation of the local stabilization (51).

We firstly recall the definition of the Aronszajn-Slobodeckij $H^{\frac{1}{2}}$ inner product over $\partial K$

[TABLE]

where $N_{s}^{K}$ denotes the number of edges of $K$ and $\{s_{i}\}_{i=1}^{N_{s}^{K}}$ denotes its set of edges. We observe that, owing to the fact that the stabilization is defined on $\ker(\Pi^{\nabla})$, it is possible to drop in (52) the contribution of the $L^{2}$ inner product.

We discuss now the evaluation of the double integral $I_{ij}$ in (52). We distinguish three different variants of the mutual locations of two edges $s_{i}$ and $s_{j}$.

$s_{i}$ and $s_{j}$ are identical ($s_{i}\equiv s_{j}$). In this case, the integrand in (52) has a removable singularity and is, in fact, a polynomial of degree $2p-2$. Such an integral is computed exactly by means of a Gauß-Lobatto quadrature formula with $p+1$ points. 2. 2.

$s_{i}$ and $s_{j}$ are distant ($\overline{s}_{i}\cap\overline{s}_{j}=\emptyset$). In this case, the integrand in (52) is an analytic function and can be efficiently approximated e.g. by a Gauß-Lobatto quadrature rule, see e.g. [15, Theorem 5.4]. 3. 3.

$s_{i}$ and $s_{j}$ share a vertex $\vec{v}$ and make an interior angle $0<\varphi<2\pi$. Then, $s_{i}$ and $s_{j}$ admit local parametrizations

[TABLE]

for some $\vec{a}$ and $\vec{b}\in\mathbb{R}^{2}$. Since the functions $u,v\in V^{\Delta}(K)$ are polynomials of degree $p$ along $s_{i}$ and $s_{j}$ and are continuous in $\vec{v}$ there holds

[TABLE]

where $f,g,q$ and $r$ are polynomials of degree $p-1$ and one has, using a change of coordinate,

[TABLE]

The integrand $F(s,t)$ is not smooth in $(0,1)^{2}$ (its derivatives blow up near the origin) and is not even defined in the origin, but it becomes regular after a coordinate transformation [19]. Having split the integral over the square $(0,1)^{2}$ into a sum of integrals over the two triangles obtained by bisecting such square with the segment of endpoints $(0,0)$ and $(1,1)$, simple algebra yields

[TABLE]

after the transformation $s=tz$ in the inner integral. The integrand admits the representation

[TABLE]

which is a rational function with a uniformly positive denominator

[TABLE]

Hence, the integrand (56) is an analytic function and can be efficiently approximated by Gauß-Lobatto quadrature, see e.g. [12].

Remark 2*.*

In [11], in the context of the approximation of a 2D Poisson problem, the possibility of using a stabilization involving only the boundary degrees of freedom was proven. More precisely, a stabilization equal to the boundary $H^{1}$ norm was employed; such norm can be related to the one introduced in (51) via $hp$ polynomial inverse estimates in one dimension. However, the analysis of [11] is not proven for the $p$ version of the method and therefore it is not clear whether the boundary stabilization therein proposed can be employed also for the $p$ analysis.

4 Exponential convergence with geometric graded polygonal meshes

In this section, we prove that, employing geometric refined towards $\mathbf{0}$ meshes and choosing appropriately a distribution of local degrees of accuracy, lead to exponential convergence of the energy error in terms of the dimension of the space, that is, in terms of the number of degrees of freedom.

We split the analysis as follows. In Section 4.1, we introduce the concept of sequences of polygonal meshes that are geometrically graded towards $\mathbf{0}$ (we recall that we are assuming that $\mathbf{0}$ is the unique “singular vertex” of $\Omega$, see (12)). In Section 4.2, we discuss the approximation results by harmonic polynomials, whereas in Section 4.3 we discuss the approximation results by functions in the harmonic virtual element space. Finally, in Section 4.4, we prove, under a proper choice of the vector of the degrees of accuracy, exponential convergence of the energy error in terms of the number of the degrees of freedom.

4.1 Geometric meshes

We describe sequences of geometrically graded meshes that we will employ for proving Theorem 4.6. Let $\sigma\in(0,1)$ be a given parameter. The sequence $\{\mathcal{T}_{n}\}$ is such that $\mathcal{T}_{n}$ consists then of $n+1$ “layers” for every $n\in\mathbb{N}_{0}$, where the “layers” are defined as follows.

We set the [math]-th layer $L_{n,0}=L_{0}$ as the set of all polygons $K\in\mathcal{T}_{n}$ abutting $\mathbf{0}$, which we recall is the unique “singular corner” of $\Omega$ by the assumption (12). The other layers are defined by induction as

[TABLE]

Next, we describe a procedure for building geometric polygonal graded meshes. Let $\mathcal{T}_{0}=\{\Omega\}$. The decomposition $\mathcal{T}_{n+1}$ is obtained by refining decomposition $\mathcal{T}_{n}$ only at the elements in the finest layer $L_{0}$. In order to have a proper geometric graded sequence of nested meshes, we demand for the following assumption.

• (D3)

[TABLE]

A consequence of the assumption (D3) is that $h_{K}\approx\sigma^{n-j}$, $j$ being the layer to which $K$ belongs. This, in addition to (60) guarantees that the distance between $K\in L_{j}$, $j=1,\dots,n$ and $\mathbf{0}$ is proportional to $\sigma^{n-j}$. Moreover, following [21, (5.6)], it can be shown that the number of elements in each layer is uniformly bounded with respect to all the geometric parameters discussed so far.

The sequence of nested meshes that we build is then characterized by very small elements near the singularity, whereas the size of the elements increases proportionally with the distance between the elements theirselves and $\mathbf{0}$.

Example 4.1**.**

In Figure 1, we depict three polygonal meshes satisfying the assumption (D3). We observe that the mesh in Figure 1 (right) does not fulfill the star-shapedness assumption (D1). We depict with different colours polygons belonging to different layers.

4.2 Approximation by harmonic polynomials

Here, we discuss approximation estimates by means of harmonic polynomials. Such results will be used for the approximation of the first term in (25), that is the best approximation in the $H^{1}$ seminorm of the solution to (11) by harmonic polynomials.

We will firstly deal with the approximation by harmonic polynomials on the polygons that are far from the singularity, see Lemma 4.2. Secondly, we will discuss the approximation estimates by harmonic polynomials on the polygons abutting the singularity, see Lemma 4.3.

Before that, we recall a (technical) auxiliary result, involving approximation on a polygon $K$ with $h_{K}=1$ by means of harmonic polynomials. The proof of this theorem can be found in [21, Theorem 4.10] and relies on the results in the pioneering works of [26, 28].

Theorem 4.1**.**

Let $\widehat{K}$ be a polygon with $h_{\widehat{K}}=1$. In particular, meas(\widehat{K})$$<1. We assume that the following parameters are given:

[TABLE]

where we recall that $\rho_{0}$ is the radius of the ball with respect to which $\widehat{K}$ is star shaped, see the assumption (D1). Let also:

[TABLE]

Then, there exists a sequence $\{\widehat{q}_{p}\}_{p=1}^{\infty}$, $\widehat{q}_{p}\in\mathbb{H}_{p}(\widehat{K})$ for all $p\in\mathbb{N}$, of harmonic polynomials such that, for any $\widehat{u}\in W^{1,\infty}(\widehat{K}_{\delta})$,

[TABLE]

We do not discuss the proof of Thorem 4.1, but we point out that in order to have this result we are using the fact that $\rho_{0}$ introduced in the assumption (D1) is such that $\rho_{0}\in(0,\frac{1}{2})$, since [21, Theorem 4.10] holds true under this hypothesis.

As a consequence of Theorem 4.1, for all the regular (in the sense of the assumptions (D1) and (D2)) polygons $\widehat{K}$ with diameter $1$ it holds that there exists an harmonic polynomial $q_{p}$ of degree $p$ such that

[TABLE]

where $c$ and $b$ are two positive constants depending uniquely on $\rho_{0}$ introduced in the assumption (D1) and the “enlargement factor” $\delta$ introduced in (61). Since both $\rho_{0}$ and $\delta$ are for the time being fixed, then $c$ and $b$ are two positive universal constants.

We assume now that the polygon $K$ belongs to $L_{j}$, $j=1,\dots,n$ and consequently has the diameter unequal to $1$. Then, a scaling argument immediately yields

[TABLE]

where $\widehat{K}$, the polygon obtained by scaling $K$, is such that $h_{\widehat{K}}=1$, where $\{\widehat{q}_{p_{\widehat{K}}}\}_{p_{\widehat{K}}=1}^{\infty}$ is the sequence validating (64), where $K_{\varepsilon}$ is defined as in (62) and where the “enlargement” factor $\varepsilon$ must be chosen in such a way that when we scale $K$ to $\widehat{K}$, then $K_{\varepsilon}$ is mapped in $\widehat{K}_{\delta}$, $\delta$ being exactly the parameter fixed in (61).

We note that sequence $\{q_{p}\}_{p=1}^{\infty}$, which is the pull-back of $\{\widehat{q}_{p}\}_{p=1}^{\infty}$, consists of harmonic polynomials since it is the composition of a sequence of harmonic polynomials with a dilatation.

What we have to check is that the size of $K_{\varepsilon}$ is not too large. In particular, we want that $K_{\varepsilon}$ is kept separated from the singularity at $\mathbf{0}$, for all $L_{j}$, $j=1,\dots,n$.

Let $u$ be the solution to problem (11). Henceforth, we assume that $\text{dist}(K,\mathbf{0})<1$ (which is always valid if one takes $\Omega$, the domain of problem (10), small enough). From Section 2, we know that $u$, the solution to problem (10), is analytic on the set $\mathcal{N}(u)$ defined in (7). In particular, $u$ is analytic on the following domain depending on $K$:

[TABLE]

since $\mathcal{N}_{K}(u)\subset\mathcal{N}(u)$. This fact has an extreme relevance in the proof of forthcoming Lemma 4.2. The important issue is that more the polygon is near the singularity, the smaller is the extended domain $\mathcal{N}_{K}(u)$, see Figure 2.

In any case, $\mathcal{N}_{K}(u)$ remains contained in the global analyticity domain $\mathcal{N}(u)$, which is fixed once and for all.

We choose $c=\frac{1}{4}$ in (66). Owing to (60) and recalling that $K\notin L_{0}$, there exist two constants $0<\alpha_{1}\leq 1\leq\alpha_{2}$ independent of $K$ such that $\alpha_{1}h_{K}\leq\text{dist}(K,\mathbf{0})\leq\alpha_{2}h_{K}$. Thus,

[TABLE]

This implies that $u$ is analytic on the following domain too:

[TABLE]

Therefore, we fix for instance $\varepsilon=\frac{1}{8}\frac{\alpha_{1}}{d_{u}}h_{K}$. In this way, we have built $K_{\varepsilon}=\widetilde{\mathcal{N}}_{K}(u)$ neighbourhood of $K$ not covering $\mathbf{0}$.

It is straightforward to note that scaling $K$ to $\widehat{K}$ with $h_{\widehat{K}}=1$, we also scale $K_{\varepsilon}$ to $\widehat{K}_{\delta}$ (see (62) for the definition of $\widehat{K}_{\delta}$), where $\delta=\frac{1}{8}\frac{\alpha_{1}}{d_{u}}$ is now independent of $K$ and only depends on $u$. Fixing such a $\delta$ in Theorem 4.1, we have that (65) holds with $\mathbf{0}\notin\overline{K_{\varepsilon}}$; in particular, the norm appearing on the right-hand side of (65) is finite for all $K\in L_{j}$, $j=1,\dots,n$.

We are now ready to state the bound on the best error with respect to harmonic polynomials on the polygons not abutting the singularity.

Lemma 4.2**.**

Let the assumptions (D1)-(D3) hold true, let $K\in L_{j}$, $j=1,\dots,n$ and let $u\in W^{1,\infty}(\widetilde{\mathcal{N}}_{K}(u))$, where $\widetilde{\mathcal{N}}_{K}(u)$ is defined in (67). Then, there exists a sequence $\{q_{p}\}_{p=1}^{\infty}\subseteq\{\mathbb{H}_{p}(K)\}_{p=1}^{\infty}$ of harmonic polynomials such that

[TABLE]

where $b$ is a constant independent of $K$.

Proof.

The proof follows from Theorem 4.1 and the subsequent discussion. In particular, the first inequality in (68) follows from a scaling argument, whereas, the second one is a consequence of computations analogous to those in (9) and the definition of $\widetilde{\mathcal{N}}_{K}(u)$ in (67). ∎

It is clear from the above discussion that we must follow a different strategy for the elements in the first layer; in fact, here, the $W^{1,\infty}$ norm of $u$ is not finite in principle.

It holds in particular the following result.

Lemma 4.3**.**

Let the assumptions (D1)-(D3) hold true. Let $K\in L_{0}$. Let $u\in H^{2,2}_{\beta}(\Omega)$. Then, there exists $q_{1}\in\mathbb{P}_{1}(K)$ such that

[TABLE]

In particular, $q_{1}$ is a harmonic polynomial.

Proof.

The polynomial $q_{1}$ is given by the linear interpolant of $u$ at, for instance, three nonaligned vertices of K. The proof follows the lines of [9, Lemma 3]. ∎

Remark 3*.*

Lemma 4.3 suggests that one could also consider harmonic virtual element spaces with nonuniform degrees of accuracy, still guaranteeing exponential convergence for the $hp$ version of the method. In particular, one could consider a distribution of degrees of accuracy which grows linearly as the distance from the singularity increases, as depicted in Figure 3.

At the interface $s$ of two nondisjoint elements $K_{0}$ and $K_{1}$ in layers $L_{0}$ and $L_{1}$ one associates $p_{s}=\max(1,p)=p$ (maximum rule) in order to define nonuniform boundary spaces $\mathbb{B}(\partial K)$ similarly to (13), as depicted in Figure 4.

In this section, we have thus built a piecewise discontinuous harmonic polynomial with certain approximation properties described in Lemmata 4.2 and 4.3. Such a discontinuous function will be used in the proof of Theorem 4.6 in the approximation of the first term on the right-hand side of (25).

4.3 Approximation by functions in the harmonic virtual element space

Here, we discuss about approximation estimates by functions in the harmonic virtual element space which will be used for the approximation of the second term in (25). As in Section 4.2, we firstly investigate approximation estimates on polygons not abutting the singularity, see Lemma 4.4; secondly, we discuss approximation estimates of polygons in the finest layer $L_{0}$, see Lemma 4.5.

Lemma 4.4**.**

Let the assumptions (D1)-(D3) hold true. Let $K\in L_{j}$, $j=1,\dots,n$ and let $\beta\in[0,1)$. Let $g$, the Dirichlet datum of problem (11), belong to space $\mathcal{B}^{\frac{3}{2}}_{\beta}(\partial\Omega)$ and let $u$, the solution to problem (11), belong to space $\mathcal{B}^{2}_{\beta}(\Omega)$, see (5). Then, there exists $u_{I}\in V^{\Delta}(K)$ such that

[TABLE]

where we recall that $\sigma$ is the geometric grading parameter of the assumption (D3).

Proof.

Before proving the result, we observe that the $H^{m+1}(s)$ seminorm exists for all edges $s$ of $K$, since $u\in\mathcal{B}_{\beta}^{2}(\Omega)$ implies that $u$ is analytic far from the singularity.

Let us consider $u_{I}\in V^{\Delta}(K)$ defined as the weak solution to the following local Laplace problem:

[TABLE]

where $u_{\text{GL}}$ is the Gauß-Lobatto interpolant of degree $p$ of $u$ on each edge $s$. Then, using the fact that $u-u_{I}$ is harmonic and using a Neumann trace inequality [33, Theorem A.33], one gets

[TABLE]

for every $c\in\mathbb{R}$.

We deduce that we must deal with the boundary error term only. We fix $c=0$ in (70) (the case $c\neq 0$ will become important in the following). Since $u$ is analytic far from the singularity, we inherit the two following results from [12, Theorems 4.2 and 4.5]:

[TABLE]

Using interpolation theory [35], recalling from the assumption (D2) that $h_{s}\approx h_{K}$ and that the number of edges of each $K\in\mathcal{T}_{n}$ is uniformly bounded, yield

[TABLE]

We apply a multiplicative trace inequality [13, Theorem 1.6.6], the fact that the maximum number of edges of $K$ is uniformly bounded, see the assumption (D2), and the trivial bound $|a||b|\leq a^{2}+b^{2}$, $a$, $b\in\mathbb{R}$, getting

[TABLE]

Recalling the definition of the weighted Sobolev seminorms (4), one obtains

[TABLE]

Combining (60), (72) and (73), we deduce

[TABLE]

Finally, recalling from the assumption (D3) that $h_{K}\approx\sigma^{n-j}$, we get the claim by inserting (74) in (71). ∎

Next, we turn our attention to the approximation in the polygons belonging to the first layer.

Lemma 4.5**.**

Let the assumptions (D1)-(D3) hold true. Let $K\in L_{0}$ and let $\beta\in[0,1)$. Let $g$, the Dirichlet datum of problem (11), belong to space $\mathcal{B}^{\frac{3}{2}}_{\beta}(\partial\Omega)$ and let $u$, the solution to problem (11), belong to space $\mathcal{B}^{2}_{\beta}(\Omega)$ (5). Then, there exists $u_{I}\in V^{\Delta}(K)$ such that

[TABLE]

where we recall that $\sigma$ is the geometric grading parameter of the assumption (D3).

Proof.

Let $u_{I}$ be defined as in (69), with $u_{\text{GL}}$ being now the linear interpolant of $u$ on each edge $s$ of $K$. Let $\widetilde{\mathcal{T}}(K)$ be the subtriangulation of $K$ obtained by joining $\mathbf{0}$ with the other vertices of $K$. Such a subtriangulation is regular, see the assumption (D1).

From (70), we have

[TABLE]

We denote by $\widetilde{u}_{\text{GL}}$ the linear interpolant of $u$ over every $T\in\widetilde{\mathcal{T}}(K)$ at the three vertices of $T$. One obviosuly has $\widetilde{u}_{\text{GL}}=u_{\text{GL}}$ on $\partial K$. Applying a trace inequality, we get

[TABLE]

By picking $c$ the average of $u-\widetilde{u}_{\text{GL}}$ over $K$, applying a Poincaré inequality and recalling that $\text{card}(\widetilde{\mathcal{T}})$ is uniformly bounded, we deduce

[TABLE]

In order to conclude, we apply [33, Lemma 4.16] and (60) obtaining

[TABLE]

which holds true owing to the fact that $u\in\mathcal{B}^{2}_{\beta}(\Omega)$. ∎

Again, for the proof of Lemma 4.5, one could have used nonuniform degrees of accuracy as discussed in Remark 3.

In order to conclude this section, we highlight that we built in Lemmata 4.4 and 4.5 a continuous approximant of $u$, which belongs to space $V_{n,g}$ (18).

The $h$ version of harmonic VEM for quasi-uniform meshes.

Although the goal of this paper is to study the $hp$ version of harmonic VEM, it is worthwhile to mention that the $h$ version of the method employing sequences of quasi-uniform meshes, see e.g. [6] for the definition of quasi-uniform meshes, easily follows by combining Lemma 3.1, [27, Theorem 2] and Lemma 4.4.

In particular, assuming that $u$, the solution to problem (11), belongs to $H^{s+1}(\Omega)$, $s\in\mathbb{R}_{+}$, and that we employ harmonic virtual element spaces with a uniform degree of accuracy $p$, one gets

[TABLE]

where the hidden constant depends on $s$, on the shape of the elements in the mesh and the choice of the stabilization, but is independent of the mesh size $h$.

4.4 Exponential convergence

Here, we discuss the main result of the work, namely the exponential convergence of the energy error in terms of the number of degrees of freedom. In order to achieve such a result, we fix as a degree of accuracy

[TABLE]

The main result of the paper follows.

Theorem 4.6**.**

Let $\{\mathcal{T}_{n}\}_{n\in\mathbb{N}_{0}}$ be a sequence of polygonal decomposition satisfying the assumptions (D1)-(D3). Let $u$ and $u_{n}$ be the solutions to problems (11) and (19), respectively. Let $g$, the Dirichlet datum introduced in (11), belong to $\mathcal{B}^{\frac{3}{2}}_{\beta}(\partial\Omega)$. Then, the following holds true:

[TABLE]

where $b$ is a constant independent of the discretization parameters and $N$ is the number of degrees of freedom of $V_{n}$ defined in (18).

Proof.

We only give the sketch of the proof. Applying Lemma 3.1, bound (50), Lemmata 4.4 to 4.5 to the first term on the right-hand side of (25) along with standard $hp$ approximation strategies [33] and Lemmata 4.2 and 4.3 to the second term of the right-hand side of (25) along with [21, Theorem 5.5], we have

[TABLE]

for some $\widetilde{b}$ independent of the discretization parameters, $n+1$ being the number of layers in $\mathcal{T}_{n}$.

In order to conclude, it suffices to find out the relation between $n$ and $N$, the number of degrees of freedom of space $V_{n}$. To this end, we recall from [21, (5.6)] that in each layer $L_{j}$ there exists a fixed maximum number of elements, see the assumption (D3). Moreover, thanks to the assumption (D2), there exists a fixed maximum number of edges per element.

If we set $N_{\text{edge}}$ the maximum number of edges per element and $N_{\text{element}}$ the maximum number of elements per layer, we conclude that

[TABLE]

In particular, $\sqrt{N}\lesssim n$. This, along with (78), gives the assertion. ∎

5 Numerical results

5.1 Numerical results: $h$ version

In this section, we present numerical results validating the algebraic rate of convergence of the $h$ version of the method stated in (75).

To this purpose, we consider the following test case. Let $\Omega$, the domain of problem (11), be the square domain

[TABLE]

and let $u$, the solution to the problem, be

[TABLE]

which is an analytic harmonic function over $\mathbb{R}^{2}$.

Moreover, we observe that since the functions in the harmonic virtual element space are known only via their degrees of freedom, we cannot explicitly compute the energy error. Therefore, we study the following normalized broken $H^{1}$ error between $u$ and the energy projection of $u_{n}$:

[TABLE]

where $\Pi^{\nabla,K}_{p}$ is defined in (28), for all $K\in\mathcal{T}_{n}$.

Importantly, the rate of convergence of the error in (79) is the same as the one of the exact $H^{1}$ error. In order to see this, we apply a triangle inequality and the stability of the $H^{1}$ projection, to get

[TABLE]

and after that we apply Lemma 3.1, [27, Theorem 2] and Lemma 4.4.

We test the method employing sequences made of three types mesh, see Figure 5, namely a squares, a Voronoi-Lloyd and an hexagonal mesh.

We also pick as uniform degrees of accuracy $p=1$, $2$, $3$ and $4$. The numerical results are collected in Figure 6 and are in accordance with the expected rate of convergence in (75).

5.2 Numerical results: $hp$ version

In this section, we present numerical experiments validating the exponential rate of convergence of the $hp$ version of the method stated in Theorem 4.6. To this end, we consider the following test case. Let $\Omega$, the domain of problem (11), be the L-shaped domain

[TABLE]

and let $u$, the solution to (11), be

[TABLE]

where $r$ and $\theta$ are the polar coordinates of the real plane. We observe that the such a function belongs to $H^{\frac{5}{3}-\varepsilon}(\Omega)$, for all $\varepsilon>0$, but not to $H^{\frac{5}{3}}(\Omega)$, and moreover that $u$ is harmonic.

5.2.1 Numerical tests on polygonal geometric graded meshes

We consider sequences of meshes as those depicted in Figure 1 and we consider two different distributions of local degrees of accuracy.

We firstly investigate in Figure 7 the performances of the harmonic VEM choosing a distribution of degrees of accuracy $p$ as in (76). Under this choice, we know that Theorem 4.6 holds true.

Secondly, we investigate in Figure 8 the performances of the harmonic VEM by taking a nonuniform distribution of degrees of accuracy. In particular, we consider the (graded) distribution given by

[TABLE]

At the interface of two polygons in different layers one associate a polynomial degree $p_{s}$ via the maximum rule as in Figure 4, thus modifying straightforwardly the definition of the space $\mathbb{B}(\partial K)$ defined in (13). It is worth to notice that under choice (83) the dimension of space $V_{n}$ is asymptotically $\frac{1}{2}n^{2}$, $n+1$ being the number of layers. Such a dimension is comparable with the one of space $V_{n}$ assuming (76), which is asymptotically $n^{2}$.

In both figures, we consider sequences of meshes with different geometric refinement parameters $\sigma$; we recall that the properties fulfilled by $\sigma$ are discussed in the assumption (D3). We fix in particular $\sigma=\frac{1}{2}$, $\sigma=\sqrt{2}-1$ and $\sigma=(\sqrt{2}-1)^{2}$.

As observed already in Section 5.1, we study the computable error in (79) in lieu of the exact one, whose convergence in terms of the number of degrees of freedom is the same as the one of the exact $H^{1}$ error.

On the $y$-axis we consider the logarithm of the error defined in (79), while in the $x$-axis we put the square root of the number of degrees of freedom.

As already stated in Example 4.1, the meshes as those in Figure 1 (right) do not satisy the assumption (D1) and then, in principle, Theorem 4.6 does not apply. The numerical experiments in Figure 7 and 8 reveal that in this case the convergence deteriorates after few $hp$ refinements, especially for small $\sigma$.

On the other hand, the other two sequences of meshes, namely those whose elements are depicted in Figure 1 (left) and (center), have the expected exponential decay.

Importantly, the exponential convergence is still observed also under choice (83) of the local degrees of accuracy. Our conjecture is that Theorem 4.6 holds under (83) as well. Nonetheless, we avoid to investigate this issue, on the one hand, in order to avoid additional technicalities, on the other, because the dimension of space $V_{n}$ under choices (76) and (83) behaves like $n^{2}$ and $\frac{1}{2}n^{2}$, respectively. This means that the exponential decay is still valid with the same exponential rate in both cases.

5.2.2 Numerical comparison between $hp$ harmonic VEM and $hp$ VEM

We also perform a numerical comparison between the performances of the harmonic VEM discussed so far and the standard $hp$ version of VEM for the approximation of corner singularities, see [9]. The main difference is that in VEM internal degrees of freedom for each element are employed in order to take care of the approximation of the right-hand side in Poisson problems. This obviously leads to a nontrivial growth of the dimension of the space of approximation and in particular it leads to a decay of the energy error of the following sort:

[TABLE]

where $b$ is a positive constant independent of the discretization parameters and $N$ is the dimension of the virtual space; see [9, Theorem 3].

In Figure 9, we compare error (79) for the two methods employing the meshes in Figure 1 (left) and in Figure 1 (center). The grading parameter is $\sigma=\frac{1}{2}$.

In both cases, we consider a distribution of local degrees of accuracy as in (76). We note that the stabilization of the VEM differs from the one introduced in (32) for the harmonic VEM. For more details concerning the construction of the $hp$ VEM we refer to [9].

From Figure 9, it is possible to observe the faster decay of error (79) when employing the $hp$ harmonic VEM, see (77), when compared to the same error employing the $hp$ VEM, see (84).

Aknowledgements

L. M. acknowledges the support of the Austrian Science Fund (FWF) project F 65. Both authors aknowledge that the major part of the research presented in this paper has been carried out at the Institute of Mathematics of the University of Oldenburg, Germany. Moreover, they wish to thank Prof. M. J. Melenk of the Technische Universität Wien for fruitful discussions on the topic.