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

Samuel Dubuis, Paride Passelli and Marco Picasso

# Abstract

The elliptic problem - div ( μ u ) = f is considered, where μ > 0 is smooth but strongly varying. Anisotropic a posteriori error estimates are derived, the effectivity index being bounded above and below by two constants independent of the data 𝑓, 𝜇, the mesh size and aspect ratio, up to higher order terms. Numerical experiments on non-adapted and adapted anisotropic meshes confirm these predictions.

MSC 2010: 65N30; 65N50

## 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

(2.1) { - div ( μ u ) = f in Ω , u = 0 on Ω .

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

0 < μ min μ ( x ) μ max , a.e. x Ω .

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

(2.2) Ω μ u h v h = Ω f v h for all v h V h .

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

(2.3) C ^ 1 Ω μ | ( u - u h ) | 2 K T h η K 2 C ^ 2 Ω μ | ( u - u h ) | 2

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

Λ K = ( λ 1 , K 0 0 λ 2 , K ) , λ 1 , K λ 2 , K > 0 , R K = ( r 1 , K T r 2 , K T ) .

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.

In particular, the second assumption excludes too distorted meshes; see for instance [29]. It guarantees that all the quantities vary smoothly in the neighborhood of every triangle 𝐾. In practice, these assumptions are fulfilled when using available anisotropic mesh generators; see Section 6.

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.

Proposition 1

## Proposition 1 (Anisotropic Clément Interpolation Error Estimate [20, 21, 29, 34])

There exists a constant C ^ > 0 depending only on the reference triangle K ^ , in particular independent of the mesh size and aspect ratio, such that, for any v H 1 ( Ω ) , for any K T h ,

(3.1) v - R h ( v ) L 2 ( K ) C ^ ω K ( v ) ,

and for any edge l i , i = 1 , 2 , 3 , of 𝐾,

(3.2) v - R h ( v ) L 2 ( l i ) C ^ ( | l i | λ 1 , K λ 2 , K ) 1 / 2 ω K ( v ) ,

where

(3.3) ω K 2 ( v ) = λ 1 , K 2 v r 1 , K L 2 ( Δ K ) 2 + λ 2 , K 2 v r 2 , K L 2 ( Δ K ) 2 = λ 1 , K 2 ( r 1 , K T G K ( v ) r 1 , K ) + λ 2 , K 2 ( r 2 , K T G K ( v ) r 2 , K ) ,

with

G K ( v ) = ( Δ K ( v x 1 ) 2 Δ K v x 1 v x 2 Δ K v x 1 v x 2 Δ K ( v x 2 ) 2 ) .

## 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

Π K f = 1 | K | K f .

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

Π l i f = 1 | l i | l i f .

We then define the local error estimator η K 2 by

(4.1) η K 2 = ( Π K f + Π K μ u h L 2 ( K ) + 1 2 i = 1 3 ( | l i | λ 1 , K λ 2 , K ) 1 / 2 [ Π l i μ u h n ] L 2 ( l i ) ) ω K ( u - u h ) .

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

( u - u h ) x i by Π h ZZ u h x i - u h x i , i = 1 , 2 ,

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

Π h ZZ v h x i ( P ) = K T h , P K | K | v h | K x i K T h , P K | K |

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

lim h 0 Π h ZZ u h - u h L 2 ( Ω ) u - u h L 2 ( Ω ) = 1 .

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.

Theorem 1

Assume that f L 2 ( Ω ) and μ W 1 , ( Ω ) ; let u H 0 1 ( Ω ) be the weak solution of (2.1), and let u h V h be the solution of (2.2). Moreover, assume that there exists C ^ > 0 depending only on the reference triangle K ^ such that, for all K T h ,

(4.2) λ 1 , K 2 ( r 1 , K T G K ( u - u h ) r 1 , K ) C ^ λ 2 , K 2 ( r 2 , K T G K ( u - u h ) r 2 , K ) .

Then there exists a constant C ^ 1 > 0 depending only on the reference triangle such that

Ω μ | ( u - u h ) | 2 C ^ 1 K T h ( η K 2 + ε K 2 μ min ) ,

where

ε K 2 = λ 2 , K 2 ( f - Π K f L 2 ( K ) 2 + ( μ - Π K μ ) u h L 2 ( K ) 2 + 1 λ 2 , K i = 1 3 [ ( μ - Π l i μ ) u h n ] L 2 ( l i ) 2 ) .

Moreover, there exists a constant C ^ 2 > 0 depending only on the reference triangle K ^ such that

K T h η K 2 C ^ 2 K T h ( K μ | ( u - u h ) | 2 ( 1 + μ - Π K μ L ( Δ K ) μ min ) + ε K 2 μ min ) .

Remark 1

Assumption (4.2) means that the error in the maximal stretching direction is smaller than or equal to the error in the minimal stretching direction. This assumption was already used in [34] to derive the equivalence between the numerical error and the estimator. One possibility to fulfill (4.2) is to equidistribute for all K T h the error in each direction r 1 , K , r 2 , K ,

(4.3) λ 1 , K 2 ( r 1 , K T G K ( u - u h ) r 1 , K ) = λ 2 , K 2 ( r 2 , K T G K ( u - u h ) r 2 , K ) .

This requirement corresponds to the matching assumption discussed in [23]. Building a mesh such that (4.3) is satisfied will be ensured by the adaptive algorithm presented in Section 6. The algorithm has also the goal to align directions r 1 , K , r 2 , K with respect to the eigenvectors of G K ( u - u h ) .

Remark 2

Assume f H 1 ( Ω ) and μ W 2 , ( Ω ) . Then, from [20, 21], we have

f - Π K f L 2 ( K ) 2 C ^ i = 1 2 λ i , K 2 f r i , K L 2 ( K ) 2 .

Similarly, we have

( μ - Π K μ ) u h L 2 ( K ) 2 C ^ ( i , j = 1 2 λ i , K λ j , K λ 2 , K r i , K T H ( μ ) r j , K L ( K ) ) 2 u h L 2 ( K ) 2 ,

where H ( ) denotes the Hessian matrix,

i = 1 3 [ ( μ - Π l i μ ) u h n ] L 2 ( l i ) 2 C ^ ( i = 1 2 λ i , K 2 μ r i , K L ( K ) 2 ) u h n L 2 ( K ) 2 , μ - Π K μ L ( K ) C ^ ( λ 1 , K μ r 1 , K L ( K ) + λ 2 , K μ r 2 , K L ( K ) ) .

In the isotropic framework, this yields a contribution

K T h ε K 2 = O ( h 3 ) , which is negligible compared to Ω μ | ( u - u h ) | 2 = O ( h 2 ) .

In the anisotropic setting, for instance if 𝑓 and 𝜇 depend only on x 2 and r 1 , K = ( 1 , 0 ) , then ε K 2 = O ( λ 2 , K 3 ) . Thus, in both cases, Theorem 1 indeed yields (2.3) up to higher order terms.

## 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

e H 1 = ( u - u h ) L 2 ( Ω ) ,

the error in the weighted H 1 seminorm

e μ , H 1 = μ 1 / 2 ( u - u h ) L 2 ( Ω ) ,

the anisotropic estimator

η A = ( K T h η K 2 ) 1 / 2 ,

the anisotropic effectivity index

ei A = η A e μ , H 1 ,

and the ZZ effectivity index

ei ZZ = u h - Π h ZZ u h L 2 ( Ω ) e H 1 .

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.

For all x R , ϵ > 0 , let

H ϵ ( x ) = { 0 x - ϵ , x + ϵ 2 ϵ + 1 2 π sin ( π x ϵ ) - ϵ x ϵ , 1 ϵ x ,

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

(5.1) u ( x ) = μ 2 sin ( π x 1 ) sin ( π x 2 ) H ϵ ( x 1 - 0.5 ) + μ 1 sin ( π x 1 ) sin ( π x 2 ) H ϵ ( 0.5 - x 1 ) .

and 𝜇 is given by

(5.2) μ ( x ) = μ 2 H ϵ ( x 1 - 0.5 ) + μ 1 ( 1 - H ϵ ( x 1 - 0.5 ) ) ,

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 .

### Table 1

True error, estimated error and effectivity index for various non-adapted meshes and various choices of μ 1 , μ 2 and 𝜖, when 𝑢 and 𝜇 are given by (5.1) and (5.2).

h 1 - h 2 η A e μ , H 1 e i A η ZZ e H 1 e i ZZ
μ 1 = 1 , μ 2 = 2 , ϵ = 0.1

0.05-0.5 9.33 3.25 2.87 1.98 2.54 0.78
0.025-0.25 6.73 2.15 3.13 1.49 1.64 0.91
0.0125-0.125 3.96 1.38 2.87 0.99 1.05 0.95
0.00625-0.0625 2.10 0.73 2.87 0.53 0.55 0.97
0.003125-0.03125 1.07 0.36 2.93 0.27 0.28 0.98
0.0015625-0.015625 0.54 0.19 2.93 0.14 0.14 0.98
0.00078125-0.0078125 0.27 0.095 2.89 0.07 0.07 0.98

μ 1 = 1 , μ 2 = 2 , ϵ = 0.01

0.05-0.5 91.07 36.24 2.51 20.67 29.83 0.69
0.025-0.25 63.39 23.47 2.70 14.44 19.72 0.73
0.0125-0.125 14.20 6.62 2.15 3.50 5.36 0.65
0.00625-0.0625 6.74 2.20 3.07 2.03 1.77 1.15
0.003125-0.03125 3.27 0.99 3.29 0.87 0.81 1.07
0.0015625-0.015625 1.67 0.51 3.28 0.42 0.42 1.02
0.00078125-0.0078125 0.85 0.26 3.30 0.21 0.21 1.00

μ 1 = 1 , μ 2 = 100 , ϵ = 0.1

0.05-0.5 3584.94 1077.18 3.33 123.56 135.07 0.91
0.025-0.25 2395.20 764.94 3.13 85.30 89.35 0.95
0.0125-0.125 1365.70 469.50 2.91 50.27 52.43 0.95
0.00625-0.0625 715.23 246.22 2.90 26.23 26.89 0.98
0.003125-0.03125 359.39 120.49 2.98 12.94 13.13 0.99
0.0015625-0.015625 182.26 61.13 2.98 6.54 6.63 0.99
0.00078125-0.0078125 92.05 31.22 2.95 3.32 3.37 0.99

μ 1 = 1 , μ 2 = 100 , ϵ = 0.01

0.005-0.05 3258.10 1183.14 2.75 162.95 183.48 0.89
0.0025-0.025 1517.05 450.98 3.36 66.70 65.89 1.01
0.00125-0.0125 777.60 231.73 3.36 32.19 32.01 1.01
0.000625-0.00625 375.00 110.71 3.39 15.60 15.61 1.00
0.0006-0.006 360.27 106.08 3.40 14.87 14.90 1.00

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.

(6.1) 0.75 TOL ( K T h η K 2 Ω μ ( x ) | u h | 2 d x ) 1 / 2 1.25 TOL .

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

η P 2 = K T h P K η K 2 so that P T h η P 2 = 3 K T h η K 2 ,

and (6.1) can be replaced by

0.75 TOL ( P T h η P 2 3 Ω μ ( x ) | u h | 2 d x ) 1 / 2 1.25 TOL .

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

3 N v 0.75 2 TOL 2 Ω μ ( x ) | u h | 2 d x η P 2 3 N v 1.25 2 TOL 2 Ω μ ( x ) | u h | 2 d x ,

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 ,

(6.2) 3 σ P 2 N v 0.75 2 TOL 2 Ω μ ( x ) | u h | 2 d x K T h P K η i , K 2 3 σ P 2 N v 1.25 2 TOL 2 Ω μ ( x ) | u h | 2 d x ,

where

η i , K 2 = ( Π K f + Π K μ u h L 2 ( K ) + 1 2 j = 1 3 ( | l j | λ 1 , K λ 2 , K ) 1 / 2 [ Π l j μ ( x ) u h n ] L 2 ( l j ) ) λ i , K 2 r i , K T G K ( u - u h ) r i , K

and

σ P = K T h , P K ( η 1 , K 2 + η 2 , K 2 ) K T h , P K ( η 1 , K 4 + η 2 , K 4 ) 1 / 2 ,

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

G P = K T h P K G K ( u - u h ) .

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

λ i , P = K T h , P K λ i , K K T h , P K 1 .

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.

### Table 2

True error, effectivity indices and aspect ratio for different values of the tolerance TOL, when 𝑢 and 𝜇 are given by (6.3) and (5.2), with μ 1 = 1 , μ 2 = 2 and ϵ = 0.01 .

TOL Vertices ei A e H 1 e i ZZ ar max ar av CPU tot CPU adapt
0.1 119 0.89 0.43 1.01 1679 355 7 6
0.05 196 0.94 0.22 0.99 3971 769 12 10
0.025 355 0.92 0.11 1.00 6643 1693 17 15
0.0125 741 0.96 0.054 1.00 17847 3055 25 21
0.00625 1388 0.99 0.027 1.00 37973 6683 36 31
0.003125 2822 0.98 0.013 1.00 77257 13292 55 47

### 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

(6.3) u ( x ) = μ 2 sin ( π x 1 ) H ϵ ( x 1 - 0.5 ) + μ 1 sin ( π x 1 ) H ϵ ( 0.5 - x 1 )

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

(6.4) μ ( x ) = μ 2 H ϵ ( x 3 - 0.05 ) + μ 1 ( 1 - H ϵ ( x 3 - 0.05 ) )

and 𝑢 by

(6.5) u ( x ) = μ 2 sin ( π x 1 ) H ϵ ( x 3 - 0.05 ) + μ 1 sin ( π x 1 ) H ϵ ( 0.05 - x 3 ) .

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.

### Table 3

True error, effectivity indices and aspect ratio for different values of tolerance TOL, when 𝑢 and 𝜇 are given by (6.5) and (6.4), with various values for and μ 1 , μ 2 and 𝜖.

TOL Vertices e i A e H 1 e i ZZ ar max ar av
μ 1 = 1 , μ 2 = 2 , ϵ = 0.1

0.5 93 2.82 0.39 0.88 121 28
0.25 397 3.22 0.18 1.00 476 63
0.125 1525 3.28 0.09 1.00 589 99
0.0625 7511 3.27 0.05 0.99 1028 124
0.03125 39512 3.29 0.02 0.99 1769 166

μ 1 = 1 , μ 2 = 2 , ϵ = 0.01

0.5 644 3.61 0.97 1.03 1220 173
0.25 2426 3.65 0.46 0.99 2533 257
0.125 13317 3.67 0.23 0.99 2552 277
0.0625 32479 3.56 0.12 0.98 12101 979
0.03125 214170 3.52 0.07 0.98 33075 1227

μ 1 = 1 , μ 2 = 100 , ϵ = 0.1

0.5 107 3.42 32.62 0.99 234 45
0.25 366 3.45 15.67 1.03 377 72
0.125 1570 3.36 8.22 1.01 645 102
0.0625 7250 3.45 3.95 1.00 2024 132
0.03125 39129 3.44 1.96 0.99 1961 184

### 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 .

### Table 4

True error, effectivity indices and aspect ratio for different values of tolerance TOL, when 𝑢 and 𝜇 are given by (6.5) and (6.4), with μ 1 = 1 , μ 2 = 100 and ϵ = 0.01 .

TOL Vertices ei A e H 1 e i ZZ ar max ar av CPU tot CPU adapt
0.25 191 1.06 195.85 1.01 864 143 23 19
0.125 452 1.09 85.17 1.04 4542 310 46 39
0.0625 1283 1.03 47.37 1.01 6684 694 94 78
0.03125 6426 1.01 26.36 0.99 15832 1098 292 237
0.015625 33745 1.10 13.60 0.99 117521 1588 2298 1402

## 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.

Funding statement: Paride Passelli is financed by Rio Tinto Aluminium LRF Research Center at Saint Jean de Maurienne (EPFL industrial grant).

## A Appendix: Proof of Theorem 1

We first prove the upper bound. Letting e = u - u h , using (2.1) and (2.2), we have, for any v h V h ,

Ω μ | e | 2 = Ω f ( e - v h ) - Ω μ u h ( e - v h ) .

By integration by parts over the triangles 𝐾, we obtain

Ω μ | e | 2 = K T h K ( f + div ( μ u h ) ) ( e - v h ) + 1 2 K [ μ u h n ] ( e - v h )
= K T h K ( Π K f + Π K μ u h ) ( e - v h ) + 1 2 i = 1 3 l i [ Π l i μ u h n ] ( e - v h )
+ K T h K ( f - Π K f ) ( e - v h ) + K T h K ( μ - Π K μ ) u h ( e - v h )
+ 1 2 K T h i = 1 3 l i [ ( μ - Π l i μ ) u h n ] ( e - v h ) .
Using the Cauchy–Schwarz inequality, choosing v h = R h e and by using interpolation error estimates (3.1) and (3.2), we obtain

Ω μ | e | 2 C ^ ( K T h η K 2 + K T h ( f - Π K f L 2 ( K ) + ( μ - Π K μ ) u h L 2 ( K ) + 1 2 λ 2 , K i = 1 3 [ ( μ - Π l i μ ) u h n ] L 2 ( l i ) ) ω K ( e ) ) .

Using assumption (4.2) and the fact that r 1 , K , r 2 , K form a basis, we have

ω K 2 ( e ) C ^ λ 2 , K 2 ( r 1 , K T G K ( e ) r 1 , K + r 2 , K T G K ( e ) r 2 , K ) = C ^ λ 2 , K 2 e L 2 ( Δ K ) 2 .

We obtain the result by using the discrete Cauchy–Schwarz inequality.

In order to prove the lower bound, we use the standard bubble functions [8, 37], adapted to the anisotropic case in [34], and modified here to account for the variations of 𝜇.

Proposition 2

There exist a function φ H 0 1 ( Ω ) and a constant C ^ > 0 (that depends only on the reference triangle K ^ ) such that, for any K T h ,

(A.1) l i [ Π l i μ u h n ] φ = ( | l i | λ 1 , K λ 2 , K ) 1 / 2 [ Π l i μ u h n ] L 2 ( l i ) ω K ( e ) , i = 1 , 2 , 3 ,
(A.2) K ( Π K f + Π K μ u h ) φ = Π K f + Π K μ u h L 2 ( K ) ω K ( e ) ,
(A.3) ( λ 1 , K 2 φ r 1 , K L 2 ( K ) 2 + λ 2 , K 2 φ r 2 , K L 2 ( K ) 2 ) 1 / 2 C ^ ω K ( e ) ,
(A.4) K μ | φ | 2 C ^ ( Π K μ + μ - Π K μ L ( K ) ) ω K 2 ( e ) λ 2 , K 2 .

## Proof

We claim that

φ | K = C ^ K Ψ K + i = 1 3 C ^ l i Ψ l i ,

where C ^ K , C ^ l i are constants that depend only on K ^ and that have to be computed, and Ψ K , Ψ l i are the usual bubble functions over 𝐾 and its edges l i , i = 1 , 2 , 3 . For (A.1) and (A.2), the proof follows the same steps as the one of [34, Proposition 2]. In order to prove (A.3), we proceed again as in [34, Proposition 2] obtaining, for i = 1 , 2 ,

K | φ r i , K | 2 C ^ ω K 2 ( e ) λ 1 , K λ 2 , K ( K | Ψ K r i , K | 2 + j = 1 3 K | Ψ l j r i , K | 2 ) ,

which together with [21, relation (2.14)] gives (A.3). To prove (A.4), we note that

K μ | φ | 2 = K Π K μ | φ | 2 + K ( μ - Π K μ ) | φ | 2

and then apply [34, Proposition 4]. ∎

We can now prove the lower bound. Using the definition of η K 2 and identities (A.2) and (A.1), one can write

K T h η K 2 = K T h K ( Π K f + Π K μ u h ) φ + 1 2 i = 1 3 l i [ Π l i μ u h n ] φ ,

where φ H 0 1 ( Ω ) is the function given by Proposition 2. Therefore, adding and subtracting the correct quantities in the right-hand side, we have

K T h η K 2 = K T h K μ e φ + K T h K ( Π K f - f ) φ + K T h K ( Π K μ - μ ) u h φ + 1 2 K T h i = 1 3 l i [ ( Π l i μ - μ ) u h n ] φ .

Using the properties of the L 2 projection Π K and Π l i , we obtain

K T h η K 2 = K T h K μ e φ + K T h K ( Π K f - f ) ( φ - Π K φ ) + K T h K ( Π K μ - μ ) u h ( φ - Π K φ ) + 1 2 K T h i = 1 3 l i [ ( Π l i μ - μ ) u h n ] ( φ - Π l i φ ) .

Integration by parts over the triangles K T h and using the Cauchy–Schwarz inequality, we get

K T h η K 2 K T h μ 1 / 2 e L 2 ( K ) μ 1 / 2 φ L 2 ( K ) + K T h ( ( Π K f - f ) L 2 ( K ) + ( Π K μ - μ ) u h L 2 ( K ) ) φ - Π K φ L 2 ( K ) + 1 2 K T h i = 1 3 [ ( Π l i μ - μ ) u h n ] L 2 ( l i ) φ - Π l i φ L 2 ( l i ) .

Now, using results from [20, 21], we have that there exists a constant C ^ depending only on the reference triangle such that, for all 𝐾 and all i = 1 , 2 , 3 ,

φ - Π K φ L 2 ( K ) C ^ ( λ 1 , K 2 φ r 1 , K L 2 ( K ) 2 + λ 2 , K 2 φ r 2 , K L 2 ( K ) 2 ) 1 / 2 , φ - Π K φ L 2 ( l i ) C ^ λ 2 , K ( λ 1 , K 2 φ r 1 , K L 2 ( K ) 2 + λ 2 , K 2 φ r 2 , K L 2 ( K ) 2 ) 1 / 2 .

Finally, using the two bounds (A.3) and (A.4), Young’s and the Cauchy–Schwarz inequalities, we obtain

(A.5) K T h η K 2 C ^ ( Ω μ | e | 2 + K T h ( Π K μ + μ - Π K μ L ( K ) ) ω K 2 ( e ) λ 2 , K 2 + K T h ( Π K f - f L 2 ( K ) + ( μ - Π K μ ) u h L 2 ( K ) + 1 2 λ 2 , K i = 1 3 [ ( μ - Π l i μ ) u h n ] L 2 ( l i ) ) ω K ( e ) ) .

In order to conclude the proof, we note that

( Π K μ + μ - Π K μ L ( K ) ) ω K 2 ( e ) λ 2 , K 2 C ^ μ 1 / 2 e L 2 ( Δ K ) 2 + C ^ μ - Π K μ L ( Δ K ) μ min μ 1 / 2 e L 2 ( Δ K ) 2

The final result is obtained inserting the above estimate in (A.5).

# Acknowledgements

The referees are acknowledged for fruitful suggestions.

### References

[1] M. Ainsworth and J. T. Oden, A posteriori error estimation in finite element analysis, Comput. Methods Appl. Mech. Engrg. 142 (1997), no. 1–2, 1–88. 10.1002/9781118032824Search in Google Scholar

[2] M. Ainsworth, J. Z. Zhu, A. W. Craig and O. C. Zienkiewicz, Analysis of the Zienkiewicz–Zhu a posteriori error estimator in the finite element method, Internat. J. Numer. Methods Engrg. 28 (1989), no. 9, 2161–2174. 10.1002/nme.1620280912Search in Google Scholar

[3] F. Alauzet, P. Frey and P. George, Anisotropic mesh adaptation for Rayleigh–Taylor instabilities, ECCOMAS 2004 - European Congress on Computational Methods in Applied Sciences and Engineering, University of Jyväskylä, Jyväskylä (2004), https://www.researchgate.net/publication/283138589. Search in Google Scholar

[4] F. Alauzet and A. Loseille, A decade of progress on anisotropic mesh adaptation for computational fluid dynamics, Comput.-Aided Des. 72 (2016), 13–39. 10.1016/j.cad.2015.09.005Search in Google Scholar

[5] T. Apel, Anisotropic Finite Elements: Local Estimates and Applications, Adv. Numer. Math., B. G. Teubner, Stuttgart, 1999. Search in Google Scholar

[6] I. Babuška and A. K. Aziz, On the angle condition in the finite element method, SIAM J. Numer. Anal. 13 (1976), no. 2, 214–226. 10.1137/0713021Search in Google Scholar

[7] I. Babuška, J. Chandra and J. E. Flaherty, Adaptive Computational Methods for Partial Differential Equations, Society for Industrial and Applied Mathematics, Philadelphia, 1983. Search in Google Scholar

[8] I. Babuška, R. Durán and R. Rodríguez, Analysis of the efficiency of an a posteriori error estimator for linear triangular finite elements, SIAM J. Numer. Anal. 29 (1992), no. 4, 947–964. 10.1137/0729058Search in Google Scholar

[9] W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations, Lectures in Math. ETH Zürich, Birkhäuser, Basel, 2003. 10.1007/978-3-0348-7605-6Search in Google Scholar

[10] E. Boey, Y. Bourgault and T. Giordano, Anisotropic space-time adaptation for reaction-diffusion problems, preprint (2017), https://arxiv.org/abs/1707.04787. Search in Google Scholar

[11] Y. Bourgault and M. Picasso, Anisotropic error estimates and space adaptivity for a semidiscrete finite element approximation of the transient transport equation, SIAM J. Sci. Comput. 35 (2013), no. 2, A1192–-A1211. 10.1137/120891320Search in Google Scholar

[12] W. Cao, Superconvergence analysis of the linear finite element method and a gradient recovery postprocessing on anisotropic meshes, Math. Comp. 84 (2015), no. 291, 89–117. 10.1090/S0025-5718-2014-02846-9Search in Google Scholar

[13] C. Carstensen, All first-order averaging techniques for a posteriori finite element error control on unstructured grids are efficient and reliable, Math. Comp. 73 (2004), no. 247, 1153–1165. 10.1090/S0025-5718-03-01580-1Search in Google Scholar

[14] P. Clément, Approximation by finite element functions using local regularization, Rev. Française Automat. Informat. Recherche Opérationnelle Sér. 9 (1975), no. R2, 77–84. 10.1051/m2an/197509R200771Search in Google Scholar

[15] T. Coupez, Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing, J. Comput. Phys. 230 (2011), no. 7, 2391–2405. 10.1016/j.jcp.2010.11.041Search in Google Scholar

[16] W. Dörfler, A convergent adaptive algorithm for Poisson’s equation, SIAM J. Numer. Anal. 33 (1996), no. 3, 1106–1124. 10.1137/0733054Search in Google Scholar

[17] S. Dubuis, Adaptive algorithms for two fluids flows with anisotropic finite elements and order two time discretizations, Ph.D. thesis, Ecole polytechnique fédérale de Lausanne, 2019. Search in Google Scholar

[18] S. Dubuis and M. Picasso, An adaptive algorithm for the time dependent transport equation with anisotropic finite elements and the Crank–Nicolson scheme, J. Sci. Comput. 75 (2018), no. 1, 350–375. 10.1007/s10915-017-0537-1Search in Google Scholar

[19] L. Formaggia, S. Micheletti and S. Perotto, Anisotropic mesh adaption in computational fluid dynamics: Application to the advection-diffusion-reaction and the Stokes problems, Appl. Numer. Math. 51 (2004), no. 4, 511–533. 10.1016/j.apnum.2004.06.007Search in Google Scholar

[20] L. Formaggia and S. Perotto, New anisotropic a priori error estimates, Numer. Math. 89 (2001), no. 4, 641–667. 10.1007/s002110100273Search in Google Scholar

[21] L. Formaggia and S. Perotto, Anisotropic error estimates for elliptic problems, Numer. Math. 94 (2003), no. 1, 67–92. 10.1007/s00211-002-0415-zSearch in Google Scholar

[22] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Class. Appl. Math. 69, Society for Industrial and Applied Mathematics, Philadelphia, 2011. 10.1137/1.9781611972030Search in Google Scholar

[23] G. Kunert, An a posteriori residual error estimator for the finite element method on anisotropic tetrahedral meshes, Numer. Math. 86 (2000), no. 3, 471–490. 10.1007/s002110000170Search in Google Scholar

[24] G. Kunert and S. Nicaise, Zienkiewicz–Zhu error estimators on anisotropic tetrahedral and triangular finite element meshes, M2AN Math. Model. Numer. Anal. 37 (2003), no. 6, 1013–1043. 10.1051/m2an:2003065Search in Google Scholar

[25] P. Laug and H. Borouchaki, The BL2D mesh generator - beginner’s guide, user’s and programmer’s manual, Technical report RT-0194, Institut National de Recherche en Informatique et Automatique (INRIA), Rocquencourt-Le Chesnay, 1996. Search in Google Scholar

[26] A. Loseille and F. Alauzet, Continuous mesh framework part I: Well-posed continuous interpolation error, SIAM J. Numer. Anal. 49 (2011), no. 1, 38–60. 10.1137/090754078Search in Google Scholar

[27] A. Lozinski, M. Picasso and V. Prachittham, An anisotropic error estimator for the Crank–Nicolson method: application to a parabolic problem, SIAM J. Sci. Comput. 31 (2009), no. 4, 2757–2783. 10.1137/080715135Search in Google Scholar

[28] S. Micheletti and S. Perotto, Reliability and efficiency of an anisotropic Zienkiewicz–Zhu error estimator, Comput. Methods Appl. Mech. Engrg. 195 (2006), no. 9–12, 799–835. 10.1016/j.cma.2005.02.009Search in Google Scholar

[29] S. Micheletti, S. Perotto and M. Picasso, Stabilized finite elements on anisotropic meshes: A priori error estimates for the advection-diffusion and the Stokes problems, SIAM J. Numer. Anal. 41 (2003), no. 3, 1131–1162. 10.1137/S0036142902403759Search in Google Scholar

[30] J.-M. Mirebeau, Optimally adapted meshes for finite elements of arbitrary order and W 1 , p norms, Numer. Math. 120 (2012), no. 2, 271–305. 10.1007/s00211-011-0412-1Search in Google Scholar

[31] M. Picasso, Adaptive finite elements for a linear parabolic problem, Comput. Methods Appl. Mech. Engrg. 167 (1998), no. 3–4, 223–237. 10.1016/S0045-7825(98)00121-2Search in Google Scholar

[32] M. Picasso, An anisotropic error indicator based on Zienkiewicz–Zhu error estimator: Application to elliptic and parabolic problems, SIAM J. Sci. Comput. 24 (2003), no. 4, 1328–1355. 10.1137/S1064827501398578Search in Google Scholar

[33] M. Picasso, Numerical study of the effectivity index for an anisotropic error indicator based on Zienkiewicz–Zhu error estimator, Comm. Numer. Methods Engrg. 19 (2003), no. 1, 13–23. 10.1002/cnm.546Search in Google Scholar

[34] M. Picasso, Adaptive finite elements with large aspect ratio based on an anisotropic error estimator involving first order derivatives, Comput. Methods Appl. Mech. Engrg. 196 (2006), no. 1–3, 14–23. 10.1016/j.cma.2005.11.018Search in Google Scholar

[35] M. Picasso, Numerical study of an anisotropic error estimator in the L 2 ( H 1 ) norm for the finite element discretization of the wave equation, SIAM J. Sci. Comput. 32 (2010), no. 4, 2213–2234. 10.1137/090778249Search in Google Scholar

[36] R. Verfürth, A posteriori error estimation and adaptive mesh-refinement techniques, J. Comput. Appl. Math. 50 (1994), 67–83. 10.1016/0377-0427(94)90290-9Search in Google Scholar

[37] R. Verfürth, A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques, Wiley-Teubner, New York, 1996. Search in Google Scholar

[38] R. Verfürth, A Posteriori Error Estimation Techniques for Finite Element Methods, Oxford University, Oxford, 2013. 10.1093/acprof:oso/9780199679423.001.0001Search in Google Scholar

[39] J. Xu and Z. Zhang, Analysis of recovery type a posteriori error estimators for mildly structured grids, Math. Comp. 73 (2004), no. 247, 1139–1152. 10.1090/S0025-5718-03-01600-4Search in Google Scholar

[40] O. C. Zienkiewicz and J. Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, Internat. J. Numer. Methods Engrg. 24 (1987), no. 2, 337–357. 10.1002/nme.1620240206Search in Google Scholar

[41] O. C. Zienkiewicz and J. Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. I. The recovery technique, Internat. J. Numer. Methods Engrg. 33 (1992), no. 7, 1331–1364. 10.1002/nme.1620330702Search in Google Scholar

[42] Spatial Corp Headquarters, Broomfield, 3D Precise Mesh, www.3ds.com. Search in Google Scholar