The elliptic problem is considered, where 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.
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 , where 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 and a function that is the solution of the elliptic equation
Throughout this paper, it is assumed that and such that there exist such that
For any , let be a conformal triangulation of Ω into triangles 𝐾 of diameter . Let be the usual finite element space of continuous, piecewise linear, functions on triangles of , with zero value on . To approximate the solution of (2.1), we are looking for such that
We will introduce an error estimator equivalent to the numerical error , that is to say there exist two constants 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 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 . 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 , we denote by one of the affine transformations mapping the reference triangle into 𝐾 defined by , where and . Since is invertible, it admits a singular value decomposition , where and are orthogonal matrices and
In the above notation, are the unit vectors corresponding to the directions of maximum and minimum stretching, respectively, and correspond to the amplitude of the maximum and minimum stretching. A geometrical interpretation is shown in Figure 1. Note that, since is not uniquely defined, then neither , , are.
Our goal is to derive residual-based, anisotropic a posteriori error estimates, the key ingredient being Clément’s interpolant . 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
for each 𝐾, the cardinality of , that is the union of triangles sharing a vertex with 𝐾, is uniformly bounded from above, independently of the mesh geometry,
for each 𝐾, the diameter of is uniformly bounded from above, independently of the mesh geometry.
Let be the Clément interpolant. Under the two above assumptions, the following interpolation error estimate can be proved.
There exists a constant depending only on the reference triangle , in particular independent of the mesh size and aspect ratio, such that, for any , for any ,
and for any edge , , of 𝐾,
4 The Anisotropic Error Estimator
We now introduce the error estimator. It is an improvement of the one presented in , adapted to the case of a non-constant diffusion coefficient 𝜇. Letting , for all , we define the projection of 𝑓 onto the set of constant functions by
For a vector-valued function , we denote , and for any edge , , of 𝐾, we define
We then define the local error estimator by
Here 𝒏 stands for the unit outer normal to 𝐾, and is given by (3.3). For any edge , , of 𝐾, we denote by the jump across the edge .
Observe that error estimator (4.1) is not standard since it involves the exact solution in the term and thus is not fully computable. However, in practice, post-processing techniques can be applied in order to approximate the quantity , contained in , for instance Zienkiewicz–Zhu (ZZ) post-processing [40, 2, 41]. More precisely, we will replace
where, for any , for any vertex 𝑃 of the mesh,
is an approximate projection of onto .
On general meshes, it was first proven that and the true error are equivalent; see for instance . More recently, the superconvergence of the ZZ gradient recovery was finally shown for unstructured anisotropic meshes . 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 , up to some higher order terms. The result is stated in the next theorem, whose proof is presented in Appendix A.
Then there exists a constant depending only on the reference triangle such that
Moreover, there exists a constant depending only on the reference triangle such that
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  to derive the equivalence between the numerical error and the estimator. One possibility to fulfill (4.2) is to equidistribute for all the error in each direction ,
This requirement corresponds to the matching assumption discussed in . 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 with respect to the eigenvectors of .
Similarly, we have
where denotes the Hessian matrix,
In the isotropic framework, this yields a contribution
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 seminorm
the error in the weighted seminorm
the anisotropic estimator
the anisotropic effectivity index
and the ZZ effectivity index
These quantities should satisfy the following properties:
independent of the solution 𝑢,
independent of the variations of 𝜇 and of the ratio ,
independent of the mesh size and aspect ratio,
close to one.
be a smoothing of the classical Heavyside function. We consider problem (2.1) in the unit square and choose 𝑓 so that 𝑢 is given by
and 𝜇 is given by
with . 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 and denote the mesh size in direction and , 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 . Moreover, the values of show that the Zienkiewicz–Zhu error estimator is asymptotically exact.
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 by 3.45. To keep the notation simple, we still write instead of 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  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 on each vertex 𝑃 by adjusting the mesh with the objective to satisfy
where is the number of vertices of . 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 ,
which satisfies .
The mesh generator BL2D requires for each vertex 𝑃 the local mesh sizes , and the angle 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 as the angle between the horizontal axis and the eigenvector corresponding to the largest eigenvalue of . Concerning the local mesh sizes , , values , , are computed at the vertices
Then, if the two inequalities in (6.2) are satisfied, we set ; if the left inequality in (6.2) is not satisfied, we set ; if the right inequality is not satisfied, we set . 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.
Consider again problem (2.1) in the unit square , with exact solution
and 𝜇 defined by (5.2). We consider , and . 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 and an initial uniform mesh. We reported also the maximum aspect ratio , the average aspect ratio , the total CPU time and the CPU time spent in the adaptation of the mesh (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 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  for details. In order to produce adapted meshes, the 3D Precise Mesh software is used . Letting , 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 , an initial uniform mesh and several values of , and 𝜖. In Figures 5 and 6, the adapted mesh and solution are shown when and , and . Then, in order to keep the effectivity index close to one, is divided by 3.45, and the adaptive algorithm is run again for several values of TOL when , and . 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.
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
Using assumption (4.2) and the fact that , form a basis, we have
We obtain the result by using the discrete Cauchy–Schwarz inequality.
There exist a function
and a constant
(that depends only on the reference triangle
) such that, for any
We claim that
where are constants that depend only on and that have to be computed, and are the usual bubble functions over 𝐾 and its edges , . 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 ,
and then apply [34, Proposition 4]. ∎
where is the function given by Proposition 2. Therefore, adding and subtracting the correct quantities in the right-hand side, we have
Using the properties of the projection and , we obtain
Integration by parts over the triangles and using the Cauchy–Schwarz inequality, we get
In order to conclude the proof, we note that
The final result is obtained inserting the above estimate in (A.5).
The referees are acknowledged for fruitful suggestions.
 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
 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
 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
 T. Apel, Anisotropic Finite Elements: Local Estimates and Applications, Adv. Numer. Math., B. G. Teubner, Stuttgart, 1999. Search in Google Scholar
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 M. Picasso, Numerical study of an anisotropic error estimator in the 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
 R. Verfürth, A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques, Wiley-Teubner, New York, 1996. Search in Google Scholar
 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
 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
 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
© 2022 the author(s), published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.