Stability and Convergence of Spectral Mixed Discontinuous Galerkin Methods for 3D Linear Elasticity on Anisotropic Geometric Meshes

TL;DR

This paper analyzes the stability and exponential convergence of spectral mixed discontinuous Galerkin methods for 3D linear elasticity problems on anisotropic meshes, especially near incompressible limits.

Contribution

It provides a computational study on the stability and convergence of mixed DG schemes for elasticity on complex meshes, including robustness near incompressibility.

Findings

Methods are stable with respect to the Poisson ratio.

Numerical evidence shows exponential convergence.

Schemes effectively handle singularities on polyhedral domains.

Abstract

We consider spectral mixed discontinuous Galerkin finite element discretizations of the Lam\'e system of linear elasticity in polyhedral domains in $\mathbb{R}^3$. In order to resolve possible corner, edge, and corner-edge singularities, anisotropic geometric edge meshes consisting of hexahedral elements are applied. We perform a computational study on the discrete inf-sup stability of these methods, and especially focus on the robustness with respect to the Poisson ratio close to the incompressible limit (i.e. the Stokes system). Furthermore, under certain realistic assumptions (for analytic data) on the regularity of the exact solution, we illustrate numerically that the proposed mixed DG schemes converge exponentially in a natural DG norm.

Equations

$$-\Delta \mathbf{u}+\nabla p$$

$$\nabla \cdot \mathbf{u}+(1-2\nu )p$$

$$\mathbf{u}$$

$$\begin{array}{rl}A(\mathbf{u},\mathbf{v})+B(\mathbf{v},p)& ={\int }_{\Omega }\mathbf{f}\cdot \mathbf{v}\mathsf{d}x,\\ -B(\mathbf{u},q)+C(p,q)& =0,\end{array}$$

$$A(\mathbf{u},\mathbf{v})$$

$$C(p,q)$$

$$a(\mathbf{u},p;\mathbf{v},q)={\int }_{\Omega }\mathbf{f}\cdot \mathbf{v}\mathsf{d}x\forall (\mathbf{v},q)\in H_0^1(\Omega {)}^3\times L_0^2(\Omega ),$$

$$a(\mathbf{u},p;\mathbf{v},q):=A(\mathbf{u},\mathbf{v})+B(\mathbf{v},p)-B(\mathbf{u},q)+C(p,q).$$

$$\underset{0\not\equiv q\in L_0^2(\Omega )}{\inf }\underset{0\not\equiv \mathbf{v}\in H_0^1(\Omega {)}^3}{\sup }\frac{-{\int }_{\Omega }q\nabla \cdot \mathbf{v}\mathsf{d}x}{\Vert \nabla \mathbf{v}{\Vert }_{L^2(\Omega )}\Vert q{\Vert }_{L^2(\Omega )}}\ge \kappa >0,$$

$$\mathcal{S}=\left(\bigcup _{\mathbf{c}\in \mathcal{C}}\mathbf{c}\right)\cup \left(\bigcup _{\mathbf{e}\in \mathcal{E}}\mathbf{e}\right)\subset \partial \Omega .$$

$$r_{\mathbf{c}}(x)=\mathrm{dist}(x,\mathbf{c}),r_{\mathbf{e}}(x)=\mathrm{dist}(x,\mathbf{e}),{\rho }_{ce}(x)=r_{\mathbf{e}}(x)/r_{\mathbf{c}}(x).$$

$$\begin{array}{rl}{\omega }_{\mathbf{c}}& =\{x\in \Omega :r_{\mathbf{c}}(x)<\epsilon \wedge {\rho }_{ce}(x)>\epsilon \forall \mathbf{e}\in {\mathcal{E}}_{\mathbf{c}}\},\\ {\omega }_{\mathbf{e}}& =\{x\in \Omega :r_{\mathbf{e}}(x)<\epsilon \wedge r_{\mathbf{c}}(x)>\epsilon \forall \mathbf{c}\in {\mathcal{C}}_{\mathbf{e}}\},\\ {\omega }_{\mathbf{c}\mathbf{e}}& =\{x\in \Omega :r_{\mathbf{c}}(x)<\epsilon \wedge {\rho }_{\mathbf{c}\mathbf{e}}(x)<\epsilon \}.\end{array}$$

$$\begin{split}|v|^{2}_{M^{m}_{\bm{\beta}}(\Omega)}=|v|^{2}_{H^{m}(\Omega_{0})}&+\sum_{{\bm{e}}\in\mathcal{E}}\sum_{\genfrac{}{}{0.0pt}{}{\bm{\alpha}\in\mathbb{N}^{3}_{0}}{|\bm{\alpha}|=m}}\big{\|}{r_{\bm{e}}^{\beta_{\bm{e}}+|\bm{\alpha}^{\perp}|}\mathsf{D}^{\bm{\alpha}}v}\big{\|}^{2}_{L^{2}(\omega_{{\bm{e}}})}\\ &+\sum_{{\bm{c}}\in\mathcal{C}}\sum_{\genfrac{}{}{0.0pt}{}{\bm{\alpha}\in\mathbb{N}^{3}_{0}}{|\bm{\alpha}|=m}}\big{\|}{r_{\bm{c}}^{\beta_{\bm{c}}+|\bm{\alpha}|}\mathsf{D}^{\bm{\alpha}}v}\big{\|}^{2}_{L^{2}(\omega_{{\bm{c}}})}\\ &+\sum_{{\bm{c}}\in\mathcal{C}}\sum_{{\bm{e}}\in\mathcal{E}_{\bm{c}}}\big{\|}{r_{\bm{c}}^{\beta_{\bm{c}}+|\bm{\alpha}|}\rho_{{\bm{c}}{\bm{e}}}^{\beta_{\bm{e}}+|\bm{\alpha}^{\perp}|}\mathsf{D}^{\bm{\alpha}}v}\big{\|}^{2}_{L^{2}(\omega_{{\bm{c}}{\bm{e}}})},\end{split}$$

$$0<{\beta }_{\mathbf{e}}<{\beta }_{\mathcal{E}},0<{\beta }_{\mathbf{c}}<{\beta }_{\mathcal{C}},\mathbf{e}\in \mathcal{E},\mathbf{c}\in \mathcal{C},$$

$$[[u]]$$

$$[[\mathbf{v}]]$$

$$[[\underline{w}]]$$

$${\int }_{\mathcal{F}}(\cdot )\mathsf{d}s:=\sum _{f\in \mathcal{F}}{\int }_f(\cdot )\mathsf{d}s,\Vert {\nabla }_h(\cdot ){\Vert }_{L^2(\Omega )}^2:=\sum _{K\in {\mathcal{T}}^{\sigma ,l}}\Vert \nabla (\cdot ){\Vert }_{L^2(K)}^2.$$

$$\begin{array}{rl}{\mathbf{V}}_h& :=\{\mathbf{v}\in L^2(\Omega {)}^3:\mathbf{v}{|}_K\in {\mathcal{Q}}_k(K{)}^3,K\in {\mathcal{T}}^{\sigma ,l}\},\\ Q_h& :=\{q\in L_0^2(\Omega ):q{|}_K\in {\mathcal{Q}}_{k-1}(K),K\in {\mathcal{T}}^{\sigma ,l}\}.\end{array}$$

$$\mathbf{V}(h)={\mathbf{V}}_h+H_0^1(\Omega {)}^3.$$

$$\mathtt{c}(x):=\gamma {\mathtt{h}}^{-1}(x)k^2,$$

$$\mathtt{h}(x):=\{\begin{array}{ll}\min \{h_{K^{♯},f}^{\perp },h_{K^{♭},f}^{\perp }\},& x\in f\subset {\mathcal{F}}_{\mathcal{I}},f=\partial K^{♯}\cap \partial K^{♭},\text{with }{K}^{♯},K^{♭}\in {\mathcal{T}}^{\sigma ,l},\\ h_{K,f}^{\perp },& x\in f\subset {\mathcal{F}}_{\mathcal{B}},f=\partial K\cap \partial \Omega ,\text{with }K\in {\mathcal{T}}^{\sigma ,l}.\end{array}$$

$$\begin{array}{rl}{A}_h({\mathbf{u}}_h,\mathbf{v})+B_h(\mathbf{v},p_h)& ={\int }_{\Omega }\mathbf{f}\cdot \mathbf{v}\mathsf{d}x,\\ -B_h({\mathbf{u}}_h,q)+C_h(p_h,q)& =0,\end{array}$$

$$A_h(\mathbf{u},\mathbf{v})$$

$$+{\int }_{\mathcal{F}}\mathtt{c}[[\mathbf{u}]]:[[\mathbf{v}]]\mathsf{d}s,$$

$$B_h(\mathbf{v},q)$$

$$C_h(p,q)$$

$$a_h({\mathbf{u}}_h,p_h;\mathbf{v},q)={\int }_{\Omega }\mathbf{f}\cdot \mathbf{v}\mathsf{d}x$$

$$a_h(\mathbf{u},p;\mathbf{v},q):=A_h(\mathbf{u},\mathbf{v})+B_h(\mathbf{v},p)-B_h(\mathbf{u},q)+C_h(p,q).$$

$$\begin{array}{rl}|||(\mathbf{v},q)||{|}_{\mathsf{DG}}^2:=& \Vert \mathbf{v}{\Vert }_h^2+(2-2\nu )\Vert q{\Vert }_{L^2(\Omega )}^2,\end{array}$$

$$\Vert \mathbf{v}{\Vert }_h^2:=\Vert {\nabla }_h\mathbf{v}{\Vert }_{L^2(\Omega )}^2+{\int }_{\mathcal{F}}\mathtt{c}|[[\mathbf{v}]]{|}^2\mathsf{d}s,\mathbf{v}\in \mathbf{V}(h).$$

$$a_h(\mathbf{u},p;\mathbf{u},p)\ge C\Vert \mathbf{u}{\Vert }_h^2+(1-2\nu )\Vert p{\Vert }_{L^2(\Omega )}^2\ge C(1-2\nu )|||(\mathbf{u},p)||{|}_{\mathsf{DG}}^2\forall (\mathbf{u},p)\in {\mathbf{V}}_h\times Q_h;$$

$${\gamma }_B:=\underset{0\not\equiv q\in Q_h}{\inf }\underset{0\not\equiv \mathbf{v}\in {\mathbf{V}}_h}{\sup }\frac{B_h(\mathbf{v},q)}{\Vert \mathbf{v}{\Vert }_h\Vert q{\Vert }_{L^2(\Omega )}}\ge \varkappa k^{-\rho },$$

$${\gamma }_a:=\underset{\genfrac{}{}{0ex}{}{(\mathbf{u},p)\in {\mathbf{V}}_h\times Q_h}{(\mathbf{u},p)\ne (0,0)}}{\inf }\underset{\genfrac{}{}{0ex}{}{(\mathbf{v},q)\in {\mathbf{V}}_h\times Q_h}{(\mathbf{v},q)\ne (0,0)}}{\sup }\frac{a_h(\mathbf{u},p;\mathbf{v},q)}{|||(\mathbf{u},p)||{|}_{\mathsf{DG}}|||(\mathbf{v},q)||{|}_{\mathsf{DG}}}\ge C\max \{k^{-2\rho },1-2\nu \},$$

Full text

1. Introduction

Consider an axi-parallel, open and bounded polyhedron $\Omega\subset\mathbb{R}^{3}$, with Lipschitz boundary $\partial\Omega$, in the three-dimensional Cartesian system. Using a spectral discontinuous Galerkin finite element method (DGFEM), we shall study the numerical approximation of the following linear elasticity problem in mixed form: Find a displacement field $\bm{u}=(u_{1},u_{2},u_{3})\in H^{1}_{0}(\Omega)^{3}$, and a pressure function $p\in L^{2}_{0}(\Omega)$ such that

[TABLE]

Here, $\nabla\cdot$ is the divergence operator, $\nu\in(0,\nicefrac{{1}}{{2}}]$ is the Poisson ratio, and $\bm{f}\in L^{2}(\Omega)^{3}$ is an external force (scaled by $\nicefrac{{2(1+\nu)}}{{E}}$, where $E>0$ is Young’s modulus). We shall include the limit case $\nu=\nicefrac{{1}}{{2}}$ which corresponds to the Stokes equations of incompressible fluid flow.

Elliptic boundary value problems in three-dimensional polyhedral domains are well-known to exhibit isotropic corner and anisotropic edge singularities, as well as a combination of the both; see, e.g., [6, 16, 17]. In a recent series of papers [20, 21, 22] on the numerical approximation of the Poisson equation in 3d polyhedra the use of $hp$-version DG methods has been proposed (see also [15] for eigenvalue problems with singular potentials). These schemes provide a convenient framework to resolve anisotropic edge singularities on (irregular) geometrically and anisotropically refined meshes, whilst using high-order spectral elements in the interior. Furthermore, supposing that the data is sufficiently smooth, it has been proved in [21, 22] that exponential convergence rates for $hp$-DG methods can be achieved.

In our previous work [24], we have employed the approach [20, 21, 22] in order to apply the high-order mixed DG methods introduced in [12] to the three-dimensional framework. More precisely, we have analyzed high-order interior penalty (IP) DG methods (of uniform but arbitrarily high polynomial degree) for the numerical approximation of (1)–(3) on geometrically refined edge meshes. They can be seen as $hp$-methods with fixed and uniform polynomial degrees, or as spectral methods on locally refined meshes; in this paper they shall simply be termed spectral DGFEM. Incidentally, in contrast to classical IPDG methods, these DG schemes feature anisotropically scaled penalty terms which account for possible element anisotropies; see also [7]. A focal point of the article [24] has been to provide an inf-sup stability analysis for mixed IPDG schemes on anisotropic meshes. Our results, which are based in parts on [19], are explicit with respect to both the (uniform) polynomial degree and the Poisson ratio $\nu$; cf. also [14] in the context of spectral methods for the Stokes equations. In particular, for fixed (but arbitrarily high) polynomial degrees, our stability analysis proves that the behaviour of the mixed DG scheme remains robust as $\nu$ tends to the critical limit of $\nicefrac{{1}}{{2}}$ of incompressible materials. Furthermore, following the techniques presented in [21, 22, 23] we showed that the proposed DG schemes are able to achieve exponential rates of convergence for the class of piecewise analytic functions in weighted Sobolev spaces studied in [6].

The goal of the present paper is to provide a computational investigation of the theoretical inf-sup stability results on anisotropic geometric edge meshes presented in [19, 24]. Furthermore, we will confirm the asserted exponential convergence of the spectral mixed DG method for a number of examples with typical edge and corner singularities in polyhedral domains. We will also look into the robustness of the DG approach with respect to the Poisson ratio as $\nu\to\nicefrac{{1}}{{2}}$. The precise outline of the paper is as follows: In Section 2, we present the mixed formulation of (1)–(3), and recall its regularity in terms of anisotropically weighted Sobolev spaces. Furthermore, in Section 3 the geometric edge meshes and the spectral mixed DG discretizations will be introduced; in addition, in Section 3.4, we revisit a discrete inf-sup stability framework together with the well-posedness of the DG scheme on anisotropic geometric edge meshes. Then, in Section 4, we discuss the practical computation of the DG solution as well as of the discrete inf-sup constants, and present some numerical results in a few typical reference situations. In Section 5 we perform a number of experiments which confirm the exponential convergence as well as the robustness of the DG method with respect to the Poisson ratio. Finally, we add a few concluding remarks in Section 6.

Notation

Throughout this article, we use the following notation: For a domain $D\subset{\mathbb{R}^{d}}$, $d\geq 1$, let $L^{2}(D)$ signify the Lebesgue space of square-integrable functions equipped with the usual norm ${\|\cdot\|_{L^{2}(D)}}$. Furthermore, we write $L^{2}_{0}(D)$ to denote the subspace of $L^{2}(D)$ of all functions with vanishing mean value on $D$. The standard Sobolev space of functions with integer regularity exponent $s\geq 0$ is signified by $H^{s}(D)$; we write $\|\cdot\|_{H^{s}(D)}$ and $|\cdot|_{H^{s}(D)}$ for the corresponding norms and semi-norms, respectively. As usual, we define $H^{1}_{0}(D)$ as the subspace of functions in $H^{1}(D)$ with zero trace on $\partial D$. For vector- and tensor-valued functions we use the standard notation $(\nabla\bm{v})_{ij}:=\partial_{j}v_{i}$, and $(\nabla\cdot\underline{\sigma})_{i}:=\sum_{j=1}^{3}\,\partial_{j}\sigma_{ij}$ and $\underline{\sigma}:\underline{\tau}:=\sum_{i,j=1}^{3}\,\sigma_{ij}\tau_{ij}$, respectively. Moreover, for vectors $\bm{v},\bm{w}\in\mathbb{R}^{3}$, let $\bm{v}\otimes\bm{w}\in\mathbb{R}^{3\times 3}$ be the matrix whose $ij$th component is $v_{i}\,w_{j}$.

2. Linear elasticity in polyhedra

We discuss the weak formulation of the mixed system of linear elasticity (1)–(3). Furthermore, we review its regularity in polyhedral domains.

2.1. Mixed formulation and well-posedness

A standard mixed formulation of (1)–(3) is to find $(\bm{u},p)\in H^{1}_{0}(\Omega)^{3}\times L^{2}_{0}(\Omega)$ such that

[TABLE]

for all $(\bm{v},q)\in H^{1}_{0}(\Omega)^{3}\times L^{2}_{0}(\Omega)$, where

[TABLE]

More compactly, we can write (4) equivalently in the form

[TABLE]

with

[TABLE]

It is straightforward to verify that $a$ is a bounded bilinear form on $H^{1}_{0}(\Omega)^{3}\times L^{2}_{0}(\Omega)$, and that, for $\nu\in(0,\nicefrac{{1}}{{2}})$, it is also coercive. In particular, by application of the Lax-Milgram theorem, we conclude that the solution $(\bm{u},p)\in H^{1}_{0}(\Omega)^{3}\times L^{2}_{0}(\Omega)$ of (5), and, thus, of (4), exists and is unique. In addition, for $\nu=\nicefrac{{1}}{{2}}$, which corresponds to the Stokes equations, problem (5) is still well-posed. Indeed, this is an immediate consequence of the inf-sup condition

[TABLE]

where $\kappa$ is the inf-sup constant, depending only on $\Omega$; see [5, 8] for details.

2.2. Regularity

Following [6] we recall the regularity of the solution of (1)–(3) in weighted Sobolev spaces; cf. also [9, 10, 11]. To this end, we denote by $\mathcal{C}$ the set of corners, and by $\mathcal{E}$ the set of edges of $\Omega$. Potential singularities of the solution are located on the skeleton of $\Omega$ given by

[TABLE]

Given a corner ${\bm{c}}\in\mathcal{C}$, an edge ${\bm{e}}\in\mathcal{E}$, and $x\in\Omega$, we define the distance functions

[TABLE]

Furthermore, for each corner ${\bm{c}}\in\mathcal{C}$, we signify by $\mathcal{E}_{{\bm{c}}}=\left\{\,{\bm{e}}\in\mathcal{E}\,:\,{\bm{c}}\cap\overline{\bm{e}}\neq\emptyset\,\right\}$ the set of all edges of $\Omega$ which meet at ${\bm{c}}$. For any ${\bm{e}}\in\mathcal{E}$, the set of corners of ${\bm{e}}$ is given by $\mathcal{C}_{\bm{e}}=\left\{\,{\bm{c}}\in\mathcal{C}\,:\,\,{\bm{c}}\cap\overline{\bm{e}}\neq\emptyset\,\right\}$. Then, for ${\bm{c}}\in\mathcal{C}$, ${\bm{e}}\in\mathcal{E}$, respectively ${\bm{e}}\in\mathcal{E}_{{\bm{c}}}$, and a sufficiently small parameter $\varepsilon>0$, we define the neighbourhoods

[TABLE]

Moreover, we define the interior part of $\Omega$ by $\Omega_{0}=\{x\in\Omega:\,\mathrm{dist}(x,\partial\Omega)>\varepsilon\}$.

Near corners ${\bm{c}}\in\mathcal{C}$ and edges ${\bm{e}}\in\mathcal{E}$, we shall use local coordinate systems in $\omega_{{\bm{e}}}$ and $\omega_{{\bm{c}}{\bm{e}}}$, which are chosen such that ${\bm{e}}$ corresponds to the direction $(0,0,1)$. Then, we denote quantities that are transversal to ${\bm{e}}$ by $(\cdot)^{\perp}$, and quantities parallel to ${\bm{e}}$ by $(\cdot)^{\|}$. For instance, if $\bm{\alpha}\in\mathbb{N}^{3}_{0}$ is a multi-index associated with the three local coordinate directions in $\omega_{{\bm{e}}}$ or $\omega_{{\bm{c}}{\bm{e}}}$, then we write $\bm{\alpha}=(\bm{\alpha}^{\perp},\alpha^{\|})$, where $\bm{\alpha}^{\perp}=(\alpha_{1},\alpha_{2})$ and $\alpha^{\|}=\alpha_{3}$. In addition, we use the notation $|\bm{\alpha}^{\perp}|=\alpha_{1}+\alpha_{2}$, and $|\bm{\alpha}|=|\bm{\alpha}^{\perp}|+\alpha^{\|}$.

Following [6, Def. 6.2 and Eq. (6.9)], we introduce anisotropically weighted Sobolev spaces. To this end, to each ${\bm{c}}\in\mathcal{C}$ and ${\bm{e}}\in\mathcal{E}$, we associate a corner and an edge exponent $\beta_{\bm{c}},\beta_{\bm{e}}\in\mathbb{R}$, respectively. We collect these quantities in the weight vector $\bm{\beta}=\{\beta_{\bm{c}}:\,{\bm{c}}\in\mathcal{C}\}\cup\{\beta_{\bm{e}}:\,{\bm{e}}\in\mathcal{E}\}\in\mathbb{R}^{|\mathcal{C}|+|\mathcal{E}|}$. Then, for $m\in\mathbb{N}_{0}$, we define the weighted semi-norm

[TABLE]

as well as the full norm $\left\|v\right\|^{2}_{M^{m}_{\bm{\beta}}(\Omega)}=\sum_{k=0}^{m}\|v\|^{2}_{M^{k}_{\bm{\beta}}(\Omega)}$; here, the operator $\mathsf{D}^{\bm{\alpha}}$ denotes the partial derivative in the local coordinate directions corresponding to the multi-index $\bm{\alpha}$. The space $M^{m}_{\bm{\beta}}(\Omega)$ is the weighted Sobolev space obtained as the closure of $C^{\infty}_{0}(\Omega)$ with respect to the norm $\left\|\cdot\right\|_{M^{m}_{\bm{\beta}}(\Omega)}$.

We notice the following regularity property of the solution of (1)–(3) in terms of the weighted Sobolev spaces defined above; see [6] (in addition, cf. [16, 17]):

Proposition 2.1**.**

There exist upper bounds $\beta_{\mathcal{E}},\beta_{\mathcal{C}}>0$ such that, if the weight vector $\bm{\beta}$ satisfies

[TABLE]

then, for every $m\in\mathbb{N}$, the solution $(\bm{u},p)\in H_{0}^{1}(\Omega)^{3}\times L^{2}_{0}(\Omega)$ of (1)–(3) with $\bm{f}\in M^{m}_{1-\bm{\beta}}(\Omega)^{3}$ fulfills $(\bm{u},p)\in M^{m}_{-1-\bm{\beta}}(\Omega)^{3}\times M^{m-1}_{-\bm{\beta}}(\Omega)$.

3. Mixed discontinuous Galerkin methods on geometric meshes

In the following section we will introduce spectral mixed DG discretizations on geometric meshes for the numerical solution of (1)–(3).

3.1. Hexahedral geometric edge meshes

In order to numerically resolve possible corner and edge singularities in the solution $(\bm{u},p)$ of (1)–(3), we employ anisotropic geometric edge meshes. To this end, we follow the construction in [19], where such meshes have been studied in the context of DGFEM for the Stokes equations; see also the earlier paper [2] on conforming $hp$-version finite element methods. Specifically, we begin from a coarse regular and shape-regular, quasi-uniform partition $\mathcal{T}^{0}=\{Q_{j}\}_{j=1}^{J}$ of $\Omega$ into $J$ convex axi-parallel hexahedra. Each of these elements $Q_{j}\in\mathcal{T}^{0}$ is the image under an affine mapping $G_{j}$ of the reference patch $\widetilde{Q}=(-1,1)^{3}$, i.e. $Q_{j}=G_{j}(\widetilde{Q})$. The mappings $G_{j}$ are compositions of (isotropic) dilations and translations.

Based on the coarse partition (macro mesh) $\mathcal{T}^{0}$ we will use three canonical geometric refinements (patches) towards corners, edges and corner-edges of $\widetilde{Q}$; see Figure 1. They feature a refinement ratio $\sigma\in(0,1)$, as well as a number of refinement levels $\ell\in\mathbb{N}_{0}$; to give an example, in Figure 1, we have selected $\sigma=\nicefrac{{1}}{{2}}$, and $\ell=3$.

Given a (fixed) refinement ratio $\sigma\in(0,1)$ as well as a refinement level value $\ell\in\mathbb{N}_{0}$, geometric meshes in $\Omega$ are now built by applying the patch mappings $G_{j}$ to transform the above canonical geometric mesh patches on the reference patch $\widetilde{Q}$ to the macro-elements $Q_{j}\in\mathcal{T}^{0}$, thereby yielding a local patch mesh $\mathcal{M}^{\sigma,\ell}_{j}$ on $Q_{j}$. The patches $Q_{j}$ away from the singular support $\mathcal{S}$ (i.e. with $\overline{Q}_{j}\cap\mathcal{S}=\emptyset$) are left unrefined, i.e. in this case we let $\mathcal{M}^{\sigma,\ell}_{j}=\{Q_{j}\}$. It is important to note that the geometric refinements in the canonical patches have to be suitably selected, oriented and combined in order to achieve a proper geometric refinement towards corners and edges of $\Omega$. Then, a $\sigma$-geometric mesh in $\Omega$ is given by $\mathcal{T}^{\sigma,l}=\bigcup_{j=1}^{J}\mathcal{M}^{\sigma,\ell}_{j}$. Furthermore, the sequence $\{\mathcal{T}^{\sigma,l}\}_{\ell\in\mathbb{N}_{0}}$ is referred to as a $\sigma$-geometric mesh family. We note that this family of meshes is anisotropic as well as irregular. For a more general construction of geometric meshes on polyhedral domains, we refer to [20].

3.2. Faces and face operators

We denote the set of all interior faces in $\mathcal{T}^{\sigma,l}$ by $\mathcal{F}_{\mathcal{I}}(\mathcal{T}^{\sigma,l})$, and the set of all boundary faces by $\mathcal{F}_{\mathcal{B}}(\mathcal{T}^{\sigma,l})$. Further, let $\mathcal{F}(\mathcal{T}^{\sigma,l})=\mathcal{F}_{\mathcal{I}}(\mathcal{T}^{\sigma,l})\cup\mathcal{F}_{\mathcal{B}}(\mathcal{T}^{\sigma,l})$ signify the set of all (smallest) faces of $\mathcal{T}^{\sigma,l}$. In addition, for an element $K\in\mathcal{T}^{\sigma,l}$, we denote the set of its faces by $\mathcal{F}_{K}=\left\{\,f\in\mathcal{F}\,:\,f\subset\partial K\,\right\}$.

Next, we recall the standard DG trace operators. For this purpose, consider an interior face $f=\partial K^{\sharp}\cap\partial K^{\flat}\in\mathcal{F}_{\mathcal{I}}(\mathcal{T}^{\sigma,l})$ shared by two neighbouring elements $K^{\sharp},K^{\flat}\in\mathcal{T}^{\sigma,l}$. Furthermore, let $u$, $\bm{v}$ and $\underline{w}$ be scalar-, vector, and tensor-valued functions, respectively, all sufficiently smooth inside the elements $K^{\sharp},K^{\flat}$. Then, we define the following trace operators along $f$:

[TABLE]

Here, for an element $K\in\mathcal{T}^{\sigma,l}$, we denote by $\bm{n}_{K}$ the outward unit normal vector on $\partial K$. Similarly, for a boundary face $f=\partial K\cap\partial\Omega\in\mathcal{F}_{\mathcal{B}}(\mathcal{T}^{\sigma,l})$, with $K\in\mathcal{T}^{\sigma,l}$, and a sufficiently smooth scalar function $u$, we let $[\![u]\!]=u|_{K}\bm{n}_{\Omega}$, and $\{\!\!\{u\}\!\!\}=u|_{K}$, where $\bm{n}_{\Omega}$ is the outward unit normal vector on $\partial\Omega$; obvious modifications are made for vector- and tensor-valued functions in accordance with the definition above.

Finally, $\nabla_{h}$ and $\nabla_{h}\cdot$ denote the element-wise gradient and divergence operators, respectively. Here and in the sequel, we use abbreviations like

[TABLE]

3.3. Spectral DG discretizations

Given a geometric edge mesh $\mathcal{T}^{\sigma,l}$ on $\Omega$ and a polynomial degree $k\geq 1$ (which is assumed uniform and isotropic on $\mathcal{T}^{\sigma,l}$), we approximate (1)–(3) by finite element functions $(\bm{u}_{h},p_{h})\in\bm{V}_{h}\times Q_{h}$, where

[TABLE]

Here, for $k\geq 0$, $K\in\mathcal{T}^{\sigma,l}$, ${\mathcal{Q}}_{k}(K)$ denotes the space of all polynomials of degree at most $k$ in each variable on $K$. In addition, we let

[TABLE]

On the space $\bm{V}_{h}$ we consider the stabilization function ${\tt c}\in L^{\infty}(\mathcal{F})$ given by

[TABLE]

where $\gamma>0$ is a penalty parameter independent of the refinement ratio $\sigma$, the number of refinement levels $\ell$, and the polynomial degree $k$. Furthermore, for $x\in f$, with $f\in\mathcal{F}$, the mesh function ${\tt h}$ is defined by

[TABLE]

In this definition, for $K\in\mathcal{T}^{\sigma,l}$ and $f\in\mathcal{F}_{K}$, we denote by $h_{K,f}^{\perp}$ the diameter of the element $K$ in the direction perpendicular to the face $f$.

Then, we consider the following mixed discontinuous Galerkin discretization of (4): Find $(\bm{u}_{h},p_{h})\in\bm{V}_{h}\times Q_{h}$ such that

[TABLE]

for all $(\bm{v},q)\in\bm{V}_{h}\times Q_{h}$. The forms $A_{h}$, $B_{h}$, and $C_{h}$ are given, respectively, by

[TABLE]

where $\theta\in[-1,1]$ is a fixed parameter. Different choices of $\theta$ refer to various types of interior penalty DG methods: for instance, the form $A_{h}$ may be chosen to correspond to the symmetric (for $\theta=1$), incomplete (for $\theta=0$), or non-symmetric (for $\theta=-1$) interior penalty DG discretization of the Laplacian; for a detailed review on a wide class of DG methods for the Poisson problem and the Stokes system, we refer to the articles [1, 18], respectively.

As in (5), the discrete DG formulation (8) is equivalent to finding $(\bm{u}_{h},p_{h})\in\bm{V}_{h}\times Q_{h}$ such that

[TABLE]

for all $(\bm{v},q)\in\bm{V}_{h}\times Q_{h}$, where

[TABLE]

3.4. Discrete inf-sup stability and well-posedness

In this section we recapitulate an inf-sup stability result from [24] for the form $a_{h}$ given in (11). We first define the DG-norm

[TABLE]

for any $(\bm{v},q)\in\bm{V}(h)\times L^{2}(\Omega)$, where

[TABLE]

If the Poisson ratio satisfies $\nu\in(0,\nicefrac{{1}}{{2}})$, and provided that the penalty parameter $\gamma$ featured in (7) is chosen sufficiently large, then the following coercivity estimate can be shown:

[TABLE]

here, $C>0$ is a constant independent of $\nu$, $k$, $l$, and the aspect ratio of the anisotropic elements.

This result can be made stronger if a discrete inf-sup condition on the form $B_{h}$ on geometric edge meshes is assumed: Let $\mathcal{T}^{\sigma,l}$ be a geometric edge mesh on $\Omega$ as defined in Section 3.1, with refinement ratio $\sigma\in(0,1)$, and $\ell\geq 1$ layers of refinement. Suppose that there exist constants $\varkappa>0$ and $\rho\geq 0$ that may depend on $\sigma$, $\gamma$, and on the macro-element mesh $\mathcal{T}^{0}$, but are independent of $k$, $\ell$, and the aspect ratio of the anisotropic elements in $\mathcal{T}^{\sigma,l}$, such that there holds

[TABLE]

as $k\to\infty$.

Remark 3.1**.**

In [19] it was proved that this assumption is fulfilled with $\rho=\nicefrac{{3}}{{2}}$ (and any $k\geq 2$). Our numerical computations in Section 4.4.1 below indicate, however, that the dependence of the right-hand side of (13) on $k$ is much weaker than $k^{-\nicefrac{{3}}{{2}}}$.

The following result, which implies the well-posedness of (8) even in the incompressible limit $\nu=\nicefrac{{1}}{{2}}$, follows immediately from [24, Theorem 5.1].

Theorem 3.2**.**

Let $\nu\in(0,\nicefrac{{1}}{{2}}]$. If (13) holds true, then we have the inf-sup condition

[TABLE]

with a constant $C>0$ that depends on the penalty parameter $\gamma$, however, is independent of $\nu$, $k$, $l$, and the aspect ratio of the anisotropic elements.

We emphasize that, for fixed $k$, the stability bound (14) does not deteriorate as $\nu\to\nicefrac{{1}}{{2}}$.

4. Computing the DG solution and the inf-sup constants

Our goal is to investigate the behavior of the inf-sup conditions from (13) and (14) numerically. We note that both of them involve the discrete space $Q_{h}$ from (6). Due to the global zero mean constraint contained in $L^{2}_{0}(\Omega)$, the construction of $Q_{h}$ in terms of standard local basis functions, as provided by most finite element packages, causes difficulties. The classical remedy is to use the full space

[TABLE]

and to impose the zero mean condition by means of a Lagrange multiplier technique. Noticing that $\dim(Q_{h})=\dim(\widetilde{Q}_{h})-1$, we emphasize that this approach is of equivalent computational cost as the original DG system (8), yet, it allows to employ standard discrete spaces. This, in turn, leads to a more convenient practical framework for the computation of the DG solution and the evaluation of inf-sup constants.

4.1. Reformulation of the mixed DG discretization

We rewrite the original system (8), which is based on the discrete DG space $\bm{V}_{h}\times Q_{h}$, on the new space $\bm{V}_{h}\times\widetilde{Q}_{h}\times\mathbb{R}$, where $\widetilde{Q}_{h}$ is the full space from (15). To this end, we introduce an auxiliary variable $\widetilde{r}\in\mathbb{R}$ which takes the role of the mean value of the pressure $p_{h}$ on $\Omega$. More precisely, let us consider the following augmented DG formulation: Find $(\widetilde{\bm{u}}_{h},\widetilde{p}_{h},\widetilde{r})\in\bm{V}_{h}\times\widetilde{Q}_{h}\times\mathbb{R}$ such that

[TABLE]

for all $(\widetilde{\bm{v}},\widetilde{q},\widetilde{s})\in\bm{V}_{h}\times\widetilde{Q}_{h}\times\mathbb{R}$. Here we use the notation

[TABLE]

to denote the mean value integral on $\Omega$. Note also that this new system may be written in a more compact way: Find $(\widetilde{\bm{u}}_{h},\widetilde{p}_{h},\widetilde{r})\in\bm{V}_{h}\times\widetilde{Q}_{h}\times\mathbb{R}$ such that

[TABLE]

where

[TABLE]

We again stress the fact that this system can be expressed in terms of standard local basis functions, and, thereby, permits to apply a straightforward implementational setting.

To show the equivalence of the two formulations (8) and (16), we require the following lemma. Here, we shall denote by $\mathcal{Q}_{0}(\Omega)\simeq\mathbb{R}$ the space of all (globally) constant functions on $\Omega$, and we note that

[TABLE]

Lemma 4.1**.**

The form $B_{h}$ from (9) satisfies

[TABLE]

Conversely, if, for given $q\in\widetilde{Q}_{h}$, there holds that $B_{h}(\bm{v},q)=0$ for any $\bm{v}\in\bm{V}_{h}$, then it follows that $q\in\mathcal{Q}_{0}(\Omega)$.

Proof.

Given $(\bm{v},q)\in\bm{V}_{h}\times\mathcal{Q}_{0}(\Omega)$. Since $\mathcal{Q}_{0}(\Omega)$ is one-dimensional we may, without loss of generality, suppose that $q\equiv 1$. Then, we have

[TABLE]

By applying the Gauss-Green theorem on each element $K\in\mathcal{T}^{\sigma,l}$, we obtain two expressions which are identical, and, thus, cancel out:

[TABLE]

Let now $q\in\widetilde{Q}_{h}$, with $q-\fint_{\Omega}q\,{\mathsf{d}}x\neq 0$, and $B_{h}(\bm{v},q)=0$ for all $\bm{v}\in\bm{V}_{h}$. Then, the inf-sup condition (13) implies that

[TABLE]

Applying (18) yields a contradiction, and, consequently, we deduce that $q-\fint_{\Omega}q\,{\mathsf{d}}x=0$. Thus, we have $q\equiv\fint_{\Omega}q\,{\mathsf{d}}x\in\mathcal{Q}_{0}$. This completes the proof. ∎

Now we can state the equivalence of the two formulations.

Proposition 4.2**.**

The augmented DG discretization from (16) has a unique solution $(\widetilde{\bm{u}}_{h},\widetilde{p}_{h},0)\in\bm{V}_{h}\times\widetilde{Q}_{h}\times\mathbb{R}$, and $(\widetilde{\bm{u}}_{h},\widetilde{p}_{h})$ is the solution of the original DG formulation (8) with $\widetilde{p}_{h}\in Q_{h}$.

Proof.

We proceed in three steps.

Step 1:

We first show that the new formulation enforces the pressure $\widetilde{p}_{h}$ to have zero mean. To this end, we choose the test variable to be $\widetilde{s}=1$. Thus, from the third equation of (16), we deduce that

[TABLE]

Furthermore, let us choose $\widetilde{q}\equiv 1\in\mathcal{Q}_{0}$. From the second equation of (16), and with the aid of (18), we infer that

[TABLE]

i.e. $\fint_{\Omega}\widetilde{p}_{h}\,{\mathsf{d}}x=\widetilde{r}=0$, since $\nu\neq 0$.

Step 2:

Next, we show that $(\widetilde{\bm{u}}_{h},\widetilde{p}_{h},\widetilde{r}):=(\bm{u}_{h},p_{h},0)\in\bm{V}_{h}\times Q_{h}\times\mathbb{R}$, where $(\bm{u}_{h},p_{h})$ is the solution from (8), solves (16). The first and last equation in (16) are clearly fulfilled with $(\widetilde{\bm{u}}_{h},\widetilde{p}_{h},\widetilde{r})=(\bm{u}_{h},p_{h},0)$. The second equation in (16) simplifies to

[TABLE]

To show that this equation does indeed hold true, we notice that

[TABLE]

Due to (18), the first term on the right-hand side is zero. For the second term, since $\widetilde{q}-\fint_{\Omega}\widetilde{q}\in Q_{h}$, it holds

[TABLE]

simply by the second equation in (8). Hence, we end up with

[TABLE]

where we have used the fact that $p_{h}$ has zero mean. Thus, all three equations from (16) are fulfilled with $(\widetilde{\bm{u}}_{h},\widetilde{p}_{h},\widetilde{r})=(\bm{u}_{h},p_{h},0)$ from above.

Step 3:

It remains to show that the solution of (16) is unique. To this end, we assume that there exist two solutions $(\widetilde{\bm{u}}_{h1},\widetilde{p}_{h1},0),(\widetilde{\bm{u}}_{h2},\widetilde{p}_{h2},0)\in\bm{V}_{h}\times Q_{h}\times\mathbb{R}$. Using again (18) as well as the fact that the pressures have zero mean, the second equation in (16) implies that

[TABLE]

Therefore, we get the following system for the difference $(\widetilde{\bm{u}}_{h1}-\widetilde{\bm{u}}_{h2},\widetilde{p}_{h1}-\widetilde{p}_{h2})\in\bm{V}_{h}\times Q_{h}$ of the two solutions:

[TABLE]

for all $(\widetilde{\bm{v}},\widetilde{q}-\fint_{\Omega}\widetilde{q}\,{\mathsf{d}}x)\in\bm{V}_{h}\times Q_{h}$. This, in turn, is just the mixed formulation (8) with the unique zero solution. ∎

4.2. Inf-sup constant of the form $B_{h}$

We will now discuss how to compute the inf-sup constant $\gamma_{B}$ from (13) numerically. To do so, we proceed along the lines of [5, Chapter II.3.2]. Let us choose two sets of basis functions $\{\bm{\phi}_{i}\}_{i=1}^{M}\subset\bm{V}_{h}$ and $\{\psi_{j}\}_{j=1}^{N}\subset\widetilde{Q}_{h}$, with $M=\dim(\bm{V}_{h})$ and $N=\dim(\widetilde{Q}_{h})$. For any $\bm{v}\in V_{h}$ and $q\in\widetilde{Q}_{h}$, we store the associated coefficients in two vectors $\bm{\mathsf{v}}:=(v_{1},\dots,v_{M})^{\top}$ and $\bm{\mathsf{q}}:=(q_{1},\dots,q_{N})^{\top}$, respectively. Counting degrees of freedom in each element, we remark that

[TABLE]

Furthermore, we define the matrix $\bm{\mathsf{B}}\in\mathbb{R}^{M\times N}$ by

[TABLE]

Due to Lemma 4.1, we notice that $q\in\mathcal{Q}_{0}$ if and only if the associated coefficient vector $\bm{\mathsf{q}}$ satisfies $\bm{\mathsf{q}}\in\ker(\bm{\mathsf{B}})$. Moreover, by virtue of (17), we conclude that $q\in Q_{h}$ if and only if $\bm{\mathsf{q}}\in\mathbb{R}^{N}/\ker(\bm{\mathsf{B}})$. Moreover, since $\dim(\mathcal{Q}_{0})=1$, it follows that $\dim(\ker(\bm{\mathsf{B}}))=1$, and, in view of (19),

[TABLE]

Let us further introduce the symmetric positive definite matrices $\bm{\mathsf{D}}\in\mathbb{R}^{M\times M}$ and $\bm{\mathsf{E}}\in\mathbb{R}^{N\times N}$ corresponding to the norms $\|\cdot\|_{h}$ and $\|\cdot\|_{L^{2}(\Omega)}$, respectively, through

[TABLE]

Taking into account our considerations above, we infer that

[TABLE]

To proceed, we define

[TABLE]

Let the singular value decomposition (SVD) of $\widetilde{\bm{\mathsf{B}}}$ be given by

[TABLE]

with the singular values $\sigma_{1}\geq\dots\geq\sigma_{N-1}>\sigma_{N}=0$, cf. (20), being contained on the diagonal of $\bm{\mathsf{\Sigma}}$, and orthogonal matrices $\widetilde{\bm{\mathsf{V}}}$ and $\widetilde{\bm{\mathsf{Q}}}$ with columns $\widetilde{\bm{\mathsf{v}}}_{1},\dots,\widetilde{\bm{\mathsf{v}}}_{M}$ and $\widetilde{\bm{\mathsf{q}}}_{1},\dots,\widetilde{\bm{\mathsf{q}}}_{N}$, respectively. In addition, we set

[TABLE]

For $i\in\{1,\dots,N\}$ we then conclude

[TABLE]

where $\bm{\mathsf{e}}_{i}$ is the $i$th standard unit vector in $\mathbb{R}^{N}$. Hence, we have

[TABLE]

where $\widetilde{\bm{\mathsf{e}}}_{i}$ is the $i$th standard unit vector in $\mathbb{R}^{M}$. Involving (22), we deduce that

[TABLE]

Given linear combinations

[TABLE]

and using (23) and (24), we obtain

[TABLE]

Moreover, employing (24) and (25), the norms are represented by

[TABLE]

We will now evaluate the inf-sup constant

[TABLE]

with $\bm{\alpha}:=(\alpha_{1},\dots,\alpha_{N-1})$ and $\bm{\beta}:=(\beta_{1},\dots,\beta_{N-1})$. Without loss of generality, we may suppose that \|\bm{v}\|_{h}=\big{(}\sum_{i=1}^{N-1}\alpha_{i}^{2}\big{)}^{\nicefrac{{1}}{{2}}}=1, and \|q\|_{L^{2}(\Omega)}=\big{(}\sum_{i=1}^{N-1}\beta_{i}^{2}\big{)}^{\nicefrac{{1}}{{2}}}=1. Then, (26) simplifies to $\gamma_{B}=\inf_{\bm{\beta}}\sup_{\bm{\alpha}}\sum_{i=1}^{N-1}\sigma_{i}\alpha_{i}\beta_{i}$. By applying the Cauchy-Schwarz inequality, that is,

[TABLE]

we observe that the supremum is attained for \alpha_{i}=\big{(}\sum_{i=1}^{N-1}\sigma_{i}^{2}\beta_{i}^{2}\big{)}^{-\nicefrac{{1}}{{2}}}\sigma_{i}\beta_{i}. This leads to

[TABLE]

Proposition 4.3**.**

The inf-sup constant $\gamma_{B}$ from (13) is given by the smallest positive singular value, $\sigma_{N-1}$, of the matrix $\widetilde{\bm{\mathsf{B}}}$ from (21).

4.3. Inf-sup constant of the form $a_{h}$

In order to compute the inf-sup constant from (14), we proceed analogously as in the previous section. To this end, we choose a basis $\{\bm{\chi}_{i}\}_{i=1}^{M+N}$ of $\bm{V}_{h}\times\widetilde{Q}_{h}$, and define the system matrix $\bm{\mathsf{M}}\in\mathbb{R}^{(M+N)\times(M+N)}$ by

[TABLE]

with $a_{h}$ from (11). For brevity, we consider only the limit case $\nu=\nicefrac{{1}}{{2}}$. Due to Lemma 4.1 and the coercivity of the form $A_{h}$ from (8), we conclude that $\bm{\mathsf{M}}$ has a one-dimensional kernel. Denoting by $\bm{\mathsf{D}}\in\mathbb{R}^{(M+N)\times(M+N)}$ the symmetric positive matrix defining the $|\!|\!|\cdot|\!|\!|_{\mathsf{DG}}$-norm through

[TABLE]

where $\bm{\mathsf{v}}\in\mathbb{R}^{M+N}$ is the coefficient vector of a given pair $(\bm{u},p)\in\bm{V}_{h}\times\widetilde{Q}_{h}$ with respect to the basis $\{\bm{\chi}_{i}\}_{i=1}^{M+N}$, we let

[TABLE]

Given the SVD of $\widetilde{\bm{\mathsf{M}}}$ by $\widetilde{\bm{\mathsf{M}}}=:\widetilde{\bm{\mathsf{X}}}\bm{\mathsf{\Sigma}}\widetilde{\bm{\mathsf{Y}}}^{\top}$, with the singular values $\sigma_{1}\geq\dots\geq\sigma_{M+N-1}>\sigma_{M+N}=0$ being contained on the diagonal of $\bm{\mathsf{\Sigma}}$, and orthogonal matrices $\widetilde{\bm{\mathsf{X}}}$ and $\widetilde{\bm{\mathsf{Y}}}$, we infer the following result.

Proposition 4.4**.**

In the incompressible case $\nu=\nicefrac{{1}}{{2}}$, the inf-sup constant from (14) satisfies $\gamma_{a}=\sigma_{M+N-1}$, where $\sigma_{M+N-1}$ is the smallest positive singular value of the matrix $\widetilde{\bm{\mathsf{M}}}$ from (27).

4.4. Numerical computation of the inf-sup constants

We shall now investigate the inf-sup constants $\gamma_{B}$ and $\gamma_{a}$ from (13) and (14), respectively, by means of a number of numerical experiments. In particular, we will investigate the dependence on the approximation degree $k$ and on the Poisson ratio $\nu$. In the sequel, we choose $\theta$ from (9) to be 1 (i.e. we use the symmetric interior penalty DG method), the mesh grading factor from Section 3.1 as $\sigma=\nicefrac{{1}}{{2}}$, and the penalty parameter from (7) is set to $\gamma=10$. As shape functions we use tensorized Lagrange polynomials in the Gauss quadrature points. All our computations are performed with the finite element library deal.II; see, e.g., [3, 4].

4.4.1. Inf-sup constant $\gamma_{B}$

We consider the canonical edge, corner, and corner-edge patch meshes presented in Section 3.1, see Figure 1, with the modification that, in case of the corner-edge mesh, we refine one corner and all adjacent neighbouring edges. Furthermore, we study the situation of geometrically refined meshes on a Fichera domain given by $(-1,1)^{3}\setminus[0,1)^{3}$, where we simultaneously refine the reentrant corner and all three adjacent edges.

Let us first fix the approximation degree $k$, and refine the meshes by increasing the number of refinement levels $\ell$ step by step. For different approximation degrees the results are depicted in Figure 2. We clearly observe that the values of $\gamma_{B}$ stabilize after some initial refinement steps, thereby underlining the robustness of the inf-sup constant of $B_{h}$ with respect to the anisotropic geometric refinements. The asymptotic values are visualized in Figure 3; it is observed that there is a mild $k$-dependence of the inf-sup constant $\gamma_{B}$, however, our results indicate that, for the given examples, the dependence is considerably more optimistic than the theoretical bound $k^{-\nicefrac{{3}}{{2}}}$ proved in [19], cf. Remark 3.1.

4.4.2. Inf-sup constant of the form ${a}_{h}$

Let us turn to the behavior of the inf-sup constant $\gamma_{a}$ from (14) with respect to $k$, with $\nu=\nicefrac{{1}}{{2}}$. From Theorem 3.2 recall the theoretical dependence $\gamma_{a}\gtrsim\max\{k^{-2\rho},1-2\nu\}$; this result holds true with a theoretical value of $\rho=\nicefrac{{3}}{{2}}$, cf. Remark 3.1. Since we set $\nu=\nicefrac{{1}}{{2}}$, we deduce $\gamma_{a}\gtrsim k^{-3}$.

We focus on the canonical geometric edge, corner, and corner-edge refinements from Section 3.1. As before, we first fix the approximation degree $k$, and refine the meshes step by step in order to monitor the inf-sup constant; see Figure 4 for the resulting plots with different approximation degrees. Again, we display the stabilized values of $\gamma_{a}$ for increasing $k$ in Figure 5. The results are qualitatively similar to the inf-sup constant $\gamma_{B}$ discussed earlier. There is a $k$-dependence of the inf-sup constant $\gamma_{a}$ which is again much weaker than $k^{-3}$. Furthermore, our results show that the inf-sup constant $\gamma_{a}$ does not deteriorate in the critical limit $\nu=\nicefrac{{1}}{{2}}$ as shown in Theorem 3.2.

5. Exponential convergence on geometric meshes

In this section we turn to the exponential convergence of the spectral mixed DG method (8). Inspired by the regularity theory from [6] for analytic data (cf., in particular, Theorem 6.9), we suppose that the solution $(\bm{u},p)$ of (1)–(3) belongs to $A_{-1-\bm{\beta}}(\Omega)^{3}\times A_{-\bm{\beta}}(\Omega)$, where, for a weight vector $\bm{\gamma}\in\mathbb{R}^{|\mathcal{C}|+|\mathcal{E}|}$, we consider the countably normed space of piecewise analytic functions

[TABLE]

with a constant $C_{v}>0$ depending on the function $v$. Under this assumption, referring to [24, Theorem 6.2], it can be shown that the DG approximation $(\bm{u}_{h},p_{h})$ from (8) converges at an exponential rate. To quantify this fact, let $\nu\in(0,\nicefrac{{1}}{{2}}]$, and consider a sequence of geometric edge meshes $\mathcal{T}^{\sigma,l}$ as in Section 3.1. Moreover, choose a uniform polynomial degree $k\geq 2$ that is proportional to the number of layers $\ell\geq 1$ in $\mathcal{T}^{\sigma,l}$. Then, the DG approximation $(\bm{u}_{h},p_{h})$ from (8) satisfies the error bound

[TABLE]

where $N:=\dim(\bm{V}_{h}\times Q_{h})$ denotes the number of degrees of freedom.

5.1. Exponential convergence on canonical patches

In our numerical examples, we use manufactured solutions, which feature typical singularities close to $\mathcal{S}$ in the displacement $\bm{u}=(u_{1},u_{2},u_{3})$, and then test the spectral DGFEM (8) with the resulting right-hand side force functions $\bm{f}$ in (1). The domain $\Omega$ is chosen to be the unit cube. To cover all possible cases, we consider one solution with an edge, one with a corner, and another one with a corner-edge singularity:

• •

Displacement with an anisotropic edge singularity along the $z$-axis:

[TABLE]

• •

Displacement with an isotropic corner singularity at the origin:

[TABLE]

• •

Displacement which combines the above edge and corner singularities:

[TABLE]

The corresponding pressures $p$ for these displacements are then given via (2):

[TABLE]

Incidentally, in order to avoid too complicated right-hand sides $\bm{f}$, we do not enforce the displacements $\bm{u}$ to vanish on the whole boundary $\partial\Omega$. More precisely, they are chosen such that \bm{u}\cdot\bm{n}\big{|}_{\partial\Omega}=0, i.e. the normal component of the displacement field is zero. Then, by the Gauss-Green theorem, notice the mean value property of the pressure,

[TABLE]

The nonzero Dirichlet boundary conditions are accounted for by means of a standard flux term which is added to the right-hand side of the DGFEM formulation  (10).

For our test examples, we use $\nu\in\left\{\nicefrac{{1}}{{8}},\nicefrac{{1}}{{2}},\nicefrac{{3}}{{8}}\right\}$. In order to solve the resulting linear systems we employ the GMRES method in combination with the SparseILU preconditioner implemented in deal.II. The iterations terminate as soon as the Euclidean norm of the (unpreconditioned) residual becomes smaller or equal to $10^{-12}$. The initial meshes consist of a single element, and an approximation degree $k=1$. In the following, we refine successively the meshes towards the singularities, and simultaneously increase the approximation degree by one in each refinement step such that $k\sim\ell$, where $\ell$ is the number of layers. Since the singularities (and thereby their location) in the examples above are known explicitly, we only refine the corresponding edge and/or corner; see Figure 6 for the corner-edge example.

In Figures 7 and 8 we display the error of the approximation in the DG-norm (12) in a semi-logarithmic coordinate system with respect to the 5th, respectively the 4th root of the number of degrees of freedom $N$; cf. (28). Indeed, on a single element, the number of degrees of freedom grows with $\mathcal{O}(k^{3})$; in addition, when resolving a corner-edge singularity, the number of elements grows with $\mathcal{O}(\ell^{2})$, while, for an edge or a corner singularity, only $\mathcal{O}(\ell)$ many elements are needed. Hence, recalling that $k\sim\ell$, in the cases of the edge and the corner singularity examples a growth of $\mathcal{O}(N^{4})$ degrees of freedom is obtained, while in the case of the corner-edge example we even have $\mathcal{O}(N^{5})$ (thus the 5th root in (28)). The graphs show that, after some initial refinement steps, we obtain nearly constant slopes in all three situations. Hence, these experiments confirm that the proposed spectral DGFEM (8) on geometric edge meshes is able to resolve isotropic as well as anisotropic singularities at exponential rates.

5.2. Robustness with respect to the Poisson ratio

The purpose of the second series of experiments is to investigate the robustness of the exponential convergence bound (28) with respect to $\nu$ as $\nu\to\nicefrac{{1}}{{2}}$. The domain $\Omega$ is again chosen to be the unit cube.

Example 5.1**.**

We first consider an example where the displacement $\bm{u}$ is smooth and divergence-free:

[TABLE]

In this case it immediately follows that $p=0$, and, hence, $-\Delta\bm{u}=\bm{f}$; in particular, the resulting right-hand side force function $\bm{f}$ is independent of $\nu$.

For this example, we use a fixed uniform mesh consisting of 64 elements, and simply vary the uniform polynomial degree $k$. For this setup, since the solution $(\bm{u},p)$ is analytic, we expect exponential convergence with respect to the 3rd root of the number of degrees of freedom. In order to study the convergence with respect $\nu$ as $\nu\to\nicefrac{{1}}{{2}}$, in Figure 9, we plot the error of the DG method with respect to different values of the Poisson ratio $\nu$. Clearly, the deviations between the different exponential convergence curves are almost negligible, and, thereby, underline the robustness of the DGFEM with respect to $\nu$ for this example.

Example 5.2**.**

In our last experiment, we choose a circular force in the $x$-$y$-plane and a linear force in the $z$-direction, i.e.

[TABLE]

Since the exact solution is not known in this example, we compute a reference solution based on refining all edges and corners of $\Omega$ with $k=\ell=5$; cf. Figure 10. The DG error for different values of $\nu$ is depicted in Figure 11; as in the previous example, we observe that the DGFEM remains stable when $\nu$ tends to the incompressible limit $\nicefrac{{1}}{{2}}$. Furthermore, the nearly straight graphs indicate that exponential convergence is also achieved in these computations.

6. Conclusions

This paper centres on spectral mixed discontinuous Galerkin discretizations for linear elasticity and Stokes flow in three dimensional polyhedral domains. In a series of numerical experiments we have validated our theoretical results from [24] for various canonical reference situations. In particular, we have performed a computational study on the inf-sup stability and the exponential convergence of this class of methods on anisotropic geometric edge meshes. For the former purpose we have derived a simple procedure to determine the discrete inf-sup constants based on a singular value decomposition approach along the lines of [5].

Following the approach [20, 21, 22], our work may be extended to variable (and possibly anisotropic) polynomial degree distributions, thereby leading to $hp$-version DGFEM. To prove stability results in that context, however, the discrete inf-sup condition (13) would need to be generalized to the corresponding $hp$-meshes. Finally, let us mention that the linear mixed DG discretizations addressed in the present work could be combined with the iterative Newton DG (NDG) approach [13] in order to approximate problems in nonlinear elasticity in three dimensions.