Non-conforming Harmonic Virtual Element Method: $h$- and $p$-Versions

Abstract

We study the $h$- and $p$-versions of non-conforming harmonic virtual element methods (VEM) for the approximation of the Dirichlet–Laplace problem on a 2D polygonal domain, providing quasi-optimal error bounds. Harmonic VEM do not make use of internal degrees of freedom. This leads to a faster convergence, in terms of the number of degrees of freedom, as compared to standard VEM. Importantly, the technical tools used in our $p$-analysis can be employed as well in the analysis of more general non-conforming finite element methods and VEM. The theoretical results are validated in a series of numerical experiments. The hp-version of the method is numerically tested, demonstrating exponential convergence with rate given by the square root of the number of degrees of freedom.

Appendix A: Details on the Implementation

In this section, we discuss some practical aspects concerning the implementation of the non-conforming harmonic VEM in 2D. We employ henceforth the notation of [16]. It is worth to underline that we present herein only the case with uniform degree of accuracy; the implementation of the $hp$ version is dealt with similarly. As a first step, we begin by fixing the notation for the various bases instrumental for the construction of the method.

Basis of $\mathbb P_{p-1}(e)$ for a given $e\in \mathcal {E}^K$. Using the same notation as in (12), we denote the basis of $\mathbb P_{p-1}(e)$, $e\in \mathcal E^K$, by $\{ m_{r}^e\}_{r=0,\ldots ,p-1}$. The choice we make is

$$\begin{aligned} m_{r}^e(\mathbf x):=\mathbb L_{r}\left( \phi _e^{-1}(\mathbf x) \right) \quad \forall r=0,\dots ,p-1, \end{aligned}$$

(63)

where $\phi _e: [-1,1] \rightarrow e$ is the linear transformation mapping the interval $[-1,1]$ to the edge $e$, and $\mathbb L_{r}$ is the Legendre polynomial of degree r over $[-1,1]$. We recall, see e.g. [56], for future use the orthogonality property

$$\begin{aligned} \left( m_r^e,m_s^e\right) _{0,e} = \frac{h_e}{2} \int _{-1}^{1} \mathbb L_r(t) \mathbb L_s(t) \, \text {d}t = \frac{h_e}{2r+1} \delta _{rs} \quad \forall r,s=0,\dots ,p-1, \end{aligned}$$

(64)

where $\delta _{rs}$ is the Kronecker delta (1 if $r=s$, 0 otherwise).

Basis of $\mathbb H_p(K)$ for a given $K\in \mathcal T_n$. We denote the basis of the space of harmonic polynomials $\mathbb H_p(K)$ by $\{ q_\alpha ^\Delta \}_{\alpha =1,\ldots ,n_p^\Delta }$, where $n_p^\Delta :=\dim \mathbb H_p(K) = 2p+1$. The choice we make for this basis is

$$\begin{aligned} \begin{aligned} q_1^\Delta (\mathbf x)&= 1;\\ q_{2l}^\Delta (\mathbf x)&= \sum _{k=1, \, k \text { odd }}^{l} (-1)^{\frac{k-1}{2}} \left( {\begin{array}{c}l\\ k\end{array}}\right) \left( \frac{x-x_K}{h_K} \right) ^{l-k} \left( \frac{y-y_K}{h_K} \right) ^{k} \quad \forall l=1,\dots ,p; \\ q_{2l+1}^\Delta (\mathbf x)&= \sum _{k=0, \, k \text { even }}^{l} (-1)^{\frac{k}{2}} \left( {\begin{array}{c}l\\ k\end{array}}\right) \left( \frac{x-x_K}{h_K} \right) ^{l-k} \left( \frac{y-y_K}{h_K} \right) ^{k} \quad \forall l=1,\dots ,p.\\ \end{aligned} \end{aligned}$$

The fact that this is actually a basis for $\mathbb H_p(K)$ is proven, e.g., in [6, Theorem 5.24].

Basis for $V^\Delta (K)$ for a given $K\in \mathcal T_n$. For this local VE space introduced in (11), we employ the canonical basis $\left\{ \varphi _{j,r} \right\} _{\begin{array}{c} j=1,\ldots , N_K\\ r=0,\ldots ,p-1 \end{array}}$ defined though (13), where we also recall that $N_K$ denotes the number of edges of K.

In the following, we derive the matrix representation of the local discrete bilinear form introduced in (23). We begin with the computation of the matrix representation of the projector $\Pi ^{\nabla ,K}_p$ acting from $V(K)$ to $\mathbb H_p(K)$ and defined in (15). To this purpose, given any basis function $\varphi _{j,r} \in V^\Delta (K)$, $j=1,\dots ,N_K$, $r=0,\dots ,p-1$, we expand $\Pi ^{\nabla ,K}_p\varphi _{j,r}$ in terms of basis $\{ q_\alpha ^\Delta \}_{\alpha =1,\ldots ,n_p^\Delta }$ of $\mathbb H_p(K)$, i.e.,

$$\begin{aligned} \Pi ^{\nabla ,K}_p\varphi _{j,r} = \sum _{\alpha =1}^{n_p^\Delta } s_\alpha ^{(j,r)} q_\alpha ^\Delta . \end{aligned}$$

(65)

Using (15) and testing (65) with functions $q_\beta ^\Delta$, $\beta =1,\dots ,n_p^\Delta$, we get that the coefficients $s_{\alpha }^{(j,r)}$ can be computed by solving for $\varvec{s}^{(j,r)}:=[s_1^{(j,r)},\ldots , s_{n_p^\Delta }^{(j,r)}]^T$ the $n_p^\Delta \times n_p^\Delta$ algebraic linear system

$$\begin{aligned} \varvec{G}\varvec{s}^{(j,r)}=\varvec{b}^{(j,r)}, \end{aligned}$$

where

$$\begin{aligned}&\varvec{G}= \begin{bmatrix} \left( q_1^\Delta ,1\right) _{0,\partial K}&\quad \left( q_2^\Delta ,1\right) _{0,\partial K}&\quad \cdots&\quad \left( q_{n_p^\Delta }^\Delta ,1\right) _{0,\partial K} \\ 0&\quad \left( \nabla q_2^\Delta ,\nabla q_2^\Delta \right) _{0,K}&\quad \cdots&\quad \left( \nabla q_{n_p^\Delta }^\Delta ,\nabla q_2^\Delta \right) _{0,K} \\ \vdots&\quad \vdots&\quad \ddots&\quad \vdots \\ 0&\quad \left( \nabla q_{n_p^\Delta }^\Delta ,\nabla q_2^\Delta \right) _{0,K}&\quad \cdots&\quad \left( \nabla q_{n_p^\Delta }^\Delta ,\nabla q_{n_p^\Delta }^\Delta \right) _{0,K} \end{bmatrix}, \nonumber \\&\quad \varvec{b}^{(j,r)}= \begin{bmatrix} \left( \varphi _{j,r},1\right) _{0,\partial K} \\ \left( \nabla \varphi _{j,r},\nabla q_2^\Delta \right) _{0,K} \\ \vdots \\ \left( \nabla \varphi _{j,r},\nabla q_{n_p^\Delta }^\Delta \right) _{0,K} \end{bmatrix}. \end{aligned}$$

Collecting all the $N_K p$ (column) vectors $\varvec{b}^{(j,r)}$ in a matrix $\varvec{B}\in \mathbb R^{n_p^\Delta \times N_K p}$, namely, setting $\varvec{B}:=[\varvec{b}^{(1,1)},\dots ,\varvec{b}^{(N_K,p)}]$, the matrix representation $\varvec{\Pi }_*$ of the projector $\Pi ^{\nabla ,K}_p$ acting from $V^\Delta (K)$ to $\mathbb H_p(K)$ is given by

$$\begin{aligned} \varvec{\Pi }_*= \varvec{G}^{-1} \varvec{B} \in \mathbb R^{n_p^\Delta \times N_K p}. \end{aligned}$$

Subsequently, we define

$$\begin{aligned} \varvec{D}:=\begin{bmatrix} \text {dof}_{1,1}\left( q_1^\Delta \right)&\quad \cdots&\quad \text {dof}_{1,1}\left( q_{n_p^\Delta }^\Delta \right) \\ \vdots&\quad \ddots&\quad \vdots \\ \text {dof}_{N_K, p}\left( q_1^\Delta \right)&\quad \cdots&\quad \text {dof}_{N_K, p}\left( q_{n_p^\Delta }^\Delta \right) \end{bmatrix} \in \mathbb R^{N_Kp \times n_p^\Delta }. \end{aligned}$$

Let $\varvec{\Pi }$ be the matrix representation of the operator $\Pi ^{\nabla ,K}_p$ seen now as a map from $V^\Delta (K)$ into $V^\Delta (K)\supseteq \mathbb H_p(K)$. Then, following [16], it is possible to show that

$$\begin{aligned} \varvec{\Pi }= \varvec{D} \varvec{G}^{-1} \varvec{B} \in \mathbb R^{N_K p \times N_K p}. \end{aligned}$$

Next, denoting by $\widetilde{\varvec{G}} \in \mathbb R^{n_p^\Delta \times n_p^\Delta }$ the matrix coinciding with $\varvec{G}$ apart from the first row which is set to zero, the matrix representation of the bilinear form in (23) is

$$\begin{aligned} \left( \varvec{\Pi }_*\right) ^T \, \widetilde{\varvec{G}} \,\left( \varvec{\Pi }_*\right) + \left( \varvec{I} - \varvec{\Pi }\right) ^T \, \varvec{S} \, \left( \varvec{I} - \varvec{\Pi }\right) . \end{aligned}$$

Here, $\varvec{S}$ denotes the matrix representation of an explicit stabilization $S^K(\cdot ,\cdot )$. For the stabilization defined in (28), we have

$$\begin{aligned}&\varvec{S}((k-1)N_K+ r,(l-1)N_K+s) = \sum _{i=1}^{N_K} \frac{p}{h_{e_i}} \left( \Pi ^{0,e_i}_{p-1}\varphi _{l,s}, \Pi ^{0,e_i}_{p-1}\varphi _{k,r}\right) _{0,e_i}\\&\quad \forall k,l=1,\dots ,N_K, \forall r,s=0,\dots ,p-1. \end{aligned}$$

By expanding $\Pi ^{0,e_i}_{p-1}\varphi _{l,s}$ and $\Pi ^{0,e_i}_{p-1}\varphi _{k,r}$ in the basis $\left\{ m_\gamma ^{e_i} \right\} _{\gamma =0,\ldots ,p-1}$ of $\mathbb P_{p-1}(e_i)$, i.e.,

$$\begin{aligned}&\Pi ^{0,e_i}_{p-1}\varphi _{l,s} = \sum _{\gamma =0}^{p-1} t_\gamma ^{(l,s),e_i} m_\gamma ^{e_i}, \quad \Pi ^{0,e_i}_{p-1}\varphi _{k,r}= \sum _{\zeta =0}^{p-1}t_\zeta ^{(k,r),{e_i}} m_\zeta ^{e_i}, \nonumber \\&\quad \forall k,l=1,\dots ,N_K,\, \forall r,s=0,\dots ,p-1, \end{aligned}$$

(66)

we can write

$$\begin{aligned} \begin{aligned}&\varvec{S}((k-1)N_K+ r,(l-1)N_K+s) = \sum _{i=1}^{N_K} \sum _{\gamma =0}^{p-1} \sum _{\zeta =0}^{p-1} t_\gamma ^{(l,s),{e_i}} t_\zeta ^{(k,r),{e_i}} \frac{p}{h_{e_i}} \left( m_\gamma ^{e_i}, m_\zeta ^{e_i}\right) _{0,{e_i}}\nonumber \\&\quad \forall k,l=1,\dots ,N_K,\, \forall r,s=0,\dots ,p-1. \end{aligned} \end{aligned}$$

For the basis defined in (63), using the orthogonality of the Legendre polynomials (64), this expression can be simplified leading to a diagonal stability matrix $\mathbf S$:

$$\begin{aligned} \begin{aligned}&\varvec{S}((k-1)N_K+ r,(k-1)N_K+r)= \sum _{i=1}^{N_K} \sum _{\zeta =0}^{p-1} \frac{p}{2r+1} \left( t_\zeta ^{(k,r),{e_i}}\right) ^2 \nonumber \\&\quad \forall k=1,\dots ,N_K,\, \forall r=0,\dots ,p-1. \end{aligned} \end{aligned}$$

For fixed $i,k \in \{1,\dots ,N_K\}$ and $r \in \{0,\dots ,p-1\}$, the coefficients $t_\zeta ^{(k,r),e_i}$ are obtained by testing $\Pi ^{0,e_i}_{p-1}\varphi _{k,r}$, defined in (66), with $m_{\zeta }^{e_i}$, $\zeta =0,\dots ,p-1$, and by taking into account the definition of $\Pi ^{0,e_i}_{p-1}$ in (14), the orthogonality relation (64) and the definition of $\varphi _{k,r}$ in (13). This gives

$$\begin{aligned} t_\zeta ^{(k,r),e_i} = \frac{2 \zeta +1}{h_{e_i}} \left( \varphi _{k,r},m_\zeta ^{e_i}\right) _{0,e_i} = (2 \zeta +1) \delta _{i k} \delta _{r \zeta } \quad \forall \zeta =0,\dots ,p-1. \end{aligned}$$

The global system of linear equations corresponding to method (10) is assembled as in the standard non-conforming FEM. Finally, one imposes in a non-conforming fashion the Dirichlet boundary datum $g$ by

$$\begin{aligned} \int _eu_nq^e_{p-1}\, d s= \int _egq^e_{p-1}\, d s\quad \forall q^e_{p-1}\in \mathbb P_{p-1}(e), \end{aligned}$$

where, in practice, g is replaced by $g_p$, see Remark 2.