An approximation of anisotropic metrics from higher order interpolation error for triangular mesh adaptation

Journal of Computational and Applied Mathematics 258 (2014) 99–115

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

An approximation of anisotropic metrics from higher order interpolation error for triangular mesh adaptation Frédéric Hecht a,b , Raphaël Kuate a,c,∗ a

UPMC Univ Paris 06, UMR 7598, Laboratoire Jacques-Louis Lions, F-75005, Paris, France

b

CNRS, UMR 7598, Laboratoire Jacques-Louis Lions, F-75005, Paris, France

c

Université Paris-Est, IFSTTAR, COSYS, F-77447, Marne-la-Vallée, France

article

info

Article history: Received 8 March 2012 Received in revised form 7 June 2013 Keywords: Mesh adaptation Metrics Interpolation error estimates Finite elements Anisotropy

abstract We propose an efficient algorithm for the numerical approximation of metrics, used for anisotropic mesh adaptation on triangular meshes with finite element computations. We derive the metrics from interpolation error estimates expressed in terms of higher order derivatives, for the Pk -Lagrange finite element, k > 1. Numerical examples of mesh adaptation done using metrics computed with our Algorithm, and derived from higher order derivatives as error estimates, show that we obtain the right directions of anisotropy. © 2013 Elsevier B.V. All rights reserved.

1. Introduction An important field of computational fluid dynamics (CFD) is the accuracy of the result obtained from a numerical scheme. One of the means used to optimize the accuracy/cost of a numerical scheme in CFD is mesh adaptation. For a given numerical scheme, the aim of a mesh adaptation process is to change locally the size of the mesh with respect to the behavior of the physics solved. Two main tools are generally used to perform mesh adaptation: the metric [1] which defines the mesh size and the error indicator which gives information on the accuracy of the solution. We do not consider meshes with curved elements (isoparametric). We call anisotropic mesh [2–8] a mesh for which in some regions, elements (triangles) are stretched in a direction. This notion of anisotropy has not exactly the usual meaning. A mesh is usually called isotropic or shape regular when the aspect ratio of its elements is uniformly bounded, see [9,10]. Some algorithms [11–15] used to generate anisotropic adapted meshes suitable for a PDE numerical solution take the Hessian matrix of the computed solution as error estimates. However, the derivation of metrics from the Hessian matrix is justified only [11,12] for the piecewise linear Lagrange finite element. So there is the problem of generating a metric when the Lagrange interpolation needed is a polynomial of degree k, k > 1. The anisotropic behavior of derivative tensors of the interpolation of a function u of dimension d can be characterized by the largest ellipse contained in the level curve of the polynomial for directional derivatives [16,17]. Algebraically, for any positive integer k, the (k + 1)th order tensor for the partial derivatives of a function u at a point has its anisotropic behavior characterized by a suitable d × d positive definite matrix Q satisfying some constraints K , which involves a nonlinear minimization with respect to the matrix Q. In [16], W. Cao proposed a method to find a suboptimal solution to the minimization problem, and gives an approximate Q for the anisotropic measure. Let λ1 , λ2 , . . . , λd be the eigenvalues (in ascending order) of Q, and v1 , v2 , . . . , vd be the

∗

Correspondence to: Universite Pierre et Marie Curie, Laboratoire Jacques-Louis Lions, Boite courrier 187, Cedex 05, 75252 Paris, France. E-mail addresses: [email protected] (F. Hecht), [email protected], [email protected], [email protected] (R. Kuate).

0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.09.002

100

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

corresponding eigenvectors. Cao approximates Q by choosing successively (λi , vi ) for i = d, d − 1, . . . , 1, such that each λi k+1

is the smallest possible to have the constraints K hold. More precisely, he first chooses λd such that |λd | 2 is the largest (k + 1)th order directional derivative of u. The corresponding eigenvector vd is chosen as the unit vector along which the (k + 1)th directional derivative is the largest. To determine the remaining d − 1 eigenpairs, he uses a set of orthonormal bases for the orthogonal complement of span{vd } in Rd and expresses recursively the initial problem in terms of a combination of the λi , i = d, d − 1, . . . already evaluated and a problem of one dimension less than the latter. In [18], J.-M. Mirebeau proposed an exact solution of the minimization problem for some cases of second order and third order derivatives. We propose an algorithm which solves the minimization problem by a numerical optimization of the parameters of the ellipse which characterizes the anisotropic behavior of the derivatives of the interpolation of u. A naive computation of the ellipse yields a cubic complexity, since an ellipse has three parameters once its center is known. We use the error estimates described in [17], which are similar to those in [16]. We discretize the level curve of the polynomial directional derivatives (error curve) into n points, with no regard to the number of vertices of the mesh. We assume that the exact largest ellipse contained in the error curve does not pass far away from the point M0 of the error curve which is the closest to the center of the ellipse. We then choose a ‘‘small’’ distance ϵ0 between the approximated ellipse and the error curve at M0 , which reduces the complexity of the problem to a quadratic complexity. We thus find an ellipse of the largest area which passes at a distance ϵ0 of M0 and at another ‘‘small’’ distance ϵ of any other point of the error curve, such that the constraints of being included into the error curve be satisfied. Then, for any other point Mi of the error curve considered (except M0 ) such that the ellipse found passes at a distance ϵ of Mi , we obtain a minimization problem in terms of one variate, of which the exact solution is known ‘‘by hand’’. The algorithm then tests the constraints for all the remaining points of the error curve, except M0 and Mi . The points Mi of which all the constraints are tested are those which lie in regions of the error curve such that the approximated ellipse be tangent to the error curve. The cost of our Algorithm is then less than a quadratic cost, with respect to n. The paper is organized as follows. In Section 2, we give a short description of the error estimates [17] used. In Section 3, we recall some notions about metrics [12] and give a formulation of the optimization problem that goes with the characterization of the anisotropic behavior of derivative tensors of the interpolation of a function u by the largest ellipse contained in the error curve. In Section 4, we explain the algebraic consequences of the equations obtained from the formulation of the optimization problem. We then present the algorithm of the resolution of the optimization problem in Section 5, and some illustrating examples on the convergence of metrics computed. In Section 6, we present numerical experiments of metrics computed, and examples of mesh adaptation for a function with metrics generated by our Algorithm with third order, fourth order and fifth order derivatives as local interpolation error estimates. We also compare the results obtained with mesh adaptation using metrics based on the Hessian matrix. Our Algorithm is coded, included into the software Freefem++ [19] as a library (MetricKuate) and the computations are done with Freefem++. 2. The error estimates We give a short description of the error estimates proposed in [17]. 2.1. Definitions We use the notations [20, Chapter VII], which are introduced in the books [21] [22, Chapter II].

• D(ℓ) u is the derivative of u of order ℓ, so if u is a real function, we have D(0) u = u. • K is a triangle, hK is the diameter of its circumscribed ball, ρK is the diameter of its inscribed ball and |K | is the area of K . • W ℓ,p (K ) is the set of functions of Lp (K ) such that all the partial derivatives of order ≤ ℓ are in Lp (K ).  1p   • |u|W ℓ,p (K ) = K α∈{1,2}ℓ |∂α u|p . Definition 2.1 (2-Simplex Finite Element of Degree k). A finite element (K , PK , ΣK ) will be called 2-simplex of degree k, if Pk (K ) ⊂ PK where Pk (K ) is the set of polynomials P [R2 ] restricted to K of degree less than or equal to k. ΣK := {li (p) | p ∈ PK , 1 ≤ i ≤ nK }, li ∈ P ′ K , 1 ≤ i ≤ nK = dim PK , are bounded linear functionals li : PK → R. Let IKk be the interpolation operator with values in PK and which keep unchanged the polynomials of Pk (K ), i.e. if u ∈ Pk (K ), then u = IKk u. Definition 2.2 (Reference Element). The reference element is a regular 2-simplex denoted Kˆ , chosen with the size 1. Definition 2.3. A constant metric M of R2 is a symmetric matrix positively definite that allows us to define a distance in the Euclidean metric space R2 . Let FK be the affine transformation such that K = FK (Kˆ ), we denote: BK the linear part of FK which corresponds to the derivate of FK , and det(BK ) the determinant of BK which is the Jacobian of FK . 1 −1 The size hK ,ξ of an element K in the direction ξ is the length of the edge of K in the direction ξ of that edge. ∥B− is K ∥ the smallest size of the element denoted h− = inf h . This element size can be represented by an ellipse associated to ξ ∈S2 K ,ξ K the ‘‘natural’’ metric of the element introduced in [23,13].

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

101

Let u be a real function of K , we define its corresponding function uˆ = u ◦ FK on Kˆ . We have x = FK (ˆx),

uˆ (ˆx) = u(x),

u = uˆ ◦ FK−1 .

A finite element K is affine-equivalent if it has the following property k  k ˆ. I K u = IKˆ u

(2.1)

Theorem 2.1 ([17]). For any finite element (K , PK , ΣK ) of Definition 2.1 satisfying Property (2.1), with an associated metric MK , let u be a function. For all integer ℓ ∈ [0..k + 1] and for all (p, q) ∈ [1, ∞]2 such that W ℓ,p (K ) be included in C 0 (K ), the set of continuous real-valued functions on K and u ∈ W ℓ,p (K ). If it exists a real constant κ ∈ R such that ℓ

− ℓ ℓ 2 ∥(D(ℓ) u)(ξ , . . . , ξ )∥Lp (K ) ≤ κ ∥ξ ∥ℓMK (h− K ) = κ (ξ , MK ξ ) (hK ) ,

(2.2)

then we have 1

1

ℓ ∥u − IKk u∥Lq (K ) ≤ κ C |K | q − p (h− K) ,

where C is a real positive constant depending only on ℓ, p, q. And if W ℓ,p (K ) ⊂ W 1,q (K ) then we have 1

1

ℓ−1 |u − IKk u|W 1,q (K ) ≤ κ C |K | q − p (h− , K)

with C positive constant depending only on ℓ, p, q. Proof. See [17] for the proof. 3. Description of the problem We first recall some notions on metrics. − → − → The scalar product of two vectors U and V under the metric space (R2 , M ) is defined as: [24]

− → − → − → − → ⟨ U , V ⟩M = t U M V ∈ R. − →

The Euclidean norm of a vector AB or the distance between the points A and B for the metric M is defined by:

− → dM (A, B) = ∥ AB ∥M =



− → − →

t AB M AB

.

There is the situation for which the metric is non constant but continuous. In this case we have a metric field M (Y ) defined

− →

in each point Y of the space. The norm of a vector AB is then given by

− →

dM (A, B) = ∥ AB ∥M =

 1

− →

t AB M

− →− → (A + t AB ) ABdt .

0

A metric M can be represented by its unit ball [12]. The unit ball under the metric M defined at the vertex P is the set of points Q which satisfy:



− → ∥PQ ∥M =

− → − →

t PQ M PQ

= 1.

Hessian based metric. We recall the derivation of a metric from the Hessian matrix by its eigendecomposition [25]. Let N be

− → − →

− →

a real symmetric and invertible matrix, (λ1 , λ2 ) its eigenvalues and (µ1 , µ2 ) its eigenvectors. Let V be a vector, such that  

− →

V =

u1 u2

− →,µ − → (µ 1 2)

. Then,

− → − → 2 2 → → |t V N V | = |u21 λ1 − µ1 + u22 λ2 − µ2 | 2 2 → → ≤ u21 − µ1 |λ1 | + u22 − µ2 |λ2 |    → − |λ1 | → − → t− = V µ1 µ2 

0

|λ2 |

0

 − → − →

− →

→ i.e., |t V N V | ≤ t V − µ1



→ − µ1



 |λ1 | − → µ2  0

−1



0

|λ2 |



 → − µ1

− → µ2 

− → V

−1 − → µ2 

− →

V .

102

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

Definition 3.1. The matrix





 |λ1 | − →  µ2

→ |N | = − µ1

0



0

|λ2 |

−1

 → − µ1

− → µ2 

previously described is called the absolute value of the matrix N. For mesh adaptation purposes, let err be the local approximation error on a triangle. Assume that the error estimate is

− →

written in terms of the Hessian matrix H , a vector V and a constant C , in a given norm ∥ · ∥:

− → − → ∥err ∥ ≤ C ⟨ V , H V ⟩. Consider δ , the maximal local error expected i.e.,

− →

− →

C ⟨ V , H V ⟩ ≤ δ. Let hmax be the maximal size of edge permitted in the mesh, let hmin be the minimal edge size permitted and let RH be the matrix of eigenvectors of H . The metric |H |, absolute value of H with respect to the mesh size bounds, the error estimates and δ is then given by −1 |H | = RH ΛH RH ,

(3.1)

where ΛH is the diagonal matrix of entries λi defined by

˜H

    C |λH 1 1 i | ˜λH , , 2 , 2 i = min max δ hmax hmin

(3.2)

where λH i is the eigenvalue associated with the ith eigenvector of H . Metric and ellipse. Consider the metric

 M=

− →

a

c /2

c /2

b

 ,

a > 0, b > 0, 4ab − c 2 > 0, P (xp , yp ), Q (x, y), xp , yp , x, y ∈ R.

∥PQ ∥M = 1 ⇐⇒ a(x − xp )2 + c (x − xp )(y − yp ) + b(y − yp )2 = 1. 2

By writing cos θ 2

a=

α2

+

sin θ 2

β2

,

b=

cos θ 2

β2

+

sin θ 2

α2

,

 c=

1

α2

−

1

β2



sin 2θ ,

(3.3)

which are bounded by: 1

α2

≤a≤

1

β2

,

1

α2

≤b≤

1

β2

,

1

α2

−

1

β2

≤c≤

1

β2

−

1

α2

,

(3.4)

we obtain:

((x − xp ) cos θ + (y − yp ) sin θ )2 ((y − yp ) cos θ − (x − xp ) sin θ )2 + = 1. α2 β2

(3.5)

(3.5) is the equation of an ellipse centered at P with semi-axes of length α > 0 and β > 0, α ≥ β , where θ is the angle between the semi-major axis of the ellipse and the x-axis of the R2 canonic base. 3.1. Introduction of the problem Thanks to Céa’s lemma, the approximation error of the finite element solution of an elliptic PDE is upper bounded by the interpolation error. So building an adapted mesh for the finite element solution of an elliptic PDE such that it satisfies locally a maximal approximation error prescribed, can be done by considering the interpolation error. Then, the problem of mesh adaptation for the finite element approximation of an elliptic PDE is solved with respect to the interpolation error as error estimates. Results described in Section 2 [17, Section 3.1], aim to bound the interpolation error of a function on a simplex with the norm of a vector inside a metric given by the D(l) derivative of the function. We characterize the anisotropic behavior of derivative tensors of the interpolation of a function u by the largest ellipse contained in the level curve of the polynomial for directional derivatives (error curve), assuming that the error curve is locally described by a curve representing a level set of the error. We suppose that we know an estimation of the local interpolation error of a function on a two-dimensional mesh

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

103

around a vertex P in each direction ξ , and we want to dispose the points Q of the triangulation around P such that the local interpolation error be upper bounded. Let Ep (ξ ) be an estimation of the interpolation error in the direction of the vector ξ . If the mesh size is less than or equal to ∥ξ ∥ in the direction of ξ then, the interpolation error will be less than Ep (ξ ). We want to have locally

|Ep (ξ )| ≤ E ,

ξ ∈ R2 .

(3.6)

Inequality (3.6) represents an area bounded by the level set

|Ep (ξ )| = E .

(3.7)

We want to dispose the points Q of the triangulation around the vertex P, inside the level set (3.7), as close as possible to that level set. This is equivalent to having the interpolation error evaluated at each vertex of the triangulation (going from the vertex P) less than or equal to its value on the level set. We also want to have a mesh of which edges are equal with respect to a metric field, i.e., all the points Q are equidistant from P in a local metric M . So they belong to the unit ball of the metric M , the biggest possible that could be contained inside the level set described by (3.7). We look for a metric which prescribes the mesh size around a vertex P such that the triangle’s vertices be as close as possible to the level set (3.7) while being inside the area delimited by that level set. How to build that metric? 3.2. Formulation We do not take care of the different norms used, the different constants of (2.2) and the bounds for element length scale in the mesh. We assume that these constants and the bounds for element length scale are inputs of the adaptation function, in the mesh adaptation procedure (it is the case in Freefem++). We assume that

|Ep (ξ )| = ∥(D(ℓ) u)(ξ , . . . , ξ )∥(K ) .

(3.8)

We look for a local metric M such that the right-hand side of inequality (2.2) be as close as possible to E in (3.7). Let F (x, y) be Ep (ξ ). Our problem becomes drawing an ellipse centered at P which has the maximum area and which is contained inside the given error curve, or looking for the parameters α, β and θ of that ellipse. Consider

α1 =

1

α2

,

α2 =

1

β2

,

we obtain from (3.3):

α2 =

 ( a − b) 2 + c 2 + ( a + b) 2

,

α1 =

 − (a − b)2 + c 2 + (a + b) 2

.

Maximizing the ellipse area, the parameters of which are α, β, θ , is the same as maximizing αβ or minimizing α1 α2 , which is equivalent to minimizing (4ab − c 2 ). We consider the problem in a numerical point of view and assume that the vertex P is at the origin (change of basis).   Consider a discretization of the error curve with n points Mi of position (xi , yi ), set Xi =

xi yi

. We have to find three reals

a, b and c that minimize (4ab − c 2 ) under the following constraints:

 b > 0, a > 0, ax2i + by2i + cxi yi ≥ 1,  4ab > c 2 .

0 ≤ i ≤ n − 1,

(3.9)

Without ‘‘good’’ properties (convexity, quasi-convexity) of the minimization function, we did not find any algorithm that guarantees a global minimum. However, the global minimum exists since the linear constraints of the second line of system (3.9) define a closed set. The resolution of the problem formulated in (3.9) needs the construction of the polyhedron of constraints where the minimum lies. As (a, b, c ) ∈ R3 , a naive resolution algorithm will be of O (n3 ) complexity, which can be an expensive part of a mesh adaptation process. To reduce that cost we propose to solve an approximated problem of which the algorithm is of O (nk ), k ≤ 2, complexity. Let M0 be the point of the error curve, the closest to the origin P of the required ellipse, and let Mi be another point of the error curve. We look for an ellipse which has the maximal area with respect to the constraints of (3.9). Let ϵ0 be the distance between the ellipse and M0 , and let ϵ be the distance between the ellipse and Mi (Fig. 1). The distance between the ellipse and a point Mj of the error curve is assumed to be the distance between Mj and the

−→

intersection point of the vector Mj P with the ellipse. Let Q be a point of the ellipse on the segment [P , Mi ] such that

−−→ ∥Mi Q ∥ = ϵ then, − →

PQ = X =

 1−

ϵ ∥Xi ∥



Xi .

104

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

Fig. 1. Parameters of the ellipse with respect to the error curve. The parameters ϵ0 and ϵ , distances between the error curve and the ellipse respectively at points M0 and Mi .

Set ri = ∥Xi ∥, the problem (3.9) becomes: Find three real numbers a, b and c that minimize (4ab − c 2 ) under the constraints

a > 0, b > 0,  2    r0    , ax20 + by20 + cx0 y0 =   r0 − ϵ0     ri

2

ax2i + by2i + cxi yi =    ri − ϵ    2 2  ax + by + cx y ≥ 1 ,  j j j   j 4ab > c 2 .

,

1 ≤ i ≤ n − 1,

(3.10)

1 ≤ j ≤ n − 1, j ̸= i,

Numerically, the constraint 4ab > c 2 needs a threshold T > 0 for which a solution of 4ab − c 2 − T will be considered as a minimizer of 4ab − c 2 . Minimization of 4ab − c 2 . 4ab − c 2 is four times the product of the eigenvalues of the matrix M . Let Rmax be the Euclidean norm of the vector defined by a point Mmax of the error curve and P, such that Mmax be the furthest away from the center P of the ellipse. As α ≤ Rmax and β ≤ r0 − ϵ0 , we obtain:

 T =

2

2 Rmax (r0 − ϵ0 )

.

4. Algebraic resolution We examine system (3.10): If x0 = 0 then, b=

1

(|y0 | − ϵ0 )2

.

• If xi = 0 then, b=

1

(|yi | − ϵ)2

,

and one must have:

ϵ = |yi | − |y0 | + ϵ0 . • If yi = 0 and yj = 0, one must have:

(4.1)

ϵ ≥ |xi | − |xj |. • If xi y0 − yi x0 = 0 and xi ̸= 0, one must have:   x0 (r0 − ϵ0 ) ϵ = ri 1 − .

(4.2)

xi r0

(4.3)

Remark 4.1. One can observe that if a coefficient of system (3.10) equals zero or if the points Mi , M0 and the origin P are aligned, the choice of the parameters ϵ and ϵ0 is not arbitrary. Those parameters must satisfy (4.1), (4.2) or (4.3). However, it will add new constraints to the problem.

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

105

Thanks to Remark 4.1, to preserve the arbitrary choice of the parameters ϵ and ϵ0 , we do the following: The points (of the error curve) for which the coefficients of system (3.10) equal zero, are up to five (the four half-axes and eventually a point on the axis (P , M0 )). So, if there exists points for which the coefficients of system (3.10) vanish, we turn these points out of their position from a small angle, with respect to the center P of the ellipse. Doing that does not sensitively modify the problem nor its solution. We assume that x0 , xi , xj , y0 , yi , yj are all not zero for all i and for all j, and also assume that xi y0 − yi x0 is not zero for all i. From the two equations of system (3.10), we extract c and a as functions of b. Set:

∆=

x0 (xi y0 − x0 yi )

2

r0

−

r0 − ϵ0



yi

γ =− τ=



xi

+

x0 (xi y0 − x0 yi ) r0 − ϵ0 yi y0 x0 yi + y0 xi xi x0

Cj = 1 −

,

λ=

,

xi x0

xj (yj xi − yi xj ) x0 (xi y0 − x0 yi )



xi ( xi y 0 − x0 y i )

2

r0



x0

ri − ϵ



y0

2

ri

xi (xi y0 − x0 yi )

,

2

ri ri − ϵ

,

βj = y2j + x2j τ − xj yj λ, 2

r0

+

r0 − ϵ0

xj (yj x0 − xj y0 )



xi (xi y0 − x0 yi )

ri ri − ϵ

2

,

then, c = ∆ − λ b,

a = γ + τ b.

(4.4)

The inequality before the last one of (3.10) becomes

βj b ≥ Cj ,

1 ≤ j ≤ n − 1, j ̸= i.

(4.5)

So, the minimization of 4ab − c 2 becomes the minimization of G(b) = (4τ − λ2 )b2 + 2(2γ + λ∆)b − ∆2 − T ,

(4.6)

the roots of which are:

 b1 =

2  

1

b2 =

2  

1

2

x0 ri

2

xi r0

2

r0 − ϵ 0

 +

ri − ϵ

xi y0 − x0 yi

 +

ri − ϵ

xi y0 − x0 yi



x0 ri

xi r0

2

r0 − ϵ 0

 +D , 

(4.7)

−D ,

where

D=

2|xi x0 |

(r0 − ϵ0 )





ri r0

2

ri − ϵ

 +

x0 yi − y0 xi Rmax

2

.

The inequalities (3.4), (4.5) and the following bounds

(r0 − ϵ0 )2 Rmax

≤ β ≤ r0 − ϵ0 ≤ α ≤ Rmax ,

give bmin and bmax such that b ∈ [bmin , bmax ]. We compare G(bmin ) and G(bmax ) with respect to b1 and b2 in (4.7), in order to determine the value of b which minimizes G in (4.6). 5. Algorithm The metric (ellipse) is computed as follows: Algorithm 1 Metric computation. 1: 2: 3:

Discretize the error curve around vertex P into n points, determine the point M0 of the curve which is the closest to the origin P and the point Mmax of the curve which is the furthest away from P. Move by a small angle the points for which the coefficients of system (3.10) vanish, with respect to the center P of the ellipse. For all fixed Mi , switch the constraints of (3.10) for all the Mj (0 ̸= j ̸= i) and determine the optimal value of b in (4.6) which by (4.4), gives the values of a and c. Conserve the triplet (a, b, c ) which maximizes the area of the ellipse.

106

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

Fig. 2. Convergence of the Algorithm. An ellipse of parameters θ = π4 , α = 3, β = 2, n = 159. Left: ϵ0 = 10−1 . Right: ϵ0 = 10−4 . We have obtained a convergence of the approximation with respect to the choice of ϵ0 .

Fig. 3. Metrics for arbitrary error curve. Left: E = 59, ϵ0 optimal = 10−5 , n = 285, F (x, y) = −0.2y3 x2 − yx4 + 0.2y3 − 9x3 − 0.5xy2 + 0.005x2 − 0.1y2 . Right: F (x, y) = 0.05x2 − 0.01y2 + 0.05xy, E = 1, ϵ0 optimal = 0.001, n = 304. The error curve of the left picture is unspecified (not a particular higher order derivative) so, one can always compute a metric for any local error estimates expressed in terms of a curve.

5.1. Algorithm order Algorithm 1 is globally quasi-linear because for each fixed Mi only one (the first) unsatisfied constraint on the Mj implies an interruption of the Mj ’s loop and switch to the next Mi+1 . The points which satisfy the constraints are only those in the regions where the required ellipse is tangent to the error curve, and their number is few in general. The unique case for which for each fixed Mi we test all the other Mj , is when the error curve is already an ellipse. In this case, the algorithm order is O (n2 ). But the error curve describing the l-derivatives is an ellipse only if it describes the error estimates given by the Hessian matrix, the case of which we do not need Algorithm 1, since the construction of the metric is known. 5.2. Convergence of the approximation It is natural to think that M0 will not be the point of the error curve that is the furthest away from the ellipse. So we can look for ϵ0 such that the solution of our approached problem converges to an optimal solution. We have an upper bound of ϵ0 : The disc area of radius r0 is less than the ellipse area of parameters r0 − ϵ0 and Rmax so,

  r0 ϵ0 < r0 1 − . Rmax

We do not have any absolute criterion for choosing ϵ0 . Nevertheless, we vary its values within a small number of iterations inside the code between 10−5 and its upper bound, and thus optimize the approximation. Figs. 2 and 3 show examples of the ellipses computed, with the obtained optimal values of ϵ0 (values of ϵ0 of which the area of the ellipse is maximal). However, the solution obtained with the value of ϵ0 = 10−3 does not differ significantly from the optimal solution obtained when we vary ϵ0 . So, we generally take ϵ0 = 10−3 . 6. Numerical experiments Computation of the error curve. Algorithm 1 needs a discretization of the error curve around a vertex P into n points. One can observe that the lth order derivatives (D(ℓ) u)(ξ , . . . , ξ ) of a function u in the direction ξ which describes the error (3.8), is

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

107

Fig. 4. Convergence of the Algorithm. E = 1, F (x, y) = 0.5x2 + 0.5y2 + 2xy, left: n = 228, right: n = 288. A better approximation of the error curve gives a better approximation of the ellipse.

Fig. 5. Case of non-invertible Hessian matrix. E = 1, n = 234 F (x, y) = x2 + y2 − 2xy. The Hessian matrix is non-invertible and we have chosen a different maximal authorized edge size of the mesh. However, the choice of the same maximal authorized edge size for both analytical and approximated ellipse gives the same result.

a homogeneous polynomial of degree l. One has to find n points (xi , yi ) coordinates of vector ξi in the Cartesian coordinates system of origin P, and which satisfy (3.7). Let r be a real number and ϕ ∈ [0, 2π ] such that

(x, y) = r (cos ϕ, sin ϕ).

Ep (x, y) = r l Ep (cos ϕ, sin ϕ),

since Ep (x, y) is a homogeneous polynomial of degree l. For each value ϕi of a uniform discretization of [0, 2π ] and by (3.7), one obtains (xi , yi ) = ri (cos ϕi , sin ϕi ), where

 ri =

E

1/l

Ep (cos ϕi , sin ϕi )

.

Many algorithms [11–14] used to generate anisotropic adapted meshes suitable for piecewise linear Lagrange finite element solution of a PDE, use the Hessian matrix of the computed solution as error estimates. As the Hessian matrix H is a symmetric matrix, its absolute value (Definition 3.1) with respect to hmin , hmax , the maximal local error expected and the error estimates defined in (3.1) can be used to define a metric. We denote by analytical metric/ellipse, the metric/ellipse given by the absolute value of the Hessian matrix defined in (3.1). In some particular cases of second order and third order derivatives, one can derive exact metrics [18]. 6.1. Some examples We apply Algorithm 1 with F (x, y) = ⟨(x, y), H · (x, y)⟩, we obtain comparative results with the analytical metric, Figs. 4 and 5. In the example of Fig. 4 for the case H invertible, one can observe that the obtained ellipse is close to the analytical one when the discretization of the error curve is ‘‘good’’, particularly when the number n of discretization points is high enough. The equation of the analytical ellipse of the example of Fig. 4 is x2 + y2 + xy = 1.

108

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

Fig. 6. Function f (x, y) = tanh(2(sin(5y) − 2x)) + yx2 + y3 .

The example of Fig. 5 shows the case H non-invertible where hmax has been arbitrarily chosen in the analytical case. The equation of the analytical ellipse of the example of Fig. 5 is 0.5



1 h2max

    1 + 2 (x2 + y2 ) + 2 2 − 2 xy = 1. hmax

We just need to vary the value of hmax to coincide with the approached solution. 6.2. Mesh adaptation with Freefem++ Lots of PDE software programs like Freefem++ [19] use metrics to build meshes, the edge sizes of which are equal with respect to the metric field. We suppose that once a metric field on a mesh is known, we can build an adapted mesh for that metric field [24, Section 21.4], [26, Section 20.3]. If we know the metric value on each vertex of the mesh for example, the method described in [26] modifies the mesh by evaluating each edge length in the metric field (control space) using an interpolation of the metrics defined on the two vertices of the edge. When an edge length exceeds the ‘‘unity’’, a subdivision of that edge is done in the control space. A loop consisting of resolution and adaptation: Algorithm 2, is executed for a small number of iterations mmax . This method is used in the software Freefem++ with the mesh generator Bamg [27]. Algorithm 1 has been coded in C++ and integrated into Freefem++as a library (MetricKuate). Algorithm 2 Mesh adaptation loop. 1: 2: 3: 4: 5: 6:

Mesh generation → T0 Solve the problem on Tm Compute error estimates and a metric field M on Tm Modify Tm with respect to M m←m+1 goto 2 until m = mmax .

6.2.1. Interpolation of a function We present the results of mesh adaptation for the following function f , Fig. 6: f (x, y) = tanh(2(sin(5y) − 2x)) + yx2 + y3 . For the Pk Lagrange finite element interpolation of f , the error function err k at the vertex P in the direction

  x0 y0

is given

by the kth order derivatives of the function f times a coefficient cf .

 errPk

(x0 , y0 ) = cf

 (k + 1)! ∂ (k+1) f k1 k2 (P )x0 y0 , k1 !k2 ! ∂ xk1 ∂ yk2 =k+1

 k1 +k2

k1 , k2 , k ∈ N, k ≥ 1.

(6.1)

The error function is written as a homogeneous polynomial of degree (k + 1) of the coordinates of the local unit ball, as is shown in [17]. We use Algorithm 1 with E = 1, to compute the metric on each vertex of the mesh, and the mesh is adapted with the obtained metric. We compute a mesh adaptation processes for the function f , on a square [−1, 1] × [−1, 1], using the kth, (2 ≤ k ≤ 5), order derivatives as error estimates and associated metrics computed with Algorithm 1.

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

109

Fig. 7. Adapted meshes for function f of Fig. 6 and the maximal error level around 10−5 , metrics computed with Algorithm 1 for the (3.8) error estimator in terms of the third order derivatives.

Fig. 8. Adapted meshes for function f of Fig. 6 and the maximal error level around 10−5 , method based on the Hessian matrix.

Comparison of the second and third order derivatives error estimates. In the code, we compute backward the real cf such that we have the same error level on each final mesh, for the Hessian error estimates and for the third order derivatives error estimates. Table 1 represents the error and the number of vertices of the adapted meshes for different iterations. N denotes the number of vertices of the mesh. The index 2 is for the mesh obtained using the third order derivatives error estimates. The index 1 is for the mesh obtained using the Hessian matrix error estimates. ‘‘error min’’ and ‘‘error max’’ are respectively the minimal error and the maximal error computed on the triangles of the mesh. The error is computed at the barycenter of each triangle, as the absolute value of the Π0 interpolation of f − Π2 f , where Πk f is the polynomial interpolation of f of degree k. The behavior of the results presented in Table 1 does not change significantly when one computes the error with f − Π1 f , or when one measures the error using classical norms. We can observe on Figs. 7, 8 and on Table 1 that for the same maximal error level around 10−5 , Algorithm 1 applied to the third order derivatives recovers correctly the directions of anisotropy and uses less elements than the classical mesh adaptation method with metrics based on the Hessian matrix. We see that the convergence is obtained after three iterations,

110

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115 Table 1 Error and number of vertices of the adapted meshes for different iterations. The index 2 is for the mesh obtained using the third order derivatives. The index 1 is for the mesh obtained using the Hessian matrix. The error is: error = |Π0 (f − Π2 f )|, taken at the barycenter of each triangle. Iteration

Error min2

Error min1

Error max2

Error max1

N2

N1

1 2 3 4 5 6 7 8 9

1.20059e−05 3.11655e−09 3.62071e−09 3.69877e−09 1.10912e−09 7.05385e−09 1.89551e−09 5.05435e−10 9.78865e−10

4.22866e−06 2.48466e−09 1.52523e−11 2.19181e−10 4.76638e−10 3.61698e−10 9.00541e−12 1.42422e−10 9.68075e−10

0.356521 0.0043082 9.67516e−05 7.95773e−05 7.79361e−05 7.77929e−05 8.01738e−05 7.75686e−05 8.02579e−05

0.358384 0.00437699 9.9182e−05 0.000103307 0.000106652 8.02886e−05 7.80091e−05 8.22038e−05 7.7833e−05

2050 5194 5198 5194 5264 5200 5200 5200 5218

3 923 17 235 15 381 15 210 15 188 15 190 15 186 15 186 15 180

Table 2 Errors of the P2 finite element solution of (6.5). Left: adapted meshes for each iteration, metrics computed with respect to the third order derivatives obtained from uh . Right: uniform and non-adapted meshes composed of right-angled and quasi-isosceles triangles. The error is significantly reduced from the uniform to the adapted meshes. Iteration

1 2 3 4 5 6

Adapted meshes for the P2 approximation

Non-adapted meshes for the P2 approximation.

NT

≀ ≀ u − uh ≀≀0

≀ ≀ ∇(u − uh ) ≀≀0

NT

≀ ≀ u − uh ≀≀0

≀ ≀ ∇(u − uh ) ≀≀0

202 153 129 125 125 125

0.00912801 0.00623116 0.0109628 0.0105778 0.0101554 0.0106713

0.274678 0.237684 0.243304 0.242451 0.229079 0.191755

200 180 162 160 128 126

0.400887 0.453795 0.503353 0.522699 0.630091 0.649742

1.9092 2.01553 2.11673 2.12949 2.32936 2.32995

Fig. 9. Adapted meshes for function f of Fig. 6 for fourth and fifth order derivatives error estimates, metrics computed with Algorithm 1. Left: errors estimates in terms of fourth order derivatives, 2894 vertices. Right: errors estimates in terms of fifth order derivatives, 2790 vertices.

Table 1. We can observe on Fig. 8 that the corners of the domain are meshed with small size elements for the Hessian based method. Those wrong effects are not observed with the third order derivatives on Fig. 7. One can also observe that the third order derivatives capture differently the middle layer of the mesh, than the second order derivatives. Directions of anisotropy. Fig. 9 shows examples of adapted meshes for function f for the fourth and for the fifth order derivatives error estimates. One can observe that we obtain the right directions of anisotropy as for the cases of second and third order derivatives error estimates, Figs. 7 and 8. Convergence rate of the interpolation error. We have also compared the interpolation errors on adapted meshes in terms of mesh complexity, for function f with metrics computed using our Algorithm. We use exact derivatives of order k + 1

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

111

Fig. 10. Convergence of the interpolation on adapted meshes in terms of the mesh complexity (number of nodes N). The meshes are adapted with metrics computed using our Algorithm with error given by the exact derivatives of order k + 1 for the Pk finite element interpolation of function f . The convergence rate of the interpolation error with respect to the mesh adaptation process increases with the degree of the polynomial interpolation, and their values are −(k+1) for the Pk finite element interpolation. approximately about the optimal bounds proposed in [18]: 2

for the Pk finite element interpolation of f and compute the metrics from the error given by those derivatives (6.1), with Algorithm 1. Fig. 10 represents the convergence rates of the interpolation error in terms of the number of nodes N of the mesh, computed on sequences of adapted meshes. One can observe that the convergence rate of the interpolation error with respect to the mesh adaptation process increases with the degree of the polynomial interpolation, and the convergence rates −(k+1) observed are approximately about the optimal bounds proposed in [18]: for the Pk finite element interpolation. 2 6.2.2. Approximation of a function We use the following numerical quadrature for the computation of the error: Let Ω be a two-dimensional polygonal domain. Let v ∈ L2 (Ω ) be a function. For any triangulation Th of Ω , we denote by ≀≀v≀≀0 the numerical quadrature of ∥v∥0 := ∥v∥L2 (Ω ) (the standard L2 norm of v on Ω ). Let xi , i = 1, 2, . . . , nt , be some nodes of a triangle T and denote vi := v(xi ).

≀≀v≀≀20;T := |T |

nt  i=1

ωi vi2 ,

≀≀v≀≀20 :=



≀≀v≀≀20;T ,

T ⊂Ω

where ωi ∈ R, xi and nt are chosen such that (6.2) be exact for a polynomial of degree five. We next give a short description of the least-square method used for the computation of derivatives.

(6.2)

112

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

Fig. 11. The patch Txn of triangles used for the computation of second and first derivatives of a function at vertex xn .

The least squares method. The approximation of Hv (xn ), Hessian matrix of function v at vertex xn is obtained from a second order Taylor expansion of v(x), in a neighborhood of xn : 1 − → − → − → → v(x) = v(xn ) + − xn x · ∇v(xn ) + ⟨xn x, Hv (xn )xn x⟩ + O(∥xn x∥3 ), 2

(6.3)

where ⟨, ⟩ denotes the scalar product of R2 . We compute simultaneously Hv (xn ) and ∇v(xn ) by writing Eq. (6.3) at as many vertices xi as Hv (xn ) and ∇v(xn ) have unknowns (#n = 5), in a neighborhood of xn . However, the number of vertices surrounding xn (first ring) may be less than #n. So, we add another ring consisting of triangles of the neighbors vertices of the previous ring, such that the patch Txn associated with vertex xn contains at least #n vertices, Fig. 11. Then, one has to solve a system with more equations than unknowns: AY = b,

(6.4)

where Y is the vector containing the unknowns (of the Hessian matrix and the gradient of v ). System (6.4) is solved in the sense of the least squares approximation, which minimizes the distance between AY and b by minimizing the square of the Euclidean norm of their difference. The minimizer is the solution of the following system: t

AAY = t Ab.

The computation of the third order derivatives is done using a third order Taylor expansion instead of (6.3). We consider the Poisson equation with Dirichlet boundary conditions in a two-dimensional polygonal domain Ω ,

− ∆u = g in Ω and u = ue on Γ := ∂ Ω .

(6.5)

We solve the weak form of problem (6.5) for the following value of the exact solution: ue (x, y) = exp(−x/ζ ) + exp(−y/ζ ),

ζ = 5 · 10−2 .

Let uh be the finite element solution of the weak form of problem (6.5). The exact solution is taken with non-homogeneous Dirichlet boundary conditions, and the appropriate g derived from ue . The domain Ω is the square [0, 5] × [0, 5]. We have computed the second and the third order derivatives with the least squares method, from uh .

• (i) The P2 finite element solution of (6.5) is computed and the mesh is adapted with respect to the third order derivatives error estimates, with metrics computed using Algorithm 1.

• (ii) The P1 finite element solution of (6.5) is computed and the mesh is adapted with respect to the second order derivatives error estimates, with associated analytical metrics. Let NT be the number of triangles of the mesh. We have done the mesh adaptation processes (i) and (ii) such that the number of triangles of the two computations be about the same. One can observe on Fig. 12 that the right direction of anisotropy is obtained with metrics computed using Algorithm 1, with respect to the third order derivatives (from uh ) computed with the least squares method. We have compared the error of the P2 finite element solution of (6.5), for adapted meshes and uniform (non-adapted) meshes composed of right-angled and quasi-isosceles triangles. One can observe on Table 2 that for about the same number of triangles, the adapted meshes reduce significantly the error obtained on the uniform mesh, which confirms the theoretical results expected. Fig. 13 shows adapted meshes obtained from the P1 and from the P2 finite element solution of (6.5) for about the same error. The P1 finite element approximation needs much more triangles for the same error (Table 3 left) or increases the error for about the same number of triangles (Table 3 right) than the P2 finite element approximation (Table 2 left). We have also compared the convergence rate of the approximation on the adapted meshes, Fig. 14. Thus, we have adapted the mesh on many cases of different mesh complexities, for both the P1 and P2 finite element solution of (6.5). The error

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

113

Fig. 12. Adapted meshes with about the same number of triangles at convergence. Left: P2 finite element solution of (6.5), metrics computed with respect to the third order derivatives obtained from the associated uh . Right: P1 finite element solution of (6.5), analytical metrics with respect to the second order derivatives obtained from the associated uh .

Fig. 13. Adapted meshes with about the same error (L2 -norm) around 10−2 at convergence. Left: P2 finite element solution of (6.5), metrics computed with respect to the third order derivatives obtained from the associated uh . Right: P1 finite element solution of (6.5), analytical metrics with respect to the second order derivatives obtained from the associated uh .

is computed using (6.2) and the Lp , 1 ≤ p ≤ 4, norms are presented, Fig. 14: The left picture represents the results for the P2 finite element solution of (6.5), and the right picture represents the results for the P1 finite element solution of (6.5). One can observe that there is a significant gain on the convergence rate by the use of third order derivatives error estimates, compared to the use of second order derivatives error estimates. 7. Conclusion We have proposed an algorithm of which the complexity is less than a quadratic complexity, for the approximation of metrics for any given error estimation represented by a curve around a mesh vertex, in dimension two of the space. The main purpose of this work is to perform anisotropic mesh adaptation for numerical simulations with Lagrange finite elements approximation of degree k, k > 1. The numerical tests of mesh adaptation done using metrics computed with our Algorithm with respect to the third order, fourth order and fifth order derivatives as interpolation error estimates, show that we obtain the right directions of anisotropy, and better results for the error estimates as theoretically expected.

114

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

Fig. 14. Comparison of the convergence rate of the approximation on adapted meshes in terms of the mesh complexity (number of nodes N). Left: errors of the P2 finite element solution of (6.5), metrics computed using our Algorithm with respect to the third order derivatives obtained from the associated uh . Right: errors of the P1 finite element solution of (6.5), analytical metrics with respect to the second order derivatives obtained from the associated uh . There is a significant increasing of the convergence rate by the use of third order derivatives error estimates, compared to the use of second order derivatives error estimates, for the mesh adaptation process. Table 3 Errors of the P1 finite element solution of (6.5). Adapted meshes for each iteration, metrics computed with respect to the second order derivatives obtained from uh . The number of triangles on the left table or the error on the right table is significantly increased compared to the P2 finite element approximation, Table 2 left. Adapted meshes for the P1 approximation.

Adapted meshes for the P1 approximation.

Iteration

NT

≀ ≀ u − uh ≀≀0

≀ ≀ ∇(u − uh ) ≀≀0

Iteration

NT

≀ ≀ u − uh ≀≀0

≀≀∇(u − uh ) ≀≀0

1 2 3 4 5 6

706 550 384 366 366 368

0.0486665 0.012564 0.0111353 0.0117546 0.0112295 0.0110409

2.42084 1.07564 0.983516 0.987868 0.966913 0.955524

1 2 3 4 5 6

196 145 132 128 126 126

0.179948 0.146021 0.104809 0.0989654 0.104739 0.100723

4.36019 3.73805 3.32904 3.31696 3.34238 3.28807

A future work could be an investigation on other types of error estimates for simulations of problems involving conforming triangular meshes, of which the level set of these error estimates may be locally expressed as a curve, since our Algorithm computes metrics for any error estimation which is locally represented by a curve. References [1] P.S. Heckbert, M. Garland, Optimal triangulation and quadric-based surface simplification, Comput. Geom. 14 (1–3) (1999) 49–65. 11. [2] W. Huang, Mathematical principles of anisotropic mesh adaptation, Commun. Comput. Phys. 1 (2006) 276–310. [3] T. Coupez, Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing, J. Comput. Phys. 230 (7) (2011) 2391–2405. [4] A. Loseille, F. Alauzet, Continuous mesh framework part II: validations and applications, SIAM J. Numer. Anal. 49 (1) (2011) 61–86. [5] J.C. Aguilar, J.B. Goodman, Anisotropic mesh refinement for finite element methods based on error reduction, J. Comput. Appl. Math. 193 (2) (2006) 497–515. [6] M. Randrianarivony, Anisotropic finite elements for the stokes problem: a posteriori error estimator and adaptive mesh, J. Comput. Appl. Math. 169 (2) (2004) 255–275. 8. [7] J.-F. Remacle, X. Li, M.S. Shephard, J.E. Flaherty, Anisotropic adaptive simulation of transient flows using discontinuous galerkin methods, Internat. J. Numer. Methods Engrg. 62 (7) (2005) 899–923. [8] X. Jiao, A. Colombi, X. Ni, J. Hart, Anisotropic mesh adaptation for evolving triangulated surfaces, Eng. Comput. 26 (2010) 363–376. [9] G. Kunert, R. Verfürth, Edge residuals dominate a posteriori error estimates for linear finite element methods on anisotropic triangular and tetrahedral meshes, Numer. Math. 86 (2000) 283–303. [10] E. Creusé, G. Kunert, S. Nicaise, A posteriori error estimation for the Stokes problem: anisotropic and isotropic discretizations, Math. Models Methods Appl. Sci. 14 (9) (2004) 1297–1341. [11] V. Dolej˜sí, Anisotropic mesh adaptation for finite volume and finite element methods on triangular meshes, Comput. Vis. Sci. 1 (1998) 165–178. [12] F. Alauzet, P. Frey, Anisotropic mesh adaptation for CFD computations, Comput. Methods Appl. Mech. Engrg. 194 (48–49) (2005) 5068–5082. [13] W. Cao, On the error of linear interpolation and the orientation, aspect ratio, and internal angles of a triangle, SIAM J. Numer. Anal. 43 (1) (2005) 19–40 (electronic). [14] E. D’Azevedo, R. Simpson, On optimal triangular meshes for minimizing the gradient error, Numer. Math. 59 (1991) 321–348. [15] M.J. Castro-Díaz, F. Hecht, B. Mohammadi, O. Pironneau, Anisotropic unstructured mesh adaption for flow simulations, Internat. J. Numer. Methods Fluids 25 (4) (1997) 475–491. [16] W. Cao, An interpolation error estimates on anisotropic meshes in Rn and optimal metrics for mesh refinement, SIAM J. Numer. Anal. 45 (6) (2007) 2368–2391.

F. Hecht, R. Kuate / Journal of Computational and Applied Mathematics 258 (2014) 99–115

115

[17] F. Hecht, Estimation d’erreur d’interpolation anisotropes pour des éléments finis de Lagrange, in: Numerical Analysis and Scientific Computing for Partial Differential Equations and their Challenging Applications, CIMNE, ISBN: 978-84-96736-41-2, 2008, pp. 108–120. [18] J.-M. Mirebeau, Optimal meshes for finite elements of arbitrary order, Constr. Approx. 32 (2010) 339–383. [19] F. Hecht, O. Pironneau, Freefem++, language for finite element method and PDEs, Laboratoire Jacques-Louis Lions, Université Pierre et Marie Curie, http://www.freefem.org/ff++/. [20] C. Bernardi, Y. Maday, F. Rapetti, Discrétisations variationnelles de problèmes aux limites elliptiques, in: Mathématiques & Applications (Berlin) (Mathematics & Applications), vol. 45, Springer-Verlag, Berlin, 2004. [21] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, in: Classics in Applied Mathematics, vol. 40, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002. Reprint of the 1978 original (North-Holland, Amsterdam; MR0520174 (58 #25001)). [22] P.G. Ciarlet, J.-L. Lions (Eds.), Handbook of Numerical Analysis, Vol. II, North-Holland, Amsterdam, 1991. Finite element methods, Part 1. [23] M. Picasso, An anisotropic error indicator based on Zienkiewicz–Zhu error estimator: application to elliptic and parabolic problems, SIAM J. Sci. Comput. 24 (4) (2003) 1328–1355 (electronic). [24] P. Frey, P.-L. George, Maillages: Applications aux Elements Finis, Hermès Science, Paris, 1999. [25] R. Kuate, Anisotropic metrics for finite element meshes using a posteriori error estimates: Poisson and stokes equations, Eng. Comput. (2012). [26] F. Hecht, P.L. Georges, Nonisotropic Grids, CRC Press, 1999. [27] F. Hecht, The mesh adapting software: bamg, INRIA report, 1998. http://www-c.inria.fr/gamma/cdrom/www/bamg/eng.htm.