ON THE ERROR OF LINEAR INTERPOLATION ... - Semantic Scholar
c 2005 Society for Industrial and Applied Mathematics
SIAM J. NUMER. ANAL. Vol. 43, No. 1, pp. 19–40
ON THE ERROR OF LINEAR INTERPOLATION AND THE ORIENTATION, ASPECT RATIO, AND INTERNAL ANGLES OF A TRIANGLE∗ WEIMING CAO† Abstract. In this paper, we attempt to reveal the precise relation between the error of linear interpolation on a general triangle and the geometric characters of the triangle. Taking the model problem of interpolating quadratic functions, we derive two exact formulas for the H 1 -seminorm and L2 -norm of the interpolation error in terms of the area, aspect ratio, orientation, and internal angles of the triangle. These formulas indicate that (1) for highly anisotropic triangular meshes the H 1 -seminorm of the interpolation error is almost a monotonically decreasing function of the angle between the orientations of the triangle and the function; (2) maximum angle condition is not |λ1 /λ2 | or essential if the mesh is aligned with the function and the aspect ratio is of magnitude less, where λ1 and λ2 are the eigenvalues of the Hessian matrix of the function. With these formulas we identify the optimal triangles, which produce the smallest H 1 -seminorm of the interpolation error, 1 to be the acute isosceles aligned with the solution and of an aspect ratio about 0.8| λ |. The L2 λ2 norm of the interpolation error depends on the orientation and the aspect ratio of the triangle, but 2 not directly on its maximum or minimum angles. The optimal triangles for the L -norm are those aligned with the solution and of an aspect ratio |λ1 /λ2 |. These formulas can be used to formulate more accurate mesh quality measures and to derive tighter error bounds for interpolations. Key words. anisotropic mesh, linear interpolation, aspect ratio, mesh alignment, maximum angle condition AMS subject classifications. 65D05, 65L50, 65N15, 65N50 DOI. 10.1137/S0036142903433492
1. Introduction. It is well known that on quasi-uniform meshes the accuracy of piecewise linear interpolation is first order in the H 1 -norm. More precisely, denote by H m (Ω) the usual Sobolev spaces of order m on a bounded domain Ω ∈ R2 . For any u ∈ H 2 (Ω), the H 1 -norm of the error between u and its piecewise linear interpolation uI can be bounded by u − uI H 1 (Ω) ≤ ch|u|H 2 (Ω) , where h is the diameter of the triangles in the mesh, and c is a constant independent of h and u. If the mesh is not quasi-uniform, then the error contributed from each triangle K can be bounded by u − uI H 1 (K) ≤ c
h2 |u|H 2 (K) , ρ
where ρ is the diameter of the largest inscribed circle in K. This error bound guarantees the H 1 -norm of the error converge to 0 as h → 0, as long as hρ remains bounded. This is equivalent to requiring that the minimum angle of K is bounded from 0 [5]. ∗ Received
by the editors August 19, 2003; accepted for publication (in revised form) May 7, 2004; published electronically April 26, 2005. This research was supported in part by the NSF under grant DMS-0209313. http://www.siam.org/journals/sinum/43-1/43349.html † Department of Applied Mathematics, University of Texas at San Antonio, San Antonio, TX 78249 ([email protected]). 19
20
WEIMING CAO
However, this error estimate is not tight when the triangle has only one small angle. Babu˘ska and Aziz [2] showed that the minimum angle condition is actually not essential for the convergence of linear interpolation. They improved the error bound as u − uI H 1 (K) ≤ Γ(α)h|u|H 2 (K) , where Γ(α) is an increasing function of the maximum angle α of K. Therefore, in order to guarantee the convergence it is only required that the maximum angle be bounded from π. This is the well-known maximum angle condition. This condition is necessary and sufficient for the convergence of the linear interpolation process over the class of functions of H 2 (K). In adaptive computation, one often knows a priori or a posteriori some information of the functions to be approximated. This information can be used to restrict the set of functions considered and to avoid the divergence case although the triangle has a maximum angle close to π. Indeed, long and thin triangles violating the maximum angle condition have been used successfully in engineering, particularly in computational fluid dynamics [1, 15]. For problems with very different length scales in different spatial directions, long and thin triangles turn out to be better choices than shape regular ones if they are properly used. This motivated an intensive study on the error analysis for anisotropic meshes in the finite element method (FEM). For instance, Apel [1] described an error estimate in terms of the length scales h1 and h2 along the x and y directions, respectively, |u − uI |W m,p (K) ≤ c
α1 +α2 =−m
1 α2 hα 1 h2 |
∂ −m u |W m,p (K) , ∂xα1 ∂y α2
m = 0, 1,
where W m,p (K) is the Sobolev space of functions whose up to mth order derivatives are Lp -integratable. If u and K are aligned and the maximum angle condition is satisfied, this estimate is asymptotically accurate, i.e., both sides of the above inequality have the same hi order. But when the maximum angle of K approaches π, the ratio of the error bound over the actual error norm goes to infinity [7]. Formaggia and Perotto [7] presented another type of estimate of the H 1 -norm for the interpolation error based on the eigendecomposition of the affine mapping from a standard element to K. Their estimate is accurate when u and K are aligned and the maximum angle of K is close to π. However, the ratio of their estimate over the actual error norm goes to infinity when two angles of K approach π2 . More recent error analyses have been based on the overall mesh properties and the behavior of the approximation functions. For instance, Berzins [4] developed a mesh quality indicator measuring the correlation between the anisotropic features of the mesh and those of the solutions. Kunert [11] proposed a so-called matching function to quantify how good overall a mesh is for a specific solution. Huang [8] introduced measures for three aspects of the mesh qualities—aspect ratio, alignment, and adaptation—and an overall quality mesh measure based on them. He formulated the error bounds in terms of these measures and proposed a variational formulation to optimize the overall mesh quality measure to control the interpolation error. A similar idea was used by Huang and Sun [9] in formulating the monitor function in variational mesh adaptation. Needless to say, the interpolation error depends on the solution and the size and shape of the elements in the mesh. Understanding this relation is crucial for the
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
21
generation of efficient and effective meshes for the FEM. However, in all the error estimates for anisotropic meshes, the relation between the error and the geometric characters of a triangle, such as the alignment, aspect ratio, and internal angles, has not been revealed explicitly. In the mesh generation community, this relation is studied more closely for the model problem of interpolating quadratic functions. This model is a reasonable simplification of the cases involving general functions, since quadratic functions are the leading terms in the local expansion of the linear interpolation errors. For instance, Nadler [12] derived an exact expression for the L2 -norm of the linear interpolation error in terms of the three sides 1 , 2 , and 3 of the triangle K, u − uI 2L2 (K) =
|K| [(d1 + d2 + d3 )2 + d1 d2 + d2 d3 + d1 d3 ], 180
where |K| is the area of the triangle, di = i · Hi with H being the Hessian matrix of u. Bank and Smith [3] gave a similar formula for the H 1 -seminorm of the linear interpolation error; see (8) in section 3. Assuming u = λ1 x2 + λ2 y 2 , D’Azevedo and Simpson [6] derived the exact formula for the maximum norm of the interpolation error u − uI L2 ∞ (K) =
D12 D23 D31 , 16λ1 λ2 |K|2
where Dij = i · diag(λ1 , λ2 )i . Based on the geometric interpretation of this formula, they proved that for a fixed area the optimal triangle, which produces the smallest maximum interpolation error, √ √ is the one obtained by compressing an equilateral triangle by factors λ1 and λ2 along the two eigenvectors of the Hessian matrix of u. Furthermore, the optimal incidence for a given set of interpolation points √ √ is the Delaunay triangulation based on the stretching map (by factors λ1 and λ2 along the two eigenvector directions) of the grid points. Rippa [13] showed that the mesh obtained this way is also optimal for the Lp -norm of the error for any 1 ≤ p ≤ ∞. Thought these formulas are exact, they do not describe explicitly the relation between the error and the geometric characters of the triangle. In this paper, we attempt to reveal this relation explicitly and precisely. Taking the model problem of linear interpolation of quadratic functions, we derive two exact expressions for the H 1 -seminorm and L2 -norm of the interpolation error in terms of the area, aspect ratio, alignment direction, and internal angles of the triangle. From these formulas the effects of the geometric characters of a triangle can be clearly identified. They indicate that (1) for highly anisotropic triangular meshes the H 1 -seminorm of the interpolation error is almost a monotonically decreasing function of the angle between the orientation of the triangle and the orientation of the function; (2) maximum angle condition is critical if the mesh is not aligned with the function or if the aspect ratio is larger than |λ1 /λ 2 |; (3) if the triangles are aligned with the function and the aspect ratio is about |λ1 /λ2 | or less, then the error is not sensitive to the maximum angle of the triangle. Also, we can easily identify the best triangles which produce the smallest H 1 - or L2 -norm of the interpolation error. It turns out that in the sense of the H 1 -seminorm, the optimal triangle with a given area is the acute isosceles aligned with the function and of the aspect ratio about 0.8| λλ12 |. The L2 -norm of the interpolation error depends on the alignment and the aspect ratio of the triangle, but not directly on its maximum or minimum angles. The optimal triangles for the L2 -norm are those aligned with the solution and of an aspect ratio |λ1 /λ2 |.
22
WEIMING CAO
The organization of this paper is as follows. In section 2 we give a precise definition of the orientation and the aspect ratio of a triangle by using the singular value decomposition (SVD) of the affine mapping from a standard element to the triangle. In section 3 we derive the exact formulas for the H 1 -seminorm and L2 -norm of the interpolation error for the quadratic functions. Then we identify in section 4 a number of special cases that are of interests to the optimal design of meshes. Optimal choices of the aspect ratio and the alignment direction are discussed here. In section 5 we elaborate in more detail the effects of various geometric characters of a triangle on the interpolation error. Finally we demonstrate by an example the accuracy of the linear interpolation errors on different meshes. 2. Mapping and the geometric characters of a triangle. Let K be a triangle in the xy-plane, and let x1 , x2 , and x3 be the vertices of K. Denote by i ˆ be the the vector of the side opposite to xi (in counterclockwise direction). Let K equilateral triangle with the vertices √ √ 3 0 − 23 2 , ξ3 = . ξ1 = , ξ2 = 1 1 −2 − 12 √ ˆ are ei = 3[cos(2(i − 1)π/3), sin(2(i − 1)π/3)]T , i = 1, 2, 3. See The three sides of K Figure 1. η
ξ
1
y x
1
1
ξ2
ξ3
ξ
x
3
x
2
x
ˆ and physical element K. Fig. 1. Standard element K
Let xc = 31 (x1 + x2 + x3 ) be the center of K. The affine mapping (which maps ˆ to K can be expressed as x = M ξ + xc with ξ i to xi ) from K 1 M = √ (x1 − x3 ), x2 − xc . 3 Denote by M = U ΣV ∗ the SVD of the 2 × 2 matrix M . Without loss of generality, we assume Σ = diag(σ1 , σ2 ) with σ1 ≥ σ2 > 0, and cos φu − sin φu cos φv − sin φv , V = Rφv = . U = Rφu = sin φu sin φv cos φu cos φv Rφu and Rφv represent the linear transform of rotation (counterclockwise) by angles φu and φv , respectively. ˆ into a triangle centered at the origin. Its effect The mapping x = M ξ maps K can be understood as the composition of three operations (see, e.g., [16]): (1) rotation clockwise by angle φv ; then (2) stretching by factors σ1 and σ2 in ξ and η directions,
23
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
respectively; and finally (3) rotation counterclockwise by angle φu . A circle (centered at (0,0)) in the ξη-plane will be mapped into an ellipse in the xy-plane with two axes σ1 and σ2 , and the longer axis forms an angle φu with the x-axis. See Figure 2. ξ1
σ
φ
ξ
ξ1
2
v
1
σ1
1
φ
u
ˆ to K under the mapping M : R∗ K, ˆ ΣR∗ K, ˆ Rφ ΣR∗ K. ˆ Fig. 2. From K u φv φv φv
We define the aspect ratio of the triangle K as r12 =
(1)
σ1 σ2
and define the orientation of K as the direction of angle φu with the x-axis. It can be seen that K is equilateral if and only if σ1 = σ2 and that K is isosceles if and only if φv = mπ 6 with integer m or σ1 = σ2 . Together with the aspect ratio, φv determines the three internal angles of the triangle. For a fixed aspect ratio, the maximum internal angle of K is a periodic even function of φv with period π3 , and it is decreasing in (0, π6 ). K is an obtuse isosceles triangle when φv = mπ 3 and an acute isosceles triangle when φv = (2m+1)π . If K is an obtuse/acute isosceles, then 6 its orientation is in parallel/perpendicular to its base. Furthermore, let b and h be the length of the base and the height over the base of an isosceles triangle; then √ √ √ σ2 = 2√ r12 = 23 hb , φv = 0 when h ≤ √23 b, σ1 = 33 b, 3 h, σ1 = 23 h, σ2 = 33 b, r12 = √23 hb , φv = π6 when h > 23 b. For general triangles, the aspect ratio r12 can be found as follows. Let σ2 1 σ1 ∈ [1, ∞). + q(K) = 2 σ2 σ1 This quantity can be calculated from the three sides and the area of K as (see formula (10) below) (2)
q(K) =
σ12 + σ22 |1 |2 + |2 |2 + |3 |2 √ = . 2σ1 σ2 4 3|K|
Therefore, r12 = q +
q 2 − 1.
Clearly, the closer r12 is to 1, the closer K is to being equilateral. The same is true for the quantity q(K). Therefore, both r12 and q(K) can be used to measure the closeness of a triangle to being equilateral. Indeed, the reciprocal of the righthand-side expression in (2) was used by Bank and Smith [3] as the “shape regularity
24
WEIMING CAO
quantity” of a triangle. q(K) can also be expressed in terms of the matrix M directly. Note that σ1 and σ2 are eigenvalues of (M T M )1/2 ; it is easy to see that q(K) =
M 2F , 2 det(M )
where · F stands for the Frobenius norm of a matrix. This formula was used by Knupp, Margolin, and Shashkov [10] in certain functionals characterizing the smoothness of the mesh. There are several other ways to define the aspect ratio and the orientation of a triangular element in the finite element analysis. For instance, the aspect ratio is usually defined as the ratio of the length of the longest side over the perpendicular distance from it to the opposite vertex, or as the ratio of the diameter of the triangle over the diameter of the largest inscribed circle in the triangle. The orientation can be defined as the direction of the longest side. It is easy to see that these definitions are equivalent to (1) up to some bounded constants. However, they are not precise enough for describing accurately the behavior of the interpolation errors. 3. Formulas for H 1 - and L2 - norms of the interpolation error. Without loss of generality, assume K is a triangle with its center at the origin and its orientation being the x-axis, i.e., xc = 0, φu = 0. We study the H 1 - and L2 -norms of the error for the linear interpolation of a quadratic function u over K. Since the first order terms of u have no contribution to the error of linear interpolation, we assume in particular that u = 21 x · Hx, where H is the Hessian of u. Since H is a 2 × 2 symmetric matrix, we may decompose it into λ1 RφT h , (3) H = Rφh λ2 where Rφh is the matrix of rotation by an angle φh . We also assume that λ1 ≥ |λ2 |. Other cases can be covered by simply considering the function −u. It is easy to see that the contour lines of u are concentrical ellipses √ (when λ √1 λ2 > 0) or hyperbolas (when λ1 λ2 < 0). Their axes are multiples of 1/ λ1 and 1/ λ2 , and the longer axis (with λ2 ) is of angle φh with the y-axis. We define this direction as the orientation of function u. When φh = π2 , u is oriented along the x-axis, i.e., in the same direction as triangle K. In this case, we say K and u are aligned. When φh = 0, u is oriented along the y-axis, i.e., perpendicular to the orientation of K. See Figure 3.
x1
Fig. 3. Triangle K and the contour lines of a function u. φh is about − π6 in this graph.
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
25
Denote by uI the linear interpolation of u at the three vertices of K. Clearly the norm of the error u − uI depends on the size and the shape of triangle K, as well as the magnitude and the orientation of function u. We derive in this section the exact formulas for the H 1 - and L2 -norms of the interpolation error in terms of these factors. Theorem 3.1. Let K be a triangle oriented along the x-axis. r12 is the aspect ratio of K, r21 = 1/r12 . Let u be a quadratic function. λ1 and λ2 are the eigenvalues of the Hessian matrix of u, and the orientation of u is of angle φh with the y-axis; then √ 3 2 ∇(u − uI )L2 (K) = |K|2 [(r12 + r21 )(λ21 + λ22 ) + (r12 − r21 )(λ21 − λ22 ) cos(2φh )] 36 1
1 − 4λ1 λ2 + ( (λ1 + λ2 )(r12 + r21 ) + (λ1 − λ2 )(r12 − r21 ) + (4) 2 4 · cos(2φh ) )2 [(r12 + r21 + (r12 − r21 ) cos(6φv + 4θ)] ,
ˆ to K, and where φv is the rotation angle in the mapping from K 1 2(λ1 − λ2 ) sin(2φh ) θ = atan (5) . 2 (λ1 + λ2 )(r12 − r21 ) + (λ1 − λ2 )(r12 + r21 ) cos(2φh ) Proof. Step 1. We start with a result established by Bank and Smith in [3]. Let ci (x, y), i = 1, 2, 3, be the barycentrical coordinates of a point (x, y) in K. The side basis functions (taking value 14 at a midpoint) are b1 (x, y) = c2 (x, y) c3 (x, y), b2 (x, y) = c3 (x, y) c1 (x, y), b3 (x, y) = c1 (x, y) c2 (x, y). It is easy to see that 1 u − uI = − (v1 b1 (x, y) + v2 b2 (x, y) + v3 b3 (x, y)), 2
(6) where
vi = i · Hi ,
(7)
i = 1, 2, 3.
Let v = [v1 , v2 , v3 ]T . It is further established in [3] that
1 (8) |∇(u − uI )|2 dxdy = v · Bv, 4 K where
∇bi · ∇bj dxdy K ⎡ |1 |2 + |2 |2 + |3 |2 = 1 ⎣ 48|K| symm.
B=
21 · 2 |1 |2 + |2 |2 + |3 |2
⎤ 21 · 3 ⎦. 22 · 3 |1 |2 + |2 |2 + |3 |2
26
WEIMING CAO
We first derive a formula for the matrix B in terms of the singular values and the rotation angle φv of the mapping M . Note that we assume φu = 0; therefore √ σ1 cos αi ∗ , i = M ei = ΣRφv ei = 3 σ2 sin αi where αi = 2(i − 1)π/3 − φv ,
(9)
i = 1, 2, 3.
With formula (A.1) in the appendix it is easy to verify that (10)
|1 | + |2 | + |3 | = 2
2
2
3σ12
3
2
cos αi +
3σ22
i=1
3
sin2 αi =
i=1
9 2 (σ + σ22 ). 2 1
Let fij =
(11)
σ12
σ12 σ2 cos αi cos αj + 2 2 2 sin αi sin αj . 2 + σ2 σ1 + σ2
Then i · j = 3(σ12 cos αi cos αj + σ22 sin αi sin αj ) = 3(σ12 + σ22 )fij , Note that |K| =
√ 3 3 4 σ1 σ2 .
We may express B as ⎡ 4 √ 2 1 2 3 f12 3 σ1 + σ2 ⎣ 1 B= 24 σ1 σ2 symm.
4 3 f13 4 3 f23
1 ≤ i, j ≤ 3.
⎤ ⎦,
1
and rewrite the norm of the error as √ 2 3 σ1 + σ22 2 (12) ∇(u − uI )L2 (K) = 96 σ1 σ2 3 8 (vi )2 + (f12 v1 v2 + f13 v1 v3 + f23 v2 v3 ) . · 3 i=1 Step 2. We next simplify the terms involving vi . Assume the Hessian matrix H = (hij ), and denote by σ1 ˜ = σ2 h11 σh12 (13) . H 2 h12 σ1 h22 ˜ is a 2 × 2 symmetric matrix, we may decompose it into Since H μ1 ˜ = Rθ RθT , (14) H μ2 where Rθ is the matrix for the counterclockwise rotation by an angle θ, and μ1 and ˜ By the facts that μ2 are the eigenvalues of H. ˜ (R∗ ei ) vi = i · Hi = (Rφ∗ v ei ) · (ΣT HΣ) (Rφ∗ v ei ) = σ1 σ2 (Rφ∗ v ei ) · H φv ∗ ∗ ∗ ∗ = σ1 σ2 (Rθ Rφv ei ) · diag(μ1 , μ2 ) (Rθ Rφv ei )
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
and that Rθ∗ Rφ∗ v ei =
√
3[cos(αi − θ), sin(αi − θ)]T ,
we have vi = 3σ1 σ2 [μ1 cos2 (αi − θ) + μ2 sin2 (αi − θ)] = 32 σ1 σ2 [(μ1 + μ2 ) + (μ1 − μ2 ) cos 2(αi − θ)]. Therefore it follows from formulas (A.2) and (A.3) that 9 3 2 2 2 (vi ) = (σ1 σ2 ) 3(μ1 + μ2 ) + (μ1 − μ2 ) 4 2 i=1 81 2 2 2 2 (σ1 σ2 ) μ1 + μ2 + μ1 μ2 . = 8 3
3
(15)
2
Expand vi vj as follows: (16)
vi vj =
9 (σ1 σ2 )2 [(μ1 + μ2 )2 + (μ21 − μ22 )(cos 2(αi − θ) + cos 2(αj − θ)) 4 + (μ1 − μ2 )2 cos 2(αi − θ) cos 2(αj − θ)].
Let β = φv + θ. Then αi − θ = 2(i − 1)π/3 − β, and ⎧ ⎨ cos(2β + π3 ) cos(2β − π3 ) cos 2(αi − θ) + cos 2(αj − θ) = ⎩ cos(2β + π)
for (i, j) = (1, 2), for (i, j) = (1, 3), for (i, j) = (2, 3)
and
⎧ 1 ⎪ ⎨ −4 − − 14 − cos 2(αi − θ) · cos 2(αj − θ) = ⎪ ⎩ −1 − 4
1 2 1 2 1 2
cos(4β − π3 ) cos(4β + π3 ) cos(4β + π)
for (i, j) = (1, 2), for (i, j) = (1, 3), for (i, j) = (2, 3).
Substituting the above formulas into (17), we have
(17)
f12 v1 v2 + f13 v1 v3 + f23 v2 v3 = 94 (σ1 σ2 )2 {[(μ1 + μ2 )2 − 14 (μ1 − μ2 )2 ] (f12 + f13 + f23 ) + (μ21 − μ22 ) [f12 cos(2β + π3 ) + f13 cos(2β − π3 ) + f23 cos(2β + π)] − 12 (μ1 − μ2 )2 [f12 cos(4β − π3 ) + f13 cos(4β + π3 ) + f23 cos(4β + π)] }.
By using the definition (11) for fij , we have from formula (A.4) that 3 f12 + f13 + f23 = − ; 4
(18)
from formulas (A.5) and (A.6) that f12 cos(2β + π3 ) + f13 cos(2β − π3 ) + f23 cos(2β + π) (19)
= − 34
σ12 − σ22 cos(2θ); σ12 + σ22
27
28
WEIMING CAO
and from formulas (A.7) and (A.8) that f12 cos(4β − π3 ) + f13 cos(4β + π3 ) + f23 cos(4β + π) (20)
= − 34
σ12 − σ22 cos(6φv + 4θ). σ12 + σ22
Put (17)–(20) into (13); we may express the norm of the error as √ 3 3 2 σ12 − σ22 2 2 2 2 2 2 ∇(u − uI )|L2 (K) = (σ + σ2 )σ1 σ2 μ1 + μ2 − 2 (μ1 − μ2 ) cos(2θ) 64 1 σ1 + σ22 1 σ12 − σ22 2 + (μ1 − μ2 ) 1 + 2 (21) cos(6φv + 4θ) . 2 σ1 + σ22 Step 3. Finally we express μ1 , μ2 , and θ in terms of the eigenvalues of the Hessian H and the aspect ratio r12 of K. According to the eigendecomposition (3) of H, λ1 cos2 φh + λ2 sin2 φh (λ1 − λ2 ) sin φh cos φh . H= (λ1 − λ2 ) sin φh cos φh λ1 sin2 φh + λ2 cos2 φh Therefore
˜ ij ) = ˜ = (h H
D sin 2φh r12 (A + D cos 2φh ) D sin 2φh r21 (A − D cos 2φh )
with A=
1 (λ1 + λ2 ), 2
D=
1 (λ1 − λ2 ). 2
Recall the eigendecomposition (14), it is not difficult to establish for the eigenvalues ˜ that μ1 and μ2 of H μ21 + μ22
˜ 11 )2 + (h ˜ 22 )2 + 2(h ˜ 12 )2 = (h 2 2 2 2 2 2 2 = (A + D cos (2φh ))(r12 + r21 ) + 2AD(r12 − r21 ) cos(2φh ) 2 2 + 2D sin (2φh ),
(μ21 − μ22 ) cos 2θ
˜ 11 )2 − (h ˜ 22 )2 = (h 2 2 2 2 2 2 = (A + D cos2 (2φh ))(r12 − r21 ) + 2AD(r12 + r21 ) cos(2φh ),
(μ1 − μ2 )2
˜ 11 − h ˜ 22 )2 + 4(h ˜ 12 )2 = (h 2 2 2 2 = 4D + A (r12 − r21 )2 + 2AD(r12 − r21 ) cos(2φh ) 2 2 2 +D (r12 − r21 ) cos (2φh ) 2 2 − r21 ) cos(2φh ) = 4(D2 − A2 ) + A2 (r12 + r21 )2 + 2AD(r12 2 2 2 +D (r12 − r21 ) cos (2φh ) = − 4λ1 λ2 + [A(r12 + r21 ) + D(r12 − r21 ) cos(2φh )]2
and that tan(2θ) =
˜ 12 2h 2D sin 2φh = . ˜ 22 ˜ 11 − h A(r12 − r21 ) + D(r12 + r21 ) cos 2φh h
Substitute the above formulas into the right-hand side of (21) and simplify, and we obtain the formula (5).
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
29
We next derive the formula for the L2 -norm of the linear interpolation error. Theorem 3.2. Let K be a triangle oriented along the x-axis. r12 is the aspect ratio of K, r21 = 1/r12 . Let u be a quadratic function. λ1 and λ2 are the eigenvalues of its Hessian, and the orientation of u is of angle φh with the y-axis. Then 16 |K|3 2 − λ1 λ2 + [(λ1 + λ2 )(r12 + r21 ) (22) u − uI L2 (K) = 160 9 + (λ1 − λ2 )(r12 − r21 ) cos(2φh )]2 . Proof. It follows from (6) that
1 (23) |u − uI |2 dxdy = v · B0 v, 4 K where
bi · bj dxdy .
B0 = K
An elementary calculation of the integrals leads to ⎡ ⎤ 2 1 1 |K| ⎣ 1 2 1 ⎦. B0 = 180 1 1 2 Therefore
(24)
|K| |u − uI |2 dxdy = 360 K
3
(vi )2 + v1 v2 + v1 v3 + v2 v3
,
i=1
where v1 , v2 , and v3 are defined as in (7). We have similar to (17) that 81 10 2 2 2 (σ1 σ2 ) μ1 + μ2 + μ1 μ2 , v1 v2 + v1 v3 + v2 v3 = 16 3 which, together with (15), yields
|K|3 4 (μ1 + μ2 )2 − μ1 μ2 . (25) |u − uI |2 dxdy = 40 9 K ˜ = det(H) = λ1 λ2 and Finally, by μ1 μ2 = det(H) ˜ 11 + h ˜ 22 = A(r12 + r21 ) + D(r12 − r21 ) cos(2φh ), μ1 + μ2 = h we prove the conclusion of this theorem. 4. Discussion about some special cases. Denote by T1 , T2 , and T3 the terms in the square brackets on the right-hand side of (5), i.e., (26) (27) (28)
T1 = (r12 + r21 )(λ21 + λ22 ) + (r12 − r21 )(λ21 − λ22 ) cos(2φh ), 1 2 T2 = −4λ1 λ2 + ((λ1 + λ2 )(r12 + r21 ) + (λ1 − λ2 )(r12 − r21 ) cos(2φh )) , 4 T3 = r12 + r21 + (r12 − r21 ) cos(6φv + 4θ).
30
WEIMING CAO
Recall that we assume r12 = σ1 /σ2 ≥ 1 and λ1 ≥ |λ2 |; therefore all the three terms above are nonnegative, and we may write √ 3 1 ∇(u − uI )2L2 (K) = (29) |K|2 T1 + T2 · T3 . 36 2 We may also restrict 0 ≤ φv < π3 and 0 ≤ φh < π. Other cases are covered by the symmetry and periodicity of the error norms. For a fixed aspect ratio, the extremum values of T1 , T2 , and T3 are min (T1 ) = 2(r21 λ21 + r12 λ22 ) when φh = π2 , (30) T1 = 2 2 when φh = 0, max (T1 ) = 2(r12 λ1 + r21 λ2 ) 2 when φh = π2 , min (T2 ) = (r21 λ1 − r12 λ2 ) T2 = (31) 2 when φh = 0, max (T2 ) = (r12 λ1 − r21 λ2 ) when 6φv + 4θ = (2m + 1)π, min (T3 ) = 2r21 T3 = (32) when 6φv + 4θ = 2mπ. max (T3 ) = 2r12 We give more details about these cases. Case 1. φh = π2 and φv = π6 , i.e., K is an acute isosceles triangle aligned with u. In this case, θ = 0 by (5). Therefore, T1 , T2 , and T3 all take their minimum values, which yields √ 3 2 |K|2 [ 2(r21 λ21 + r12 λ22 ) + r21 (r21 λ1 − r12 λ2 )2 ]. (33) ∇(u − uI )L2 (K) = 36 This is the smallest H 1 -seminorm of the interpolation error for all the triangles with a fixed aspect ratio and a fixed area. We may consider the optimal choice of the aspect ratio in this case. Let 2 f (r) = 2 1r λ21 + rλ22 + 1r 1r λ1 − rλ2 = 3λ22 r + 2(2λ21 − 4λ1 λ2 ) 1r + λ21 13 . r (1)
(1)
It can be seen that f is decreasing in [1, r∗ ) and increasing in (r∗ , ∞), where ! λ (λ − λ ) + λ (λ − λ )2 + 9λ2 1 1 2 1 1 2 2 (1) (34) . r∗ = 2 3λ2 (1)
Therefore r12 = r∗ is the best aspect ratio. Moreover, since T1 , T2 , and T3 all take minimum values in this case, we conclude that the acute isosceles triangle, which is (1) aligned with u and of the aspect ratio r∗ , is the best one that produces the smallest H 1 -seminorm of the interpolation error among all the triangles with the same area! (1) It is noted that for u with λ1 = λ2 (isotropic u), the best aspect ratio ∗ =1 ! is r√ (1) (with isotropic K). For u with λ2 = −λ1 , the best aspect ratio is r∗ = (2 + 13)/3 ≈ 1.367. When λ1 >> |λ2 |, we have " # # # # # λ1 # 2 ## λ1 ## (1) ≈ 0.816 ## ## , r∗ ≈ # # 3 λ2 λ2
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
31
with which the interpolation error is of the smallest possible magnitude √ 2 2 |K|2 |λ1 λ2 | ≈ 0.235|K|2 |λ1 λ2 |. (35) ∇(u − uI )L2 (K) ≈ 6 We may also compare the above optimal choice of the aspect ratio with some intuitive choices. For instance, if we choose r12 = 1, then √ 3 2 |K|2 (3(λ21 + λ22 ) − 2λ1 λ2 ), (36) ∇(u − uI )L2 (K) ≈ 36 which is between 0.128|K|2 λ21 and 0.193|K|2 λ21 . If we choose r12 = | λ1 |, then in the λ2 case λ1 >> |λ2 |, √ 5 3 2 |K|2 |λ1 λ2 | ≈ 0.240|K|2 |λ1 λ2 |; ∇(u − uI )L2 (K) ≈ (37) 36 if we choose r12 = |λ1 /λ2 |, then we have √ ! ! 3 2 ∇(u − uI )L2 (K) = |K|2 |λ31 λ2 | ≈ 0.0962|K|2 |λ31 λ2 |. (38) 18 Note that when u and K are aligned, the inverse mapping ξ = M −1 x with r12 = |λ1 /λ2 | transforms u(x, y) into const.(ξ 2 ± η 2 ) and K into an equilateral triangle. It was shown by D’Azevedo and Simpson [6] and Rippa [13] that triangles with such an aspect ratio lead to the smallest maximum-norm and Lp -norm of the interpolation error among all triangles of the same area. However, this aspect ratio is not the optimal for the interpolation error in the sense of H 1 -seminorm. Case 2. φh = π2 and φv = 0, i.e., K is an obtuse isosceles triangle aligned with u. We have in this case √ 3 2 |K|2 [ 2(r21 λ21 + r12 λ22 ) + r12 (r21 λ1 − r12 λ2 )2 ]. (39) ∇(u − uI )L2 (K) = 36 This error formula differs from that of Case 1 only in the coefficient of the second term. To study the optimal choice of the aspect ratio in this case, let 2 g(r) = 2 1r λ21 + rλ22 + r 1r λ1 − rλ2 = 3λ21 1r + 2(λ22 − λ1 λ2 )r + λ22 r3 . It can be shown that the minimum of g is attained at ! λ (λ − λ ) + |λ | (λ − λ )2 + 9λ2 2 1 2 2 1 2 2 (2) r∗ = (40) . 2 3λ2 (2)
When λ1 = λ2! , the best aspect ratio is r∗ = 1. When λ2 = −λ2 , the best aspect √ (1) ratio is r∗ = (2 + 13)/3 ≈ 1.367. When |λ1 | >> |λ2 |, we have ⎧ "√ " ⎪ λ 10 + 1 1 ⎪ | | ≈ 1.178 | λ1 | when λ1 λ2 > 0, ⎨ 3 λ2 λ (2) " " 2 r∗ ≈ √ ⎪ ⎪ 10 − 1 | λ1 | ≈ 0.849 | λ1 | ⎩ when λ1 λ2 < 0. 3 λ2 λ2
32
WEIMING CAO
With this choice of aspect ratio, the interpolation error is of the magnitude when λ1 λ2 > 0, 0.0878|K|2 |λ31 λ2 | 2 ∇(u − uI )L2 (K) ≈ (41) 2 3 when λ1 λ2 < 0. 0.281|K| |λ1 λ2 | We may compare this case (obtuse isosceles) with the previous one (acute isosceles). When the triangle is aligned with u and the aspect ratio r12 ≈ |λ1 /λ2 |, the interpolation error is nearly the same for both acute and obtuse triangles if λ1 λ2 > 0. For u with λ1 λ2 < 0, ∇(u − uI )L2 (K) is about 1.7 times smaller with the acute isosceles triangle than that with the obtuse one. Furthermore, because for general triangles with arbitrary φv , the error norm ∇(u − uI )L2 (K) is between those of Case 1 and Case 2, we may conclude that the H 1 -seminorm of the interpolation error is not sensitive to the maximum angle of the triangle, as long as the triangle is aligned with the function u and the aspect ratio r12 ≈ |λ1 /λ2 | is used. Case 3. φh = 0. In this case the orientation of u is perpendicular to the triangle. We have √ 3 2 |K|2 [2(r12 λ21 + r21 λ22 ) + (r12 λ1 − r21 λ2 )2 ∇(u − uI )L2 (K) = 18 (42) · (r12 + r21 + (r12 − r21 ) cos(6φv ))] For given |K| and φv , it can be shown the error norm is an increasing function of r12 . Therefore, the best aspect ratio is r12 = 1. This is not surprising, since when the triangle is perpendicular to the orientation of u, increasing the aspect ratio leads to lower resolution in the needy direction and larger interpolation error, no matter what value φv takes. We should emphasize that in the above discussion on the optimal aspect ratios the area and the orientation of the triangle are assumed fixed. In adaptive mesh generation there may be some other constraints on the triangles in a mesh, e.g., one side of a triangle or the vertices of the triangles are fixed. Optimal choices of the aspect ratio subject to those constraints can be quite different from the above values. 5. General discussion about orientation, aspect ratio, and angle condition. In this section we present a general discussion about the effects of the orientation, aspect ratio, and the angle φv . We first study the H 1 -seminorm of the interpolation error. Mesh alignment. Given a triangle K (fixed |K|, r12 , and φv ), it is commonly believed that the interpolation error is the smallest when the triangle is aligned with the function, i.e., when φh = π2 . We present a justification of this viewpoint. We consider the case where both the solutions and the triangles are highly anisotropic, i.e., we assume λ1 >> |λ2 | and r12 >> 1. In this case we have for the three terms (26)–(28) in the error norm ∇(u − uI )2L2 (K) that T1 = r12 λ21 (1 + cos 2φh ) + LOT, T2 = 14 (r12 λ1 )2 (1 + cos 2φh )2 + LOT, T3 = r12 (1 + cos(6φv + 4θ)) + LOT and for the angle θ in (5) that tan(2θ) ≈
2λ1 sin 2φh = 2r21 tan φh , λ1 r12 (1 + cos 2φh )
33
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
where the lower order terms LOT include r12 λ22 , r12 λ21 , r21 λ22 , etc. In this case the error norm is dominated by ∇(u − uI )2L2 (K) ≈
√
3 2 36 |K|
3 2 [r12 λ21 + 14 r12 λ1 (1 + cos 2φh )(1 + cos(6φv + 4θ))] · (1 + cos 2φh ) + LOT.
Note that by (5) θ is an increasing function of φh , and cos(2φh ) is a decreasing function of φh in (0, π2 ). Therefore, we may conclude that the H 1 -seminorm of the error is almost a monotonically decreasing function of the angle φh in (0, π2 ). We plot in Figure 4 the graph of ∇(u − uI )2L2 (K) versus φh for various λ1 , λ2 , φv and r12 . It is observed that for all aspect ratios and φv , the error is smallest when φh is near π2 . Also seen in Figure 4 is a small peak of the error norm around φh = π2 . This is because when φh is very close to π2 , the leading terms in T1 and T2 go to 0 and the lower order terms become the major components in the error. 4
3
10
10
2
(K)
L
||∇(u-u )||2
v
2
10
φv=0 φv=π/8 φ =π/6
2
10
||∇(u-u L)||(K)
φv=0 φv=π/8 φ =π/6
3
10
v
1
10
1
10
0
10
0
0
0.2 0.4 0.6 0.8 with λ1=10, λ2=1, r12=8
10
1
10
0.2 0.4 0.6 0.8 with λ1=10, λ2=1, r12=4
1
8
10
10
φv=0 φv=π/8 φv=π/6
φv=0 φv=π/8 φv=π/6
6
10 2
||∇(u-u L)||(K)
2
||∇(u-u L)||(K)
0
5
10
4
10
2
10 0
10
0
0
0.2 0.4 0.6 0.8 with λ1=100, λ2=1, r12=80
Fig. 4. u − uI 2 2
L (K)
1
10
0
0.2 0.4 0.6 0.8 with λ1=100, λ2=1, r12=40
versus the alignment angle φh /π with |K| =
3
√ 3 4
1
fixed.
Aspect ratio. We are interested here in the question of, given a function u and a fixed triangle orientation φh , how one chooses the aspect ratio of the triangle. We plot in Figure 5 the graph of ∇(u − uI )L2 (K) versus r12 with various φh . It is observed that when φh = π2 (K aligned with u), the best aspect ratio r12 is in between (2) (1) r∗ = |λ1 /λ2 | and r∗ ≈ 0.8|λ1 /λ2 |, while when φh ≈ 0 (K perpendicular to u), the best aspect ratio is 1. For other φv and φh , the optimal choice of r12 lies between alignment is ensured, the 1 and |λ1 /λ2 |. Therefore, in practice, when a good mesh proper aspect ratio should be between the magnitude of |λ1 /λ2 | and 0.8|λ1 /λ2 |, with the lower end taken for the mesh with mostly obtuse triangles and the upper
34
WEIMING CAO 4
5
10
10
φv=0 φv=π/8 φ =π/6
2
||∇(u-u L)||(K)
v
2
10
v
3
10
1
2
10
10
0
10
φv=0 φv=π/8 φ =π/6
4
10 ||∇(u-u L)||2(K)
3
10
1
0
10
20 with φ =π/2
30
40
10
0
10
h
6
6
φ =0 v φv=π/8 φv=π/6
4
φ =0 v φv=π/8 φv=π/6
4
10
||∇(u-u L)||2(K)
2
40
10
10
||∇(u-u L)||(K)
30
h
10
2
10
2
10
0
10
20 with φ =π/3
0
0
10
20 with φh=π/6
Fig. 5. ∇(u − uI )2 2
L (K)
30
40
10
0
10
20 with φh=0
30
versus the aspect ratio r12 with λ1 = 10, λ2 = 1, and |K| =
40
3
√ 3 . 4
end taken for the mesh with mostly acute triangles. The less the mesh alignment, the smaller r12 should be. Maximum angle condition. It is well known that a good mesh for finite element approximation should respect the so-called maximum angle condition, i.e., the maximum angles of the triangular elements should be bounded from π. In [2], Babu˘ska and Aziz gave an example showing that when the maximum angle condition is violated, the H 1 -seminorm of the error of linear interpolation can go to infinity, although the element area approaches 0. Their example is in our nota√ tion the case with λ1 = 1, λ2 = 0, |K| = 2 , r12 = 2 3 , φh = 0, φv = 0, and thus 9 3 ∇(u − uI )2L2 (K) = 256 ( + 8
). We discuss here the impact of the maximum angle condition on the accuracy of linear interpolation. Taking as an example with λ1 = 100 and λ2 = 1, we plot in Figure 6 the graph of ∇(u − uI )2L2 (K) versus φv with given aspect ratios and alignment directions. It is easy to see the following: (a) When λ1 /λ2 and r12 are high, triangles with φv ≈ π6 or π2 produce the smallest error, and triangles with φv ≈ 0 or π3 produce the largest error, whether they are aligned with u or not. This can be justified by looking into formula (5) for the error norm. When λ1 >> |λ2 | and r12 >> 1, we have θ ≈ 0. Hence, the term T3 , and therefore ∇(u − uI )L2 (K) , attains its minimum when φv ≈ π6 or π2 . Note that the maximum internal angle of K is an even periodic function of φv with the period π/3. It is decreasing in (0, π/6). Hence φv = 0, π/3 corresponds to the case where the maximum internal angle of K is the largest among all the triangles of the same aspect ratio, and φv = π/6, π/2 corresponds to the case where
35
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
the maximum internal angle is the smallest among all triangles. Therefore, in this case it is important to use triangles with smaller maximum internal angles, regardless of whether or not the mesh alignment is satisfied. (b) When the triangle is aligned with u and the aspect ratio r12 ≈ |λ1 /λ2 |, the error is insensitive to φv and is therefore insensitive to the maximum internal angle of the triangle. In particular, if φh = π2 , λ1 λ2 > 0, and r12 = λ1 /λ2 , then T2 = 0, and by (29) ∇(u − uI )L2 (K) is independent of φv . See the lower left graph of Figure 6 and the analysis in Cases 1 and 2 in section 4. (c) When r12 ≈ 1, the best φv may be different from π6 and π2 . However, in this case the difference between the maximum and the minimum of the error is not significant. 10
8
10
φ =0 h φh=π/6 φh=π/3 φh=π/2
5
10
φ =0 h φh=π/6 φh=π/3 φh=π/2
6
10 2
||∇(u-u L)||(K)
2
||∇(u-u L)||(K)
10
4
10
2
10 0
10
0
0
0.2 0.4 0.6 0.8 with λ =100, λ =1, r =80 1
2
10
1
6
2
1
12
4
10
2
2
4
10
φh=0 φh=π/6 φ =π/3 h φh=π/2
||∇(u-u L)||(K)
φh=0 φh=π/6 φ =π/3 h φh=π/2
5
10 L
0.2 0.4 0.6 0.8 with λ =100, λ =1, r =20 1
10
||∇(u-u )||(K)
0
12
3
10
2
10
3
0
0.2 0.4 0.6 0.8 with λ1=100, λ2=1, r12=10
Fig. 6. ∇(u − uI )2 2
L (K)
1
10
0
0.2 0.4 0.6 0.8 with λ1=100, λ2=1, r12=1.5
versus the angle φv /π with |K| =
3
√ 3 4
1
fixed.
In summary, we conclude that ifthe triangle is aligned with the function and the aspect ratio is of an magnitude |λ1 /λ2 | or smaller, then the maximum angle condition is not essential to the H 1 -seminorm of the interpolation error. Otherwise, the maximum angle is critical to the magnitude of the error. L2 -norm. From (23) we can conclude the following for the L2 -norm of the linear interpolation error: (1) u−uI L2 (K) does not depend on the angle φv . Therefore, for a given function, the L2 -norm of its linear interpolation error depends only on the area, aspect ratio, and orientation of the triangle. Maximum/minimum angle conditions are irrelevant to the L2 -norm of the interpolation error. (2) u − uI L2 (K) is a π-periodic even function of φh . It is decreasing in φh in
36
WEIMING CAO
(0, π2 ). When φh =
π 2
(K aligned with u),
u − uI 2L2 (K) =
|K|3 40
4 (λ1 r21 + λ2 r12 )2 − λ1 λ2 9
is the smallest, and when φh = 0 (K perpendicular to u) |K|3 4 2 2 (λ1 r12 + λ2 r21 ) − λ1 λ2 u − uI L2 (K) = 40 9 is the largest. (3) With the same orientation, it can be shown that when φh is in (0, π4 ) ∪ ( 3π 4 , π) (i.e., K is aligned more to the perpendicular direction of u), the best aspect ratio is r12 = 1; when π4 ≤ φh ≤ 3π 4 (i.e., K is aligned more to the orientation of u), the best aspect ratio is $# # # λ1 + λ2 − (λ1 − λ2 ) cos(2φh ) # #. (43) r12 = ## λ1 + λ2 + (λ1 − λ2 ) cos(2φh ) # We plot in Figure 7 the L2 -norm of the error versus the aspect ratio for the case of √ λ1 = 400, λ2 = 1, and |K| = 3 4 3 . 8
10
φ =0 h φh=π/6 φ =π/4 h φ =π/3 h φ =2π/5 h φh=π/2
7
10
6
10
5
||u-u ||2
L (K)
10
4
10
3
10
2
10
1
10
Fig. 7. u − uI 2 2
0
L (K)
5
10
15
20 25 aspect ratio r12
30
35
40
45
versus the aspect ratio r12 with λ1 = 400, λ2 = 1 and |K| =
3
√ 3 . 4
(4) For a fixed triangle area, the minimum L2 -norm of the linear interpolation π error is attained at φh = 2 and r12 = |λ1 /λ2 |. The minimum value is 8 3 when λ1 λ2 > 0, 2 90 |K| |λ1 λ2 | u − uI L2 (K) = 1 3 |K| |λ λ | when λ1 λ2 < 0. 1 2 90 6. An example of linear interpolation on different meshes. In this section, we present the results of piecewise linear interpolation of a quadratic function on different types of mesh. We choose u(x, y) = x2 + 100y 2 ,
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
(a)
(b)
37
(c)
Fig. 8. Three types of mesh used for piecewise linear interpolation.
which corresponds to λ1 = 100, λ2 = 1, and the orientation along the x-axis. We consider the following three types of mesh on [0, 1] × [0, 1]: (a) The first type of mesh is as shown in Figure 8(a). There√ are M and N equal partitions in x and y directions, respectively. When M/N < 2/ 3, all the triangles in this mesh, except those on the top and the bottom boundaries, are aligned exactly with the x-axis. Their aspect ratio is 2N r12 = √ 3M and their angle φv is π6 . (b) The second type of mesh is as shown in Figure 8(b)√with M and N equal partitions in x and y directions, respectively. When M/N < 3/2, all the triangles, except those on the left and the right boundaries, are aligned with the x-axis. Their aspect ratio is √ 3N r12 = 2M and their angle φv is 0. (c) The third type is the unstructured meshes generated in the following way. First we create a Delaunay triangulation by using the package Triangle [14] on a rectangle [0, 1] × [0, rs ], where rs ≥ 1 is a real number. Then we compress the rectangle in y direction by a factor rs . The desired mesh on the unit square is the image of the Delaunay triangulation under the compression. See Figure 8(c) for a typical mesh of this type. By Delaunay triangulation the minimum internal angle of the triangles in the mesh is maximized, and all the triangles in the mesh over [0, 1] × [0, rs ] are approximately of the same size and close to unilateral. Therefore, most triangles in the mesh over the unit square are roughly aligned with the x-axis (when rs > 1) and are of the aspect ratio r12 ≈ rs . However, in this type of mesh the angle φv varies for different triangles (approximately uniformly distributed between 0 and π/3). We are interested in the accuracy of the linear interpolation of u on different types of mesh and with various aspect ratios, in particular with the following three (1) ratios: (1) r∗ = 81.6, which is the best aspect ratio (for H 1 -seminorm) calculated
38
WEIMING CAO (2)
according to (34) for the case of acute isosceles triangles; (2) r∗ = 11.79, which is the best aspect ratio (in H 1 -seminorm) calculated according to (40) for the case ! (∞)
of obtuse isosceles triangles; and (3) r∗ = | λλ12 | = 10, which is the best aspect ratio for the Lp -norm, 1 ≤ p ≤ ∞. We also report the results with the aspect ratios r12 = λ1 /λ2 = 10 and r12 = 1. We list in Table 1 the various norms of the interpolation error with different choices of M , N , and rs . Note that for all the meshes, the total number of elements is around 4000 and the total number of nodes is around 2200. The smallest H 1 seminorm for type (a) meshes is with r12 = 92.37; the smallest H 1 -seminorm for type (b) meshes is with r12 = 11.54. For type (c) meshes, the smallest H 1 -seminorm is (2) achieved at rs = 11.78 ≈ r∗ . This is because for this type of mesh, the error from the obtuse triangles (with φv ≈ 0) dominates in the global error. Therefore, the aspect (2) ratio close to r∗ is the best choice for the global error norm. For the L2 -, L1 -, and maximum norms, the smallest interpolation error is obtained (∞) with r12 ≈ r∗ for all the three meshes. In summary, we conclude that when the mesh is in good alignment with the solution and most of the triangles are acute isosceles, the aspect ratio should be (1) chosen around r∗ ≈ 0.8| λλ12 |. If the mesh is in good alignment but with mostly obtuse triangles, or with varied maximum internal angles, then the aspect ratio should be (2) chosen around r∗ , which is approximately 1.178 |λ1 /λ2 | for the case λ1 λ1 > 0, and 0.849 |λ1 /λ2 | for the case λ1 λ1 < 0. Table 1 The L1 - and L2 -norms, H 1 -seminorm, and maximum norm of the interpolation error over the entire domain. (∗1) and (∗2) indicate the aspect ratio close to the best values for the H 1 seminorm in each case. M ×N
Node #
Elem #
r12
u−uI L1
u−uI L2
|u−uI |H 1
u−uI ∞
Type (a) meshes 4×500
2507
4004
144.34
5.23e−3
6.18e−3
8.28e−2
7.82e−3
5×400
2409
4005
92.37(∗1)
3.37e−3
3.97e−3
7.68e−2
5.01e−3
6×333
2341
4002
64.08
2.37e−3
2.78e−3
7.76e−2
3.50e−3
14×143
2167
4018
11.79
7.30e−4
8.01e−4
1.43e−1
9.74e−4
15×133
2152
4005
10.24
7.22e−4
7.92e−4
1.54e−1
9.64e−4
48×42
2131
4080
1.01
3.55e−3
4.20e−3
5.88e−1
7.08e−3
Type (b) meshes 3×572
2578
4004
165.11
6.19e−3
7.60e−3
7.22e+0
1.38e−2
4×445
2543
4005
96.34
3.62e−3
4.39e−3
3.21e+0
7.81e−3
5×364
2372
4004
63.04
2.39e−3
2.87e−3
1.68e+0
5.00e−3
12×160
2173
4000
11.54(∗2)
7.47e−4
8.23e−4
1.53e−1
1.01e−3
13×149
2175
4023
9.92
7.35e−4
8.08e−4
1.60e−1
9.93e−4
41×49
2125
4067
1.03
3.50e−3
4.13e−3
5.87e−1
5.22e−3
Type (c) meshes 2474
4101
100
4.24e−3
5.34e−3
2.95e+0
1.91e−2
2426
4079
81.6
3.55e−3
4.57e−3
2.48e+0
1.67e−2
2163
4036
11.78
8.38e−4
9.62e−4
1.84e−1
2.42e−3
2127
4005
10
8.34e−4
9.53e−4
1.87e−1
3.05e−3
2092
4036
1
4.15e−3
5.21e−3
6.75e−1
2.01e−2
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
39
Appendix. We list here some basic trigonometric formulas used to derive the results in the previous sections. Let α be any real number. Denote by αi = 2(i − 1)π/3 − α, i = 1, 2, 3.
(A.2)
3 , 2 cos(2α1 ) + cos(2α2 ) + cos(2α3 ) = 0,
(A.3)
cos2 (2α1 ) + cos2 (2α2 ) + cos2 (2α3 ) =
(A.4)
3 cos(α1 ) cos(α2 ) + cos(α1 ) cos(α3 ) + cos(α2 ) cos(α3 ) = − . 4
(A.1)
cos2 (α1 ) + cos2 (α2 ) + cos2 (α3 ) =
3 , 2
Similar relations hold for sine functions, too. Let θ be any real number, and let β = α + θ. Then π π + cos(α1 ) cos(α3 ) cos 2β − (A.5) cos(α1 ) cos(α2 ) cos 2β + 3 3 3 + cos(α2 ) cos(α3 ) cos(2β + π) = − cos(2θ), 4 π π + sin(α1 ) sin(α3 ) cos 2β − (A.6) sin(α1 ) sin(α2 ) cos 2β + 3 3 3 + sin(α2 ) sin(α3 ) cos(2β + π) = cos(2θ), 4 π π + cos(α1 ) cos(α3 ) cos 4β + (A.7) cos(α1 ) cos(α2 ) cos 4β − 3 3 3 + cos(α2 ) cos(α3 ) cos(4β + π) = − cos(6α + 4θ), 4 π π + sin(α1 ) sin(α3 ) cos 4β + (A.8) sin(α1 ) sin(α2 ) cos 4β − 3 3 3 + sin(α2 ) sin(α3 ) cos(4β + π) = cos(6α + 4θ). 4 REFERENCES [1] T. Apel, Anisotropic Finite Elements: Local Estimates and Applications, Adv. Numer. Math., Teubner, Stuttgart, 1999. [2] I. Babu˘ ska and A. K. Aziz, On the angle condition in the finite element method, SIAM J. Numer. Anal., 13 (1976), pp. 214–226. [3] R. E. Bank and R. K. Smith, Mesh smoothing using a posteriori error estimates, SIAM J. Numer. Anal., 34 (1997), pp. 979–997. [4] M. Berzins, A solution-based triangular and tetrahedral mesh quality indicator, SIAM J. Sci. Comput., 19 (1998), pp. 2051–2060. ´mal, Triangular elements in the finite element method, Math. Comp., [5] J. Bramble and M. Zla 24 (1970), pp. 809–820. [6] E. F. D’Azevedo and R. B. Simpson, On optimal interpolation triangle incidences, SIAM J. Sci. Statist. Comput., 10 (1989), pp. 1063–1075. [7] L. Formaggia and S. Perotto, New anisotropic a priori error estimates, Numer. Math., 89 (2001), pp. 641–667. [8] W. Huang, Measuring mesh qualities and application to variational mesh adaptation, SIAM J. Sci. Comput., 26 (2005), pp. 1643–1666. [9] W. Huang and W. Sun, Variational mesh adaptation II: Error estimates and monitor functions, J. Comput. Phys., 184 (2003), pp. 619–648. [10] P. Knupp, L. G. Margolin, and M. Shashkov, Reference Jacobian optimization-based rezone strategies for arbitrary Lagrange Eulerian methods, J. Comput. Phys., 176 (2002), pp. 93– 128.
40
WEIMING CAO
[11] G. Kunert, A Posteriori Error Estimation for Anisotropic Tetrahedral and Triangular Finite Element Meshes, Ph.D. dissertation, TU Chemnitz, Chemnitz, Germany, 1999. [12] E. J. Nadler, Piecewise Linear Approximation on Triangulations of a Planar Region, Ph.D. thesis, Division of Applied Mathematics, Brown University, Providence, RI, 1985. [13] S. Rippa, Long and thin triangles can be good for linear interpolation, SIAM J. Numer. Anal., 29 (1992), pp. 257–270. [14] J. R. Shewchuk, Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator, Technical report, School of Computer Science, Carnegie Mellon University, Pittsburgh, PA, 1996. [15] R. B. Simpson, Anisotropic mesh transformations and optimal error control, Appl. Numer. Math., 14 (1994), pp. 183–198. [16] L. N. Trefethen and D. Bau, Numerical Linear Algebra, SIAM, Philadelphia, 1997.
SIAM J. NUMER. ANAL. Vol. 43, No. 1, pp. 19–40
ON THE ERROR OF LINEAR INTERPOLATION AND THE ORIENTATION, ASPECT RATIO, AND INTERNAL ANGLES OF A TRIANGLE∗ WEIMING CAO† Abstract. In this paper, we attempt to reveal the precise relation between the error of linear interpolation on a general triangle and the geometric characters of the triangle. Taking the model problem of interpolating quadratic functions, we derive two exact formulas for the H 1 -seminorm and L2 -norm of the interpolation error in terms of the area, aspect ratio, orientation, and internal angles of the triangle. These formulas indicate that (1) for highly anisotropic triangular meshes the H 1 -seminorm of the interpolation error is almost a monotonically decreasing function of the angle between the orientations of the triangle and the function; (2) maximum angle condition is not |λ1 /λ2 | or essential if the mesh is aligned with the function and the aspect ratio is of magnitude less, where λ1 and λ2 are the eigenvalues of the Hessian matrix of the function. With these formulas we identify the optimal triangles, which produce the smallest H 1 -seminorm of the interpolation error, 1 to be the acute isosceles aligned with the solution and of an aspect ratio about 0.8| λ |. The L2 λ2 norm of the interpolation error depends on the orientation and the aspect ratio of the triangle, but 2 not directly on its maximum or minimum angles. The optimal triangles for the L -norm are those aligned with the solution and of an aspect ratio |λ1 /λ2 |. These formulas can be used to formulate more accurate mesh quality measures and to derive tighter error bounds for interpolations. Key words. anisotropic mesh, linear interpolation, aspect ratio, mesh alignment, maximum angle condition AMS subject classifications. 65D05, 65L50, 65N15, 65N50 DOI. 10.1137/S0036142903433492
1. Introduction. It is well known that on quasi-uniform meshes the accuracy of piecewise linear interpolation is first order in the H 1 -norm. More precisely, denote by H m (Ω) the usual Sobolev spaces of order m on a bounded domain Ω ∈ R2 . For any u ∈ H 2 (Ω), the H 1 -norm of the error between u and its piecewise linear interpolation uI can be bounded by u − uI H 1 (Ω) ≤ ch|u|H 2 (Ω) , where h is the diameter of the triangles in the mesh, and c is a constant independent of h and u. If the mesh is not quasi-uniform, then the error contributed from each triangle K can be bounded by u − uI H 1 (K) ≤ c
h2 |u|H 2 (K) , ρ
where ρ is the diameter of the largest inscribed circle in K. This error bound guarantees the H 1 -norm of the error converge to 0 as h → 0, as long as hρ remains bounded. This is equivalent to requiring that the minimum angle of K is bounded from 0 [5]. ∗ Received
by the editors August 19, 2003; accepted for publication (in revised form) May 7, 2004; published electronically April 26, 2005. This research was supported in part by the NSF under grant DMS-0209313. http://www.siam.org/journals/sinum/43-1/43349.html † Department of Applied Mathematics, University of Texas at San Antonio, San Antonio, TX 78249 ([email protected]). 19
20
WEIMING CAO
However, this error estimate is not tight when the triangle has only one small angle. Babu˘ska and Aziz [2] showed that the minimum angle condition is actually not essential for the convergence of linear interpolation. They improved the error bound as u − uI H 1 (K) ≤ Γ(α)h|u|H 2 (K) , where Γ(α) is an increasing function of the maximum angle α of K. Therefore, in order to guarantee the convergence it is only required that the maximum angle be bounded from π. This is the well-known maximum angle condition. This condition is necessary and sufficient for the convergence of the linear interpolation process over the class of functions of H 2 (K). In adaptive computation, one often knows a priori or a posteriori some information of the functions to be approximated. This information can be used to restrict the set of functions considered and to avoid the divergence case although the triangle has a maximum angle close to π. Indeed, long and thin triangles violating the maximum angle condition have been used successfully in engineering, particularly in computational fluid dynamics [1, 15]. For problems with very different length scales in different spatial directions, long and thin triangles turn out to be better choices than shape regular ones if they are properly used. This motivated an intensive study on the error analysis for anisotropic meshes in the finite element method (FEM). For instance, Apel [1] described an error estimate in terms of the length scales h1 and h2 along the x and y directions, respectively, |u − uI |W m,p (K) ≤ c
α1 +α2 =−m
1 α2 hα 1 h2 |
∂ −m u |W m,p (K) , ∂xα1 ∂y α2
m = 0, 1,
where W m,p (K) is the Sobolev space of functions whose up to mth order derivatives are Lp -integratable. If u and K are aligned and the maximum angle condition is satisfied, this estimate is asymptotically accurate, i.e., both sides of the above inequality have the same hi order. But when the maximum angle of K approaches π, the ratio of the error bound over the actual error norm goes to infinity [7]. Formaggia and Perotto [7] presented another type of estimate of the H 1 -norm for the interpolation error based on the eigendecomposition of the affine mapping from a standard element to K. Their estimate is accurate when u and K are aligned and the maximum angle of K is close to π. However, the ratio of their estimate over the actual error norm goes to infinity when two angles of K approach π2 . More recent error analyses have been based on the overall mesh properties and the behavior of the approximation functions. For instance, Berzins [4] developed a mesh quality indicator measuring the correlation between the anisotropic features of the mesh and those of the solutions. Kunert [11] proposed a so-called matching function to quantify how good overall a mesh is for a specific solution. Huang [8] introduced measures for three aspects of the mesh qualities—aspect ratio, alignment, and adaptation—and an overall quality mesh measure based on them. He formulated the error bounds in terms of these measures and proposed a variational formulation to optimize the overall mesh quality measure to control the interpolation error. A similar idea was used by Huang and Sun [9] in formulating the monitor function in variational mesh adaptation. Needless to say, the interpolation error depends on the solution and the size and shape of the elements in the mesh. Understanding this relation is crucial for the
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
21
generation of efficient and effective meshes for the FEM. However, in all the error estimates for anisotropic meshes, the relation between the error and the geometric characters of a triangle, such as the alignment, aspect ratio, and internal angles, has not been revealed explicitly. In the mesh generation community, this relation is studied more closely for the model problem of interpolating quadratic functions. This model is a reasonable simplification of the cases involving general functions, since quadratic functions are the leading terms in the local expansion of the linear interpolation errors. For instance, Nadler [12] derived an exact expression for the L2 -norm of the linear interpolation error in terms of the three sides 1 , 2 , and 3 of the triangle K, u − uI 2L2 (K) =
|K| [(d1 + d2 + d3 )2 + d1 d2 + d2 d3 + d1 d3 ], 180
where |K| is the area of the triangle, di = i · Hi with H being the Hessian matrix of u. Bank and Smith [3] gave a similar formula for the H 1 -seminorm of the linear interpolation error; see (8) in section 3. Assuming u = λ1 x2 + λ2 y 2 , D’Azevedo and Simpson [6] derived the exact formula for the maximum norm of the interpolation error u − uI L2 ∞ (K) =
D12 D23 D31 , 16λ1 λ2 |K|2
where Dij = i · diag(λ1 , λ2 )i . Based on the geometric interpretation of this formula, they proved that for a fixed area the optimal triangle, which produces the smallest maximum interpolation error, √ √ is the one obtained by compressing an equilateral triangle by factors λ1 and λ2 along the two eigenvectors of the Hessian matrix of u. Furthermore, the optimal incidence for a given set of interpolation points √ √ is the Delaunay triangulation based on the stretching map (by factors λ1 and λ2 along the two eigenvector directions) of the grid points. Rippa [13] showed that the mesh obtained this way is also optimal for the Lp -norm of the error for any 1 ≤ p ≤ ∞. Thought these formulas are exact, they do not describe explicitly the relation between the error and the geometric characters of the triangle. In this paper, we attempt to reveal this relation explicitly and precisely. Taking the model problem of linear interpolation of quadratic functions, we derive two exact expressions for the H 1 -seminorm and L2 -norm of the interpolation error in terms of the area, aspect ratio, alignment direction, and internal angles of the triangle. From these formulas the effects of the geometric characters of a triangle can be clearly identified. They indicate that (1) for highly anisotropic triangular meshes the H 1 -seminorm of the interpolation error is almost a monotonically decreasing function of the angle between the orientation of the triangle and the orientation of the function; (2) maximum angle condition is critical if the mesh is not aligned with the function or if the aspect ratio is larger than |λ1 /λ 2 |; (3) if the triangles are aligned with the function and the aspect ratio is about |λ1 /λ2 | or less, then the error is not sensitive to the maximum angle of the triangle. Also, we can easily identify the best triangles which produce the smallest H 1 - or L2 -norm of the interpolation error. It turns out that in the sense of the H 1 -seminorm, the optimal triangle with a given area is the acute isosceles aligned with the function and of the aspect ratio about 0.8| λλ12 |. The L2 -norm of the interpolation error depends on the alignment and the aspect ratio of the triangle, but not directly on its maximum or minimum angles. The optimal triangles for the L2 -norm are those aligned with the solution and of an aspect ratio |λ1 /λ2 |.
22
WEIMING CAO
The organization of this paper is as follows. In section 2 we give a precise definition of the orientation and the aspect ratio of a triangle by using the singular value decomposition (SVD) of the affine mapping from a standard element to the triangle. In section 3 we derive the exact formulas for the H 1 -seminorm and L2 -norm of the interpolation error for the quadratic functions. Then we identify in section 4 a number of special cases that are of interests to the optimal design of meshes. Optimal choices of the aspect ratio and the alignment direction are discussed here. In section 5 we elaborate in more detail the effects of various geometric characters of a triangle on the interpolation error. Finally we demonstrate by an example the accuracy of the linear interpolation errors on different meshes. 2. Mapping and the geometric characters of a triangle. Let K be a triangle in the xy-plane, and let x1 , x2 , and x3 be the vertices of K. Denote by i ˆ be the the vector of the side opposite to xi (in counterclockwise direction). Let K equilateral triangle with the vertices √ √ 3 0 − 23 2 , ξ3 = . ξ1 = , ξ2 = 1 1 −2 − 12 √ ˆ are ei = 3[cos(2(i − 1)π/3), sin(2(i − 1)π/3)]T , i = 1, 2, 3. See The three sides of K Figure 1. η
ξ
1
y x
1
1
ξ2
ξ3
ξ
x
3
x
2
x
ˆ and physical element K. Fig. 1. Standard element K
Let xc = 31 (x1 + x2 + x3 ) be the center of K. The affine mapping (which maps ˆ to K can be expressed as x = M ξ + xc with ξ i to xi ) from K 1 M = √ (x1 − x3 ), x2 − xc . 3 Denote by M = U ΣV ∗ the SVD of the 2 × 2 matrix M . Without loss of generality, we assume Σ = diag(σ1 , σ2 ) with σ1 ≥ σ2 > 0, and cos φu − sin φu cos φv − sin φv , V = Rφv = . U = Rφu = sin φu sin φv cos φu cos φv Rφu and Rφv represent the linear transform of rotation (counterclockwise) by angles φu and φv , respectively. ˆ into a triangle centered at the origin. Its effect The mapping x = M ξ maps K can be understood as the composition of three operations (see, e.g., [16]): (1) rotation clockwise by angle φv ; then (2) stretching by factors σ1 and σ2 in ξ and η directions,
23
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
respectively; and finally (3) rotation counterclockwise by angle φu . A circle (centered at (0,0)) in the ξη-plane will be mapped into an ellipse in the xy-plane with two axes σ1 and σ2 , and the longer axis forms an angle φu with the x-axis. See Figure 2. ξ1
σ
φ
ξ
ξ1
2
v
1
σ1
1
φ
u
ˆ to K under the mapping M : R∗ K, ˆ ΣR∗ K, ˆ Rφ ΣR∗ K. ˆ Fig. 2. From K u φv φv φv
We define the aspect ratio of the triangle K as r12 =
(1)
σ1 σ2
and define the orientation of K as the direction of angle φu with the x-axis. It can be seen that K is equilateral if and only if σ1 = σ2 and that K is isosceles if and only if φv = mπ 6 with integer m or σ1 = σ2 . Together with the aspect ratio, φv determines the three internal angles of the triangle. For a fixed aspect ratio, the maximum internal angle of K is a periodic even function of φv with period π3 , and it is decreasing in (0, π6 ). K is an obtuse isosceles triangle when φv = mπ 3 and an acute isosceles triangle when φv = (2m+1)π . If K is an obtuse/acute isosceles, then 6 its orientation is in parallel/perpendicular to its base. Furthermore, let b and h be the length of the base and the height over the base of an isosceles triangle; then √ √ √ σ2 = 2√ r12 = 23 hb , φv = 0 when h ≤ √23 b, σ1 = 33 b, 3 h, σ1 = 23 h, σ2 = 33 b, r12 = √23 hb , φv = π6 when h > 23 b. For general triangles, the aspect ratio r12 can be found as follows. Let σ2 1 σ1 ∈ [1, ∞). + q(K) = 2 σ2 σ1 This quantity can be calculated from the three sides and the area of K as (see formula (10) below) (2)
q(K) =
σ12 + σ22 |1 |2 + |2 |2 + |3 |2 √ = . 2σ1 σ2 4 3|K|
Therefore, r12 = q +
q 2 − 1.
Clearly, the closer r12 is to 1, the closer K is to being equilateral. The same is true for the quantity q(K). Therefore, both r12 and q(K) can be used to measure the closeness of a triangle to being equilateral. Indeed, the reciprocal of the righthand-side expression in (2) was used by Bank and Smith [3] as the “shape regularity
24
WEIMING CAO
quantity” of a triangle. q(K) can also be expressed in terms of the matrix M directly. Note that σ1 and σ2 are eigenvalues of (M T M )1/2 ; it is easy to see that q(K) =
M 2F , 2 det(M )
where · F stands for the Frobenius norm of a matrix. This formula was used by Knupp, Margolin, and Shashkov [10] in certain functionals characterizing the smoothness of the mesh. There are several other ways to define the aspect ratio and the orientation of a triangular element in the finite element analysis. For instance, the aspect ratio is usually defined as the ratio of the length of the longest side over the perpendicular distance from it to the opposite vertex, or as the ratio of the diameter of the triangle over the diameter of the largest inscribed circle in the triangle. The orientation can be defined as the direction of the longest side. It is easy to see that these definitions are equivalent to (1) up to some bounded constants. However, they are not precise enough for describing accurately the behavior of the interpolation errors. 3. Formulas for H 1 - and L2 - norms of the interpolation error. Without loss of generality, assume K is a triangle with its center at the origin and its orientation being the x-axis, i.e., xc = 0, φu = 0. We study the H 1 - and L2 -norms of the error for the linear interpolation of a quadratic function u over K. Since the first order terms of u have no contribution to the error of linear interpolation, we assume in particular that u = 21 x · Hx, where H is the Hessian of u. Since H is a 2 × 2 symmetric matrix, we may decompose it into λ1 RφT h , (3) H = Rφh λ2 where Rφh is the matrix of rotation by an angle φh . We also assume that λ1 ≥ |λ2 |. Other cases can be covered by simply considering the function −u. It is easy to see that the contour lines of u are concentrical ellipses √ (when λ √1 λ2 > 0) or hyperbolas (when λ1 λ2 < 0). Their axes are multiples of 1/ λ1 and 1/ λ2 , and the longer axis (with λ2 ) is of angle φh with the y-axis. We define this direction as the orientation of function u. When φh = π2 , u is oriented along the x-axis, i.e., in the same direction as triangle K. In this case, we say K and u are aligned. When φh = 0, u is oriented along the y-axis, i.e., perpendicular to the orientation of K. See Figure 3.
x1
Fig. 3. Triangle K and the contour lines of a function u. φh is about − π6 in this graph.
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
25
Denote by uI the linear interpolation of u at the three vertices of K. Clearly the norm of the error u − uI depends on the size and the shape of triangle K, as well as the magnitude and the orientation of function u. We derive in this section the exact formulas for the H 1 - and L2 -norms of the interpolation error in terms of these factors. Theorem 3.1. Let K be a triangle oriented along the x-axis. r12 is the aspect ratio of K, r21 = 1/r12 . Let u be a quadratic function. λ1 and λ2 are the eigenvalues of the Hessian matrix of u, and the orientation of u is of angle φh with the y-axis; then √ 3 2 ∇(u − uI )L2 (K) = |K|2 [(r12 + r21 )(λ21 + λ22 ) + (r12 − r21 )(λ21 − λ22 ) cos(2φh )] 36 1
1 − 4λ1 λ2 + ( (λ1 + λ2 )(r12 + r21 ) + (λ1 − λ2 )(r12 − r21 ) + (4) 2 4 · cos(2φh ) )2 [(r12 + r21 + (r12 − r21 ) cos(6φv + 4θ)] ,
ˆ to K, and where φv is the rotation angle in the mapping from K 1 2(λ1 − λ2 ) sin(2φh ) θ = atan (5) . 2 (λ1 + λ2 )(r12 − r21 ) + (λ1 − λ2 )(r12 + r21 ) cos(2φh ) Proof. Step 1. We start with a result established by Bank and Smith in [3]. Let ci (x, y), i = 1, 2, 3, be the barycentrical coordinates of a point (x, y) in K. The side basis functions (taking value 14 at a midpoint) are b1 (x, y) = c2 (x, y) c3 (x, y), b2 (x, y) = c3 (x, y) c1 (x, y), b3 (x, y) = c1 (x, y) c2 (x, y). It is easy to see that 1 u − uI = − (v1 b1 (x, y) + v2 b2 (x, y) + v3 b3 (x, y)), 2
(6) where
vi = i · Hi ,
(7)
i = 1, 2, 3.
Let v = [v1 , v2 , v3 ]T . It is further established in [3] that
1 (8) |∇(u − uI )|2 dxdy = v · Bv, 4 K where
∇bi · ∇bj dxdy K ⎡ |1 |2 + |2 |2 + |3 |2 = 1 ⎣ 48|K| symm.
B=
21 · 2 |1 |2 + |2 |2 + |3 |2
⎤ 21 · 3 ⎦. 22 · 3 |1 |2 + |2 |2 + |3 |2
26
WEIMING CAO
We first derive a formula for the matrix B in terms of the singular values and the rotation angle φv of the mapping M . Note that we assume φu = 0; therefore √ σ1 cos αi ∗ , i = M ei = ΣRφv ei = 3 σ2 sin αi where αi = 2(i − 1)π/3 − φv ,
(9)
i = 1, 2, 3.
With formula (A.1) in the appendix it is easy to verify that (10)
|1 | + |2 | + |3 | = 2
2
2
3σ12
3
2
cos αi +
3σ22
i=1
3
sin2 αi =
i=1
9 2 (σ + σ22 ). 2 1
Let fij =
(11)
σ12
σ12 σ2 cos αi cos αj + 2 2 2 sin αi sin αj . 2 + σ2 σ1 + σ2
Then i · j = 3(σ12 cos αi cos αj + σ22 sin αi sin αj ) = 3(σ12 + σ22 )fij , Note that |K| =
√ 3 3 4 σ1 σ2 .
We may express B as ⎡ 4 √ 2 1 2 3 f12 3 σ1 + σ2 ⎣ 1 B= 24 σ1 σ2 symm.
4 3 f13 4 3 f23
1 ≤ i, j ≤ 3.
⎤ ⎦,
1
and rewrite the norm of the error as √ 2 3 σ1 + σ22 2 (12) ∇(u − uI )L2 (K) = 96 σ1 σ2 3 8 (vi )2 + (f12 v1 v2 + f13 v1 v3 + f23 v2 v3 ) . · 3 i=1 Step 2. We next simplify the terms involving vi . Assume the Hessian matrix H = (hij ), and denote by σ1 ˜ = σ2 h11 σh12 (13) . H 2 h12 σ1 h22 ˜ is a 2 × 2 symmetric matrix, we may decompose it into Since H μ1 ˜ = Rθ RθT , (14) H μ2 where Rθ is the matrix for the counterclockwise rotation by an angle θ, and μ1 and ˜ By the facts that μ2 are the eigenvalues of H. ˜ (R∗ ei ) vi = i · Hi = (Rφ∗ v ei ) · (ΣT HΣ) (Rφ∗ v ei ) = σ1 σ2 (Rφ∗ v ei ) · H φv ∗ ∗ ∗ ∗ = σ1 σ2 (Rθ Rφv ei ) · diag(μ1 , μ2 ) (Rθ Rφv ei )
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
and that Rθ∗ Rφ∗ v ei =
√
3[cos(αi − θ), sin(αi − θ)]T ,
we have vi = 3σ1 σ2 [μ1 cos2 (αi − θ) + μ2 sin2 (αi − θ)] = 32 σ1 σ2 [(μ1 + μ2 ) + (μ1 − μ2 ) cos 2(αi − θ)]. Therefore it follows from formulas (A.2) and (A.3) that 9 3 2 2 2 (vi ) = (σ1 σ2 ) 3(μ1 + μ2 ) + (μ1 − μ2 ) 4 2 i=1 81 2 2 2 2 (σ1 σ2 ) μ1 + μ2 + μ1 μ2 . = 8 3
3
(15)
2
Expand vi vj as follows: (16)
vi vj =
9 (σ1 σ2 )2 [(μ1 + μ2 )2 + (μ21 − μ22 )(cos 2(αi − θ) + cos 2(αj − θ)) 4 + (μ1 − μ2 )2 cos 2(αi − θ) cos 2(αj − θ)].
Let β = φv + θ. Then αi − θ = 2(i − 1)π/3 − β, and ⎧ ⎨ cos(2β + π3 ) cos(2β − π3 ) cos 2(αi − θ) + cos 2(αj − θ) = ⎩ cos(2β + π)
for (i, j) = (1, 2), for (i, j) = (1, 3), for (i, j) = (2, 3)
and
⎧ 1 ⎪ ⎨ −4 − − 14 − cos 2(αi − θ) · cos 2(αj − θ) = ⎪ ⎩ −1 − 4
1 2 1 2 1 2
cos(4β − π3 ) cos(4β + π3 ) cos(4β + π)
for (i, j) = (1, 2), for (i, j) = (1, 3), for (i, j) = (2, 3).
Substituting the above formulas into (17), we have
(17)
f12 v1 v2 + f13 v1 v3 + f23 v2 v3 = 94 (σ1 σ2 )2 {[(μ1 + μ2 )2 − 14 (μ1 − μ2 )2 ] (f12 + f13 + f23 ) + (μ21 − μ22 ) [f12 cos(2β + π3 ) + f13 cos(2β − π3 ) + f23 cos(2β + π)] − 12 (μ1 − μ2 )2 [f12 cos(4β − π3 ) + f13 cos(4β + π3 ) + f23 cos(4β + π)] }.
By using the definition (11) for fij , we have from formula (A.4) that 3 f12 + f13 + f23 = − ; 4
(18)
from formulas (A.5) and (A.6) that f12 cos(2β + π3 ) + f13 cos(2β − π3 ) + f23 cos(2β + π) (19)
= − 34
σ12 − σ22 cos(2θ); σ12 + σ22
27
28
WEIMING CAO
and from formulas (A.7) and (A.8) that f12 cos(4β − π3 ) + f13 cos(4β + π3 ) + f23 cos(4β + π) (20)
= − 34
σ12 − σ22 cos(6φv + 4θ). σ12 + σ22
Put (17)–(20) into (13); we may express the norm of the error as √ 3 3 2 σ12 − σ22 2 2 2 2 2 2 ∇(u − uI )|L2 (K) = (σ + σ2 )σ1 σ2 μ1 + μ2 − 2 (μ1 − μ2 ) cos(2θ) 64 1 σ1 + σ22 1 σ12 − σ22 2 + (μ1 − μ2 ) 1 + 2 (21) cos(6φv + 4θ) . 2 σ1 + σ22 Step 3. Finally we express μ1 , μ2 , and θ in terms of the eigenvalues of the Hessian H and the aspect ratio r12 of K. According to the eigendecomposition (3) of H, λ1 cos2 φh + λ2 sin2 φh (λ1 − λ2 ) sin φh cos φh . H= (λ1 − λ2 ) sin φh cos φh λ1 sin2 φh + λ2 cos2 φh Therefore
˜ ij ) = ˜ = (h H
D sin 2φh r12 (A + D cos 2φh ) D sin 2φh r21 (A − D cos 2φh )
with A=
1 (λ1 + λ2 ), 2
D=
1 (λ1 − λ2 ). 2
Recall the eigendecomposition (14), it is not difficult to establish for the eigenvalues ˜ that μ1 and μ2 of H μ21 + μ22
˜ 11 )2 + (h ˜ 22 )2 + 2(h ˜ 12 )2 = (h 2 2 2 2 2 2 2 = (A + D cos (2φh ))(r12 + r21 ) + 2AD(r12 − r21 ) cos(2φh ) 2 2 + 2D sin (2φh ),
(μ21 − μ22 ) cos 2θ
˜ 11 )2 − (h ˜ 22 )2 = (h 2 2 2 2 2 2 = (A + D cos2 (2φh ))(r12 − r21 ) + 2AD(r12 + r21 ) cos(2φh ),
(μ1 − μ2 )2
˜ 11 − h ˜ 22 )2 + 4(h ˜ 12 )2 = (h 2 2 2 2 = 4D + A (r12 − r21 )2 + 2AD(r12 − r21 ) cos(2φh ) 2 2 2 +D (r12 − r21 ) cos (2φh ) 2 2 − r21 ) cos(2φh ) = 4(D2 − A2 ) + A2 (r12 + r21 )2 + 2AD(r12 2 2 2 +D (r12 − r21 ) cos (2φh ) = − 4λ1 λ2 + [A(r12 + r21 ) + D(r12 − r21 ) cos(2φh )]2
and that tan(2θ) =
˜ 12 2h 2D sin 2φh = . ˜ 22 ˜ 11 − h A(r12 − r21 ) + D(r12 + r21 ) cos 2φh h
Substitute the above formulas into the right-hand side of (21) and simplify, and we obtain the formula (5).
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
29
We next derive the formula for the L2 -norm of the linear interpolation error. Theorem 3.2. Let K be a triangle oriented along the x-axis. r12 is the aspect ratio of K, r21 = 1/r12 . Let u be a quadratic function. λ1 and λ2 are the eigenvalues of its Hessian, and the orientation of u is of angle φh with the y-axis. Then 16 |K|3 2 − λ1 λ2 + [(λ1 + λ2 )(r12 + r21 ) (22) u − uI L2 (K) = 160 9 + (λ1 − λ2 )(r12 − r21 ) cos(2φh )]2 . Proof. It follows from (6) that
1 (23) |u − uI |2 dxdy = v · B0 v, 4 K where
bi · bj dxdy .
B0 = K
An elementary calculation of the integrals leads to ⎡ ⎤ 2 1 1 |K| ⎣ 1 2 1 ⎦. B0 = 180 1 1 2 Therefore
(24)
|K| |u − uI |2 dxdy = 360 K
3
(vi )2 + v1 v2 + v1 v3 + v2 v3
,
i=1
where v1 , v2 , and v3 are defined as in (7). We have similar to (17) that 81 10 2 2 2 (σ1 σ2 ) μ1 + μ2 + μ1 μ2 , v1 v2 + v1 v3 + v2 v3 = 16 3 which, together with (15), yields
|K|3 4 (μ1 + μ2 )2 − μ1 μ2 . (25) |u − uI |2 dxdy = 40 9 K ˜ = det(H) = λ1 λ2 and Finally, by μ1 μ2 = det(H) ˜ 11 + h ˜ 22 = A(r12 + r21 ) + D(r12 − r21 ) cos(2φh ), μ1 + μ2 = h we prove the conclusion of this theorem. 4. Discussion about some special cases. Denote by T1 , T2 , and T3 the terms in the square brackets on the right-hand side of (5), i.e., (26) (27) (28)
T1 = (r12 + r21 )(λ21 + λ22 ) + (r12 − r21 )(λ21 − λ22 ) cos(2φh ), 1 2 T2 = −4λ1 λ2 + ((λ1 + λ2 )(r12 + r21 ) + (λ1 − λ2 )(r12 − r21 ) cos(2φh )) , 4 T3 = r12 + r21 + (r12 − r21 ) cos(6φv + 4θ).
30
WEIMING CAO
Recall that we assume r12 = σ1 /σ2 ≥ 1 and λ1 ≥ |λ2 |; therefore all the three terms above are nonnegative, and we may write √ 3 1 ∇(u − uI )2L2 (K) = (29) |K|2 T1 + T2 · T3 . 36 2 We may also restrict 0 ≤ φv < π3 and 0 ≤ φh < π. Other cases are covered by the symmetry and periodicity of the error norms. For a fixed aspect ratio, the extremum values of T1 , T2 , and T3 are min (T1 ) = 2(r21 λ21 + r12 λ22 ) when φh = π2 , (30) T1 = 2 2 when φh = 0, max (T1 ) = 2(r12 λ1 + r21 λ2 ) 2 when φh = π2 , min (T2 ) = (r21 λ1 − r12 λ2 ) T2 = (31) 2 when φh = 0, max (T2 ) = (r12 λ1 − r21 λ2 ) when 6φv + 4θ = (2m + 1)π, min (T3 ) = 2r21 T3 = (32) when 6φv + 4θ = 2mπ. max (T3 ) = 2r12 We give more details about these cases. Case 1. φh = π2 and φv = π6 , i.e., K is an acute isosceles triangle aligned with u. In this case, θ = 0 by (5). Therefore, T1 , T2 , and T3 all take their minimum values, which yields √ 3 2 |K|2 [ 2(r21 λ21 + r12 λ22 ) + r21 (r21 λ1 − r12 λ2 )2 ]. (33) ∇(u − uI )L2 (K) = 36 This is the smallest H 1 -seminorm of the interpolation error for all the triangles with a fixed aspect ratio and a fixed area. We may consider the optimal choice of the aspect ratio in this case. Let 2 f (r) = 2 1r λ21 + rλ22 + 1r 1r λ1 − rλ2 = 3λ22 r + 2(2λ21 − 4λ1 λ2 ) 1r + λ21 13 . r (1)
(1)
It can be seen that f is decreasing in [1, r∗ ) and increasing in (r∗ , ∞), where ! λ (λ − λ ) + λ (λ − λ )2 + 9λ2 1 1 2 1 1 2 2 (1) (34) . r∗ = 2 3λ2 (1)
Therefore r12 = r∗ is the best aspect ratio. Moreover, since T1 , T2 , and T3 all take minimum values in this case, we conclude that the acute isosceles triangle, which is (1) aligned with u and of the aspect ratio r∗ , is the best one that produces the smallest H 1 -seminorm of the interpolation error among all the triangles with the same area! (1) It is noted that for u with λ1 = λ2 (isotropic u), the best aspect ratio ∗ =1 ! is r√ (1) (with isotropic K). For u with λ2 = −λ1 , the best aspect ratio is r∗ = (2 + 13)/3 ≈ 1.367. When λ1 >> |λ2 |, we have " # # # # # λ1 # 2 ## λ1 ## (1) ≈ 0.816 ## ## , r∗ ≈ # # 3 λ2 λ2
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
31
with which the interpolation error is of the smallest possible magnitude √ 2 2 |K|2 |λ1 λ2 | ≈ 0.235|K|2 |λ1 λ2 |. (35) ∇(u − uI )L2 (K) ≈ 6 We may also compare the above optimal choice of the aspect ratio with some intuitive choices. For instance, if we choose r12 = 1, then √ 3 2 |K|2 (3(λ21 + λ22 ) − 2λ1 λ2 ), (36) ∇(u − uI )L2 (K) ≈ 36 which is between 0.128|K|2 λ21 and 0.193|K|2 λ21 . If we choose r12 = | λ1 |, then in the λ2 case λ1 >> |λ2 |, √ 5 3 2 |K|2 |λ1 λ2 | ≈ 0.240|K|2 |λ1 λ2 |; ∇(u − uI )L2 (K) ≈ (37) 36 if we choose r12 = |λ1 /λ2 |, then we have √ ! ! 3 2 ∇(u − uI )L2 (K) = |K|2 |λ31 λ2 | ≈ 0.0962|K|2 |λ31 λ2 |. (38) 18 Note that when u and K are aligned, the inverse mapping ξ = M −1 x with r12 = |λ1 /λ2 | transforms u(x, y) into const.(ξ 2 ± η 2 ) and K into an equilateral triangle. It was shown by D’Azevedo and Simpson [6] and Rippa [13] that triangles with such an aspect ratio lead to the smallest maximum-norm and Lp -norm of the interpolation error among all triangles of the same area. However, this aspect ratio is not the optimal for the interpolation error in the sense of H 1 -seminorm. Case 2. φh = π2 and φv = 0, i.e., K is an obtuse isosceles triangle aligned with u. We have in this case √ 3 2 |K|2 [ 2(r21 λ21 + r12 λ22 ) + r12 (r21 λ1 − r12 λ2 )2 ]. (39) ∇(u − uI )L2 (K) = 36 This error formula differs from that of Case 1 only in the coefficient of the second term. To study the optimal choice of the aspect ratio in this case, let 2 g(r) = 2 1r λ21 + rλ22 + r 1r λ1 − rλ2 = 3λ21 1r + 2(λ22 − λ1 λ2 )r + λ22 r3 . It can be shown that the minimum of g is attained at ! λ (λ − λ ) + |λ | (λ − λ )2 + 9λ2 2 1 2 2 1 2 2 (2) r∗ = (40) . 2 3λ2 (2)
When λ1 = λ2! , the best aspect ratio is r∗ = 1. When λ2 = −λ2 , the best aspect √ (1) ratio is r∗ = (2 + 13)/3 ≈ 1.367. When |λ1 | >> |λ2 |, we have ⎧ "√ " ⎪ λ 10 + 1 1 ⎪ | | ≈ 1.178 | λ1 | when λ1 λ2 > 0, ⎨ 3 λ2 λ (2) " " 2 r∗ ≈ √ ⎪ ⎪ 10 − 1 | λ1 | ≈ 0.849 | λ1 | ⎩ when λ1 λ2 < 0. 3 λ2 λ2
32
WEIMING CAO
With this choice of aspect ratio, the interpolation error is of the magnitude when λ1 λ2 > 0, 0.0878|K|2 |λ31 λ2 | 2 ∇(u − uI )L2 (K) ≈ (41) 2 3 when λ1 λ2 < 0. 0.281|K| |λ1 λ2 | We may compare this case (obtuse isosceles) with the previous one (acute isosceles). When the triangle is aligned with u and the aspect ratio r12 ≈ |λ1 /λ2 |, the interpolation error is nearly the same for both acute and obtuse triangles if λ1 λ2 > 0. For u with λ1 λ2 < 0, ∇(u − uI )L2 (K) is about 1.7 times smaller with the acute isosceles triangle than that with the obtuse one. Furthermore, because for general triangles with arbitrary φv , the error norm ∇(u − uI )L2 (K) is between those of Case 1 and Case 2, we may conclude that the H 1 -seminorm of the interpolation error is not sensitive to the maximum angle of the triangle, as long as the triangle is aligned with the function u and the aspect ratio r12 ≈ |λ1 /λ2 | is used. Case 3. φh = 0. In this case the orientation of u is perpendicular to the triangle. We have √ 3 2 |K|2 [2(r12 λ21 + r21 λ22 ) + (r12 λ1 − r21 λ2 )2 ∇(u − uI )L2 (K) = 18 (42) · (r12 + r21 + (r12 − r21 ) cos(6φv ))] For given |K| and φv , it can be shown the error norm is an increasing function of r12 . Therefore, the best aspect ratio is r12 = 1. This is not surprising, since when the triangle is perpendicular to the orientation of u, increasing the aspect ratio leads to lower resolution in the needy direction and larger interpolation error, no matter what value φv takes. We should emphasize that in the above discussion on the optimal aspect ratios the area and the orientation of the triangle are assumed fixed. In adaptive mesh generation there may be some other constraints on the triangles in a mesh, e.g., one side of a triangle or the vertices of the triangles are fixed. Optimal choices of the aspect ratio subject to those constraints can be quite different from the above values. 5. General discussion about orientation, aspect ratio, and angle condition. In this section we present a general discussion about the effects of the orientation, aspect ratio, and the angle φv . We first study the H 1 -seminorm of the interpolation error. Mesh alignment. Given a triangle K (fixed |K|, r12 , and φv ), it is commonly believed that the interpolation error is the smallest when the triangle is aligned with the function, i.e., when φh = π2 . We present a justification of this viewpoint. We consider the case where both the solutions and the triangles are highly anisotropic, i.e., we assume λ1 >> |λ2 | and r12 >> 1. In this case we have for the three terms (26)–(28) in the error norm ∇(u − uI )2L2 (K) that T1 = r12 λ21 (1 + cos 2φh ) + LOT, T2 = 14 (r12 λ1 )2 (1 + cos 2φh )2 + LOT, T3 = r12 (1 + cos(6φv + 4θ)) + LOT and for the angle θ in (5) that tan(2θ) ≈
2λ1 sin 2φh = 2r21 tan φh , λ1 r12 (1 + cos 2φh )
33
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
where the lower order terms LOT include r12 λ22 , r12 λ21 , r21 λ22 , etc. In this case the error norm is dominated by ∇(u − uI )2L2 (K) ≈
√
3 2 36 |K|
3 2 [r12 λ21 + 14 r12 λ1 (1 + cos 2φh )(1 + cos(6φv + 4θ))] · (1 + cos 2φh ) + LOT.
Note that by (5) θ is an increasing function of φh , and cos(2φh ) is a decreasing function of φh in (0, π2 ). Therefore, we may conclude that the H 1 -seminorm of the error is almost a monotonically decreasing function of the angle φh in (0, π2 ). We plot in Figure 4 the graph of ∇(u − uI )2L2 (K) versus φh for various λ1 , λ2 , φv and r12 . It is observed that for all aspect ratios and φv , the error is smallest when φh is near π2 . Also seen in Figure 4 is a small peak of the error norm around φh = π2 . This is because when φh is very close to π2 , the leading terms in T1 and T2 go to 0 and the lower order terms become the major components in the error. 4
3
10
10
2
(K)
L
||∇(u-u )||2
v
2
10
φv=0 φv=π/8 φ =π/6
2
10
||∇(u-u L)||(K)
φv=0 φv=π/8 φ =π/6
3
10
v
1
10
1
10
0
10
0
0
0.2 0.4 0.6 0.8 with λ1=10, λ2=1, r12=8
10
1
10
0.2 0.4 0.6 0.8 with λ1=10, λ2=1, r12=4
1
8
10
10
φv=0 φv=π/8 φv=π/6
φv=0 φv=π/8 φv=π/6
6
10 2
||∇(u-u L)||(K)
2
||∇(u-u L)||(K)
0
5
10
4
10
2
10 0
10
0
0
0.2 0.4 0.6 0.8 with λ1=100, λ2=1, r12=80
Fig. 4. u − uI 2 2
L (K)
1
10
0
0.2 0.4 0.6 0.8 with λ1=100, λ2=1, r12=40
versus the alignment angle φh /π with |K| =
3
√ 3 4
1
fixed.
Aspect ratio. We are interested here in the question of, given a function u and a fixed triangle orientation φh , how one chooses the aspect ratio of the triangle. We plot in Figure 5 the graph of ∇(u − uI )L2 (K) versus r12 with various φh . It is observed that when φh = π2 (K aligned with u), the best aspect ratio r12 is in between (2) (1) r∗ = |λ1 /λ2 | and r∗ ≈ 0.8|λ1 /λ2 |, while when φh ≈ 0 (K perpendicular to u), the best aspect ratio is 1. For other φv and φh , the optimal choice of r12 lies between alignment is ensured, the 1 and |λ1 /λ2 |. Therefore, in practice, when a good mesh proper aspect ratio should be between the magnitude of |λ1 /λ2 | and 0.8|λ1 /λ2 |, with the lower end taken for the mesh with mostly obtuse triangles and the upper
34
WEIMING CAO 4
5
10
10
φv=0 φv=π/8 φ =π/6
2
||∇(u-u L)||(K)
v
2
10
v
3
10
1
2
10
10
0
10
φv=0 φv=π/8 φ =π/6
4
10 ||∇(u-u L)||2(K)
3
10
1
0
10
20 with φ =π/2
30
40
10
0
10
h
6
6
φ =0 v φv=π/8 φv=π/6
4
φ =0 v φv=π/8 φv=π/6
4
10
||∇(u-u L)||2(K)
2
40
10
10
||∇(u-u L)||(K)
30
h
10
2
10
2
10
0
10
20 with φ =π/3
0
0
10
20 with φh=π/6
Fig. 5. ∇(u − uI )2 2
L (K)
30
40
10
0
10
20 with φh=0
30
versus the aspect ratio r12 with λ1 = 10, λ2 = 1, and |K| =
40
3
√ 3 . 4
end taken for the mesh with mostly acute triangles. The less the mesh alignment, the smaller r12 should be. Maximum angle condition. It is well known that a good mesh for finite element approximation should respect the so-called maximum angle condition, i.e., the maximum angles of the triangular elements should be bounded from π. In [2], Babu˘ska and Aziz gave an example showing that when the maximum angle condition is violated, the H 1 -seminorm of the error of linear interpolation can go to infinity, although the element area approaches 0. Their example is in our nota√ tion the case with λ1 = 1, λ2 = 0, |K| = 2 , r12 = 2 3 , φh = 0, φv = 0, and thus 9 3 ∇(u − uI )2L2 (K) = 256 ( + 8
). We discuss here the impact of the maximum angle condition on the accuracy of linear interpolation. Taking as an example with λ1 = 100 and λ2 = 1, we plot in Figure 6 the graph of ∇(u − uI )2L2 (K) versus φv with given aspect ratios and alignment directions. It is easy to see the following: (a) When λ1 /λ2 and r12 are high, triangles with φv ≈ π6 or π2 produce the smallest error, and triangles with φv ≈ 0 or π3 produce the largest error, whether they are aligned with u or not. This can be justified by looking into formula (5) for the error norm. When λ1 >> |λ2 | and r12 >> 1, we have θ ≈ 0. Hence, the term T3 , and therefore ∇(u − uI )L2 (K) , attains its minimum when φv ≈ π6 or π2 . Note that the maximum internal angle of K is an even periodic function of φv with the period π/3. It is decreasing in (0, π/6). Hence φv = 0, π/3 corresponds to the case where the maximum internal angle of K is the largest among all the triangles of the same aspect ratio, and φv = π/6, π/2 corresponds to the case where
35
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
the maximum internal angle is the smallest among all triangles. Therefore, in this case it is important to use triangles with smaller maximum internal angles, regardless of whether or not the mesh alignment is satisfied. (b) When the triangle is aligned with u and the aspect ratio r12 ≈ |λ1 /λ2 |, the error is insensitive to φv and is therefore insensitive to the maximum internal angle of the triangle. In particular, if φh = π2 , λ1 λ2 > 0, and r12 = λ1 /λ2 , then T2 = 0, and by (29) ∇(u − uI )L2 (K) is independent of φv . See the lower left graph of Figure 6 and the analysis in Cases 1 and 2 in section 4. (c) When r12 ≈ 1, the best φv may be different from π6 and π2 . However, in this case the difference between the maximum and the minimum of the error is not significant. 10
8
10
φ =0 h φh=π/6 φh=π/3 φh=π/2
5
10
φ =0 h φh=π/6 φh=π/3 φh=π/2
6
10 2
||∇(u-u L)||(K)
2
||∇(u-u L)||(K)
10
4
10
2
10 0
10
0
0
0.2 0.4 0.6 0.8 with λ =100, λ =1, r =80 1
2
10
1
6
2
1
12
4
10
2
2
4
10
φh=0 φh=π/6 φ =π/3 h φh=π/2
||∇(u-u L)||(K)
φh=0 φh=π/6 φ =π/3 h φh=π/2
5
10 L
0.2 0.4 0.6 0.8 with λ =100, λ =1, r =20 1
10
||∇(u-u )||(K)
0
12
3
10
2
10
3
0
0.2 0.4 0.6 0.8 with λ1=100, λ2=1, r12=10
Fig. 6. ∇(u − uI )2 2
L (K)
1
10
0
0.2 0.4 0.6 0.8 with λ1=100, λ2=1, r12=1.5
versus the angle φv /π with |K| =
3
√ 3 4
1
fixed.
In summary, we conclude that ifthe triangle is aligned with the function and the aspect ratio is of an magnitude |λ1 /λ2 | or smaller, then the maximum angle condition is not essential to the H 1 -seminorm of the interpolation error. Otherwise, the maximum angle is critical to the magnitude of the error. L2 -norm. From (23) we can conclude the following for the L2 -norm of the linear interpolation error: (1) u−uI L2 (K) does not depend on the angle φv . Therefore, for a given function, the L2 -norm of its linear interpolation error depends only on the area, aspect ratio, and orientation of the triangle. Maximum/minimum angle conditions are irrelevant to the L2 -norm of the interpolation error. (2) u − uI L2 (K) is a π-periodic even function of φh . It is decreasing in φh in
36
WEIMING CAO
(0, π2 ). When φh =
π 2
(K aligned with u),
u − uI 2L2 (K) =
|K|3 40
4 (λ1 r21 + λ2 r12 )2 − λ1 λ2 9
is the smallest, and when φh = 0 (K perpendicular to u) |K|3 4 2 2 (λ1 r12 + λ2 r21 ) − λ1 λ2 u − uI L2 (K) = 40 9 is the largest. (3) With the same orientation, it can be shown that when φh is in (0, π4 ) ∪ ( 3π 4 , π) (i.e., K is aligned more to the perpendicular direction of u), the best aspect ratio is r12 = 1; when π4 ≤ φh ≤ 3π 4 (i.e., K is aligned more to the orientation of u), the best aspect ratio is $# # # λ1 + λ2 − (λ1 − λ2 ) cos(2φh ) # #. (43) r12 = ## λ1 + λ2 + (λ1 − λ2 ) cos(2φh ) # We plot in Figure 7 the L2 -norm of the error versus the aspect ratio for the case of √ λ1 = 400, λ2 = 1, and |K| = 3 4 3 . 8
10
φ =0 h φh=π/6 φ =π/4 h φ =π/3 h φ =2π/5 h φh=π/2
7
10
6
10
5
||u-u ||2
L (K)
10
4
10
3
10
2
10
1
10
Fig. 7. u − uI 2 2
0
L (K)
5
10
15
20 25 aspect ratio r12
30
35
40
45
versus the aspect ratio r12 with λ1 = 400, λ2 = 1 and |K| =
3
√ 3 . 4
(4) For a fixed triangle area, the minimum L2 -norm of the linear interpolation π error is attained at φh = 2 and r12 = |λ1 /λ2 |. The minimum value is 8 3 when λ1 λ2 > 0, 2 90 |K| |λ1 λ2 | u − uI L2 (K) = 1 3 |K| |λ λ | when λ1 λ2 < 0. 1 2 90 6. An example of linear interpolation on different meshes. In this section, we present the results of piecewise linear interpolation of a quadratic function on different types of mesh. We choose u(x, y) = x2 + 100y 2 ,
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
(a)
(b)
37
(c)
Fig. 8. Three types of mesh used for piecewise linear interpolation.
which corresponds to λ1 = 100, λ2 = 1, and the orientation along the x-axis. We consider the following three types of mesh on [0, 1] × [0, 1]: (a) The first type of mesh is as shown in Figure 8(a). There√ are M and N equal partitions in x and y directions, respectively. When M/N < 2/ 3, all the triangles in this mesh, except those on the top and the bottom boundaries, are aligned exactly with the x-axis. Their aspect ratio is 2N r12 = √ 3M and their angle φv is π6 . (b) The second type of mesh is as shown in Figure 8(b)√with M and N equal partitions in x and y directions, respectively. When M/N < 3/2, all the triangles, except those on the left and the right boundaries, are aligned with the x-axis. Their aspect ratio is √ 3N r12 = 2M and their angle φv is 0. (c) The third type is the unstructured meshes generated in the following way. First we create a Delaunay triangulation by using the package Triangle [14] on a rectangle [0, 1] × [0, rs ], where rs ≥ 1 is a real number. Then we compress the rectangle in y direction by a factor rs . The desired mesh on the unit square is the image of the Delaunay triangulation under the compression. See Figure 8(c) for a typical mesh of this type. By Delaunay triangulation the minimum internal angle of the triangles in the mesh is maximized, and all the triangles in the mesh over [0, 1] × [0, rs ] are approximately of the same size and close to unilateral. Therefore, most triangles in the mesh over the unit square are roughly aligned with the x-axis (when rs > 1) and are of the aspect ratio r12 ≈ rs . However, in this type of mesh the angle φv varies for different triangles (approximately uniformly distributed between 0 and π/3). We are interested in the accuracy of the linear interpolation of u on different types of mesh and with various aspect ratios, in particular with the following three (1) ratios: (1) r∗ = 81.6, which is the best aspect ratio (for H 1 -seminorm) calculated
38
WEIMING CAO (2)
according to (34) for the case of acute isosceles triangles; (2) r∗ = 11.79, which is the best aspect ratio (in H 1 -seminorm) calculated according to (40) for the case ! (∞)
of obtuse isosceles triangles; and (3) r∗ = | λλ12 | = 10, which is the best aspect ratio for the Lp -norm, 1 ≤ p ≤ ∞. We also report the results with the aspect ratios r12 = λ1 /λ2 = 10 and r12 = 1. We list in Table 1 the various norms of the interpolation error with different choices of M , N , and rs . Note that for all the meshes, the total number of elements is around 4000 and the total number of nodes is around 2200. The smallest H 1 seminorm for type (a) meshes is with r12 = 92.37; the smallest H 1 -seminorm for type (b) meshes is with r12 = 11.54. For type (c) meshes, the smallest H 1 -seminorm is (2) achieved at rs = 11.78 ≈ r∗ . This is because for this type of mesh, the error from the obtuse triangles (with φv ≈ 0) dominates in the global error. Therefore, the aspect (2) ratio close to r∗ is the best choice for the global error norm. For the L2 -, L1 -, and maximum norms, the smallest interpolation error is obtained (∞) with r12 ≈ r∗ for all the three meshes. In summary, we conclude that when the mesh is in good alignment with the solution and most of the triangles are acute isosceles, the aspect ratio should be (1) chosen around r∗ ≈ 0.8| λλ12 |. If the mesh is in good alignment but with mostly obtuse triangles, or with varied maximum internal angles, then the aspect ratio should be (2) chosen around r∗ , which is approximately 1.178 |λ1 /λ2 | for the case λ1 λ1 > 0, and 0.849 |λ1 /λ2 | for the case λ1 λ1 < 0. Table 1 The L1 - and L2 -norms, H 1 -seminorm, and maximum norm of the interpolation error over the entire domain. (∗1) and (∗2) indicate the aspect ratio close to the best values for the H 1 seminorm in each case. M ×N
Node #
Elem #
r12
u−uI L1
u−uI L2
|u−uI |H 1
u−uI ∞
Type (a) meshes 4×500
2507
4004
144.34
5.23e−3
6.18e−3
8.28e−2
7.82e−3
5×400
2409
4005
92.37(∗1)
3.37e−3
3.97e−3
7.68e−2
5.01e−3
6×333
2341
4002
64.08
2.37e−3
2.78e−3
7.76e−2
3.50e−3
14×143
2167
4018
11.79
7.30e−4
8.01e−4
1.43e−1
9.74e−4
15×133
2152
4005
10.24
7.22e−4
7.92e−4
1.54e−1
9.64e−4
48×42
2131
4080
1.01
3.55e−3
4.20e−3
5.88e−1
7.08e−3
Type (b) meshes 3×572
2578
4004
165.11
6.19e−3
7.60e−3
7.22e+0
1.38e−2
4×445
2543
4005
96.34
3.62e−3
4.39e−3
3.21e+0
7.81e−3
5×364
2372
4004
63.04
2.39e−3
2.87e−3
1.68e+0
5.00e−3
12×160
2173
4000
11.54(∗2)
7.47e−4
8.23e−4
1.53e−1
1.01e−3
13×149
2175
4023
9.92
7.35e−4
8.08e−4
1.60e−1
9.93e−4
41×49
2125
4067
1.03
3.50e−3
4.13e−3
5.87e−1
5.22e−3
Type (c) meshes 2474
4101
100
4.24e−3
5.34e−3
2.95e+0
1.91e−2
2426
4079
81.6
3.55e−3
4.57e−3
2.48e+0
1.67e−2
2163
4036
11.78
8.38e−4
9.62e−4
1.84e−1
2.42e−3
2127
4005
10
8.34e−4
9.53e−4
1.87e−1
3.05e−3
2092
4036
1
4.15e−3
5.21e−3
6.75e−1
2.01e−2
LINEAR INTERPOLATION ERROR AND TRIANGLE GEOMETRY
39
Appendix. We list here some basic trigonometric formulas used to derive the results in the previous sections. Let α be any real number. Denote by αi = 2(i − 1)π/3 − α, i = 1, 2, 3.
(A.2)
3 , 2 cos(2α1 ) + cos(2α2 ) + cos(2α3 ) = 0,
(A.3)
cos2 (2α1 ) + cos2 (2α2 ) + cos2 (2α3 ) =
(A.4)
3 cos(α1 ) cos(α2 ) + cos(α1 ) cos(α3 ) + cos(α2 ) cos(α3 ) = − . 4
(A.1)
cos2 (α1 ) + cos2 (α2 ) + cos2 (α3 ) =
3 , 2
Similar relations hold for sine functions, too. Let θ be any real number, and let β = α + θ. Then π π + cos(α1 ) cos(α3 ) cos 2β − (A.5) cos(α1 ) cos(α2 ) cos 2β + 3 3 3 + cos(α2 ) cos(α3 ) cos(2β + π) = − cos(2θ), 4 π π + sin(α1 ) sin(α3 ) cos 2β − (A.6) sin(α1 ) sin(α2 ) cos 2β + 3 3 3 + sin(α2 ) sin(α3 ) cos(2β + π) = cos(2θ), 4 π π + cos(α1 ) cos(α3 ) cos 4β + (A.7) cos(α1 ) cos(α2 ) cos 4β − 3 3 3 + cos(α2 ) cos(α3 ) cos(4β + π) = − cos(6α + 4θ), 4 π π + sin(α1 ) sin(α3 ) cos 4β + (A.8) sin(α1 ) sin(α2 ) cos 4β − 3 3 3 + sin(α2 ) sin(α3 ) cos(4β + π) = cos(6α + 4θ). 4 REFERENCES [1] T. Apel, Anisotropic Finite Elements: Local Estimates and Applications, Adv. Numer. Math., Teubner, Stuttgart, 1999. [2] I. Babu˘ ska and A. K. Aziz, On the angle condition in the finite element method, SIAM J. Numer. Anal., 13 (1976), pp. 214–226. [3] R. E. Bank and R. K. Smith, Mesh smoothing using a posteriori error estimates, SIAM J. Numer. Anal., 34 (1997), pp. 979–997. [4] M. Berzins, A solution-based triangular and tetrahedral mesh quality indicator, SIAM J. Sci. Comput., 19 (1998), pp. 2051–2060. ´mal, Triangular elements in the finite element method, Math. Comp., [5] J. Bramble and M. Zla 24 (1970), pp. 809–820. [6] E. F. D’Azevedo and R. B. Simpson, On optimal interpolation triangle incidences, SIAM J. Sci. Statist. Comput., 10 (1989), pp. 1063–1075. [7] L. Formaggia and S. Perotto, New anisotropic a priori error estimates, Numer. Math., 89 (2001), pp. 641–667. [8] W. Huang, Measuring mesh qualities and application to variational mesh adaptation, SIAM J. Sci. Comput., 26 (2005), pp. 1643–1666. [9] W. Huang and W. Sun, Variational mesh adaptation II: Error estimates and monitor functions, J. Comput. Phys., 184 (2003), pp. 619–648. [10] P. Knupp, L. G. Margolin, and M. Shashkov, Reference Jacobian optimization-based rezone strategies for arbitrary Lagrange Eulerian methods, J. Comput. Phys., 176 (2002), pp. 93– 128.
40
WEIMING CAO
[11] G. Kunert, A Posteriori Error Estimation for Anisotropic Tetrahedral and Triangular Finite Element Meshes, Ph.D. dissertation, TU Chemnitz, Chemnitz, Germany, 1999. [12] E. J. Nadler, Piecewise Linear Approximation on Triangulations of a Planar Region, Ph.D. thesis, Division of Applied Mathematics, Brown University, Providence, RI, 1985. [13] S. Rippa, Long and thin triangles can be good for linear interpolation, SIAM J. Numer. Anal., 29 (1992), pp. 257–270. [14] J. R. Shewchuk, Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator, Technical report, School of Computer Science, Carnegie Mellon University, Pittsburgh, PA, 1996. [15] R. B. Simpson, Anisotropic mesh transformations and optimal error control, Appl. Numer. Math., 14 (1994), pp. 183–198. [16] L. N. Trefethen and D. Bau, Numerical Linear Algebra, SIAM, Philadelphia, 1997.
