The optimal aspect ratio for piecewise quadratic anisotropic finite element approximation

0.29 MB·English

by  Jean-Marie Mirebeau

THEOPTIMALASPECTRATIOFORPIECEWISEQUADRATIC
ANISOTROPICFINITEELEMENTAPPROXIMATION
Jean-MarieMirebeau
LaboratoireJacquesLouisLions,Universite´ PierreetMarieCurie,Paris,France
E-mail: [email protected]
ABSTRACT forsomeµ ∈IP andπ ∈IH . IfT isasufficientlysmall
z m−1 z m
1 triangle,wethushaveatleastheuristicallyonT
Mesh adaptation for finite element approximation is a proce-
1
0 dure used in numerous applications. The use of thin and long ∇(f −Im−1f)(cid:39)∇(π −Im−1π ), (3)
2 anisotropictrianglesimprovestheefficiencyoftheprocedure. T z T z
n When piecewise linear finite elements are used, the aspect sincetheLagrangeinterpolationoperatorImT−1 onthetriangle
a ratio for mesh adaptation is generally dictated by the absolute T reproducestheelementsofIPm−1.
J value of the (estimated) hessian matrix of the approximated ForanytriangleT andanyf ∈ H1(T)∩C0(T),wedefine

3 function. We give in this paper the corresponding aspect ratio theaveragedH1interpolationerrore (f) asfollows
T m
forpiecewisequadraticfiniteelements.
] 1 (cid:90)
A Keywords— Anisotropic finite elements, Adaptive meshes, eT(f)2m := |T| |∇(f −ImT−1f)|2.
N Interpolation,Nonlinearapproximation. T
. Thelocalcounterpartof(1)istheproblemoftheoptimaltrian-
h
t 1. INTRODUCTION gle: findforallπ ∈IHm
a
m Consider a bounded polygonal domain Ω ⊂ IR2, a sufficiently sup{|T|; T s.t. eT(π)m ≤1}. (4)
[ smooth function f : Ω → IR, and an integer m ≥ 2. We
Indeed the cardinality of a triangulation is inversely propor-
1 introducetheproblemofoptimalmeshadaptation
tional to the area of its elements. This approach is developed
v
0 min{#(T); T s.t. (cid:107)∇(f −Im−1f)(cid:107) ≤ε}, (1) in Chapter 2 of [1] and leads to asymptotically optimal error
1 T L2(Ω) estimates of (1) as ε → 0 (or more precisely estimates of ε
6 where T stands for an arbitrary triangulation of Ω, and #(T) as #(T) → ∞, which is equivalent). Unfortunately these
0 for its cardinality. Here Im−1 denotes the Lagrange interpola- estimates are not completely realistic for applications because
. T
1 tionoperatorontofiniteelementsofdegreem−1onT. currentlyavailablenumericalanisotropicmeshgeneratorsonly
0 In practical applications, the problem (1) is generally in- control the aspect ratio and orientation of the generated trian-
1 tractable for at least three reasons. 1: The function f may gles.
1
: have complicated local features, difficult to analyze. We thus ForeachtriangleT,ofverticesv1,v2 andv3,wedenoteby
v first make a local analysis based on Taylor developments. 2: zT := (v1 +v2 +v3)/3 its barycenter. We denote by S2+ the
Xi ThecollectionoftriangularmeshesofΩisacombinatorialset collectionof2×2symmetricpositivedefinitematrices,andwe
r andproblemssuchas(1)aretypicallyNP-complete(afterdis- defineamatrixHT ∈S2+bytheequality
a cretization).Weavoidthisproblembyfirstconsideringthecase
2 (cid:88)
of a single triangle. 3: Currently available anisotropic mesh HT−1 := 3 (vi−zT)(vi−zT)T.
generationalgorithmsonlygivecontrolontheaspectratioand 1≤i≤3
orientationofthegeneratedtriangles,butnotontheirotherfea-
IfAisaninvertible2×2matrixandifT(cid:48)ismappedontoT by
tures. Wethusonlyoptimizethisaspectratio.
thelinearmapz (cid:55)→Az,thenoneeasilychecksthat
2. ANOPTIMIZATIONPROBLEM H =ATH A. (5)
T(cid:48) T
We denote by IP the space of bivariate polynomials of de- By construction the triangle T of vertices
m−1 eq
gree≤m−1,andbyIH thespaceofhomogeneouspolynomi- (cos(2kπ/3),sin(2kπ/3)) satisfies H = Id.
m 0≤k≤2 Teq
alsofdegreem. Iff ∈Cm(Ω),ifz ∈Ωisfixedandifh∈IR2 Combining these two properties, Proposition 5.1.3 in [1]
issmall,thenlocally establishesthatforanytriangleT
(cid:112)
f(z+h)=µ (h)+π (h)+o(|h|m), (2) |T| detH =|T |,
z z T eq
wherediscπ := 4(ac−b2)(bd−c2)−(ad−bc)2 andλ :=
+
max{λ,0}.
Theorem. For m ∈ {2,3} the map π ∈ IH → M (π) is
m m
a near-minimizer of the problem (7) in the following sense. If
π is non-univariate then M (π) is non-degenerate. Further-
m
more there exists a constant C, independent of π, such that
e (π) ≤C and
Fig.1. Lagrange interpolationpointsfor IP andIP finiteele- Mm(π) m
1 2
ments(left),triangleT andassociatedellipseET (right). detM (π)≤Cinf{detM; M ∈S+s.t. e (π) ≤1}.
m 2 M m
andthatthereexistsarotationU (dependingonT)suchthat Proof: Theintegerm∈{2,3}isfixed,andwedenoteforeach
π ∈IH
m
1
z (cid:55)→UH2(z−z ) (6) (cid:107)∇π(cid:107):= sup |∇π(u)|.
T T
|u|≤1
mapsT ontoT (thepowerαofasymmetricpositivedefinite
eq Foreach2×2matrixAwedenotebyπ◦AtheelementofIH
m
matrixisobtainedbyelevatingtheeigenvaluestothepowerαin definedby(π◦A)(u):=π(A(u)), u∈IR2.Werecallthat
adiagonalization). Furthermoretheellipseofminimalvolume
containingT isET := {z; (z−zT)THT(z−zT)},seeFig1. ∇(π◦A)(u)=AT∇π(A(u)), u∈IR2,
The matrix H thus encodes the area, the aspect ratio and the
T
orientationofT. whichimpliesforanyrotationU
ForeachM ∈S+andeachπ ∈IH wedefine
2 m (cid:107)∇π(cid:107)=(cid:107)∇(π◦U)(cid:107). (8)
e (π) :=sup{e (π) ; T s.t. H =M}.
M m T m T
The main difficulty of this proof is to show that there exists a
Wefinallyintroduceforeachπ ∈ IH theproblemoftheopti- constantC ≥ 1suchthatforallπ ∈ IH andallM ∈ S+ one
m m 2
malaspectratioforIP interpolation has
m−1
inf{detM; M ∈S2+s.t. eM(π)m ≤1}. (7) C−1eM(π)m ≤(cid:107)M(cid:107)12(cid:107)∇(π◦M−21)(cid:107)≤CeM(π)m. (9)
Assume that this point is established. Proposition 6.5.4 in [1],
3. MAINRESULT
states that the map π (cid:55)→ M (π) is a near-minimizer for the
m
optimizationproblem
Ourmainresultisthesolutionoftheoptimizationproblem(7)
ienmtehnetsc.asTehoefppiieecceewwiisseeqluinaedarrataincdcapsieeciesweinsteirqeulyadnreawticafindnigteiveels- inf{detM; M ∈S2+s.t. (cid:107)M(cid:107)12(cid:107)∇(π◦M−12)(cid:107)},
a well founded answer to a long standing question: which as- in the same sense as in the statement of this theorem. Com-
pect ratio, depending on the third derivatives of the approxi- biningthisresultwiththeequivalence(9),andusingthehomo-
mated function, should be used in finite element software that geneity of π, we immediately conclude the proof of this theo-
combineanisotropyandIP2elements? rem.
Wefirstintroducesomenotation. Weequipthevectorspace Wethusturntotheproofof(9). Ourfirstobservationisthat
IHmwiththenorm thereexistsaconstantC0suchthatforallπ ∈IHm
(cid:107)π(cid:107):= sup |π(u)|. (cid:115) 1 (cid:90)
|u|≤1 e (π) := |∇(π−Im−1π)|2 ≤C (cid:107)∇π(cid:107),
Foreachπ ∈IH ,π =ax2+2bxy+cy2,wedefine Teq m |Teq| Teq Teq 0
2 (10)
(cid:18) a b (cid:19) indeedtheleftandrighthandsidearenormsonIHm.
[π]= b c . ConsiderasymmetricmatrixM ∈S2+andatriangleT such
thatH = M. Accordingto(6)thereexistsarotationU such
T
Theabsolutevalueofasymmetricmatrix(resp. thesquareroot that the image of T by the map z (cid:55)→ UM12(z − zT) is the
of a non negative symmetric matrix) is obtained by taking the triangleT .Injectingthischangeofvariablesin(10)weobtain
eq
absolutevalue(resp. squareroot)oftheeigenvaluesinadiago-
(cid:115)
nalization. Foreachπ ∈IH2weset |T1|(cid:90) |UM−12∇(π−ImT−1π)|2 ≤C0(cid:107)∇(π◦(M−21U−1))(cid:107).
M2(π):=(cid:107)π(cid:107)|[π]|. T
Foreachπ ∈IH ,π =ax3+3bx2y+3cxy2+dy3,weset Observing that (cid:107)Av(cid:107) ≥ (cid:107)A−1(cid:107)−1|v| for any invertible 2×2
3
matrixAandvectorv ∈IR2,andrecalling(8),weobtain
(cid:113) (cid:18)−discπ(cid:19)13
M (π):= [∂ π]2+[∂ π]2+ Id,
3 x y (cid:107)π(cid:107) + (cid:107)M(cid:107)−21eT(π)m ≤C0(cid:107)∇(π◦M−21)(cid:107).
Taking thesupremum of the lefthand side amongall triangles
T suchthatH =M weestablishtheleftpartof(9),provided
T
thatC ≥C .
0
WenowremarkthatthereexistsaconstantC suchthatfor
1
allπ ∈IH
m
(cid:115)
1 (cid:90)
(cid:107)∇π(cid:107)≤C |∂ (π−Im−1π)|2. (11)
1 |Teq| Teq x Teq
Indeed assume that the right hand side vanishes. Then µ :=
π−Im−1πisapolynomialofdegreemdependingonlyonthe
Teq
variable y, and which vanishes on the Lagrange interpolation
√
pointsofT ,seeFig1. Henceµvanishesfory = ± 3/2and
eq √ √
y = 0ifm = 2(resp. y = ± 3/2,y = ± 3/4andy = 0if
m=3). Thereforeµ=0whichimpliesthatπ =0. Bothsides
Fig. 2. Interpolation of (13) with IP elements on a uniform,
1
of(11)arethusequivalentnormsonthevectorspaceIH .
m isotropicoranisotropicmeshofcardinality500.
We consider a diagonalization of a symmetric matrix M12,
M ∈S2,
+ whereλ > 0isaconstant(thisexpressionneedstobeslightly
(cid:18) (cid:19)
M21 =UTDU, D = α0 β0 , imnoodridfieerdtiof πenzsuvraenitshhaetsHor i∈s uCni0v(aΩri,aSte+f)o).r sSoommeevamlueesshogfenz-,
2
erators such as [2] can, at least heuristically, and provided H
whereU isarotationandα = (cid:107)M(cid:107)12. ConsiderthetriangleT has sufficient regularity, produce a mesh T of Ω such that
whichismappedontoT bythechangeofcoordinates C−1H(z) ≤ H ≤ CH(z) for each T ∈ T each z ∈ T,
eq T
where C is a constant not too large. In other words the aspect
z (cid:55)→UM12z =DUz, ratio of the elements of T is dictated by the metric H. Some
rigorous results in this direction can be found in Chapter 5 of
and thus satisfies HT = M according to (5). Injecting this [1].
changeofvariablesinto(11)weobtain In the expression (12) the matrix M (π ) ensures that the
m z
elementsofT havetheoptimalaspectratio,whilethescalarfac-
(cid:115)
(cid:107)∇(π◦(M−21U−1))(cid:107)≤C1 |T1|(cid:90) |α−1v·∇(π−ImT−1π)|2, teoqru(iddiesttrMibumte(dπza)m)−on2g1mthgeuaerlaenmteeenststhoaftTthe(aingteenrperoallatpiroinnceirprloeriins
T
adaptiveapproximation).
wherev :=U−1e ,e :=(1,0),andwhereweusedforthe∂ Weconductedsomenumericalexperimentsusing[2]andfor
x x x
derivativethatU−1D−1ex =α−1v.Recallingthatα=(cid:107)M(cid:107)12, thesyntheticfunction
|v|=1,andusing(8)weobtain f(x,y):=tanh(10(sin(5y)−2x))+x2y+y3 (13)
(cid:115) 1 (cid:90) onthedomainΩ := (−1,1)2. Theyillustratetheimprovement
(cid:107)M(cid:107)21(cid:107)∇(π◦M−12)(cid:107) ≤ C1 |T| |∇(π−ImT−1π)|2 offeredbyanisotropicmeshadaptation, bothinthecaseofIP1
T andIP elements,foratriangulationofcardinality500.
2
= C e (π) ≤C e (π) .
1 T m 1 M m
#(T)=500 Uniform Isotropic anisotropic
Thisconcludestheproofof(9)withC :=max{C0,C1},hence (cid:107)∇(f −I1T f)(cid:107)L2 110 51 11
theproofofthistheorem. (cid:5) (cid:107)∇(f −I2 f)(cid:107) 79 14 0.88
T L2
Ournextobjectiveistocombineouranalysiswithanadaptive
4. APPLICATIONSANDCONCLUSION anisotropicmeshrefinementprocedure,forapartialdifferential
equationsolvedwithIP finiteelements.Theoptimizationprob-
2
Consider a function f for which one desires to solve, at least
lem(1)isparticularlyrelevantinthecaseofellipticequations.
heuristically, the optimization problem (1). Assume that some
estimateof 5. REFERENCES
(cid:88) ∂mf xk yl
πz := ∂kx∂ly(z)k! l! [1] Jean-Marie Mirebeau, Ph.D Thesis, Adaptive and
k+l=m anisotropic finite element approximation : Theory
isknownateachpointz ∈ Ω, anddefineariemannianmetric and algorithms, tel.archives-ouvertes.fr/
H onΩasfollows tel-00544243/en/
[2] FreeFem++ software, developped by Frederic Hecht,
H(z):=λ(detMm(πz))−21mMm(πz), (12) www.freefem.org/ff++/