Jean-Daniel Boissonnat

Camille Wormser

Mariette Yvinec

Anisotropic Delaunay Mesh Generation

Anisotropic meshes are triangulations of a given domain in the plane or in higher dimensions, with elements elongated along prescribed directions. Anisotropic triangulations are known to be well suited for interpolation of functions or solving PDEs. Assuming that the anisotropic shape requirements for mesh elements are given through a metric field varying over the domain, we propose a new approach to anisotropic mesh generation, relying on the notion of anisotropic Delaunay meshes. An anisotropic Delaunay mesh is defined as a mesh in which the star of each vertex v consists of simplices that are Delaunay for the metric associated to vertex v.

This definition works in any dimension and allows to define a simple refinement algorithm. The algorithm takes as input a domain and a metric field and provides, after completion, an anisotropic mesh whose elements are sized and shaped according to the metric field.

Preliminaries

Anisotropic Metric

An anisotropic metric in R d is defined by a symmetric positive definite quadratic form represented, in some vector basis, by a d × d matrix M . The distance between two points a and b as measured by metric M , is defined as

This definition provides a definition for M -lengths and, by integration, for higher dimensional M -volume measures.

Introduction

Anisotropic meshes are triangulations of a given domain in the plane or in higher dimensions, with elements elongated along prescribed directions. Anisotropic triangulations have been shown to be particularly well suited for interpolation of functions [START_REF] D'azevedo | Optimal triangular mesh generation by coordinate transformation[END_REF][START_REF] Richard | What is a good linear finite element? Interpolation, conditioning, anisotropy, and quality measures[END_REF] and for solving PDEs [START_REF] Babuska | Error estimates for adaptive finite element computations[END_REF]. They allow to minimize the number of elements in the mesh while retaining a good accuracy in computations.

The required anisotropy is generally described through a metric field defined over the domain to be meshed. The directions along which the elements should be elongated are usually given, at each point of the domain, as a quadratic form. The eigenvectors and eigenvalues of the quadratic form describe the preferred directions and their anisotropic ratios.

Two main issues arise in this context. The first is to define the metric field. The second one is to generate a mesh whose elements are shaped according to the chosen metric field.

Defining good metric fields and error estimates is still an active research area. Alauzet et al. introduced the notion of continuous metrics and continuous meshes to minimize interpolation error [START_REF] Alauzet | Multi-dimensional continuous metric for mesh adaptation[END_REF][START_REF] Loseille | Optimal 3d highly anisotropic mesh adaptation based on the continuous mesh framework[END_REF][START_REF] Alauzet | Size gradation control of anisotropic meshes[END_REF]. Loseille et al. [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3d steady Euler equations[END_REF] applied this notion to a posteriori error estimates in order to minimize the approximation error during the process of solving some PDEs. Chen et al [START_REF] Chen | Optimal anisotropic meshes for minimizing interpolation errors in l p -norm[END_REF] considered anisotropic finite element approximation of functions in the L p norm. Their result reveals that the accuracy of the approximation is governed by a quantity that depends non-linearly on the hessian of the function. In his thesis, Mirebeau [START_REF] Mirebeau | Anisotropic finite element approximation[END_REF] extends this result to finite elements of arbitrary degree and to Sobolev norms, and provides sharp asymptotic error estimates for the approximation of functions of two variables.

Various methods have been proposed to generate anisotropic meshes whose elements are shaped according to a given metric field. In their early work on 2D meshes, Bossen and Heckbert [START_REF] Bossen | A pliant method for anisotropic mesh generation[END_REF] proposed to adapt their pliant method for mesh generation to the anisotropic setting. Starting from a constrained Delaunay triangulation, the pliant method performs local optimization operations including centroidal smoothing and retriangulation, and possibly insertion or removal of vertices. Li et al. [START_REF] Li | Biting ellipses to generate anisotropic mesh[END_REF] and Shimada et al. [START_REF] Shimada | Anisotropic Triangulation of Parametric Surfaces via Close Packing of Ellipsoids[END_REF] have proposed to place the mesh vertices close to the centers of ellipsoid bubbles optimally packed in the domain. Borouchaki et al. [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications. part I algorithms[END_REF] proposed to adapt the standard Delaunay incremental construction to the anisotropic context. This construction is then combined with an anisotropic version of the unit mesh approach that aims at producing meshes whose edges have unit length. Lengths, in the anisotropic case, are measured in the Riemanian metric provided by the metric field. The efficiency of the method has been demonstrated in various contexts [START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF][START_REF] Dobrzynski | Anisotropic Delaunay mesh adaptation for unsteady simulations[END_REF]. Following a different line of research, some attempts have been done recently to define anisotropic Delaunay triangulation and meshes as the duals of some Voronoi diagrams derived from the metric field. Labelle and Shewchuk [START_REF] Labelle | Anisotropic Voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF] have defined an anisotropic mesh as the dual of the so-called anisotropic Voronoi diagram. The sites of this diagram are the mesh vertices and the distance to a site is computed with respect to the metric attached to this site. In the 2-dimensional case, Labelle and Shewchuk have proposed a refinement algorithm that can provably produce anisotropic meshes. Their approach has somehow been simplified in [START_REF] Boissonnat | Anisotropic diagrams: Labelle and Shewchuk approach revisited[END_REF], leading to a direct computation of the dual mesh, and extended by Cheng et al. [START_REF] Cheng | Anisotropic surface meshing[END_REF] to produce anisotropic meshes of surfaces embedded in 3D. Extending Labelle and Shewchuk's approach to higher dimensions seems however difficult due to the presence of flat tetrahedra called slivers [START_REF] Schoen | Robust, guaranteed-quality anisotropic mesh generation[END_REF]. Du and Wang [START_REF] Du | Anisotropic centroidal voronoi tessellations and their applications[END_REF] have proposed to use a definition of anisotropic Voronoi diagrams which is somewhat symmetric to the definition used by Labelle and Shewchuk. The Voronoi regions are based on distances from points to sites that are computed with respect to the metric of the point. Du and Wang compute centroidal Voronoi diagrams using this definition and show experimentally that the dual structures are generally anisotropic meshes of high quality. However they could not provide theoretical guarantees or conditions that ensure that the dual structure is a valid triangulation.

In this paper, we introduce a new notion of anisotropic mesh which extends nicely in any dimension. The resulting meshes can be computed using standard Delaunay algorithms. As in the previous approaches, we assume that the anisotropy is prescribed by a metric field that associates to each point p of the domain a symmetric positive definite square matrix M p , describing the metric at point p. Given a set of points V called sites, we consider, for each site v ∈ V , the Delaunay triangulation Del v (V ) of V , computed for the metric M v attached to location of v. Each triangulation Del v (V ) is well defined and can be computed using the standard Euclidean Delaunay triangulation on affinely transformed input points. For each site v ∈ V , we keep the star S v of v in Del v (V ), i.e. the set of simplices of Del v (V ) that are incident to v. The collection of stars is called the star set of V . In general, there are inconsistencies among the stars : a simplex s may appear in the stars of some of its vertices without appearing in the stars of all of them. As a result, the simplices in the star set of V do not form a triangulation of V . However, we show in this paper that, given a compact domain of R d and a smooth metric field, one can insert new sites in V at carefully chosen locations so that all inconsistencies are removed. The simplices in the star set then form a d-dimensional triangulation that we call an anisotropic Delaunay mesh. When the domain has smooth boundaries, a faithfull representation of those boundaries may be obtained using the method of restricted Delaunay triangulations. The refinement algorithm is then extended to achieve also consistency between surface stars which are defined as the restrictions of the stars to the boundary surfaces. The algorithm produces then a mesh whose vertices lie within the input domain and whose boundary is within a controlled Hausdorff distance from the input domain boundary. Sharp features could possibly be handled using protecting balls but this issue is not handled in the present paper.

The idea of maintaining independent stars for each vertex of a mesh has been first proposed by Shewchuk [START_REF] Shewchuk | Star splaying: an algorithm for repairing Delaunay triangulations and convex hulls[END_REF] for maintaining triangulations of moving points. The star set has also been used [START_REF] Schoen | Robust, guaranteed-quality anisotropic mesh generation[END_REF] to build the dual of an anisotropic Voronoi diagram as defined by Labelle and Shewchuk [START_REF] Labelle | Anisotropic Voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF]. The method we use to ensure consistency among the stars is inspired by the work of Li and Teng [START_REF] Li | Generating well-shaped Delaunay meshed in 3d[END_REF][START_REF] Li | Generating well-shaped d-dimensional Delaunay meshes[END_REF] for removing slivers in isotropic meshes. In our context, the method is extended so as to take into account the metric distortion between neighboring stars and also to avoid, in addition to slivers, more general quasi-cospherical configurations that may prevent the termination of the algorithm.

In addition to conforming to the given anisotropic metric field, this mesh generation method has several notable advantages.

-It is not limited to the plane and works in any dimension; -It is easy to implement. Through a stretching transform, the star of each vertex in the mesh can be computed as part of an Euclidean Delaunay triangulation. Therefore the algorithm relies only on the usual Delaunay predicates (applied in some stretched spaces);

-The star of each vertex in the output mesh is formed with simplices that are Delaunay with respect to the metric of the central vertex. This provides a neat characterization of the output mesh from its set of vertices.

-The method provides some theoretical guarantees on the size and shape of the output mesh elements. Each element is guaranteed to be sized and shaped according to the metrics of all its vertices.

The main downside of this anisotropic Delaunay mesh approach is that no consistent mesh is obtained before reaching the very end of the refinement algorithm. This may leads to over dense meshes when the metric field is highly distorted. In such cases, the only way out consists in somehow smoothing the input metric field.

This paper is an extension of a preliminary work limited to the 3-dimensional case, that has been presented at the Symposium on Computational Geometry [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF].

The paper is organized as follows. Section 2 recalls basic facts about anisotropic metrics, metric fields, metric distortion and sizing fields. Section 3 introduces the main notions underlying our refinement strategy: restricted Delaunay triangulation, star sets, inconsistencies, slivers and quasi-cosphericities. Section 4 presents the anisotropic mesh generation algorithm. For pedagogical reason, we focus in this section on the generation of a mesh covering a given domain and conforming to a varying field of anisotropic metrics defined on this domain. We postpone to section 7 the additionnal problem to get into the mesh a faithfull representation of the domain boundaries and internal subdivisions. Sections 5 and 6 detail the proof that the refinement algorithm terminates. Section 7 explains how to handle domain boundaries and sharp features in the anisotropic setting. At last Section 8 provides concluding remarks and some insights on on-going and future work.

In the following, we often use the same notation, M , for a metric and the associated matrix in a given basis. Given the symmetric positive definite matrix M , we denote by F M any matrix such that det(F M ) > 0 and F T M F M = M . Note that F M is not unique. The Cholesky decomposition provides an upper triangular F M , while a symmetric F M can be obtained by diagonalizing the quadratic form M and computing the quadratic form with the same eigenvectors and the square root of each eigenvalue.

Note that

d M (a, b) = (a -b) T F T M F M (a -b) = F M (a -b) (1)
where the notation . stands for the Euclidean norm. Equation (1) proves that d M is a distance and, in particular, enjoys the standard triangular inequality. In the following we call such a F M matrix a stretching transform of M .

Given some metric M , an M -sphere C M (c, r), with center c and radius r, is defined as the set of points p such that d M (c, p) = r, and likewise an M -ball B M (c, r), is defined as the set of points p such that d M (c, p) ≤ r. Note that an M -sphere is an ellipsoid in the Euclidean space, with its axes aligned along the eigenvectors of M .

Given a k-simplex s in R d and a metric M , we define the M -circumsphere C M (s) as the circumscribing M -sphere of s with smallest radius. The M -circumball B M (s) is the M -ball bounded by C M (s) and the M -circumradius r M (s) of a simplex s is the radius of C M (s).

Equation (1) shows that C M (s) is the reciprocal image F -1 M (C(F M (s))) of the Euclidean circumscribing sphere of the simplex F M (s).

Let M be a metric and V be a set of points, called sites. The Delaunay triangulation of V for metric M , denoted Del M (V ), is the triangulation of V such that the interior of the M -circumball of each d-simplex is empty, i.e. contains no site of V . Owing to Equation (1), the Delaunay triangulation Del M (V ) of a finite set of points V for metric M is simply obtained by computing the Euclidean Delaunay triangulation of the stretched image F (V ) = {F M v, v ∈ V }, and stretching the result back with F -1 M . The Delaunay triangulation Del M (V ) is thus viewed as the dual of a stretched Voronoi diagram.

Metric Field and Distortion

In the rest of the paper, we consider a compact domain D ⊂ R d and assume that we are given a metric field defined over D, i.e. a metric M x is given at each point x ∈ D.

In the following, to avoid double subscripts, we replace subscript M x by x and simply write Y x for Y Mx . Hence, we will write for instance

F x for F Mx and d x (a, b) for d Mx (a, b).
We recall some definitions due to Labelle and Shewchuk [START_REF] Labelle | Anisotropic Voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF].

Given two metrics M and N , with stretching transforms F M and F N respectively, the distortion γ(M, N ) between M and N is defined as

γ(M, N ) = max{ F -1 M F N , F -1 N F M },

Ax

x . Note that the distortion γ(M, N ) does not depend the streching matrices F M and F N choosen for the metrics M and N . In the context of a metric field, the relative distortion between two points p and q of the domain D is defined as γ(p, q) = γ(M p , M q ). Observe that γ ≥ 1 and is equal to 1 iff M p = M q .

A fundamental property of γ(p, q) is that it bounds the ratio between d p and d q :

∀x, y, 1 γ(p, q) d q (x, y) ≤ d p (x, y) ≤ γ(p, q) d q (x, y). (2)
Let s = p 0 p 1 . . . p d be d-simplex, let M i be the metric attached the vertex p i , for i = 0, ..., d and let B i (s) be the M i -circumball of s. The distortion γ(B M ) of a M -ball B M is defined as the maximal distortion between any pairs of points of B M ∩ D. We define the distortion γ(s) of a simplex s as the maximum of the distortion of its circumballs:

γ(s) = max{γ(B i (s)), i = 0, . . . , d}.

Sizing field

In this paper, we will assume that the metric field is smooth over the domain D. The distortion γ(p, q) is then a continuous function and the maximum distortion over D, Γ = sup x,y∈D γ(x, y), is finite since D is compact.

We now consider a local view of the distortion. Given a constant γ 0 > 1, called the distortion bound, we define for each point p ∈ D the bounded distortion radius, bdr(p, γ 0 ), as the upper bound on distances such that, for all q and r in D, max(d p (p, q), d p (p, r)) ≤ ⇒ γ(q, r) ≤ γ 0 .

Lemma 2.1 (Bounded distortion radius lemma) The bounded distortion radius bdr(p, γ 0 ) enjoys the following property for any p, q in D:

1 γ(p, q) [bdr(p, γ 0 ) -d p (p, q)] ≤ bdr(q, γ 0 ) ≤ γ(p, q) [bdr(p, γ 0 ) + d p (p, q)] .
Proof Let x, y be any two points in D so that:

d q (q, x) ≤ 1 γ(p, q) (bdr(p, γ 0 ) -d p (p, q)) , (3)
d q (q, y) ≤ 1 γ(p, q) (bdr(p, γ 0 ) -d p (p, q)) . (4)
Then, we have, using the triangular inequality,

d p (p, x) ≤ d p (p, q) + d p (q, x) ≤ d p (p, q) + γ(p, q)d q (q, x) ≤ bdr(p, γ 0 )
and, similarly,

d p (p, y) ≤ bdr(p, γ 0 ).
Then, by definition of the bounded distortion radius, γ(x, y) ≤ γ 0 . Because the last inequality is true for any pair of points x, y satisfying inequalities (3) and ( 4), we conclude that 1 γ(p, q)

[bdr(p, γ 0 )d p (p, q))] ≤ bdr(q, γ 0 ).

(

) 5
To prove the second inequality of Lemma 2.1, we simply write inequality (5) for the pair (q, p), which yields:

1 γ(p, q)
[bdr(q, γ 0 )d q (p, q)] ≤ bdr(p, γ 0 ) from which we deduce bdr(q, γ 0 ) ≤ γ(p, q) bdr(p, γ 0 ) + d q (p, q) ≤ γ(p, q) [bdr(p, γ 0 ) + d p (p, q)] .

We will further assume that the bounded distortion radius has a strictly positive lower bound on domain D: bdr 0 = min p∈D bdr(p, γ 0 ) > 0.

Our algorithm uses a sizing field to control the size of mesh elements according to the local metric. The most basic sizing field we use is the bounded distorsion radius bdr(p, γ 0 ) which will enforce the mesh density to adapt to the metric distorsion. However, our algorithm can take into account additionnal user defined sizing criteria.

Definition 2.2 (Sizing field) Let γ 0 ≥ 1 be a given distortion bound. We call sizing field and denote by sf(p, γ 0 ) (or sf(p) for short if γ 0 is understood), any function defined over the domain D, that satisfies the three following conditions:

positiveness ∃ sf 0 > 0 such that ∀x ∈ D, sf(x, γ 0 ) ≥ sf 0 (6) distortion ∀x ∈ D, sf(x, γ 0 ) ≤ bdr(x, γ 0 ) (7) continuity ∀x, y ∈ D, (8) 1 γ(x, y) [sf(x, γ 0 ) -d x (x, y)] ≤ sf(y, γ 0 ) ≤ γ(x, y) [sf(x, γ 0 ) + d x (x, y)]

Stars and Refinement

We now define the local structures that are built and refined by our algorithm. These definitions rely on the notion of restricted Delaunay widely used in reconstruction area, see e.g [START_REF] Chew | Guaranteed-quality mesh generation for curved surfaces[END_REF][START_REF] Edelsbrunner | Triangulating topological spaces[END_REF][START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF][START_REF] Tamal | Curve and surface reconstruction[END_REF].

Let D be as before a domain of R d and let V be a finite set of points of D that are called hereafter sites or vertices.

The restriction to D of the Delaunay triangulation Del(V ) of V is the subcomplex of Del(V ) whose maximal faces are the d-simplices of Del(V ) that have their dual Voronoi vertices inside the domain D. The natural extension of this definition in the anisotropic setting is to define the restriction of Del M (V ) to D as the subcomplex of Del M (V ) whose maximal faces are the d-simplices s of Del M (V ) that have their M -circumcenter inside D.

Stars and Inconsistencies

For each site v in V , we consider the Delaunay triangulation Del v (V ) of V for the metric M v . We define the restricted star S v of site v as the set of d-simplices incident to v in the restriction of Del v (V ) to D.

The collection of all restricted stars S(V

) = {S v , v ∈ V }, is called the restricted star set of V .
We say that a simplex is inconsistent if it appears in the star of at least one of its vertices but does not appear in the stars of all of them. For instance, in Figure 1, edge vw or facet xvw are inconsistent because they appear in the triangulation Del v (V ) computed with metric M v but not in the Delaunay triangulations Del w (V ) computed with metrics M w .

Our algorithm incrementally inserts new sites in V and updates the restricted star set S(V ) until it contains no more inconsistent simplices. As shown below, when the mesh is dense enough with respect to the variation of the metric field, inconsistencies are related to the occurence of special configurations of subsets of sites that are called quasi-cospherical configurations or QC-configurations for short. The algorithm will therefore aim at avoiding those QC-configurations. As it turns out, QC-configurations can be avoided but only when even more special configurations called slivers do not occur. Both notions are now defined.

Slivers

In the following we use the definition of slivers provided by Cheng et al [START_REF] Cheng | Silver exudation[END_REF] and we extend it in the anisotropic setting. Let s be a k-simplex. We denote by C M (s) the M -circumsphere of s, by r M (s) the Mcircumradius of s, by e M (s) the M -length of the shortest edge of s for the metric M and by Vol M (s) the M -volume of s. We define two quality measures of s for metric M . The M -radius-edge ratio is defined as the ratio ρ M (s) = r M (s)/e M (s). The sliverity ratio σ M (s) is the ratio Vol M (s)/e k M (s)

1 k .
Definition 3.1 (Sliver) Let ρ 0 and σ 0 be two positive constants and let M be a metric. A k-simplex s is said to be

• well-shaped for M , if ρ M (s) ≤ ρ 0 and σ M (s) ≥ σ 0 • a sliver for M , if ρ M (s) ≤ ρ 0 , σ M (s) < σ 0
• a k-sliver for M , if it is a sliver and all its (k -1)-dimensional faces are well-shaped.

It is easily shown that any k-dimensional simplex that is a sliver is either a k-sliver or include as a subface a k -sliver for some k < k.

The following lemma is known for slivers in dimension 3, see e.g. [START_REF] Edelsbrunner | Smoothing and cleaning up slivers[END_REF]. It has been extended to higher dimensions [START_REF] Li | Generating well-shaped d-dimensional Delaunay meshes[END_REF] and extends naturally to anisotropic metrics as proved in the appendix. If s is a k-sliver wih respect to M , the M -distance from v to aff(s(v)) is at most 2kσ 0 r M (s(v)) and the M -distance from v to C M (s(v)) ∩ aff(s(v)) is at most 4πkρ 0 σ 0 r M (s(v)).

Quasi-Cosphericity

Let γ 0 > 1 be a bound on the distortion and M be a metric. We now introduce the notion of (γ 0 , M )-cosphericity and show its link with inconsistent simplices.

Definition 3.3 (QC-configuration) A subset U of d + 2 sites {p 0 , p 1 , . . . , p d+1 } is said to be a (γ 0 , M )-cospherical configuration if there exist two metrics N and N such that :

• γ(M, N ) ≤ γ 0 , γ(M, N ) ≤ γ 0 and γ(N, N ) ≤ γ 0 ;

• the triangulations Del N (U ) and Del N (U ) are different.

Metrics N and N are said to witness the (γ 0 , M )-cosphericity of U . If M is clear from the context, we simply say that U is a γ 0 -cospherical configuration and if both M and γ 0 are understood, we say that U is a quasi-cospherical configuration or a QC-configuration for short.

See Figure 2 for an illustration in the plane. Note that the d+2 points in U play symmetric roles in the above definition. In the sequel, U will often consist of the set of vertices of a d-simplex s belonging to the star S v of some site v ∈ V , together with an additionnal site p of V . In such a case, we write U = (s, p).

From Radon theorem, there are only two distinct triangulations of U = (s, p) and any d-simplex with vertices in U belongs to exactly one of them [START_REF] Charles L Lawson | Properties of n-dimensional triangulations[END_REF][START_REF] Richard | Updating and constructing constrained delaunay and constrained regular triangulations by flips[END_REF]. Therefore, we have the following easy lemma.

Lemma 3.4 A configuration (s, p) is (γ 0 , M )-cospherical iff there exist two metrics N and N such that

• γ(M, N ) ≤ γ 0 , γ(M, N ) ≤ γ 0 and γ(N, N ) ≤ γ 0 ;
• p belongs to the interior of exactly one of the two circumballs B N (s) and B N (s).

The following lemma relates QC-configurations and inconsistencies.

Lemma 3.5 Let s be an inconsistent simplex of star S v . If γ(s) < γ 0 , then there exists a site q ∈ V such that the configuration (s, q) is (γ 0 , M v )-cospherical. Proof By definition s appears in the Delaunay triangulation Del v (V ) computed with the metric of vertex v, but not in the triangulation Del w (V ) computed with the metric of some other vertex w of s. Take N = M v and N = M w . Because the distortion of s is less than γ 0 , we have γ(v, w) = γ(M v , M w ) ≤ γ 0 . Since s is a d-simplex in S v but not in S w , it belongs to Del v (V ) and not to Del w (V ). Hence, there is a site q ∈ V such that q is inside B w (s) and not inside B v (s). It then follows from Lemma 3.4 that (s, q) is a (γ 0 , M v )-cospherical configuration witnessed by the metrics N = M v and N = M w .

Given a metric M and a (γ 0 , M )-cospherical configuration U , the M -radius r M (U ) of U is defined as the minimum of the M -circumradii of all the d-simplices with vertices in U .

Definition 3.6 (Well-shaped QC-configuration) A (γ 0 , M )-cospherical configuration U is said to be well-shaped if any d-simplex formed with vertices in U is well-shaped with respect to M .

4 Meshing Algorithm

Algorithm Outline

To mesh a given compact domain D, the algorithm constructs the set of sites from a small set of initial sites by inserting new sites in a greedy way. The algorithm maintains the star set for the current set of sites and the addition of new sites is steered by this star set : while there remain bad simplices in the star set, the algorithm selects one bad simplex and kills this simplex by inserting a new site. Bad simplices are d-simplices that have a high distortion or are oversized with respect to a user defined sizing field, those that are badly shaped (high radius-edge ratio or small sliverity ratio), and those that are inconsistent. To kill a bad simplex s appearing in a star S v , a new site p, called the refinement point, is chosen in the M v -circumscribing ball of s and inserted in the star set. The maintainance of the star set upon the insertion of the refinement point p involves the creation of a new star S p for p and the insertion of p in the star S v and in any star S w where p will appear as a vertex. Note that a new site p has to be inserted in star S w iff p is in conflict with some d-simplex of S w , where point p is said to be in conflict with the d-simplex s of S w if p is included in the M w -circumball of s. Upon each insertion, the algorithm maintains the star set by calling the following Insert procedure:

Algorithm 1 Insert(p)
1. insert p in all the stars S w that contain a simplex in conflict with p;

2. create the new star S p .

As noticed in Section 3.3, once the set of vertices is dense enough with respect to the variation of the metric field so that all simplices in the star set have a distortion smaller then γ 0 , inconsistenties arise only from QC-configurations. The refinement algorithm therefore aims at avoiding those configurations. However, as will be clear from the proof of Lemma 5.3, it is not possible to avoid QC-configurations involving slivers. The algorithm thus needs to remove slivers before removing inconsistent simplices.

Recall that for a d-simplex s in some star S v , we write B v (s) or B v (c v (s), r v (s)) for the M v -circumball of s with center c v (s) and radius r v (s), ρ v (s) for the M v -radius-edge ratio of s and σ v (s) for its M v -sliverity ratio. The refinement algorithm (see Algorithm 2) applies four rules in turn. The rules are applied with a priority order : rule (i) is applied only if no rule (j) with j < i can be applied. The algorithm ends when no rule applies any more. The algorithm relies on two procedures: procedure Insert inserts a new site in the data structures, and procedure Pick valid chooses the location of the new site (see the next section).

The refinement algorithm depends on parameters α 0 , ρ 0 , σ 0 , and γ 0 while the Pick valid procedure depends on two more parameters β and δ. The values of constants α 0 , ρ 0 , σ 0 , γ 0 control the quality of the mesh elements and their adaptation to the metric field. Parameters β and δ influence the behaviour of the algorithm and their values are chosen in Section 5 in order to ensure the termination of the refinement algorithm.

Remark. Note that the sizing field sf(p) used in Rule (1) takes care of the distortion bound γ 0 and also possibly of a user defined sizing field (see Definition 2.2). Parameter α 0 is always chosen less than 1. Therefore, when Rule (1) does not apply anymore, the distortion of any d-simplex in any star is bounded by γ 0 .

Algorithm 2 Refinement algorithm

Rule (1) Size:

If ∃ a d-simplex s in star S v such that r v (s) ≥ α 0 sf(c v (s)), Insert(c v (s)); Rule (2) Radius-edge ratio: If ∃ a d-simplex s in star S v such that ρ v (s) > ρ 0 , Insert(Pick valid(s, M v )); Rule (3) Sliver removal: If a d-simplex s in star S v is a M v sliver (i.e. ρ v (s) ≤ ρ 0 and σ v (s) < σ 0 ), Insert(Pick valid(s, M v )); Rule (4) Inconsistency: If a d-simplex s in some star S v is inconsistent, Insert(Pick valid(s, M v ));
Sections 5 and 6 will prove that the algorithm terminates. Before that, Subsection 4.2 describes the procedure Pick valid while Subsection 4.3 analyses the properties of the resulting mesh.

Picking Region and Hitting Sets

In this section, we describe in more detail procedure Pick valid. The simplest idea to kill a simplex would be to insert a refinement point at its circumcenter. However, with this simple strategy, the algorithm may loop, creating cascading configurations of slivers and QCconfigurations, and is not guaranteed to terminate. To avoid slivers and QC-configurations, the algorithm resorts to a strategy analog to the one used by Li and Teng [START_REF] Li | Generating well-shaped Delaunay meshed in 3d[END_REF][START_REF] Li | Generating well-shaped d-dimensional Delaunay meshes[END_REF] to avoid slivers in isotropic meshes. The basic idea is to relax the choice of the refinement point of a bad simplex. Instead of using systematically the circumcenter, the refinement point of a bad simplex is picked from a small region around the circumcenter, called the picking region. The refinement point is carefully chosen in the picking region so as to avoid the formation of new slivers and new QC-configurations.

Definition 4.1 (Picking region) Let δ < 1 be a constant called the picking ratio. If s is a bad simplex in star S v , with M v -circumball B v (c v (s), r v (s)), the M v -picking region of s, noted P v (s), is the M v -ball B v (c v (s), δr v (s)).
In fact, it is not possible, when choosing a refinement point in the picking region P v (s) of a simplex s of S v to completely avoid the formation of new slivers and new QC-configurations.

The Pick valid procedure will only avoid the creation of small slivers and small QCconfigurations where the meaning of small, precisely defined below, is relative to the radius r v (s) and controlled by a parameter β.

Definition 4.2 (Hitting set) Let p be a point in the M v -picking region of a simplex s.

Let r v (s) be the M v -circumradius of s and β be a constant. A subset t of the current set of sites V is said to hit p if one of the two following conditions is satisfied:

• t consists of k ≤ d sites and, for some metric M such that γ(M p , M ) ≤ γ 0 , the k-simplex s = (t, p) is a k-sliver with M -circumradius r M (s ) ≤ βr v (s).
• t consists of d + 1 sites and, for some metric

M such that γ(M p , M ) ≤ γ 0 , U = (t, p) is a well-shaped (γ 0 , M )-cospherical configuration with M -radius r M (U ) ≤ βr v (s).
A point p in P v (s) is said to be a valid refinement point if it is not hit by any subset of V . Each subset t of sites in V induces a forbidden region where the refinement point should not lie in order to avoid being hit by t. A subset t of sites in V that hits some point in the picking region P v (s) is said to be a hitting set for P v (s). A point p in P v (s) is therefore valid if it avoids the forbidden region of any hitting set of P v (s).

v P v (s) q r u p s Figure 3: {q, r, u} is a hitting set for the picking region P v (s). It defines a forbidden region (dashed area) to be avoided by the refinement point p of simplex s.

Note that the definition of valid refinement points depends on the constants δ and β: δ defines the size of the picking regions and β bounds from below the size of acceptable new slivers and new QC-configurations, with respect to the circumradius of the simplex being refined. The definition of valid refinement points also depends on the constants ρ 0 and σ 0 that define well-shaped simplices and slivers, and on the constant γ 0 that defines QC-configurations. We prove in the following sections (Section 5 and Section 6) that it is possible to choose the algorithm parameters β, δ, ρ 0 , σ 0 , γ 0 and α 0 so that valid refinement points do exist in any picking region considered by the refinement algorithm.

To find a valid refinement point in the M v -picking region P v (s) of some bad d-simplex s, the insertion algorithm calls the following Pick valid procedure. This procedure randomly chooses a point in the picking region P v (s) until it finds one that avoids all forbidden regions. This procedure depends on constants ρ 0 , σ 0 , γ 0 , δ and β, to be fixed later in Section 5.

Algorithm 3 Pick valid(s, M v )
Step 1 Pick randomly a point p in the picking region P v (s)

Step 2 Avoid small slivers For k = 3 to d, if there exists a subset of k sites in V that hits p, then discard p and go back to step 1.

Step 3 Avoid small QC-configurations If there exists a subset of d + 1 sites in V that hits p, then discard p and go back to step 1.

Step 4 Return p.

Quality of the final mesh

Upon termination of the algorithm, all stars are consistent. Therefore they can be merged together to form a triangulation T of the domain. Each simplex s in T is well-shaped with respect to the metric of all its vertices, i.e. ρ p (s) ≤ ρ 0 and σ p (s) ≥ σ 0 for any vertex p of s. Moreover each simplex s in T complies to the sizing field sf(γ 0 ) which implies that s has a distortion smaller than γ 0 . If needed, the sizing field sf(γ 0 ) may also take into account a user defined sizing field.

Complexity of the meshing algorithm

The complexity of algorithm 2 is roughly linear with respect to the number of vertices of the output mesh. This might appear suprising since the algorithm maintains a distinct 3D triangulation for each vertex in the mesh. In fact, these triangulations are quite small, since they are designed to maintain the star of a single central vertex. In each star, the generated sites form a well spaced set of points with respect to the local metric, and therefore the number of simplices in the star is bounded by a constant depending only on the dimension d. Hence the star set is a data structure whose size is linear with respect to the number of vertices of the output mesh.

Each new vertex is inserted in a constant number of stars and its own stars is initialized by the insertion of a constant number of the current vertices. The insertion of a point in a star takes a constant time. One of the main concerns of our implementation is to efficiently filter out the stars in which a new vertex has to be inserted, and the subset of current vertices that have to be inserted in the star of the new vertex. These problems are handled using additionnal standard data structures for range queries among bounding boxes.

Finally, the most costly part of the algorithm is the computation of a valid refinement point using the pick valid procedure. This amounts to randomly choose candidates and check if they are valid. The validity check amounts to simulating the insertion of the point and its cost is similar to the cost of an insertion. The proof of the Picking Lemma below, yields that the expected number of performed trials is constant. Therefore, altogether the expected complexity of the algorithm is linear with respect to the number of vertices of the ouput mesh.

It now remains to prove that the algorithm terminates, which will be done in the two next sections.

Termination of the Algorithm

The refinement algorithm (Algorithm 2) depends on parameters α 0 , γ 0 , ρ 0 , σ 0 that respectively control the size, the distortion, the radius-edge ratio and the sliverity ratio of the simplices, and on parameters δ and β that define the picking-region and valid refinement points. In this section and in the following one, we prove that for a suitable choice of those parameters, Algorithm 2 terminates providing as claimed a consistent mesh that is an anisotropic Delaunay mesh.

Let us first notice that our algorithm will never refine a star element that is not part of the restricted Delaunay triangulation of the domain to be meshed. As a consequence, the Steiner vertices inserted by the algorithm are within the domain, or very close to it when they are chosen through a call to the Pick valid procedure. The proof of termination then relies on a volume argument based on a minimum separation distance between any two vertices of the mesh.

For any vertex p in V , we define the separation distance and insertion-radius as follows.

Definition 5.1 If p is a vertex in V , the separation distance sd(p) of p is defined as:

sd(p) = min q∈V d p (p, q)
Definition 5.2 If p is a vertex of V and V (p) ⊂ V the subset of vertices inserted before p, the insertion radius r(p) of p is defined as:

r(p) = min q∈V (p) d p (p, q).
We will mainly show that there is a constant Λ (depending on parameters α 0 , γ 0 , ρ 0 , σ 0 , δ and β and on properties of the domain and of the metric field) such that the separation bound sd(p) ≥ Λ sf(p) holds for any site in V .

We first need to ensure that valid refinement points will be found in any picking region considered by the algorithm. This is the goal of the next lemma whose proof is deferred to Section 6).

Lemma 5.3 (Picking lemma) For any values of parameters β, δ, α 0 and ρ 0 , it is possible to choose σ 0 small enough and the distortion bound γ 0 close enough to 1, so that, if the separation bound sd(p) ≥ Λ sf(p) holds for any site in the current set V , valid refinement points do exist in the M v -picking region of any bad simplex s in star S v .

To prove the separation bound on vertices, we first prove a lower bound on the insertion radii of vertices, considering in turn each of the refinement rules. We begin with a technical lemma relating the circumradius of a simplex with the insertion radius of its refinement point.

Lemma 5.4 (Insertion radius lemma) Let s be a d-simplex of star S v with the M vcircumball B v (c v (s), r v (s)). Assume that s is a bad simplex. Let p be the refinement point of s and r(p) the insertion radius of p.

• If Rule (1) applies, the refinement point p of s is the circumcenter c v (s), and

r(p) ≥ r v (s) Γ , ( 9
)
where Γ is the maximal distortion over D: Γ = max x,y∈D γ(x, y).

• If one of Rule (2), (3) or ( 4) is applied, the refinement point p is taken from the picking region P v (s) and :

r(p) ≥ 1 -δ γ 0 r v (s). ( 10
)
Proof In the first case, p = c v (s), and therefore min

q∈V (p) d v (p, q) ≥ r v (s) r(p) = min q∈V (p) d p (p, q) ≥ r v (s) Γ .
In the second case, p belongs to the picking region P v (s), and we know that the distorsion γ(s), hence also the distorsion γ(v, p), are at most γ 0 min q∈V (p)

d v (p, q) ≥ (1 -δ)r v (s) r(p) = min q∈V (p) d p (p, q) ≥ 1 -δ γ 0 r v (s).
Lemma 5.5 (Rule (1) lemma) When Rule ( 1) is applied, the insertion radius r(p) of the inserted site p is at least:

r(p) ≥ Λ 1 sf(p) with Λ 1 = α 0 Γ . (11)
Proof Rule ( 1) is applied to a simplex s in star S v when the M v -circumradius r v (s) of s is greater than α 0 sf(c v (s)). The refinement point p is then c v (s) and we get from Lemma 5.4

r(p) ≥ r v (s) Γ ≥ α 0 sf(c v (s)) Γ = α 0 Γ sf(p). ( 12
)
Lemma 5.5 proves a lower bound on the insertion radius of any vertex p introduced by application of Rule (1) . The next lemmas aim at finding a constant Λ 2 , and some conditions on α 0 , ρ 0 , γ 0 , β and δ so that Rules (2)-(4) will maintain the invariant that the insertion radius of any inserted point is at least Λ 2 sf(p).

Lemma 5.6 (Rule (2) lemma) Let s ∈ S v be a simplex to be refined by application of Rule (2) and let p be the refinement point of s. If, for any vertex q inserted before p, r(q) ≥ Λ 2 sf(q) then we have r(p) ≥ Λ 2 sf(p), provided that the following conditions hold

(1 -δ)ρ 0 γ 3 0 ≥ 2, ( 13
) (1 + δ)ρ 0 γ 0 Λ 2 ≤ 1. ( 14
)
Proof First, observe that γ(s) ≤ γ 0 since Rule (1) does not apply. Then, because p is inserted by application of Rule ( 2), the M v -circumradius of s, r v (s), is such that r v (s) ≥ ρ 0 e v (s), where e v (s) is the M v -length of the M v -shortest edge of s, which is the shortest edge of s according to metric M v . Let qq be the M v -shortest edge of s and q be the last inserted vertex of qq .

r v (s) ≥ ρ 0 e v (s) = ρ 0 d v (q, q ) ≥ ρ 0 γ 0 d q (q, q ) ≥ ρ 0 γ 0 r(q)
Using the induction hypothesis and the continuity condition of the sizing field (see Definition 2.2), we have

r v (s) ≥ ρ 0 γ 0 Λ 2 sf(q) ( 15
)
≥ ρ 0 γ 0 Λ 2 [sf(p) -d p (p, q)] γ(p, q) ≥ ρ 0 γ 2 0 Λ 2 [sf(p) -d p (p, q)] ≥ ρ 0 γ 2 0 Λ 2 [sf(p) -γ 0 d v (p, q)] (since v, p ∈ B v (s)). (16)
Now, because q is a vertex of s and p is chosen in the picking region P v (s), d v (p, q) ≤ (1 + δ)r v (s) which, together with inequality (16), gives:

r v (s) ≥ ρ 0 γ 2 0 Λ 2 [sf(p) -γ 0 (1 + δ)r v (s)] r v (s) ≥ ρ 0 γ 2 0 Λ 2 sf(p) 1 + ρ 0 γ 0 (1 + δ) Λ 2 . ( 17
)
Then, using the insertion radius lemma (Lemma 5.4), we get:

r(p) ≥ (1-δ)ρ 0 γ 3 0 Λ 2 sf(p) 1 + ρ 0 γ 0 (1 + δ) Λ 2 , ( 18
)
which proves that r(p) ≥ Λ 2 sf(p) when conditions ( 13) and ( 14) are fulfilled.

Lemma 5.7 (Rule (3)-( 4) lemma) Let s ∈ S v be a simplex to be refined by application of Rule (3) or Rule (4), and let p be the refinement point of s. If, for any vertex q inserted before p, r(q) ≥ Λ 2 sf(q) then we have r(p) ≥ Λ 2 sf(p), provided that the following conditions hold

β(1 -δ) γ 3 0 (1 + δ) ≥ 2 (19) Λ 2 ≤ (1 -δ)Λ 1 2γ 3 0 + γ 2 0 (1 + δ)Λ 1 ( 20
)
βΛ 2 γ 0 ≤ 1. ( 21
)
Proof Assume first that s was created by application of Rule (1). Then, if q is the last inserted vertex of s, we have r(q) ≥ Λ 1 sf(q) by Lemma 5.5. Furthermore, r v (s) is at least half the M v -length of any edge of s and, in particular, of any edge of s that is incident to q. Therefore, using the fact that v and q belong to s which has a low distorsion, we get

r v (s) ≥ min q ∈s d v (q, q ) 2 ≥ min q ∈V d q (q, q ) 2γ 0 ≥ r(q) 2γ 0 ≥ Λ 1 2γ 0 sf(q)
Now, repeating the calculations performed to deduce inequality ( 17) from ( 15), we obtain

r v (s) ≥ Λ 1 2γ 2 0 [sf(p) -γ 0 (1 + δ)r v (s)] r v (s) ≥ Λ 1 2γ 2 0 sf(p) 1 + Λ 1 (1+δ) 2γ 0 = Λ 1 sf(p) 2γ 2 0 + γ 0 (1 + δ)Λ 1 . (22)
Then, using the insertion radius lemma (Lemma 5.4), we get:

r(p) ≥ (1 -δ)Λ 1 sf(p) 2γ 3 0 + γ 2 0 (1 + δ)Λ 1 . ( 23
)
It follows that the bound r(p) ≥ Λ 2 sf(p) holds, provided condition (20) is satisfied. Now consider the case where s was created by application of Rule (2), (3) or (4). Assume that s has been created when inserting the refinement point q of a simplex s in some star S w (see Figure 4). The refinement point q was chosen by the procedure Pick valid(s , M w ) and therefore, r v (s) ≥ βr w (s ). Let us bound r w (s ) from below. Vertex q is the last inserted vertex of s. It has been chosen in the picking region of s and therefore the vertices of s are at M w -distance at most (1 + δ)r w (s ) from q. Hence, since q and w belong to s and γ(s ) ≤ γ 0 , r(q) ≤ γ 0 (1 + δ)r w (s ).

Therefore: Repeating the calculations performed to deduce inequality ( 17) from ( 15), we obtain

r v (s) ≥ βr w (s ) ≥ β γ 0 (1 + δ) r(q). r v (s) ≥ βΛ 2 γ 0 (1 + δ) sf(q).
r v (s) ≥ βΛ 2 γ 2 0 (1 + δ) [sf(p) -γ 0 (1 + δ)r v (s)]
Hence,

r v (s) ≥ βΛ 2 γ 2 0 (1+δ) sf(p) 1 + βΛ 2 γ 0 , (24)
and, from the insertion radius lemma (Lemma 5.4),

r(p) ≥ β(1-δ)Λ 2 γ 3 0 (1+δ) 1 + βΛ 2 γ 0 sf(p). (25)
Conditions ( 21) and [START_REF] Tamal | Curve and surface reconstruction[END_REF] ensure that r(p) ≥ Λ 2 sf(p).

Lemma 5.8 (Separation bound) Assume that each vertex p has been inserted in the set V with an insertion radius r(p) such that r(p) ≥ Λ 2 sf(p) where Λ 2 is a constant. Then the set V admits the following separation bound :

sd(p) ≥ Λ sf(p) with Λ = Λ 2 /Γ 2 1 + Λ 2 /Γ 2 .
Proof Observe that for any pair of vertices p, q ∈ V , we have either d p (p, q) ≥ r(p) or d q (p, q) ≥ r(q), where the first is true if p has been inserted after q and the second is true otherwise. In the first case, we have

d p (p, q) ≥ r(p) ≥ Λ 2 sf(p) ≥ Λ sf(p)
In the second case, we have

d p (p, q) ≥ d q (p, q) Γ ≥ Λ 2 sf(q) Γ ≥ Λ 2 Γ 2 (sf(p) -d p (p, q)) (using (8) ) d p (p, q) ≥ Λ 2 /Γ 2 1 + Λ 2 /Γ 2 sf(p) = Λ sf(p)
which proves the separation bound.

We can now give the main theorem of this section that proves the separation bound on the set of vertices and ensures that the refinement algorithm terminates. Theorem 5.9 Given a compact domain D and a sizing field over D satisfying conditions ( 6)-( 8), assume that the constants ρ 0 , γ 0 , β and δ are chosen in such a way that:

(1 -δ)ρ 0 γ 3 0 ≥ 2 (26) (1 -δ)β γ 3 0 (1 + δ) ≥ 2. ( 27
)
Assume furthermore that σ 0 is small enough and γ 0 close enough to 1 so that the picking lemma (Lemma 5.3) holds. Then the refinement algorithm (Algorithm 2) terminates.

Proof Observe that the inequalities 26 and 27 are just conditions ( 13) and ( 19) of Lemma 5.6 and 5.7. Assume that these inequalities hold.

We choose the value of Λ 2 small enough so that Λ 2 ≤ Λ 1 = α 0 Γ and condition (14) of Lemma 5.6, and conditions ( 20) and ( 21) of Lemma 5.7 hold. For further reference, let us notice here that it suffices to choose

Λ 2 = min 1 (1 + δ)ρ 0 , (1 -δ)α 0 16Γ + 4(1 + δ)α 0 , 1 β (28)
to get a value of Λ 2 , independant of σ 0 and γ 0 , and suitable for any value of σ 0 and any

γ 0 ∈ [1, 2].
We first prove by induction that any vertex p is inserted in V with an insertion radius r(p) ≥ Λ 2 sf(p). First notice that the induction hypothesis can be enforced on the set of initial points. Then assume that the hypothesis is true up to a given stage. From the separation bound lemma (Lemma 5.8), any vertex in the current set has a separation bound sd(p) ≥ Λ sf(p). From the picking lemma (Lemma 5.3), we know that the algorithm will find a valid refinement point if the next vertex is to be searched in a picking region. Now, from Lemmas 5.5, 5.6, and 5.7, we know that the insertion radius of the next vertex p is still going to be bounded by r(p) ≥ Λ 2 sf(p) which achieves the inductive proof.

Finally, we can set up the volume argument for the algorithm termination in the metric M y of any point y ∈ D. Indeed, for any pair p, q of vertices in V , we have either

d p (p, q) ≥ r(p) ≥ Λ 2 sf(p) ≥ Λ 2 sf 0 or d q (p, q) ≥ r(q) ≥ Λ 2 sf(q) ≥ Λ 2 sf 0 .
In both cases, we may conclude that

d y (p, q) ≥ Λ 2 Γ sf 0 .
Since D is a compact domain and has therefore a bounded M y -volume, this proves that the algorithm can only insert a finite number of vertices and therefore terminates.

Proof of Lemma 5.3 (Picking lemma)

To complete the proof of termination of the algorithm, it remains to prove the picking lemma (Lemma 5.3), which is done in this section.

Let us recall briefly the context. Assume that the algorithm needs to refine a simplex s in star S v , with circumball B v (c v (s), r v (s)). The picking lemma states that a valid refinement point can always be found provided that the bound on the sliverity ratio σ 0 is small enough and that the bound on the distortion γ 0 is sufficiently close to 1. The refinement point is searched in the picking region P v (s), a M v -ball with radius δr v (s) centered at the circumcenter c v (s). The refinement point is valid when it does not belong to the so-called forbidden regions. Each forbidden region is associated to a hitting set and consists of the points in the picking region that form with the hitting set either a small sliver or a small well-shaped QC-configuration. Small is here relative to the circumradius r v (s) and controlled by parameter β (see Definition 4.2).

The proof shows that the union of the forbidden regions does not cover the picking region.

In a first step, we show that the volume of each forbidden region is bounded and in fact can be made as small as required with a good choice of the parameters σ 0 and γ 0 (Lemmas 6.5 and 6.7). In a second step, we show that the number of hitting sets, or equivalently of forbidden regions to be avoided, is bounded (Lemma 6.8).

We begin with two technical lemmas. The first one bounds the difference between the two circumspheres of a well-shaped simplex (see Definition 3.1) with respect to two metrics with a bounded distortion.

Lemma 6.1 (Circumsphere lemma) Let M v and M w be two metrics with a distortion γ(M v , M w ) ≤ γ 0 for some γ 0 > 1. Let s be a k-simplex that is well shaped with respect to metric M v . We write c v and r v for the M v -circumcenter and M v -circumradius of s respectively, and we write c w and r w for the M w -circumcenter and M w -circumradius of s, respectively.

• The M v -distance d v (c v , c w ) between the circumcenters satisfies

d v (c v , c w ) ≤ f k (ρ 0 , σ 0 , γ 0 ) r v ( 29
)
where

f k (ρ 0 , σ 0 , γ 0 ) = 1 + 2 k k γ 2 0 ρ k 0 σ k 0 γ 2 0 -1 .
• The circumradius r w is bounded as follows

r w ∈ h - k (ρ 0 , σ 0 , γ 0 )r v , h + k (ρ 0 , σ 0 , γ 0 ) r v where h - k (ρ 0 , σ 0 , γ 0 ) = 1 γ 0 (1 -f k (ρ 0 , σ 0 , γ 0 )) , h + k (ρ 0 , σ 0 , γ 0 ) = γ 0 (1 + f k (ρ 0 , σ 0 , γ 0 )) .
Proof The proof is given in the appendix.

Observe that f k (ρ 0 , σ 0 , γ 0 ) tends to zero when σ 0 tends to 0 and γ 0 tends to 1 in such a way that (γ 0 -1)/σ k 0 tends to 0. We give a name to such functions for further reference.

Definition 6.2 (k-regularly vanishing function)

In the following, a function f of σ 0 and γ 0 is said to be k-regularly vanishing if f tends to zero when σ 0 tends to 0 and γ 0 tends to 1 in such a way that (γ 0 -1)/σ k 0 tends to 0.

The second lemma considers a well-shaped (γ 0 , M )-cospherical configuration U and shows that all the d simplices with vertices in U have nearly the same M -circumradii.

Lemma 6.3 (Circumradii in QC-configurations) Let U be a well shaped (γ 0 , M )-cospherical configuration. The M -circumradii of the d-simplices whose vertices belong to U satisfy:

max s⊂U r M (s) ≤ [1 + η(γ 0 , ρ 0 , σ 0 )] min s⊂U r M (s),
where η(γ 0 , ρ 0 , σ 0 ) is a d-regularly vanishing function.

Note that min s⊂U r M (s) has been defined in Subsection 3.3 as the M -radius r M (U ) of the configuration U .

Proof The proof is given in the appendix.

In the following lemmas, we bound the volume of forbidden regions induced by slivers and QC-configurations by enclosing those regions within spherical shells. Definition 6.4 (Spherical shells) The M v -spherical shell S v (c, r + , r -) with center c and radii r + > r -is the difference between the two M v -balls B v (c, r + ) and B v (c, r -).

The M v -volume of the spherical shell S v (c, r + , r -) is upper bounded by:

V ol v (S v (c, r + , r -) ≤ φ d (r + ) d-1 r + -r -, (30)
where φ d is the measure of the unit (d -1)-sphere and therefore depends only on the dimension d.

Avoiding slivers

Let s be a k-simplex of a star S v . We again write r v for its M v -circumradius. Consider a refinement point p to be taken from the M v -picking region P v (s) of s. Point p is required to lie outside all forbidden regions. We first consider the case of the forbidden region Y v (t) associated to a hitting set t formed by k ≤ d sites such that the k-simplex s = (t, p) is a small k-sliver with respect to a metric M close to M p . More precisely (see Definition 4.2), by a metric M close to M p , we mean a metric M such that γ(M, M p ) ≤ γ 0 , and by a small k-sliver we mean a k-sliver whose M -circumradius is smaller than βr v . Here, as in the rest of the paper, we use the same notation for a subset of sites and the simplex formed by the convex hull of the subset. Note that t is not required to be a simplex appearing in some current star.

Lemma 6.5 (Forbidden regions due to slivers) The M v -volume of the region Y v (t) forbidden by the hitting set t is bounded from above as follows

Vol v (Y v (t)) ≤ µ k (ρ 0 , σ 0 , γ 0 ) β d r d v ,
where µ k (ρ 0 , σ 0 , γ 0 ) is a k-regularly vanishing function.

Proof By the definition of a hitting set, there is a metric M satisfying γ(M p , M ) ≤ γ 0 such that s = (t, p) is a small k-sliver with respect to M . Now, since p belongs to the M v -picking region P v (s) of s, we have γ(M v , M p ) ≤ γ 0 . It follows that γ(M v , M ) ≤ γ 2 0 . Let C(c , r ) and C v (c v , r v ) denotes respectively the M -circumscribing sphere of t and the M v -circumsphere of t. Since s = (t, p) is a small k-sliver with respect to M , t is a wellshaped (k -1)-simplex and its M -circumradius r smaller than βr v . From the sliver lemma (Lemma 3.2), we know that p is at M -distance at most 4πkρ 0 σ 0 r from C(c , r ) ∩ aff(t), where aff(t) is the affine hull of t. Applying the circumsphere lemma (Lemma 6.1) to the well-shaped (k -1)-simplex t, we get

d v (c v , p) ≤ d v (c v , c ) + d v (c , p) ≤ γ 2 0 d M (c v , c ) + γ 2 0 d M (c , p) ≤ γ 2 0 f k-1 (ρ 0 , σ 0 , γ 2 0 )r + γ 2 0 [1 + 4πkρ 0 σ 0 ] r ≤ γ 2 0 1 + f k-1 (ρ 0 , σ 0 , γ 2 0 ) + 4πkρ 0 σ 0 r def = λ + r , writing λ + = γ 2 0 1 + f k-1 (ρ 0 , σ 0 , γ 2
0 ) + 4πkρ 0 σ 0 . In the same way, we have:

d v (c v , p) ≥ d v (c , p) -d v (c v , c ) ≥ 1 γ 2 0 d M (c , p) -γ 2 0 d M (c v , c ) ≥ 1 γ 2 0 [1 -4πkρ 0 σ 0 ] r -γ 2 0 f k-1 (ρ 0 , σ 0 , γ 2 0 ) r def = λ -r , writing λ -= 1 γ 2 0 [1 -4πkρ 0 σ 0 ] -γ 2 0 f k-1 (ρ 0 , σ 0 , γ 2 0 ).
It follows that the forbidden region Y v (s ) associated to s is included in the M v -spherical shell S v (c v , r + v , r - v ) centered at c v and with radii r + v = λ + r and r - v = λ -r . Then, we get from Equation 30 an upperbound for the volume of Y v (s ). The M v -volume of the spherical shell

S v (c v , r + v , r - v ) and thus the M v -volume of Y v (s ) is at most Vol v (Y v (s )) ≤ φ d λ + d-1 (λ + -λ -)r d ,
where

λ + -λ -= γ 2 0 - 1 γ 2 0 + 2γ 2 0 f k-1 (ρ 0 , σ 0 , γ 2 0 ) + γ 2 0 + 1 γ 2 0 4πkρ 0 σ 0 (31)
Since the M -circumradius r of t is smaller than βr v , we get :

Vol v (Y v (s )) ≤ µ k (ρ 0 , σ 0 , γ 0 ) β d r d v , with µ k (ρ 0 , σ 0 , γ 0 ) = φ d λ + d-1 (λ + -λ -)β d r d v .
From the definitions of λ + , Equation 31 and Lemma 6.1 (Circumsphere lemma), the function µ k (ρ 0 , σ 0 , γ 0 ) is k-regularly vanishing.

Avoiding QC-configurations

In Lemma 6.5, we bounded the volume of a forbidden region associated to a sliver. We will now bound the volume of a forbidden region associated to a QC-configuration. We first prove a technical lemma.

Lemma 6.6 (QC-configuration lemma) Given the following:

1. a metric M and a distortion bound γ 0 > 1,

2. a d-simplex s that is well shaped with respect to M . We denote by c and r the M -circumcenter and the M -circumradius of s.

3. a point p such that the configuration (p, s) is a (γ 0 , M )-cospherical configuration.

Then p belongs to the M -spherical shell S M (c, g - d (ρ 0 , σ 0 , γ 0 ) r, g + d (ρ 0 , σ 0 , γ 0 ) r) enclosed between two M -spheres centered at c, with repective radii g - d (ρ 0 , σ 0 , γ 0 )r and g + d (ρ 0 , σ 0 , γ 0 )r, where:

g + d (ρ 0 , σ 0 , γ 0 ) = γ 2 0 + (1 + γ 2 0 )f d (ρ 0 , σ 0 , γ 0 ) g - d (ρ 0 , σ 0 , γ 0 ) = 1 γ 2 0 -(1 + 1 γ 2 0 )f d (ρ 0 , σ 0 , γ 0 ) .
Proof Let N and N be two metrics that witness the (γ 0 , M )-cospherical configuration (s, p), such that p belongs to the interior of the N -circumball B N (s) while p does not belong to the interior of the N -circumball B N (s). Let c N , c N denote respectively the N and N -circumcenters of s. Then, using Lemma 6.1,

d M (p, c) ≤ d M (p, c N ) + d M (c N , c) ≤ γ 0 d N (p, c N ) + f d (ρ 0 , σ 0 , γ 0 )r ≤ γ 0 h + d (ρ 0 , σ 0 , γ 0 )r + f d (ρ 0 , σ 0 , γ 0 )r ≤ γ 2 0 + 1 + γ 2 0 f d (ρ 0 , σ 0 , γ 0 ) r = g + d (ρ 0 , σ 0 , γ 0 ) (32)
and

d M (p, c) ≥ d M (p, c N ) -d M (c N , c) ≥ 1 γ 0 d N (p, c N ) -f d (ρ 0 , σ 0 , γ 0 )r ≥ 1 γ 0 h - d (ρ 0 , σ 0 , γ 0 )r -f d (ρ 0 , σ 0 , γ 0 )r ≥ 1 γ 2 0 -1 + 1 γ 2 0 f d (ρ 0 , σ 0 , γ 0 ) r = g - d (ρ 0 , σ 0 , γ 0 ) (33)
Inequalities ( 32) and ( 33) are just another way to state Lemma 6.6.

Let s be a k-simplex of a star S v , write r v for its M v -circumradius, and consider a refinement point p to be taken from the M v -picking region P v (s) of s. Point p is required to lie outside all forbidden regions. After considering the case of a forbidden region associated to a sliver in the previous section, we consider now the case of a forbidden region W v (t) associated to a hitting set t of d + 1 sites that form with p a small well-shaped M -cospherical configuration for some metric M close to M p . Again (see Definition 4.2), by a metric M close to M p , we mean such that γ(M, M p ) ≤ γ 0 , and by a small configuration, we mean a configuration whose M -circumradius is smaller than βr v . For convenience, as before we denote by t either a subset of sites or the simplex formed by the convex hull of these sites. Note that t is not required to be a simplex appearing in some current star.

Lemma 6.7 (Forbidden regions due to QC-configurations) The M v -volume of the region W v (t) forbidden by the hitting set t is upper bounded as follows

Vol v (W v (t)) ≤ ω(ρ 0 , σ 0 , γ 0 )β d r d v ,
where ω(ρ 0 , σ 0 , γ 0 ) is a d-regularly vanishing function.

Proof As for the proof of Lemma 6.5, we prove that the forbidden region

W v (t) is included in a M v -spherical shell S v (c v , r + v , r - v ) enclosed between two M v -spheres centered at c v , the M v -circumcenter of t.
For the same reason as in the proof of Lemma 6.5, there exists a metric M satisfying γ(M, M v ) ≤ γ 2 0 such that t forms with p a (γ 0 , M )-cospherical configuration. Let c , r be respectively the M -circumcenter and the M -circumradius of t. Applying Lemmas 6.1 and 6.6 to t, we get:

d v (p, c v ) ≤ d v (p, c ) + d v (c , c v ) ≤ γ 2 0 d M (p, c ) + γ 2 0 d M (c , c v ) ≤ γ 2 0 g + d (ρ 0 , σ 0 , γ 0 ) r + γ 2 0 f d (ρ 0 , σ 0 , γ 2 0 ) r def = λ + r , (34)
where

λ + = γ 2 0 g + d (ρ 0 , σ 0 , γ 0 ) + γ 2 0 f d (ρ 0 , σ 0 , γ 2 0 ). Similarly, d v (p, c v ) ≥ d v (p, c ) -d v (c , c v ) ≥ 1 γ 2 0 d M (p, c ) -γ 2 0 d M (c , c v ) ≥ 1 γ 2 0 g - d (ρ 0 , σ 0 , γ 0 ) r -γ 2 0 f d (ρ 0 , σ 0 , γ 2 0 )r def = λ -r , (35)
where

λ -= 1 γ 2 0 g - d (ρ 0 , σ 0 , γ 0 ) -γ 2 0 f d (ρ 0 , σ 0 , γ 2 0 ). It follows that the forbidden region W v (t) associated to t is included in the M v -spherical shell S v (c v , r + v , r - v ) enclosed by the two M v -spheres centered at c v of radii r + v = λ + r and r - v = λ -r .
Therefore, using Equation 30, the volume Vol v (W v (t)) is upper bounded as follows:

Vol

v (W v (t)) ≤ φ d λ + d-1 (λ + -λ -)r d , with λ + -λ -= γ 4 0 - 1 γ 4 0 + γ 2 0 1 + γ 2 0 + 1 γ 2 0 1 + 1 γ 2 0 f d (ρ 0 , σ 0 , γ 0 ) + 2γ 2 0 f d (ρ 0 , σ 0 , γ 2 0 ) (36)
By definition of a hitting set, the QC-configuration (t, p) is required to be small. Specifically, (t, p) has a circumradius that is at most βr v , which, owing to Lemma 6.3, implies that the M -circumradius r of t is at most (1 + η(ρ 0 , σ 0 , γ 0 ))βr v Hence, we can write Vol v (W v (t)) ≤ ω(ρ 0 , σ 0 , γ 0 )

β d r d v with ω(ρ 0 , σ 0 , γ 0 ) = φ d λ + d-1 (λ + -λ -)(1 + η(ρ 0 , σ 0 , γ 0 )) d .
Owing to the definition of λ + (Equation 34), Equation 36, Lemmas 6.1 (Circumsphere lemma) and 6.6 (QC-configuration lemma), function ω(ρ 0 , σ 0 , γ 0 ) is a d-regularly vanishing function.

Bounding the number of forbidden regions Lemma 6.8 Assume that the separation bound sf(p) ≥ Λ sf(p) holds for the current set of vertices and that the algorithm parameters α 0 , β, δ, ρ 0 , σ 0 , and γ 0 satisfy the relation

γ 0 α 0 δ + 2γ 2 0 β (1 + η(ρ 0 , σ 0 , γ 0 )) < 1 ( 37
)
Assume that a refinement point is searched in the M v -picking region P v (s) of the d-simplex s in the star S v , and write K v (s) for the set of hitting subsets of P v (s). The cardinality of K v (s) is bounded by a constant K that depends on α 0 , β, δ, ρ 0 , σ 0 , and γ 0 and remains bounded when σ 0 tends to 0 and γ 0 tends to 1 in such a way that γ 0 -1

σ d 0 tends to 0.
Proof First observe that the cardinality of each hitting subset t in K v (s) is at most d + 1. To bound the cardinality of K v (s), we first bound the cardinality of the set Q v (s) of vertices that may be part of a hitting set t. For this, we use a volume argument based on an upper bound on the distance d v (c v , q) for each q ∈ Q v (s) and a lower bound on the distance d v (q, q ) for any two sites (q, q ) in Q v (s).

Let q be a vertex of Q v (s). The slivers or QC-configurations to avoid are required to have respectively M -circumradii and M -radii smaller than βr v for some metric M such that γ(M, M v ) ≤ γ 2 0 . If q belongs to a hitting set corresponding to a sliver, the M -distance from q to p is at most 2βr v . If q belongs to a hitting set corresponding to a QC-configuration Lemma 6.3, implies that the M -distance from q to p is at most 2β (1 + η(ρ 0 , σ 0 , γ 0 ))r v . In any case, the M v -distance d v (p, q) is therefore at most 2γ 2 0 β (1 + η(ρ 0 , σ 0 , γ 0 ))r v . Moreover, if c v denotes as usual the M v -circumcenter of the simplex s to be refined,

d v (c v , q) ≤ d v (c v , p) + d v (p, q) ≤ δ + 2γ 2 0 β (1 + η(ρ 0 , σ 0 , γ 0 )) r v
We have r v ≤ α 0 sf(c v ) since, when a point is searched in the picking region of a simplex, Rule (1) does not apply anymore. Hence, the inequality above becomes

d v (c v , q) ≤ l 1 sf(c v ), (38)
with

l 1 = α 0 δ + 2γ 2 0 β (1 + η(ρ 0 , σ 0 , γ 0 )) . ( 39
)
We need now to bound the M v -distance d v (q, q ) between two sites in Q v (s) as a function of sf(c v ). Starting from the separation bound hypothesis, we get sd(q) = min q ∈V d q (q, q ) ≥ Λ sf(q).

Then,

min q ∈V d v (q, q ) ≥ Λ Γ sf(q) ≥ Λ Γγ(q, c v ) (sf(c v ) -d cv (c v , q)) (from (8) ) ≥ Λ Γ 2 (sf(c v ) -γ 0 d v (c v , q)) (from (2) ) ≥ Λ Γ 2 (1 -γ 0 l 1 ) sf(c v ) (from (38) ).
Therefore, we have

d v (q, q ) ≥ l 2 sf(c v ) (40) with l 2 = Λ Γ 2 (1 -γ 0 l 1 ) (41)
Observe that l 2 is positive when condition (37) is satisfied. Inequality [START_REF] Tournois | Interleaving Delaunay refinement and optimization for practical isotropic tetrahedron mesh generation[END_REF] shows that the M v -balls centered at the vertices of Q v (s) and with radii l 2 sf(c v )/2 are disjoint and inequality [START_REF] Shewchuk | Star splaying: an algorithm for repairing Delaunay triangulations and convex hulls[END_REF] shows that those balls are contained in the

M v -ball B(c v , (l 1 + l 2 /2) sf(c v )).
A volume argument then proves that the cardinality of

Q v (s) is bounded by (1 + 2l 1 /l 2 ) d .
By considering all possible simplices with vertices in Q v (s), we get a bound on the number |K v (s)| of forbidden regions we need to avoid when picking a refinement point in P v (s)

|K v (s)| ≤ |Q v (s)| d+1 ≤ (1 + 2l 1 /l 2 )) d(d+1) .
Therefore,

|K v (s)| ≤ K (42) with K = (1 + 2l 1 /l 2 ) d(d+1) (43) = 1 + Γ 2 Λ l 1 1 -γ 0 l 1 d(d+1) . ( 44
)
Assume that we choose Λ 2 as in Equation 28. Then Λ is independant of σ 0 and γ 0 . Lemma 6.3 says that η(ρ 0 , σ 0 , γ 0 ) tends to 0 when σ 0 tends to 0 and γ 0 tends to 1 in such a way that γ 0 -1 σ d 0 tends to 0. Therefore l 1 and K remain bounded in the same conditions, which achieves the proof of Lemma 5.3.

Proof of the picking lemma

Proof When a refinement point p has to be picked in the picking region P v (s) of some d-simplex s in star S v , the M v -volume of the picking region P v (s) is δ d r d v (s)u d where u d is the volume of the unit Euclidean ball of dimension d.

To be valid, the refinement point has to lie outside the forbidden regions. In the previous lemmas, we have bounded the M v -volume of the forbidden regions. More precisely, in Lemma 6.5, we gave a bound on the volume of the forbidden region associated to a small ksliver and, in Lemma 6.7, we gave a bound on the volume of the forbidden region associated to a small QC-configuration. Lemma 6.8 bounds the total number of forbidden regions to avoid.

A valid refinement point exists in P v (s) if the volume of the picking region exceeds the total volume of the forbidden regions which is guaranteed if the two following conditions hold:

K µ k (ρ 0 , σ 0 , γ 0 )β d ≤ δ d u d , k = 1, . . . , d (45) K ω(ρ 0 , σ 0 , γ 0 )β d ≤ δ d u d (46)
Assume that α 0 , β, δ, and ρ 0 have been choosen in such a way that equations ( 26), ( 27) and ( 37) are satisfied for any value of γ 0 in [1,[START_REF] Alauzet | Size gradation control of anisotropic meshes[END_REF]. We may for example start with some δ ∈]0, 1[, then choose ρ 0 and β such that equations ( 26) and ( 27) are satisfied for γ = 2. We then choose α 0 so that equation ( 37) is satisfied for γ = 2. Note that these three inequations will remain satisfied for any value of γ 0 in [1,[START_REF] Alauzet | Size gradation control of anisotropic meshes[END_REF]. From Lemma 6.5 and Lemma 6.7, we know that µ k (ρ 0 , σ 0 , γ 0 ) and ω(ρ 0 , σ 0 , γ 0 ) can be made arbitrarily small when σ 0 tends to 0 and γ 0 tends to 1 in such a way that (γ 0 -1)/σ d 0 tends to 0, while lemma 6.8 guarantees that K remains bounded under the same circumstances. It is therefore it is possible to choose σ 0 and γ 0 so as to satisfy Equations ( 45) and (46).

Boundaries and sharp features

Up to this section, we have focused on generating anisotropic meshes that cover a given d-dimensional domain D and conform to a varying anisotropic metric field defined on D. By restricting the stars to the domain D, we have ensured to insert Steiner vertices only within the domain D and we got meshes that roughly cover D. Still, no special attention was paid to get into the final mesh a faithful representation of the domain boundary. This is the purpose of this section. Though the following algorithm could be put at work in any dimension, to be more concrete, we assume in the following that we work in R 3 : the domain is a 3-dimensional domain that is bounded by smooth or piecewise smooth surfaces. The domain may also be subdivided in subdomains by smooth or piecewise smooth surfaces. In the following, we call boundary surface any surface that bounds the domain or one of the subdomains and has to be faithfully represented in the mesh. We denote the domain by D and the set of boundary surfaces by ∂D.

Domains bounded by smooth surfaces

We handle first the case where boundary surfaces are smooth surfaces. In the isotropic setting, the problem of meshing a 3-dimensional domain bounded by smooth surfaces may be solved by a Delaunay refinement algorithm [START_REF] Oudot | Meshing volumes bounded by smooth surfaces[END_REF], based on the notion of restricted Delaunay triangulation. The algorithm refines a set of sites V and its Delaunay triangulation Del(V ), the refinement being guided not only by the restriction of Del(V ) to the domain D but also by its restriction to the set of boundary surfaces ∂D. The restriction of the Delaunay triangulation Del(V ) to a surface is the subcomplex of Del(V ) formed by the facets whose dual Voronoi edges intersect the surface. In the isotropic case, the Delaunay refinement algorithm is known to provide a mesh whose boundary is a a faithfull approximation of the domain surface [START_REF] Oudot | Meshing volumes bounded by smooth surfaces[END_REF]. The algorithm we propose here combines the Delaunay refinement algorithm with the star set system of an anisotropic Delaunay mesh.

The algorithm is summarized in Algorithm 4 below. For each vertex v ∈ V , it maintains two restrictions of the star of v in the Delaunay triangulation Del v (V ) computed using the metric M v of vertex v. The first one is the restricted star S v formed by the tetrahedra of Del v (V ) incident to v and whose M v -circumcenter belongs to the domain D. The second one is the restricted surface star T v formed with the facets of Del v (V ) incident to v and whose M v -dual edges intersect ∂D.

We note S(V ) = {S v , v ∈ V } and T (V ) = {T v , v ∈ V }
those restricted star sets. Facets in T (V ) are also sometimes called surface facets hereafter. The refinement algorithm will insert new Steiner points in V applying refinement rules that aim to get rid of bad facets in T (V ) and bad tetrahedra in S(V ).

A facet in T (V ) is considered bad if either some of its vertices do not belong to a surface in ∂D (topological defect) or if it is oversized, overdistorded, badly shaped or inconsistent. By definition, each facet t in the restricted surface star T v admits an M v -circumball B v (c v (t), r v (t)) centered on a surface in ∂D and empty of vertices of V . Such a ball is called an M v -surface Delaunay ball. The size condition for t is an upper bound on the radius r v (t) and in addition the sizing field used to upper bound r v (t) takes care of the distorsion condition. The shape condition for t is an upper bound on the radius-edge ratio ρ v (t) where ρ v (t) is the ratio from r v (t) to the M v -length of the M v -shortest edge of t. At last, a facet t in T v is considered inconsistent iff it does not belong to the restricted surface star sets of all its vertices.

As in Algorithm 2 above, tetrahedra in S(V ) are considered bad if they are oversized, overdistorded, badly shaped (radius-edge ratio and sliverity conditions) or inconsistent.

Rules are applied with a priority order: Rule (i) is applied only if no Rule (j) with j < i can be applied. Rules for facets have a higher priority than rules for tetrahedra except the rule for inconsistent facets that we postpone at the before last position. Algorithm 4 uses the Pick valid procedure to choose refinement points inserted to get rid of bad surface facest or bad tetrahedra. The refinement point computed for a bad tetrahedron s in star S v is not the M v circumcenter c v (s) of s but a point chosen in the picking region P v (s). The refinement point inserted to get rid of a bad surface facet t in the surface star

T v is not the center c v (t) of a M v -surface Delaunay ball B v (c v (t), r v (t))
of t, but a point chosen in the picking region P v (t) defined as the intersection of the ball B v (c v (t), δr v (t)) with the surface in ∂D including c v (t). Note that the refinement point of a surface facet is always a point of the surface.

Let c be the refinement point computed for a tetrahedron s in some restricted star S v .

Point c is said to encroach a facet t in the restricted surface star T w if it is included in the M w -surface Delaunay ball of t. When a refinement rule is applied to a tetrahedron, the computed refinement point c is first tested for encroachment against the current surface star set and inserted only if no encroachment occurs. Otherwise, the refinement point c is rejected and one of the facets encroached by c is refined instead of the tetrahedron. Algorithm 5 (Insert or snap valid) given below takes care of this behavior. As a result, a refinement point of a tetrahedron is never inserted in the star system if some surface facet is encroached. This ensures that only points in D are inserted as vertices of the star system.

The refinement algorithm uses the constants α 0 , γ 0 , ρ 0 , σ 0 introduced in Algorithm 2, and an additionnal constant α 1 to tune the density of mesh vertices on the domain boundary. The Pick valid procedure still depends on constant β and δ.

At the end of the algorithm, the stars are consistent and the star sets S(V ) and T (V ) can be merged into a consistent mesh. Each simplex in the resulting mesh is well shaped with respect to the metrics of its vertices. The boundary of the mesh is a two-manifold triangulated surfaces whose Hausdorff distance to the domain boundary is controlled. If the sizing field is dense enough, the mesh includes a faithful approximation of all the domain boundary surfaces.

The proof of termination of Algorithm 4 is still based on a volume argument.

First, we notice that the Picking Lemma (Lemma 5.3) is still valid for a refinement point of a surface facet. Indeed, the number of hitting configurations and the volume of forbidden regions can still be bounded as in section 6.

Then, the following theorem whose proof is given in the appendix, provides a lower bound on the insertion radius of each mesh vertex.

Theorem 7.1 Assume that the constant δ is chosen in ]0, 0.5[, that the constants ρ 0 , β, are chosen so as to satisfy: and that α 0 and α 1 are chosen in ]0, 1[ so as to satisfy

ρ 0 ≥ 5 (1 + δ) (1 -δ) 2 Γ 2 γ 8 0 (47) β ≥ 5 1 + δ 1 -δ 2 Γ 2 γ 8 0 , (48)
α 1 ≤ 1 7Γγ 7 0 α 0 1 + α 0 2γ 0 (49) α 1 ≤ 1 8 1 (1 -δ) Γ 2 γ 0 ρ 0 (50) α 1 ≤ 1 8 (1 + δ) (1 -δ) Γ 2 γ 0 β . (51)
Then, there are constants Λ 2 > Λ 4 > Λ 3 such that :

-for each mesh vertex p that is on the boundary surface, the insertion radius r(p) is such that r(p) ≥ Λ 3 sf(p), -for each mesh vertex p that is not a boundary vertex, the insertion radius r(p) is such that r(p) ≥ Λ 2 sf(p) and furthermore the M p -distance δ(p) = d p (p, ∂D) from p to the boundary surface is such that such that : δ(p) ≥ Λ 4 sf(p).

The constants Λ 2 , Λ 3 and Λ 4 depend on the algorithm parameters α 0 , α 1 , γ 0 , ρ 0 , and β and on the global distorsion Γ. Note that conditions (50) and (51) together with conditions (47) and (48) imply that

α 1 ≤ 1 40 1 -δ 1 + δ 1 γ 7 0 . (52)
Equations ( 49) and (52) imply a very dense mesh at least on boundary surfaces. Note however that bounds given in Theorem 7.1 are not tight but largely reflect our will to keep the proof relatively simple.

From the lower bound Λ 3 sf(p) on the insertion radius of each mesh vertex p, we establish a separation bound on mesh vertices as in Lemma 5.8 and conclude the proof of termination by a volume argument as in Section 5.

Domain bounded by piecewise smooth surfaces

In the case of domains bounded by piecewise smooth surfaces, the meshes are also required to include a faithfull representation of the sharp edges (creases) of the bounding surfaces.

A first idea to handle piecewise smooth boundary surfaces is to generalize to sharp edges the notion of restricted Delaunay triangulation and to add to the Delaunay refinement process a refinement level for sharp edges. This has been attempted for the generation of isotropic meshes [START_REF] Rineau | Meshing 3D domains bounded by piecewise smooth surfaces[END_REF] and could be generalized to the star set system. The generation of anisotropic meshes for domains bounded by polyhedral input surfaces is handled this way in [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF]. However this approach has to cope with the problem of small angles subtended by sharp edges and surface patches incident on sharp edges. Indeed input angles smaller than π/2 are known to jeopardize the termination of a Delaunay refinement process. The termination of the refinement process is therefore only granted under severe unrealistic restrictions on the angles formed by boundary surface patches and sharp edges. Such restrictions on input angles are even more stringent in the anisotropic setting where the angular condition on input features has to be respected in the local metric of every point around the feature.

A more promising approach is the method of protecting balls proposed by Cheng et al. [START_REF] Cheng | Theory of a practical Delaunay meshing algorithm for a large class of domains[END_REF][START_REF] Cheng | A practical Delaunay meshing algorithm for a large class of domains*[END_REF]. In this approach, sharp edges are first covered by a set of protecting balls whose centers belong to the sharp edges and define a subdivision of these edges into smaller edges. Protecting balls are considered as weighted points and included as initial points in a weighted Delaunay triangulation. The Delaunay refinement is then performed using this weighted Delaunay triangulation where every additionnal Steiner vertex is inserted with a null weight. Such a weighting scheme ensures the preservation in the final mesh of the initial subdivision of sharp edges induced by the centers of the protecting balls. A solution to generate anisotropic meshes respecting sharp edges would be to transpose the protecting balls approach to the star system of anisotropic Delaunay meshes. Obviously protecting balls in the star system should be turned into balls for the local metric. The possible occurence of metric discontinuities on a sharp edge could be handled by taking for the metric at each point on the sharp edge the intersection of the metrics of both incident patches. The implementation and full study of such an approach will be reported elsewhere.

Conclusion

We have proposed a new class of anisotropic meshes, the so-called anisotropic Delaunay meshes. These meshes conform to a given metric field, can be defined in any dimension, and keep locally the nice properties of Delaunay meshes. We also described an algorithm to generate such meshes in any dimension d. Differently from other methods that have been proposed in dimensions higher than 2, our algorithm produces meshes with a precise characterization and theoretical guarantees.

The algorithm is simple and has been implemented for d = 2 and 3 using the CGAL library [1]. We have also implemented a variant of Algorithm 4 using only Rules (1), ( 3) and [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF] to generate anisotropic surface meshes. Results appear in [START_REF] Boissonnat | Anisotropic Delaunay meshes of surfaces[END_REF]. Figure 5 shows the output of the algorithm on a 3-dimensional ball where the metric is stretched horizontally in the left part and vertically in the right part. The metric field varies slowly on the figure on the left and rapidly on the figure on the right. In this example, we did not enforce any size bound, so that the refinement is only governed by the need to remove inconsistencies. As expected, the mesh density depends on the distortion of the metric. The line where the eigenvectors exchange their eigenvalues is clearly visible on the figure on the right. Further experimental results will be reported elsewhere.

By placing anisotropic meshes in the realm of Delaunay meshes, our framework allows to benefit from recent advances in isotropic mesh generation. In particular, our approach can benefit from local optimization techniques that greatly improve the quality of Delaunay meshes generated by refinement [START_REF] Tournois | Interleaving Delaunay refinement and optimization for practical isotropic tetrahedron mesh generation[END_REF]. For example, since generated meshes are locally Delaunay, ODT methods (optimal Delaunay triangulations) [START_REF] Chen | Optimal anisotropic meshes for minimizing interpolation errors in l p -norm[END_REF][START_REF] Alliez | Variational tetrahedral meshing[END_REF] can be applied in our anisotropic framework.

At last, since our algorithm computes the stars independently and then look for inconsistencies among neighboring stars, it is naturally amenable to parallel computation.

Appendix

9.1 Proof of Lemma 3.2 (Sliver lemma)

Proof (Sliver lemma) In this proof, all lengths, volumes and angles are measured with respect to metric M . We denote by r and r(v) the circumradii of s and s(v) respectively, by V and V (v) their respective volumes, and by e and e(v) the lengths of their respective shortest edges. Let a be the distance from v to the affine hull aff(s(v)) of s(v) and let a be the distance from v to the sphere aff(s(v)) ∩ C(v).

Using the fact that s is a sliver, we have

V = 1 k a V (v) < σ k 0 e k , which yields a < kσ k 0 e k V (v)
.

As s(v) is a face of s, we have e ≤ e(v), and, since s(v) is not a sliver,

V (v) ≥ σ k-1 0 e(v) k-1 . Then, a < kσ k 0 e k σ k-1 0 e(v) k-1 , ≤ kσ 0 e(v) ≤ 2kσ 0 r(v),
which proves the first part of the lemma.

To bound the distance a , we consider the 2-plane through v and the centers c and c of the circumspheres C and C(v) of s and s(v) respectively. See Figure 6. Let p be the projection of v on the affine hull aff(s(v)) and let p be the projection of v on the sphere aff(s(v)) ∩ C(v). Thus a = vp and a = vp . Let q be the point where the ray issued from c that passes through c intersects C. Let ϕ = pp v and θ = qcp . Observe that a a = sin ϕ and sin θ = r(v) r ≥ r(v)/e(v) r/e ≥ 1 2ρ 0 , because r(v) ≥ e(v)/2 and the radius-edge ratio r e of s is smaller than ρ 0 . We distinguish two cases depending on the position of c and v with respect to the affine hull aff(s(v)) of s(v).

In the first case (Figure 6, left part), c and v are on different sides of aff(s(v)). We have ϕ ≥ θ 2 and therefore 9.2 Proof of Lemma 6.1 (Circumsphere lemma)

a = a sin ϕ ≤ a sin( θ 2 ) ≤ 2kσ 0 r(v) sin( 1 2 arcsin 1 2ρ 0 ) ≤ πkσ 0 r(v)
Proof We first prove the circumsphere lemma when s is a d-simplex. The case of a k-simplex, which easily follows,will be considered in subsection 9.2.4.

Computing the circumcenters

Let s = (p 0 , . . . , p d ) be a d-simplex. Since the M v -circumcenter c v of s is at equal M vdistance from all the vertices of s, we have

d 2 v (c v , p i ) = r 2 v for i = 0, . . . , d. Therefore, (p i -c v ) T F T v F v (p i -c v ) = (p 0 -c v ) T F T v F v (p 0 -c v ) i = 1, . . . , d.
Equivalently, we have for i = 0, . . . , d

Inequalities (59), (60), ( 61) and (62) yield (58) which achieves to prove claim 9.1 and inequality (56).

We finally get from (56), ( 57) and (58) :

d v (c v , c w ) ≤ (γ 2 0 -1)r v + 1 2 γ 0 2 γ 2 0 -1 2 d+1 d ρ d 0 σ d 0 r v . ≤ 1 + 2 d d γ 2 0 ρ d 0 σ d 0 γ 2 0 -1 r v
This ends the proof of the first part of Lemma 6.1 in the case of a d-simplex.

Bounding the circumradius r w

Let p be a vertex of s. We have r v = d v (c v , p) and r w = d w (c w , p). Since metric M v satisfies the triangular inequality,

d v (c w , p) -d v (c v , c w ) ≤ d v (c v , p) ≤ d v (c w , p) + d v (c v , c w ).
Then, using the fact γ(M v , M w ) ≤ γ 0 and the first part of Lemma 6.1,

d w (c w , p) γ 0 -f d (ρ 0 , σ 0 , γ 0 )r v ≤ r v ≤ γ 0 d w (c w , p) + f d (ρ 0 , σ 0 , γ 0 )r v r w γ 0 -f d (ρ 0 , σ 0 , γ 0 )r v ≤ r v ≤ γ 0 r w + f d (ρ 0 , σ 0 , γ 0 )r v .
Therefore,

r v γ 0 (1 -f d (ρ 0 , σ 0 , γ 0 )) ≤ r w ≤ r v γ 0 (1 + f d (ρ 0 , σ 0 , γ 0 )) , (63)
which proves the second part of Lemma 6.1 in the case of a d-simplex.

The case of a k-simplex

In the case of a k-simplex s, the circumcenters c v and c w belong to the k-dimensional subspace that is the affine hull, aff(s), of s. If C(v) and C(w) are respectively the M v and M w circumspheres of s, the above proof applies verbatim to the spheres aff(s) ∩ C(v) and aff(s) ∩ C(w) that are the circumspheres of s in the subspace aff(s). This yields the proof of Lemma 6.1 in the case of a k-simplex.

Furthermore, because c encroaches t, we have for any vertex q of t:

r w (t) ≥ 1 2 d w (c, q) ≥ 1 2γ 0 d c (c, q) ≥ 1 2γ 2 0 d v (c, q) ≥ (1 -δ) 2γ 2 0 r v (s) (66)
and, using Equations ( 65) and (66),

r(p) ≥ (1 -δ) 2 2γ 3 0 r v (s).
We complete now the proof of Theorem 7.1 by induction. Let's assume that there exists constants Λ 2 > Λ 4 > Λ 3 such that up to a given stage of the algorithm : -for any mesh vertex q inserted on the boundary surface, the insertion radius r(q) is such that: r(q) ≥ Λ 3 sf(q), -for any mesh vertex q that is not on the boundary surface, the insertion radius r(q) is such that r(q) ≥ Λ 2 sf(q) and the M q -distance δ(q) = d q (q, ∂D) from q to the boundary surface satisfies : δ(q) ≥ Λ 4 sf(q). Performing a case analysis on the rule that triggers the insertion of the next vertex p, we compute a lower bound on the insertion radius r(p) and a lower bound on the distance δ(p) if p does not belong to the boundary surface.

Lower bound on the insertion radius r(p)

Rule (1). Assume p is inserted by Rule (1) applied on facet t of T v . Then we have p = Pick valid(t, M v ) and, using Lemma 9.2, r(p)

≥ (1-δ) Γ r v (t). Then, r v (t) ≥ α 1 sf(c v (t)) ≥ α 1 Γ [sf(p) -d p (p, c v (t))] ≥ α 1 Γ [sf(p) -Γd v (p, c v (t))] ≥ α 1 Γ [sf(p) -Γδr v (t)] (67)
Hence,

r v (t) ≥ α 1 Γ(1 + δα 1 ) sf(p) and r(p) ≥ (1 -δ) Γ 2 α 1 (1 + δα 1 ) sf(p). ( 68
)
The induction hypothesis is fulfilled for p, if we have:

Λ 3 ≤ (1 -δ) Γ 2 α 1 (1 + δα 1 ) . ( 69
)
Rule [START_REF] Alauzet | Size gradation control of anisotropic meshes[END_REF]. Assume now that Rule 2 is applied. Then p = Pick valid(t, M v ) where t is a facet in T v with a vertex q that does not belong to ∂D. From the induction hypothesis δ(q) ≥ Λ 4 sf(q). We have:

d v (p, q) ≤ (1 + δ)r v (t)
and

r v (t) ≥ 1 (1 + δ) d v (p, q) ≥ 1 γ 0 (1 + δ) d q (p, q) ≥ 1 γ 0 (1 + δ) δ(q) ≥ 1 γ 0 (1 + δ) Λ 4 sf(q), ≥ Λ 4 γ 2 0 (1 + δ) [sf(p) -d p (p, q)] ≥ Λ 4 γ 2 0 (1 + δ) [sf(p) -γ 0 (1 + δ)r v (t)] r v (t) ≥ Λ 4 γ 2 0 (1+δ) 1 + Λ 4 γ 0 sf(p)
Then, using Lemma 9.2:

r(p) ≥ 1 -δ 1 + δ Λ 4 γ 2 0 (γ 0 + Λ 4 ) sf(p).
The induction hypothesis is fulfilled for p, if we have:

Λ 3 ≤ 1 -δ 1 + δ Λ 4 γ 2 0 (γ 0 + Λ 4 ) . ( 70
)
Rule (3). Assume p is inserted by Rule (3) applied on facet t of T v . Then, p = Pick valid(t, M v ) and, from Lemma 9.2, we have r(p) ≥ (1-δ) γ 0 r v (t) To get a lower bound for r v (t), we argue as in the proof of Lemma 5.6 and, replacing Λ 2 by Λ 3 in Equation ( 17)), we get: where the last equation makes use of the fact that, since p = Pick valid(t, M w ), r(p) ≥

r v (t) ≥
(1-δ) γ 0 r w (t). Then, we get:

r v (s) ≥ α 0 Γ sf(p) -Γγ 0 (1 + δ) (1 -δ) r(p) -Γδr v (s) r v (s) ≥ α 0 Γ(1 + α 0 δ) sf(p) -Γγ 0 (1 + δ) (1 -δ) r(p) (77)
from which we deduce:

r(p) ≥ (1 -δ) 2 2Γ 2 γ 0 α 0 1 + α 0 δ sf(p) -Γγ 0 (1 + δ) (1 -δ) r(p) r(p) ≥ (1 -δ) 2 2Γ 2 γ 0 α 0 1 + α 0 δ 1 1 + (1-δ 2 ) 2Γγ 0 α 0 (1+α 0 δ) sf(p). ( 78
)
It follows that the inductive hypothesis on the insertion radii is still fulfilled in this case if (87)

Λ 3 ≤ (1 -δ) 2 2Γ 2 γ 0 α 0 (1 + α 0 δ) 1 1 + (1-δ 2 )
The inductive hypothesis on the insertion radii is fulfilled in this case provided that Rule ( 6) and [START_REF] Boissonnat | Anisotropic Delaunay meshes of surfaces[END_REF]. Assume p is inserted by Rule [START_REF] Boissonnat | Anisotropic diagrams: Labelle and Shewchuk approach revisited[END_REF] or [START_REF] Boissonnat | Anisotropic Delaunay meshes of surfaces[END_REF] applied to a tetrahedron s of S v .

No snapping

If no encroachement occurs, p = Pick valid(s, M v ) and from Lemma 9.2 r(p) ≥ (1-δ) γ 0 r v (s). Let q be the last inserted vertex of simplex s. Vertex q has been inserted as the refinement point of a simplex s in some star S w and we have r v (s) ≥ βr w (s ). Then we argue as in the proof of Lemma 5.

Snapping

If encroachment occurs, p = Pick valid(t, M w ) for some facet t in T w encroached by p and, according to Lemma 9.2, we have r(p) ≥ (1-δ) 2 2γ 3 0 r v (s). Let q be the last inserted vertex of q. We argue as in the proof of Lemma 5.7:

r v (s) ≥ βr w (s ) ≥ β γ 0 (1 + δ) r (q)
≥ βΛ 3 γ 0 (1 + δ) sf(q), which is r v (s) ≥ A sf(q) with A = βΛ 3 γ 0 (1 + δ) .

Repeating the calculation leading from Equation (82) to Equation (87), we get:

r(p) ≥ (1 -δ) 2 2γ 5 0 A (1 + (1 + δ)A) 1 1 + (1-δ 2 ) 2γ 3 0 A (1+(1+δ)A) sf(p). (91)
The inductive hypothesis on the insertion radii is fulfilled in this case provided that

Let q be the last inserted vertex of t. Vertex q is the refinement point of a simplex t and we have:

r v (t) ≥ βr w (t ) ≥ β γ 0 (1 + δ) r(q). r v (t) ≥ βΛ 3 γ 0 (1 + δ) sf(q). ( 94
)
Observe that we can drop Equations ( 72) and (96) because they are implied respectively by Equations ( 88) and (92). From Equation (69), we know that we will have Λ 3 ≤ (1δ) α 1 Γ 2 , In fact, for technical reasons, we will choose:

Λ 3 = 0.8 (1 -δ) (1 + δ) α 1 Γ 2 . ( 104
)
Then from hypothesis in Equations ( 50) and (51) we will have

(1 + δ)ρ 0 Λ 3 γ 0 ≤ 0.1 ( 105
)
βΛ 3 γ 0 ≤ 0.1, (106)
which allows to simplify Equations ( 81), ( 88), ( 90), (92) leading to the system:

Λ 3 ≤ (1 -δ) Γ 2 α 1 (1 + δα 1 ) (69) Λ 3 ≤ 1 -δ 1 + δ Λ 4 γ 2 0 (γ 0 + Λ 4 ) (70)
Λ 3 ≤ γ 0 ρ 0 (1 + δ) ρ 0 (1 -δ) γ 3 0 -1 (73)
Λ 2 ≥ 2γ 5 0 α 1 (1 -γ 2 0 α 1 ) (102)
Λ 4 ≤ 1 1 + 2γ 5 0 α 1 1-γ 2 0 α 1 Λ 2 - 2γ 5 0 α 1 1 -γ 2 0 α 1 (103)
Let us then choose Λ 2 from Equations ( 75) and (102s). From Equation (52), we have γ 2 0 α 1 ≤ γ 7 0 α 1 ≤ 0.1, so that 2γ 5 0 α 1 (1γ 2 0 α 1 ) ≤ 2γ 5 0 α 1 0.9 ≤ 2.23γ 5 0 α 1 (107)

Then Equation (102s) is satisfied if we choose:

Λ 2 = 3.5γ 5 0 α 1 (108)
Condition (49) ensures that this choice satisfies Equation (75). We then choose Λ 4 from Equation (103). Let

A = 1 1 + 2γ 5 0 α 1 1-γ 2 0 α 1 Λ 2 - 2γ 5 0 α 1 1 -γ 2 0 α 1
Since from Equation (107), 2γ 5 0 α 1 (1-γ 2 0 α 1 ) ≤ 2.3 γ 5 0 α 1 and from Equation (52), γ 5 0 α 1 ≤ γ 7 0 α 1 ≤ 0.1, we have:

A ≥ 1 1 + 2.3 γ 5 0 α 1 Λ 2 - 2γ 5 0 α 1 1 -γ 2 0 α 1 ≥ 1 1.3
1.2 γ 5 0 α 1 ≥ 0.9 γ 5 0 α 1 , so that Equation ( 103) is satisfied if we choose Λ 4 = 0.9 γ 5 0 α 1

.

It remains to check that Equations (70), (73), (79), (81s), (88s), (90s) and (92s) are satisfied.

Equation ( 70) is satisfied since 1δ 1 + δ Λ 4 γ 2 0 (γ 0 + Λ 4 ) = 1δ 1 + δ 0.9 γ 5 0 α 1 γ 2 0 (γ 0 + 0.9 γ 5 0 α 1 ) ≥ 1δ 1 + δ 0.9 γ 2 0 α 1 (1 + 0.9 γ 4 0 α 1 )

≥ 1 -δ 1 + δ 0.8 γ 2 0 α 1 ≥ 1 -δ 1 + δ 0.8 α 1 Γ 2 = Λ 3
Since Equation (105) is granted, Equation (73) is satisfied if

ρ 0 (1 -δ) γ 3 0 -1 ≥ 0.1 ⇐⇒ ρ ≥ 1.1 γ 3 0 (1 -δ) ,
which is implied by Condition (47).

Equation ( 79) is equivalent to:

0.8 (1 -δ) (1 + δ) α 1 Γ 2 ≤ (1 -δ) 2 4Γ 2 γ 0 α 0 (1 + α 0 δ) ⇐⇒ α 1 ≤ (1 -δ 2 ) 3.2 γ 0 α 0 (1 + α 0 δ) Since (1-δ 2 )
3.2 γ 0 α 0 (1+α 0 δ) ≥ 0.75 3.2 γ 0 α 0 1.5 ≥ α 0 7γ 0 , this is granted by Condition (49). Equation (81s) is equivalent to:

ρ 0 ≥ 1 0.9(1 -δ) γ 3 0 Λ 2 Λ 3
From our choice for Λ 2 and Λ 3 ( Equation ( 104) and ( 108) ),

Λ

Figure 1 :

 1 Figure 1: Example of inconsistent stars in 2D: stars S v and S w are inconsistent because edge [vw] belongs to S v but not to S w .

Lemma 3 . 2 (

 32 Sliver lemma) Let s be a k-simplex and M a metric. If v is a vertex of s, we denote by s(v) the (k -1)-face of s opposite to vertex v, by aff(s(v)) the affine hull of s(v), i.e. the (k -1)-flat spanned by the vertices of s(v), by C M (s(v)) the M -circumsphere of s(v), and by r M (s(v)) the M -radius of C M (s(v)).

Figure 2 :

 2 Figure 2: (s = abc, d) is a QC-configuration because d is inside B N (s) but outside B N (s)

Figure 4 :

 4 Figure 4: For the proof of Lemma 5.7.

Figure 6 :

 6 Figure 6: For the proof of the sliver lemma.

  p (p, c v (s)) ≤ d p (p, c) + d p (c, c v (s)) ≤ Γd w (p, c) + Γd v (c, c v (s)) ≤ Γ(1 + δ)r w (t) + Γδr v (s) ≤ Γγ 0 (1 + δ) (1δ) r(p) + Γδr v (s)

.( 1 -Rule 5 . 2 0( 1

 1521 Because δ and α 1 belongs to [0, 1] while Γ and γ 0 are greater than 1, Assume p is inserted by Rule (5) applied to the tetrahedron s of S v . andd p (p, q) ≤ d p (p, c) + d p (c, q) ≤ γ 0 d w (p, c) + γ 0 d c (c, q) ≤ γ 0 ((1 + δ)r w (t) + γ 2 0 d v (c, q) ≤ γ 0 (1 + δ)r w (t) + γ 2 0 (1 + δ)r v (s) δ)r v (s) r v (s) ≥ A γ 2 0 (1 + (1 + δ)A) sf(p)γ

( 1 -

 1

2

 2 7 and, replacing Λ 2 by Λ 3 in Equation (25)) induction hypothesis on the insertion radii is still verified in this case if Λ

( 1 -Rule ( 7 )

 17 Assume p is inserted by Rule[START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF] applied to the facet t of T v . We have p = Pick valid(t, M v ) from Lemma 9.2,r(p) ≥ (1δ) γ 0 r v (t).

  where . is the matrix norm operator associated with the Euclidean norm, i.e. for a d × d square matrix A, A = sup x∈R d

  Algorithm 4 Refinement algorithm for domain with smooth boundary Facet radius-edge ratio If a facet t in star T v is such that ρ v (t) > ρ 0 , Insert(Pick valid(t, M v )); Tet size and distorsion If a tetrahedron s in some star S v , is such that r v (s) ≥ α 0 sf(c v (s)), Insert or snap valid(s, M v ); Tet radius-edge ratio: If a tetrahedron s in some star S v is such that ρ v (s) > ρ 0 , Insert or snap valid(s, M v );

	Rule (1) Facet size and distortion
	If there is a facet t in star T v
	such that r v (t) ≥ α 1 sf(c v (t)), Insert(Pick valid(t, M v ));
	Rule (2) Facet without topological defect
	If there is a facet t in star T v
	with some vertex ∈ ∂D Insert(Pick valid(t, M v ));
	Rule (3) Rule (4) Rule (5) Rule (6) Sliver removal:
	If a tetrahedron s in star S v is a M v -sliver
	Insert or snap valid(s, M v );
	Rule (7) Facet consistency:
	If a facet t in some star T v is inconsistent
	Insert(Pick valid(t, M v ));
	Rule (8) Tetrahedron consistency
	If a tetrahedron s in some star S v is inconsistent
	Insert or snap valid(s, M v );
	Algorithm 5 Insert or snap valid(s, M v )
	c = Pick valid(s, M v )
	If c encroaches some facet t in some restricted surface star T w
	insert(Pick valid(t, M w ))
	else insert(c)

Acknowledgments

The authors would like to thank Tamal Dey for fruitful discussions. This work has been partially supported by the Agence Nationale de la Recherche (ANR) under project GIGA (Geometric Inference and Geometric Approximation). The research leading to these results has also received funding from the European Research Council (ERC) under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement No. 339025 GUDHI (Algorithmic Foundations of Geometry Understanding in Higher Dimensions) Figure 5: Two examples of anisotropic meshes produced by our algorithm.

Writing P = (p 1p 0 , . . . , p dp 0 ) for the square matrix whose columns are the vectors p ip 0 , i = 1, . . . , d, and Diag(A) for the column matrix whose elements are the elements of the main diagonal of a square matrix A, the last equation becomes

, from which we get the position of c v with respect to the position of the vertices of s

An equivalent formula gives the M w -circumcenter c w of s.

9.2.2

Bounding the distance between c v and c w

In the following, we choose a coordinate system in which M v = F T v F v and F v are identity matrices. (Equivalently, we could assume without loss of generality that M v is the Euclidean metric since the distance d v (c v , c w ) is only related to the relative distortion between M v and M w .) Then, we deduce from (53) :

where q = Diag(P T P ), q = Diag(P T F T w F w P ). We further write

where I is the identity matrix. By our choice of the coordinate system, the M v -norm of a vector x is just the Euclidean norm x of its coordinates in this reference system. Therefore,

The following claim provides bounds for the two terms on the right hand side of (56).

Claim 9.1

Proof Writing A for the Euclidean norm of a matrix A, ( A = sup x =1 Ax ), we have,

. The absolute values of the eigenvalues of matrix I -(F T w F w ) -1 are thus at most γ 2 0 -1. Moreover, from (54), P -T q = 2d v (c v , p 0 ) is just twice the M v -circumradius of s, which proves inequality (57).

To prove (58), we write

We will bound the three terms on the right hand side of (59). We first note that

Then, for P -T , we use the fact that P -T ≤ P -T ∞ where P -T ∞ is the maximum absolute value of any entry in P -T . Each entry in P -T is a cofactor of matrix P T divided by the determinant of P T . The determinant of P T is d! times the M v -volume of s. Each entry in P T is a coordinate of some p ip 0 and therefore less than p ip 0 ≤ 2r v , which implies that each cofactor of P T is at most (d -1)!(2r v ) d-1 . Therefore,

where e v is the M v -length of the shortest (for M v ) edge of s. We now bound qq :

9.3 Proof of Lemma 6.3 (Circumradii in QC-configurations lemma)

Proof Let s min and s max be the simplices with vertices in U having respectively the minimum and maximum M -circumradius.

Let N and N be the two metrics witnessing the quasi-cosphericity of U . We consider the set of metrics with distortion less than γ 0 from M and a continuous path joining N to N within this set, for instance the linear interpolation between N and N . Since the Delaunay triangulations Del N (U ) and Del N (U ) are different and since the metric evolve continously along the path, there is at least a metric M on the path, with γ(M, M ) ≤ γ 0 and such that U is M -cospherical which means that all d-simplices with vertices in U have the same M -circumradius.

Then, applying twice the Circumsphere lemma 6.1 respectively to s min and s max , we get:

which proves lemma 6.3, setting

Proof of Theorem 7.1

Proof We begin by a lemma relating the insertion radius of each mesh vertex to the radius of the bad simplex triggering the insertion. This generalizes Lemma 5.4 by taking care of the surface facet refinement rules and of the effect of snapping mesh vertices to the surface when some encroachment occurs (see the Insert or snap valid procedure).

Let p be a mesh vertex. The vertex p is inserted by application of one of the refinement rules 1-8 to either a surface facet t of some star T v whose M v -surface Delaunay ball radius is denoted by r v (t), or to a tetrahedron s of some star S v whose M v -circumradius is denoted by r v (s).

Lemma 9.2 (Second Insertion radius lemma) The insertion radius r(p) of the mesh vertex p is such that: [START_REF] Alauzet | Multi-dimensional continuous metric for mesh adaptation[END_REF] or Rule [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF] applies.

• When Rule (4) applies, r(p) ≥ (1-δ) Γ r v (s) if no encroachment occurs and r(p) ≥

(1-δ) 2 2Γγ 0 r v (s) otherwise.

• When one of Rule ( 5),( 6) or (8) applies, r(p)

Proof If Rule(1) applies, p = Pick valid(t, M v ), and we have

The cases of Rule ( 2), ( 3) and ( 7) is analogous to the case of Rule (1) except that now the distorsion γ(p, v) is known to be less than γ 0 .

When no snapping occurs, the cases of Rule ( 4), ( 5), ( 6) and ( 8) are analogous to the cases of Rule ( 2), ( 3) and [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF].

Assume that Rule (4) applies to a tetrahedron s and that snapping occurs. Then the point c output by Pick valid(s, M v ) encroaches a facet t in some surface star T w , and p = Pick valid(t, M w ). We have:

Furthermore, because c encroaches facet t in surface star T w , we have, for any q that is a vertex of t,

and therefore

Assume now that one of Rule ( 5), ( 6) or ( 8) applies, and that snapping occurs. As above c = Pick valid(s, M v ) encroaches a facet t in some surface star T w , and p = Pick valid(t, M w ). We have:

which implies

and

Rule ( 4) Assume p is inserted by Rule (4) applied to the tetrahedron s of the star S v .

No snapping

If no snapping occurs, p = Pick valid(s, M v ), and according to lemma 9.2: r(p) ≥

(1-δ)

Γ r v (s). Then,

and we perform a computation analogous to the computation leading to Equation (67) to get:

Then

and the induction hypothesis is fulfilled for p, if we have:

Snapping

If snapping occurs, p = Pick valid(t, M w ) where t is a facet of some star T w that is encroached by the computed refinement point c = Pick valid(s, M v ).

From Lemma 9.2, we have r(p)

No snapping If no encroachment occurs, we have p = Pick valid(s, M v ) and from Lemma 9.2 r(p) ≥ (1-δ) γ 0 r v (s). We argue as in the proof of Lemma 5.6 and, replacing Λ 2 by Λ 3 in Equation ( 17)), we get:

The induction hypothesis on the insertion radii is satisfied in that case provided that:

Snapping

If encroachment occurs, p = Pick valid(t, M w ) for some facet t in T w encroached by c = Pick valid(s, M v ) and, according to Lemma 9.2, we have r(p)

where e v (s) is the M v -length of the M v -shortest edge of s and q is the last inserted vertex of e v (s). Then,

For future reference, we set r v (s) ≥ A sf(q) (82)

Then we have:

Furthermore,

Then, since p = Pick valid(t, M v ), we have

Combining Equation ( 94) and (95) we get:

and, from Equation (93),

Therefore, the inductive hypothesis on the insertion radii is still satisfied in this case if:

Lower bound on δ(p).

It remains to establish a lower bound on δ(p) in case one of Rules (4), ( 5), ( 6) or ( 8) is applied and no snapping on the surface occurs. Let p be the vertex inserted and x be the point of ∂D closest to p according to the metric M p :

At the time p is inserted, surface Delaunay balls of surface facets cover ∂D and point x belongs to the Delaunay surface ball B w (c w (t), r w (t)) of some facet t in the star S w of some vertex w. Therefore,

Then,

57 Furthermore, we may assume that γ(x, p) ≤ γ 0 . Indeed, otherwise δ(p) = d p (p, x) ≥ bdr 0 and we are done. Therefore γ(p, w) ≤ γ(p, x)γ(x, w) ≤ γ 2 0 . Let us now consider δ(p) = d p (p, x). We have:

≥ r(p)γ 2 0 d w (w, x) ≤ r(p) -2γ 2 0 r w (t)

≥ r(p) -

where Equation (99) makes use of Equation (98).

To get a lower bound for δ(p), we now have to get an upper bound for sf(x). We have: and using the induction hypothesis for the insertion radii of internal mesh vertices we get:

so that the inductive hypothesis on distance to the surface is satisfied if

and

Close up

The inductive hypothesis is fulfilled if we can satisfy Equations 69 , 70, 72, 73, 75, 79, 81, 88, 90, 92, 96, 102 and 103.