Table Of Content

Greedy bisection generates optimally adapted triangulations Jean-Marie Mirebeau and Albert Cohen January 10, 2011 1 1 0 2 Abstract n Westudythepropertiesofasimplegreedyalgorithmintroducedin[8]forthegenerationofdata- a adaptedanisotropictriangulations. Givenafunctionf,thealgorithmproducesnestedtriangulations J TN andcorrespondingpiecewisepolynomialapproximationsfN off. Therefinementprocedurepicks 7 thetriangle which maximizes thelocal Lp approximation error, and bisectsit in adirection which is chosen sotominimize thiserrorat thenextstep. Westudytheapproximation errorin theLp norm ] A when the algorithm is applied to C2 functions with piecewise linear approximations. We prove that asthealgorithm progresses, thetriangles tendtoadopt an optimalaspect ratio whichis dictatedby N the local hessian of f. For convex functions, we also prove that the adaptive triangulations satisfy h. the convergence bound kf −fNkLp ≤CN−1kpdet(d2f)kLτ with τ1 := 1p +1, which is known to be t asymptotically optimal among all possible triangulations. a m [ 1 Introduction 1 v In finite element approximation, a classical and important distinction is made between uniform and 2 adaptive methods. In the first case all the elements which constitute the mesh have comparable shape 5 and size, while these attributes are allowed to vary strongly in the second case. An important feature of 4 adaptive methods is the fact that the mesh is not fixed in advance but rather tailored to the properties 1 of the function f to be approximated. Since the function approximating f is not picked from a fixed . 1 linearspace,adaptivefiniteelementscanbeconsideredasaninstanceofnon-linearapproximation. Other 0 instances include approximation by rational functions, or by N-term linear combinations of a basis or 1 dictionary. We refer to [9] for a general survey on non-linear approximation. 1 : In this paper, we focus our interest on piecewise linear finite element functions defined over tri- v angulations of a bidimensional polygonal domain Ω IR2. Given a triangulation we denote by i ⊂ T X V := v s.t. v Π , T the associated finite element space. The norm in which we measure the T 1 r apTprox{imation|err∈oristhe∈LpTn}ormfor1 p andwethereforedonotrequirethatthetriangulations a ≤ ≤∞ areconformingandthat the functions of V arecontinuousbetweentriangles. For a givenfunctionf we T define e (f) := inf inf f g , N Lp Lp #( ) Ng VT k − k T ≤ ∈ the best approximation error of f when using at most N elements. In adaptive finite element approximation, critical questions are: 1. Given a function f and a number N >0, how can we characterize the optimal mesh for f with N elements corresponding to the above defined best approximation error. 2. What quantitative estimates are available for the best approximation error e (f) ? Such esti- N Lp mates should involve the derivatives of f in a different way than for non-adaptive meshes. 3. Can we build by a simple algorithmic procedure a mesh of cardinality N and a finite element N T function f V such that f f is comparable to e (f) ? N ∈ TN k − NkLp N Lp 1 While the optimal mesh is usually difficult to characterize exactly, it should satisfy two intuitively desirable features: (i) the triangulationshould equidistribute the local approximationerror between each triangle and (ii) the aspect ratio of a triangle T should be isotropic with respect to a distorted metric inducedbythelocalvalueofthehessiand2f onT (andthereforeanisotropicinthesenseoftheeuclidean metric). Under such prescriptions on the mesh, quantitative errorestimates have recently been obtained in [7, 1] when f is a C2 function. These estimates are of the form 1 1 limsupNe (f) C det(d2f) , = +1, (1.1) N Lp Lτ ≤ k | |k τ p N →∞ p where det(d2f) is the determinant of the 2 2 hessian matrix. For a convex C2 function f this estimate × has been proved to be asymptotically optimal in [7], in the following sense liminf Ne (f) c det(d2f) . (1.2) N Lp Lτ N + ≥ k | |k → ∞ p Theconvexityassumptioncanactuallybereplacedbyamildassumptiononthesequenceoftriangulations which is used for the approximation of f: a sequence ( ) is said to be admissible if #( ) N TN N≥N0 TN ≤ and sup N1/2 max diam(T) < . N≥N0(cid:18) T∈TN (cid:19) ∞ Then it is proved in [10], that for any admissible sequence and any C2 function f, one has liminf N inf f g c det(d2f) . (1.3) Lp Lτ N→+∞ g∈VTN k − k ≥ k | |k p The admissibility assumption is not a severe limitation for an upper estimate of the error since it is also proved that for all ε>0, there exist an admissible sequence such that limsupN inf f g C det(d2f) +ε. (1.4) Lp Lτ N→+∞ g∈VTN k − k ≤ kp| |k We also refer to [10] for a generalization of such upper and lower estimates to higher order elements. Fromthe computationalviewpoint,acommonlyusedstrategyfordesigninganoptimalmeshconsists therefore in evaluating the hessian d2f and imposing that each triangle of the mesh is isotropic with respect to a metric which is properly related to its local value. We refer in particular to [3] where this programis executed using Delaunay mesh generation techniques. While these algorithms fastly produce anisotropic meshes which are naturally adapted to the approximated function, they suffer from two intrinsic limitations: 1. They use the data of d2f, and therefore do not apply to non-smooth or noisy functions. 2. They are non-hierarchical: for N >M, the triangulation is not a refinement of . N M T T In [8], an alternate strategy was proposed for the design of adaptive hierarchical meshes, based on a simple greedy algorithm: starting from an initial triangulation , the algorithm picks the triangle TN0 T with the largest local Lp error. This triangle is then bisected from the mid-point of one of its N ∈ T edgestotheoppositevertex. Thechoiceoftheedgeamongthethreeoptionsistheonethatminimizesthe newapproximationerrorafterbisection. ThealgorithmcanbeappliedtoanyLpfunction,smoothornot, in the context of piecewise polynomial approximation of any given order. In the case of piecewise linear approximation, numerical experiments in [8] indicate that this elementary strategy generates triangles with an optimal aspect ratio and approximations f V such that f f satisfies the same N ∈ TN k − NkLp estimate as e (f) in (1.1). N Lp Thegoalofthispaperistosupporttheseexperimentalobservationsbyarigorousanalysis. Ourpaper is organized as follows: In 2,weintroducenotationswhichareusedthroughoutthepaperandcollectsomeavailableapprox- § imation theory results for piecewise linear finite elements, making the distinction between (i) uniform, (ii)adaptiveisotropicand(iii)adaptiveanisotropictriangulations. Inthe lastcase,whichis inthe scope 2 ofthis paper,weintroduceameasureofnon-degeneracyofatriangleT withrespecttoaquadraticform. Weshowthattheoptimalerrorestimate(1.1)ismetwheneachtriangleisnon-degenerateinthesenseof the above measure with respect to the quadratic form given by the local hessian d2f. We end by briefly recalling the greedy algorithm which was introduced in [8]. In 3,we study the behaviorofthe refinementprocedurewhenappliedto aquadraticfunctionq such § that its associated quadratic form q is of positive or negative sign. A key observation is that the edge which is bisected is the longest with respect to the metric induced by q. This allows us to prove that the triangles generated by the refinement procedure adopt an optimal aspect ratio in the sense of the non-degeneracy measure introduced in 2. § In 4, we study the behavior of the algorithm when applied to a general C2 function f which is § assumedto be strictly convex(or strictly concave). We firstestablisha perturbationresult, whichshows that whenf is locally close to a quadratic function q the algorithmbehavesin a similar manner as when applied to q. We then prove that the diameters of the triangles produced by the algorithm tend to zero so that the perturbation result can be applied. This allows us to show that the optimal convergence estimate limsupN f f C det(d2f) (1.5) N Lp Lτ k − k ≤ k | |k N →∞ p is met by the sequence of approximations f V generated by the algorithm. The extension of this result to an arbitNrar∈y CT2Nfunction f remains an open problem. It is possible to proceed to an analysis similar to 3 in the case where the quadratic form q is of mixed sign, also § provingthat the trianglesadopt anoptimal aspectratio as they get refined. We describe this analysis in 9.1 of [11]. However, it seems difficult to extend the perturbation analysis of 4 to this new setting. In § § particular the diameters of the triangles are no more ensured to tend to zero, and one can even exhibit examplesofnon-convexC2 functionsf forwhichtheapproximationf failstoconvergetowardsf dueto N thisphenomenon. Suchexamplesarediscussedin[8]whichalsoproposesamodificationofthealgorithm for which convergence is always ensured. However, we do not know if the optimal convergence estimate (1.5)holdsfor any f 2 with this modified algorithm,althoughthis seems plausible fromthe numerical ∈C experiments. 2 Adaptive finite element approximation 2.1 Notations We shall make use of a linear approximation operator that maps continuous functions defined on T T A onto Π . For an arbitrary but fixed 1 p , we define the local Lp approximationerror 1 ≤ ≤∞ e (f) := f f . T p T Lp(T) k −A k The critical assumptions in our analysis for the operator will be the following: T A 1. is continuous in the L norm. T ∞ A 2. AT commutes with affine changes of variables: AT(f)◦φ=Aφ−1(T)(f ◦φ) for all affine φ. 3. reproduces Π : (π)=π, for any π Π . T 1 T 1 A A ∈ Note that the commutation assumption implies that for any function f and any affine transformation φ:x x +Lx we have 0 7→ e (f) = det(L)1/pe (f φ) , (2.6) φ(T) p T p | | ◦ Two particularly simple admissible choices of approximation operators are the following: =P , the L2(T)-orthogonalprojection operator: (f P f)π =0 for all π Π . • AT T T − T ∈ 1 =I , the local interpolation operator: I f(v )=fR(v ) with v ,v ,v the vertices of T. T T T i i 0 1 2 • A { } 3 All our results are simultaneously valid when is either P or I , or any linear operator that fulfills T T T A the three above assumptions. Given a function f and a triangulation with N = #( ), we can associate a finite element N N T T approximationf defined on eachT by f (x)= f(x). The globalapproximationerroris given N N N T ∈T A by 1 kf −fNkLp = eT(f)pp p, (cid:16)TX∈TN (cid:17) with the usual modification when p= . ∞ Remark 2.1 The operator of best Lp(T) approximation which is defined by T B f f = min f π , T Lp(T) Lp(T) k −B k π∈Πmk − k does not fall in the above category of operators, since it is non-linear (and not easy to compute) when p=2. However, it is clear that any estimate on f f with f defined as f on each T implies N Lp N T 6 k − k A a similar estimate when f is defined as f on each T. N T B Here and throughout the paper, when q(x,y)=a x2+2a xy+a y2+a x+a y+a 2,0 1,1 0,2 1,0 0,1 0,0 we denote by q the associated quadratic form : if u=(x,y) q(u)=a x2+2a xy+a y2. 2,0 1,1 0,2 a a Note that q(u)= Qu,u where Q= 2,0 1,1 . We define h i a1,1 a0,2 (cid:18) (cid:19) det(q):=det(Q). If q is a positive or negative quadratic form, we define the q-metric v := q(v) (2.7) q | | | | which coincides with the euclidean norm when q(v)p= x2 +y2 for v = (x,y). If q is a quadratic form of mixed sign, we define the associated positive form q which corresponds to the symmetric matrix Q | | | | that has same eigenvectors as Q with eigenvalues (λ, µ) if (λ,µ) are the eigenvalues of Q. Note that | | | | generally q(u)= q(u) and that one always has q(u) q(u). | | 6 | | | |≤| | Remark 2.2 If detQ>0, then there exists a 2 2 matrix L and ε +1, 1 such that × ∈{ − } 1 0 LtQL=ε . 0 1 (cid:18) (cid:19) The linear change of coordinates φ(u) := Lu, where u = (x,y) R2, therefore satisfies q φ(u) = ∈ ◦ ε(x2+y2). On the other hand, if detQ<0 then there exists a 2 2 matrix L such that × 1 0 LtQL= . 0 1 (cid:18) − (cid:19) Defining again φ(u):=Lu we obtain in this case q φ(u)=x2 y2. ◦ − 4 2.2 From uniform to adaptive isotropic triangulations A standard estimate in finite element approximation states that if f W2,p(Ω) then ∈ inf f g Lp Ch2 d2f Lp, g∈Vhk − k ≤ k k where V is the piecewise linear finite element space associated with a triangulation of mesh size h h T h:=max diam(T). If we restrict our attention to uniform triangulations, we have T∈Th N :=#( ) h 2. h − T ∼ Therefore, denoting by eunif(f) the Lp approximation error by a uniform triangulation of cardinality N Lp N, we can re-express the above estimate as eunif(f) CN 1 d2f . (2.8) N Lp ≤ − k kLp This estimate canbe significantly improvedwhenusing adaptive partitions. We give here some heuristic arguments, which are based on the assumption that on each triangle T the relative variation of d2f is small so that it can be considered as constant over T (which means that f is replaced by a quadratic function on each T), and we also indicate the available results which are proved more rigorously. First consider isotropic triangulations, i.e. such that all triangles satisfy a uniform estimate h T ρ = A, (2.9) T r ≤ T where h := diam(T) denotes the size of the longest edge of T, and r is the radius of the largest disc T T contained in T. In such a case we start from the local approximationestimate on any T e (f) Ch2 d2f , T p ≤ Tk kLp(T) and notice that h2 d2f T d2f = d2f , Tk kLp(T) ∼| |k kLp(T) k kLτ(T) with 1 := 1+1and T theareaofT,wherewehaveusedtheisotropyassumption(2.9)intheequivalence τ p | | and the fact that d2f is constant over T in the equality. It follows that 1 1 e (f) C d2f , := +1. T p ≤ k kLτ(T) τ p Assumenowthatwecanconstructadaptiveisotropictriangulations withN :=#( )whichequidis- N N T T tributes the local error in the sense that for some prescribed ε>0 cε e (f) ε, (2.10) T p ≤ ≤ with c > 0 a fixed constant independent of T and N. Then defining f as (f) on each T , we N T N A ∈ T have on the one hand f f N1/pε, N Lp k − k ≤ and on the other hand, with 1 := 1 +1, τ p N(cε)τ f f τ Cτ d2f τ Cτ d2f τ . ≤ k − NkLp(T) ≤ k kLτ(T) ≤ k kLτ TX∈TN TX∈TN Combining both, one obtains for eiso(f) := f f the estimate N Lp k − NkLp eiso(f) CN 1 d2f . (2.11) N Lp ≤ − k kLτ This estimate improves upon (2.8) since the rate N 1 is now obtained with the weaker smoothness − condition d2f Lτ and since, even for smooth f, the quantity d2f Lτ might be significantly smaller ∈ k k 5 than d2f Lp. This type of result is classical in non-linear approximation and also occurs when we k k consider best N-term approximation in a wavelet basis. Theprincipleoferrorequidistributionsuggestsasimplegreedyalgorithmtobuildanadaptiveisotropic triangulation for a given f, similar to our algorithm but where the bisection of the triangle T that maximizes the local error e (f) is systematically done from its most recently created vertex in order to T p preservethe estimate (2.9). Such an algorithmcannotexactly equilibrate the errorin the sense of (2.10) and therefore does not lead to the same the optimal estimate as in (2.11). However, it was proved in [2] that it satisfies f fN Lp C f B2 N−1, k − k ≤ | | τ,τ for all τ such that 1 < 1 +1, provided that the local approximation operator is bounded in the τ p AT Lp norm. Here B2 denotes the usual Besov space which is a natural substitute for W2,τ when τ < 1. τ,τ Therefore this estimate is not far from (2.11). 2.3 Anisotropic triangulations: the optimal aspect ratio We now turn to anisotropic adaptive triangulations, and start by discussing the optimal shape of a triangle T for a given function f at a given point. For this purpose, we again replace f by a quadratic function assuming that d2f is constant over T. For such a q Π and its associated quadratic form q, 2 ∈ wefirstderiveanequivalentquantity forthe localapproximationerror. Hereandaswellasin 3and 4, § § weconsideratriangleT andwe denote by(a,b,c)its edgevectorsorientedinclockwiseoranti-clockwise direction so that a+b+c=0. Proposition 2.3 The local Lp-approximation error satisfies eT(q)p =eT(q)p T p1 max q(a), q(b), q(c) , ∼| | {| | | | | |} where the constant in the equivalence is independent of q, T and p. Proof: The first equality is trivial since q and q differ by an affine function. Let T be an equilateral eq triangle of area T = 1, and edges a,b,c. Let E be the 3-dimensional vector space of all quadratic eq | | forms. Then the following quantities are norms on E, and thus equivalent: e (q) max q(a), q(b), q(c) . (2.12) Teq p ∼ {| | | | | |} Note that the constants in this equivalence are independent of p since all Lp(T) norms are uniformly equivalent on E. If T is an arbitrary triangle, there exists an affine transform φ : x x +Lx such that T = φ(T ). 0 eq 7→ For any quadratic function q, we thus obtain from (2.6) eT(q)=eT(q)=eφ(Teq)(q)=|detL|p1eTeq(q◦φ)=|detL|p1eTeq(q◦L) since q L is the homogeneous part of q φ. By (2.12), we thus have ◦ ◦ eT(q) detL p1 max q(La), q(Lb), q(Lc) , ∼| | {| | | | | |} where a,b,c are againthe edge vectorsof T . Remarking that T = detL and that La,Lb,Lc are eq { } | | | | { } the edge vectors of T, this concludes the proof of this proposition. ⋄ In order to describe the optimal shape of a triangle T for the quadratic function q, we fix the area of T and try to minimize the error e (q) or equivalently max q(a), q(b), q(c) . The solution to T p | | {| | | | | |} this problem can be found by introducing for any q such that det(q) = 0 the following measure of 6 non-degeneracy for T: max q(a), q(b), q(c) ρq(T):= {| | | | | |}. (2.13) T det(q) | | | | p 6 Letφbealinearchangeofvariables,qaquadraticformandT atriangleofedgesa,b,c. Thendet(q φ)= ◦ (detφ)2det(q), the edges of φ(T) are φ(a),φ(b),φ(c) and φ(T) = detφ T . Hence we obtain | | | || | max q φ(a), q φ(b), q φ(c) max q(φ(a)), q(φ(b)), q(φ(c)) ρq φ(T)= {| ◦ | | ◦ | | ◦ |} = {| | | | | |} =ρq(φ(T)). ◦ T det(q φ) detφ T det(q) | | | ◦ | | || | | | (2.14) p p The last equation, combined with Remark 2.2, allows to reduce the study of ρ (T) to two elementary q cases by change of variable: 1. The case where det(q) > 0 is reduced to q(x,y) = x2 +y2. Recall that for any triangle T with edges a,b,c we define h :=diam(T)=max a, b, c , with the euclidean norm. In this case T {| | | | | |} |·| we therefore haveρq(T)= hT2T, which correspondsto a standardmeasure of shape regularityin the | | sense that its boundedness is equivalent to a property such as (2.9). This quantity is minimized when the triangle T is equilateral,with minimal value 4 (in fact it was alsoprovedin [5] that the √3 minimum of the interpolation error q I q among all triangles of area T =1 is attained T Lp(T) k − k | | when T is equilateral). For a general quadratic form q of positive sign, we obtain by change of variable that the minimal value 4 is obtained for triangles which are equilateral with respect to √3 the metric . More generally triangles with a good aspect ratio, i.e. a small value of ρ (T), are q q |·| those which are isotropic with respect to this metric. Of course, a similar conclusion holds for a quadratic form of negative sign. 2. Thecasewheredet(q)<0isreducedtoq(x,y)=x2 y2. Inthiscase,theanalysispresentedin[4] − shows that the quantity ρ (T) is minimized when T is a half of a square with sides parallel to the q x andy axes,withminimalvalue 2. Butusing(2.14)wealsonotice thatρ (T)=ρ (L(T))for any q q linear transformation L such that q = q L. This holds if L has eigenvalues (λ, 1), where λ = 0, ◦ λ 6 andeigenvectors(1,1)and( 1,1). Therefore,allimagesofthehalfsquarebysuchtransformations − L are also optimal triangles. Note that such triangles can be highly anisotropic. For a general quadratic form q of mixed sign, we notice that ρ (T) ρ (T), and therefore triangles which are q q ≤ | | equilateral with respect to the metric have a good aspect ratio, i.e. a small value of ρ (T). q q |·|| | In addition, by similar arguments, we find that all images of such triangles by linear transforms L with eigenvalues (λ, 1) and eigenvectors (u,v) such that q(u) = q(v) = 0 also have a good aspect λ ratio, since q=q L for such transforms. ◦ We leaveaside the special case where det(q)=0. In sucha case, the trianglesminimizing the errorfor a givenareadegenerateinthe sensethatthey shouldbe infinitely longandthin, alignedwiththe direction of the null eigenvalue of q. Summing up, we find that triangleswith a goodaspectratio arecharacterizedby the fact that ρ (T) q is small. In addition, from Proposition 2.3 and the definition of ρ (T), we have q 1 1 eT(q)p ∼|T|1+p1 |det(q)|ρq(T)=k |det(q)|kLτ(T)ρq(T), τ := p +1. (2.15) p p We nowreturnto afunctionf suchthatd2f is assumedtobe constantoneveryT . Assuming that N ∈T all triangles have a good aspect ratio in the sense that ρ (T) C q ≤ for some fixed constant C and with q the value of d2f over T, we find up to a change in C that e (f) C det(d2f) (2.16) T p Lτ(T) ≤ k | |k p By a similar reasoning as with isotropic triangulations, we now obtain that if the triangulation equidis- tributes the error in the sense of (2.10) f fN Lp CN−1 det(d2f) Lτ, (2.17) k − k ≤ k | |k p 7 andtherefore(1.1)holds. This estimateimprovesupon(2.11)sincethe quantity det(d2f) Lτ might k | |k be significantly smaller than d2f , in particular when f has some anisotropic features, such as sharp k kLτ p gradients along curved edges. The above derivation of (1.1) is heuristic and non-rigorous. Clearly, this estimate cannot be valid as suchsincedet(d2f)mayvanishwhiletheapproximationerrordoesnot(considerforinstancef depending only on a single variable). More rigorous versions were derived in [7] and [1]. In these results d2f is | | typically replaced by a majorant d2f +εI, avoiding that its determinant vanishes. The estimate (1.1) | | canthenberigorouslyprovedbutholdsforN N(ε,f)largeenough. Thislimitationisunavoidableand ≥ reflects the fact that enough resolution is needed so that the hessian can be viewed as locally constant over each optimized triangle. Another formulation, which is rigorously proved in [10], reads as follows. Proposition 2.4 There exists an absolute constant C > 0 such that for any polygonal domain Ω and any function f C2(Ω), one has ∈ limsupNe (f) C det(d2f) . N Lp Lτ ≤ k | |k N + → ∞ p 2.4 The greedy algorithm Givena targetfunction f, our algorithmiteratively builds triangulations with N =#( ) and finite N N T T element approximations f . The starting point is a coarse triangulation . Given , the algorithm N TN0 TN selects the triangle T which maximizes the local error e (f) among all triangles of , and bisects it T p N T from the mid-point of one of its edges towards the opposite vertex. This gives the new triangulation . N+1 T The critical part of the algorithm lies in the choice of the edge e a,b,c from which T is bisected. ∈{ } Denoting by T1 and T2 the two resulting triangles, we choose e as the minimizer of a decision function e e d (e,f), which role is to drive the generated triangles towards an optimal aspect ratio. While the most T natural choice for d (e,f) corresponds to the split that minimizes the error after bisection, namely T d (e,f)=e (f)p+e (f)p, T Te1 p Te2 p we shall instead focus our attention on a decision function which is defined as the L1 norm of the interpolation error d (e,f)= f I f + f I f . (2.18) T T1 L1(T1) T2 L1(T2) k − e k e k − e k e Forthisdecision,theanalysisofthealgorithmismadesimpler,duetothefactthatwecanderiveexplicit expressions of f I f when f = q is a quadratic polynomial with a positive homogeneous part T L1(T) k − k q. We prove in 3 that this choice leads to triangles with an optimal aspect ratio in the sense of a small § ρ (T). This leadsus in 4to aproofthatthe algorithmsatisfiesthe optimalconvergenceestimate (2.17) q § in the case where f is C2 and strictly convex. Remark 2.5 It should be well understood that while the decision function is based on the L1 norm, the selection of the triangle to be bisected is done by maximizing e (f) . The algorithm remains therefore T p governedbytheLp norminwhichwewishtominimizetheerror f f foragivennumberoftriangles. N p k − k Intuitively, this means that the Lp-norm influences the size of the triangles which have to equidistribute the error, but not their optimal shape. Remark 2.6 It was pointed out to us that the L1 norm of the interpolation error to a suitable convex function is also used to improve the mesh in the context of moving grid techniques, see [6]. We define a variant of the decision function as follows D (e,f):= f I f d (e,f). T T L1(T) T k − k − Note that D (e,f) is the reduction of the L1 interpolationerror resulting from the bisection of the edge T e,andthattheselectededgethatminimizesd (,f)isalsotheonethatmaximizesD(,f). Thefunction T · · D has a simple expression in the case where f is a convex function. T 8 Lemma 2.7 LetT beatriangleandletf beaconvexfunctiononT. LetebeanedgeofT withendpoints z and z . Then 0 1 T f(z )+f(z ) z +z 0 1 0 1 D (e,f)= | | f . (2.19) T 3 2 − 2 (cid:18) (cid:18) (cid:19)(cid:19) If in addition f has C2 smoothness, we also have T 1 D (e,f)= | | d2f(z )e,e min t,1 t dt, where z :=(1 t)z +tz . (2.20) T t t 0 1 6 h i { − } − Z0 Proof: Since f is convex, we have I f f on T, hence T ≥ f I f = (I f f). T L1(T) T k − k − ZT Similarly I f f on T1 and I f f on T2, hence Te1 ≥ e Te2 ≥ e D (e,f)= I f I f I f. T T T1 T2 − e − e ZT ZTe1 ZTe2 Let z be the vertex of T opposite the edge e. Since the function f is convex, it follows the previous 2 expression that D (e,f) is the volume of the tetrahedron of vertices T z +z z +z z +z 0,x 1,x 0,y 1,y 0 1 , ,f and (z ,z ,f(z )) for i=0,1,2. i,x i,y i 2 2 2 (cid:18) (cid:18) (cid:19)(cid:19) where (z ,z ) are the coordinates of z . Let u = z z and v = z z . We thus have D (e,f) = i,x i,y i 0 2 1 2 T − − 1 det(M) where 6| | u v ux+vx x x 2 M := u v uy+vy .  y y 2  f(z ) f(z ) f(z ) f(z ) f z0+z1 f(z ) 0 − 2 1 − 2 2 − 2   Subtracting the half of the first two columns to the third one w(cid:0)e find(cid:1)that M has the same determinant as u v 0 x x M˜ := uy vy 0 .   f(z ) f(z ) f(z ) f(z ) f z0+z1 f(z0)+f(z1) 0 − 2 1 − 2 2 − 2   Recalling that 2T = det(u,v) we therefore obtain (2.19). I(cid:0)n orde(cid:1)r to establish (2.20), we observe that | | | | we have in the distribution sense ∂2(min t,1 t +) = δ 2δ +δ , where δ is the one-dimensional t { − } 0− 1/2 1 t Dirac function at a point t. Hence for any univariate function h C2([0,1]), we have ∈ 1 h (t)min t,1 t dt=h(0) 2h(1/2)+h(1). ′′ { − } − Z0 Combining this result with (2.19) we obtain (2.20). ⋄ 3 Positive quadratic functions In this section, we study the algorithm when applied to a quadratic polynomial q such that det(q) >0. We shall assume without loss of generality that q is positive definite, since all our results extend in a trivial manner to the negative definite case. Our first observation is that the refinement procedure based on the decision function (2.18) always selects for bisection the longest edge in the sense of the q-metric defined by (2.7). q |·| 9 Lemma 3.1 An edge e of T maximizes D (e,q) among all edges of T if and only if it maximizes e T q | | among all edges of T. Proof: The hessian d2q is constant and for all e R2 one has ∈ d2qe,e =2q(e). h i If e is an edge of a triangle T, and if q is a convex quadratic function, equation (2.20) therefore gives T 1 T D (e,q)= | |q(e) min t,1 t dt= | | e2. (3.21) T 3 { − } 12| |q Z0 This concludes the proof. ⋄ It follows from this lemma that the longest edge of T in the sense of the q-metric is selected for bisection by the decision function. In the remainder of this section, we use this fact to prove that the refinement procedure produces triangles which tend to adopt an optimal aspect ratio in the sense that ρ (T) becomes small in an average sense. q For this purpose, it is convenient to introduce a close variant to ρ (T): if T is a triangle with edges q a,b,c, such that a b c , we define q q q | | ≥| | ≥| | q(b)+q(c) b2 + c2 σ (T):= = | |q | |q . (3.22) q 4T √detq 4T √detq | | | | Using the inequalities b2 + c2 2a2 and a2 2(b2 + c2), we obtain the equivalence | |q | |q ≤ | |q | |q ≤ | |q | |q ρ (T) ρ (T) q q σ (T) . (3.23) q 8 ≤ ≤ 2 Similar to ρ , this quantity is invariant under a linear coordinate changes φ, in the sense that q σ (T)=σ (φ(T)), q φ q ◦ From (2.15) and (3.23) we can relate σ to the local approximation error. q Proposition 3.2 There exists a constant C , which depends only on the choice of , such that for any 0 T A triangle T, quadratic function q and exponent 1 p , the local Lp-approximation error satisfies ≤ ≤∞ C0−1eT(q)p ≤σq(T)k detqkLτ(T) ≤C0eT(q)p. (3.24) where 1 := 1 +1. p τ p Our next result shows that σ (T) is always reduced by the refinement procedure. q Proposition 3.3 If T is a triangle with children T and T obtained by the refinement procedure for the 1 2 quadratic function q, then max σ (T ),σ (T ) σ (T). q 1 q 2 q { }≤ Proof: Assuming that a b c , we know that the edge a is cut andthat the childrenhave area q q q | | ≥| | ≥| | T /2 and edges a/2,b,(c b)/2 and a/2,(b c)/2,c (recall that a+b+c=0). We then have | | − − a b c 2T detq σ (T ) q +q − (3.25) q i | | ≤ 2 2 p (cid:16) (cid:17) (cid:18) (cid:19) b+c b c = q +q − (3.26) 2 2 (cid:18) (cid:19) (cid:18) (cid:19) q(b)+q(c) = (3.27) 2 = 2T detq σ (T). (3.28) q | | p ⋄ 10