A robust discontinuous Galerkin scheme on anisotropic meshes

Discontinuous Galerkin (DG) methods are extensions of the usual Galerkin finite element methods. Although there are vast amount of studies on DG methods, most of them have assumed shape-regularity conditions on meshes for both theoretical error analysis and practical computations. In this paper, we present a new symmetric interior penalty DG scheme with a modified penalty term. We show that, without imposing the shape-regularity condition on the meshes, the new DG scheme inherits all of the good properties of standard DG methods, and is thus robust on anisotropic meshes. Numerical experiments confirm the theoretical error estimates obtained.

The above mentioned papers and textbooks confirm that a shape-regularity condition has always been imposed on meshes for theoretical error analysis of DG methods. The shape-regularity condition requires that elements in the meshes must be neither very "flat" nor "degenerated" (see Definition 1). If a mesh contains very flat elements, it is said to be anisotropic. The simple numerical experiment described in Section 3.3 shows that the standard SIP-DG method is not robust on anisotropic meshes, and that the shape-regularity condition is crucial for practical computations.
The purpose of this paper is to introduce a new SIP-DG scheme that is robust on anisotropic meshes. To this end, we use the general trace inequality to define a new penalty term for the proposed SIP-DG scheme in Section 4. In Section 5, we show that the new scheme inherits all of the good properties of standard SIP-DG methods. That is, if the penalty parameter is sufficiently large, the new SIP-DG scheme is consistent, coercive, stable, and bounded on arbitrary (proper) meshes (Lemma 3 new and Lemma 8). Those properties immediately yield error estimations of the new SIP-DG scheme without imposing the shape-regularity condition (see Corollary 11). An immediate consequence is the error estimate of order O(h k ) under the maximum angle condition on meshes (see Corollary 12). In Section 6, we present the results of numerical experiments to confirm the theoretical results obtained. From the results, we conclude that the newly presented SIP-DG scheme is robust on anisotropic meshes.

The model problem
Let Ω ⊂ R d , (d = 2, 3) be a bounded polyhedral domain. Let L 2 (Ω), H 1 (Ω), and H 1 0 (Ω) be the usual Lebesgue and Sobolev spaces (see Section 2.3 for notation of their (semi)norms and inner products). We consider the following Poisson problem: find u ∈ H 1 (Ω) such that −∆u = φ in Ω, u = 0 on ∂Ω, (2.1) where φ ∈ L 2 (Ω) is a given function. The weak form of the model problem is as follows: a(u, v) := Ω ∇u · ∇v dx = (φ, v) Ω , ∀v ∈ H 1 0 (Ω). (2.2) The model problem (2.1) is said to satisfy elliptic regularity if there exists a positive constant C ell such that the following a priori estimate holds for the exact solution u: It is well known that if Ω is convex, then the model problem (2.1) satisfies elliptic regularity [8].

Meshes of Ω
Although DG methods allow elements with a variety of geometries, we consider only simplicial elements in this paper. Thus, we suppose that the domain Ω is divided into a finite number of triangles (d = 2) or tetrahedrons (d = 3), which are assumed to be closed sets. A mesh (or triangulation) of Ω is denoted by T h . That is, T h is a finite set of triangles or tetrahedrons that has the following properties: where int T i is the interior of T i . For T ∈ T h , let n T be the unit outer normal vector on ∂T .
In this paper, we assume that meshes are proper (or face-to-face). This means that, for T 1 , For a d-simplex T and its facet f , their Lebesgue and Hausdorff measures are denoted by |T | and | f |, respectively.
For T , we define h T := diam T . If d = 2, R T is the circumradius of T . Note that where l 1 ≤ l 2 ≤ · · · ≤ h T are the lengths of the edges of T . As has been seen in [10,11,12,13,14,15], R T is an important parameter in measuring interpolation errors on d-simplices. For example, the errors of Lagrange interpolation on T are bounded in terms of R T , as presented by (5.10) and (5.11).
Remark. The best definition of R T for tetrahedrons remains an open problem. A simple example immediately rejects the idea that R T for a tetrahedron might be the radius of its circumsphere [14, p. 3]. The definition (2.5) is given in [10]. In [10], another parameter, denoted by H T , is introduced, and it is shown that R T and H T are equivalent (see also [15]). In [14], the projected circumradius of T is defined for a tetrahedrons. It is conjectured that R T and the projected circumradius are equivalent.
Suppose that we consider a (possibly infinite) family of meshes {T h } h>0 with h → 0.
Definition 1 (1) The family of meshes {T h } h>0 is said to satisfy the shape-regularity condition with respect to σ if there exists a positive constant σ such that where ρ T is the diameter of the inscribed ball of T . We call σ the shape-regular constant in this paper.
(2) The family of meshes {T h } h>0 is said to satisfy the maximum angle condition with respect to a constant C max (π/3 < C max < π) if the following hold for all T ∈ T h and for all T h : -An arbitrary inner angle θ of T is θ ≤ C max (d = 2), or -An arbitrary inner angle θ of any facet of T is θ ≤ C max , and an arbitrary dihedral angle η of T is η ≤ C max (d = 3).
(3) The family of meshes {T h } h>0 is said to satisfy the circumradius condition if the family satisfies In this paper, we always assume that the family {T h } h>0 of meshes satisfies the circumradius condition.
Remark. If we deal with only a finite family of meshes {T h } and we take a sufficiently large σ > 0, then the family satisfies the shape-regularity condition because {T h } contains only a finite number of d-simplices. However, if σ is too large (say, σ ≥ 10), we commonly say that {T h } is not shape-regular. In such a case, as mentioned in Section 1, {T h } is said to be anisotropic.

Function spaces
Let k ≥ 1 be a positive integer. We use the notation L 2 (Ω), L 2 ( f ) ( f ∈ F h ), H k (Ω), H 1 0 (Ω) for the usual Lebesgue and Sobolev spaces. We denote their norms and seminorms by, for example, · 0, f , | · | 1,Ω , and their inner products by (·, ·) f , (·, ·) Ω . Let P k (K) be the set of all polynomials defined on the closed set K with degree less than or equal to a positive integer k.

Jump and mean of functions on f ∈ F h
For each interior facet f ∈ F o h , there are two simplices that share f . We number those simplices as T i f ∈ T h (i = 1, 2) and fix the numbering once the mesh is obtained.
Recalling that n T is the unit outer nomal vector on ∂T , define n f := n T 1 f . For v ∈ H 2 (T h ), we set v 1 := v| T 1 f and v 2 := v| T 2 f . We denote the trace operator on T i f to f by γ i f (i = 1, 2). Define The jump [∇v] f and average {∇v} f are defined in a similar way. If g ∈ F ∂ h , then g ⊂ ∂Ω. Let g ⊂ ∂T g with T g ∈ T h . Then, define [v] = {v} := γ g (v| T g ).
3 Standard SIP-DG scheme

Definition of SIP-DG scheme
In the SIP-DG scheme, the bilinear form a(u, v) in (2.2) is discretized as Here, η is a penalty parameter that is taken to be sufficiently large. To make the notation concise, we set The terms J h (v, w h ) and P std h (v, w h ) are called the jump term and penalty term, respectively. The discretized bilinear form a std h (v, w h ) is written as

Definition 2
The SIP-DG scheme for the model problem is defined as follows: find

Properties of SIP-DG scheme and error analysis
In the following, we summarize the properties of the SIP-DG method. For their proofs, readers are referred to the standard textbooks [4], [6], [7], [17]. In this section, we mainly refer to [6].
Lemma 3 (Consistency) [6, Lemma 4.8] The exact solution u ∈ V * of the model problem (2.1) is consistent: Therefore, the solution u h ∈ V h of the SIP-DG method (3.2) satisfies the Galerkin orthogonality: We define the norms associated with the bilinear form a std h as:

Lemma 4
Suppose that the mesh T h is shape-regular with respect to a constant σ > 0 and the penalty parameter η is sufficiently large. Then, (2) (Discrete stability) The following inequality holds: (3) (Boundedness) [6, Lemma 4.16] The following inequality holds: where the constant C := C(η, σ) is independent of h.
Theorem 5 [6, Theorem 4.17] Suppose that the mesh T h is shape-regular with respect to a constant σ > 0 and the penalty parameter η is sufficiently large. Then, there exists a unique SIP-DG solution u h ∈ V h of (3.2), and the following error estimate holds: where the constant C depends only on the penalty parameter η and σ.
where the constant C depends on η and σ, but is independent of h.

Numerical experiments (part 1)
We consider a numerical experiment to examine how the shape-regular constant σ affects the practical computations involved in the standard SIP-DG scheme.
Set Ω := (0, 1) × (0, 1) and φ(x, y) := π 2 sin(πx) sin(πy) in the model problem (2.1). Then, the exact solution is u(x, y) = sin(πx) sin(πy)/2. Let n and m be positive integers. We divide the horizontal and vertical sides of Ω into n and m equal segments, respectively. We then draw diagonal lines in each small rectangle to define the mesh, as depicted in Figure 1. We fix n = 40 and the penalty parameter η = 10. We apply the standard SIP-DG method to the model problem with various m. The conjugate gradient method with the incomplete Cholesky decomposition preconditioner (ICCG) is used for the linear solver. The successive over-relaxation (SOR) method is also used occasionally to check whether the obtained u h is reasonable. The results are summarized in Table 1. Table 1 Errors produced by the SIP-DG method. 1/2 , and the "DG-error" is u − u h DG . We employ the 4-point Gauss quadrature of degree 3 on triangles to compute those errors. We see that the errors given by the SIP-DG method decrease as m increases until m = 100, which is consistent with the theoretical error estimates. However, the errors increase as m increases from m = 120. The ICCG iterations do not converge for m = 200, while the SOR iterations give almost the same results until m = 160. For m = 180, 200, the SOR iterations converge quickly but the obtained u h are not reasonable.
We also examined the case m = 400. In this case, we required η = 30 to obtain reasonable u h .

New penalty term and SIP-DG scheme
From the numerical experiments described in the previous section, we can conclude that the shape-regularity condition is crucial for the standard SIP-DG method with a fixed penalty parameter η. It is natural to wonder why this is the case.
The most important term in the SIP-DG scheme is the penalty term P std h (v, w h ) defined by (3.1). The penalty term originates from the trace inequalities where f is an arbitrary facet of T ∈ T h . Note that the constants C tr i (i = 1, 2) strongly depend on the shape-regular constant σ.
To avoid imposing the shape-regularity condition, we adopt the general trace inequalities v 0, f ≤ C tr which are valid on an arbitrary d-simplex T . Warburton and Hesthaven [18] presented explicit forms of the constants C tr i (i = 3, 4) that are independent of the geometry of T . Note that (4.1) is a "simplified" version of (4.2) under the shape-regularity condition. We introduce a quantity on each f ∈ F h below. Let Let T i f be the d-simplex whose vertices are those of f and the barycenter of T i Let g ⊂ ∂T g with T g ∈ T h and T g be the d-simplex whose vertices are those of g and the barycenter of T g . Define Note that See Figure 2. We remark that the only information we need for the simplex T i f is its g ⊂ ∂Ω With the quantity | f | | T f | defined for f ∈ F h , we (re)define the penalty term and SIP-DG bilinear form a std h as Note that, if a mesh satisfies the shape-regularity condition, | f | | T f | ≈ 1 h f for any facet f ∈ F h , and P new h becomes equivalent to the standard penalty term P std h .

Definition 7
The SIP-DG scheme for the model problem is redefined as follows: find Remark. To see the meaning of the new penalty terms (4.4), let us consider a small part of the mesh of Figure 1, as shown in Figure 3. The mesh consists of congruent right triangles. Let T be one of the right triangles, and f 1 , f 2 be its edges as depicted in Figure 3. Set h = | f 1 |. Then, | f 2 | = αh, where 0 < α < 1. Elementary geometry tells us that f 2 | T f 2 | = 12/h, Note that if T is becoming degenerated, the coefficient 12η/α will increase. This means that we can interpret the coefficient η | f | | T f | as an "improved" version of the standard penalty term coefficient η h f with the "adaptive" parameter α.
We redefine the norms associated with the SIP-DG scheme: The following lemma holds. Its proof is quite similar to the standard proofs. Here, we loosely follow the proofs given by Di Pietro and Ern [6].
(2) (Discrete stability) The following inequality holds: (3) (Boundedness) The following inequalities hold: where the constant C := C(η, C tr 4 ) is independent of h and the geometry of elements in T h .
Hence, we have We also obtain a similar inequality for the case Thus, it follows from the arithmetic-geometric mean that for some constant δ > 0. Set δ := 1/2, and let η be sufficiently large so that η ≥ (C tr 4 ) 2 . Then, the following coercivity holds: (2) For arbitrary v h ∈ V h , we have because of the coercivity.
(3) By the Cauchy-Schwarz inequality, we see that

It follows from (5.5) and (4.3) that
Therefore, we obtain Gathering these inequalities together, we conclude that (5.3) and (5.4) hold.
Theorem 9 For the exact solution u ∈ V * h of the model problem (2.1) and its SIP-DG solution u h ∈ V h of (4.6), the following error estimate holds: where the constant C = C(η, C tr 4 ) is independent of h and the geometry of elements in T h .
Proof By the consistency of a new h , we have Thus, the discrete coercivity and boundedness yield for an arbitrary y h ∈ V h . Therefore, we obtain Taking the infimum for y h ∈ V h and rewriting C, we conclude that (5.7) holds.
To derive a more practical error estimation, we prepare another general trace inequality. Let T be a d-simplex and f be a facet of T . The following inequality holds [17, p.24]: ∇v · n 0, f ≤ C tr 5 | f | 1/2 |T | 1/2 |v| 1,T + h T |v| 2,T , ∀v ∈ H 2 (T ). (5.8) Now, let I k h u ∈ P k (T h ) be an interpolation of u that satisfies Note that the usual Lagrange interpolation satisfies (5.9). Insert I k h u into y h in (5.7), and set U := u − I k h u. Then, [U] f = 0 on any f ∈ F h , and P new h (U, U) = 0. Therefore, we see that , and thus That is, the following theorem has been proved.
Theorem 10 Let u ∈ V * h be the exact solution of the model problem (2.1), and u h ∈ V h be its SIP-DG solution (4.6). Then, we have the following error estimation: where I k h u is an interpolation that satisfies (5.9), and the constant C = C(C tr 4 , C tr 5 , η) is independent of h and the geometry of elements in T h .
Let I 1 h u be the usual Lagrange interpolation of u and let R T be the quantity defined in Section 2.2. Suppose that k = 1, d = 2 and u ∈ H 2 (Ω). Then, |u − I 1 h u| 2,T = |u| 2,T , and, from the results in [11,13], |u − I 1 h u| 1,T ≤ C L1 R T |u| 2,T , ∀u ∈ H 2 (T ), (5.10) where the constant C L1 is independent of the geometry of T . Thus, we find that where R := max T ∈T h R T and h := max T ∈T h h T . Suppose that k ≥ 2, d = 2, 3, and u ∈ H k+1 (Ω). Then, from [13,14,10], and where the constant C Lk is independent of the geometry of T . Therefore, we have obtained the following corollary.
Corollary 11 Let u ∈ V * h be the exact solution of the model problem (2.1), and u h ∈ V h be its SIP-DG solution (4.6). Suppose that u ∈ H k+1 (Ω). Then, we have the following error estimations: where the constant C = C(C tr 4 , C tr 5 , C Lk , η) is independent of h, R, and the geometry of elements in T h .
Let d = 2. In this case, elementary geometry (the law of sines) tells us that a mesh T h satisfies the maximum angle condition if and only if there exists a constant σ 2 such that Recently, it is reported that the same situation holds for tetrahedrons. That is, a tetrahedron satisfies the maximum angle condition if and only if (5.12) holds [9]. See also [15]. Hence, we have the following corollary.

Corollary 12
Suppose that k ≥ 1 if d = 2, and k ≥ 2 if d = 3. Let u ∈ V * h be the exact solution of the model problem (2.1), and u h ∈ V h be its SIP-DG solution (4.6). Suppose that u ∈ H k+1 (Ω) and T h satisfies the maximum angle condition. Then, we have the following error estimation: where the constant C depends only on C tr 4 , C tr 5 , C Lk , η, and σ 2 .
For an L 2 error estimate, we have the following theorem, which is quite similar to the Aubin-Nitsche lemma.
Theorem 13 Let d = 2. Let u ∈ V * h be the exact solution of the model problem (2.1), and u h ∈ V h be its SIP-DG solution (4.6). Suppose that the model problem satisfies the elliptic regularity (2.3). Then, we have the following error estimation: where C = C(C tr 4 , C tr 5 , η, C ell , C L1 ).
Proof. Setting φ := u − u h in (2.2), we consider an auxiliary problem: find ζ ∈ H 1 0 (Ω) such that Because of the elliptic regularity (2.3), we have ζ 2,Ω ≤ C ell u − u h 0,Ω . Because ζ ∈ H 2 (Ω), we have [∇ζ] · n f = 0 on any f ∈ F o h , and [ζ] = 0 on any f ∈ F h . Hence, the definition of a new h implies that Let I 1 h ζ be the piecewise linear Lagrange interpolation of ζ on T h . Because of the Galerkin orthogonality in Lemma 3 new , we infer that Therefore, it follows from (5.4) and (5.10) that This completes the proof.
Corollary 14 Let d = 2. Let u ∈ V * h be the exact solution of the model problem (2.1), and u h ∈ V h be its SIP-DG solution (4.6). Suppose that elliptic regularity holds, and u ∈ H k+1 (Ω). Then, we have the following error estimation: where the constant C = C(C tr 4 , C tr 5 , η, C ell , C L1 ) is independent of h, R, and the geometry of elements in T h .
Remark. It is conjectured that Theorem 13 and Corollary 14 hold for d = 3. To show this conjecture, we need to use a different interpolation, such as the Crouzeix-Raviart interpolation, and analyze its error in terms of the · new DG * norm.

Numerical experiments for the new SIP-DG scheme
In this section, we report the results of numerical experiments to confirm the theoretical results obtained in the previous section. First, we consider the same numerical experiment as in Section 3.3 (the same domain Ω, same function φ, and same meshes T h ) with the new SIP-DG scheme. We fix n = 40 and the penalty parameter η = 0.8. 1 Then, the new SIP-DG method is applied to the model problem (2.1) with various m.
The results are summarized in Table 2. Note that the discretized solution u h is stable for all cases, and the errors decrease as m increases, which is consistent with the error estimations.
For a given positive integer N and α > 1, we consider the isosceles triangle with base length h := 2/N and height 2/ 2/h α ≈ h α (then, R T ≈ h α /2 + h 2−α /8), as depicted in Figure 4. We triangulate Ω with this triangle as shown in Figure 5. Note that if h → 0, these meshes satisfy neither the shape-regularity condition nor the maximum angle condition. Furthermore, if α ≥ 2.0 and h → 0, the meshes do not satisfy the circumradius condition. See Kobayashi-Tsuchiya [12].  We compute the SIP-DG solutions on this mesh with α = 1.2, 1.5, 1.8, 2.0, 2.1 and various N, and measure the errors. Figures 6 and 7 list the numerical results. Tables 3 and 4 contain concrete values of the numerical results for α = 1.5 and α = 2.0, respectively.
It seems that both errors behave similarly as α varies. Observing Figure 6 and Figure 7, we can infer that , that is, the errors are governed by the parameter R, not h. For example, when α = 2.1, the errors a (0) h (u − u h , u − u h ) 1/2 and P (0) h (u − u h , u − u h ) 1/2 increase as h decreases. However, the ratio a (0) h (u − u h , u − u h ) 1/2 /R and P (0) h (u − u h , u − u h ) 1/2 /R seem to be bounded by a constant.  Fig. 6 Behavior of the error a (0) h (u − u h , u − u h ) 1/2 with respect to α. The horizontal axis represents the maximum diameter (left) and the circumradius of the triangles (right), and the vertical axis represents the error. The legend indicates the value of α in each case.

Conclusion
Using the general trace inequality (4.2), we have presented a new penalty term (4.4) and a new SIP-DG scheme (4.5). We have established its error estimates without imposing the shape-regularity condition. From numerical experiments, we have confirmed that the new SIP-DG scheme is robust on anisotropic meshes.
The main idea in this paper is to use | f | | T f | instead of the standard term 1 h f in the penalty term. Because this idea is very simple, we expect that it can be extended in  Fig. 7 Behavior of the error P new h (u − u h , u − u h ) 1/2 with respect to α. The horizontal axis represents the maximum diameter (left) and the circumradius of the triangles (right), and the vertical axis represents the error. The legend indicates the value of α in each case.  many directions and will be used in many DG schemes. In the following, we mention some directions that are immediately apparent.
-Application to many variants of DG methods.
-Application of the new SIP-DG scheme to non-proper meshes, and the establishment of error estimations for these cases. -Investigation of the possibility of applying the proposed approach to hybridized DG methods.