Manuscriptsubmittedto Website: http://AIMsciences.org AIMS’Journals VolumeX,Number0X,XX200X pp. X–XX LEVELSETS AND ANISOTROPIC MESH ADAPTATION Alexandra Claisse and Vincent Ducrot UPMCUnivParis06,UMR7598,LaboratoireJ.L.Lions,F-75005Paris,France Pascal Frey UPMCUnivParis06,UMR7598,LaboratoireJ.L.Lions,F-75005Paris,France UniversidaddeChile,UMI2807,CentrodeModelamientoMatem´atica,Santiago,Chile (Communicated by the associate editor name) Dedicated to Prof. Li Ta-tsien on the occasion of his 70th birthday Abstract. In this paper, we focus on the problem of adapting dynamic tri- angulationsduringnumericalsimulationstoreducetheapproximationerrors. Dynamically evolving interfaces arise in many applications, such as free sur- faces in multiphase flows and moving surfaces in fluid-structure interactions. In such simulations, it is often required to preserve a high quality interface discretization thus posing significant challenges in adapting the triangulation inthevicinityoftheinterface,especiallyifitsgeometryoritstopologychanges dramaticallyduringthesimulation. Ourapproachcombinesanefficientlevelset formulationtorepresenttheinterfaceintheflowequationswithananisotropic mesh adaptation scheme based on a Riemannian metric tensor to prescribe size,shapeandorientationoftheelements. Experimentalresultsareprovided toemphasizetheeffectivenessofthistechniquefordynamicallyevolvinginter- facesinflowsimulations. 1. Introduction. Nowadays, many computational simulations of PDE problems involve adaptive triangulations. In such triangulations, the mesh points are locally adapted in order to achieve high concentration of nodes in regions where the gradi- ent of the solution is important. In general, adaptive methods attempt to equidis- tribute the approximation errors by adapting the mesh density to a discrete metric tensor field based on the Hessian of the solution [22]. In recent years, anisotropy allowed to significantly improve the accuracy of the solution and of the surface dis- cretization as well as to maintain the order of convergence of the computational schemes[5]. These techniques have successfully proved their ability in reducing the number of unknowns while preserving the accuracy of the numerical solutions, at least in steady-state simulations. However, when dealing with unsteady phenom- ena involving dynamically changing interfaces, such as free surfaces in multiphase flows or moving boundaries in fluid-structure interactions, the challenge becomes to keep updated with the arbitrary evolution of the interfaces, including changes of topology. Usually, an interface can be characterized by a discontinuity of physicial quan- tities like pressure, viscosity or density. Motion can be described in Eulerian or 2000 Mathematics Subject Classification. Primary: 65N60;Secondary: 28A75. Key words and phrases. levesets,meshadaptation,anisotropy,metrictensor. ThefirstauthorissupportedbyagrantfromtheR´egionIledeFrance. 1 2 A. CLAISSE, V. DUCROT AND P. FREY Lagrangiancoordinatesincontinuummechanics. InLagrangiannumericalschemes, each mesh vertex is associated with a particle of material and shall follow the dis- placement of this particle during time. Hence, flows with significant deformations mayleadtolargemeshdistortionsandsubsequentaccuracyproblemsrelatedtonu- mericalerrorsthatareverytedioustoremoveusingcomplexremeshingprocedures. In this spirit, vorticity is often a cause of mesh tangling, when i.e., overlapping, inverted or non-convex elements are created, thus leading to aborting the ongo- ing computation. Eulerian schemes, on the other hand, the computational mesh remains fixed as the medium moves with respect to the mesh. Hence, these meth- ods intrinsically avoid any mesh-related difficulty although the interfaces and the boundaries must be accurately discretized to resolve flow details. Moreover, Euler- ian methods may lead to poor performances in solid mechanics simulations. It is thus important to rely on robust and efficient meshing techniques to discretize the computational domain. In this paper, we investigate the problem of adapting dynamic triangulations during numerical flow simulations to reduce approximation errors. More precisely, given an interface of codimension one between two fluids, we are interested in com- puting its subsequent motion under a velocity field v related to the physics of the problem. Following the ideas suggested by Dervieux-Thomasset [17] and Osher- Sethian [33], we introduce a levelset function that represents the interface in the incompressibleflowequationsandthemotionistrackedbyconvectingthelevelsets with the velocity field v using a first order Hamilton-Jacobi equation. We propose a dynamic anisotropic adaptation of the triangulation using a Riemannian metric tensor to prescribe the size and the stretching of the elements in the vicinity of the interfaceincombinationwithmeshmodificationoperations. Thelatterareusedfor maintainingthemeshqualityaswellastheaccuracyofthecalculation. Weprovide a theoretical bound on the geometric approximation error when discretizing the interface with anisotropic elements. As we have complex flow simulations in mind, it is very important to have an accurate representation of the computational domain boundaries (curves or sur- faces). Very few approaches have been proposed for anisotropic (re)meshing of curve or surface triangulations [23] and dynamic surfaces introduce significant dif- ficulties and constraints to mesh adaptation. A desirable procedure should be able to deal with complex geometry and topology as well as with non uniform data. We will show that the control of the geometric representation of an interface using an anisotropic metric tensor field can be used to the purpose of reconstructing a curve orasurfacefromanunorganizeddataset, i.e., apointcloud. FollowingZhaoet al. [41],weproposeavariationalapproachbasedonadistancefunctiontothedataset. The regularization method is related to the local density of points and provide a smoother surface reconstruction than piecewise affine. The differential properties of the distance function will allow to define an anisotropic metric tensor used to enhance the boundary representation of the domain. The remainder of this paper is structured as follows. Section 2 presents some background material on anisotropic mesh generation and metric-based mesh adap- tation. Section3introducesalevelsetformulationfortime-dependentandinterface tracking problems. Section 4 presents numerical results and discussions for static and dynamic interface capturing. Finally, Section 5 concludes the paper with an introduction on open problems in dynamic simulations. LEVELSETS AND ANISOTROPIC ADAPTATION 3 2. Metric-based mesh adaptation. In this section, we present some basic no- tions and definitions on triangulations and meshes, metric tensor fields and aniso- tropic mesh adaptation. Throughout this paper, Ω denotes a simply connected, open, bounded domain in Rn; Ω is the closure of Ω and |Ω| represents its n-dimensional measure or its volume. 2.1. Basicnotions. LetconsiderapolyhedraldomainΩinRnandsupposegivena familyoftriangulationsT onΩ,wheretheparameterh>0representstheelement h size characterizing this family. We assume that each element K in T is a closed h (cid:83) sub-domain of Ω and that Ω⊂ K and also that the intersection of any two K∈Th different elements is either empty or reduced to a n−1-simplex (elements cannot overlap). Under these assumptions, let define a uniform mesh T of Ω as a mesh in h which all elements are equally sized and regular (equilateral). In such case, (cid:18) (cid:19)1 |Ω| n h ≈ ∀K ∈T , (1) K |T | h h where |T | denotes the number of mesh elements in T . Following [30], a regu- h h lar family of meshes is a family which satisfies the conditions: (i) there exists a constant τ such that: hK ≤τ, ∀K ∈(cid:91)T , ρ h K h where h = diam(K) and ρ = sup{diam(B),B is a ball contained in K} K K are the diameter and in-diameter of K, respectively, (ii) the size is such that: h= max h →0. (2) K K∈Th A quasi-uniform mesh is a mesh for which (i) is satisfied and the variation of its element size is bounded by a constant. Hence, Relation (1) is true for a quasi- uniform meshes. In the context of finite element approximation spaces, it may be usefull to introduce the following notion. An affine family of triangulations is a family such that any element K in the family can be mapped from a reference element Kˆ through an affine mapping. This is equivalent to state that there exists an invertible affine mapping F : Kˆ (cid:55)→ K, such that K = F(Kˆ), and in addition, the reference element is chosen to be equilateral with unitary size |Kˆ|=1. 2.2. Mesh adaptation. Mesh adaptation aims at improving the efficiency and theaccuracyofnumericalsolutionsbyconcentratingmorenodesinregionsoflarge solution variations than in other regions of the computational domain. As a con- sequence, the number of mesh nodes required to achieve a given accuracy can be dramatically reduced thus resulting in a reduction of the computational cost. Tra- ditionally, isotropic mesh adaptation has received much attention, where almost regular (equilateral) mesh elements are only ajusted in size based on an error es- timate. However, in regions of large solution gradient, adaptive isotropic meshes usuallycontaintoomanyelements. Moreover,iftheproblemathandexhibitsstrong anisotropicfeatures, likeshockwavesorboundarylayersinfluiddynamics, itisde- sirable to adjust the element size as well as the element shape and orientation to betterrepresentforthesolutionvariations. Inmostcases,elementswithhighaspect ratio are needed, thus resulting in an anisotropic mesh. In recent years, research 4 A. CLAISSE, V. DUCROT AND P. FREY on anisotropic mesh adaptation has been intensified and consequently mathemat- ically sound error estimates have been developed, see for instance [8], [11], [21], [19]. These works have led to the concept of optimal triangles [10] and the affine interpolation error and its high-order derivatives, gradient and Hessian, for regular (quadratic) functions provided a convenient way of dealing with anisotropy using the so-called metric tensor field [3] [25]. Fromapracticalpointofview,classicalmeshadaptationalgorithmsfallintotwo categories,h-methodsorr-methods. Ontheonehand,h-methodsinvolvelocalmesh modification operations such as edge flipping, edge-vertex collapsing, vertex addi- tion and vertex relocation. Anisotropic meshes are generated using various strate- gies, like the Delaunay-based triangulation method, the advancing-front method, the quadtree-octree-based method and methods combining local mesh modifica- tions; see [22] for an overview. A metric tensor prescribes the size, stretching and orientation of mesh elements at any place in the domain and is usually associated with mesh vertices. On the other hand, r-methods employ only vertex relocation procedures to adapt the mesh. Typically, they are based on solving elliptic or par- abolic PDE systems [14] or use a functional to account for desired mesh properties like orthogonality or smoothness [39]. For unsteady problems, the nodes are moved according to the solution of an equation involving velocities of nodes in order to preserve the concentration of nodes in regions of high gradient of the numerical solution. Moving mesh methods can be based on a variational mesh generator [29]. 2.3. Metric tensor. Asmentioned,anisotropicmeshadaptationreliesonanerror estimate or an error indicator to prescribe the element size, shape and orientation using a matrix-valued field. In the continuous geometric viewpoint, the metric defininganelementisrepresentedbyanellipsoid. Hence,size,shapeandorientation notionsareassociatedwithitsvolume,lengthratiosbetweenthelengthsofitssemi- axes and its principal axis vectors, respectively. Commonly in anisotropic mesh generation, the metric tensor denoted as M(x) allowstocreateaquasi-uniformmeshinthemetricrelatedtoM. Moreprecisely,the meshelementsareequilateralinthemetric,thiscorrespondingtothecondition[30]: tr(cid:0)(F(cid:48))tM F(cid:48)(cid:1)=ndet(cid:0)(F(cid:48))tM F(cid:48)(cid:1)n1 , ∀K ∈T , (3) K K h where M is an average of M(x) on element K and F(cid:48) is the Jacobian matrix of K the affine mapping F. Moreover, the volume of the element in T is unitary: h (cid:90) (cid:112) det(M(x))dx=1, ∀K ∈T , (4) h K or, in a discrete formulation: (cid:112) |K| det(M )=1, ∀K ∈T . (5) K h By extension, in the metric given by M(x) for any x in Ω, the length of a curve γ is defined as follows: (cid:90) 1(cid:113) (cid:90) 1(cid:112) |γ| = (γ(t),γ(t)) dt= (γ(t),M(γ(t))γ(t))dt (6) M M(γ(t)) 0 0 where γ(t):[0,1]→Rn is a parametrization of γ. In Rn, the metric M associated with an element K can be represented by its K unit sphere, an ellipsoid B(K) defined as [25]: (cid:110) (cid:112) (cid:111) B(K)= P ∈Rn,(cid:107)OP(cid:107) = (OP,M OP)=1 MK K LEVELSETS AND ANISOTROPIC ADAPTATION 5 where O denotes the center of B(K). We assume that the metric tensor M(x) is a symmetricpositivedefiniten×nmatrix. Then,thespectraldecompositiontheorem leads to express matrix M in terms of its n eigenvalue-eigenvector pairs (λ ,e ) as: i i n (cid:88) M =P ΛPt = λ e et i i i i=1 where the normalized eigenvectors of M are the columns of matrix P =[e ,...,e ] 1 n such that P Pt =PtP =I and Λ is the diagonal matrix with eigenvalues λ on the i diagonal. Thegeometricmeaningofthisdecompositioniseasytounderstand. Asa matteroffact,thematrixP determinestheorientationandthematrixΛprescribes the size and shape of any element K in T . h To describe the mesh modifications used to create anisotropic meshes, we need to introduce primarily the notion of mesh quality. 2.4. Mesh quality. Itisimportanttointroducequalitymeasuresinordertoeval- uatehowcloselythecontrolsonsize,shapeandorientationprescribedbythemetric M(x) at any x in Ω are satisfied. Mesh quality assessment has been an area of ex- tensive research in the context of finite element methods, essentially for isotropic meshes (see [9] for a review). Mesh quality measures can be defined to assess individual elements as well as the entire mesh T . The geometric (shape) quality measure of an element K can h be defined based on the Jacobian matrix as:  tr(cid:0)(F(cid:48) )tF(cid:48) (cid:1) 2(nn−1) (cid:32) (cid:107)F(cid:48) (cid:107) (cid:33)n−n1 Qgeo(K)=ndet(cid:0)(FKK(cid:48) )tFKK(cid:48) (cid:1)n1  = √ndeKt(FFK(cid:48) )n1 (7) where (cid:107)·(cid:107) denotes the Frobenius matrix norm. Notice that Q (K) ≥ 1 and F geo that the best quality Q (K) = 1 is achieved for an equilateral element in T . geo h Following [30], it is easy to prove that for all K ∈T : h  (cid:32)(cid:18) (cid:19) 1 (cid:33)2 2(nn−1) n(n1−1) µµmax 2n −1 +1 ≤ Qgeo(K) ≤ µµmax , (8) min min where µ and µ are the minimal and maximal singular values of F(cid:48) , respec- min max K tively. These inequalities mean that the geometric quality of K is equivalent to the aspect ratio of K, as the dimension n is small. By extension, the geometric mesh quality is defined as: Q (T )= max Q (K) (9) geo h geo K∈Th Similarqualitymeasurescanbedefinedtoassessthealignmentandequidistribution conditions. However, for practical reasons, we advocate the use of a single measure to assess the quality of an element, like for instance: (cid:18) k (cid:19)n (cid:80) (e ,M e ) i K i Q (K)=α i=1 (10) ani (cid:112) |K| det(M ) K where e denotes here any of the k edges of K and α is a normalizaton coefficient i such that Q (K) = 1 for an equilateral element. Again, Q (K) ≥ 1 for all K ani ani in T . The larger max Q (K) is, the farther the mesh deviates from satisfying h K ani the shape and size and orientation conditions. 6 A. CLAISSE, V. DUCROT AND P. FREY 2.5. Error estimates. Among the various meshing strategies developed over the years,thedefinitionofthemetrictensorbasedontheHessianofasolutionvariable seems nowadays commonly generalized in the meshing community. This choice is largelymotivatedbytheinterestingnumericalresultsobtainedbytheseminalwork of D’Azevedo [10] on linear interpolation for quadratic functions on triangles. In this approach, the metric tensor is defined as: M =|H(v)|=P |Λ|Pt, |Λ|=diag{|λ |} , (11) i i=1,n giventhedecompositionoftheHessianoffunctionv. ThemetricM isfurthermod- ified by imposing a minimal and maximal edge lengths to avoid unrealistic metric. The local sizes of an element with respect to the main directions prescribed by the √ eigenvectors are related to the eigenvalues: h =1/ λ , for i=1,...,n. Moreover, i i the metric indicates the local density (distribution) of mesh nodes at any position x in the domain directly from the local sizes h . Let suppose the mesh density √ i function d be defined as d = (cid:81)N λ , where N denotes the current number of i=1 i vertices,thentheoptimalnumberofmeshverticesorthemeshcomplexity,hereafter denoted as C, will be given by the formula, for all x in Ω: (cid:90) (cid:90) (cid:89)N 1 C = d(x)dx= dx. (12) h (x) Ω Ωi=1 i Given a metric tensor M defined on Ω and considering a solution variable u, we proposed the following error model for the local interpolation error on a mesh T , h for all x in Ω [4], : e (x)= max (cid:107)u−P u(x)(cid:107) , (13) M h ∞,K K∈B(x) (cid:112) where B(x) = {y ∈ Ω, (y,M(x)x) ≤ 1} denotes the unit ball centered at x. A practical bound on the interpolation error on an element K ∈ T in L∞ norm is h given by (see [25] for more details on the estimation): e (K)=(cid:107)u−P u(x)(cid:107) ≤ c max max(v,|H(u(y))|v), (14) M h ∞,K n y∈K v⊂K where v ⊂ K represents any vector inscribed in K and c is a constant related to n the dimension of the space Rn. It is possible to show that the optimal directions of the metric coincide with the main directions of the Hessian, thus leading to the following error model: eM(x)=(cid:88)N h2i (cid:12)(cid:12)(cid:12)(cid:12)∂∂α2u2(cid:12)(cid:12)(cid:12)(cid:12) , (15) i=1 i under the assuption that the Hessian of a C2 continuous function u can be decom- posed as follows: (cid:18)∂2u(cid:19) H(u)=R(u)diag R(u)−1 (16) ∂α2 i where the eigenvectors of H(u) are the columns of the matrix R(u) = [e ,...,e ]. 1 n ThenextstepistosearchforafunctionM(x)minimizingtheLp normoftheerror, for a given number N of vertices. This means solving the problem: mMin(cid:90)Ω(eM(x))pdx=mhiin(cid:90)Ω(cid:32)(cid:88)i=N1h2i(x)(cid:12)(cid:12)(cid:12)(cid:12)∂∂α2ui2(x)(cid:12)(cid:12)(cid:12)(cid:12)(cid:33)p dx, (17) LEVELSETS AND ANISOTROPIC ADAPTATION 7 under the constraint: (cid:90) N (cid:90) (cid:89) N = h−1(x)dx = d(x)dx. (18) i Ωi=1 Ω This optimization problem is solved in two steps; at first the anisotropic quotients (cid:18) hn (cid:19)1/n r = i i (cid:81)n h k=1 k are computed an then the density d(x) is evaluated. As a result, the optimal continuous metric is defined using the following decomposition: M =D R(u)−1diag((cid:107)λ (cid:107))R(u), (19) Lp Lp i where the matrix D and the eigenvalues λ are given in the following table: Lp i Norm D (cid:107)λ (cid:107) Lp i Lp N23 (cid:90)Ω(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)i(cid:89)=n1∂∂α2ui2(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)2pp+3−23 (cid:12)(cid:12)(cid:12)(cid:12)∂∂α2ui2(cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)i(cid:89)=n1∂∂α2ui2(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)−2p1+3 L1 N23 (cid:90)Ω(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)i(cid:89)=n1∂∂α2ui2(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)15−32 (cid:12)(cid:12)(cid:12)(cid:12)∂∂α2ui2(cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)i(cid:89)=n1∂∂α2ui2(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)−51 L2 N23 (cid:90)Ω(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)i(cid:89)=n1∂∂α2ui2(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)27−32 (cid:12)(cid:12)(cid:12)(cid:12)∂∂α2ui2(cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)i(cid:89)=n1∂∂α2ui2(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)−71 L∞ N23 (cid:90)Ω(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)i(cid:89)=n1∂∂α2ui2(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)12−32 (cid:12)(cid:12)(cid:12)(cid:12)∂∂α2ui2(cid:12)(cid:12)(cid:12)(cid:12) Definition 2.1. AssumingasequenceofmetricsM correspondingtoagivennum- i ber of vertices N, the order of convergence k for the norm (cid:107)e (x)(cid:107) is attained Mi Lp if the following inequality holds: (cid:107)e (x)(cid:107) ≤cstN−k/n, (20) Mi Lp According to this definition, all the metrics previously defined lead to a conver- gence order of two (k =2). Inanumericaladaptivecomputation,theunknownHessianofasolutionvariable H(u(x))atanyxinΩ,mustbereplacedbyadiscreteapproximationH (x). In[2], h sufficientconditionsaregivenfortheapproximationH underwhichthemeshquasi- h uniform in the metric |H | is also quasi-uniform in the metric |H(x)|. Obviously, h if the initial mesh used to recover the Hessian H is far from the optimal one, the h discrete hessian may not approximate the continuous one. A few methods allow to recover the hessian H from the piecewise linear function u defined at mesh h h vertices. For instance, a weak definition of the discrete Hessian is established using Green formula as follows: (cid:90) (cid:90) ∂u ∂ϕ Hijϕ dx=− h x dx, (21) h x ∂x ∂x B(x) B(x) i j where B(x) denotes the stencil of vertex x and ϕ represents the piecewise linear x P finite element function associated with the interior node x. At a boundary 1 node, the Hessian H (x) can be extrapolated from its neighboring interior nodes. h 8 A. CLAISSE, V. DUCROT AND P. FREY OtherHessianrecoverymethodsinvolveaprojectiontechnique[42]oraleast-square approximation of the Hessian matrix based on a second order Taylor expansion of a C2 continuous function u∗ sufficiently close to the numerical solution u [25]. All h h these techniques have linear complexity with respect to the number of mesh nodes. It turns out that the choice of the recovery method depends strongly on the norm in which the error has to be minimized. 2.6. Anisotropic mesh adaptation. For the generation of anisotropic meshes, h-adaptation has been the most widely used approach, mainly because of it is ro- bust and reliable as well as conceptually simple (see [22] for a survey of h and p-methods). Typically, a local error estimate is designed and then local mesh mod- ification operations such as vertex deletion, vertex addition, edge swapping and vertex relocation, are iteratively used to produce the desired adapted mesh. The main idea consists in generating quasi-uniform meshes in the metric supplied by the error estimate that prescribes the size, shape and orientation of the elements throughout the computational domain. Severalmeshingstrategieshavebeendevelopedforgeneratinganisotropicmeshes according to a metric prescription. For example, such techniques include the De- launay triangulation method [27], [24], [18], the advancing-front method [26], the bubblemeshmethod[40]andmethodscombininglocalmodificationswithsmooth- ing [13]. On the other hand, variational methods have received much attention in the recent years typically for generating structured meshes as they are especially well suited for finite difference schemes. In these approaches, an adaptive mesh is considered as the image of a computational mesh under a coordinate transforma- tiondeterminedbyafunctional,forinstancebasedongeometricconsiderations[14]. The key issue is then to design the proper convex functionals, usually not directly provided by an error estimate. In this section, attention is focused on the anisotropic Delaunay-based method developed in [18]. Since the early paper of Delaunay [16], the properties of Delau- nay triangulations have been extensively analyzed, especially with respect to their computational geometry features. Here, we consider the extension of the Delaunay measure to generate anisotropic simplices. At an elementary level, constructing a Delaunay triangulation requires finding all simplices having a circumsphere con- taining a given point p in Rn. This operation can be performed using the following test, the Delaunay measure, for all K in T : h d(O,p) α(p,K)= <1, (22) ρ(K) where ρ(K) is the circumradius of K and O is the center of the circumsphere. The anisotropic measure of α is derived from the previous formula by computing all distance terms with respect to a prescribed metric tensor field M [25]. Considering a matrix-valued field leads however to solving a system of non linear equations and is then more computationally expensive [18]. Since mesh quality is related to the quality of the worst element in the mesh, it isimportanttoattemptremovingbadlyshapedorincorrectlysizedelementsinthe mesh[22]. Tothisend,localmeshmodificationsarecarriedoutinaniterativeman- ner, basedonanedgelengthanalysiswithrespecttothemetricM. Hence, smaller elements are removed using an vertex removal or edge contraction procedure while larger elements are further refined using an edge splitting technique. Additionally, edgeflippingandnoderelocationoperationsareusedtoimprovetheoverallquality LEVELSETS AND ANISOTROPIC ADAPTATION 9 ofthemesh. Anefficiency index allowstomeasuretheadequacyofthemeshedges with the supplied metric tensor field: (cid:32) N (cid:33) 1 (cid:88) τ(T )=exp min(l ,l−1)−1 (23) h N i i i=1 where l denotes the length of edge i in the current mesh. i Regarding the adaptation of curves and surfaces, specific procedures must be developed to account for the differential geometric properties of the manifolds [25]. To generalize the above formulation to surface meshes, we prefer to define a three- dimensionalmetrictensoronalocalparametricspace,forinstanceusingaquadric- based surface analysis [23]. In this approach, the curvature tensors at each ver- tex of a given triangulation are estimated to approximate the Hessian of a locally parametrized surface. Aswehaveseen,metric-basedadaptivemeshingtechniquesprovideaconvenient yet powerful framework to numerical simulations. The numerical accuracy can be largely improved by adapting the mesh node density to the metric tensor field provided by an error estimate ased on the Hessian of the solution. In the next section,wewillprovidetwoexamplesofanisotropicadaptationcoupledwithlevelset techniques for dynamic simulations. 3. Levelset methods. Introduced by Dervieux [17] more than three decades ago and rediscovered a few years later by Osher [33], levelsets have since proven to be successful for dealing with moving fronts and interfaces. In this section, we will reviewthebasicconceptsbehindthelevelsetmethodsandproposetwoapplications of levelsets to emphasize their effectiveness in tracking evolving interfaces using anisotropic triangulations. 3.1. Basic concepts. The underlying concept behind the levelset method is quite simple. Given an interface in Rn of codimension one bounding an open domain Ω, the objective is to compute its motion under a velocity field v(t,x). To this end, a Lipschitz continuous function φ(t,x), x ∈ Rn, is introduced and the interface is chosen to coincide with the set where φ(t,x) = 0. The function φ is assumed to have the following properties:  φ(t,x) < 0, for x∈Ω(t)  φ(t,x) > 0, for x∈/ Ω¯(t) (24)  φ(t,x) = 0, for x∈∂Ω(t)=Γ(t). Hence, at all time the interface coincide with the set Γ(t) for which φ vanishes. In other words, if the velocity v depends on the surface geometry and if the domain Ω(t) evolves in time at a speed v(t,x), for instance v = (α−βH)n, with n the normal to the surface and H its mean curvature, then φ(t,x(t)) = 0 for all x in ∂Ω(t). Using the derivation rule will lead to the subsequent equation: ∂φ ∂φ +x(cid:48)(t).∇φ= +vn.∇φ=0, ∂t ∂t and as one can write: ∇φ |∇φ|2 n.∇φ= .∇φ= =|∇φ|, (25) |∇φ| |∇φ| 10 A. CLAISSE, V. DUCROT AND P. FREY the following first order Hamilton-Jacobi equation is then posed in Rn: ∂φ +v|∇φ|=0. (26) ∂t As n and ∇φ are colinear vectors then T .∇φ=0 for all tangent vector T. In two dimensions, v = v n+v T, the equation φ +(v n+v T).∇φ = 0 shows that n t t n t only the normal component of the velocity is involved and so it is easy to deduce the characteristic equation of levelsets: ∂φ (t,x)+v (t,x)|∇φ(t,x)|=0. (27) ∂t n By embedding the interface as a levelset of a higher dimensional function, the dimension of the advection problem is augmented by one. But, even if they are computationally more expensive to discretize than other representation methods, levelsets offer a convenient way to handle topological changes. Unfortunately, singularities are typical of the solutions of the Hamilton-Jacobi equation (27). From the mathematical point of view, viscosity solutions have been proposed to this problem by [15] to find a physically meaningfull solution to this problem. From the numerical point of view, the existence and appareance of singu- larities in the solution implies that adequate numerical methods have to be used to obtain the unique viscosity solution. This key point has been largely discussed on uniformgridsorCartesianmeshesandad-hocschemeshavebeenproposedinvolving monotonicity, upwinding, essentially non-oscillatory (ENO) schemes and weighted essentially non-oscillatory (WENO) schemes, see for instance [34], [31]. The exten- sion of these schemes on unstructured isotropic meshes has also been investigated; for instance in [1], [28]. Inalltheseapproaches,theapproximationerrorisstronglyrelatedtothequality of the discretisation as well as to the mesh density in the vicinity of the interface. To this end, we propose the construction of a Riemannian metric tensor to help controlling the generation of anisotropic elements close to the interface. 3.2. Geometric approximation of the interface. While analytical (implicit) descriptions of the interface may be convenient, complex manifolds do not usually have such simple representations. The implicit representation can be stored with a discretization of the domain Ω ⊂ Rn in a finite number of sample points at which function φ is approximated. Since only the points near the manifold φ = 0 are important to accurately represent the interface, clustering the points near the interface will be the main challenge of approximation techniques. Intwodimensions,wehaveshown[20]thatitispossibletoconstructanaccurate piecewise linear approximation Γ of the manifold Γ(t) describing the interface h φ(x,t)=0. Let suppose the levelset function φ is known at the vertices of a given triangulation T of Ω. A discrete anisotropic metric tensor field can be defined h in order and the triangulation T can be adapted to create quasi-uniform meshes h in this metric. We consider as the approximation space of φ the space of affine Lagrange finite elements P (the basis functions being denoted as ϕ) and write the 1 approximation φ of φ for all x in Ω as (omitting the time variable): h np (cid:88) φ (x)= φ(x )ϕ (x), (28) h j j j=1 where np denotes the number of mesh vertices in T . Considering the set E of all h elements intersected by Γ(t), E ={K ∈T ,K∩Γ(t)(cid:54)=0}, and denoting by λ the h j