Keywords: Unstructured mesh, continuous mesh, Riemannian metric space, interpolation error, linear interpolate Maillage non structuré, maillage continu, espaces métriques riemanniens, erreur d'interpolation, interpolé linéaire

In the context of mesh adaptation, Riemannian metric spaces have been used to prescribe orientation, density and stretching of anisotropic meshes. Such structures are used to compute lengths in adaptive mesh generators. In this report, a Riemannian metric space is shown to be more than a way to compute a distance. It is proven to be a reliable continuous mesh model. In particular, we demonstrate that the linear interpolation error can be derived continuously for a continuous mesh.

In its tangent space, a Riemannian metric space reduces to a constant metric tensor so that it simply spans a metric space. Metric tensors are then used to continuously model discrete elements. On this basis, geometric invariants have been extracted. They connect a metric tensor to the set of all the discrete elements which can be represented by this metric. As the behavior of a Riemannian metric space is obtained by patching together the behavior of each of its tangent spaces, the global mesh model arises from gathering together continuous element models. We complete the continuous-discrete analogy by providing a continuous interpolation error estimate and a well-posed definition of the continuous linear interpolate. The later is based on an exact relation connecting the discrete error to the continuous one.

From one hand, this new continuous framework freed the analysis of the topological mesh constraints. On the other hand, powerful mathematical tools are available and well defined on the space of continuous meshes: calculus of variations, differentiation, optimization, . . . , whereas these tools are not defined on the space of discrete meshes.

Introduction

In this article, a continuous mesh concept is introduced using the notion of Riemannian metric space of the differential geometry. The notion of continuous mesh aims at modeling unstructured discrete computational meshes classically used in numerical simulations solved by the finite element, the finite volume, the discontinuous Galerkin,... numerical methods. The final purpose of this model is twice. From a theoretical point of view, it enables us to prove the consistency of classical metric-based mesh adaptation. To this end, it provides exact relations between the linear interpolation error and the mesh prescription. From a practical point of view, the objective is to derive a continuous interpolation error estimate that can be used to automatically perform anisotropic mesh adaptation. The proposed continuous formalism permits a complete abstraction of the notion of mesh. In particular the geometric data (e.g. the vertex coordinates) and the topology description (e.g. the mesh entities) do not exist anymore. As regards interpolation error estimates in the continuous mesh framework, they are free of any a priori hypothesis on the mesh such as alignment or density requirements. As a consequence, this report aims at proving that Riemannian metric space is more than a way to define distance computation within anisotropic mesh generators. It is a reliable mesh model that is well suited for metric-based mesh adaptation.

Fundamentals of anisotropic mesh adaptation. When dealing with mesh adaptation, the prescription and the generation of adapted meshes are crucial issues. There exists a large class of methods to prescribe and to generate adapted meshes depending on the problem at hand along with the mesh specificity: uniform, isotropic, anisotropic, . . . The simplest algorithms consist in refining or coarsening the current mesh according to patterns. However, such strategies encounter several bottlenecks. In particular, mesh coarsening can only be applied to regions already refined by patterns, i.e., only added patterns can be removed, and they do not allow anisotropic meshes to be generated. Moreover, the quality of the sequence of refined meshes is strongly related to the quality of the initial mesh and this quality can only decrease during the refinement process. A generic and elegant way to generate anisotropic meshes is to use the notion of metric and Riemannian metric space. An adapted anisotropic mesh in this framework is simply the image in the Euclidean space of a uniform mesh in a Riemannian metric space.

To this end, the distance in the adaptive mesh generator is computed in the Riemannian metric space instead of the Euclidean one. This technique generalizes the case of uniform meshes where distances are computed in the classical Euclidean space. Then, the main idea to generate an anisotropic adapted mesh is to complete a unit mesh in the given Riemannian metric space. It consists in generating a mesh where the edges length are equal to one with respect to the prescribed Riemannian metric space. This approach was initiated by Hecht, see developments in [START_REF] Hecht | Mesh adaptation by metric control for multi-scale phenomena and turbulence[END_REF] and in Vallet's thesis [START_REF] Vallet | Génération de maillages éléments finis anisotropes et adaptatifs[END_REF]. This method is commonly called metric-based mesh adaptation. There are actually a lot of available codes that are based on the metric concept. Let us mention Bamg [START_REF] Hecht | BAMG: bidimensional anisotropic mesh generator[END_REF] and BL2D [START_REF] Laug | BL2D-V2, mailleur bidimensionnel adaptatif[END_REF] in 2D, Yams [START_REF] Frey | Yams, a fully automatic adaptive isotropic surface remeshing procedure[END_REF] for discrete surface mesh adaptation and Forge3d [START_REF] Coupez | Génération de maillages et adaptation de maillage par optimisation locale[END_REF], Fun3d [START_REF] Jones | Validation of 3D adjoint based error estimation and mesh adaptation for sonic boom prediction[END_REF], Gamic [START_REF] George | Gamhic3d, a fully automatic adaptive mesh generation method in three dimensions[END_REF], MeshAdap [START_REF] Li | 3D anisotropic mesh adaptation by mesh modification[END_REF],

Mmg3d [START_REF] Dobrzynski | Anisotropic delaunay mesh adaptation for unsteady simulations[END_REF], Mom3d [START_REF] Tam | Anisotropic mesh adaptation for 3D flows on structured and unstructured grids[END_REF], Tango [START_REF] Bottasso | Anisotropic mesh adaption by metric-driven optimization[END_REF] and [START_REF] Pain | Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations[END_REF] in 3D. It is worth mentioning that all these codes have arisen from different mesh generation methods. The method in [START_REF] George | Gamhic3d, a fully automatic adaptive mesh generation method in three dimensions[END_REF][START_REF] Hecht | BAMG: bidimensional anisotropic mesh generator[END_REF] is based on a global constrained Delaunay kernel. In [START_REF] Laug | BL2D-V2, mailleur bidimensionnel adaptatif[END_REF], the Delaunay method and the frontal approaches are coupled. [START_REF] Frey | Yams, a fully automatic adaptive isotropic surface remeshing procedure[END_REF] is based on local mesh modifications. And, [START_REF] Coupez | Génération de maillages et adaptation de maillage par optimisation locale[END_REF] is based on the minimal volume principle. Another benefit is that this mesh prescription approach is naturally derived from anisotropic interpolation error estimate. Indeed, interpolation error involved a Hessian matrix from which a metric is easily derived. If the use of a Hessian matrix to defined a metric tensor is now classical for generating anisotropic meshes, it remains to evaluate the impact of the practical algorithm used to generate the mesh. In other words, is it a relevant choice to use metric fields for the control of the interpolation error ? A lot of numerical examples for real life problems [START_REF] Bottasso | Anisotropic mesh adaption by metric-driven optimization[END_REF][START_REF] Dompierre | Anisotropic mesh adaptation: towards a solver and user independent cfd[END_REF][START_REF] Formaggia | Anisotropic mesh adaptation in computational fluid dynamics: Application to the advection-diffusion-reaction and the Stokes problems[END_REF][START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF][START_REF] Loseille | Achievement of second order mesh convergence for discontinuous flows with adapted unstructured mesh adaptation[END_REF][START_REF] Pain | Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations[END_REF][START_REF] Schall | Mesh adaptation as a tool for certified computational aerodynamics[END_REF][START_REF] Tam | Anisotropic mesh adaptation for 3D flows on structured and unstructured grids[END_REF]] tend to answer affirmatively to this question. In this work, this question is theoretically studied by using a continuous mesh model.

The proposed continuous framework. In this paper, a fully continuous mesh model is introduced based on the notion of Riemannian metric space. The continuous mesh is a function that prescribes at each point a density, anisotropic quotients and orientations. This model is completely generic and geometric. Several geometric relations connect the continuous mesh to the set of its discrete representatives. The main pro of this approach is to be independent of the mesh generation algorithm used to generate adapted meshes.

This model is then used to study the interpolation error. Some mathematical developments lead to a single relation that connects an infinite set of discrete elements to a unique continuous estimate. From this estimate, the interpolation error for a given continuous mesh can be accurately predicted whatever the considered smooth function. Contrary to classical approaches inherited from [START_REF] Castro-Díaz | Anisotropic unstructured mesh adaptation for flow simulations[END_REF], there are optimal conditions leading to alignment requirements between the mesh and the Hessian function. The derivation of a general interpolation error estimate has already been studied for a single element in the case of a quadratic function [START_REF] Cao | On the error of linear interpolation and the orientation, aspect ratio, and internal angles of a triangle[END_REF]. These studies enable the best shaped element for a given norm to be derived. However, it is, from a practical point of view, impossible to generate a mesh only composed of optimal elements. Moreover, the variation of the function needs to be taken into account when the function is no more quadratic. In this work, estimates are derived locally both for the continuous mesh and the functions thanks to Taylor expansion. If the Taylor expansion of a function is common, the Taylor expansion of a continuous mesh consists in working in the tangent space of the Riemannian metric space where the metric tensor is constant. Then both the variations of the mesh and the function are integrated on the whole computational domain. As a consequence, estimates are still very accurate for non quadratic functions and non uniform continuous meshes. This point is demonstrated in the numerical examples. Theoretically, the use of a discrete support is no more mandatory to compute the interpolation error. Several fully continuous and analytical examples are given to exemplify this feature. More practically, only a background mesh that supports a discrete representation of the continuous mesh is needed when dealing with real-life appli-INRIA cations. There is no need to generate a discrete mesh that corresponds to the continuous one.

Overview. In a first section, we recall the differential geometry notions that are recurrently used in the continuous approach. Next, the continuous element and mesh models are introduced. Geometrical invariants, that connect continuous elements and discrete elements, are proved. In a third section, these results are used to exactly predict the linear interpolation error for discrete and continuous elements. Then, the notion of continuous linear interpolation is well established and exactly estimated. Finally, numerical experiments emphasize the possibility to compute, for a given function and a given continuous mesh, the continuous interpolation error without any discrete supports. The correlation between continuous and discrete estimations of the interpolation error is finally shown.

Metric notion for mesh adaptation

In order to be self-contained, we recall the differential geometry notions that are used in the sequel. It mainly concerns the computation of lengths for different kinds of metric spaces: Euclidean, Euclidean metric space and Riemannian metric space. A complete review and the mathematical study of these spaces are available in [START_REF] Berger | Geometry 1 & 2[END_REF][START_REF] Berger | A panoramic view of Riemannian geometry[END_REF][START_REF] Carmo | Differential geometry of curves and surfaces[END_REF].

Notations. Bolded symbols, as a, b, u, v, x, . . ., denote vectors or points of R n . Vector coordinates are denoted by x = (x 1 , . . . , x n ). And, the natural dot product between two vectors u and v of R n is:

u, v = n i=1 u i v i .

Euclidean metric space

An Euclidean metric space (R n , M) of dimension n is a finite vector space where the dot product is defined by means of a symmetric definite positive form M :

u, v M = u, Mv = t uMv , for (u, v) ∈ R n × R n .
The form M is usually written as a n × n matrix that is:

(i) (symmetric) ∀(u, v) ∈ R n × R n , u, Mv = v, Mu (ii) (positive) ∀u ∈ R n , u, Mu ≥ 0 (iii) (definite) u, Mu = 0 =⇒ u = 0.
These properties ensure that M defines a dot-product. In the following, the matrix M is simply called a metric tensor or a metric. The simplest example of an Euclidean metric space is given by the identity matrix I n which spans the canonical Euclidean space R n .

The dot product defined by M spans a normed vector space (R n , . M ) and a metric vector space (R n , d M (., .)) supplied by the following norm and distance definition:

• ∀u ∈ R n , u M = u, Mu • ∀(u, v) ∈ R n × R n , d M (u, v) = u -v M .
In these spaces, the length ℓ M of a segment ab = [a, b] is given by the distance between its extremities:

ℓ M (ab) = d M (a, b).
Note that this property is generally wrong for a general Riemannian metric space defined hereafter.

In an Euclidean metric space, angles and volumes are still well posed. These features are of main interest when dealing with mesh generation. These quantities are generally introduced to define quality functions. Given a bounded subset K of R n , the volume of K computed with respect to metric tensor M is:

|K| M = K √ det M dK = √ det M |K| In . (1)
The angle between two non-zero vectors u and v is defined by the unique real-value θ ∈ [0, π] verifying:

cos(θ) = u, v M u M v M .
Geometric interpretation. We will often refer to the geometric interpretation of a metric tensor. This geometric view plays an important role in the continuous mesh model. In the vicinity V(a) of a point a, the set of points, that are at a distance ε of a, is given by:

Φ M (ε) = x ∈ V(a) | t (x -a) M (x -a) ≤ ε 2 .
We note that it is sufficient to describe Φ M (1) as Φ M (ε) can be deduced from Φ M (1) for all ε by homogeneity:

Φ M (ε) = ε -1 x ∈ V(a) | x ∈ Φ M (1) .
To describe Φ M (1), the spectral decomposition

M = R Λ t R is used. R is an orthonormal matrix verifying t RR = R t R = I n .
It is composed of the eigenvectors of M. Λ is a diagonal matrix composed of the eigenvalues of M. Eigenvalues (λ i ) i=1,n are strictly positive. In the eigenvectors frame, the initial quadratic form

t (x -a) M (x -a) becomes t (x -ã) Λ (x -ã).
Consequently, Φ M (1) can be rewritten in this basis:

Φ M (1) = x ∈ V(ã) | n i=1 λ i (x i -ãi ) 2 ≤ 1 = x ∈ V(ã) | n i=1 xi -ãi h i 2 ≤ 1 .

INRIA

The last relation defines an ellipsoid centered at a with its axes aligned with the principal directions of M. Sizes along these directions are given by

h i = λ -1 2
i . We denote by E M this ellipsoid. Figure 1 depicts E M . In the sequel, the set Φ M (1) is called the unit ball of metric M and it is denoted by B M . Natural metric mapping. The last information handled by a metric tensor M is the definition of an application that maps the unit ball B In of identity metric I n onto the unit ball B M of metric M. This mapping is given by the application M -1 2 : R n → R n . M -1 2 is defined by the spectral decomposition M -1 2 = R Λ -1 2 t R, where Λ -1 2 is the diagonal matrix composed of the inverse of the square root of the eigenvalues of M. This mapping provides another description of the ellipsoid E M :

v 1 v 2 v 3
E M = M -1 2 x | x 2 2 = 1 .

Riemannian metric space

When a metric tensor field is varying smoothly in the whole domain Ω, a Riemannian metric space is defined. We denote this space by M = (M(x)) x∈Ω . The continuous mesh model is based on such a space. To give a practical visualization of a Riemannian metric space, the unit ball of the metric at some points of the domain has been drawn, see Figure 1 (right).

The main operation performed in this space is the computation of edges lengths. It is important to note that, in a Riemannian metric space, computing the length of a segment (i.e., an edge) differs from evaluating the distance between the extremities of this segment. Indeed, the straight line is no more the shortest path between two points which is given by a geodesic. To take into account the variation of the metric along the edge, the edge length is evaluated with an integral formula: Definition 1 (Edge length computation). In a Riemannian metric space M = (M(x)) x∈Ω , the length of edge ab is computed using the straight line parameterization γ(t) = a + t ab, where t ∈ [0, 1]:

ℓ M (ab) = 1 0 γ ′ (t) M dt = 1 0 t ab M(a + t ab) ab dt. ( 2
)
Figure 2 depicts iso-values of segment length from the origin for different Riemannian metric spaces. The plotted function is f (x) = ℓ M (ox) where o is the origin of the plane. The isovalues are isotropic for the Euclidean space. They are anisotropic in the case of an Euclidean metric space defined by M. The two principal directions of M clearly appear. In the case of a Riemannian metric space M, all previous symmetries are lost.

Continuous mesh model

In this section, the definition of a continuous element is first introduced. It is defined locally by considering a constant metric tensor. Then, a continuous mesh model is defined based on the notion of Riemannian metric space. In both cases, we give the set of all discrete elements or all discrete meshes that are well represented by these continuous models. These classes of equivalence are based on the notion of unit element and of unit mesh.

INRIA

Continuous element and class of unit elements

In the continuous framework, a metric tensor M is a continuous element.

Remark 1. A continuous element M is a metric tensor and it also defines an Euclidean metric space (R n , M). In the sequel, we may use one these equivalent definitions of a continuous element.

The class of all discrete elements, that are represented by this continuous model, is given by the following definition: Definition 2 (Unit element). An element K is unit with respect to a continuous element M if the length of all its edges is unit in the metric M. If K is given by its list of edges (e i ) i=1..n(n+1)/2 , then :

∀i = 1, ..., n(n + 1) 2 , ℓ M (e i ) = 1 .
The volume of K is given by:

|K| M = √ n + 1 2 n/2 n! and |K| In = √ n + 1 2 n/2 n! det(M).
Volume and length values are deduced from the Euclidean space example where the metric M is the identity matrix I n . In this case, the element is the regular n-simplex. It is composed of n + 1 vertices connected by n(n+1) 2 edges of unit length. We can prove by induction that the volume V n of the regular n-simplex is given by:

V n = √ n + 1 2 n/2 n! .
To this end, we use the following induction formulas:

t 2 n + R 2 n-1 = 1 and V n = 1 n V n-1 t n
, where t n is the length from one vertex to the center of the opposite face. The formulas are initialized by the 2D area or the 3D volume of the regular triangle or tetrahedron:

V 2 = √ 3 4 , t 2 = √ 3 2 , V 3 = √ 2 12
and t 3 = 2 3 .

Figure 3 gives two examples of unit elements for two different metric tensors. More generally, the relationships between unit discrete elements and a continuous element are stated in the following proposition:

Proposition 1 (Equivalence classes). Let M be a continuous element, there exists a nonempty infinite set of unit elements with respect to M. Conversely, given an element K = (e i ) i=1..n(n+1)/2 such that |K| In = 0, then there is a unique continuous element M for which element K is unit with respect to M.

The relation unit with respect to M defines a class of equivalence among the set of all discrete elements. Proof. We first examine the uniform case where M = I 3 . The general case is deduced from it by using the mapping M -1 2 . Let K 0 be a regular tetrahedron, thereby K 0 is unit with respect to I 3 . Whatever the rotation matrix R verifying t RR = R t R = I 3 , the tetrahedron R K 0 is still unit for I 3 . Consequently, the class of all unit elements for the continuous element I 3 is:

K = {K | ∀R ∈ O 3 : K = RK 0 } with O 3 = R | t RR = R t R = I 3 .
The equivalence class of the unit elements with respect to M is then given by the set:

M -1 2 K ∀K ∈ K .
Conversely, given a non-degenerated discrete element K = (e i ) i=1..n(n+1)/2 , such that |K| I3 = 0, let us demonstrate that there exists a unique metric for which K is unit. It is sufficient to solve the following linear system:

(S)      ℓ 2 M (e 1 ) = 1 . . . ℓ 2 M (e n(n+1)/2 ) = 1 .
The determinant of (S) is equal to |K| I3 = 0. Consequently, (S) admits a unique solution.

Figure 4 depicts some unit elements with respect to a continuous element.

INRIA

Geometric invariants

So far, only Definition 2 has established relationships between unit elements and continuous elements. Other properties exist. They connect the geometric properties of unit elements to the linear algebra properties of metric tensors. The following proposition gives geometric invariants that hold for all unit elements with respect to a continuous element.

Proposition 2 (Geometric invariants). Let M be a continuous element and K be a unit element with respect to M. We denote by (e i ) i its edges list, see conventions in Figure 5, and by |K| its Euclidean volume. Then, the following invariants hold:

• standard invariants:

∀ (e i , e j ), t e i M e i = 1,

2 t e i M e j + 1 = 0 if i = j. (3)
Figure 4: Several unit elements with respect to a continuous element in 2D and 3D. • invariant related to the Euclidean volume |K|:

|K| = √ 3 4 det(M -1 2 ) in 2D and |K| = √ 2 12 det(M -1 2 ) in 3D. (4)
• invariant related to the square length of the edges for all symmetric matrix H:

3 i=1 t e i He i = 3 2 trace(M -1 2 HM -1 2 ) in 2D, 6 i=1 t e i He i = 2 trace(M -1 2 HM -1 2 ) in 3D.
(

) 5
Proof. The first invariant of Relation (3) comes from the definition of a unit element.

The second invariant of Relation [START_REF] Berger | A panoramic view of Riemannian geometry[END_REF] states that the angle between two edges of a unit element face is constant in the metric. Let (e i , e j ) be a couple of edges of element K, this couple defines a face. We denote by e k the third edge of this face. According to the conventions depicted in Figure 5, these edges verify: e i + e je k = 0. Expanding the following relation

t (e i + e j + e k ) M (e i + e je k ) = 0, leads to the second invariant of Relation (3). Invariant ( 4) is proved by a direct integration. Given a unit element K for M, there exists a unique regular tetrahedron K 0 , which is unit with respect to the identity matrix I 3 , such that K = M -1 2 K 0 . The volume of K is then given by:

|K| = K 1 dx = K0 det(M -1 2 ) dx = det(M -1 2 ) |K 0 |, where |K 0 | = √ 2
12 . The same proof applies in 2D. Invariant [START_REF] Bottasso | Anisotropic mesh adaption by metric-driven optimization[END_REF] is first proved in the simpler case where H = I 3 and M = I 3 . The general case will be deduced from this proof. Let us consider the regular tetrahedron K 0 = (v 1 , v 2 , v 3 , v 4 ) unit for I 3 defined by the list of vertices:

v 1 = (0, 0, 0) , v 2 = (1, 0, 0) , v 3 = 1 2 , √ 3 2 , 0 and v 4 = 1 2 , √ 3 6 , 2 3 .
The proof does not depend on this specific choice of coordinates. We first demonstrate the following preliminary result: For all lines (D) passing through one of the vertices of K 0 , the sum of the square lengths of the edges projected on (D) is invariant. Without loss of generality, we assume that line (D) passes through the vertex v 1 of K 0 . If line (D) is defined by the vector n = (cos(u) cos(v), cos(u) sin(v), sin(u)) ,

INRIA

with (u, v) ∈ R 2 , then the length of the first three edges of K 0 projected on (D) are given by: a

= e 1 . n = cos(u) cos(v), b = e 2 . n = 1 2 cos(u) cos(v) + √ 3 2 cos(u) sin(v), c = e 3 . n = 1 2 cos(u) cos(v) + √ 3 6 cos(u) sin(v) + 2 3 sin(u),
with conventions of Figure 5. A direct trigonometric calculus shows that the sum of the square length of all the edges projected on (D) is equal to 2. Indeed, it comes:

Σ = a 2 + b 2 + c 2 + (b -c) 2 + (c -a) 2 + (a -b) 2 = 3 a 2 + 3 b 2 + 3 c 2 -2ab -2ac -2bc.
After expanding and factorizing, Σ reads:

Σ = 2 cos(u) 2 cos(v) 2 + 2 cos(u) 2 sin(v) 2 + 2 sin(u) 2 = 2.
When M is different from I 3 , we use the mapping M -1 2 that maps the unit ball of I 3 onto the unit ball of M. As regards line (D), we select the specific line which has for direction vector one the main direction of M, e.g. u j , and which is passing through v 1 . The lengths a, b and c are thus multiplied by h j which is the size prescribed by M in the direction u j . Consequently, the square length of the edges projected on (D) are multiplied by h 2 j . It comes:

Σ j = 6 i=1 |e i . u j | 2 = h 2 j Σ = 2 h 2 j .
Considering the previous relation for all the principal directions of M and summing the Σ j complete the proof:

6 i=1 e i 2 2 = 3 j=1 6 i=1 |e i . u j | 2 = 2 h 2 1 + h 2 2 + h 2 3 = 2 trace(M -1 ). ( 6
)
We now consider the more general case where H is symmetric definite positive. The matrices H 1 2 and H -1 2 are well defined and are symmetric. We first prove that if K is a unit element for M then H

1 2 K is a unit element with respect to M 1 2 H -1 M 1 2 . Indeed, if we consider the edge ẽj = H 1 2 e j of H 1 2
K where e j is a unit length edge of K with respect to M, it comes:

t ẽj M 1 2 H -1 M 1 2 ẽj = t (H 1 2 e j ) M 1 2 H -1 M 1 2 (H 1 2 e j ) = 1.
Then, from Relation (6) we get the last invariant:

6 i=1 t e i He i = 6 i=1 ẽi 2 2 = 2 trace(M -1 2 HM -1 2 ) .

RR n°6846

Other geometric invariants can be found in [START_REF] Loseille | Adaptation de maillage 3D anisotrope multi-échelles et ciblé à une fonctionnelle. Application à la prédiction haute-fidélité du bang sonique[END_REF].

Continuous mesh and class of unit meshes

A continuous mesh of a domain Ω ⊂ R n is a Riemannian metric space M = (M(x)) x∈Ω .

We recall the spectral decomposition of M(x):

M : x ∈ Ω → M(x) = R(x)    λ 1 (x) . . . λ n (x)    t R(x) = R(x)    h -2 1 (x) . . . h -2 n (x)    t R(x) ,
where R(x) is an orthonormal matrix providing the local orientation, (λ i (x)) i=1,n are the local eigenvalues and (h i (x)) i=1,n are the local sizes along the principal directions of the continuous mesh. Practically, another decomposition is used that points out the local characteristics of the continuous mesh. This decomposition is given by the following proposition.

Proposition 3 (Continuous mesh).

A continuous mesh M = (M(x)) x∈Ω locally writes:

M(x) = d 2 n (x) R(x)     r -2 n 1 (x) . . . r -2 n n (x)     t R(x),
where

• the density d is equal to: d = n k=1 h k -1 = n k=1 λ k 1 2
,

• the n anisotropic quotients r i are equal to:

r i = h n i n k=1 h k -1
.

Proof. The proof consists in computing d

2 n r -2 n i : d 2 n r -2 n i = n k=1 h k -2 n h -2 i n k=1 h k 2 n = h -2 i = λ i .

INRIA

The density d controls only the local level of accuracy of the continuous mesh. Increasing or decreasing d does not change the anisotropic properties or the orientation, see Figure 6 (left). In 3D, anisotropic quotients arises from the quotient of different parallelepipeds, see Figure 6 (right).

We also define the complexity C of a continuous mesh:

C(M) = Ω d(x) d x = Ω det(M(x)) d x.
This real-value parameter is useful to quantify the global level of accuracy of the continuous mesh (M(x)) x∈Ω . It can also be interpreted as the continuous counterpart of the number of vertices of a discrete mesh. This quantity also leads to the definition of sequence of continuous embedded meshes. Two continuous embedded meshes have the same anisotropic ratios and orientations. They only differ from their complexity: Definition 3 (Embedded continuous meshes). Two continuous meshes (M(x)) x∈Ω and (N (x)) x∈Ω are embedded if a constant c exists such that:

∀x ∈ Ω, N (x) = c M(x).
Conversely, from a continuous mesh M = (M(x)) x∈Ω , we can deduce a continuous mesh N = (N (x)) x∈Ω of complexity N having the same anisotropic properties (anisotropic orientations and ratios) by considering:

N (x) = N C(M) 2 n M(x).
In the context of error estimation, this notion enables the study of the order of convergence with respect to an increasing complexity N . Consequently, the complexity C(M) is also the continuous counterpart of the classical parameter h used for uniform meshes while studying convergence. In the continuous mesh framework, the uniform refinement consisting in dividing by two each edge of a uniform mesh of size h writes:

M i = 4 i   1 h 2 1 h 2 1 h 2   ,
where i is the level of refinement. (M i ) i=1...k defines a sequence of continuous embedded meshes. Consequently, this simple practical adaptive strategy has a simple continuous interpretation in term of embedded continuous meshes. Such a uniform refinement is exemplified on Figure 7. However, when dealing with anisotropic meshes, a unique size h is no more sufficient to give a quantitive information on the accuracy. The size h is then replaced by the continuous mesh complexity. Figure 7: In a sequence of uniformly refined uniform meshes, a single size information h is sufficient to describe the accuracy of the current mesh in the whole domain. On the contrary when the mesh involves strong differences in sizes and orientations, this size is replaced by another measure, the complexity.

h 1 h 2 h 3

INRIA

Unit mesh. The notion of unit mesh is far more complicated than the notion of unit element as the existence of a mesh composed only of unit regular simplexes with respect to a given continuous mesh is not guaranteed. For instance, if the continuous mesh is not compatible with the domain size, then it clearly does not exist such discrete mesh. To avoid this problem, let us look at the existence of a discrete mesh composed only with unit regular simplexes with respect to a continuous mesh in R n . To simplify even more the problem, we first consider the continuous mesh (

I n (x)) x∈R n .
It is well known that R3 cannot be filled only with the regular tetrahedron while it is possible to fill R 2 with the equilateral triangle. Consequently, even for the simplest continuous mesh (I 3 (x)) x∈R 3 , there is no discrete mesh composed only of the unit regular tetrahedron. Therefore, the notion of unit mesh has to be released:

Definition 4 (Unit mesh). A discrete mesh H of a domain Ω ⊂ R n is unit for a continuous mesh (M(x)) x∈Ω if all its elements are quasi-unit.
Now, let us give a meaning to quasi-unit in three dimensions. A first way to release the definition of unity is to take into account technical constraints imposed by mesh generators. To converge (and to avoid cycling) while analyzing edges length, the meshing algorithm considers an admissible edge length interval of the form [ 1 α , α] with α > 0 [START_REF] Frey | Yams, a fully automatic adaptive isotropic surface remeshing procedure[END_REF]. If the symmetry property is required, i.e., α 2 = 1 α , then we obtain α = √ 2. Therefore, as regards the meshing requirement, a tetrahedron K defined by its list of edges

(e i ) i=1...6 is said quasi-unit if ∀i ∈ [1, 6], ℓ M (e i ) ∈ [ 1 √ 2 , √ 2
]. Nevertheless, we do not know if this definition provide the existence of a unit mesh for the continuous mesh (I 3 (x)) x∈R 3 . In the following, this question of existence is studied by means of the space filling tetrahedra.

Non-regular space filling tetrahedra. The study of space filling tetrahedra is an old geometrical question [START_REF] Naylor | Filling space with tetrahedra[END_REF][START_REF]Which tetrahedra fill space?[END_REF]. In the past, it has been demonstrated that there exist sets of non-regular space filling tetrahedra: the Sommerville tetrahedra [START_REF] Sommerville | Space-filling tetrahedra in euclidean space[END_REF] and the Goldberg tetrahedra family [START_REF] Goldberg | Three infinite families of tetrahedral space-fillers[END_REF].

The Sommerville tetrahedra are based on particular splittings of the unit cube, see Figure 8. We recall these tetrahedra thanks to their vertices coordinates, only the last vertex distinguishes them.

K is denoted (v 1 , v 2 , v 3 , v 4 ) with v 1 = (0, 0, 0), v 2 = ( 1 2 , -1 2 , 1 2 ), v 3 = ( 1 2 , 1 2 , 1 2 ) and • v 4 = ( 1 2 , 0, 0) for the Sommeville tetrahedron 1 • v 4 = (1, 0, 0) for the Sommeville tetrahedron 2 • v 4 = ( 1 2 , -1 2 , - 1
2 ) for the Sommeville tetrahedron 3

• v 4 = ( 1 2 , 0, - 1 4
) for the Sommeville tetrahedron 4. The Goldberg tetrahedra are based on the splitting of a prism the basis of which is the equilateral triangle, see Figure 9. Their coordinates depend on an initial choice of two lengths a and e. We specify one of the Goldberg tetrahedra for the specific choice a = 1 • v 1 = (0, 0, 0), v 2 = (0, 0, 1), v 3 = (0, 1, 1

3 ) and v 4 = (

√ 3 2 , 1 2 , 2
3 ) for the Goldberg tetrahedron. We propose now to compare these space filling tetrahedra to the unit regular tetrahedron. To this end, all these tetrahedra are scaled such that their volumes are equal to

√ 2
12 . The resulting edges lengths for each tetrahedron are specified in Table 1. We notice that the proposed notion of quasi-unit element is only verified for the Sommerville tetrahedra 1 and 2, and the Goldberg tetrahedron. Therefore, there exists space filling tetrahedra that are quasi-units for the metric I 3 in the sense proposed above.

INRIA

Now, the case of a constant anisotropic metric M is studied. We consider the pattern around a vertex, i.e., the vertex ball, composed only with the second Sommerville tetrahedron. This pattern exists as it fills (R 3 , M). The vertex ball is mapped back in the natural Euclidean space thanks to the application M 1 2 . Then, we notice that the non-regularity of the second Sommerville tetrahedron leads necessarily to the creation (in the Euclidean space) of several different anisotropic tetrahedra. However, all these different tetrahedra have the same edges lengths and the same volume in the metric M. Consequently, filling space with only one tetrahedra is possible for all isotropic metrics of the form αI 3 , but a set of tetrahedra is required to fill the Euclidean space anisotropically.

Controlling the volume. Unfortunately, the weaker constraint on the edges length can lead to the generation of quasi-unit elements with a null volume. For instance in (R 3 , I 3 ), the regular tetrahedron with edges length equal to √ 2 is quasi-unit for I 3 . However, if one of its vertex is projected orthogonally on the opposite face, then a quasi-unit element of null volume is obtained. Indeed, three edges are of length √ 2 and the three other are of length

√ 3 6 ≈ 0.816 ∈ 1 √ 2 , √ 2 .
In consequence, controlling only the edges length is not sufficient, the volume must also be controlled to release the notion of unit element. Practically, this is achieved by prescribing a quality function:

Q M (K) = 36 3 1 3 |K| 2 3 M 6 i=1 ℓ 2 M (e i ) ∈ [0, 1] . (7)
For the perfect regular tetrahedron, whatever its edges length, the quality function is equal to 1. For a null volume tetrahedron, Q M is 0. The qualities of the space filling tetrahedra are given in Table 1. Notice that Q M only quantifies the gap to the regular tetrahedron shape.

We deduce the following definition of quasi-unit element, which is also practically used by mesh generators,

Tetrahedron

Coeff. The edges length and the tetrahedron quality Q I3 given by Formula (7) are provided.

Edges length

Definition 5 (Quasi-unit element). A tetrahedron K defined by its list of edges

(e i ) i=1...6 is said quasi-unit for M if ∀i ∈ [1, 6], ℓ M (e i ) ∈ 1 √ 2 , √ 2 and Q M (K) ∈ [α, 1] with α > 0 .
In our case, α = 0.8 is an acceptable value as it enables the Sommerville tetrahedra 1 and 2, and the Goldberg tetrahedron to be generated.

Remark 2. Instead of considering Q M , the quality function 1 Q M can be considered. As the variation range becomes [1, ∞[, the discrimination of bad elements is made easier.

Continuous linear interpolation error

In the previous section, a continuous framework has been introduced to model elements and meshes. Now, we aim at applying this framework in the context of error estimation. Let (M(x)) x∈Ω be a continuous mesh of a domain Ω and let u be a non linear function which is assumed to be only twice continuously differentiable. We seek a well-posed definition of the continuous linear interpolation error uπ M u L 1 (Ω) related to a continuous mesh (M(x)) x∈Ω which implies a well-posed definition of a linear continuous interpolate π M u. More precisely, we would like the continuous linear interpolation error to be a reliable mathematical model of u -Π h u L 1 (Ω h ) where Π h is defined by a mesh H of a discretized domain Ω h which is a unit mesh with respect to (M(x)) x∈Ω . This means that considering

u -π M u L 1 (Ω) is equivalent to consider u -Π h u L 1 (Ω h ) .
The error analysis is first done locally, i.e, in a tangent space of (M(x)) x∈Ω at a given point a. In the tangent space, the continuous mesh (M(x)) x∈Ω reduces to the continuous element M(a), i.e., the analysis is performed locally at the element level. The function u is approximated by its local quadratic taylor expansion. Indeed, terms of order greater than two can be neglected while studying the linear interpolation error. Then, a global error estimate is derived by taking into account the variation of the continuous mesh and of the function.

Interpolation error in L 1 norm for quadratic functions

In this section, we consider a quadratic function u defined on a domain Ω ⊂ R 3 . The function is given by its matrix representation:

∀x ∈ Ω, u(x) = 1 2 t x H x,
where H is a symmetric matrix representing the Hessian of u. For every symmetric matrix H, |H| denotes the positive symmetric matrix deduced from H by taking the absolute values of its eigenvalues. The function u is linearly interpolated on an tetrahedron K defined by its vertices list:

K = (v 1 , v 2 , v 3 , v 4 ).
|K| denotes the Euclidean volume of K. We denote by Π h u the linear interpolate of u on K. We can now state the following result:

INRIA Proposition 4.
For every quadratic function u, its linear interpolation error in L 1 norm on an element K verifies:

u -Π h u L 1 (K) ≤ |K| 40 6 i=1 t e i |H|e i ,
where {e i } i=1,6 is the set of edges of K.

The previous inequality becomes an equality when u is elliptic or parabolic.

Proof. The proof consists in deriving an exact error estimate of the point-wise interpolation error within the element K: e(x) = (u -Π h u)(x) for x ∈ K. This error is then integrated over K.

To derive e, we use the standard reference element technique. Reference element K ref is defined by its four vertices coordinates:

v 1 = t (0, 0, 0), v 2 = t (1, 0, 0), v 3 = t (0, 1, 0) and v 4 = t (0, 0, 1).
All the computations are done on K ref and the result is then mapped onto the current element K by using the following affine mapping:

x = v 1 + B K x with B K = (v 2 -v 1 , v 3 -v 1 , v 4 -v 1 ), x ∈ K, x ∈ K ref .
The matrix B k is given as a function of the following edges:

e 1 = v 2 -v 1 , e 2 = v 3 -v 1 and e 3 = v 4 -v 1 ,
so that B K = (e 1 , e 2 , e 3 ). The quadratic function u reads in the frame of K ref :

u(x( x)) = 1 2 t v 1 H v 1 + 1 2 t v 1 H B K x + 1 2 t x t B K H v 1 + 1 2 t x t B K H B K x.
As we consider the linear interpolation, linear and constant terms of u(x( x)) are exactly interpolated. Without loss of generality, these terms are neglected and only quadratic terms are kept. Indeed, if we consider ũ(x) = 1 2 t x t B K H B K x, then it comes:

e(x) = (u -Π h u)(x) = (ũ -Π h ũ)(x).
We can now consider ũ instead of u. However for the sake of clarity, we keep on writing u and not ũ. We rewrite u in a matrix form:

u(x( x)) = 1 2 t   x ŷ ẑ     t e
v i ) = u(v i ) for all i ∈ [1, 4]:              Π h u(v 1 ) = d = u(x((0, 0, 0)) = 0, Π h u(v 2 ) = a = u(x((1, 0, 0)) = 1 2 ( t e 1 He 1 ), Π h u(v 3 ) = b = u(x((0, 1, 0)) = 1 2 ( t e 2 He 2 ), Π h u(v 4 ) = c = u(x((0, 0, 1)) = 1 2 ( t e 3 He 3
). The solution of the previous linear system gives the final expression of Π h u:

Π h u(x( x)) = 1 2
( t e 1 He 1 ) x + ( t e 2 He 2 ) ŷ + ( t e 3 He 3 ) ẑ .

The exact point-wise interpolation error e(x) is then given by:

e(x( x)) = 1 2 [ ( t e 1 He 1 ) (x 2 -x) + ( t e 2 He 2 ) (ŷ 2 -ŷ) + ( t e 3 He 3 ) (ẑ 2 -ẑ) + 2 ( t e 1 He 2 ) xŷ + 2 ( t e 1 He 3 ) xẑ + 2 ( t e 2 He 3 ) ŷẑ ].
From this equality, estimate in L 1 , L 2 or H 1 can be deduced by considering the change of variables given by the mapping B K . Indeed, for every function F , its integration over K can be computed through its expression in K ref :

K F (x) dxdydz = K ref F (x( x)) | det(B K )| dxdŷdẑ,
As 6 |K| = det(B k ), previous equality becomes:

K F (x) dxdydz = 6|K| K ref F (x( x)) dxdŷdẑ.
Consequently, the interpolation error in L 1 norm is evaluated by a direct integration of |e(x)|. When u is concave or convex, we have:

|(u -Π h u)(x)| = (u -Π h u)(x) in the convex case and |(u -Π h u)(x)| = -(u -Π h u)(x)
in the concave case. The error reads: Using the conventions of Figure 5, the crossed terms can be expressed only in terms of e i H e i for i = 1, .., 6: If u is hyperbolic, the following inequality is used:

u -Π h u L 1 (K) = |K| 40
       2 t e
1 2 |x H x| ≤ 1 2 x |H| x,
to conclude the proof in the general case.

The same proof applies in 2D. For a quadratic function u, the linear interpolation error on a triangle K is given by:

u -Π h u L 1 (K) ≤ |K| 24
3 i=1 t e i |H|e i .
Error estimates in L 2 norm and in H 1 norm can also be derived. We refer to [START_REF] Berzins | Solution-based mesh quality for triangular and tetrahedral meshes[END_REF][START_REF] Nadler | Piecewise linear approximation on triangulations of a planar regions[END_REF] and references therein for their evaluations. These error estimates are classically used to exhibit mesh quality functions and to obtain the best element shape minimizing the interpolation error. In the continuous mesh framework, the interpolation error estimate in L 1 norm is used to prove some exactness properties of the continuous linear interpolate.

Remark 3 (Safety principle). Even if it is possible to define exactly the linear interpolation error in L 1 norm for hyperbolic functions, we do not consider these expressions from a practical point of view. We prefer to consider |H| instead of H transforming the function into an elliptic or a parabolic one. It comes to in over-estimating the error for hyperbolic functions. Indeed, it seems that we do not take any advantages of considering the null error directions. However, the maximal error directions are those of the gradient of u. These directions correspond to the eigenvectors direction of H. All these choices are illustrated for the 2D example where the function x 2y 2 is considered. In that case, the following elliptic approximates can be used:

|x 2 -y 2 | ≤ x 2 + y 2 |x 2 -y 2 | ≤ r 2 (x -y) 2 + 2 r (x + y) 2 for all r > 0 |x 2 -y 2 | ≤ r 2 (x + y) 2 + 2 r (x -y) 2 for all r > 0 .
The first one is the approximation retained in this paper. The others are elliptic form aligned with null directions of the hyperbolic function. The interesting point to note is that, whatever the considered bound, the unit balls of the three elliptic bounds have the same area. Consequently, even if the safety principle consists in over-estimating the interpolation error in the hyperbolic case, it does not result in over-meshing.

Linear interpolation on a continuous element

Let M be a continuous element and u be a quadratic positive function (see Remark 3).

We study the interpolation error for the class of all unit discrete elements in M, given by Definition 2. Figure 4 depicts for a given metric tensor M some unit elements. We can now state the main result:

Theorem 1. For all unit element K with respect to M, the interpolation error of u in L 1 norm does not depend on the element shape and is only a function of the Hessian H of u and of the metric M.

• In 3D, for all unit elements K in M, the following equality holds:

u -Π h u L 1 (K) = √ 2 240 det(M -1 2 ) trace(M -1 2 H M -1 2 ). ( 8
)
• In 2D, for all unit elements K in M, the following equality holds:

u -Π h u L 1 (K) = √ 3 64 det(M -1 2 ) trace(M -1 2 H M -1 2 ).
Proof. According to Proposition 4, the interpolation error in L 1 norm of a quadratic positive function u on an element K is:

u -Π h u L 1 (K) = |K| 40 6 i=1 t e i He i .
Then, if K is unit with respect to M, the previous interpolation error is expressed by:

u -Π h u L 1 (K) = √ 2 240 det(M -1 2 ) trace(M -1 2 H M -1 2 ).
thanks to the geometric invariants related to the volume, Relation (4), and to the square lengths of the edges, Relation [START_REF] Bottasso | Anisotropic mesh adaption by metric-driven optimization[END_REF].

We note the strong analogy with classical interpolation error estimate for Lagrange interpolation [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]:

• The term det(M - 1
2 ) stands for the Jacobian of the affine transformation from the reference element K onto the current element K. In our continuous framework, it is the Jacobian of the affine mapping between the reference continuous element unit ball B I3 onto the current continuous element unit ball B M .

• The term trace(M -1 2 H M -1 2 ) stands for the semi-norm involved in classical error estimates. Generally, this semi-norm contains the anisotropic behavior of the estimate. In the continuous framework, the trace-term gives the alignment correlation between the principal directions of the Hessian H and the principal directions of the metric M.

INRIA

Relation [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] shows that the infinite set of discrete elements that are unit for a given metric M achieves the same interpolation error, and moreover, shows that this interpolation error is only expressed with continuous quantities: the continuous element M and the Hessian of the function u. Consequently, Theorem 1 points out that the metric alone contains enough information to describe completely the linear interpolation error in L 1 norm. In other words, this theorem confirms that the use of metric-based mesh adaptation is particularly well suited to control anisotropically the interpolation error. In the past, this efficiency has been observed practically on real life problems, see for instance [START_REF] Bottasso | Anisotropic mesh adaption by metric-driven optimization[END_REF][START_REF] Dompierre | Anisotropic mesh adaptation: towards a solver and user independent cfd[END_REF][START_REF] Formaggia | Anisotropic mesh adaptation in computational fluid dynamics: Application to the advection-diffusion-reaction and the Stokes problems[END_REF][START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF][START_REF] Loseille | Achievement of second order mesh convergence for discontinuous flows with adapted unstructured mesh adaptation[END_REF][START_REF] Pain | Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations[END_REF][START_REF] Schall | Mesh adaptation as a tool for certified computational aerodynamics[END_REF][START_REF] Tam | Anisotropic mesh adaptation for 3D flows on structured and unstructured grids[END_REF].

Continuous linear interpolate

The main difficulty in defining the continuous linear interpolate is to connect a discrete error computed on an element to a local continuous error that is defined point-wise. Indeed, the discrete interpolation error in norm L 1 is integrated on the element K. On the contrary, a continuous mesh is a function x → M(x) defined at each point x of Ω.

Suppose now that the continuous mesh (M(x)) x∈Ω is varying and that the function u is no more quadratic but only twice continuously differentiable. If Equality (8) of Theorem 1 does not hold anymore, all the terms of the right-hand-side M and H are well defined continuously. The definition of a continuous interpolate follows up from this consideration.

We denote by u Q the quadratic approximation of a smooth function u. At point a, u Q is defined in the vicinity of a as the truncated second order Taylor expansion of u:

∀x ∈ V(a) , u Q (a; x) = u(a) + ∇u(a)(x -a) + 1 2 (x -a), H(a)(x -a) .
When no confusion is possible, the notation u Q (a; x) is replaced by u Q (x). u Q is a complete quadratic form composed of a constant term, a linear term and finally a quadratic term.

The first result of this section provides an equivalence formula between discrete and continuous views locally around a point a of the domain. In the vicinity of a, u Q approximates u and (M(x)) x∈Ω reduces to M(a) in the tangent space. We can now state the main result: Theorem 2 (Discrete-continuous equivalence). Let u be a twice continuously differentiable fonction of a domain Ω and (M(x)) x∈Ω be a continuous mesh of Ω. Then, there exists a unique continuous linear interpolate function π M such that:

∀a ∈ Ω , |u -π M u|(a) = 2 u Q -Π h u Q L 1 (K) |K| ,
for every K unit element with respect to M(a).

The proof is given hereafter, we first look at the consequences of this theorem.

RR n°6846

Corollary 1. Let u be a twice continuously differentiable fonction of a domain Ω and (M(x)) x∈Ω be a continuous mesh of Ω. Then, the following continuous linear interpolation estimate holds in 3D:

∀a ∈ Ω , |u -π M u|(a) = 1 10 trace M(a) -1 2 |H(a)| M(a) -1 2 = 1 10 d(a) -2 3 3 i=1 r i (a) 2 3 t u i (a) |H(a)| u i (a) .
In 2D, the estimate is:

∀a ∈ Ω , |u -π M u|(a) = 1 8 trace M(a) -1 2 |H(a)| M(a) -1 2 = 1 8 d(a) -1 2 i=1 r i (a) t u i (a) |H(a)| u i (a) .
Proof. In 3D, for all unit elements K with respect to M(a), the error estimation (8) can be rewritten as follow for the quadratic function u Q approximating u in the vicinity of a:

u Q -Π h u Q L 1 (K) |K| = 1 20 trace(M(a) -1 2 |H(a)| M(a) -1 2 ).
Then, expressing M(a) as a function of the continuous mesh parameters given by the decomposition of Proposition 3 leads to:

u Q -Π h u Q L 1 (K) |K| = 1 20 d(a) -2 3 3 i=1 r i (a) 2 3 t u i (a) |H(a)| u i (a)
where the (u i (a)) i=1,3 stand for the eigenvectors of M(a).

This result shows that the continuous point-wise linear interpolation can be decomposed into the product of two terms:

• a first term that control the accuracy, this density term is directly connected to the size of the continuous element,

• a second term that measures alignment deviation between the continuous element orientation and the anisotropy features of the function u.

It is possible to give a geometric interpretation of this estimate. This interpretation illustrates the impact of a continuous mesh on the error iso-values and, consequently, gives some clue toward the control of the error by means of a continuous mesh. The term M -1 2 (a) H(a) M -1 2 (a) corresponds to the frame change related to the continuous mesh local orientation. Given a symmetric matrix |H(a)|, the corresponding quadratic form is:

f = 1 2 t x |H(a)| x.

INRIA

The matrix M -

f (x) = t x |H(a)| x = t (M -1 2 (a) x) |H(a)| M -1 2 (a) x = t x M -1 2 (a) |H(a)| M -1 2 (a) x.
Iso-values of f and f are different when they are seen in the canonical Euclidean space. In fact, viewing a quadratic form for different continuous meshes changes its iso-values in the Euclidean space, i.e., the real physical space. It is then possible to control the error isovalues by modifying M. This is the main principle of mesh adaptation, but here formulated in a continuous framework. Classical metric-based mesh adaptation consists in finding a metric field that provides isotropic iso-values of the error function, see pionneer work [START_REF] Castro-Díaz | Anisotropic unstructured mesh adaptation for flow simulations[END_REF].

These results demonstrate that both the interpolation error and the linear interpolate Π h have continuous counterparts. It is then a step forward in finding a complete analogy between the discrete and the continuous views. From a practical point of view, we deduce the following analogy. Given a unit mesh H of a domain Ω h with respect to a continuous mesh (M(x)) x∈Ω , the global interpolation error is:

u -Π h u L 1 (Ω h ) = K∈H u -Π h u L 1 (K) . (9)
In the continuous case, the discrete summation becomes an integral:

u -π M u L 1 (Ω) = Ω |u -π M u|(x) dx. ( 10
)
Note that there is no global guarantee on the continuous interpolation error reliability given by Relation [START_REF] Coupez | Génération de maillages et adaptation de maillage par optimisation locale[END_REF]. For instance, there is no a priori relationship between ( 9) and [START_REF] Coupez | Génération de maillages et adaptation de maillage par optimisation locale[END_REF]. The only guarantee is the local equivalence given by Theorem 2. However, the local guarantee becomes global when the mesh is unit with respect to a constant metric tensor (this does not necessary implied that the mesh is uniform) and when the function is quadratic. In this specific case, by neglecting error due to the boundary discretization, we have the equality:

2 u -Π h u L 1 (Ω h ) = u -π M u L 1 (Ω) ,
for all unit meshes H with respect to (M(x)) x∈Ω . The numerical examples of Section 4 will numerically demonstrate the efficiency of the continuous model. In particular, we will observe that:

• the model is accurate and the equivalence (9)≈( 10) is observed for non quadratic functions and non-constant continuous meshes,

• the error due to the fact that mesh generator generates edges with length not stricly equal to one is negligible. In particular, edges length range given in Definition 5 ensures reliable numerical results.

The constant 2 involved in Theorem 2 arrises in the perfect case where the mesh is only composed of perfect unit elements. Note this is only possible in 2D whereas it is always false in 3D due the impossibility to tessel the space only with the regular tetrahedron. In a more practical situation, this constant needs to be evaluated from one unit mesh in order to estimate the deviation between the continuous mesh complexity with respect to the number of vertices of the unit mesh. In other words, the number of vertices N v of a unit mesh verifies the following function:

N v = C N,
where N is the continuous mesh complexity and C a constant. The constant C depends on the domain shape, the mesh generator used and the unit mesh resulting quality. Consequently, it reflects how far the current mesh is from the perfect unity. Several examples are given in the numerical examples section and illustrates this relation.

To conclude this section, the proof of Theorem 2 is given. This proof is based on the exact expression of the continuous linear interpolate: Proposition 5. The continuous interpolate π M u evaluated at a ∈ Ω for a continuous mesh (M(x)) x∈Ω and for a smooth function u is given by:

π M u(a) = p * (0),
where p * is the unique linear polynomial solution of:

p * = min p∈P 1 u Q -p L 2 (B M ) ,
where u Q is the quadratic model of u at a and B M is the unit ball of M at a. π M is given by

π M u(a) = u(a) + ∇u(a) + 1 c n trace(M -1 2 (a) H(a) M -1 2 (a)),
where c n is a constant that depends only on the space dimension:

c 2 = 1 8 and c 3 = 1 10 .
Proof. The quadratic model u Q of u at point a defined by:

u Q (a; x) = u(a) + ∇u(a)(x -a) + 1 2 (x -a), H(a)(x -a) , becomes after the translation x → x + a u Q (x) = 1 2 t x H(a) x + ∇u(a) x + u(a).

INRIA

The linear polynomial p * is given by:

p * : x ∈ B M → t g x + c,
where c ∈ R and g ∈ R 3 . As the space of linear polynomials P 1 and the unit ball B M are convex, p * exists and is unique. We seek for p * verifying the following condition:

∀p ∈ P 1 , B M (u Q (x) -p * (x)) p(x) dx = 0.
In particular, it is true for the following basis of P 1 :

x → 1, x → x 1 , x → x 2 and x → x 3 .
The previous condition leads to:

B M (u Q (x) -p * (x)) dx = 0, B M (u Q (x) -p * (x)) x i dx = 0, (11)
for i = {1, 2, 3}. The initial integration domain B M is mapped onto the unit sphere B I3 by using the following one-to-one change of variables:

B M -→ B I3 x -→ y = M(a) 1 2 x.
Functions u Q and p * becomes:

u Q (x) = ũQ (y) = 1 2 t y M(a) -1 2 H(a) M(a) -1 2 y + t y M(a) -1 2 ∇u(a) + u(a), p * (x) = p * (y) = t y M(a) -1 2 g + c.
We now consider the following basis:

y → 1, y → y 1 , y → y 2 and y → y 3 , Equations ( 11) become:

B I 3 (ũ Q (y) -p * (y)) det(M(a) -1 2 ) dy = 0, B I 3 (ũ Q (y) -p * (y)) y i det(M(a) - 1
2 ) dy = 0, for i = {1, 2, 3}. By using integration formula of Annexe A, it comes in 3D:

2 15 trace(M(a) -1 2 H(a) M(a) -1 2 ) + 4 3 (u(a) -c) det(M(a) -1 2 ) = 0, 4 15 M(a) -1 2 (∇u(a) -g) det(M(a) -1 2 ) = 0,
from which the expression of p * is deduced:

g = ∇u(a), c = u(a) + 1 10 trace(M(a) -1 2 H(a) M(a) -1 2 ).
Finally, for the 2D case, integration formula of Annexe A provides:

1 8 trace(M(a) -1 2 H(a) M(a) -1 2 ) + (u(a) -c) det(M(a) -1 2 ) = 0, 1 4 M(a) -1 2 (∇u(a) -g) det(M(a) -1 2 ) = 0.
and the expression of p * is given by:

g = ∇u(a), c = u(a) + 1 8 trace(M(a) -1 2 H(a) M(a) -1 2 ).
The proof of Theorem 2 is deduced from the definition of π M u(a).

In the case where the continuous mesh is constant and the function u quadratic, we verify (9)= [START_REF] Coupez | Génération de maillages et adaptation de maillage par optimisation locale[END_REF].

Notice that using the L 2 projection of the quadratic model u Q of u is necessary to ensure the specific equivalence between the discrete linear interpolate Π h and the continuous linear interpolate π M of Theorem 2. However, the continuous linear interpolate is still well defined if we use the function u instead of u Q . It seems then possible to define the continuous linear interpolate for far less regular functions. For instance, one may consider only functions that are locally L 2 . An open problem is then to find a discrete linear interpolation operator which enables discrete and discontinuous approaches to be linked. Some works involving other interpolation operators have already been considered, we can cite the developments in [START_REF] Formaggia | New anisotropic a prioiri error estimate[END_REF] that derive interpolation error estimate based on the Clément's interpolate [START_REF] Clément | Approximation by finite element functions using local regularization[END_REF].

Numerical examples

Using the continuous framework introduced above, the continuous interpolation error in L 1 norm of a function over a domain Ω can be computed analytically for any function u and any continuous mesh (M(x)) x∈Ω that are defined by analytic functions. We exemplify this continuous view on examples in 2D and 3D. This calculus does not require any discrete support, e.g. any mesh. To validate the approach, each continuous error is compared to the discrete interpolation error computed on a unit mesh with respect to the continuous one.

Note. In this section, the analytical integrations are computed by using symbolic calculus, when possible, with Maple.

INRIA

Continuous interpolation error computation

Embedded continuous meshes. We consider the set of continuous embedded meshes M(α) = (M α (x)) x∈Ω1 . M(α) is defined on the square domain Ω 1 = [0, 1] × [0, 1] and is given by:

M α (x, y) = α h -2 1 (x, y) 0 0 h -2 2 (x, y)
, where h 1 (x, y) = 0.1(x + 1) + 0.05(x -1) and h 2 (x, y) = 0.2. The parameter α is used to control the level of accuracy of the mesh. The continuous mesh becomes coarser when α decreases but anisotropic quotients and orientation remain constant. This trend is given by the computation of the complexity C(M(α)):

C(M(α)) = N (α) = Ω1 1 h 1 h 2 (x, y) dxdy = 200 3 ln(2)α.
The parameter α defines embedded continuous meshes accordingly to Definition 3. The continuous interpolation error on M(α) is computed for two analytical functions u 1 and u 2 :

• u 1 is a quadratic function given by:

u 1 (x, y) = 6x 2 + 2xy + 4y 2 ,
• u 2 is a non quadratic function given by: u 2 (x, y) = e (2x 2 +y) .

As regards the function u 1 , the point-wise continuous interpolation error on M(α) is:

(u 1 -π M u 1 )(x, y) = 1 8 trace(M -1 2 α (x, y) |H u1 |(x, y) M -1 2 α (x, y)) = 3 (0.15x + 0.05) 2 2 α + 0.04 α = 27 x 2 + 18 x + 35 800 α .
The previous expression is then integrated over Ω 1 :

Ω1 |u 1 -π M u 1 |(x, y) dxdy = 53 800 α = 53 21 ln(2) N (α) .
For function u 2 , the point-wise continuous interpolation error on M(α) is:

(u 2 -π M u 2 )(x, y) = 1 8 trace(M -1 2 α (x, y) |H u2 |(x, y) M -1 2 α (x, y)) = e 4x 2 +y 8 α (0.15x + 0.05) 2 (4 + 16x 2 ) + 0.05 .
By a direct integration over Ω 1 , it comes:

Ω1 |u 2 -π M u 2 |(x, y) dxdy ≈ 0.2050950191 α ≈ 13.673 ln(2) N (α) .
Before comparing the continuous evaluation of interpolation error to the discrete one, some remarks can be made. According to the continuous mesh M(α), the continuous interpolation for a smooth function u writes:

u -π M u L 1 (Ω1) = C u N (α) ,
where the constant C u depends on u and N (α) is the complexity of M(α). The previous expression gives a quantitative information on the order of convergence of the interpolation error on a sequence of continuous embedded meshes issued from M(1). Indeed, with a simple analogy with uniform meshes, we have N (α) = O(h -2 (α)), so that:

u -π M u L 1 (Ω1) = C u h 2 (α).
Consequently, the continuous interpolation error model predicts an order of convergence of two on a sequence of embedded continuous meshes. From the discrete view, it is well known that a uniform refinement leads to a second order of convergence for the linear interpolation error with respect to a smooth function. This fact is also given by the continuous analysis.

Non-embedded continuous mesh. Now, let us give a more complex example in order to demonstrate the ability of the continuous mesh model to predict the order of convergence.

In this example, a discrete study of the prediction of the interpolation error is impossible whereas a clear convergence order is exhibited with the continuous mesh model.

We consider the following set of continuous meshes I(α) = (M α (x)) x∈Ω2 . I(α) is defined on the square domain Ω 2 = [0.5, 1] × [0.5, 1] and is given by:

M α (x, y) = F (x, y) α 2 h -2 1 (x, y) 0 0 α h -2 2 (x, y) t F (x, y),
where

F (x, y) = 1 x 2 + y 2 x -y y x , and h -2 1 (x, y) = 4(x 2 + y 2 ) and h 2 (x, y) -2 = 1 2 x 2 + y 2 .
Note that this set is no more embedded accordingly to Definition 3. On the contrary, the continuous meshes are rather non uniformly embedded as one size is scaled by α and the other is scaled by α 2 . The equivalent discrete refinement process is no more homogeneous,

INRIA

i.e., the factor of division of each edge while increasing α depends on the edge coordinates. Consequently, the order of convergence seems unpredictable a priori contrary to the embedded case. However, we show that we are able to predict the asymptotic convergence order of the continuous interpolation error for the set spanned by I(α) by using the continuous analysis. The complexity of I(α) is given by:

C(I(α)) = N (α) = Ω2 1 h 1 h 2 (x, y) dxdy ≈ 0.364 α 3 2 .
We consider the interpolation error of the quadratic function u 3 :

u 3 (x, y) = x 2 + y 2 .
The point-wise continuous interpolation error on I(α) is given by:

(u 3 -π M u 3 )(x, y) = 1 8 trace(M -1 2 α (x, y) |H u3 |(x, y) M -1 2 α (x, y)) = 1 16 (x 2 + y 2 ) α 2 + x 2 + y 2 2 α .
The previous expression is then integrated over Ω 2 , it results:

Ω2 |u 3 -π M u 3 |(x, y) dxdy = 0.133 N (α) 4 3 + 0.014 N (α) 2 3
.

The inhomogeneity in the scaling of the sizes leads to two terms with different order of convergence: 8 3 and 4 3 for the first and the second terms, respectively. Consequently, the asymptotic order of convergence for the continuous interpolation error of u 3 on I(α) is only O(N -2

3 ) leading to a convergence order of 4 3 ≈ 1.33. This is less than the second order reached on the previous set of meshes defined by M(α). Indeed, for a sufficiently large value of N , the term of order O(N -4

3 ) becomes negligible with respect to the low order term O(N -2 3 ). However, this approximation is only true asymptotically. Practically, the complexity allowing this simplification depends on the constant 0.014 and 0.133. According to Figure 10 (left), as soon as the complexity becomes greater than 1000, the asymptotic order of convergence is fully represented by 0.014

N (α) 2 3
. Note that this value is reachable in practice on discrete meshes. Remark 4. Depending on the constants involved in the estimation of the order of convergence, a convergence order different from the asymptotic one can be observed. If we suppose that the error for a function u on I(α) now writes:

Ω2 |u -π M u|(x, y) dxdy = 133 N (α) 4 3
+ 0.014

N (α) 2 3
.

Then, the asymptotic complexity such that the term of order 4 3 becomes negligible is difficult to reach in practice with discrete meshes (greater than 10 6 ), see Figure 10 3 ) as predominant term in the reachable range of complexity [1, 10 6 ] whereas the asymptotic convergence order is O(N -2 3 ).

of the order of O(N -4 3 ). This example shows that predicting both the order of convergence and the magnitude of the constants is crucial to get a reliable asymptotic prediction of the interpolation error.

In the next section, the analytic evaluation of the constant C u along with the convergence order are compared to discrete estimations obtained by generating unit discrete meshes with respect to M(α) and I(α) for different values of α.

Numerical interpolation error computation

Sequences of unit meshes are generated with respect to M(α) and I(α). Interpolation errors in L 1 norm are first computed on these meshes and then are compared to continuous error estimations.

Unit meshes. To validate the previous continuous evaluation of the interpolation error, a series of discrete unit meshes with respect to M(α) = (M α (x)) x∈Ω is generated. These meshes are considered for α = {1, 2, 4, 8, 16, 32}. We denote by (H α ) α∈ [1...32] this sequence of discrete meshes. They have been generated using Yams [START_REF] Frey | Yams, a fully automatic adaptive isotropic surface remeshing procedure[END_REF]. Figure 11 depicts all these meshes. The histograms reporting the meshes edges length and elements quality are given in Table 2. These histograms point out the gap between the generated unit meshes and a perfect unit mesh. We notice that an almost perfect quality is reached for each mesh and that more than 80% of the edges length lie in the range [ meshes sequence for the continuous mesh I(α) are depicted in Figure 12. Similar conclusions arise.

0.20 < L < 0.
Complexity vs. number of vertices. We first study the correlation between the discrete number of vertices N v and the continuous complexity N . This preliminary study is necessary in the continuous discrete comparison. Indeed, the continuous interpolation error estimate involves the continuous complexity whereas the discrete number of points is used in the discrete error computation. As regards the sequence of discrete meshes (H α ) α with respect to (M(α)) α , the discrete number of nodes N v is plotted as a function of the complexity N in Figure 13 (top right). This function is linear and the slope evaluated numerically gives:

∀H α , N v (α) = 1.54 N (α).
Note that this constant handles the discrepancy between the perfect unit mesh and the generated constrained discrete mesh with respect to the continuous one. Constraints arise from the domain boundary, the used mesh generator and the smoothness of the continuous mesh. This constant is in the theoretical framework equal to 2, see Theorem 2. When the constant is equal to 2, it ensures the following equality between discrete and continuous errors:

C u N = 2 C u N v .
In the case of the set of inhomogeneous continuous meshes defined by I(α), the constant is equal to 3.37.

Discrete interpolation error. The discrete interpolation error of a function u on a mesh H α reads: where K stands for an element of H α . The element-wise interpolation error u -Π h u L 1 (K) is computed by means of a 5-order Gauss quadrature numerical integration. The discrete interpolation error is compared to the continuous one by considering N v instead of N in the continuous estimates. Figure 13 plots discrete and continuous interpolation errors for functions u 1 and u 2 on M, and for function u 3 on I. An excellent correlation is obtained between the continuous and the discrete evaluations. The slight differences between the continuous and the discrete expressions come from the practical difficulty to generate a perfect unit mesh. Indeed, mesh generators generate unit meshes in the sense of Definition 5 implying that elements are quasi-unit. In consequence, most of the edges length lie in the range [

u -Π h u L 1 (Hα) = K∈Hα u -Π h u L 1 (K) ,
√ 2 2 , √ 2
]. Bounds of the interpolation error can be deduced from this edge length tolerance. For every function u and for every continuous mesh, we have :

u -π Mmin u L 1 ≤ u -π M u L 1 ≤ u -π Mmax u L 1 ,
where M min = 2 M and M max = 1 2 M. Unit meshes composed of edges length that lie in the acceptable range will verify:

INRIA

1 2 u -π M u L 1 ≤ u -Π h u L 1 ≤ 2 u -π M u L 1 .
These bounds are plotted in Figure 13. Note that previous bounds impact only the constant evaluation C u and not the order of convergence of the interpolation error.

3D examples

In this section, a 3D analytical and a 3D numerical examples are presented.

A fully analytical example. We first consider the set of 3D continuous embedded meshes ] which are given by: where h 1 (x, y, z) = 0.1(x + 1) + 0.05(x -1), h 2 (x, y, z) = 0.2 and h 3 (x, y, z) = 0.2(z + 2). The parameter α is used to control the level of accuracy of the mesh. The continuous mesh becomes coarser when α decreases. This trend is given by the computation of the complexity C(M(α)):

M(α) = (M α (x)) x∈Ω3 defined on the domain Ω 3 = [0, 1] × [0, 1] × [0, 1
M α (x, y, z) = α   h -2 1 (x, y, z) 0 0 0 h -2 2 (x, y, z) 0 0 0 h -2 3 (x, y, z)   , RR n°6846
C(M(α)) = N (α) = Ω3 1 h 1 h 2 h 3 (x, y, z) dxdydz = 1000 3 ln(2)(ln(3) -ln(2))α
We consider the interpolation error of the function u 4 given by:

u 4 (x, y, z) = e 2x+y+z .
The point-wise continuous linear interpolation error is

(u 4 -π M u 4 )(x, y, z) = 1 10 5
e 2 x+y+z 441 x 2 + 798 x + 2361 + 400 z 2 + 1600 z α .

Integration over Ω 3 leads to:

Z Z Z Ω 3 |u4 -πMu4|(x,
N (α) 2 3
.

Sequence of unit meshes with respect to M(α) for α ∈ {2, 4, 8, 16, 32} have been generated using Yams [START_REF] Frey | Yams, a fully automatic adaptive isotropic surface remeshing procedure[END_REF] and Gamanic [START_REF] George | Gamhic3d, a fully automatic adaptive mesh generation method in three dimensions[END_REF], see Figure 14. In this example, the constant linking the number of vertices and the continuous complexity is 2.07, see Figure 15 (left). Consequently, the perfect theoretical case is almost reached. The comparison between discrete and continuous interpolation error is depicted in Figure 15 (right). We observe an excellent correlation between the continuous model and the numerical computation.

A fully numerical example. The theory is now applied to a more realistic 3D example where the continuous interpolation error results from numerical computations and nor from an analytical evaluation which is not possible in that case. More precisely, the function and the continuous mesh (i.e.,the metric field) are linearly interpolated on a discrete mesh. Then, the continuous interpolation error is evaluated by using these discrete data. This case is more representative of the practical use of the continuous mesh model. The aim of this example is to point out that, even with numerical computations, the continuous mesh model turns out to be a reliable model to predict interpolation error.

If u is the function of interest and I(α) = (M α (x)) x∈Ω is the continuous mesh, the continuous interpolation error is computed on the discrete mesh H 0 as follow:

Ω |u(x) -π M u(x)| dx ≈ 1 10 K∈H0 trace(M α (K) -1 2 |H u (K)|M α (K) -1 2 ) |K|, (12)
where M α (K) and |H u (K)| are the linear discrete representation on the element K of the continuous mesh and of the Hessian of u, and where a simple 1-order Gauss quadrature numerical integration has been considered.

We first define the continuous mesh

I(α) = (M α (x)) on the domain Ω = [-1, 1] × [-1, 1] × [-1, 1] ⊂ R 3 .
The analytical expression of M α is inherited from the following quadratic form: q(x, y, z) = 2x 2yx + 6y 2 + z 2 . by F = (n, t 1 , t 2 ), where the vectors t 1 , t 2 lies in the plane orthogonal to n = g/ g 2 . We assume that (t 1 , t 2 ) is an orthonormal basis of R 2 . We also consider the Hessian matrix H q of q:

H q =   4 -1 0 -1 12 0 0 0 2   .
From this matrix, we define the sub-matrix:

H q = t t 1 H q t 1 t t 1 H q t 2 t t 1 H q t 2 t t 2 H q t 2 , INRIA 10
4 -Π h u 4 L 1 (Ω3) .
which is the projection of H q onto the plane defined by (t 1 , t 2 ). We denote by λ 1 and λ 2 the eigenvalues of H q . The continuous mesh I(α) parametrized by α is then given by:

M α = F       α 2 g 2 2 α |λ 1 | 2 g 2 α |λ 2 | 2 g 2       t F.
The complexity of I(α) is:

C(I(α)) = α 2 2 Ω g 2 2 |λ 1 | |λ 2 | ,
which can be evaluated numerically. This continuous mesh has the physical feature to follow the iso-surfaces of the function q and to adapt the sizes to the local curvature of the isosurfaces. The iso-surfaces of q are represented in Figure 16 (left). A sequence of discrete meshes (H α ) α for a complexity of [4 000, 8 000, 16 000, 32 000, 64 000, 128 000] ,

have been generated. Resulting unit meshes are depicted in Figures 17 and18. In this case, the number of vertices is almost 2.3 times the continuous complexity.

The interpolation error of the function

u 5 (x, y, z) = x 2 + y 2 + z 2 .
is computed on this sequence of discrete meshes. These discrete errors are compared to continuous interpolation errors evaluated on a fine uniform mesh of constant size h = 0.004 using Relation [START_REF] Dobrzynski | Anisotropic delaunay mesh adaptation for unsteady simulations[END_REF]. This size corresponds to the minimal size prescribed by the continuous mesh I(128 000). The correlation between the continuous and the discrete interpolation error is shown in Figure 16 (right). As expected, the order of convergence is well predicted but the numerical computations induces an over error estimation. Consequently, if the slope (i.e., the convergence order) of the error is well approximated, the perfect fitting due the approximated evaluations of C u results in a slightly shifted curve.

Some conclusions

These examples reveal that the interpolation error can be computed continuously without any discrete support. From this point of view, a discrete mesh is simply the projection of the continuous one. Discrepancies between continuous and discrete interpolation errors are due to projection errors. Practically, this gap depends on the mesh generator used and the difficulty to generate a unit mesh.

When the continuous interpolation error cannot be anymore evaluated analytically, we have shown that a numerical evaluation using a discrete mesh still provides accurate predictions of either the convergence order or the error magnitude. Figure 16: Left, iso-surfaces of q. Right, comparison between the continuous and the discrete interpolation error for u 5 on I(α). In that case, the continuous interpolation error is computed numerically by means of the linear approximation of the continuous mesh on a fine uniform discrete mesh.

INRIA

Perspective uses

We give some possible extensions related to the use of continuous meshes. In particular, we emphasize the available mathematical tools and their possible use with the concept of continuous mesh. These extensions are related to the wish of describing the optimal mesh for a given functional.

Optimal mesh adaptation. Given a smooth function u and a continuous mesh M, the linear continuous interpolation error reads, at least asymptotically:

u -π M u L 1 (Ω) = C u N k n ,
where n is the space dimension, k the expected order of convergence and C u a constant revealing alignment correlation between M and the Hessian of u. We have seen previously that the previous relation leads to a second order of convergence on a sequence of embedded continuous meshes having an increasing complexity. In the case of a sequence of non embedded continuous meshes, the order of convergence is related to the specific refinement strategy. Deriving an optimal mesh adaptation strategy consists in increasing the convergence order while reducing the constant C u . As regards the convergence order, it is fixed by the chosen refinement (embedded or inhomogeneous) along with the order of the interpolation. Consequently, the remaining possibility to reduce the error is to diminish the constant C u . This problem is rewritten as an optimization problem set up on the space of continuous meshes: min

M F (M) (13)
In the case of the linear interpolation error, the functional is equal to:

F (M) = Ω |u -π M u|.
To avoid the trivial solution giving a null error, a constraint is imposed on the complexity:

C(M) = N. (14)
In the following, two methods to solve ( 13)-( 14) are briefly presented. The first one is based on calculus of variations while the second uses differentiable optimization. Note that these methods are not defined on the space of discrete meshes.

Calculus of variations. This approach consists in deriving analytically a necessary condition for the optimal solution of ( 13)-( 14). For instance, the necessary Euler-Lagrange condition states that there exists λ such that: ∀δM, δF (M; δM) = λδC(M; δM). This method needs to be able to exhibit the variation of the functionals F and C with respect to

INRIA

M. This can be done, at least formally, on the space of continuous meshes. The variation of F is:

δF (M; δM) = lim ε→0 F (M + ε δM) -F (M) ε .
where

M + ε δM = (M(x) + εδM(x)) x∈Ω .
The previous variation is well defined on the space of continuous meshes for ε ≥ 0.

Numerical optimization. In the case of a non-linear functional F , deriving an optimality condition can be unpractical. Consequently, a numerical optimization approach is advised. Such an approach is still possible on the set of continuous meshes. Let us consider a classic descent step: from a continuous mesh M i , a new continuous mesh M i+1 is found:

M i+1 = M i -αd i ,
where d i is a descent direction ensuring that F (M i ) > F (M i+1 ). It generally involves derivatives information or approximation of F with respect to M i . The difficulty in this approach is to ensure that at each step the new continuous mesh M i+1 is well defined, i.e., whatever the point x ∈ Ω, M i+1 (x) is a metric tensor. At first, it seems that a lot of additional non linear constraints are necessary to ensure this property. For instance, we may impose all principal determinants of M i+1 (x) to be strictly positive. An elegant solution to avoid additional constraints is the use the log-Euclidean framework [START_REF] Arsigny | Log-Euclidean metrics for fast and simple calculus on diffusion tensors[END_REF]. This framework ensures that a classic descent direction step will result in a valid metric tensors field.

Conclusion

In this paper, we have proposed a continuous framework to model a mesh and its elements.

We have investigated how a discrete mesh can be viewed as the discrete projection of a continuous mesh. Far to be able to give a general continuous mesh model for any discrete meshes, we only consider a model based on the notion of Riemannian metric space. This model corresponds to a wide class of existing 2D and 3D mesh generators [START_REF] Bottasso | Anisotropic mesh adaption by metric-driven optimization[END_REF][START_REF] Coupez | Génération de maillages et adaptation de maillage par optimisation locale[END_REF][START_REF] Dobrzynski | Anisotropic delaunay mesh adaptation for unsteady simulations[END_REF][START_REF] Frey | Yams, a fully automatic adaptive isotropic surface remeshing procedure[END_REF][START_REF] George | Gamhic3d, a fully automatic adaptive mesh generation method in three dimensions[END_REF][START_REF] Hecht | BAMG: bidimensional anisotropic mesh generator[END_REF][START_REF] Jones | Validation of 3D adjoint based error estimation and mesh adaptation for sonic boom prediction[END_REF][START_REF] Laug | BL2D-V2, mailleur bidimensionnel adaptatif[END_REF][START_REF] Li | 3D anisotropic mesh adaptation by mesh modification[END_REF][START_REF] Pain | Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations[END_REF][START_REF] Tam | Anisotropic mesh adaptation for 3D flows on structured and unstructured grids[END_REF]. Practically, these mesh generators compute edges length in a Riemannian metric space and aim at generating unit meshes, i.e., meshes with only unit length edges and unit volume elements with respect to the prescribed Riemannian metric field.

The continuous mesh model is mathematically given by rewriting the symmetric definite positive matrix representing the metric tensor at each point of the computational domain. Continuous mesh features as density, orientation and stretching ratio are then exhibited. Geometric invariants that connect discrete elements and continuous elements have been provided. These invariants point out all the information contained in a metric. The discrete projection is then based on the notion of unit mesh. A unit mesh with respect to a continuous mesh is provided by the use of an adaptive mesh generator. Consequently, the discrete image depends on the method used to generate the mesh. This model has then been considered to study the interpolation error. We have demonstrated that the interpolation error estimate can be formulated in this continuous framework. The main result states that the linear interpolation error on the class of unit elements is only a function of the associated continuous element and of the considered function. In other words, the notion of continuous mesh handles sufficient information to completely model the interpolation error which is a priori only discrete. The numerical examples section has illustrated that, for a given continuous mesh (known analytically), we are able to compute the interpolation error of a given function on this continuous mesh without any discrete support. Consequently, the continuous mesh model demonstrates the well-foundedness of the metric-based mesh adaptation. The numerical examples have also shown that the continuous framework is able to predict perfectly the convergence order of the interpolation error on a sequence of embedded continuous meshes.

The continuous mesh concept is both a theoretical and a practical tool. Theoretically, it can be used to predict order of convergence for complex refinement strategies. The correlation between the continuous and the discrete estimates depends on the correlation between the continuous complexity and the discrete number of nodes. This dependance is linear. Consequently, one unit mesh is sufficient practically to fit the continuous model on the discrete one. This constant contains non linear features induced by the geometry of the domain or the specific used mesh generator. The numerical examples have demonstrated that the continuous mesh concept can also be used numerically and still gives reliable results. This feature is mandatory to deal with and to predict interpolation error on real-life problem involving complex geometries. In that case, symbolic calculus and analytical integration are impractical.

The on-going work is to use this model to propose numerical and analytical mesh adaptation strategies. Indeed, as a mathematical meaning of the continuous interpolation error uπ M u for a given continuous mesh (M(x)) x∈Ω of a domain Ω and a given function u has been exhibited in this paper, it is then possible to tackle the following global optimization problem:

min M Ω |u -π M u| p 1 p
, subject to the constraint C(M) = N . This problem is simply the continuous equivalent of seeking the optimal discrete mesh minimizing globally the L p norm of the interpolation error of the function u. The main interest of adopting a continuous view is to be able to use well defined mathematical tools, such as the calculus of variations or differential optimization that are not available on the space of discrete meshes.

A Integration formula on ellipsoids

All the results provided here can be found in [START_REF] Sykora | Quantum theory and the bayesian inference problems[END_REF]. We briefly recall the main results and the way to find them.

Notations. x is a vector R n , its coordinates are denoted x = (x 1 , . . . , x n ). We recall that B In is the unit ball of the identity matrix metric of R n . It is fully described by the set of points verifying: {x = (x 1 , . . . , x n ) | x 2 ≤ 1} .

We denote by V n (R) the volume of a sphere of R n with a radius R. With this notation, it comes V n (1) = |B In |. Given an index vector p = (p 1 , . . . , p n ) and a polynomial function

E(x, p) = i (x 2 i ) pi ,
the expression of the following integral can be computed:

W n (p) = B In E(x, p) dx. ( 15
)
We propose an induction on n to compute [START_REF] Formaggia | New anisotropic a prioiri error estimate[END_REF]. We first consider the following one-to-one change of variables: (x 1 , . . . , x n ) -→ (x 1 , . . . , x n-1 , cos(θ)), 0 ≤ θ ≤ π, the Jacobian of which is | sin(θ)| = sin(θ). If we write: p = (p 1 , . . . , p n-1 ) and x = (x 1 , . . . , x n-1 ), using the previous change of variables leads to

W n (p) = π 0
Vn-1(sin(θ))

E(x, p) dx (cos 2 (θ)) pn sin(θ) dθ.

The domain V n-1 (sin(θ)) is mapped onto V n-1 (1) by simply dividing each coordinate by sin(θ). This new change of variables is given by: (x 1 , . . . , x n-1 ) -→ 1 sin(θ) (x 1 , . . . , x n-1 ).

Its Jacobian equals sin(θ) n-1 . It comes:

Vn-1(sin(θ))

E(x, p) dx = (sin(θ)) 2pi+1 (cos 2 (θ)) pn sin(θ) dθ.

RR n°6846

We finally get: Examples. Given a quadratic form f defined on R n :

W n (p) = 2 W n-1 (p)
f (x) = 1 2 t x H x + t xG + c,
where H ∈ L(R n , R n ) is a symmetric matrix, G = (g i ) i=1...n ∈ R n a vector and c ∈ R a constant, then the following integration results hold in 2D

x 2 2 ≤1

f (x) = π 1 8 trace(H) + c ,

x 2 2 ≤1 f (x) x 1 = 1 4 π g 1 , x 2 2 ≤1 f (x) x 2 = 1 4 π g 2 .
In 3D, integration formula become:

x 2 2 ≤1 f (x) = π 2

Figure 1 : 2 i

 12 Figure 1: Left, geometric interpretation of B M = Φ M (1). v i are the eigenvectors of M and h -2 i are the eigenvalues of M. Right, geometric visualization of a Riemannian metric space (M(x)) x∈[0,1]×[0,1] . At each point x of the domain, the unit ball of metric M(x) is drawn.

Figure 3 :

 3 Figure 3: 3D examples of a unit element with respect to the identity metric (left) and to an anisotropic metric tensor (right). In each case, the unit ball of the metric is drawn at each vertex of the unit element.

Figure 5 :

 5 Figure 5: Conventions used to enumerate the edges and the faces of a triangle and of a tetrahedron.

Figure 6 :

 6 Figure 6: Left, different unit elements where only the density increases from left to right. Right, the geometric interpretation of anisotropic quotients as quotients of parallelepipeds volumes.

1 AFigure 8 :

 18 Figure 8: From left to right, Sommerville tetrahedra 1, 2, 3 and 4.

Figure 9 :

 9 Figure 9: Goldberg's tetrahedra family. They are parameterized by the prescription of the lengths a and e, with b 2 = a 2 + e 2 and c 2 = 4a 2 + e 2 . Left, gathering tetrahedra together fills a prism which basis is the regular triangle of side length e.

1 2 ( 1 2

 21 a) |H(a)| M -1 2 (a) corresponds to a new quadratic form observed in the space deformed by M(a). Indeed, if we consider the change of coordinates x = M (a) x, we define a new quadratic form f :

  (right). Consequently, the observed order of convergence for acceptable complexities (C(I(α)) ∈ [1, 10 6 ]) is

Figure 10 :

 10 Figure 10: Left, convergence history obtained for the function u 3 on the set of inhomogeneous continuous meshes spans by I(α) defined on Ω 2 . Right, convergence history for the error function of Remark 4. In that case, the coefficients involved in the error expression lead to observe O(N -43 ) as predominant term in the reachable range of complexity[1, 10 6 ] whereas the asymptotic convergence order is O(N -2 3 ).

Figure 11 :

 11 Figure 11: Unit meshes with respect to M(α) on Ω 1 for α = {1, 2, 4, 8, 16, 32} from left to right and from top to bottom.

Figure 12 :

 12 Figure 12: Unit meshes with respect to I(α) on Ω 2 for α = {16, 32, 64, 128, 256, 512} from left to right and from top to bottom.

Figure 13 :

 13 Figure 13: Top left, comparison between the continuous complexity and the number of vertices for a sequence of discrete meshes (H α ) α with respect to (M(α)) α . Top right and bottom, comparison between continuous interpolation error uπ M u L 1 (Ω) and discrete interpolation error u -Π h u L 1 (Ω) for function u 1 (top right), u 2 (bottom left) on M(α) and u 3 (bottom right) on I(α). Black plain lines represent interpolation error when considering M max and M min .

Figure 14 :

 14 Figure 14: 3D discrete meshes with respect to M(α) for α = {1, 2, 4, 8, 16, 32} from left to right and from top to bottom.

Figure 15 :

 15 Figure 15: Left, comparison between the continuous complexity and the number of vertices for a sequence of unit meshes with respect to M(α) on Ω 3 . Right, comparison between the continuous interpolation error u 4π M u 4 L 1 (Ω3) and the discrete interpolation error u 4 -Π h u 4 L 1 (Ω3) .

Figure 17 :

 17 Figure 17: A cut through the plane y = 0 of the discrete meshes with respect to I(α) of complexity 4000, 8000, 16000, 32000, 64000, 128000 from top to bottom, from left to right.

Figure 18 :

 18 Figure 18: A cut through the plane x = 0 of the discrete meshes with respect to I(α) of complexity 4000, 8000, 16000, 32000, 64000, 128000 from top to bottom, from left to right.

n- 1 i=1(

 1 sin(θ)) 2pi+1 W n-1 (p).W n (p) becomes:W n (p) = W n-1 (p)

π 2 0n 1 x i x i x j x 2 i

 22 (sin(θ)) 2α-1 (cos(θ)) 2β-1 dθ, The last integral is the Beta function. The following tables sums up examples of integrals both in 2D and 3D. These formula are used to integrate quadratic forms on V n (1):|x i | |x i x j |

Table 1 :

 1 Space filling tetrahedra characteristics. Coeff. is the coefficient that scales the tetrahedron onto a unit volume tetrahedron, i.e., |K| = √ 2/12.

	Quality

1

  He 1 t e 1 He 2 t e 1 He 3 t e 2 He 1 t e 2 He 2 t e 2 He 3 t e 3 He 1 t e 3 He 2 t e 3 He 3 He 1 ) x2 + ( t e 2 He 2 ) ŷ2 + ( t e 3 He 3 ) ẑ2 + 2( t e 1 He 2 ) xŷ + 2( t e 1 He 3 ) xẑ + 2( t e 2 He 3 ) ŷẑ ]. is now linearly interpolated on K ref . Its linear interpolate Π h u( x) writes a x + b y + c z + d, where coefficients (a, b, c, d) ∈ R 4 satisfies the following linear system ensuring the P 1 exactness, i.e., Π h u(

	u
	 	  x ŷ ẑ   .
	u in K ref reads:
	u(x( x)) = [ ( t e 1 RR n°6846 1 2

  2( t e 1 H e 2 + t e 1 H e 3 + t e 2 H e 3 ) -3( t e 1 H e 1 + t e 2 H e 2 + t e 3 H e 3 ) .

  1 H e 2 = t e 1 H e 1 + t e 2 H e 2 -t e 4 H e 4 , 2 t e 1 H e 3 = t e 1 H e 1 + t e 3 H e 3 -t e 5 H e 5 , 2 t e 2 H e 3 = t e 2 H e 2 + t e 3 H e 3 -t e 6 H e 6 . e 1 H e 2 + t e 1 H e 3 + t e 2 H e 3 ) -3( t e 1 H e 1 + t e 2 H e 2 + t e 3 H e 3 ) =

	6
	t e
	i=1

INRIA

We deduce:

2( t i H e i .

Table 2 :

 2 Quality 1/Q M and edges length for the unit meshes with respect to M(α) in 2D for α = {1, 2, 4, 6, 8, 16, 32}.

	50	5	1.94 %	0.20 < L < 0.50		1	0.19 %	0.20 < L < 0.50		6	0.71 %
	0.50 < L < 0.71	53	20.54 %	0.50 < L < 0.71		130	25.15 %	0.50 < L < 0.71		82	9.69 %
	0.71 < L < 0.90	116	44.96 %	0.71 < L < 0.90		204	39.46 %	0.71 < L < 0.90		352	41.61 %
	0.90 < L < 1.11	68	26.36 %	0.90 < L < 1.11		133	25.73 %	0.90 < L < 1.11		279	32.98 %
	1.11 < L < 1.41	16	6.20 %	1.11 < L < 1.41		49	9.48 %	1.11 < L < 1.41		127	15.01 %
	1 < Q < 2	160	99.38 %	1 < Q < 2		328	99.39 %	1 < Q < 2		541	99.27 %
	2 < Q < 3	1	0.62 %	2 < Q < 3		2	0.61 %	2 < Q < 3		4	0.73 %
	M1			M2 0.20 < L < 0.50	9	0.28 %	M4 0.20 < L < 0.50	8	0.12 %
	0.50 < L < 0.71	281	16.6 %	0.50 < L < 0.71		276	8.44 %	0.50 < L < 0.71		905	13.69 %
	0.71 < L < 0.90	564	33.31 %	0.71 < L < 0.90	1436	43.89 %	0.71 < L < 0.90	2491	37.67 %
	0.90 < L < 1.11	542	32.01 %	0.90 < L < 1.11	1071	32.73 %	0.90 < L < 1.11	2081	31.47 %
	1.11 < L < 1.41	306	18.07 %	1.11 < L < 1.41		480	14.67 %	1.11 < L < 1.41	1127	17.04 %
	1 < Q < 2	328	99.18 %	1 < Q < 2		328	99.77 %	1 < Q < 2		328	99.70 %
	2 < Q < 3	2	0.82 %	2 < Q < 3		2	0.23 %	2 < Q < 3		2	0.30 %
	M8			M16			M32

  y, z) dxdydz = `-2133 -1772 e 4 + 7466 e 2 -10322 e 6 + 6761 e 8 ´e-4 4. 10 5 α

	≈	0.73 α	≈	126.215

RR n°6846

and e = 1:RR n°6846

, √ 2] as soon as α ≥ 4. Unit INRIA

.INRIA

Acknowledgment. The authors are greatly indebted to M. Mehrenberger for his fruitful discussions and advice for proving geometric invariants and to G. Olivier for her careful reading and comments on this work.

INRIA