Anisotropic Adaptive Finite Elements for an Elliptic Problem with Strongly Varying Diffusion Coefficient

1 Introduction

Adaptive meshes have proved to be very efficient for solving PDEs at reduced computational cost given a prescribed level of accuracy. In the framework of finite elements methods, the adaptive criteria is often based on a posteriori error estimates, which were first derived for isotropic meshes [9, 38, 31, 16, 36, 7, 1] and more recently for anisotropic meshes [23, 30, 39, 26, 15, 20, 21, 19, 34, 33, 27, 11, 10], that is to say meshes with large aspect ratio. Anisotropic finite elements have demonstrated their efficiency when boundary or internal layers are involved, for instance in computational fluid dynamics [3, 4].

Here we focus on the elliptic problem - div ⁡ ( μ ⁢ ∇ ⁡ u ) = f , where μ > 0 is smooth but strongly varying. We will derive anisotropic a posteriori error estimates and show that the effectivity index (the ratio between the estimated error and the true error) is bounded above and below by two constants independent of the data 𝑓, 𝜇, the mesh size and aspect ratio, up to higher order terms. These theoretical findings will be confirmed when using adapted meshes with large aspect ratio.

The outline of the paper is the following. In Section 2, we present the model equation and the numerical method. In Section 3, the anisotropic setting is introduced. In Section 4, we introduce the error estimator and state Theorem 1 which shows the equivalence between the true and the estimated error. In Section 5, numerical experiments with non-adapted meshes are performed. In Section 6, an adaptive algorithm is introduced, and its efficiency is demonstrated in both 2D and 3D situations. Finally, in Appendix A, the proof of Theorem 1 is presented.

2 Problem Statement and Numerical Method

We consider a polygonal domain Ω ⊂ R 2 and a function u : Ω → R that is the solution of the elliptic equation

Throughout this paper, it is assumed that f ∈ L 2 ⁢ ( Ω ) and μ ∈ L ∞ ⁢ ( Ω ) such that there exist μ min , μ max such that

Under these assumptions, there exists a unique weak solution u ∈ H 0 1 ⁢ ( Ω ) of (2.1). Moreover, if μ ∈ W 1 , ∞ ⁢ ( Ω ) and Ω is convex, then u ∈ H 2 ⁢ ( Ω ) ; see [22].

For any h > 0 , let T h be a conformal triangulation of Ω into triangles 𝐾 of diameter h K ≤ h . Let V h be the usual finite element space of continuous, piecewise linear, functions on triangles of T h , with zero value on ∂ ⁡ Ω . To approximate the solution of (2.1), we are looking for u h ∈ V h such that

We will introduce an error estimator ∑ K ∈ T h η K 2 equivalent to the numerical error u - u h , that is to say there exist two constants C ^ 1 , C ^ 2 > 0 that are independent of the data 𝜇, 𝑓, the mesh size and aspect ratio such that

up to higher order terms. It is important to note that (i) (2.3) can be proved provided some local error term is equidistributed in the direction of maximum and minimum stretching – see assumption (4.2) hereafter – which is one of the goals of our adaptive algorithm; (ii) the error estimator η K involves the exact solution, but post-processing techniques can be used to obtain an accurate approximation; (iii) the notion of higher order terms will be addressed more precisely in Remark 2.

3 Anisotropic Estimates for Clément’s Interpolant

In the standard finite element theory, interpolation estimates involve constants that may depend on the aspect ratio and thus yield a posteriori error estimates that are not useful when using anisotropic finite elements. Recently, interpolation estimates that are sharp for meshes with large aspect ratio have been derived. We briefly recall the theoretical framework developed in [20, 21, 29]; similar results can be found in [23]. It should be noted that the maximum angle condition discussed in [6, 5] is not required; see [20, Section 2.1] for details. For simplicity, the presentation is done in a two-dimensional framework; it can easily be extended to the three dimensions.

For any K ∈ T h , we denote by T K : K ^ → K one of the affine transformations mapping the reference triangle K ^ into 𝐾 defined by x = T K ⁢ ( x ^ ) = M K ⁢ x ^ + t K , where M K ∈ R 2 × 2 and t K ∈ R 2 . Since M K is invertible, it admits a singular value decomposition M K = R K T ⁢ Λ K ⁢ P K , where R K and P K are orthogonal matrices and

In the above notation, r 1 , K , r 2 , K are the unit vectors corresponding to the directions of maximum and minimum stretching, respectively, and λ 1 , K , λ 2 , K correspond to the amplitude of the maximum and minimum stretching. A geometrical interpretation is shown in Figure 1. Note that, since T K is not uniquely defined, then neither λ i , K , r i , K , i = 1 , 2 , are.

Figure 1

Transformation T K mapping the usual reference element K ^ into a right triangle 𝐾. The reference triangle is stretched in the direction r 1 , K = ( 1 , 0 ) T (resp. r 2 , K = ( 0 , 1 ) T ), with amplitude λ 1 , K (resp. λ 2 , K ).

Our goal is to derive residual-based, anisotropic a posteriori error estimates, the key ingredient being Clément’s interpolant [14]. When using anisotropic meshes, some additional geometrical assumptions must be made in order to ensure that the constants involved in the interpolation estimates will not depend on the mesh aspect ratio. From now, it is assumed that

1. for each 𝐾, the cardinality of Δ ⁢ K , that is the union of triangles sharing a vertex with 𝐾, is uniformly bounded from above, independently of the mesh geometry,

2. for each 𝐾, the diameter of Δ ⁢ K ^ = T K - 1 ⁢ ( Δ ⁢ K ) is uniformly bounded from above, independently of the mesh geometry.

Let R h : H 1 ⁢ ( Ω ) → V h be the Clément interpolant. Under the two above assumptions, the following interpolation error estimate can be proved.

4 The Anisotropic Error Estimator

We now introduce the error estimator. It is an improvement of the one presented in [34], adapted to the case of a non-constant diffusion coefficient 𝜇. Letting f ∈ L 2 ⁢ ( Ω ) , for all K ∈ T h , we define the L 2 ⁢ ( Ω ) projection of 𝑓 onto the set of constant functions by

For a vector-valued function f = ( f 1 , f 2 ) , we denote Π K f = : ( Π K f 1 , Π K f 2 ) , and for any edge l i , i = 1 , 2 , 3 , of 𝐾, we define

We then define the local error estimator η K 2 by

Here 𝒏 stands for the unit outer normal to 𝐾, and ω K is given by (3.3). For any edge l i , i = 1 , 2 , 3 , of 𝐾, we denote by [ ⋅ ] the jump across the edge l i .

Observe that error estimator (4.1) is not standard since it involves the exact solution in the term ω K ⁢ ( u - u h ) and thus is not fully computable. However, in practice, post-processing techniques can be applied in order to approximate the quantity G K ⁢ ( u - u h ) , contained in ω K ⁢ ( u - u h ) , for instance Zienkiewicz–Zhu (ZZ) post-processing [40, 2, 41]. More precisely, we will replace

where, for any v h ∈ V h , for any vertex 𝑃 of the mesh,

is an approximate L 2 ⁢ ( Ω ) projection of ∂ ⁡ v h ∂ ⁡ x i onto V h .

It is well known [41, 13] that, for elliptic equations and structured meshes, superconvergence of the ZZ recovery occurs, implying that the post-processing is asymptotically exact, that is to say

On general meshes, it was first proven that ∥ Π h ZZ ⁢ ∇ ⁡ u h - ∇ ⁡ u h ∥ L 2 ⁢ ( Ω ) and the true error ∥ ∇ ⁡ u - ∇ ⁡ u h ∥ L 2 ⁢ ( Ω ) are equivalent; see for instance [24]. More recently, the superconvergence of the ZZ gradient recovery was finally shown for unstructured anisotropic meshes [12]. In general, the accuracy of the ZZ post-processing is better than the theoretical predictions, even on anisotropic meshes. We refer for instance to the numerical results presented in [18, 33, 34, 32, 27, 35, 11, 28] in the case of anisotropic meshes for elliptic, parabolic and hyperbolic equations.

Under extra assumptions on 𝑓, 𝜇 and the mesh, it is possible to prove that the error estimator is equivalent to the true numerical error u - u h , up to some higher order terms. The result is stated in the next theorem, whose proof is presented in Appendix A.

5 Numerical Experiments with Non-adapted Meshes

The goal of this section is to numerically verify the equivalence between the true error and the error estimator derived in Section 4. To this end, we consider problems with a known exact solution, and we define the following quantities: the error in H 1 seminorm

the error in the weighted H 1 seminorm

the anisotropic estimator

the anisotropic effectivity index

and the ZZ effectivity index

These quantities should satisfy the following properties:

• ei A independent of the solution 𝑢,

• ei A independent of the variations of 𝜇 and of the ratio μ max μ min ,

• ei A independent of the mesh size and aspect ratio,

• ei ZZ close to one.

be a smoothing of the classical Heavyside function. We consider problem (2.1) in the unit square Ω = ( 0 , 1 ) 2 and choose 𝑓 so that 𝑢 is given by

and 𝜇 is given by

with μ 1 , μ 2 > 0 . Thus 𝜇 is constant except in a thin boundary layer of width 𝜖 where a strong gradient can be observed. The results are reported in Table 1 with unstructured non-adapted meshes of various size, where h 1 and h 2 denote the mesh size in direction x 1 and x 2 , respectively; see Figure 2. We observe that the error estimator is equivalent to the true error uniformly in the mesh size, the mesh aspect ratio and the ratio μ max μ min . Moreover, the values of ei ZZ show that the Zienkiewicz–Zhu error estimator is asymptotically exact.

Figure 2

Example of non-adapted mesh with h 1 = 0.05 and h 2 = 0.5 .

6 An Adaptive Algorithm

We present an adaptive algorithm similar to the one presented in [33, 34]. For the sake of simplicity, we present the two-dimensional version of the algorithm, which can be easily extended to the three-dimensional case.

In all numerical experiments described hereafter, we observe that the anisotropic effectivity index converges to a number close to 3.45 as ℎ goes to zero. In order to keep the effectivity index close to one, we will divide η K by 3.45. To keep the notation simple, we still write η K instead of η K / 3.45 in the sequel.

The goal of the adaptive algorithm is to build a sequence of meshes with possibly large aspect ratio such that the relative estimated error is close to a prescribed tolerance TOL, i.e.

We use the mesh generator BL2D [25] which requires information at the vertices rather than at the triangles. We define the error indicator for each vertex 𝑃 as

and (6.1) can be replaced by

We then equidistribute η P 2 on each vertex 𝑃 by adjusting the mesh with the objective to satisfy

where N v is the number of vertices of T h . In order to satisfy (4.2), we want to equidistribute the error in the directions of minimum and maximum stretchings, and therefore modify the mesh with the goal to satisfy, for i = 1 , 2 ,

where

and

which satisfies 1 2 ≤ σ P ≤ 1 .

The mesh generator BL2D requires for each vertex 𝑃 the local mesh sizes h 1 , P , h 2 , P and the angle θ P between the direction of maximum stretching and the horizontal axis. We align the directions of maximum and minimum stretching with the eigenvectors of the post-processed gradient

More precisely, we set θ P as the angle between the horizontal axis and the eigenvector corresponding to the largest eigenvalue of G P . Concerning the local mesh sizes h 1 , P , h 2 , P , values λ i , P , i = 1 , 2 , are computed at the vertices

Then, if the two inequalities in (6.2) are satisfied, we set h i , P = λ i , P ; if the left inequality in (6.2) is not satisfied, we set h i , P = 1.5 ⁢ λ i , P ; if the right inequality is not satisfied, we set h i , P = λ i , P / 1.5 . Starting from an initial coarse mesh and a coarse value of TOL, a sequence of 40 adapted meshes is generated. Then TOL is divided by two, 40 adapted meshes are generated, and so on. The results are reported hereafter.

Figure 3

Solution on adapted mesh when TOL = 0.025 , when 𝑢 and 𝜇 are given by (6.3) and (5.2), with μ 1 = 1 , μ 2 = 2 and ϵ = 0.01 .

Figure 4

In blue: number of conjugate gradient iterations needed to reach a tolerance 10 - 8 for each iteration of the adaptive algorithm. In red: number of vertices of each adapted mesh at each iteration of the adaptive algorithm. The tolerance TOL is halved every 40 iterations of the adaptive algorithm.

Consider again problem (2.1) in the unit square Ω = ( 0 , 1 ) 2 , with exact solution

and 𝜇 defined by (5.2). We consider μ 1 = 1 , μ 2 = 2 and ϵ = 0.01 . Note that the exact solution is not zero on the boundary anymore; thus, in principle, an extra-term should be added in the error estimator. However, numerical results indicate that this is not needed. Moreover, since the solution is one-dimensional, the adapted meshes will have very large aspect ratio. In Table 2, the results are reported when running the adaptive algorithm starting with TOL = 0.1 and an initial 10 × 10 uniform mesh. We reported also the maximum aspect ratio ar max = max K ∈ T h ⁡ λ 1 , K λ 2 , K , the average aspect ratio ar av , the total CPU time CPU tot and the CPU time spent in the adaptation of the mesh CPU adapt (seconds, Intel Core 2.80 GHz). Numerical results show the sharpness of our anisotropic error indicator, although the mesh aspect ratio is very large. Moreover, most of the CPU time is spent for building adapted meshes rather than solving the linear systems. To illustrate this fact, the number of conjugate gradient iterations needed to solve the linear system together with the number of vertices are reported for all iterations of the adaptive algorithm in Figure 4. Each time the tolerance is halved, first the number of vertices increases and then decreases to a value which is twice the one obtained with the previous tolerance. This is due to the fact that vertices are first added in an isotropic manner and then removed. Moreover, apart from the 40 first mesh iterations, the number of conjugate gradient iterations does not increase with the number of vertices due to the fact that, with mesh iterations, the initial guess of the conjugate gradient method is closer and closer to the solution of the linear system. In Figure 3, the mesh and solution obtained when TOL = 0.025 are shown.

We then present a 3D numerical experiment, the adaptive algorithm and the error estimator described above being adapted to the 3D case; see [17] for details. In order to produce adapted meshes, the 3D Precise Mesh software is used [42]. Letting Ω = ( 0 , 1 ) × ( 0 , 1 ) × ( 0 , 0.1 ) , we choose 𝑓 such that 𝜇 is given by

and 𝑢 by

In Table 3, the results are reported when running the adaptive algorithm starting with an initial tolerance TOL = 0.5 , an initial 10 × 10 × 2 uniform mesh and several values of μ 1 , μ 2 and 𝜖. In Figures 5 and 6, the adapted mesh and solution are shown when TOL = 0.25 and μ 1 = 1 , μ 2 = 2 and ϵ = 0.01 . Then, in order to keep the effectivity index close to one, η K is divided by 3.45, and the adaptive algorithm is run again for several values of TOL when μ 1 = 1 , μ 2 = 100 and ϵ = 0.01 . The results are reported in Table 4. As expected, the effectivity index is now close to one. Again, a non-negligible fraction of the CPU time is needed to adapt the mesh.

Figure 5

Solution on adapted mesh when TOL = 0.25 , 𝑢 and 𝜇 are given by (6.5) and (6.4), with μ 1 = 1 , μ 2 = 2 and ϵ = 0.01 .

Figure 6

Cut at x = 0.5 and y = 0.5 of the adapted mesh obtained when TOL = 0.25 , when 𝑢 and 𝜇 are given by (6.5) and (6.4), with μ 1 = 1 , μ 2 = 2 and ϵ = 0.01 .

7 Conclusion

We introduced an error estimator for the approximation of elliptic problems with strongly varying diffusion coefficient. The error estimator is shown to be equivalent to the true error, up to high order terms. An adaptive algorithm is presented, and its efficiency and accuracy is demonstrated with various numerical experiments.