Superconvergence analysis of the linear finite ... - Semantic Scholar

Comment

Report 4 Downloads 71 Views

MATHEMATICS OF COMPUTATION Volume 84, Number 291, January 2015, Pages 89–117 S 0025-5718(2014)02846-9 Article electronically published on May 28, 2014

SUPERCONVERGENCE ANALYSIS OF THE LINEAR FINITE ELEMENT METHOD AND A GRADIENT RECOVERY POSTPROCESSING ON ANISOTROPIC MESHES WEIMING CAO Abstract. For the linear ﬁnite element method based on general unstructured anisotropic meshes in two dimensions, we establish the superconvergence in energy norm of the ﬁnite element solution to the interpolation of the exact solution for elliptic problems. We also prove the superconvergence of the postprocessing process based on the global L2 -projection of the gradient of the ﬁnite element solution. Our basic assumptions are: (i) the mesh is quasiuniform under a Riemannian metric and (ii) each adjacent element pair forms an approximate (anisotropic) parallelogram. The analysis follows the same methodology developed by Bank and Xu in 2003 for the case of quasi-uniform meshes, and the results can be considered as an extension of their conclusion to the adaptive anisotropic meshes. Numerical examples involving both internal and boundary layers are presented in support of the theoretical analysis.

1. Introduction Recovery type error estimators are widely used in computations based on the ﬁnite element method, particularly in various engineering applications, due to their simplicity and robustness. They are often asymptotically exact, i.e., the ratio of the estimated error over the true error tends to 1 as the mesh gets ﬁner and ﬁner. These desirable properties are closely associated with the superconvergence of the ﬁnite element solution and its postprocessing. There have been extensive studies on the superconvergence and the recovery type error estimators; see, e.g., monographs by Wahlbin [28], Ainsworth and Oden [1], Babuˇska and Strouboulis [2], and Zhu [35], for an overview on these topics. Since the superconvergence is generally resulted from the cancellation of certain lower order terms in the ﬁnite element approximation due to local mesh symmetry, its mathematical analysis typically begins with the FE approximations based on uniform meshes, then extends to discretizations on mildly structured meshes [30,31], unstructured quasi-uniform meshes [3, 4, 15], and more recently on adaptively reﬁned meshes [29]. For problems exhibiting anisotropic features, e.g., internal and boundary layers, FE solutions based on anisotropic meshes can be much more effective than those based on isotropic ones. Recovery type error estimators have also been applied to this type of computation [12, 14, 19, 25]. Numerical examples demonstrate that these estimators are still robust and reliable at the presence of Received by the editor June 30, 2012 and, in revised form, April 29, 2013. 2010 Mathematics Subject Classiﬁcation. Primary 65N30, 65N15, 65N50. Key words and phrases. Superconvergence, anisotropic mesh, recovery type error estimates, post processing, linear ﬁnite element. This work was supported in part by NSF grant DMS-0811232. c 2014 American Mathematical Society Reverts to public domain 28 years from publication

89

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

90

WEIMING CAO

highly anisotropic elements [20, 24]. There have been some theoretical studies on the superconvergence and recovery techniques on structured anisotropic meshes. For instance, Li and Wheeler [16] proved the superconvergence of the FE solution to singular perturbed equations on anisotropic tensor product meshes. Shi, Mao, and Chen [18, 27] published a number of papers on the superconvergence analysis on uniform meshes that stretch anisotropically along coordinate axes. Lin and Lin [17] showed the superconvergence of bilinear elements on anisotropic rectangular meshes. Zhang [32], Xu and Zhang [30] demonstrated that the superconvergence of the FE solution and the asymptotic exactness of their polynomial preserving recovery estimators hold on a class of structured anisotropic meshes. However, we have seen little theoretical work on the case of general adaptively reﬁned unstructured anisotropic meshes. The purpose of this paper is to provide an analysis of the superconvergence property for the linear ﬁnite element approximation of second order problems on adaptively reﬁned anisotropic meshes and an analysis of the recovery type error estimator based on the global L2 -projection. More speciﬁcally, we consider a class of triangulations that are quasi-uniform under a problem dependent Riemannian metric. By assuming each pair of adjacent elements in the partition forms an approximate parallelogram, we demonstrate that ∇uN − ∇uI is superconvergent, where uN and uI are the ﬁnite element solution and the linear interpolation of the exact solution, respectively, and N is the total number of elements. Furthermore, we prove that ∇u − P(∇uN ) is superconvergent, where P is the global L2 -projection onto the space of continuous piecewise linear functions. These conclusions imply immediately the asymptotic exactness of the recovery based error estimator ∇uN − P(∇uN ). This work is an extension of Bank and Xu’s results [3] on the recovery technique based on L2 -projection on quasi-uniform meshes. Needless to say, the convergence and superconvergence of the ﬁnite element solution and the error estimator based on its recovery relies intimately on the interplay between the anisotropic behaviors of the PDE solution and the anisotropic features of the meshes. Here we measure the anisotropic behavior of the PDE solution by using some non-negative matrices determined by its higher order directional derivatives introduced in [6, 7, 21, 22]. Then the error bounds can be expressed as a multiple of certain powers of N with the multiplicative constants being some integrals involving the mesh metrics and the anisotropic measures. A sketch of the paper is as follows: We ﬁrst extend the concept of approximate parallelograms to the anisotropic meshes. Then we describe the model problem, the anisotropic meshes, and the assumption about the anisotropic behavior of the PDE solutions. We prove the superconvergence of uN − uI 1 in Section 4 and the superconvergence of uN − P(∇uN ) in Section 5. Numerical examples are presented in Section 6. We conclude the paper with some discussions in Section 7.

2. Approximate parallelograms on anisotropic meshes Local mesh symmetry plays an essential role in the superconvergence analysis of the ﬁnite element approximations. For linear elements, this symmetry typically requires each pair of elements sharing a common edge to form a parallelogram or an approximate parallelogram. We generalize here the concept of O(h1+α )approximate parallelogram for quasi-uniform meshes introduced in [3, 15, 26].

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

SUPERCONVERGENCE ON ANISOTROPIC MESHES

91

First, recall that a parallelogram can be characterized in several ways, e.g., the ratio of opposite side lengths is 1, or the distance between the midpoints of two diagonals is 0, or each pair of opposite sides are in parallel. For a family {Th } of quasi-uniform triangulations, with h being the characteristic diameter of elements in Th , τ1 ∪ τ2 forms an O(h1+α ) approximate parallelogram, if the ratio of its opposite side lengths is of order O(hα ), or if the distance between the midpoints of two diagonals is O(h1+α ), or the angle between the outward normals on each pair of opposite sides are within O(hα ) range of π; see [3, 15, 26]. Instead of using the above geometric characters to deﬁne an approximate parallelogram, we use here the Jacobians of the aﬃne mapping from an equilateral standard element τˆ to τ1 and τ2 . Suppose τ1 is with vertices i, j, and k, and τ2 with vertices j, i, and k . Let J1 and J2 be the Jacobians of the aﬃne mappings from τˆ to τ1 and τ2 with vertices k and k being the images of the same vertex on τˆ, respectively. Then τ1 ∪ τ2 forms a parallelogram iﬀ J1 = −J2 , and τ1 ∪ τ2 forms an O(h1+α ) approximate parallelogram, iﬀ I + J1−1 J2 I + J1 J2−1 = O(hα ). In addition, {Th } is a regular family of partition iﬀ cond(Jτ ) = O(1) for all τ ∈ Th . For discretization involving anisotropic meshes, the element diameter is not a proper parameter for describing the asymptotic behaviors, since an element can have much diﬀerent length scales in diﬀerent directions. Instead, the total number N of elements is essential. Thus we introduce the following concept of the O(N −(1+α)/2 )-approximate parallelogram. Deﬁnition 2.1. Let τ1 and τ2 be a pair of elements sharing a common edge, and let J1 and J2 be the Jacobians of the aﬃne mappings from equilateral standard element τˆ to τ1 and τ2 . τ1 ∪ τ2 is called forming an O(N −(1+α)/2 )-approximate parallelogram, if I + J1−1 J2 = O(N −α/2 ).

(1)

When {TN } is quasi-uniform, we have h = diam(τ ) = O(N −1/2 ), and the above deﬁnition is reduced to the O(h1+α )-approximate parallelogram for shape regular elements described in [3, 15]. It should be noted that this deﬁnition is independent of the ordering of τ1 and τ2 . It is easy to see that condition (1) is equivalent to I + J2−1 J1 = O(N −α/2 ).

(2)

Indeed, we have from (1) that J2−1 J1 = (I − (I + J1−1 J2 ))−1 ≤

1 = O(1) 1 − I + J1−1 J2

Hence I + J2−1 J1 = J2−1 J1 (I + J1−1 J2 ) ≤ J2−1 J1 · I + J1−1 J2 = O(N −α/2 ). However, this deﬁnition is not equivalent to requiring that (3)

I + J1 J2−1 = O(N −α/2 ).

We give here two examples of the anisotropic element pairs, one forms an approximate parallelogram by satisfying (1), and another does not, even though it satisﬁes (3). We choose condition (1) over (3) since the former measures the deviation

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

92

WEIMING CAO

of τ1 ∪ τ2 from an exact parallelogram in the coordinate system for the standard element, which is ﬁxed for all elements and all mesh sizes. (a)

2ε

(b)

(1,ε)

ε

1 ε

1 ε

0

Figure 1. Left: a pair of elements that forms an approximate parallelogram; Right: a pair of elements that does not form an approximate parallelogram. Example 1. Let = N −α/2 → 0, and let the standard element be the equilateral triangle with vertices (cos θm , sin θm ), θm = π2 + 2m 3 π, m = 0, 1, 2. Choose x1 = (0, 0), x2 = (1 − , 0), x3 = (1, ), and x4 = (0, ). Elements τ1 = Δx4 x1 x3 , and τ2 = Δx2 x3 x1 ; see Figure 1(a). The center of τ1 is xc = ( 31 , 2 3 ). Then the aﬃne mapping from τˆ to τ1 (with (0, 1) to x4 ) is x = J1 ξ + xc with the Jacobian J1 for the mapping being 1 J1 = √ (x3 − x1 ), x4 − xc = 3

√1 3 √ 3

− 13

.

3

Similarly, the Jacobian for the aﬃne mapping from τˆ to τ2 (with (0, 1) to x2 ) is − √13 1−2 0 − 2 3 3 J2 = = −J1 + . − √3 − 3 0 0 It is elementary to verify that I+

J1−1 J2

=

J1−1

and that I+

J2 J1−1

=

·

0 − 2 3 0 0

0 − 2 3 0 0

·

0 − √3 0

=

J1−1

=

−1 0 0

, .

Thus (1) holds for τ1 ∪τ2 , and it forms an O(N −(1+α)/2 )-approximate parallelogram, even though the left and√right sides always form a 45◦ angle, and their lengths diﬀer by a constant multiple 2. Note that (3) is not satisﬁed in this case. Example 2. Let x1 = (−1, 0), x2 = (0, −), x3 = (1, 0), and x4 = (0, 2). Elements τ1 = Δx1 x2 x4 , and τ2 = Δx3 x4 x2 ; see Figure 1(b). The Jacobians for the mapping from τˆ to τ1 and τ2 with (1, 0) to x1 and x3 , respectively, are 2 2 0 0 3 √ 0 − 3 √ 0 J1 = + , J2 = = −J . 1 0 − 2 3 − 3 − 3 − 3 3 Clearly, I+

J1−1 J2

=

2 0 − 3√ 3 0 0

,

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

SUPERCONVERGENCE ON ANISOTROPIC MESHES

and I + J2 J1−1 =

0 0 0

93

.

Thus (1) does not hold for τ1 ∪τ2 , and it does not form an O(N −(1+α)/2 )-approximate parallelogram. Nevertheless, (3) is satisﬁed. To further elaborate the geometric meaning of the approximate parallelograms, we list in the following several lemmas. Lemma 2.1. If τ1 ∪ τ2 forms an O(N −(1+α)/2 )-approximate parallelogram, then |τ2 | = 1 + O(N −α/2 ). |τ1 | This conclusion comes from the fact that (1) implies |τ2 | = det(J1−1 J2 ) = det(−I + (I + J1−1 J2 )) = 1 + O(N −α/2 ). |τ1 | Lemma 2.2. A pair of adjacent elements τ1 and τ2 forms an O(N −(1+α)/2 )approximate parallelogram, if and only if there exists an aﬃne mapping F such that F (τ1 ) ∪ F (τ2 ) forms a shape regular O(N −(1+α)/2 )-approximate parallelogram. Proof. Note that the Jacobians for the aﬃne mapping from τˆ to F (τ1 ) and F (τ2 ) are J˜1 = F J1 and J˜2 = F J2 , resp. If F (τ1 )∪F (τ2 ) forms a shape regular O(N −(1+α)/2 )approximate parallelogram, then I + J˜−1 J˜2 = O(N −α/2 ), 1

which implies (1), and thus τ1 ∪ τ2 forms an O(N −(1+α)/2 )-approximate parallelogram. On the other hand, if τ1 ∪ τ2 forms an O(N −(1+α)/2 )-approximate parallelogram, choose F = J1−1 . Then F (τ1 ) is equilateral. F (τ2 ) is shape regular since (1) implies cond(F J2 ) ≤ C. Furthermore, I + J˜−1 J˜2 = I + J −1 J2 = O(N −α/2 ). 1

1

Thus F (τ1 )∪F (τ2 ) forms a shape regular O(N −(1+α)/2 )-approximate parallelogram. The above lemma indicates that an anisotropic O(N −(1+α)/2 )-approximate parallelogram has to be the aﬃne image of a shape regular O(N −(1+α)/2 )-approximate parallelogram. Note that for a pair of shape regular elements, if they form an O(h1+α )-approximate parallelogram, then another element pair obtained by swapping the diagonals of this approximate parallelogram also forms an O(h1+α )-approximate parallelogram. Namely, edge swapping in a mesh does not aﬀect whether two adjacent elements form an approximate parallelogram or not. This prompts us to introduce yet another equivalent condition for approximate parallelogram by using only the location of four vertices of the adjacent elements. ˆ be the standard quadrilateral element, Lemma 2.3. Let Q = τ1 ∪ τ2 , and let Q 2 ˆ to Q, and d is e.g., [0, 1] . JQ is the Jacobian of the bilinear mapping from Q the vector connecting two midpoints of the diagonals of Q. Then τ1 ∪ τ2 forms an O(N −(1+α)/2 )-approximate parallelogram, iﬀ ˆ | JQ−1 (ξ)d | = O(N −α/2 ), (4) ∀ξ ∈ Q.

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

94

WEIMING CAO

Proof. First we note that both conditions (1) and (4) are invariant under any aﬃne transform of τ1 ∪ τ2 . Thus we only have to show the above lemma in the cases τ1 and τ2 are shape regular. In particular, we may assume, without loss of generality, that Q is with vertices x1 = (0, 0), x2 = (1, 0), x3 = (1 + δ1 , 1 + δ2 ), and x4 = (0, 1). It is ready to verify that d = 12 (x1 + x3 − x2 − x4 ) = 12 [δ1 , δ2 ]T and that the bilinear ˆ to Q is mapping x = B(ξ) from Q (5)

B(ξ) = x1 + ξ1 (x2 − x1 ) + ξ2 (x4 − x1 ) + 2ξ1 ξ2 d

and ∂x = ∂ξ An elementary calculation leads to

JQ (ξ) =

JQ−1 (ξ)d =

1 + δ1 ξ 2 δ2 ξ 2

δ1 ξ 1 1 + δ2 ξ 1

1 2(1 + δ1 ξ2 + δ2 ξ1 )

δ1 δ2

. .

Since for all ξ ∈ [0, 1]2 , we have √ √ 1 − 2|d| ≤ 1 + δ1 ξ2 + δ2 ξ1 ≤ 1 + 2|d|. Thus (4) is satisﬁed iﬀ |d| = O(N −α/2 ), i.e., τ1 ∪ τ2 forms an O(N −(1+α)/2 )approximate parallelogram. We point out that to study the superconvergence of the ﬁnite element approximation on quadrilateral meshes, Zhang [33] showed that (4) is a necessary condition for a shape regular approximate parallelogram. The above lemma indicates that (4) is also a suﬃcient condition and can be used even in the case of anisotropic meshes. 3. Model problems and their FE approximation We consider the following homogeneous Dirichlet problem of a second order elliptic equation: −∇ · (A∇u + b u) + d u = f, in Ω, (6) u|∂Ω = 0. We assume in this paper A is a positive deﬁnite constant matrix, b, d are suitably smooth functions and equation (6) is strongly elliptic. FE approximation: Let {TN } be a family of regular partitions of Ω [9]. Here, N represents the total number of elements in TN . Deﬁne SN be the space of continuous piecewise linear functions based on partition TN . VN = SN ∩H01 (Ω). Then the ﬁnite element method for solving (6) is to ﬁnd the approximate solution uN ∈ VN such that (7)

a(uN , vN ) = (f, vN ),

∀vN ∈ VN ,

where bilinear form a(·, ·) is deﬁned as a(u, v) = (A∇u, ∇v) + (u, b · ∇v) + (du, v). Anisotropic meshes: Since the convergence and superconvergence of the ﬁnite element solution depend closely on the properties of meshes, we restrict our analysis on a class of anisotropic meshes that are quasi-uniform under a given metric. More

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

SUPERCONVERGENCE ON ANISOTROPIC MESHES

95

¯ For each element τ ∈ TN , precisely, Let M be a continuous Riemannian metric on Ω. let Mτ be the average of M over τ , and its eigen-decomposition is of the form Mτ = Tτ · Λτ · Tτ ,

(8)

where Λτ is diagonal and Tτ is orthonormal. Deﬁne −1

Fτ = Tτ Λτ 2 .

(9)

A family of triangulations {TN } is called quasi-uniform under metric M , if for all τ ∈ TN , τ˜ = Fτ−1 τ are shape regular and of about the same size; see [6–8]. Let Jτ be the Jacobian of the aﬃne mapping from standard element τˆ to element τ . Then {TN } is quasi-uniform under metric M if and only if (10)

Fτ−1 Jτ (Jτ−1 Fτ )−1 (CM /N )1/2 ,

where (11)

∀τ ∈ TN ,

| det(M )|1/2 .

CM = Ω

Clearly, the error for the ﬁnite element approximation depends on the interplay between the anisotropic features of the meshes and the anisotropic behavior of the solutions. In order to derive the convergence and superconvergence for the linear ﬁnite element approximation, we need the following assumption to characterize the anisotropic behavior of the second and third derivatives of solution u. Assumption on D2 u and D3 u: Let Q2 (x) and Q3 (x) be 2 × 2 symmetric non-negative deﬁnite matrices. For m = 2, 3, we assume at each x ∈ Ω that (12)

|(s · ∇)m u(x)| ≤ (s · Qm (x)s)m/2 ,

∀s ∈ R2 .

Such matrices always exist as long as D2 u and D3 u exist. For instance, Q2 (x) can be chosen as |∇2 u| I, where I is the identity matrix and |∇2 u| = maxs=1 |(s · ∇)2 u|. It can also be chosen as abs(∇2 u(x)), which is the matrix of the same eigenvectors as ∇2 u(x) but with the eigenvalues equal to the absolute values of ∇2 u(x)’s. Note that the left-hand side of (12) is the directional derivative along s, the “smaller” the matrix Qm , the more precise it describes the anisotropic behaviors of Dm u. A geometric explanation of (12) and a numerical procedure to ﬁnd these matrices can be found in [6, 7]. More recently, Mirebeau studied systematically how to ﬁnd these matrices and applied them in anisotropic mesh generation. He also derived an explicit formula to deﬁne Q3 ; see [21, 22] for details. The eﬀects of the anisotropic mesh for approximating a solution of anisotropic derivatives can be appreciated by the following inequalities. Let u ˆ be the pull-back of u|τ under the aﬃne mapping from the standard element τˆ to element τ . Then ˆ = JτT ∇, and we have ∇ (13)

ˆ mu |(s · ∇) ˆ| ≤ (s · JτT Qm Jτ s)m/2 ,

∀s ∈ R2

which implies (14)

ˆ mu |D ˆ| ≤ JτT Qm Jτ m/2 .

When Dm u is highly anisotropic, the eigenvalues of Qm are of very diﬀerent magˆ mu nitudes. A proper choice of Jτ may make |D ˆ| signiﬁcantly smaller than |Dm u|. a-priori error estimate: By the ellipticity of equation (6), it follows from the interpolation error estimates (see Theorem 2.1 in [7]) the following a priori error estimate for the linear ﬁnite element approximation of (6).

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

96

WEIMING CAO

Theorem 3.1. Suppose {TN } is quasi-uniform under metric M , and uN is the linear ﬁnite element approximation of the solution u of (6). Furthermore, assume the second derivative D2 u satisﬁes assumption (12) about its anisotropic behavior. Then there exists a constant c independent of u and N such that

1/2 (15) F −1 2 F T Q2 F 2 , u − uN 1,Ω ≤ cN −1/2 · CM Ω

where F and CM are determined by M as in (9) and (11). Remark 3.1. When ∇2 u is positive deﬁnite on Ω, we may choose Q2 = ∇2 u in (12), and choose M = c(λmax /λmin )1/4 · ∇2 u, where λmax and λmin are the largest and smallest eigenvalues of ∇2 u, respectively; see [5, 7]. In this case, (15) is reduced to u − uN 1,Ω ≤ cN −1/2 λmin λ3/4 max L1 (Ω) . 1/4

In contrast, if we allow only shape regular elements, then the best isotropic mesh metric to minimize the error bound in (15) would be M = λmax I, and the error bound for u − uN 1,Ω would be cN −1/2 λmax L1 (Ω) , which is clearly much larger than for the anisotropic case when λmax and λmin are of much diﬀerent magnitudes. 4. Superconvergence of FE solution First we recall a basic lemma developed by Bank and Xu [3] for analyzing the “weak convergence” of linear interpolation on a triangular element. Let τ be an element with vertices xk , k = 1, 2, 3. Denote by θk the internal angle of τ at xk . Let k = xk−1 − xk+1 be the edge (directed counterclockwise) opposite to xk , and

k = |k | is its length. tk = k / k and nk are the unit tangential and outward normal directions on k , respectively, and ξk = −nk+1 · Ank−1 . Furthermore, let φk be the nodal basis function associated with xk , qk = φk+1 φk−1 , and ψk = φk (1 − φk ). Lemma 4.1. Let τ be an element in TN , and let uI be the linear interpolation of u at the vertices of τ . Then (16) 3 ξk qk ∂2u ∂2u ∇(u − uI ) · A∇vN = · {( 2k+1 − 2k−1 ) 2 + 4|τ | } 2 sin θk ∂tk ∂nk ∂tk τ k=1 ek 3 ξk k ∂3u ∂3u ∂vN · − {

ψ +

ψ }. k+1 k−1 k−1 k+1 2 ∂tk ∂ 2 tk+1 ∂tk−1 ∂ 2 tk−1 ∂tk+1 τ k=1 2 sin θk This lemma plays a fundamental role in the superconvergence analysis [3, 30]. It separates not only the error contribution from each edge of the elements, but more importantly, the contributions from normal and tangential derivatives. It makes the cancellation of lower order error terms much easier. Moreover, (16) is an identity valid on any triangle. We shall use it on the standard element τˆ. Theorem 4.1. Let uI and uN be the linear interpolation and the ﬁnite element approximation of u, respectively. Suppose D2 u and D3 u satisfy assumption (12),

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

SUPERCONVERGENCE ON ANISOTROPIC MESHES

97

and each pair of adjacent elements in TN forms an O(N −(1+α)/2 )-approximate parallelogram. Then

uN − uI 1,Ω ≤ c (17) (1 + N −α Jτ−1 4 Jτ 2 ) · JτT Q2 Jτ 2 τ

τ

+ Jτ−1 4 Jτ 2 · JτT Q3 Jτ 3

1/2 .

In addition, if the partition TN is quasi-uniform under metric M , then we have

−1/2 uN − uI 1,Ω ≤ cN (18) ((CM )2 N −1 + CM N −α F −1 4 F 2 ) Ω

· F T Q2 F 2 + (CM )2 N −1 F −1 4 F 2 · F T Q3 F 3

1/2 ,

where F and CM are determined by M as in (9) and (11). Proof. By the ellipticity of a(·, ·) and the orthogonality of the ﬁnite element approximation, uN − uI 1,Ω ≤ c sup vN

|a(uN − uI , vN )| |a(u − uI , vN )| = c sup . vN 1,Ω vN 1,Ω vN

To prove the theorem, we only have to bound a(u − uI , vN ) by a certain multiple of vN 1 for arbitrary vN ∈ SN . Note that (19) a(u − uI , vN ) = ∇(u − uI ) · A∇vN + (u − uI ) (b · ∇vN + dvN ). τ

τ

Ω

First we estimate the second term on the right-hand side of the above equation. On each element τ ∈ TN , let Fτ be the aﬃne mapping from the equilateral standard ˆ 1u element τˆ to τ , u ˆτ = u ◦ Fτ , and let u ˆI = Π ˆτ be the linear interpolation of u ˆτ on τˆ. It follows from assumption (12) on the anisotropic behavior of D2 u that | (u − uI ) · (b · ∇vN + dvN ) | ≤ c u − uI vN 1 Ω

1/2 ≤ c vN 1 u − uI 2 τ τ

1/2 ≤ c vN 1 |τ | · ˆ u−u ˆI 2 τ τˆ

1/2 ˆ 2u ≤ c vN 1 |τ | · |D ˆ |2 τ τˆ

1/2 T ≤ c vN 1 Jτ Q2 Jτ 2 . τ

τ

Now we deal with the ﬁrst term on the right-hand side of (19). On each element τ , let Jτ be the Jacobian of the mapping Fτ . Let Aˆτ = Jτ−1 AJτ−T and ξˆk = ˆ = JτT ∇ and −ˆ nk+1 · Aˆ nk−1 . Then by a change of variables, we have ∇ ˆ uτ − u ˆ vN . (20) ∇(u − uI ) · A∇vN = |τ | · ∇(ˆ ˆI ) · Aˆτ ∇ˆ τ ∈TN

τ

τ ∈TN

τˆ

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

98

WEIMING CAO

Applying Lemma 4.1 to the integrals on τˆ for each τ , and noting that |ˆ τ| = √ π ˆ |k | = 3 and θk = 3 for all k = 1, 2, 3, we have 3 ˆτ ∂2u ∂ˆ vN ˆ uτ − u ˆ vN = ∇(ˆ (21) ˆI ) · Aˆτ ∇ˆ · ξˆk qk · ∂ tˆk ∂ n ˆ k ∂ tˆk τˆ k=1 eˆk 3 ∂3u ˆτ ∂3u ˆτ ∂ˆ vN −2 + ψk−1 2 )· . ξˆk (ψk+1 2 ˆ ˆ ˆ ˆ ∂ tk+1 ∂ tk−1 ∂ tk−1 ∂ tk+1 ∂ tˆk k=1 τˆ

√ 3 3 4 ,

Summing up the above integrals for all τ ∈ TN and grouping together the two line integrals associated with a common edge e = τ ∩ τ , we have ∇(u − uI ) · A∇vN = I1 + I2 , τ ∈TN

τ

where

3

ˆτ ∂2u ∂ˆ vN · ξˆk qk · ˆ ∂ t ∂ n ˆ ∂ tˆk k k τ ∈TN k=1 eˆk ∂2u ˆτ ∂ˆ ∂2u ˆτ ∂ˆ vN |τ vN |τ ˆ ˆ = |τ | q ξτ · + |τ | q ξτ · , ∂ tˆ∂ n ˆ| ∂ tˆ ∂ tˆ∂ n ˆ ∂ tˆ eˆ eˆ e=τ ∩τ 3 3 3 ∂ u ˆ ∂ u ˆ ∂ˆ v τ τ N I2 = −2 + ψk−1 2 ) · |τ | . ξˆk ( ψk+1 2 ∂ tˆk+1 ∂ tˆk−1 ∂ tˆk−1 ∂ tˆk+1 ∂ tˆk τ ∈T k=1 τˆ

I1 =

|τ | ·

N

We ﬁrst bound the integrals in I2 . Note that |ψk | ≤ 14 , |qk | ≤ 1, |ξˆk | ≤ Aˆτ ≤ Jτ−1 2 A and ∂ˆ vN ˆ vN | ≤ Jτ |∇vN |. | ≤ |∇ˆ ∂ tˆk Furthermore, by assumption (12), we have |

| Thus |I2 | ≤ c

∂3u ˆτ | 2ˆ ∂ tk+1 ∂ tˆk−1

τ ∈TN

≤c

ˆ 3u ≤ D ˆτ ≤ JτT Q3 Jτ 3/2 .

Jτ−1 2 Jτ · JτT Q3 Jτ 3/2 |∇vN |

τ

Jτ−1 4 Jτ 2 JτT Q3 Jτ 3

1/2

· ∇vN L2 (Ω) .

τ ∈TN τ

Next, we deal with I1 . Note that for each term in the sum up to I1 we may assume eˆ = eˆ , i.e., the common edge e of τ and τ is the images of the same edge in the equilateral standard element under Fτ and Fτ . Otherwise, a rotation of the standard element by angle ±120◦ for one of the integrals will convert it to such a case. With this in mind, we have ∂ˆ vN |τ ∂ˆ vN |τ =− ∂ tˆ ∂ tˆ

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

SUPERCONVERGENCE ON ANISOTROPIC MESHES

and

99

2 ∂2u ˆτ ˆτ ∂ˆ vN |τ ˆ ∂ u ˆ I1 = − |τ |ξτ ·q· . |τ |ξτ ˆ ˆ ∂ t∂ n ˆ ∂ t∂ n ˆ ∂ tˆ e=τ ∩τ eˆ

We split the integrand into (22)

2 ˆτ − |τ |ξˆ ∂ 2 u ˆτ = |τ |(ξˆ − ξˆ ) ∂ 2 u ˆτ |τ |ξˆτ ∂ u τ τ τ ˆ ˆ ˆ∂ n ∂ t∂ n ˆ ∂ t ∂ n ˆ ∂ t ˆ 2 ˆτ − ∂ 2 u ˆτ ) + (|τ | − |τ |)ξˆ ∂ 2 u ˆτ . + |τ |ξˆτ ( ∂ u τ ∂ tˆ∂ n ˆ ∂ tˆ∂ n ˆ ∂ tˆ∂ n ˆ

It follows from the assumption that τ ∪ τ forms an O(N −(1+α)/2 )-approximate parallelogram such that |ξˆτ − ξˆτ | = |ˆ nk+1 · (Aˆτ − Aˆτ )ˆ nk−1 | ≤ Aˆτ − Aˆτ −T = Jτ−1 AJτ−T − Jτ−1 AJτ

−T −T + Jτ−1 + Jτ− T ) ≤ (Jτ−1 + Jτ−1 )AJτ A(Jτ −1 + Jτ−1 ≤ Jτ−1 + Jτ−1 · A · (Jτ )

−1 · A · Jτ−1 (1 + Jτ−1 ≤ (I + Jτ−1 Jτ )Jτ Jτ )

≤ CN −α/2 Jτ−1 2 . On the other hand, ∂2u ˆτ ˆ n · ∇)ˆ ˆ uτ = (Jτ tˆ · ∇)(Jτ n ˆ · ∇)u = (tˆ · ∇)(ˆ ∂ tˆ∂ n ˆ and

∂2u ˆτ ˆ n · ∇)ˆ ˆ uτ = (Jτ tˆ · ∇)(Jτ n ˆ · ∇)u. = (tˆ · ∇)(ˆ ˆ ∂ t∂ n ˆ By the fact that Jτ tˆ = −Jτ tˆ, we have 2 ˆτ − ∂ 2 u ˆτ | = |(J tˆ · ∇)((J + J )ˆ | ∂ˆu τ τ τ n · ∇)u| ∂ t∂ n ˆ ∂ tˆ∂ n ˆ ˆ τ−1 (Jτ + Jτ )ˆ ˆ uτ | = |(tˆ · ∇)(J n · ∇)ˆ −1 2 ˆ u ≤ |(Jτ (Jτ + Jτ )ˆ n| · |D ˆτ |

ˆ 2u ≤ I + Jτ−1 Jτ · |D ˆτ | −α/2 ˆ 2 ≤ cN |D u ˆτ |. Put the above two bounds and the bound from Lemma 2.1 into (22), we have ˆ vN | · |D ˆ 2u |τ | · Jτ−1 2 · |∇ˆ ˆτ | |I1 | ≤ cN −α/2 e ˆ e=τ ∪τ

1/2 −α/2 −1 2 ˆ 2u ˆ 3u |τ | · Jτ · Jτ · |∇vN | · (|D ˆτ |2 + |D ˆ τ |2 ) ≤ cN τˆ τ ∈TN −α/2 1/2 −1 2 ≤ cN |τ | · Jτ · Jτ · |∇vN | τ ∈T N

1/2 · (JτT Q2 Jτ 2 + JτT Q3 Jτ 3 ) τ

1/2 −α/2 ≤ cN Jτ−1 4 Jτ 2 (JτT Q2 Jτ 2 +JτT Q3 Jτ 3 ) ·∇vN L2 (Ω) , τ ∈TN τ

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

100

WEIMING CAO

where in the second step we used the trace theorem to convert the integral on edge eˆ into the integral on element τˆ. Combining the estimates for I1 , I2 , and (4), we have for any vN ∈ SN ,

(23) |a(u − uI , vN )| ≤ c ∇vN · (1 + N −α Jτ−1 4 Jτ 2 ) τ

τ

· JτT Q2 Jτ 2 + Jτ−1 4 Jτ 2 · JτT Q3 Jτ 3

1/2 ,

which completes the proof of (17). Finally, by the assumption that TN is quasi-uniform under metric M , we have from (10) that Jτ (CM N −1 )1/2 Fτ . By putting it into (17) and use the continuity of M , we have the error bound (18) in integral form. Remark 4.1. Similar to Remark 3.1, if ∇2 u is positive deﬁnite on Ω, and the triangulation is quasi-uniform under metric M = c(λmax /λmin )1/4 · ∇2 u, we may also derive an error bound for uN − uI 1,Ω in terms of the eigenvalues of ∇2 u and Q3 . In particular, if α < 1, then the leading term on the right-hand side in (18) is of 1/4 3/4 −3/4 7/4 magnitude cN −(1+α)/2 (λmin λmax L1 (Ω) · λmin λmax L1 (Ω) )1/2 . 5. Superconvergence of postprocessing Let P be the global L2 -projection onto SN . We discuss in this section the superconvergence of P(∇uN ) to ∇u. A key step leading to this result is the superconvergent error estimate for ∇u − P(∇uI ). To obtain this estimate, we ﬁrst derive an identity for the diﬀerence ∇uq − ∇uI , where uq is the piecewise quadratic “projection based interpolation” of u deﬁned as 3 6 αk qk with αk = (u − uI ), uq = uI +

k ek k=1 where qk = φk+1 φk−1 is the side basis associated with edge ek . Clearly, ek (u−uq ) = 0 on each edge ek . In addition, uq is continuous on Ω, and vanishes on ∂Ω if u does. Lemma 5.1. For any u with u|∂Ω = 0 and any w ∈ (SN )2 , we have

1 1 (∇uq − ∇uI , w)Ω = (24) (w(mk ) · nk ) k +

k+1 ek+1 k−1 ek−1 ek ∈Ω 1 1 − − (u − uI )

k +1 ek +1 k −1 ek −1

1 1 (w(mk ) · nk ) k + + (u − uI ),

k+1 ek+1 k−1 ek−1 ek ∈∂Ω

where the summations are over all interior edges ek ∈ Ω and boundary edges ek ∈ ∂Ω in the partition, respectively. w(mk ) is the value of w at the midpoint mk of each ek , refer to Figure 2.

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

SUPERCONVERGENCE ON ANISOTROPIC MESHES

101

e k+1

k

τ

e k–1

mk

ek

τ’

’ –1

e k +1

k’

’

Figure 2. Element pair around an edge ek = τ ∪ τ .

Proof. By the facts that u|∂Ω = 0 and that w is linear on each element, we have (25)

(∇uq − ∇uI , w)Ω = −(uq − uI , ∇ · w)Ω = −

(∇ · w|τ ) (uq − uI , 1)τ . τ

Note that (qk , 1)τ = (26)

1 12 |τ |,

(uq − uI , 1)τ = (

we have

3

αk qk , 1)τ =

k=1

On the other hand, let w|τ =

3

|τ | |τ | 1 αk = 12 2

k 3

3

k=1

k=1

(u − uI ). ek

wk φk . By ∇φk = −(∇φk−1 + ∇φk+1 ), we have

k=1

∇ · w|τ =

3

w k · ∇φk = −

k=1

3

w k · (∇φk−1 + ∇φk+1 ) = −

k=1

3

(wk−1 + w k+1 ) · ∇φk .

k=1

Furthermore, by ∇φk = − k nk and w(mk ) = 12 (w k+1 + w k−1 ), 2|τ | 1

k w(mk ) · nk . |τ | 3

∇ · w|τ = −

(27)

k=1

By putting the above equation into (25), we have 3 1 (uq − uI , 1)τ ·

k (w(mk ) · nk )|τ |τ | τ k=1 1 1 (uq − uI , 1)τ − (uq − uI , 1)τ = (w(mk ) · nk |τ ) k |τ | |τ |

(∇uq − ∇uI ), w)Ω =

ek ∈Ω

+

ek ∈∂Ω

k (w(mk ) · nk |τ ) (uq − uI , 1)τ , |τ |

where in the second step we used the fact that nk |τ = −nk |τ , and w is continuous across element edges. Finally, by putting (26) into the right-hand side of the above equation, we have formula (24).

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

102

WEIMING CAO

Theorem 5.1. Suppose that D2 u and D3 u satisfy assumption (12), and that each pair of adjacent elements in {TN } forms an O(N −(1+α)/2 )-approximate parallelogram. Then ∇u − P∇uI ≤ c

τ

(28)

Jτ−1 2 · (N −α JτT Q2 Jτ 2 + JτT Q3 Jτ 3 ) τ 1/2

+c Jτ−1 2 · JτT Q2 Jτ 2 , τ ∈Ωbh

1/2

τ

where Ωbh = τ ∩∂Ω =∅ τ is the union of elements next to the boundary ∂Ω. In addition, if partition TN is quasi-uniform under metric M , then we have (29)

∇u − P∇uI ≤ cN

−1/2

F −1 2 · (CM N −α F T Q2 F 2

Ω

1/2 + (CM )2 N −1 F T Q3 F 3 )

1/2 −1/2 CM + cN F −1 2 · F T Q2 F 2 , Ωbh

where F and CM are determined by M as in (9) and (11). Proof. By the triangle inequality, we have (30)

∇u − P∇uI ≤ ∇u − P∇u + P(∇u − ∇uq ) + P(∇uq − ∇uI ).

The ﬁrst term on the right-hand side can be bounded by the interpolation error estimate. Note on each element τ , ˆ u) = Jτ−1 (Π ˆ u). ˆ 1 (Jτ−1 ∇ˆ ˆ 1 ∇ˆ Π1 (∇u) = Π Hence

∇u − P∇u ≤ ∇u − Π1 (∇u) = |∇u − Π1 (∇u)|2 τ τ

ˆ u − Jτ−1 (Π ˆ u)|2 ˆ 1 ∇ˆ = |τ | |Jτ−1 ∇ˆ τ τˆ

ˆ u − Πˆ1 (∇ˆ ˆ u)|2 ≤ |τ | Jτ−1 2 · |∇ˆ τ τˆ

ˆ 3u ≤ |τ | Jτ−1 2 · |D ˆ |2 τ τˆ

≤ |τ | Jτ−1 2 · JτT Q3 Jτ 3 τ τˆ

= Jτ−1 2 · JτT Q3 Jτ 3 , 2

2

τ

τ

where in the second last step, we used the assumption (12) about the anisotropic behavior of D3 u.

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

SUPERCONVERGENCE ON ANISOTROPIC MESHES

103

The second term on the right-hand side of (30) can also be bounded by using the error estimates for quadratic interpolation,

2 2 P(∇u − ∇uq ) ≤ ∇u − ∇uq = |∇u − ∇uq |2 τ τ

ˆ u − ∇ˆ ˆ uq |2 ≤ |τ | Jτ−1 2 |∇ˆ τ τˆ

ˆ 3u ≤ |τ | Jτ−1 2 |D ˆ |2 τ τˆ

≤ Jτ−1 2 JτT Q3 Jτ 3 . τ

τ

Now we estimate the third term P(∇uq − ∇uI ). By the fact that (31)

P(∇uq − ∇uI ) = =

sup

w∈(SN )2 sup

w∈(SN

)2

|(P(∇uq − ∇uI ), w)| w |(∇uq − ∇uI , w)| , w

we only have to bound |(∇uq − ∇uI , w)| by a certain multiple of w for any w ∈ (SN )2 . It follows from Lemma 5.1 that |(∇uq − ∇uI , w)|

1 1 |w(mk )| k (u − uI ) − (u − uI ) ≤

k+1 ek+1

k +1 ek +1 ek ∈Ω 1 1 (32) + (u − uI ) − (u − uI )

k−1 ek−1

k −1 ek −1

1 1 + |w(mk )| k (u − uI ) + (u − uI ) .

k+1 ek+1

k−1 ek−1 ek ∈∂Ω

In order to estimate the diﬀerence between the averages of u − uI on opposite edges ˆ = [0, 1]2 under a bilinear ek+1 and ek +1 , we consider Q = τ ∪ τ as the image of Q transform x = B(ξ). Then ek+1 and ek +1 are the images of two opposite sides of ˆ say η = 1 and η = 0. Furthermore, let u Q, ˜ = u ◦ B and JQ = ∂x . Thus ∂ξ 1 1 (u − uI ) − (u − uI )

k+1 ek+1

k +1 ek +1 1 1 = (˜ u − Π1 u ˜)(ξ, 1) dξ − (˜ u − Π1 u ˜)(ξ, 0) dξ 0 0 (33) 1 1 ˜ ˜ ∂2u ∂2u ξ(ξ − 1) 2 (ξ, 1) dξ − ξ(ξ − 1) 2 (ξ, 0) dξ = ∂ξ ∂ξ 0 0 ˜ ∂3u = ξ(ξ − 1) 2 (ξ, η) dηdξ. ˆ ∂ξ ∂η Q ˜| ≤ Now we bound the third order derivative in the above integral. Note that |D3 u 3 ∂ u ˜ sup|s|=1 | 3 |. By using the tensor form of the chain rule for derivatives (see, e.g., ∂s Theorem 37.1 in [9]), we have for any unit vector s = [s1 , s2 ]T , ˜ ∂3u = D3 u(xs , xs , xs ) + 3 D2 u(xs , xss ) + Du(xsss ). ∂s3

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

104

WEIMING CAO

Here xs =

∂x ∂x = · s = JQ s ∂s ∂ξ

3 and xsss = ∂ x = 0, since B is bilinear. By formula (5) for B(ξ), ∂s3

∂2x = s1 s2 d, ∂s2 where d is the vector connecting the midpoints of two diagonals of Q. Therefore, we have xss =

|D3 u(xs , xs , xs )| ≤ |D3 u(JQ s, JQ s, JQ s)| ≤ [JQ s · Q3 JQ s]3/2 ≤ JQT Q3 JQ 3/2 and |D2 u(xs , xss )| ≤ |JQ s · D2 u d| ≤ JQT D2 JQ |JQ−1 d| ≤ cN −α/2 JQT Q2 JQ , where in the last step we used Lemma 2.3 under the assumption that Q = τ ∪ τ forms an O(N −(1+α)/2 )-approximate parallelogram. Combining them together, we have ˜| ≤ cN −α/2 JQT Q2 JQ + JQT Q3 JQ 3/2 . |D3 u Therefore, we may derive from (33) that 1 1 (u − uI ) − (u − uI )

k+1 ek+1

k +1 ek +1 ≤ (cN −α/2 JQT Q2 JQ + JQT Q3 JQ 3/2 ) ˆ Q ≤ (cN −α/2 JQT Q2 JQ + JQT Q3 JQ 3/2 ) · det(JQ−1 ). Q

Note that Jτ−1 JQ (ξ) ≤ const,

ˆ ∀ξ ∈ Q,

which is obviously true in the case when τ and Q are shape regular. It is also true when they are anisotropic, because Jτ−1 JQ is invariant under any aﬃne transformation of Q. Therefore, JQT Q2 JQ ≤ c JτT Q2 Jτ ,

on τ.

The term JQT Q3 JQ can be bounded similarly. In addition, c . det(JQ−1 ) ≤ c det(Jτ−1 ) ≤ |τ | Hence we have 1

k+1

(u − uI ) − ek+1

c ≤ |τ |

1

k +1

(u − uI )

ek +1

(N −α/2 JτT Q2 Jτ + JτT Q3 Jτ 3/2 ).

τ

Put it into (32) and note that 4 4 Jτ 4 |Jτ eˆk |

k = √ · ≤ · ≤ Jτ−1 , |τ | 3 det(Jτ ) 3 3 3 det(Jτ )

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

SUPERCONVERGENCE ON ANISOTROPIC MESHES

105

then we have |(∇uq − ∇uI , w)| ≤ c

k=1

· (34)

3

k |w(mk )|) ( |τ | τ

(N

−α/2

Jτ Q2 Jτ + JτT Q3 Jτ 3/2 ) T

τ

3 ≤ c ( |w(mk )|) · Jτ−1 · (N −α/2 JτT Q2 Jτ + JτT Q3 Jτ 3/2 ) τ

τ k=1

≤ c w ·

Jτ−1 2 · (N −α JτT Q2 Jτ 2 + JτT Q3 Jτ 3 )

1/2 ,

τ

τ

where we have used the fact that 3 2 ( |w(mk )|) ≤ c |w|2 . τ

τ

k=1

This completes the estimate for the ﬁrst sum on the right-hand side of (32). For the second sum on the right-hand side of (32) involving edges on the domain boundary, note that u − uI = 0 on ∂Ω, thus we have for any linear function p1 , 1 1 (u − uI ) + k−1 (u − uI ) = (ˆ u−u ˆI ) k+1 ek+1 ek−1 ∂ τˆ ˆ 1 (ˆ = (ˆ u + p1 − Π u + p1 )) ∂ τˆ

u + p1 H 2 (ˆτ ) . ≤ cˆ u + p1 C(ˆτ ) ≤ cˆ Since p1 is arbitrary, it follows from Bramble-Hilbert Lemma [9] that 1 1 (u − u ) + (u − uI ) ≤ c inf ˆ u + p1 H 2 (ˆτ ) ≤ c |ˆ u|H 2 (ˆτ ) I k+1 k−1 ∀p1 ek+1 ek−1

1/2 ˆ 2u = c |D ˆ |2 1/2

τˆ 1 J T Q J 2 ≤ c . 2 τ |τ | τ τ Therefore, we have

1 |w(mk )| k k+1

(u − uI ) + k−1 (u − uI ) ek+1 ek−1 ek ∈∂Ω 1/2

2

k J T Q J 2 ≤ c |w(mk )| |τ |1/2 · 2 τ τ 2 τ |τ | ek ∈∂Ω 1/2

(35) 1/2 2 ≤ c |w(mk )| |τ | · Jτ−1 2 JτT Q2 Jτ 2 ek ∈∂Ω

≤ c w ·

τ ∈Ωbh

τ ∈Ωbh

1

Jτ−1 2 JτT Q2 Jτ 2

τ

1/2 .

τ

By combining (34) and (35), we complete the estimate of the right hand side of (32), and hence the proof of (28). The proof of (29) is similar to that of (18).

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

106

WEIMING CAO

By using Theorems 4.1 and 5.1, and the triangle inequality (36)

∇u − P∇uN ≤ ∇u − P∇uI + P(∇uN − ∇uI ) ≤ ∇u − P∇uI + ∇uN − ∇uI ,

we have immediately the following result on the superconvergence of the postprocessed ﬁnite element solution. Theorem 5.2. Under the same assumption of Theorem 5.1, we have

(1 + N −α Jτ−1 4 Jτ 2 ) · JτT Q2 Jτ 2 ∇u − P∇uN ≤ c τ

(37)

τ

1/2 + Jτ−1 4 Jτ 2 · JτT Q3 Jτ 3 1/2

+ c Jτ−1 2 · JτT Q2 Jτ 2 . τ ∈Ωbh

τ

In addition, if the partition TN is quasi-uniform under metric M , then

−1/2 ∇u − P∇uN ≤ c N ((CM )2 N −1 + CM N −α/2 F −1 4 F 2 ) Ω

(38)

· F T Q2 F 2 + (CM )2 N −1 F −1 4 F 2 · F T Q3 F 3

1/2 −1/2 CM + cN F −1 2 · F T Q2 F 2 .

1/2

Ωbh

Remark 5.1. In practical computation, the mesh metric M is typically chosen to be a certain scalar multiple of Q2 . In this case, the integrals in (37) involving elements next to the domain boundary are only an O(N −1/2 ) portion of the same integral over the entire Ω, and the second integral in (38) is of order O(N −1/2 ), which implies ∇u − P∇uN converges at the rate of O(N −(1+min(1,α))/2 ). It also implies the asymptotic exactness (as N → ∞) of the error estimator ∇uN − P∇uN for solutions satisfying condition ∇u − ∇uN ≥ cN −1/2 . Remark 5.2. Comparing the estimates for ∇uN − ∇uI in Theorem 4.1 and ∇u − P(∇uI ) in Theorem 5.1, we note that there is an extra factor cond(Jτ ) = Jτ Jτ−1 for the bound of the former. This factor measures the aspect ratio of τ . Thus for highly anisotropic meshes, ∇uN − ∇uI can be much larger than ∇u − P(∇uI ), even though both of them are of the same convergence rate. This is a new feature not existing in the superconvergence estimates on isotropic meshes. We believe this extra factor is not due to the technical treatment in our analysis, but rather the fact that ∇uN −∇uI is basically the H 1 -norm of the H 1 -projection of u − uI , while ∇u − P(∇uI ) is mainly the L2 -error of an L2 -projection. While L2 -projection is insensitive to the aspect ratio of elements, H 1 -projection is not. Numerical results in the next section seem to support this observation. 6. Numerical results We present in this section some numerical results for the linear ﬁnite element solution and its error estimates of the model problem. We choose in (6) the domain Ω = [0, 1]2 , and the coeﬃcients A = I2 , b = 0, d = 1. The right-hand side function f and the Dirichlet boundary condition are selected so that (6) admits the following exact solutions:

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

SUPERCONVERGENCE ON ANISOTROPIC MESHES

107

Example (1). u(x, y) = tanh[−50(r − 0.3)], where r is the distance to the center of Ω. Example (2). u(x, y) = [cosh(kx) + cosh(ky)] / [2 cosh(k)] with k = 50. The ﬁrst example involves an internal layer along a circle in the domain, while the second one a boundary layer along the top and right sides of the domain. To obtain the adaptive meshes and the ﬁnite element solutions, we start with a Delaunay triangulation of Ω, and calculate the ﬁnite element solution based on it. Then we recover the second order derivatives of the solution, and deﬁne a mesh metric on Ω. We use a two-dimensional anisotropic generator, bamg [13], to create an adaptive mesh that is quasi-uniform under the mesh metric, and calculate again the ﬁnite element solution based on the adaptive mesh. This process is repeated 20 times, which is hand-picked but seems suﬃcient to obtain a convergent mesh and the ﬁnite element solution. In practical computation, more judicious criterions, e.g., a threshold on the diﬀerences between the meshes or solutions in consecutive iterates, should be applied to control this iterative process. In order to deﬁned the mesh metrics, we ﬁrst approximate the second derivatives ∇2 u by (39)

˜ = ∇(P (∇uN )), H

where P is the L2 -projection onto SN . We choose here the L2 -projection of uN because it is the quantity involved in our analysis; see Theorem 5.2. Other types ˜ is piecewise of postprocessing, e.g., local averaging, are also possible. Since H constant, it is further approximated by a continuous piecewise linear function, still ˜ with the nodal values equal to the weighted average over the element denoted by H, patch around each node. The anisotropic behavior of ∇2 u is then characterized by ˜ + δ I, Q2 = abs(H) ˜ is the matrix of the same eigenvectors as H, ˜ but with eigenvalues where abs(H) ˜ the absolute value of H’s. δ ≥ 0 is a user deﬁned parameter to avoid the mesh metric based on it being singular in case ∇2 u becomes singular. A smaller δ makes Q2 a more precise characterization of the anisotropic behavior of ∇2 u and the mesh metric more eﬀective, while a larger δ results in a less aggressive anisotropic reﬁnement and is easier for reﬁnement. An optimal value of δ is problem dependent and should be selected dynamically. However, this is not our main objective in this paper. We choose for simplicity in all our numerical examples δ = 1 (a small value compared to the magnitude of the second derivatives of u). The mesh metric M is then deﬁned as in [5, 7] for the minimization of the H 1 error for the linear interpolation of u on anisotropic meshes, (40)

M = c [λmax /λmin ]1/4 Q2 ,

where λmax and λmin are the largest and smallest eigenvalues of Q2 , and c is a constant to control the total number of elements. In order to check the O(N −(1+α)/2 )-approximate parallelogram property for the adaptive meshes used in our computation, we introduce a quantity Map to measure how close a pair of adjacent elements is to a parallelogram. Suppose τ1 and τ2 share a common edge e. J1 and J2 are the Jacobian of the aﬃne mapping from an equilateral standard element τˆ to τ1 and τ2 , respectively, with the two vertices

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

108

WEIMING CAO

opposite to e being the images of the same vertex of τˆ. Deﬁne 1 ( I + J1−1 J2 + I + J2−1 J1 ), Map (e) = (41) 2 where the matrix norm is chosen to be Frobenius norm for easy calculation. Note that when τ1 and τ2 form an O(N −(1+α)/2 )-approximate parallelogram, both terms on the right-hand side of (41) are of magnitude O(N −α/2 ) (see Section 2), thus Map = O(N −α/2 ). We use their average in (41) to produce a measurement symmetric for both sides of an edge. As noted in Lemma 2.2, this measure is invariant under any aﬃne mapping of the element pair. When τ1 and τ2 form an exact parallelogram, Map (ek ) = 0. When they are congruent but form a triangle, √ 2 3 Map = 3 ≈ 1.155. If τ1 is shape regular but τ2 is of a high aspect ratio, then Map >> 1. ˜ (with some We notice that Bank and Xu established in [4] the convergence of H 2 additional smoothing steps to improve the convergence rate) to ∇ u for quasiuniform meshes. In the case of adaptive anisotropic meshes, there is no such theory ˜ for mesh reﬁnement has been in practice for yet, even though the idea of using H ˜ is a good approximation of ∇2 u, the a long time; see, e.g., [8, 10–12]. When H above deﬁned mesh metric will be optimal for minimizing the H 1 semi-norm of the interpolation error, and thus for minimizing the a priori bound of the ﬁnite element solution error as stated in Theorem 3.1. We examine the following quantities involved at various stages of our analysis: EuN EpuI EN I

= ∇u − ∇uN , = ∇u − P (∇uI ), = ∇uN − ∇uI ,

EN pN EpuN EpN I

= ∇uN − P (∇uN ), = ∇u − P (∇uN ), = P (∇uN − ∇uI ).

In addition, we check the eﬀectivity factor η of the error estimator EN pN : η = EN pN / EuN . For each of the above errors E and √ |η − 1|, we estimate its rate of convergence β by linear regression of ln E vs. ln(1/ N ), i.e., E ≈ cN −β/2 . For comparison, we also report the results for solutions based on uniform meshes obtained from bisecting uniform square meshes on Ω with the 135◦ diagonals. Note that these uniform meshes can be considered as quasi-uniform under mesh metric I and each pair of the adjacent elements forms an O(N −(1+α)/2 )-approximate parallelogram with α = ∞. In this case, EN I will be of order N −1 according to Theorem 4.1, and EpuI and EpuN will also be of order N −1 in the absence of boundary errors according to Theorems 5.1 and 5.2. We list in Tables 1 and 2 the results for Example (1) with adaptive and uniform meshes, respectively, and Tables 3 and 4 for Example (2) similarly. Clearly, in all the cases, EuN converges linearly, i.e., at the rate of O(N −1/2 ), and the eﬀectivity factor η → 1, which indicates that the error estimator EN pN is asymptotically exact. In addition, EN I converges quadratically on uniform meshes as predicted in Theorem 4.1 with α = ∞. On adaptive meshes, the superconvergence of EN I is much weaker, with α ≈ 0.12 only. This result is similar to that reported in [3] for the case of adaptive meshes composed of shape regular elements. The slower convergence rate with adaptive meshes results mainly from the fact that there are elements in the triangulation that do not form approximate parallelograms due to

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

SUPERCONVERGENCE ON ANISOTROPIC MESHES

109

Table 1. Example (1) based on adaptive meshes generated from bamg with mesh metrics deﬁned by (40). N 4070 8024 16082 32114 64112 128762 257114 Order β

EuN 6.62e-1 4.45e-1 3.01e-1 2.10e-1 1.46e-1 1.04e-1 7.40e-2 1.053

EN pN 6.57e-1 4.41e-1 2.99e-1 2.09e-1 1.45e-1 1.04e-1 7.36e-2 1.051

EpuI 1.95e-1 1.09e-1 6.18e-2 3.64e-2 2.26e-2 1.41e-2 8.82e-3 1.485

EpuN 1.90e-1 1.10e-1 6.32e-2 3.82e-2 2.37e-2 1.53e-2 9.81e-3 1.425

EN I 1.36e-1 9.11e-2 6.16e-2 4.07e-2 2.78e-2 1.99e-2 1.30e-2 1.126

EpN I 8.06e-2 5.36e-2 3.63e-2 2.39e-2 1.62e-2 1.14e-2 7.55e-3 1.137

|η − 1| 7.81e-3 8.74e-3 7.78e-3 7.09e-3 6.78e-3 6.65e-3 6.11e-3 0.147

Table 2. Example (1) based on uniform meshes. N 4050 7938 16200 32258 64082 128018 256328 Order β

EuN 3.33e+0 2.45e+0 1.76e+0 1.26e+0 9.01e-1 6.39e-1 4.52e-1 0.965

EN pN EpuI EpuN EN I 3.01e+0 2.04e+0 2.22e+0 6.33e-1 2.35e+0 1.09e+0 1.26e+0 3.64e-1 1.72e+0 5.28e-1 6.44e-1 1.91e-1 1.25e+0 2.61e-1 3.29e-1 9.90e-2 8.95e-1 1.30e-1 1.67e-1 5.07e-2 6.37e-1 6.46e-2 8.38e-2 2.56e-2 4.52e-1 3.22e-2 4.20e-2 1.29e-2 0.924 2.013 1.927 1.891

EpN I 3.74e-1 2.63e-1 1.61e-1 9.11e-2 4.87e-2 2.51e-2 1.27e-2 1.655

|η − 1| 9.63e-2 4.35e-2 2.23e-2 1.15e-2 5.92e-3 3.00e-3 1.51e-3 1.976

Table 3. Example (2) based on adaptive meshes generated from bamg with mesh metrics deﬁned by (40). N 4030 8086 15926 32064 63867 129809 256928 Order β

EuN 8.06e-2 5.45e-2 3.79e-2 2.67e-2 1.89e-2 1.32e-2 9.35e-3 1.031

EN pN 7.88e-2 5.37e-2 3.74e-2 2.65e-2 1.88e-2 1.31e-2 9.31e-3 1.023

EpuI 1.72e-2 9.25e-3 5.51e-3 3.21e-3 2.01e-3 1.21e-3 7.67e-4 1.485

EpuN 1.75e-2 9.52e-3 5.92e-3 3.39e-3 2.20e-3 1.32e-3 8.70e-4 1.438

EN I 1.21e-2 7.16e-3 5.38e-3 3.49e-3 2.27e-3 1.56e-3 1.16e-3 1.129

EpN I 7.39e-3 4.46e-3 3.26e-3 2.00e-3 1.34e-3 8.78e-4 6.51e-4 1.178

|η − 1| 2.28e-2 1.48e-2 1.20e-2 7.93e-3 6.74e-3 5.01e-3 4.32e-3 0.797

geometric constraints, or the approximate parallelograms are rather poorly shaped; see Figure 3 and Figure 6 for a close-up look of typical meshes. Strictly speaking, these meshes belong to a more general type satisfying the so-called condition (α, σ) introduced in [3, 30] for shape regular meshes; namely, most element pairs in the partition form O(N −(1+α)/2 )-approximate parallelograms and the measure of those that do not is at most O(N −σ/2 ). The extension of our analysis to this type of meshes is straightforward. We choose not to include it here for the purpose of

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

110

WEIMING CAO

Table 4. Example (2) based on uniform meshes. N 4050 7938 16200 32258 64082 128018 256328 Order β

EuN 1.07e+00 7.86e-01 5.58e-01 3.99e-01 2.84e-01 2.01e-01 1.42e-01 0.975

EN pN 7.60e-01 6.16e-01 4.72e-01 3.54e-01 2.61e-01 1.90e-01 1.37e-01 0.835

EpuI 6.91e-01 4.54e-01 2.83e-01 1.76e-01 1.08e-01 6.53e-02 3.93e-02 1.387

EpuN 6.91e-01 4.54e-01 2.83e-01 1.76e-01 1.08e-01 6.53e-02 3.93e-02 1.387

EN I 1.37e-04 6.99e-05 3.42e-05 1.72e-05 8.66e-06 4.33e-06 2.16e-06 2.000

EpN I 1.33e-04 6.86e-05 3.39e-05 1.71e-05 8.63e-06 4.33e-06 2.16e-06 1.986

|η − 1| 2.90e-01 2.16e-01 1.55e-01 1.11e-01 7.96e-02 5.66e-02 4.01e-02 0.957

simplifying our presentation. Instead, we report some statistics about the measurement Map for how close each pair of elements is to a parallelogram. We display in Figure 4 and Figure 7 the distribution of Map for a typical adaptive mesh generated by bamg in Examples (1) and (2), respectively. To further quantify this distribution, we plot in Figure 5(a) and Figure 8(a) the percentile of the element pairs whose Map are below a given value for each of the adaptive meshes reported in Table 1 and Table 3. It is clear from these two ﬁgures that as the mesh gets reﬁned by bamg, Map decreases in general, and more and more element pairs become closer and closer to exact parallelograms. We plot in Figure 5(b) and Figure 8(b) the average value, M ap , of Map over the entire mesh for each of the adaptive meshes. ¯ ) with Based on linear regression of ln M ap vs. ln N , we ﬁnd that M ap ≈ O(N −α/2 α ¯ = 0.208 and 0.278 for the adaptive meshes generated by bamg in Example (1) and Example (2), respectively. It is also noted that EpuI exhibits superconvergence for both examples and both types of meshes as described in Theorem 5.1. However, with uniform meshes, the convergence rate is about O(N −1 ) for Example (1) and O(N −1.387/2 ) for Example (2). We believe the loss of the convergence rate in Example (2) is due to the boundary term in (28). Note that the solution of Example (1) is nearly 0 around the boundary of Ω, while it it not the case for Example (2). With adaptive meshes, EpuI converges at a rate O(N −1.485/2 ) for both examples, and the eﬀect of boundary terms for Example (2) seems to be masked by that of non-ideal approximate parallelograms. By comparing EpuI to EN I , we note that on uniform meshes EN I is always smaller than EpuI , while on adaptive meshes it is the other way around. This phenomenon reﬂects their dependence on the element aspect ratio stated in Remark 5.2; namely, higher element aspect ratios have a stronger negative impact on the accuracy of EN I than on EpuI . On the other hand, for both examples on adaptive meshes, EpuI seems to behave very similarly to EpuN , even though the bound for the latter involves an extra factor of the element aspect ratio cond(Fτ ); see Theorems 5.1 and 5.2. It raises the question of whether this extra factor in the bound for EpuN is indeed necessary. In our analysis this extra factor is introduced due to the use of triangle inequality (36), which relies on the estimate (17) of EN I . Unfortunately, we still do not know if it is possible not to use EN I to provide EpuN with a similar bound as for EpuI .

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

SUPERCONVERGENCE ON ANISOTROPIC MESHES 1

111

0.45

0.9 0.4

0.8 0.7

0.35

0.6 0.5

0.3

0.4 0.25

0.3 0.2

0.2

0.1 0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.15 0.15

0.2

0.25

0.3

0.35

0.4

Figure 3. Example (1): Left: An adaptive mesh generated by bamg based on metric (40) with the total number of elements N = 4070. Right: Its close-up look.

Figure 4. Example (1): Left: The distribution of measure Map deﬁned in (41) for the closeness of a pair of adjacent elements to a parallelogram, on the adaptive mesh with N = 4070. The color around an edge e represents the magnitude of Map (e) for the pair sharing e. Right: Its close-up look.

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

0.45

112

WEIMING CAO

90

0.36

Average value of Map in each mesh

0.38

80 70

Percentage of element pairs whose M

ap

are smaller than μ

100

60 50 40 30

0.34 0.32 0. 3 0.28 0.26

20

0.24 10 0

0

0.5

1

0.22 3 10 10 4 10 5 10 6 Number of elements for meshes generated by BAMG

1.5

μ

Figure 5. Example (1): Left: Percentile curves for the number of element pairs whose Map are below a given value speciﬁed by the horizontal axis. The 7 curves from the bottom up represent the 7 adaptive meshes reported in Table 1 in increasing order of N . Right: The average value of Map vs. the number N of elements for each of the adaptive meshes.

1

1

0.9

0.98

0.8

0.96

0.7

0.94

0.6

0.92

0.5

0.9

0.4

0.88

0.3

0.86

0.2

0.84

0.1

0.82

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.8 0.8

0.82

0.84

0.86

0.88

0.9

0.92

0.94

0.96

Figure 6. Example (2): Left: An adaptive mesh generated by bamg based on metric (40) with the total number of elements N = 4030. Right: Its close-up look.

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

0.98

1

SUPERCONVERGENCE ON ANISOTROPIC MESHES

113

0.26

100 90

0.24 Average value of M ap in each mesh

Percentage of element pairs whose M ap are smaller than

Figure 7. Example (2): Left: The distribution of measure Map deﬁned in (41) for the closeness of a pair of adjacent elements to a parallelogram, on the adaptive mesh with N = 4050. The color around an edge e represents the magnitude of Map (e) for the pair sharing e. Right: Its close-up look.

80 70 60 50 40 30 20

0.22 0. 2 0.18 0.16 0.14

10 0

0

0.5

1

1.5

0.12 10 3

10 4

10 5

Number of elements for meshes generated by BAMG

Figure 8. Example (2): Left: Percentile curves for the number of element pairs whose Map are below a given value speciﬁed by the horizontal axis. The 7 curves from the bottom up represent the 7 adaptive meshes reported in Table 3 in increasing order of N . Right: The average value of Map vs. the number N of elements for each of the adaptive meshes.

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

10 6

114

WEIMING CAO

7. Conclusion and discussions For the linear ﬁnite element method based on adaptive anisotropic meshes in two dimensions, we established the superconvergence of the ﬁnite element solutions to the interpolation of exact solutions in energy norm. We also proved the superconvergence of the recovery process using the global L2 -projection of the gradients of the ﬁnite element solutions. Our analysis is based on the notion that the mesh is quasi-uniform under a Riemannian metric and the notion that each adjacent element pair forms an approximate parallelogram. These notions are extensions of the corresponding concepts in classical superconvergence studies on meshes composed of shape regular elements. In particular, our analysis follows the same methodology developed in Bank and Xu [3] for the case of quasi-uniform meshes, and can be viewed as an extension of their conclusion to the case of unstructured anisotropic meshes. Numerical examples involving both internal and boundary layers are in support of these superconvergence properties. It should be noted that our analysis is concerned with the asymptotic behavior of the solution and its post-processing. Therefore, the conclusion is meaningful only when the mesh is suﬃciently ﬁne for the underlying problem. In the initial stages of an adaptive ﬁnite element solution where mesh resolution is far from enough, the error from data approximation can be more pronounced and the error estimates based on post-processing may not be accurate. In this case, except the anisotropic behavior of the PDEs is known a-priori, one should be extremely careful about the use of anisotropic coarse meshes, since a highly anisotropic element may introduce a much larger error than an isotropic element does if used inappropriately. Therefore, it is generally advised to begin an adaptive process with isotropic coarse meshes and other types of error estimators should be considered, too. Several issues remain unresolved in this paper. First, it is observed in numerical experiments that the convergence of ∇u − P(∇uN ) is similar to ∇u − P(∇uI ) and both are less sensitive to high element aspect ratio than ∇uN −∇uI . However, the estimate for ∇u − P(∇uN ) involves an extra factor of the element aspect ratio than the one for ∇u − P(∇uI ). It seems a sharper analysis tool is needed to remove this extra factor if it is indeed non-essential. Second, with unstructured adaptive meshes the superconvergence for both ∇uN − ∇uI and ∇u − P(∇uN ) is rather weak due to the non-uniformity of the elements. Bank and Xu proposed in [4] to reduce the dependence on mesh uniformity and strengthen the superconvergence of the gradient recovery by employing several cycles of multigrid type smoothing in addition to the L2 -projection. Their numerical results showed such a treatment may improve the superconvergence substantially. This idea is yet to be tested for the case of anisotropically reﬁned meshes. Furthermore, we hope that with the improved superconvergence property for the ˜ in (39) can recovered gradients the convergence for the approximation of ∇2 u by H be established in the case of anisotropic meshes. There are several other types of gradient recovery techniques used in practice. The most popular one is probably Zienkiewicz-Zhu’s superconvergence patch recovery (SPR) based on local averaging [36, 37], and a more recent variation, the polynomial preserving recovery (PPR), by Zhang [23, 32, 34] based on higher order local approximation of the solution. There has been extensive work on the analysis of these techniques on uniform meshes, structured and unstructured isotropic

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

SUPERCONVERGENCE ON ANISOTROPIC MESHES

115

meshes [15, 30], and recently on structured anisotropic meshes [32]. With the establishment of the superconvergence of ∇uN −∇uI , it is now possible to establish the superconvergence of these local gradient recovery techniques on general adaptively reﬁned anisotropic meshes. We will report a detailed analysis and numerical results for these techniques in a separate paper.

References [1] Mark Ainsworth and J. Tinsley Oden, A Posteriori Error Estimation in Finite Element Analysis, Pure and Applied Mathematics (New York), Wiley-Interscience [John Wiley & Sons], New York, 2000. MR1885308 (2003b:65001) [2] Ivo Babuˇska and Theofanis Strouboulis, The Finite Element Method and its Reliability, Numerical Mathematics and Scientiﬁc Computation, The Clarendon Press Oxford University Press, New York, 2001. MR1857191 (2002k:65001) [3] Randolph E. Bank and Jinchao Xu, Asymptotically exact a posteriori error estimators. I. Grids with superconvergence, SIAM J. Numer. Anal. 41 (2003), no. 6, 2294–2312 (electronic), DOI 10.1137/S003614290139874X. MR2034616 (2004k:65194) [4] Randolph E. Bank and Jinchao Xu, Asymptotically exact a posteriori error estimators. II. General unstructured grids, SIAM J. Numer. Anal. 41 (2003), no. 6, 2313–2332 (electronic), DOI 10.1137/S0036142901398751. MR2034617 (2004m:65212) [5] Weiming Cao, On the error of linear interpolation and the orientation, aspect ratio, and internal angles of a triangle, SIAM J. Numer. Anal. 43 (2005), no. 1, 19–40 (electronic), DOI 10.1137/S0036142903433492. MR2177954 (2006k:65023) [6] Weiming Cao, An interpolation error estimate in R2 based on the anisotropic measures of higher order derivatives, Math. Comp. 77 (2008), no. 261, 265–286 (electronic), DOI 10.1090/S0025-5718-07-01981-3. MR2353953 (2008j:65218) [7] Weiming Cao, An interpolation error estimate on anisotropic meshes in Rn and optimal metrics for mesh reﬁnement, SIAM J. Numer. Anal. 45 (2007), no. 6, 2368–2391 (electronic), DOI 10.1137/060667992. MR2361894 (2008k:65264) [8] Long Chen, Pengtao Sun, and Jinchao Xu, Optimal anisotropic meshes for minimizing interpolation errors in Lp -norm, Math. Comp. 76 (2007), no. 257, 179–204, DOI 10.1090/S00255718-06-01896-5. MR2261017 (2008e:65385) [9] P. G. Ciarlet, The Finite Element Methods for Elliptic Problems, SIAM Classics in Applied Mathematics, SIAM, Philadelphia, 2002. [10] L. Formaggia and S. Perotto, New anisotropic a priori error estimates, Numer. Math. 89 (2001), no. 4, 641–667, DOI 10.1007/s002110100273. MR1865506 (2002j:65110) [11] L. Formaggia and S. Perotto, Anisotropic error estimates for elliptic problems, Numer. Math. 94 (2003), no. 1, 67–92, DOI 10.1007/s00211-002-0415-z. MR1971213 (2004c:65127) [12] Wagdi G. Habashi, Julien Dompierre, Yves Bourgault, Djaﬀar Ait-Ali-Yahia, Michel Fortin, and Marie-Gabrielle Vallet, Anisotropic mesh adaptation: towards user-independent, meshindependent and solver-independent CFD. I. General principles, Internat. J. Numer. Methods Fluids 32 (2000), no. 6, 725–744, DOI 10.1002/(SICI)1097-0363(20000330)32:6725::AIDFLD9353.0.CO;2-4. MR1911901 (2003d:76106) [13] F. Hecht, Bidimensional anisotropic mesh generator (User’s Manual), INRIA, Rocquencourt, 1997. http://www-rocqq.inria.fr/gamma/cdrom/www/bamg/eng.htm [14] Gerd Kunert and Serge Nicaise, Zienkiewicz-Zhu error estimators on anisotropic tetrahedral and triangular ﬁnite element meshes, M2AN Math. Model. Numer. Anal. 37 (2003), no. 6, 1013–1043, DOI 10.1051/m2an:2003065. MR2026406 (2004k:65196) [15] A. M. Lakhany, I. Marek, and J. R. Whiteman, Superconvergence results on mildly structured triangulations, Comput. Methods Appl. Mech. Engrg. 189 (2000), no. 1, 1–75, DOI 10.1016/S0045-7825(99)00281-9. MR1779678 (2001f:74041) [16] Jichun Li and Mary F. Wheeler, Uniform convergence and superconvergence of mixed ﬁnite element methods on anisotropically reﬁned grids, SIAM J. Numer. Anal. 38 (2000), no. 3, 770–798, DOI 10.1137/S0036142999351212. MR1781203 (2001f:65137)

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

116

WEIMING CAO

[17] Qun Lin and Jia-Fu Lin, Superconvergence on anisotropic meshes, Int. J. Inf. Syst. Sci. 3 (2007), no. 2, 261–266. MR2297053 (2008c:65337) [18] Shipeng Mao, Shaochun Chen, and Dongyang Shi, Convergence and superconvergence of a nonconforming ﬁnite element on anisotropic meshes, Int. J. Numer. Anal. Model. 4 (2007), no. 1, 16–38. MR2289730 (2007k:65188) [19] F. Marcuzzi, M. Morandi Cecchi, and M. Venturin, An anisotropic unstructured triangular adaptive mesh algorithm based on error and error gradient information, Math. Comput. Simulation 78 (2008), no. 5-6, 645–652, DOI 10.1016/j.matcom.2008.04.006. MR2435592 (2009e:65197) [20] Stefano Micheletti and Simona Perotto, Reliability and eﬃciency of an anisotropic Zienkiewicz-Zhu error estimator, Comput. Methods Appl. Mech. Engrg. 195 (2006), no. 9-12, 799–835, DOI 10.1016/j.cma.2005.02.009. MR2195291 (2006i:65179) [21] Jean-Marie Mirebeau, Optimal meshes for ﬁnite elements of arbitrary order, Constr. Approx. 32 (2010), no. 2, 339–383, DOI 10.1007/s00365-010-9090-y. MR2677884 (2011g:65278) [22] J. -M. Mirebeau, Approximation adaptative et anisotrope par ´ el´ eements ﬁnis: Th´ eorie et Algorithmes, Ph.D. dissertation, Laboratoire Jacques-Louis Lions, Universit´ e Pierre et Marie Curie, 2010. [23] Ahmed Naga and Zhimin Zhang, A posteriori error estimates based on the polynomial preserving recovery, SIAM J. Numer. Anal. 42 (2004), no. 4, 1780–1800 (electronic), DOI 10.1137/S0036142903413002. MR2114301 (2005k:65231) [24] M. Picasso, Numerical study of the eﬀectivity index for an anisotropic error indicator based on Zienkiewicz-Zhu error estimator, Comm. Numer. Methods Engrg. 19 (2003), no. 1, 13–23, DOI 10.1002/cnm.546. MR1952014 (2004c:65131) [25] M. Picasso, An anisotropic error indicator based on Zienkiewicz-Zhu error estimator: application to elliptic and parabolic problems, SIAM J. Sci. Comput. 24 (2003), no. 4, 1328–1355 (electronic), DOI 10.1137/S1064827501398578. MR1976219 (2004e:65124) [26] R. Rannacher and S. Turek, Simple nonconforming quadrilateral Stokes element, Numer. Methods Partial Diﬀerential Equations 8 (1992), no. 2, 97–111, DOI 10.1002/num.1690080202. MR1148797 (92i:65170) [27] Dong-yang Shi, Shi-peng Mao, and Shao-chun Chen, An anisotropic nonconforming ﬁnite element with some superconvergence results, J. Comput. Math. 23 (2005), no. 3, 261–274. MR2133919 (2006h:65194) [28] Lars B. Wahlbin, Superconvergence in Galerkin Finite Element Methods, Lecture Notes in Mathematics, vol. 1605, Springer-Verlag, Berlin, 1995. MR1439050 (98j:65083) [29] Haijun Wu and Zhimin Zhang, Can we have superconvergent gradient recovery under adaptive meshes?, SIAM J. Numer. Anal. 45 (2007), no. 4, 1701–1722, DOI 10.1137/060661430. MR2338406 (2008k:65252) [30] Jinchao Xu and Zhimin Zhang, Analysis of recovery type a posteriori error estimators for mildly structured grids, Math. Comp. 73 (2004), no. 247, 1139–1152 (electronic), DOI 10.1090/S0025-5718-03-01600-4. MR2047081 (2005f:65141) [31] Ningning Yan and Aihui Zhou, Gradient recovery type a posteriori error estimates for ﬁnite element approximations on irregular meshes, Comput. Methods Appl. Mech. Engrg. 190 (2001), no. 32-33, 4289–4299, DOI 10.1016/S0045-7825(00)00319-4. MR1832657 (2002c:65189) [32] Zhimin Zhang, Polynomial preserving recovery for meshes from Delaunay triangulation or with high aspect ratio, Numer. Methods Partial Diﬀerential Equations 24 (2008), no. 3, 960– 971, DOI 10.1002/num.20300. MR2402584 (2009c:65327) [33] Zhimin Zhang, Polynomial preserving gradient recovery and a posteriori estimate for bilinear element on irregular quadrilaterals, Int. J. Numer. Anal. Model. 1 (2004), no. 1, 1–24. MR2052728 (2005b:65119) [34] Zhimin Zhang and Ahmed Naga, A new ﬁnite element gradient recovery method: superconvergence property, SIAM J. Sci. Comput. 26 (2005), no. 4, 1192–1213 (electronic), DOI 10.1137/S1064827503402837. MR2143481 (2006d:65137) [35] Q.D. Zhu, Theory on High Accuracy Post-processing in the Finite Element Method, Science Press, Beijing, 2008 (in Chinese).

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

SUPERCONVERGENCE ON ANISOTROPIC MESHES

117

[36] O. C. Zienkiewicz and J. Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. I. The recovery technique, Internat. J. Numer. Methods Engrg. 33 (1992), no. 7, 1331–1364, DOI 10.1002/nme.1620330702. MR1161557 (93c:73098) [37] O. C. Zienkiewicz and J. Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. II. Error estimates and adaptivity, Internat. J. Numer. Methods Engrg. 33 (1992), no. 7, 1365–1382, DOI 10.1002/nme.1620330703. MR1161558 (93c:73099) Department of Mathematics, University of Texas at San Antonio, San Antonio, Texas 78249 E-mail address: [email protected]

Licensed to Univ of Texas at San Antonio. Prepared on Mon Jan 12 00:54:52 EST 2015 for download from IP 129.115.103.99. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use