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

We consider spectral mixed discontinuous Galerkin finite element discretizations of the Lamé system of linear elasticity in polyhedral domains in 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.

## Authors

• 8 publications
• 1 publication
04/23/2020

### Analytic regularity for the incompressible Navier-Stokes equations in polygons

In a plane polygon P with straight sides, we prove analytic regularity o...
12/28/2017

### Spectral Methods in the Presence of Discontinuities

Spectral methods provide an elegant and efficient way of numerically sol...
12/21/2020

### A p-robust polygonal discontinuous Galerkin method with minus one stabilization

We introduce a new stabilization for discontinuous Galerkin methods for ...
10/16/2019

### Convergence in the incompressible limit of new discontinuous Galerkin methods with general quadrilateral and hexahedral elements

Standard low-order finite elements, which perform well for problems invo...
12/11/2020

### A Split-Form, Stable CG/DG-SEM for Wave Propagation Modeled by Linear Hyperbolic Systems

We present a hybrid continuous and discontinuous Galerkin spectral eleme...
03/14/2020

### Parameter-robust Stochastic Galerkin mixed approximation for linear poroelasticity with uncertain inputs

Linear poroelasticity models have a number of important applications in ...
11/23/2020

### Stability of Discontinuous Galerkin Spectral Element Schemes for Wave Propagation when the Coefficient Matrices have Jumps

We use the behavior of the L_2 norm of the solutions of linear hyperboli...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

Consider an axi-parallel, open and bounded polyhedron , with Lipschitz boundary , 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 , and a pressure function  such that

 (1) −Δu+∇p =f in Ω, (2) ∇⋅u+(1−2ν)p =0 in Ω, (3) u =0 on ∂Ω.

Here, is the divergence operator, is the Poisson ratio, and is an external force (scaled by , where is Young’s modulus). We shall include the limit case  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 [21, 20, 22] on the numerical approximation of the Poisson equation in 3d polyhedra the use of -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

[20, 22] that exponential convergence rates for -DG methods can be achieved.

In our previous work [24], we have employed the approach [21, 20, 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 -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 ; 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  tends to the critical limit of  of incompressible materials. Furthermore, following the techniques presented in [20, 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 . 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 , , let signify the Lebesgue space of square-integrable functions equipped with the usual norm . Furthermore, we write  to denote the subspace of  of all functions with vanishing mean value on . The standard Sobolev space of functions with integer regularity exponent is signified by ; we write and for the corresponding norms and semi-norms, respectively. As usual, we define as the subspace of functions in with zero trace on

. For vector- and tensor-valued functions we use the standard notation

, and and , respectively. Moreover, for vectors , let be the matrix whose th component is .

## 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 such that

 (4) A(u,v)+B(v,p)=∫Ωf⋅vdx,−B(u,q)+C(p,q)=0,

for all , where

 A(u,v) :=∫Ω∇u:∇vdx, B(v,q) :=−∫Ωq∇⋅vdx, C(p,q) :=(1−2ν)∫Ωpqdx.

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

 (5) a(u,p;v,q)=∫Ωf⋅vdx∀(v,q)∈H10(Ω)3×L20(Ω),

with

 a(u,p;v,q):=A(u,v)+B(v,p)−B(u,q)+C(p,q).

It is straightforward to verify that  is a bounded bilinear form on , and that, for , it is also coercive. In particular, by application of the Lax-Milgram theorem, we conclude that the solution  of (5), and, thus, of (4), exists and is unique. In addition, for , which corresponds to the Stokes equations, problem (5) is still well-posed. Indeed, this is an immediate consequence of the inf-sup condition

 inf0≢q∈L20(Ω)sup0≢v∈H10(Ω)3−∫Ωq∇⋅vdx∥∇v∥L2(Ω)∥q∥L2(Ω)≥κ>0,

where is the inf-sup constant, depending only on ; 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 [11, 9, 10]. To this end, we denote by  the set of corners, and by  the set of edges of . Potential singularities of the solution are located on the skeleton of  given by

 S=(⋃c∈Cc)∪(⋃e∈Ee)⊂∂Ω.

Given a corner , an edge , and , we define the distance functions

 rc(x)=dist(x,c),re(x)=dist(x,e),ρce(x)=re(x)/rc(x).

Furthermore, for each corner , we signify by the set of all edges of  which meet at . For any , the set of corners of  is given by . Then, for , , respectively , and a sufficiently small parameter , we define the neighbourhoods

 ωc={x∈Ω:rc(x)<ε∧ρce(x)>ε∀e∈Ec},ωe={x∈Ω:re(x)<ε∧rc(x)>ε∀c∈Ce},ωce={x∈Ω:rc(x)<ε∧ρce(x)<ε}.

Moreover, we define the interior part of  by .

Near corners and edges , we shall use local coordinate systems in  and , which are chosen such that corresponds to the direction . Then, we denote quantities that are transversal to  by , and quantities parallel to  by . For instance, if  is a multi-index associated with the three local coordinate directions in  or , then we write , where and . In addition, we use the notation , and .

Following [6, Def. 6.2 and Eq. (6.9)], we introduce anisotropically weighted Sobolev spaces. To this end, to each and , we associate a corner and an edge exponent , respectively. We collect these quantities in the weight vector . Then, for , we define the weighted semi-norm

 |v|2Mmβ(Ω)=|v|2Hm(Ω0)+∑e∈E∑α∈N30|α|=m∥∥rβe+|α⊥|eDαv∥∥2L2(ωe)+∑c∈C∑α∈N30|α|=m∥∥rβc+|α|cDαv∥∥2L2(ωc)+∑c∈C∑e∈Ec∥∥rβc+|α|cρβe+|α⊥|ceDαv∥∥2L2(ωce),

as well as the full norm ; here, the operator denotes the partial derivative in the local coordinate directions corresponding to the multi-index . The space is the weighted Sobolev space obtained as the closure of  with respect to the norm .

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 such that, if the weight vector  satisfies

 0<βe<βE,0<βc<βC,e∈E, c∈C,

then, for every , the solution  of (1)–(3) with  fulfills .

## 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  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 -version finite element methods. Specifically, we begin from a coarse regular and shape-regular, quasi-uniform partition of  into  convex axi-parallel hexahedra. Each of these elements  is the image under an affine mapping  of the reference patch , i.e. . The mappings  are compositions of (isotropic) dilations and translations.

Based on the coarse partition (macro mesh)  we will use three canonical geometric refinements (patches) towards corners, edges and corner-edges of ; see Figure 1. They feature a refinement ratio , as well as a number of refinement levels ; to give an example, in Figure 1, we have selected , and .

Given a (fixed) refinement ratio  as well as a refinement level value , geometric meshes in  are now built by applying the patch mappings  to transform the above canonical geometric mesh patches on the reference patch  to the macro-elements , thereby yielding a local patch mesh  on . The patches  away from the singular support  (i.e. with ) are left unrefined, i.e. in this case we let . 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 . Then, a -geometric mesh in  is given by . Furthermore, the sequence  is referred to as a -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 [21].

### 3.2. Faces and face operators

We denote the set of all interior faces in  by , and the set of all boundary faces by . Further, let  signify the set of all (smallest) faces of . In addition, for an element , we denote the set of its faces by .

Next, we recall the standard DG trace operators. For this purpose, consider an interior face shared by two neighbouring elements . Furthermore, let , and be scalar-, vector, and tensor-valued functions, respectively, all sufficiently smooth inside the elements . Then, we define the following trace operators along :

 [[u]] =u|K♯nK♯+u|K♭nK♭, {{u}} =\nicefrac12(u|K♯+u|K♭), [[v]] =v|K♯⋅nK♯+v|K♭⋅nK♭, {{v}} =\nicefrac12(v|K♯+v|K♭), [[w––]] =w––|K♯⊗nK♯+–w|K♭⊗nK♭, {{w––}} =\nicefrac12(w––|K♯+w––|K♭).

Here, for an element , we denote by  the outward unit normal vector on . Similarly, for a boundary face , with , and a sufficiently smooth scalar function , we let , and , where is the outward unit normal vector on ; obvious modifications are made for vector- and tensor-valued functions in accordance with the definition above.

Finally, and denote the element-wise gradient and divergence operators, respectively. Here and in the sequel, we use abbreviations like

 ∫F(⋅)ds:=∑f∈F∫f(⋅)ds,∥∇h(⋅)∥2L2(Ω):=∑K∈Tσ,l∥∇(⋅)∥2L2(K).

### 3.3. Spectral DG discretizations

Given a geometric edge mesh on and a polynomial degree (which is assumed uniform and isotropic on ), we approximate (1)–(3) by finite element functions , where

 (6) Vh:={v∈L2(Ω)3:v|K∈Qk(K)3, K∈Tσ,l},Qh:={q∈L20(Ω):q|K∈Qk−1(K), K∈Tσ,l}.

Here, for , , denotes the space of all polynomials of degree at most in each variable on . In addition, we let

 V(h)=Vh+H10(Ω)3.

On the space  we consider the stabilization function given by

 (7) c(x):=γh−1(x)k2,

where is a penalty parameter independent of the refinement ratio , the number of refinement levels , and the polynomial degree . Furthermore, for , with , the mesh function  is defined by

 h(x):=⎧⎨⎩min{h⊥K♯,f,h⊥K♭,f},x∈f⊂FI,f=∂K♯∩∂K♭,with K♯,K♭∈Tσ,l,h⊥K,f,x∈f⊂FB,f=∂K∩∂Ω,with K∈Tσ,l.

In this definition, for  and , we denote by the diameter of the element  in the direction perpendicular to the face .

Then, we consider the following mixed discontinuous Galerkin discretization of (4): Find such that

 (8) Ah(uh,v)+Bh(v,ph)=∫Ωf⋅vdx,−Bh(uh,q)+Ch(ph,q)=0,

for all . The forms , , and are given, respectively, by

 Ah(u,v) :=∫Ω∇hu:∇hvdx−∫F(θ{{∇hv}}:[[u]]+{{∇hu}}:[[v]])ds +∫Fc[[u]]:[[v]]ds, (9) Bh(v,q) :=−∫Ωq∇h⋅vdx+∫F{{q}}[[v]]ds, Ch(p,q) :=(1−2ν)∫Ωpqdx,

where is a fixed parameter. Different choices of  refer to various types of interior penalty DG methods: for instance, the form may be chosen to correspond to the symmetric (for ), incomplete (for ), or non-symmetric (for ) 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 such that

 (10) ah(uh,ph;v,q)=∫Ωf⋅vdx

for all , where

 (11) ah(u,p;v,q):=Ah(u,v)+Bh(v,p)−Bh(u,q)+Ch(p,q).

### 3.4. Discrete inf-sup stability and well-posedness

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

 (12) |||(v,q)|||2DG:=∥v∥2h+(2−2ν)∥q∥2L2(Ω),

for any , where

 ∥v∥2h:=∥∇hv∥2L2(Ω)+∫Fc|[[v]]|2ds,v∈V(h).

If the Poisson ratio satisfies , and provided that the penalty parameter  featured in (7

) is chosen sufficiently large, then the following coercivity estimate can be shown:

 ah(u,p;u,p)≥C∥u∥2h+(1−2ν)∥p∥2L2(Ω)≥C(1−2ν)|||(u,p)|||2DG∀(u,p)∈Vh×Qh;

here, is a constant independent of , , , and the aspect ratio of the anisotropic elements.

This result can be made stronger if a discrete inf-sup condition on the form on geometric edge meshes is assumed: Let be a geometric edge mesh on as defined in Section 3.1, with refinement ratio , and layers of refinement. Suppose that there exist constants and  that may depend on , , and on the macro-element mesh , but are independent of , , and the aspect ratio of the anisotropic elements in , such that there holds

 (13) γB:=inf0≢q∈Qhsup0≢v∈VhBh(v,q)∥v∥h∥q∥L2(Ω)≥ϰk−ρ,

as .

###### Remark 3.1.

In [19] it was proved that this assumption is fulfilled with (and any ). Our numerical computations in Section 4.4.1 below indicate, however, that the dependence of the right-hand side of (13) on is much weaker than .

The following result, which implies the well-posedness of (8) even in the incompressible limit , follows immediately from [24, Theorem 5.1].

###### Theorem 3.2.

Let . If (13) holds true, then we have the inf-sup condition

 (14) γa:=inf(u,p)∈Vh×Qh(u,p)≠(0,0)sup(v,q)∈Vh×Qh(v,q)≠(0,0)ah(u,p;v,q)|||(u,p)|||DG|||(v,q)|||DG≥Cmax{k−2ρ,1−2ν},

with a constant that depends on the penalty parameter , however, is independent of , , , and the aspect ratio of the anisotropic elements.

We emphasize that, for fixed , the stability bound (14) does not deteriorate as .

## 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  from (6). Due to the global zero mean constraint contained in , the construction of  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

 (15) ˜Qh:={q∈L2(Ω):q|K∈Qk−1(K), K∈Tσ,l},

and to impose the zero mean condition by means of a Lagrange multiplier technique. Noticing that , 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 , on the new space , where  is the full space from (15). To this end, we introduce an auxiliary variable which takes the role of the mean value of the pressure  on . More precisely, let us consider the following augmented DG formulation: Find such that

 (16) Ah(˜uh,˜v)+Bh(˜v,˜ph)=∫Ωf⋅˜vdx,−Bh(˜uh,˜q)+Ch(˜ph,˜q)−˜r∫Ω˜qdx=0,˜s\fintΩ˜phdx−˜r˜s=0,

for all . Here we use the notation

 \fintΩ(⋅)dx:=1|Ω|∫Ω(⋅)dx

to denote the mean value integral on . Note also that this new system may be written in a more compact way: Find such that

 ˜ah(˜uh,˜ph,˜r;˜v,˜q,˜s)=∫Ωf⋅˜vdx∀(˜v,˜q,˜s)∈Vh×˜Qh×R,

where

 ˜ah(˜uh,˜ph,˜r;˜v,˜q,˜s) :=Ah(˜uh,˜v)+Bh(˜v,˜ph)−Bh(˜uh,˜q)+Ch(˜ph,˜q) −˜r∫Ω˜qdx+˜s\fintΩ˜phdx−˜r˜s.

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 the space of all (globally) constant functions on , and we note that

 (17) ˜Qh=Qh⊕Q0.
###### Lemma 4.1.

The form from (9) satisfies

 (18) Bh(v,q)=0∀(v,q)∈Vh×Q0(Ω).

Conversely, if, for given , there holds that  for any , then it follows that .

###### Proof.

Given . Since  is one-dimensional we may, without loss of generality, suppose that . Then, we have

 Bh(v,q)=−∑K∈Tσ,l∫K∇h⋅vdx+∫F[[v]]ds.

By applying the Gauss-Green theorem on each element , we obtain two expressions which are identical, and, thus, cancel out:

 Bh(v,q)=−∑K∈Tσ,l∫∂Kv⋅nKds+∫F[[v]]ds=0.

Let now , with , and  for all . Then, the inf-sup condition (13) implies that

 0<γB≤sup0≠v∈VhBh(v,q−\fintΩqdx)∥v∥h∥∥q−\fintΩqdx∥∥L2(Ω)=sup0≠v∈Vh−Bh(v,\fintΩqdx)∥v∥h∥∥q−\fintΩqdx∥∥L2(Ω).

Applying (18) yields a contradiction, and, consequently, we deduce that . Thus, we have . 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 , and is the solution of the original DG formulation (8) with .

###### Proof.

We proceed in three steps.

#### Step 1:

We first show that the new formulation enforces the pressure to have zero mean. To this end, we choose the test variable to be . Thus, from the third equation of (16), we deduce that

 \fintΩ˜phdx=˜r.

Furthermore, let us choose . From the second equation of (16), and with the aid of (18), we infer that

 0=(1−2ν)∫Ω˜phdx−(\fintΩ˜phdx)(∫Ω1dx)=−2ν∫Ω˜phdx,

i.e. , since .

#### Step 2:

Next, we show that , where is the solution from (8), solves (16). The first and last equation in (16) are clearly fulfilled with . The second equation in (16) simplifies to

 −Bh(uh,˜q)+Ch(ph,˜q)=0,∀˜q∈˜Qh.

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

 −Bh(uh,˜q)+Ch(ph,˜q) =−Bh(uh,\fintΩ˜