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
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 
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 , we have employed the approach [21, 20, 22] in order to apply the high-order mixed DG methods introduced in  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 . A focal point of the article  has been to provide an inf-sup stability analysis for mixed IPDG schemes on anisotropic meshes. Our results, which are based in parts on , are explicit with respect to both the (uniform) polynomial degree and the Poisson ratio ; cf. also  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 .
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.
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, and and , respectively. Moreover, for vectors , let be the matrix whose th component is .
2. Linear elasticity in polyhedra
2.1. Mixed formulation and well-posedness
for all , where
More compactly, we can write (4) equivalently in the form
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
Following  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
Given a corner , an edge , and , we define the distance functions
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
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
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 .
3. Mixed discontinuous Galerkin methods on geometric meshes
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 , where such meshes have been studied in the context of DGFEM for the Stokes equations; see also the earlier paper  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 .
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 :
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
3.3. Spectral DG discretizations
Here, for , , denotes the space of all polynomials of degree at most in each variable on . In addition, we let
On the space we consider the stabilization function given by
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
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
for all . The forms , , and are given, respectively, by
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.
3.4. Discrete inf-sup stability and well-posedness
for any , where
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:
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
Let . If (13) holds true, then we have the inf-sup condition
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
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
for all . Here we use the notation
to denote the mean value integral on . Note also that this new system may be written in a more compact way: Find such that
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.
The form from (9) satisfies
Conversely, if, for given , there holds that for any , then it follows that .
Given . Since is one-dimensional we may, without loss of generality, suppose that . Then, we have
By applying the Gauss-Green theorem on each element , we obtain two expressions which are identical, and, thus, cancel out:
Now we can state the equivalence of the two formulations.
We proceed in three steps.
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
i.e. , since .