Nicolas Kowalski
email: kowalskin@ocre.cea.fr

Franck Ledoux
email: franck.ledoux@cea.fr

Pascal Frey
email: frey@ann.jussieu.fr

A PDE based

Keywords: mesh generation, finite element methods, partial differential equations, elliptic, quadrilateral, partitioning, multi-domain

Introduction

In numerous computational engineering applications, such as automobile crash simulations, structural mechanics, neutronics or fluid-structure interactions, quadrilateral meshes are preferred over triangular meshes [START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF]. Moreover, some numerical simulations can require several materials, gaz or fluids. The meshes are then referred as being multi domain.

Problem statement

The most desirable properties of quadrilateral meshes usually are:

1. to minimize the number of discrete singularities, defined here as the non 4-valent vertices (internal vertices that do not have exactly 4 neighbors);

2. to align elements along the boundary, i.e. to have several layers of quadrilaterals aligned along the boundary in a parallel fashion (e.g. boundary layers);

3. to achieve high-quality elements, i.e. to have each quadrilateral as close to a square or a rectangle as possible;

4. and to have the elements abide by a prescribed size map.

In case of multi domain simulations, meshes may have to comply wit a fifth desirable property: the mesh may be conform along the interface lines between incident domains. Then the line interface must be explicitly discretized in the mesh 3 . Generating meshes fulfilling all those properties is a difficult task where the optimization of one property may be in contradiction with the optimization of another one. Hence, a compromise between the various constraints and requirements is often required. Several numerical and algorithmic solutions have been proposed over the last decades to generate quadrilateral meshes abiding by those properties. Those methods can be categorized as either structured [START_REF] Faux | Computational geometry for design and manufacture[END_REF][START_REF] Whiteley | Two and Three-Quarter Dimensional Meshing Facilitators[END_REF] or unstructured [START_REF] Schneiders | A grid-based algorithm for the generation of hexahedral element meshes[END_REF][START_REF] Staten | Unconstrained Paving and Plastering: A New Idea for All Hexahedral Mesh Generation[END_REF][START_REF] Borouchaki | Adaptive triangular-quadrilateral mesh generation[END_REF][START_REF] Owen | Q-morph: An indirect approach to advancing front quad meshing[END_REF][START_REF] Lévy | Lp Centroidal Voronoi Tesselation and its Applications[END_REF][START_REF] Kowalski | Fun sheet matching: towards automatic block decomposition for hexahedral meshes[END_REF], depending on the importance they give to the underlying mesh structure. On the one hand, structured methods are unable to mesh most geometries and to deal with boundary discretization constraints along interface lines. On the other hand, unstructured methods usually yield meshes lacking some structure and having many undesired singularities. A good compromise between the geometrical flexibility of fully unstructured meshes and the numerical efficiency achieved on globally structured meshes can be achieved through block-structured meshes. Such meshes can be considered as unstructured arrangements of quadrilateral regions meshed in a structured way. Block-structured meshes have one main advantage over other kinds of meshes: they can theoretically adapt to any given domain while maintaining the good element quality associated with structured meshes. Moreover, block-structured meshes have two important added benefits that no unstructured mesh generation method would dispute. First of all, efficient solvers can be employed within individual regions due to their logical rectangular structure [START_REF] Lindquist | A comparison of numerical schemes on triangular and quadrilateral meshes[END_REF]. Second of all, block-structured meshes have a hierarchical memory structure, allowing processors to have dedicated memories with quick access containing the whole information they need.

However, to the best of our knowledge, no algorithm has been developed to provide a good enough partitioning of an arbitrary domain into four-sided regions. Most of the attempts to provide such an algorithm have more or less revolved around the use of the medial axis of the domain [START_REF] Tam | 2D finite element mesh generation by medial axis subdivision[END_REF]. However, the medial axis is numerically unstable, in that small changes of the domain geometry can greatly perturb the obtained partitioning. Moreover, the resulting regions are generally not four-sided, and thus are not optimal for quadrilateral meshing. Even for a domain as simple as a square, a medial axis decomposition provides triangular parts, and not the grid one would expect.

The approach we propose here is an attempt to get a block-structured partitioning that overcomes the main drawback of the previous developed techniques.

3. In the case of numerical simulations where contact and sliding effects between materials are modeled it is preferable to have duplicated vertices and edges in the contact areas.

Overview of our approach

Given a 2d domain Ω of arbitrary shape, we introduce a new automatic domain partitioning method where the geometric features of ∂Ω (see Figure 1a) are propagated in Ω by solving a PDE (see Figure 1-b). By this way, we ensure the three first properties given previously to be satisfied. Our process assumes a uniform triangulation T h of Ω, that is used to solve a PDE problem that yields a directionality field (see Figure 1-c and Section 2). Singularities of this field are then identified numerically and lines connecting singularities are deduced from the directionality field. Using these lines, our approach achieves a partition of Ω into four-sided regions (see Figure 1-c and Section 3). This algebraic approach to domain partitioning can be used on arbitrary 2d domains, as well as with domains presenting interfaces or multiple closed curves delimiting boundaries. Moreover it preserves the symmetry of Ω and minimizes the number of singularities. Once such a partitioning is defined, each region can then be easily meshed through a bilinear transfinite interpolation. The resulting block-structured meshes are usually of very high quality (see Figure 1-d and Section 4). Our method, however, cannot create meshes with a prescribed varying size map. It creates meshes with constant size elements, thus lacking the fourth properties a mesh shall exhibit. One way to overcome this drawback may to partition the domains and then to mesh each quadrilateral region using an unstructured method dealing with varying size maps [START_REF] Remacle | Blossom-quad: a non-uniform quadrilateral mesh generator using a minimum cost perfect matching algorithm[END_REF].

Creation of a directionality field

The main idea of our method consists in using a unit vector field prescribed on ∂Ω that we propagate on the whole Ω. This propagation is performed by solving a PDE with Dirichlet conditions defined on ∂Ω. Numerically, this PDE is solved by applying a P1 finite element method onto a triangulation T h of Ω with the Dirichlet conditions given on the vertices of ∂T h . Note that the interfaces between different domains are considered as being part of the geometry boundary in the following. A first intuition would be to use the outer normal vector field defined on the boundary. However, such a field is discontinuous at a corner with an angle of π 2 , whereas the desired mesh at such a corner is perfectly regular. To overcome this problem, we use the concept of directionality [START_REF] Ray | N-symmetry direction field design[END_REF][START_REF] Palacios | Rotational symmetry field design on surfaces[END_REF], which is continuous on such corners and thus better fits the expected structure of quadrilateral meshes.

Definition of the directionality of Ω

A regular vertex P of a quadrilateral mesh of Ω, that does not belong to ∂Ω, is 4-valent, and the 4 edges connected to P can be seen as two sets of opposite edges. Each of these sets, in turn, can be seen as the local discretization of a curve, with P being the intersection of two curves. Our aim being to generate quadrilateral elements, we assume that the the two directions of a cross must be orthogonal one to the other. The two tangent directions of these curves at P can thus be described by a cross, that we call the directionality of the mesh at P (see Figure 2). This notion of cross field is often used in computer visualization nowadays [START_REF] Ray | N-symmetry direction field design[END_REF], and is usually defined as follows : An important point to notice is that the cross field is not related to the reference axis used to compute angles. Indeed, choosing as reference axis a direction that makes an angle θ 0 with the X-axis is equivalent to rotating Ω by an angle -θ 0 . In that case, when each cross rotates by an angle -θ 0 , each representation vector rotates by an identical angle of -4θ 0 , thus leaving the problem globally unchanged.

C θ = v k = cos(θ + kπ 2 ), sin(θ + kπ 2 ) T , 0 ≤ k ≤ 3 , with θ ∈ [0, π 2 [. (1)

Definition of the representation vector field as the solution of a non-linear PDE

Our goal is to propagate an information known on the boundary of the domain inside the domain in a smooth fashion. The solution of this problem can be considered as the steady-state of a conduction problem

∂u ∂t -div(k∇u) = 0 (2)
With k constant, we get the Laplace equation. Solving this problem consists here in minimizing a functional J defined as:

J(u) = Ω |∇u| 2 dx, u(x) = u 0 (x) ∀x ∈ ∂Ω, (3)
where x is a vector function and the Dirichlet condition is defined at a point P of ∂Ω as the representation vector corresponding to the cross defined by tangents and normals to ∂Ω at P . This guarantees that the cross field is locally "aligned" with each boundary in its vicinity. To initialize this vector at every C 0 corner C of ∂Ω, we arbitrary define it at C as the average of the representation vectors of the two geometric edges incident to C. Solving this problem is easy using FEM or other numerical methods. It is an elliptic problem, well-posed with good properties as ensuring the maximum principle on u. However, in order to keep vectors unitary in the context of this work, we add a non-linear constraint. Indeed, allowing norms to vary introduce a bias leading to non-constant variations of the crosses orientations, which is not desirable as only directions are meaningful in our problem. Thus, given a geometrical domain Ω, we minimize the functional J defined as4 :

   J(u) = Ω |∇u| 2 dx, u(x) = u 0 (x) ∀x ∈ ∂Ω, |u(x)| = 1 ∀x ∈ Ω. (4)
As such, unlike Problem 2, Problem 4 is ill-posed, as a solution may not exists. Indeed, the solution of the same problem without considering the nonlinear constraint may result in some cases in a vector field having one or more zeroes; but as the non-linear constraint prevents the vector field from having any zero, the resulting vector field is not defined at some points. Instead of solving Problem (4), we find an approximate solution defined on Ω deprived of a finite number of points by proceeding in two steps:

-first, we solve problem (4) while neglecting the non-linear constraint; -then, this approximation is refined by linearizing of the norm constraint.

Solving the linear approximation of the problem

By neglecting the non-linear constraint in the third equation, Problem (4) becomes a simple steady-state heat problem. Classically, we introduce a weak formulation of this problem:

∀v ∈ V, Ω ∇ũ∇v dx = - Ω ∇u 0 ∇v dx (5)
with V = H 1 0 (Ω) being the space of H 1 functions vanishing on ∂Ω, ũ ∈ V and u = ũ + u 0 . The Lax-Milgram theorem guarantees the existence and the uniqueness of a solution to this problem.

Numerically, we solve this problem using a Galerkin finite element approach on a triangulation T h of Ω with P 1 -Lagrange elements. This yields a linear problem A × x = b. Since matrix A is symmetric and positive-definite, the vector x is computed using a conjugate gradient method [START_REF] Saad | Iterative methods for sparse linear systems[END_REF].

Linearization of the norm constraint

The non-linear constraint |u| = 1 is then considered. Due to the nonlinearity of this constraint, functional J is no longer convex and its argmin is zero, and, as it has a lower bound of zero, it means that it has several local energy minima. In our approach, we find the local minima closest to our first approximated solution using a linear approximation of the norm constraint [START_REF] Alouges | A new algorithm for computing liquid crystal stable configurations: the harmonic mapping case[END_REF] based on Lagrange multipliers. The algorithm is based onto the fact that normalizing vectors with a norm greater than one reduces the value of the functional J. In practice, we use an iterative process starting with the solution of Problem 3 as initial solution. Each step of this process is divided in two parts, a PDE resolution and a normalization process, each of them decreasing the value of J. In the following, the term v n (p) denotes the value of the vector v at point p during the step n. Considering v as a solution to Problem 5, we have |v n (p)| = 1 at each point p of the triangulation T h . In order to compute the solution at step n + 1, we add a linear constraint to each point p ∈ T h :

v n+1 (p) • v n (p) = 1. (6)
and, using the conjugate residuals method, we solve the previous steady-state heat problem 3 under constraints given by Equation 6at each point p ∈ T h . We introduce Lagrange multipliers [START_REF] Bramble | The Lagrange multiplier method for Dirichlet's problem[END_REF] to add the linear constraint at each point leading to the following problem:

A C C 0 u λ = b 1 , (7)
with 1 the column matrix having all its elements equal to 1, λ the column matrix, which coefficients λ p , p ∈ T h , corresponds to the Lagrange multiplier associated to the constraint |v(p)| = 1, and C a diagonal matrix, such that C(2p, 2p) = α and C(2p + 1, 2p + 1) = β, for v n (p) = (α, β). Adding those constraints changes the initial problem into an optimization problem under constraint. Eventually, we normalize all obtained vectors, thus decreasing the value of the functional J as the constraint ensures us that |v n+1 (p)| ≥ 1.

Remark 1. Except the traditional numerical errors that are due to the discretization, the triangulation T h that we use to solve the Problem 4 has no significative influence on the structure of the generated directionality field with regard to the convergence and the stability of the procedure.

Domains partitioning

At completion of the first step, we obtained a piecewise affine, continuous unit vector field defined at the vertices of T h . It can be extended by interpolation to define a continuous vector field over Ω. This vector field satisfies a prescribed Dirichlet boundary condition on ∂T h . Each of these representation vectors corresponds to a cross, leading to a smooth cross field CF defined over all T h . The general idea of our approach consists in analyzing the topology of this cross field. This means identifying its singularities, or zeroes, and connecting them, thus emphasizing its underlying structure. We will also consider the geometric singularities of T h , which are the C 0 corners of the domains, so as to obtain a partition of T h into quadrilateral regions.

The notion of topological skeleton of a vector field is well-known [START_REF] Helman | Vizualizing Vector Field Topology in Fluid Flows[END_REF], and can be understood as a set of curves representing the defining features of the vector flow, for example, the set of curves integrated from saddles points along the directions of the eigenvectors. We extend it here to cross fields, in a similar manner to what has been suggested in other works [START_REF] Ray | N-symmetry direction field design[END_REF]. A semantic difference is that we use this topological structure to partition the geometric domain.

Singularities of a directionality field

Definition 1. Given a vector field v defined on Ω, a point x ∈ Ω such that v(x) = 0 is a singularity. By extension, in our context, the singularities of the directionality field are detected as the zeroes of the piecewise linear interpolation of the representation vector field defined at the vertices of T h . Moreover, we include the discontinuities of the domain boundaries as singularities. Considering these intrinsic geometric singularities allows us to obtain a fully quadrilateral partitioning of the domain (see Section 3.3). We recall here some definitions and propositions [START_REF] Tricoche | Continuous Topology Simplification of Planar Vector Fields[END_REF] related to the Poincaré index of the zeroes of a vector field. Definition 2. The Poincaré index of an isolated zero x of a 2d vector field v = (v 1 , v 2 ) is the number of rotations of the vector field while traveling in positive direction along any closed curve enclosing x but no other singularity. Put another way, for γ a closed curve containing x and no other singularity, the Poincaré index i x of P 0 is defined as:

i x = 1 2π γ dφ, with φ = arctan v 1 v 2 . ( 8
)
Proposition 1. Let v be a vector field on a 2d geometrical domain Ω. Suppose v is defined as a linear interpolation of non-zero values given at the vertices of a triangulation, T h of Ω, all singularities are necessarily of first order, i.e. their index is either 1 or -1.

Proof: The result is exhibited by enumerating the different possibilities and summing the angular rotations of the field along the edges of the triangle around the singularity.

Streamlines and separatrices

A streamline of a vector field is a special kind of field line that can be defined as follows: Definition 3. Let v be a C 1 -continuous vector field defined over a domain Ω, and let γ(s), s ∈ [0, 1] be the parametrization of a C 1 curve defined over Ω, with ∂γ(s) its derivative. Then γ is a streamline of v if

∀ s ∈ [0, 1], ∂γ(s) × v(γ(s)) = 0, (9)
where "×" denotes here the cross product of two vectors. A streamline is thus a C 1 -continuous curve for which the tangent to the curve at any given point has the same direction as the vector field value at this point. It is to be noted that there is exactly one streamline going through any regular point of a vector field, and that if γ is not a closed curve then both end points of γ are either a singularity of the representation vector field or lie on ∂Ω. Now we extend this concept to cross fields. Around any regular point of a cross field, this field can locally be described as 4 vector fields. At singular points, such a separation of the components of the crosses is impossible. Thus, in regular areas (containing no singularity), locally to a small neighborhood of each point, it can be described as the orthogonal intersection of two streamlines. Those streamlines can be computed, starting from this point, by selecting one of the 4 vectors defined at this point and, doing so, by selecting one of the 4 smooth vector fields and then using the same numerical scheme one would use to compute a vector streamline. By repeating this process, one can compute whole streamlines on the field as long as they do not go in the vicinity of field singularities. We can now give a suitable definition for streamlines in cross fields. Definition 4. Let CF = (v 1 , v 2 , v 3 , v 4 ) be a continuous cross field defined over a domain Ω, except for a finite number of singularities 5 . Let γ(s), s ∈ [0, 1] be a parametrization of a C 1 curve defined over Ω, with ∂γ(s) its derivative. Then γ is a streamline of CF if

∀ s ∈ [0, 1], ∃i ∈ [1..4], ∂γ(s) × v i (γ(s)) = 0. ( 10
)
For each singularity S, we build the separatrices, which are the streamlines of the cross field ending on S (the red streamlines in Figure 3-a, for example). Definition 5. A streamline γ is a separatrix if there is a singularity S 0 such that ∃ t ∈ [0, 1], γ(t) = S 0 . In that case, we say that γ is a separatrix of S 0 . a b Fig. 3. In a and b, streamlines in a triangle containing a singularity or not.

To compute the separatrices at a singularity S 0 , we have to find a streamline γ such that γ(t) = S 0 for t ∈ [0, 1]. Numerically, in the triangulation T h of Ω, we apply the following process. As a singularity S 0 is necessarily located in a triangle T 0 of T h or on an edge separating two triangles T 1 and T 2 of T h , any separatrix of S 0 crosses the edges of T 0 in the first case and the edges of T 1 and T 2 in the latter case. Let us consider the first case, the latter case being a straightforward extension. Considering a singularity S 0 located into a triangle T 0 ∈ T h defined by its 3 vertices S 1 , S 2 and S 3 . We traverse all the edges of T 0 to find the points where separatrices of S 0 intersect the edges of T 0 . Let P be a point of the edge [S1, S2]. Considering the cross field CF , we compute the value CF (P ) by linearly interpolating the value of CF (S 1 ) and CF (S 2 ). The value of CF (P ) corresponds to a cross or to 4 vectors v k , 0 ≤ k ≤ 3 representing the directions of the streamlines at P . Let u(P ) be the vector S 0 P. Hence, a streamline going through P is a separatrix of S 0 iff 5. we remind the reader here that this means that the representation vector field of F is C 1 -continuous, and not that any of the vi field is

C 1 -continuous ∃k ∈ [0..3], v k |v k | = u(P ) |u(P )| . ( 11
)
In this case, the tangent of the streamline at P is v k (see Figure 4). Applying this process for every singularities of CF , we define all the intersection points between the triangles containing a singularity and the associated separatrices. All these points are on the boundary of the triangulation T ′ h = T h -{T 0 , . . . , T n } with {T 0 , . . . , T n } the sets of triangles containing the singularities of CF . By construction, CF is regular on T ′ h . Then the expected separatrices can be approximated on T ′ h using the following Runge-Kutta numerical scheme until either the boundary of the geometry or another triangle containing a singularity is met. Considering a triangle T ∈ T ′ h , a point X i ∈ ∂T and a direction d i at X i , our aim is to find the intersection point X i+1 between ∂T and the line defined by the direction d i and the point X i . Hence the direction d i+1 in X i+1 should be computed (see Figure 5-a and 5-b).

To this end, a representation vector in X i is linearly interpolated from the representation vectors of the vertices S i of T . The vector v di of the corresponding cross C i that is closest to the input direction d i is selected. The selection of this vector is propagated onto all crosses interpolated over T , thus defining a vector field Y = f (X) over T that is consistent with d i (see Figure 5-a). Then we use the Heun's method [START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF], an easy-to-use variation of a Runge-Kutta method, to integrate the ODE Y = f (X) over T , starting from X i . A line can thus be built using segments, by iterating this integration process over the successive encountered triangles. It is to be noted that, even though the singularities of the representation vector field and the ones of the cross field are identical, their separatrices are totally different.

Another interesting point is that the Poincaré index of a singularity S of the representation vector field and the number of separatrices of S are closely related. Each singularity of the vector field having an index of +1 corresponds

Ci di v d i = f (Xi) Xi T S1 S2 S3 f (S1) f (S2) f (S3) Xi X ′ i+1 f (Xi) di+1 f (X ′ i+1 ) Xi+1 a b
Fig. 5. Illustration of the integration process over a triangle.

to a singularity of the cross field having 3 separatrices, and each singularity of the vector field having an index of -1 corresponds to a singularity of the cross field having 5 separatrices, while regular points of the vector field correspond to points where a cross is defined and therefore are the intersection of two streamlines in the cross field (see Figure 6).

Definition of the domain partitioning.

Considering a domain Ω and a cross field defined over Ω, the separatrices of this cross field along with their corresponding singularities allow us to partition Ω. Every separatrix ends either on a singularity or on ∂Ω, there is no guarantee that the obtained regions are four-sided. In fact, it only provides regions that can be meshed through a submapping technique [START_REF] White | Automated Hexahedral Mesh Generation by Virtual Decomposition[END_REF], i.e. regions without any inner singularities. In order to get four-sided regions, we extend the set of of representation vector field singularities by adding a set of geometric singularities, which are defined as the C 0 corners on ∂Ω. We then build separatrices for those geometric singularities by using the same process than for the representation vector field singularities. Proposition 2. The set of separatrices of the geometric and field singularities leads to a partition of Ω into 3 or 4-sided regions. The 3-sided regions necessarily contain a singularity.

Proof: By construction, all the regions we obtained contains no singularity except eventually on their corners. This means that a parametrization f (u, v) of the inside of a region is always possible. Moreover, the boundaries of the considered region can always be included inside the parametrization by continuity, except if one edge of the region is of null length, which corresponds to a triangular region. In this case, one of the two sets of parallel streamlines of the region has all its streamlines that converge to a same point, which is the degenerated edge of the region.

v1 v2 v3 v1 v2 v3 ∆12 > 0 ∆23 > 0 ∆31 > 0 ∆12 + ∆23 + ∆31 = 2π v1 v2 v3 v1 v2 v3 ∆12 < 0 ∆23 < 0 ∆31 < 0 ∆12 + ∆23 + ∆31 = -2π
Fig. 6. A 3-valent singularity, on top, and a 5-valent singularity, on bottom. In blue, the local cross field, and in green the corresponding representation vector field. Separatrices are in red, along with the shape of surrounding streamlines in blue dots. On the right, in light green dots, the linear interpolation of the representation vectors along the edges is shown. The path along which we compute the index is the triangle around the singularity.

For all the streamlines to behave like this, the angle at this corner of the region must be at maximum π 4 , by continuity of the cross field inside the region. As the directions of the separatrices going from an inner singularity are uniformly distributed, due to the piecewise linear nature of the cross field, and as we only have 3 and 5-valent singularities, the minimum angle for a region corner that is a singularity inside Ω is 2π 5 . It follows that the degenerated edge of the considered region is a geometric corner featuring a small angle. In the absence of such a degenerated edge, the region can be parametrized by a one to one mapping function f (u, v), with u ∈ [0, 1] and v ∈ [0, 1].⊓ ⊔ Note that, while the computational cost of the method is strictly related to the number of elements in the background triangulation, the time needed to partition the domain for the different examples presented in this paper varied from about one second, for most of them, to twenty seconds for the most complicated ones.

Results

Once each of the domain has been partitioned into four-sided parts, each of these parts is meshed using a transfinite bilinear interpolation. This meshing process creates a parametrization of the region that coincides with the boundary parametrization. It leads to an immediate grid meshing process by considering any discretization of the parameter spaces of u and v. It shall be noted that the computational cost of this step, even for millions of elements, is negligible with respect to the cost of the domains partitioning step.

Quality comparison

While it is not easy to define a quality measurement for a domain partitioning, it is a lot more convenient to define one for quadrilateral meshes. In this section, we compare meshes obtained with our approach with meshes obtained using an unstructured algorithm provided by the software GMSH [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. Both the minimum angle measure and the scaled Jacobian measure are used to compare the results. Both meshes have about 9000 quadrilaterals. On Figure 7, we present the minimum angle measure as color maps for each mesh. On Figure 8, we do the same for the scaled Jacobian measure and, on Figure 9, we presents histograms of those results. For both measures, we get better results with -our approach. It is mainly due to the weak number of singularities that arises in the partitioning.

Single domain examples

The first example of Figure 10 illustrate the symmetry preservation of our approach. We can also note the weak number of singularities: four 3-valent singularities and three 4-valent ones.

In Figure 11, we show how the domain partitioning evolves when the geometry of the domain is slightly modified. In this case, the domain in a rectangular area where two quarter of circles have been removed. In Fig. 11-a the two quarters of circles are symmetric and then we obtain a symmetric quadrilateral mesh. In Fig. 11-b, the vertical side of the right quarter of circle is greater and the two singularities that are inside the domain are not more connected. In Fig. 11-c, the difference between the two quarters of circle is more important but the global topology of the domain partitioning is similar.

The last example of Figures 12 and 13 is a two-dimensional aerofoil section consisting in three elements, a leading-edge flat, a main element and a trailingedge flap. In order to compute flows around these three elements, it is better to get boundary layers of quadrilateral cells along the the aerofoil. It is what we get in Fig. 12-b where just few singularities appear around the aerofoil. In Figure 13, we highlight two different regions of the final mesh. The boundary layers we generate allows to get cells orthogonal to the boundary.

A multi-domain example

In Figure 14, two incident domains are partitioned separately. Even if the obtained meshes have a similar element size, the common curve is not discretized in the same way (see Fig. 14-b). We can also see that the partitioning lines of each domain do not cross the common curve (see Fig. 14-a). On the contrary, in Figure 15, the two domains were partitioned and discretized together. Thus, we obtain a contact line that was discretized only once and a conformal mesh all over both domains. It may be noted that unlike in Figure 14, the separatrices emanating from a singularity of a domain can spread through the other domain.

Conclusion

In this paper, we presented a novel approach to domain partitioning, which can be used as an efficient block-structured quadrilateral mesh generation algorithm. This fully automatic algorithm uses the solution of a PDE to extend the information contained by normals and tangents at the boundary inside the domain. It generates a partition of the domain into four-sided regions with curvilinear edges that can be easily meshed into quadrilateral elements using a transfinite bilinear interpolation. The number of singular vertices in the resulting mesh is minimal or close to minimal. Being an algebraic method at its core, our method provides results that are independent of the initial orientation of the domain and that respects its geometrical symmetries. An added benefit of our method is that the main computational cost is due to the domain partitioning. As a result, the additional cost for having a more refined mesh is very low, making our method quicker than most other methods for refined meshes. Should the need arise, both the partitioning and the meshing processes have data structures and implementation allowing for parallelization. Although the meshing part of the process is currently unable to deal with varying prescribed element sizes, future work is scheduled to improve this by using known methods to refine the mesh inside desired regions.

Fig. 1 .

 1 Fig.1. Overview of the main steps of our algorithm. In a, the background mesh used. The geometry has two domains with distinct mesh colors. In b, the unit representation vector field obtained by solving a PDE. In c, the corresponding cross field and the domain partitioning into four-sided regions. In d, a full quadrilateral mesh is obtained by using a bilinear transfinite interpolation over each region.

Fig. 2 .

 2 Fig. 2. Examples of crosses along the corresponding representation vectors. Notice that the crosses in (a) and (e) have the same representation vector, and that the representation vector may be different from any vector of the cross.

Fig. 4 .

 4 Fig. 4. The intersection of a separatrix (in red) of the singularity S0 with the triangle around S0, at a point P .

Fig. 7 .

 7 Fig. 7. Color map of the minimum angle quality measure on our mesh (in a) and on the unstructured mesh (in b).

Fig. 8 .Fig. 9 .

 89 Fig. 8. Color map of the Jacobian quality measure on our mesh (in a) and on the unstructured mesh (in b) as a color map.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. A symmetric domain partitioning on the left and two focuses on the obtained quadrilateral mesh

Fig. 12 .Fig. 13 .

 1213 Fig. 12. A 2d space around an aerofoil in a and and the corresponding domain partitioning in b

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Two domains sharing a curve are discretized separately: domain partitioning in (a) and quadrilateral meshes in (b).

here, the used norm is the Froebenius norm.


Pierre Meye
email: pierre.meye@orange.com

Philippe Raipin
email: philippe.raipin@orange.com

Frédéric Tronel
email: frederic.tronel@supelec.fr

Emmanuelle Anceaume
email: anceaume@irisa.fr

Mistore: A distributed storage system leveraging the DSL infrastructure of an ISP

Internet Service Providers furnishing cloud storage services usually rely on big data centers. These centralized architectures induce many drawbacks in terms of scalability, reliability, and high access latency as data centers are single points of failure and are not necessarily located close to the users. This paper introduces Mistore, a distributed storage system aiming at guaranteeing data availability, durability, low access latency by leveraging the Digital Subscriber Line infrastructure of an ISP. Mistore uses the available storage resources of a large number of home gateways and points of presence for content storage and caching facilities reducing the role of the data center to a load balancer. Mistore also targets data consistency by providing multiple types of consistency criteria on content and a versioning system allowing users to get access to any prior versions of their contents. Mistore validation has been achieved through extensive simulations

I. INTRODUCTION

Most of the existing architectures for data storage of Internet Service Providers (ISPs) are based on very large centralized datacenters that store and manage all the information related to their clients and their data. With the ever growing number of clients and the amount of data, centralized architectures reach their limit in terms of scalability, reliability, and access latencies. These drawbacks are commonly addressed by decentralizing the system over multiple geo-distributed nodes. Data are then placed on nodes close to users which reduces data access latencies and improves the scalability and reliability of the data by eliminating the single point of failure issue. Most of the storage systems addressing these issues either rely on peer to peer (P2P) technologies where peers are responsible for storing a subset of the data, or take advantage of the presence of large datacenters. In the former case, transient peers availability requires expensive maintenance procedures to reach an acceptable overall availability, while in the latter case, clients suffer from high data access latencies due to the remoteness of datacenters from users. This has recently led to the design of hybrid approaches relying on both peers and datacenters, where peers reduce datacenters workload and data access latencies, and improve system scalability while datacenters compensate peers instability.

In this paper, we investigate the design of a hybrid distributed storage system by leveraging highly available nodes in the network infrastructure of an ISP. Specifically, our approach combines the use of datacenters with home gateways and points of presence (POPs) equipments. Let us first briefly describe the Digital Subscriber Line (DSL) of an ISP. Fig. 1 illustrates a simple but common network infrastructure providing access to Internet through a DSL technology. Each user is equipped with a home gateway that provides access to multiple services like telephone, television and Internet. Each line of subscribers is first aggregated by a Digital Subscriber Line Access Multiplexer (DSLAM) which can aggregate thousands of lines. At a second level, the subscriber lines are aggregated in a high-capacity Asynchronous Transfer Mode (ATM) or Gigabit Ethernet link from the DSLAM to the IP network of the ISP. The flow coming from a DSLAM enters the ISP network via a POP that can handle links from a large number of DSLAMs. The ISP network is directly connected to the Internet backbone, and uses large datacenters to store information related to their subscribers and their data when a storage service is provided.

Based on this classical infrastructure, we aim at leveraging the Home Gateways (HGs) to take advantage of their native resources (i.e. computing, memory, and storage) a , and to the fact that users let most of the time their HGs powered on [START_REF] Valancius | Greening the internet with nano data centers[END_REF]. Thus, exploiting resources provided by HGs should yield a large number of high available and intelligent storage nodes located very close to the users. Second, we also aim at leveraging the POPs to benefit from their natural geographic repartition and their position at the edge of the ISP network, and the fact that all the traffic within and across ISPs goes through the POPs. We plan to use the POPs to bring the intelligence of the storage system close to users, as pioneered in [START_REF] Pallis | Insight and perspectives for content delivery networks[END_REF] for content caching and distribution in Content Delivery Networks (CDNs). This paper presents Mistore, a distributed storage system dedicated to users who access Internet via a DSL technology that ensures data availability, durability, and low access latency by fully exploiting HGs storage capacities, POPs physical localizations, and datacenters. The focus of the paper is on how data consistency and replication are managed in Mistore.

a For instance the livebox play (http://liveboxplay.orange.fr/) of the french ISP Orange contains a 1.2GHz Intel Atom CE4257 processor, 2GB of RAM, and a hard drive of 320GB. Data consistency: System designers usually face the CAP theorem (tradeoffs about Consistency, Availability and network Partition tolerance) when designing a distributed storage system [START_REF] Vogels | Eventually consistent[END_REF], [START_REF] Abadi | Consistency Tradeoffs in Modern Distributed Database System Design: CAP is Only Part of the Story[END_REF]. We choose to serialize the write requests in order to avoid the complexity of conflicts resolution to application developers. On the other hand multiple consistency criteria to parametrize the read requests are provided, namely a readers/writer mutual exclusion, an atomic consistency, and an eventual monotonic-read consistency criteria. We also provide a versioning system on data. This mechanism creates a new version of a piece of data when an update (i.e. a write request) is performed on it. This mechanism naturally provides a linear history about data and their updates.

Data replication:

We ensure data availability and durability through replication. Users send their write requests to the closest POP that will perform the data replication on behalf to them. It aims to mitigate the bottleneck that users low upload bandwidths may cause by moving the replication overhead to the IP core network of the ISP where the bandwidth is much larger. To reduce latencies data are replicated synchronously on POPs and asynchronously on HGs and stripping techniques can be used to store and retrieve the data fragments.

The remainder of the paper is organized as follows. Section II presents related works. Section III presents the system model and the assumptions considered. The architecture is presented in section IV. The different consistency criteria are defined in Section V and the way we implement it as well as read and write operations is described in Section VI. Section VII presents the evaluation of our design choices and we conclude in Section VIII.

II. RELATED WORK

Architecture: Recently, some works have investigated new designs of distributed storage systems to face dynamicity issues of pure P2P systems and the cost on large data centers [START_REF] Greenberg | The cost of a cloud: research problems in data center networks[END_REF]. The hybrid approaches, sometimes called peerassisted approaches [6] are promising and are developed even for content distribution purpose [START_REF] Sun | FS2You: Peer-assisted semipersistent online storage at a large scale[END_REF]. Some works investigate the use of distributed and stable small datacenters or home gateways [START_REF] Valancius | Greening the internet with nano data centers[END_REF]. In contrast to those approaches, Mistore is dedicated to an ISP which is a clear advantage over other providers because they master all the network infrastructure and thus they can push contents closer to the users [START_REF] Huang | Understanding hybrid CDN-P2P: why limelight needs its own Red Swoosh[END_REF].

Multiple consistency criteria: It is well known that the strong consistency models have lower performance than the weak consistency models in terms of availability, latency, and scalability [START_REF] Vogels | Eventually consistent[END_REF], [START_REF] Abadi | Consistency Tradeoffs in Modern Distributed Database System Design: CAP is Only Part of the Story[END_REF], [START_REF] Saito | Optimistic replication[END_REF]. Most of cloud storage provide the eventual consistency forms [START_REF] Vogels | Eventually consistent[END_REF]. However, even if the eventual consistency is sufficient for a large panel of applications some others require strong consistency criteria [START_REF] Decandia | Dynamo: amazon's highly available key-value store[END_REF]. Moreover, a strong consistency is a desirable property for application programmers as it is easy to reason about. As a consequence, several systems target strong consistency while others aims at adding strong consistency features to provide multiple consistency criteria in their system [START_REF] Lloyd | Don't settle for eventual: scalable causal consistency for wide-area storage with COPS[END_REF], [START_REF] Patterson | Serializability, not serial: concurrency control and availability in multidatacenter datastores[END_REF]. Mistore follows this vision by giving the possibility to choose the consistency criteria adapted to the application needs.

III. SYSTEM MODEL

Mistore targets personal data (e.g., documents, pictures, videos, music, etc). In the following, these data will be called objects. Communications are assumed to be reliable, and messages between two nodes are assumed to be received in the order they have been sent. POPs are interconnected via highly available and redundant channels provided by the ISP. Network partitions at this level occur with a negligible probability as the IP core network can recover and restore its state to avoid service disruption due to links/nodes failures [START_REF] Francois | Achieving sub-second igp convergence in large ip networks[END_REF], [START_REF] Iannaccone | Feasibility of ip restoration in a tier 1 backbone[END_REF]. So we assume a system with no network partitions between the POPs. We assume different failures models for POPs, HGs, and the datacenter. The model of failure for POPs and datacenter is crash recovery, i.e., we suppose that they both recover from a consistent state by using stable storage. Conversely, the model of failures exhibited by HGs is fail stop. We assume that crash are eventually detected by an eventual perfect detector failure [START_REF] Chandra | Unreliable failure detectors for reliable distributed systems[END_REF]. We leave for future work more aggressive failures, Byzantine failures and objects corruption.

IV. MISTORE ARCHITECTURE

We consider a DSL infrastructure operated by a single ISP (See Fig. 1), whose main components are as follows.

The home gateways are connected in a client/server mode to the POP aggregating their xDSL lines and form its region. They communicate with the storage system via this POP. A fraction of the storage capacity of HGs is dedicated to the storage system. One part of this fraction stores some of the data that have been warehoused within the storage system. Note that these data do not necessarily belong to the owner of the HG. The other part of the fraction caches the data recently or frequently accessed by the HG owner. The Points of Presence provide several functionalities. They cache some data that transit through them, implement the data consistency and replication strategy of the storage system, and monitor the HGs they aggregate. Monitoring data (e.g., storage capacity, data stored, availability state, etc.) are periodically transferred to datacenters. The datacenter main role is to collect and store metadata and usage patterns from the HGs and POPs. It also offers a backup functionality for objects when their primary POP is unavailable (See Section V).

We leave for future work our vision to leverage the DCs knowledge on the system to implement multi-and auto-tiered storage features which can refer to an automated placement of data on nodes to optimize performance, availability, and recovery [START_REF] Aguilera | Improving recoverability in multi-tier storage systems[END_REF]. Indeed, we have nodes of different types and locations that could be distinguished in three tiers for storing cold, warm, and hot data while allowing data to move from one tier to another depending on data access patterns. Cold data are infrequently or never accessed, and thus could be efficiently and inexpensively handled by datacenters. Warm data are data that have been recently accessed in the storage system and present a higher probability to be accessed than cold data. They could be stored in regions from where their users access the storage system to reduce access latencies. Finally, Hot data are frequently accessed and could be stored in POPs and distributed like in a CDN. It would allow to answer users requests while avoiding to overload HGs that have low upload bandwidths or datacenters that exhibit large access latencies. Looking in this perspective, we recall that this paper focuses on the way Mistore deals data consistency and data replication on the HGs and the POPs.

V. OBJECT CONSISTENCY Mistore replicates user objects to guarantee both their durability and availability. Now, as objects can be updated, consistency issues must be taken into account. Before diving into the different consistency criteria provided by Mistore, we briefly recall the main concepts and definitions that we will use. Traditionally consistency criteria are classified according to the fact that, given a distributed execution satisfying the considered criteria, it is always possible or not to build an history of this execution that could have been executed on a single processor and produce the same visible result. When this property is satisfied, the considered consistency criteria is said to be strong, otherwise it is weak. A strong consistency criteria requires stronger synchronization between processes that participates to the computation. This impacts the performance of the protocol but make the life of application developers easier. One can expect better performance of weak consistency criteria at the cost of more complex situations to deal with when manipulating data. Mistore implements both strong and weak consistency criteria. We rapidly expose the theory behind strong consistency criteria [START_REF] Raynal | A suite of formal definitions for consistency criteria in distributed shared memories[END_REF], and give a formal treatment of the weak consistency criteria implemented by Mistore.

A. Processes and operations

We consider a concurrent system composed of a set Π of n processes denoted p 1 , p 2 , . . . , p n that cooperate through a finite set of shared objects X (e.g., files in the case of Mistore).

Each shared object x ∈ X supports two types of operations: a read operation, and a write operation. The execution of writing value v into object x by process p i is denoted by w i (x)v. Conversely, the execution of reading value v from object x by process p i is denoted r i (x)v. We may omit the identity of the process that performs the operation when it is not relevant. To simplify we assume that each value written in an object is unique. b

B. History, legality and linear extension

Let h i| W denotes the set of write operations performed by process p i . Conversely let h i| R denotes the set of read operations performed by process p i . Let h i = h i| R ∪ h i| W denotes the set of operations (both read and write) performed by process p i . We assume that on a given process p i only a single operation can occur at a given time. Hence operations in h i are naturally totally ordered by an order relation that we denote → i . We call local history of p i the set h i ordered by → i . It is denoted by

h i = (h i , → i ).

Definition 1 (Global history):

A global history H = (H, → H ) of a concurrent system (Π, X) is a set H partially ordered by → H , with H = i h i and → H a partial order relation containing → rf (that is → rf ⊆→ H ),
∈ H if op 1 → H op 2 then op 1 → op 2 .
A linear extension of a global history H can be seen as a topological sort of the directed acyclic graph c induced by the partial order relation → H .

We assume that an initial value has been written (by a fictitious write operation) in each variable. With this assumption, a linear history is legal if the following holds.

Definition 3 (Legality of a read): Let H = (H, →) be a linear extension. A read operation op r = r(x)v ∈ H is legal if and only if:

1) there exists a corresponding write operation op w = w(x)v ∈ H such that op w → op r , b This hypothesis can be easily implemented by assuming that each value written in an object is a triplet of the form of (p i , v, count i ) where count i is a counter maintained by process p i counting its write operations.

c This graph G = (V, E) is defined by its vertexes V = H and there is a edge between two vertices op 1 ∈ V and op 2 ∈ V if and only if op 1 → H op 2 .

2) and for all operations op such that op w → op → op r op = w(x)⋆.

Definition 4 (Legality of a linear extension):

A linear extension H = (H, →) is legal if and only if all of its read operations are legal.

A linear history is legal when both its read and write operations respect the expected semantics of variables on a single threaded processor: a read operation on a variable returns the last value written in this variable.

C. Strong consistency criteria

Definition 5 (Sequential consistency): A global history H = (H, → H ) is sequentially consistent if it admits a linear extension which is legal. This consistency criterion is the simplest and weakest example of strong consistency criteria. It means that the concurrent execution could have happened on a single processor. This criteria has been first proposed by Lamport in [START_REF] Lamport | Concurrent reading and writing[END_REF].

We now present a stronger consistency criteria which brings into play physical real time. Each operation op happening on a process starts at a given real time denoted by op.start and ends at a given real time denoted by op.end (we have op.start < op.end). With these notations in hand, we can specify the partial order relation denoted by → rt as follows.

Definition 6 (Real time precedence): Let H be a set of operations, and let op 1 , op 2 ∈ H be any two operations. We say that op 1 precedes op 2 with respect to real time, which is denoted by op 1) H is sequentially consistent, 2) and → rt ⊆→ (i.e → is compatible with → rt ).

Definition 9 (Read/write order): Let H be a set of operations, and let op 1 , op 2 ∈ H be any two operations. We say that op 1 precedes op 2 with respect to the read/write order, which will be denoted by op 1 → rw op 2 , if and only if:

• either op 1 → rt op 2 , • or op 1 rt op 2 , op 1 = r(x)⋆ and op 2 = w(x)⋆, then 1) op 1 → rw op 2 ⇔ op 1 .start ≤ op 2 .start, 2) op 2 → rw op 1 ⇔ op 2 .start < op 1 .start.
It is clear that → rt ⊆→ rw .

Definition 10 (Read/write mutual exclusion): A global history H = (H, → H ) is compatible with the read-write mutual exclusion consistency criteria if it admits a linear extension H = (H, →) such that → rw ⊆→.

D. Weak consistency criterion

From these criteria, we propose to define a weak consistency criteria, that we will call monotonic read consistency. To the best of our knowledge it is the first time that a weak criteria is defined in the same formalism framework as strong ones.

Definition 11 (Local history): Let H = (H, → H ) be a global history. A local history (H i , → Hi ) with respect to process p i is defined as follows :

• H i = h i ∪ j =i h j | W , • → Hi is compatible

VI. IMPLEMENTATION

Mistore implements the different levels of consistency described in Section V through a lock service. Two types of locks exist. A read lock that ensures the mutual exclusion consistency criteria, and a write lock that ensures the SWMR model. The lock service follows a primary-backup scheme [START_REF] Alsberg | A principle for resilient sharing of distributed resources[END_REF]. A primary and backup lock server per object is implemented and are executed by the POPs. Fig. 2 and Fig. 3 show the pseudocode executed by HGs and POPs to acquire a lock from respectively their associated POP and the primary POP of the concerned object. The lock on an object is created by function GRANTLOCK shown in Fig. 4. Function ISAVAILABLE checks if a requested lock on a specific object can be granted to a HG. To deal with nodes failures, a lease is associated with each granted lock [START_REF] Thekkath | Frangipani: a scalable distributed file system[END_REF], [START_REF] Gray | Leases: an efficient fault-tolerant mechanism for distributed file cache consistency[END_REF]. At lease time, that is, at the expiration time of the lock, a client has to renew the lease if the current read/write operation is not completed. This is achieved by contacting the primary in charge of the object. Note that for communication costs reasons, the lease time is only activated on the primary, not on the backup. Thus, the primary does not have to acknowledge the backup each time a lock lease is renewed. The backup activates the lease time of a lock only when it takes over this role upon detection of the primary failure.

A. Write operation

A write operation on an object creates a new version of this object in Mistore. When a client wants to write (update) an object o, it sends a request to the POP p associated to his HG as presented in Fig. 5 and Fig 6 . If o is written for the first time, a primary and a backup POP are affected to it (See Fig 7). These POPs are responsible for the replication of that object and for its consistency management. The primary is the POP associated with the client HG issuing the creation request. The backup is randomly determined by the primary among the other POPs of the system (invocation of the GETBACKUP() function in Fig. 7). The create operation returns o id the identifier of an object o, which is the concatenation of the identifiers of o primary, o backup, o HG, and a creation counter that represents the number of objects created by the HG that owns o. In the following the primary and the backup of object o are respectively denoted by o p and o b . When the primary of an object is unavailable, its backup becomes the primary and the datacenter becomes the new backup. When the primary recovers, it takes over its role. Our solution only makes use of one backup because of the high availability of these nodes [START_REF] Gill | Understanding network failures in data centers: measurement, analysis, and implications[END_REF]. A write operation on object on o that has already been created in the system must contain the write lock for o to respect the SWMR model. Note that as a client can crash or goes off before the end of a write operation, the write lock on o is maintained by POP p as long as the operation is not completed to ensure the SWMR model. POP p locally caches a copy of o and asks both the primary o p and the backup o b to replicate o. When o p and o b receive this request, they both locally cache o, acknowledge POP p (as described in Fig. 8), and trigger the replication operation on their HGs (see Lines 3-7 in Fig. 8, and Fig. 9). Upon receipt of both acknowledgements, POP p releases the write lock and notifies the client that the write operation is completed. Object o is then available for subsequent read and write operations, which makes the replication process on HGs transparent to clients.

B. Read operation

Mistore gives the opportunity to specify the consistency criteria that the read object must guarantee, or a specific version of desired object. Both choices are specified in the version flag parameter of the read operation. In addition, in the quest of reducing the latency of read operations, Mistore favors readings from local caches prior to propagating the request to both the primary and the backup of the object. Read operations handle two parameters, the identifier o id of the requested object, and the version flag which indicates both the consistency criteria and the versioning view. The different types of consistency criteria handled by Mistore are the read/write mutal exclusion criteria, the atomicity, and the eventual monotonic-read consistency one. If the read operation must be handled by the primary or the backup POPs then a single one is chosen (through a random choice). Figure 10 shows the code executed by HGs to handle a read operation, Figure 11 the code executed by a POP when it receives a read operation from a HG located in its region. Figure 12 shows the code executed by a POP to retrieve a requested object from its cache or from the HGs in its region and finally Figure 13 describes the code executed by a HG to retrieve an object from its local disk. In the following we present how the different consistency criteria are handled.

Read/write mutual exclusion: When a read operation is invoked with the read/write mutual exclusion flag, a read lock must be provided. When a POP receives this request from a HG, it maintains the read lock active until the requested object has been received by the HG that invoked the read operation.

A read lock contains the last version number of an object o and ensures that as long as the read lock is active, no concurrent write operation on o will be executed. Thus, the read operation can read the object from the cache of any POP if it contains the last version of o, otherwise the request must be forwarded to the primary or the backup of the object.

Atomic consistency: When a read operation is invoked with the atomic consistency flag, no read lock must be provided. The read operation is forwarded to the primary or the backup of the object.

Eventual monotonic-read: When a read operation is invoked with the eventual monotonic-read flag, no read lock must be provided but the version flag must contain the most recent version number known by the client. The read operation can read the object from the cache of any POP if it contains a version equal to or newer than the version flag, otherwise the request must be forwarded to the primary or the backup of the object.

Versioning view: An object may be read from the cache of any POP or HG if it contains a version of the object equal to the one requested with the version flag, otherwise the request must be forwarded to the primary or the backup of the object.

VII. EVALUATION

We have implemented Mistore on the Peersim simulator and performed all the simulations based on network parameters observed in an ISP infrastructure (see Table I). Fig. 14 illustrates different write requests schemes. We observe that latencies are the lowest when objects are written synchronously on both their primary and backup POPs, the best case being when the write requests originate from the primary or the backup region. This gives us the insight that the primary and backup POPs of an object should be the POPs of the regions from which originate most of the write requests on the object. We also observe that writing synchronously an object in a remote HG shows the worst performance due to the low users network bandwidth. Both observations validate our choice to replicate objects on their primary and backup POPs in a synchronous way, and to replicate them asynchronously in the HGs of their regions. However the performance of writing objects in HGs can be considerably improved when a data stripping method is applied on objects (i.e. the object is fragmented in several blocks that are stored in parallel in different HGs). Nevertheless we observe that gains obtained by using data stripping are more substantial on big size objects than on smaller ones. Fig. 15 illustrates performances of different read requests scenarios on remote nodes. First, latencies are the lowest when objects are read on a POP, the best case being the one when this POP is the one of the region from where the read request comes. Reading an object entirely on a remote HG is in general slow due to the low upload bandwidth of users. However, the latencies can be reduced when objects are read from several HGs in parallel. This motivates to use data stripping methods when storing objects such that they can be retrieved by aggregating the bandwidth of several HGs. Nevertheless, as for write requests, data stripping is more efficient on big size data objects. These results give us some insights for cache replacement policies. Actually, differences, regarding latencies, between storing small objects on their primary and backup or synchronously in remote HGs is not very large so we believe that big objects should have the priority to be kept in cache. Moreover, small objects may also be stored synchronously in remote HGs without a big impact on performance. We have also evaluated the impact of HGs and POPs failures on write and read operation latencies. Prior to describing simulation results, we briefly present how the system reacts to failures.

Crash of HGs: When a HG fails by crashing, the list of objects it stores is retrieved from the index maintained by the POP in its region. To keep a given degree of redundancy, these objects may need to be recovered from other replicas and replicated in another HG. The results of the pending requests of a lost HG are cached in the POP in its region. Thus, when the HG recovers, acknowledgments are sent to the HG for write requests, and read requests are satisfied by accessing the cache of the POP in the region.

Crash of POPs:

The failure of a POP means no Internet connection for the clients in its region. Thus we cannot not evaluate read and write latencies for users in a region where the POP is unavailable. Moreover, the crash of a POP means the crash of the primary or the backup of several objects. Without the primary of an object, its consistency can not be ensured so a new primary must be elected. As soon as the crash of a primary of an object is detected, its backup becomes the new primary and the datacenter becomes the new backup when an operation requiring a lock is activated. The backup can then trigger the lease of the locks it holds and takes over the service to the point where the old primary crashed. A crash of a POP is considered to be transient so a primary takes over the service when it recovers. Fig. 16 shows that during a write operation, if either the primary or the backup POP of the object becomes available, the operation can finish without incurring a large latency overhead. This is due to the fact that the replication occurs in the IP core network, and thus the low user bandwidth is not involved when an object has to be replicated to a datacenter due to the unavailability of one of its primary or backup POP. Fig. 17 shows the results of read operations during POPs and HGs unavailability. As argued before, we do not plot scenarios where the closest POP is unavailable -as it cuts users Internet connection -nor scenarios where the primary or the backup is unavailable -it does not incur any overhead as the object may be read on the available primary or backup POP. Similarly, we do not plot scenarios where HGs are unavailable since the POP knows which HGs are available and then will choose the available HGs to retrieve the requested objects. Now, when the POP (primary or backup) becomes unavailable after retrieving objects from the HGs and before forwarding them to their requester, we can observe that the overhead regarding latencies is very large when objects are not stripped over several HGs.

VIII. CONCLUSION AND FUTURE WORK

In this paper we have presented the main principles of Mistore, a distributed storage system to backup and share content dedicated to users who access Internet via a DSL technology. For future work, we plan to bring more storage intelligence close to the users (i.e. in POPs or/and HGs). We also seek to integrate efficient placement algorithms to minimize the amount of metadata to maintain in order to place/locate objects in Mistore. We are convinced that our architecture could also be suitable for user-generated content distribution and storage costs reduction via a multi-tiered data management. Evaluation of those points are other directions for future work.

Figure 1 .

 1 Figure 1. Simple overview of a DSL infrastructure of an ISP.

Figure 4 .

 4 Figure 4. Point of Presence : Function GRANTLOCK.

Figure 9 .

 9 Figure 9. Home Gateway: Function STORE.

Figure 14 .

 14 Figure 14.Write latencies. HG-DC: writing on the datacenter. HG-POP-P/B: writing an object on both its primary and backup POPs from a region different from these POPs region. HG-POP-P/B-HG: writing an object on HGs of both its primary and backup POPs from a region different from the primary and backup regions. HG-P/B: writing an object on both its primary and backup POPs from one these POPs region. HG-P/B-HG): writing on HGs of both the primary and backup POPs of an object from one of these POPs regions. No stripping: an object replica is written entirely in one HG. Stripping x20: an object replica is fragmented in 20 blocks stored in parallel in 20 different HGs.

Figure 15 .

 15 Figure 15.Read latencies. HG-DC: reading on the datacenter. HG1-POP1-POP2: reading on the POP of a region different from the one where the request is issued. HG1-POP1-POP2-HG2: reading on one HG connected to the POP of a region different from the one where the request is issued. HG1-POP1: reading on the POP of a region from where the request is issued. No stripping: the object is read entirely in 1 HG. Stripping x20: the 20 blocks stored in different HGs composing the object are read in parallel.

Figure 16 .Figure 17 .

 1617 Figure 16. Write latencies in presence of crashes. See Fig.14for meanings of HG-P/B and HG-POP-P/B. HG-[P/B]: writing from the primary POP region of an object and the primary or the backup POP is unavailable. Notation "[.]" refers to the failed element. HG-POP-[P/B]: writing of an object from a region different from its primary and backup regions while either the primary or the backup POP is unavailable.

  where → rf is a partial order relation called read from order and defined as follows. For any two operations op 1 ∈ H and op 2 ∈ H, we say that op 1 → rf op 2 if and only if 1) either ∃i such that op 1 ∈ h i and op 2 ∈ h i and op 1 → i op 2 , 2) or ∃x, v such that op 1 = w(x)v and op 2 = r(x)v, 3) or ∃op 3 ∈ H such that op 1 → rf op 3 and op 3 → rf op 2 . Given a global history H = (H, → H ), a linear extension of H and denoted H = (H, →) is a total order → on H that is compatible with → H , that is for any two operations op 1 ∈ H and op 2

	Definition 2 (Linear extension of a global history):

  1 → rt op 2 , if and only if op 1 .end ≤ op 2 .start. We say that two operations op 1 and op 2 are real-time concurrent, which will be denoted by op 1 rt op 2 , if and only if neither op 1 → rt op 2 , nor op 2 → rt op 1 .

	Definition 7 (Real time concurrency): Definition 8 (Atomicity): A global history H = (H, → H )
	is atomically consistent if it admits a linear extension H =
	(H, →) such that:

  with → H : for any operations op 1 and op 2 in H i , we have op 1 → i op 2 ⇔ op 1 → op 2 . Definition 12 (Read monotonic local history): Let H = (H, → H ) be a global history. H is consistent with the read monotonic consistency criteria if and only if for all processes p i ∈ Π, the local history (H i , → Hi ) is sequentially consistent.