Anisotropic meshes and stabilized parameters for the stabilized finite element methods

# Anisotropic meshes and stabilized parameters for the stabilized finite element methods

Yana Di  Hehu Xie  Xiaobo Yin LSEC, NCMIS, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China (yndi@lsec.cc.ac.cn). The author is supported in part by National Natural Science Foundation of China (11271358) and the National Center for Mathematics and Interdisciplinary Science, CAS. LSEC, NCMIS, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China (hhxie@lsec.cc.ac.cn). The author is supported in part by National Natural Science Foundation of China (91330202, 11371026, 11001259, 2011CB309703) and the National Center for Mathematics and Interdisciplinary Science, CAS. Corresponding author, School of Mathematics and Statistics & Hubei Key Laboratory of Mathematical Sciences, Central China Normal University, Wuhan 430079, China, (yinxb@mail.ccnu.edu.cn). The author is supported in part by the National Natural Science Foundation of China (11201167).
###### Abstract

We propose a numerical strategy to generate the anisotropic meshes and select the appropriate stabilized parameters simultaneously for two dimensional convection-dominated convection-diffusion equations by stabilized continuous linear finite elements. Since the discretized error in a suitable norm can be bounded by the sum of interpolation error and its variants in different norms, we replace them by some terms which contain the Hessian matrix of the true solution, convective fields, and the geometric properties such as directed edges and the area of the triangle. Based on this observation, the shape, size and equidistribution requirements are used to derive the corresponding metric tensor and the stabilized parameters. It is easily found from our derivation that the optimal stabilized parameter is coupled with the optimal metric tensor on each element. Some numerical results are also provided to validate the stability and efficiency of the proposed numerical strategy.

Key words. Metric tensor, anisotropic, convection-diffusion equation, stabilized parameter, stabilized finite element methods

AMS subject classifications. 65N30, 65N50

## 1 Introduction

This paper is concerned with the numerical solution of the following scalar convection-diffusion equation

 {−εΔu+b⋅∇u=f,     in Ω,u=g,   on ∂Ω, (1.1)

where is a bounded polygonal domain with boundary , is the constant diffusivity, is the given convective field satisfying the incompressibility condition in , is the source function, and represents the Dirichlet boundary condition.

Despite the apparent simplicity of problem (1.1), its numerical solution become particularly challenging when convection dominates diffusion (i.e., when ). In such cases, the solution usually exhibits very thin layers across which the derivatives of the solution are large. The widths of these layers are usually significantly smaller than the mesh size and hence the layers can be hardly resolved. As a result of this, on meshes which do not resolve the layers, standard Galerkin finite element methods have poor stability and accuracy properties.

To enhance the stability and accuracy of the Galerkin discretization in convection-dominated regime, various stabilized strategies have been developed. Examples include upwind scheme [30], streamline diffusion finite element methods (SDFEM), also known as streamline-upwind Petrov-Galerkin formulation (SUPG) [14, 33], the Galerkin/least squares methods (GLS) [34], residual-free-bubbles (RFB) methods [10, 11, 17, 24], local projection schemes [5, 8, 38, 43, 44], exponential fitting [3, 12, 52], discontinuous Galerkin methods [4, 31], subgrid-scale techniques [27, 26], continuous interior penalty methods [9, 15, 16, 19, 21, 22], and spurious oscillations at layers diminishing (SOLD) methods (also known as shock capturing methods) [36, 37, 39], interested readers are referred to [49] for an extensive survey of the literature.

However, if uniform meshes are used for stabilized finite element methods, oscillations still exist near the layers in some cases although very fine meshes are used [35]. Hence, it is more appropriate to generate adaptively anisotropic meshes to capture the layers. There are some recent efforts directed at constructing adaptive anisotropic meshes which combine a stabilized scheme and some mesh modification strategies. For example, the resolution of boundary layers occurring in the singularly perturbed case is achieved using anisotropic mesh refinement in boundary layer regions [2], where the actual choice of the element diameters in the refinement zone and the determination of the numerical damping parameters is addressed. In [48] an adaptive meshing algorithm is designed by combining SUPG method, an adapted metric tensor and an anisotropic centroidal Voronoi tessellation algorithm, which is shown to be robust in detecting layers and efficient in avoiding non-physical oscillations in the numerical approximation. Sun et al. [50] develop a multilevel-homotopic-adaptive finite element method (MHAFEM) by combining SDFEM, anisotropic mesh adaptation, and the homotopy of the diffusion coefficient. The authors use numerical experiments to demonstrate that MHAFEM can efficiently capture boundary or interior layers and produce accurate solutions.

This list is by no means exhausted, and there are many efforts to construct adaptive anisotropic meshes by combining a stabilized scheme and some mesh modification strategies. However, so far as we know, there are still two key problems which haven’t been solved in rigorous ways.

First, although there are many results on optimal anisotropic meshes for minimizing the interpolation error and also the discretized error of finite element methods for solving the Laplace equation due to céa’s lemma, its extension to the discretized error of stabilized finite element methods for the convection-dominated convection-diffusion equation is not clear. So far as we know most results on optimal anisotropic meshes for the stabilized finite element methods applied to the convection-dominated convection-diffusion equation is just the same with that for the interpolation error. In fact, this strategy is not always optimal, which will be illustrated experimentally later in this paper.

Second, there is a crucial factor for most SFEMs, that is, the proper selection of the stabilization parameter on the element (there are some attempts avoiding explicitly mesh-dependent stabilization parameter, e.g. [6, 7, 25]). For example, the standard choice for quasi-uniform triangulations for SUPG is ([49, P. 305-306])

 αK={α0hKPeK>1,(% convection-dominated case)α1h2K/εPeK≤1,(diffusion-dominated case)

with appropriate positive constants and . Here is the mesh Peclét number and is the diameter of the element . A more sophisticated choice is to replace the diameter of the element in the above definition for by its streamline diameter which is the maximal length of any characteristic running through . How to extend the strategy for quasi-uniform triangulations to the case of anisotropic meshes? There are some attempts which basically use the analog of isotopic case to get the following form of stabilization parameters

 (1.2)

For example, Nguyen et al. [48] use stabilization parameters as the form of (1.2) by setting as the length of the longest edge of the element projected onto the convective field , this choice of together with (1.2) is denoted by “LEP” in this paper. Another similar choice of is the length of the projection of the longest edge of the element onto the convective field , which together with (1.2) is denoted by “PLE”. As pointed in isotropic case, a more sophisticated choice is to choose as the maximal length of any characteristic running through , i.e., the diameter of in the direction of the convection , which together with (1.2) is denoted by “DDC”. There are also some other choices of stabilized parameters derived by relatively rigorous theory on anisotropic meshes. For example, in [2] an anisotropic a priori error analysis is provided for the advection-diffusion-reaction problem. It is shown that the diameter of each element (together with (1.2), it is denoted by “DEE”), should be used as in (1.2) for the design of the stability parameters in the case of external boundary layers. In [40] an alternative approach is proposed showing that the diameter of each element is again the correct choice. Micheletti et al. [45] consider the GLS methods for the scalar advection-diffusion and the Stokes problems with approximations based on continuous piecewise linear finite elements on anisotropic meshes, where new definitions of the stability parameters are proposed. Cangiani and Süli [17] use the stabilizing term derived from the RFB method to redefine the mesh Péclet number and propose a new choice of the streamline-diffusion parameters (This well-known fact that the RFB method and the SDFEM are equivalent under certain conditions was first observed by Brezzi and Russo [13]. The similar idea was also used in [41]) that is suitable for use on anisotropic partitions. Although there are so many strategies on selection of the stabilization parameters, it is still hard to show which is optimal. Besides, there is little result on the relationship between the strategy to generate the anisotropic meshes and the selection of stabilized parameters.

In this paper, we propose a strategy to generate the anisotropic meshes and select the appropriate stabilized parameters simultaneously of stabilized continuous linear finite elements for two dimensional convection-dominated problems. As in [23, 34, 45], the discretized error (the difference between the true solution and the stabilized finite element solution) in a suitable norm can be bounded by the sum of interpolation error and its variants in different norms. Based on this result, we use the idea in our recent work [51] to replace these norms of interpolation error by some terms which contain the Hessian matrix of the true solution, convective fields , and the geometric properties such as directed edges and the area of the triangle. After that, we use the shape, size and equidistribution requirements to derive the correspond metric tensor and the stabilized parameters. From our derivation it is easily found that the optimal stabilized parameter is coupled with the optimal metric tensor on each element. Specifically, the relationship between the optimal metric tensor and the optimal stabilized parameter on each element is given approximately by (4.7).

The rest of the paper is organized as follows. In Section 2, we state the GLS stabilized finite element method for the problem (1.1). Section 3 is devoted to obtaining the estimate for the discretized error in a suitable norm via the anisotropic framework similar to that used in [51]. The optimal choice of the metric tensor and the stabilized parameters for the stabilized linear finite element methods are then derived in Section 4. Some numerical examples are provided in Section 5 to demonstrate the stability and efficiency of the proposed numerical strategy. Some concluding remarks will be given in the last section.

## 2 Stabilized finite element discretization

We shall use the standard notations in for the Sobolev spaces and their associated inner products , norms , and seminorms for . The Sobolev space coincides with , in which case the norm and inner product are denoted by and , respectively. Let and . The variational formulation of (1.1) reads as follows: find which satisfies

 A(u,v)=F(v),∀v∈H10(Ω), (2.1)

where and define the bilinear and linear forms

 A(u,v)=(ε∇u,∇v)+(b⋅∇u,v),

and

 F(v)=(f,v),

respectively.

Given a triangulation of , we denote the piecewise linear and continuous finite element space by , i.e.,

 Vh={v∈H1(Ω),v|K∈P1(K),∀K∈Th},

where is linear polynomial space in one element . We then define and . The standard finite element discretization of (2.1) is to find such that

 A(uh,vh)=f(vh),∀vh∈Vh0. (2.2)

For convection-dominated problems () , (2.2) using standard grid sizes are not able to capture steep layers without introducing non-physical oscillations. To enhance the stability and accuracy in the convection dominated regime, various stabilization strategies have been developed. Here we take the GLS stabilized finite element method as an example which reads as follows: find such that

 Ah(uh,vh)=F(vh),∀vh∈Vh0, (2.3)

with

 Ah(uh,vh)=A(uh,vh)+∑K∈ThαK(−εΔuh+b⋅∇uh,−εΔvh+b⋅∇vh)K,

and

 Fh(vh)=F(vh)+∑K∈ThαK(f,−εΔvh+b⋅∇vh)K.

In this paper we use linear finite element method, so the terms and in the two above equations are identically equal to zero. At this time, the GLS approach is the same as the SUPG method, which also enjoys the result in this paper if the linear finite element method is used. We endow the space with the discrete norm defined, for any , by

 ∥w∥2h:=ε∥∇w∥2L2(Ω)+∑K∈ThαK∥b⋅∇w∥2L2(K). (2.4)
###### Lemma 2.1.

The stabilized finite element approximation defined by (2.3) has the following error estimate

 ∥u−uh∥2h ≤C∑K∈Th(α−1K∥u−Πhu∥2L2(K)+ε∥∇(u−Πhu)∥2L2(K) +αK∥b⋅∇(u−Πhu)∥2L2(K)+αKε2∥Δ(u−Πhu)∥2L2(K)), (2.5)

where denotes the standard continuous piecewise linear interpolation operator.

###### Proof.

See, for example, [23, 34, 45]. ∎

## 3 Estimates for the interpolation error and its variants

As stated in Lemma 2.1, the discretized error in norm is bounded by four terms of interpolation error and its variants in different norms. In fact the interpolation error depends on the solution, the size and shape of the elements in the mesh. Understanding this relation is crucial for generating efficient meshes. In the mesh generation fields, this relation is often studied for the model problem of interpolating quadratic functions. For instance, Nadler [47] derived an exact expression for the -norm of the linear interpolation error in terms of the three sides , , and of the triangle :

 ∥u−Πhu∥2L2(K)=|K|180[(d11+d22+d33)2+d211+d222+d233], (3.1)

where is the area of the triangle, with being the Hessian of .

Three element-wise error estimates in different norms are derived by the following lemmas, which, together with [47], are fundamental for further discussion. Suppose is a quadratic function on a triangle . The function is given by its matrix representation:

 ∀x∈K,u(x)=12xtH(u)x.
###### Lemma 3.1.

Let be a quadratic function on a triangle , and be the Lagrangian linear interpolation of on . The following relationship holds:

 ∥∇(u−Πhu)∥2L2(K)=148|K|∑i,j=1,2i≤jDijℓtiℓj, (3.2)

where

 D11=d212+d223,D22=d212+d213,D12=2d212. (3.3)
###### Proof.

The proof is similar to but easier than that of Lemma 3.2. ∎

###### Remark 3.1.

Lemma 3.1 is also used in [51].

###### Lemma 3.2.

Let be a quadratic function on a triangle , and be the Lagrangian linear interpolation of on . Assume is a constant vector on , the following relationship holds:

 ∥b⋅∇(u−Πhu)∥2L2(K)=148|K|∑i,j=1,2i≤jDijkikj, (3.4)

where

 k1=(b2,−b1)ℓ1,k2=(b2,−b1)ℓ2. (3.5)
###### Proof.

Following [42], we first define the reference element by its three vertices coordinates:

 ^x1=(0,0)t,^x2=(1,0)t,and^x3=(0,1)t.

All the terms are computed on and then converted onto the element at hand by using the following affine mapping:

 x=x1+BK^xwithBK=[ℓ1,−ℓ2],x∈K,^x∈Kr,

where

 ℓ1=x2−x1,andℓ2=x1−x3.

In the frame of , the quadratic function turns into:

 u(x(^x))=12xt1H(u)x1+12xt1H(u)BK^x+12^xtBtKH(u)x1+12^xtBtKH(u)BK^x.

Since the linear interpolation is concerned, linear and constant terms of are exactly interpolated, these terms are neglected and only quadratic terms are kept. So we could set , with a matrix form:

 u(x(^x))=12^xtBtKH(u)BK^x=12(^x^y)t[d11−d12−d21d22](^x^y).

Then the exact point-wise interpolation error of the function in reads:

 e(x(^x))=12[d11(^x2−^x)+d22(^y2−^y)−2d12^x^y]. (3.6)

It is obvious that the following formulas hold

 BK=[ℓ1,−ℓ2]=[x2−x1x3−x1y2−y1y3−y1],B−1K=1det(BK)[y3−y1x1−x3y1−y2x2−x1].

After that it is easily to obtain

 det(BK)⋅btB−tk=bt[[]ccy3−y1y1−y2x1−x3x2−x1]=(b2,−b1)[ℓ2,ℓ1]=(k2,k1).

Then we have

 ∫K(b⋅∇xe(x))2dxdy =det(BK)∫Kr(btB−tK∇^xe(x(^x)))2d^xd^y =1det(BK)[k22∫Kr(∂e(x(^x))∂^x)2d^xd^y+k21∫Kr(∂e(x(^x))∂^y)2d^xd^y +2k1k2∫Kr∂e(x(^x))∂^x∂e(x(^x))∂^yd^xd^y]. (3.7)

Due to (3.6), we can easily obtain

 ∇^xe(x(^x))=(∂e(x(^x))/∂^x∂e(x(^x))/∂^y)=12(d11(2^x−1)−2d12^yd22(2^y−1)−2d12^x).

After simple calculation, the following results hold:

 ∫Kr^x2d^xd^y=∫Kr^y2d^xd^y=112,   ∫Kr^x^yd^xd^y=124, ∫Kr^xd^xd^y=∫Kr^yd^xd^y=16,   ∫Kr1d^xd^y=12.

Then we have

 24∫Kr(∂e(x(^x))∂^x)2d^xd^y=d212+d213=D22, (3.8)
 24∫Kr(∂e(x(^x))∂^y)2d^xd^y=d212+d223=D11, (3.9)
 48∫Kr∂e(x(^x))∂^x∂e(x(^x))∂^yd^xd^y=2d212=D12. (3.10)

Substituting (3.8), (3.9), (3.10) into (3) we get the desired results (3.4) due to the fact . ∎

###### Lemma 3.3.

Let be a quadratic function on a triangle . The following relationship holds:

 ∥Δu∥2L2(K)=(|ℓ2|2d11−2ℓt1ℓ2d12+|ℓ1|2d22)216|K|3. (3.11)
###### Proof.

Similar to that of Lemma 3.2 we can easily obtain

 ∫K(Δu)2dxdy =det(BK)∫Kr(|ℓ2|2det(BK)2d11−2ℓt1ℓ2det(BK)2d12+|ℓ1|2det(BK)2d22)2d^xd^y =(|ℓ2|2d11−2ℓt1ℓ2d12+|ℓ1|2d22)22det(BK)3=(|ℓ2|2d11−2ℓt1ℓ2d12+|ℓ1|2d22)216|K|3.

This is the desired result (3.11) and the proof is complete. ∎

###### Remark 3.2.

Since is a constant under our assumption, there is a rather direct and easy way to prove Lemma 3.3. However, we still use the frame of proof for Lemma 3.2 to make the error expression be a consistent manner.

###### Theorem 3.1.

Assume the exact solution is quadratic on each element , the error of the stabilized finite element approximation has the following estimate

 ∥u−uh∥2h≤C∑K∈ThEK, (3.12)

with

 EK =|K|180αK[(3∑i=1dii)2+3∑i=1d2ii]+ε48|K|∑i,j=1,2i≤jDijℓtiℓj+αK48|K|∑i,j=1,2i≤jDijkikj +αKε216|K|3(|ℓ2|2d11−2ℓt1ℓ2d12+|ℓ1|2d22)2,

where and are defined by (3.3) and (3.5), respectively.

###### Proof.

Together with Lemmas 2.1, 3.1, 3.2 and 3.3, the conclusion is obtained directly. ∎

Even the error estimate (3.12) is only valid for those piecewise quadratic functions, however, it could catch the main properties of the errors for general functions. In fact, the treatment to replace the general solution by its second order Taylor expansion yields a reliable and efficient estimator of the interpolation error for general functions provided a saturation assumption is valid [1, 20].

For simplicity of notation in the following discussion, , we denote by

 Q1,K=|K|180[(3∑i=1dii)2+3∑i=1d2ii],Q2,K=148|K|∑i,j=1,2i≤jDijℓtiℓj,

And then in (3.12) can be recast into

 EK=Q1,K⋅α−1K+Q2,K⋅ε+˜Q2,K⋅αK+Q3,K⋅αKε2.

## 4 Metric tensors for anisotropic mesh adaptation

We now use the error estimates obtained in Section 3 to develop the metric tensor for the discretized error in norm and give a new definition of the stability parameters which are optimal in a certain sense. In the field of mesh generation, the metric tensor, , is commonly used such that an anisotropic mesh is generated as a quasi-uniform mesh in the metric space determined by . Mathematically, this can be interpreted as the following three requirements ([32]). 1. The shape requirement. The elements of the new mesh, , are (or are close to being) equilateral in the metric. 2. The size requirement. The elements of the new mesh have a unitary volume in the metric. 3. The equidistribution requirement. The anisotropic mesh is required to minimize the error for a given number of mesh points (or equidistribute the error on every element).

### 4.1 Optimal metric tensor and stabilized parameters

We derive the monitor function first. At this time, we just need the shape and equidistribution requirements. Assume is a symmetric positive definite matrix on every point and this restriction can be explained by Remark 2 in [42]. Set . Denoted by and the projection of and to the constant space on , and . Since is a symmetric positive definite matrix, we do the singular value decomposition , where is the diagonal matrix of the corresponding eigenvalues () and is the orthogonal matrix with rows being the eigenvectors of . Denote by and the matrix and the vector defining the invertible affine map from the generic element to the reference triangle (see Figure 1). Here we take as an equilateral triangle with one edge which has angle with the horizontal line. Let , then . Mathematically, the shape requirement can be expressed as

 |^ℓ1|=|^ℓ2|=|^ℓ3|=L,^ℓ1⋅^ℓ3|^ℓ1|⋅|^ℓ3|=^ℓ2⋅^ℓ3|^ℓ2|⋅|^ℓ3|=^ℓ1⋅^ℓ2|^ℓ1|⋅|^ℓ2|=cos(2π/3), (4.1)

where is a constant.

###### Theorem 4.1.

Under the shape requirement, the following results hold:

 Q1,K=L4|K|15C2K,Q2,K=L4tr(HK)32√3C2Kdet(HK)12,Q3,K=3L4tr(HK)216C2K|K|det(HK),

and

 ˜Q2,K=L4det(HK)1224C2K(A21λ1,K+A22λ2,K)=L432√3C2K⋅btHKbdet(HK)12,

where

 A=(A1A2)=R[01−10]b.
###### Proof.

From the definition for and the relation between and , we have

 Q1,K=|K|180C2K[(3∑i=1|^ℓi|2)2+3∑i=1|^ℓi|4]=L4|K|15C2K.

Similarly, using (4.1), the following reduction is straight-forward:

 Q2,K=148|K|∑i,j=1,2i≤jDijℓtiℓj=148|K|(|ℓ1|2(d212+d223)+|ℓ2|2(d212+d213)+2(ℓ1⋅ℓ2)d212) =148C2K|K|(|ℓ1|2((^ℓ1⋅^ℓ2)2+(^ℓ2⋅^ℓ3)2)+|ℓ2|2((^ℓ1⋅^ℓ2)2+(^ℓ1⋅^ℓ3)2)+2(ℓ1⋅ℓ2)(^ℓ1⋅^ℓ2)2) =L4det(HK)1296CK|^K|∑i,j=1,2i≤j(C−12KR−1Λ−12^ℓi)⋅(C−12KR−1Λ−12^ℓj)by |K|=|^K|CK√det(HK) =L4det(HK)1296C2K|^K|∑i,j=1,2i≤j(Λ−12^ℓi)⋅(Λ−12^ℓj), (4.2)

Since and , after simple calculation,

 Λ−12^ℓ1=L(λ−121,Kcosθ,λ−122,Ksinθ)t,Λ−12^ℓ2=−L(λ−121,Kcos(π3+θ),λ−122,Ksin(π3+θ))t,

and then the following three equalities hold:

 =L2(λ−11,Kcos2θ+λ−12,Ksin2θ), =L2(λ−11,Kcos2(π3+θ)+λ−12,Ksin2(π3+θ)), =−L2(λ−11,Kcosθcos(π3+θ)+λ−12,Ksinθsin(π3+θ)).

Thus, we get

which gives to

 L4det(HK)1296C2K|^K|∑i,j=1,2i≤j(Λ−12^ℓi)⋅(Λ−12^ℓj)=L6det(HK)12128C2K|^K|tr(HK)det(HK)=L4tr(HK)32√3C2Kdet(HK)12. (4.3)

Insert (4.3) into (4.1) and then the formula for is proved.

Similar calculation can produce the corresponding formula for :

 Q3,K =116|K|3(|ℓ2|2d11−2ℓt1ℓ2d12+|ℓ1|2d22)2=L416C2K|K|3(∑i,j=1,2i≤jℓtiℓj)2 =L416C4K|K|3⋅916(λ−11,K+λ−12,K)2=9L8256C4K|K|3tr(HK)2det(HK)2=3L4tr(HK)216C2K|K|det(HK).

To analyze the term , more patience should be paid. First, similar to that of ,

 ˜Q2,K=L496C2K|K|∑i,j=1,2i≤jkikj. (4.4)

Second, using , we have

 k1=(b2,−b1)ℓ1=C−12K(b2,−b1)R−1Λ−12^ℓ1=C−12KL(A1λ−121,Kcosθ+A2λ−122,Ksinθ).

Similarly, we have

Thus,

 ∑i,j=1,2i≤jkikj=34C−1KL2(A21λ−11,K+A22λ−12,K). (4.5)

Finally, insert (4.5) into (4.4), it is easily got that

 ˜Q2,K =L6