Frédéric Alauzet
email: frederic.alauzet@inria.fr

A parallel matrix-free conservative solution interpolation on unstructured tetrahedral meshes

Keywords: Solution interpolation, matrix-free conservative interpolation, parallel interpolation, unstructured mesh, mesh adaptation, conservation laws Interpolation de solution, interpolation conservative sans matrice, interpolation parallèle, maillage non-structuré, adaptation de maillage, adaptation de maillage, loi de conservation

This document presents an interpolation operator on unstructured tetrahedral meshes that satisfies the properties of mass conservation, P 1 -exactness (order 2) and maximum principle. Interpolation operators are important for many applications in scientific computing. For instance, in the context of anisotropic mesh adaptation for time-dependent problems, the interpolation stage becomes crucial as the error due to solution transfer accumulates throughout the simulation. This error can eventually spoil the overall solution accuracy. When dealing with conservation laws in CFD, solution accuracy requires enforcement of mass preservation throughout the computation, in particular in long time scale computations. In the proposed approach, the conservation property is achieved by local mesh intersection and quadrature formulae. Derivatives reconstruction is used to obtain a second order method. Algorithmically, our goal is to design a method which is robust and efficient. The robustness is mandatory to obtain a reliable method on real-life applications and to apply the operator to highly anisotropic meshes. The efficiency is achieved by designing a matrix-free operator which is highly parallel. A multi-thread parallelization is given in this work.

Several numerical examples are presented to illustrate the efficiency of the proposed approach.

Introduction

Solution interpolation or solution transfer is an important stage for several applications in scientific computing. For instance, this stage is important for coupled problems, e.g. multiphysics simulations or fluid-structure interaction (FSI) problems, where specific meshes are considered for each sub-problem, see [START_REF] Geuzaine | A Galerkin projection method for mixed finite elements[END_REF][START_REF] Jiao | Common-refinement-based data transfer between non-matching meshes in multiphysics simulations[END_REF]. The preservation of the conservation property during the interpolation stage is also crucial for the accuracy in long time scale simulations, see [START_REF] Ringler | A potential enstrophy and energy conserving numerical scheme for solution of the shallow-water equations on a geodesic grid[END_REF][START_REF] Thuburn | Some conservation issues for the dynamical cores of NWP and climate models[END_REF]. It is an essential component of Arbitrary Lagrangian-Eulerian (ALE) methods as well. An accurate remapping must satisfy several properties such as conservation, high order accuracy, bound preserving, etc. This is also a key point in the context of time-accurate mesh adaptation. Indeed, it links the mesh generation and the numerical flow solver allowing the simulation to be restarted from the previous state. More precisely, after generating a new adapted mesh, called current mesh, the aim is to recover the solution field defined on the previous mesh, called background mesh, on this new mesh to pursue the computation. This recurrent stage in adaptive simulations is crucial for time-accurate unsteady problems as errors introduced by the interpolation procedure accumulate throughout the computations. The negative impact of such errors on solution accuracy was pointed out in [START_REF] Alauzet | 3D transient fixed point mesh adaptation for time-dependent problems: Application to CFD simulations[END_REF] where standard linear interpolation is applied. In [START_REF] Alauzet | P1-conservative solution interpolation on unstructured triangular meshes[END_REF], it has been demonstrated -for the 2D case -the importance of a conservative interpolation operator on the accuracy of anisotropic timeaccurate mesh adaptation for conservation laws.

A conservative interpolation based on a Galerkin projection has been proposed in [START_REF] Geuzaine | A Galerkin projection method for mixed finite elements[END_REF][START_REF] Farrell | Conservative interpolation between unstructured meshes via supermesh construction[END_REF][START_REF] Menon | Conservative interpolation on unstructured polyhedral meshes: An extension of the supermesh approach to cell-centered finite-volume variables[END_REF]. The Galerkin projection requires the assembling and the resolution of a linear system. This method can as well suffer from oscillations. To bound the projected quantities specific treatments have to be done [START_REF] Farrell | Conservative interpolation between unstructured meshes via supermesh construction[END_REF]. A global supermesh construction is used in [START_REF] Farrell | Conservative interpolation between unstructured meshes via supermesh construction[END_REF][START_REF] Menon | Conservative interpolation on unstructured polyhedral meshes: An extension of the supermesh approach to cell-centered finite-volume variables[END_REF] but it limits considerably the efficiency of the method and the maximal size of the problem in 3D due F. Alauzet to memory usage. A local approach has been proposed in [START_REF] Alauzet | P1-conservative solution interpolation on unstructured triangular meshes[END_REF] which improves the algorithm efficiency. It is based on a local mesh intersection procedure and local quadrature. This local approach has been followed in [START_REF] Farrell | Conservative interpolation between volume meshes by local galerkin projection[END_REF] coupled with local Galerkin projection. In [START_REF] Alauzet | P1-conservative solution interpolation on unstructured triangular meshes[END_REF], the maximum principle is strictly enforced using L 2 -optimal gradient correction, the method is thus free from oscillations.

In this paper, we consider the 3D solution interpolation for anisotropic adapted tetrahedral meshes where the background and the current meshes are distinct, in the sense that the number of entities and the connectivities can be completely different. Flows are modeled by the conservative compressible Euler equations and resolved by a second order finite volume scheme. Therefore, to obtain a consistent mesh adaptation loop, the proposed interpolation scheme must satisfy the following properties:

• mass conservation • P 1 exactness implying an order 2 for the method • maximum principle.

Moreover, this method has to be algorithmically very robust as we deal with highly stretched elements and it has to be very efficient to be applied to real-life applications. The word efficient signifies that it requires low memory storage and that the additional required CPU time over that for standard linear interpolation is acceptable. In consequence, we propose a matrix-free approach based on local mesh intersections and appropriate local reconstructions.

The mass conservation property of the interpolation operator is achieved by local mesh intersection, i.e., intersections are performed at the element level. The use of mesh intersection for conservative interpolation seems natural for unconnected meshes and has already been alluded in [START_REF] George | Delaunay triangulation and meshing[END_REF] or applied in [START_REF] Grandy | Conservative remapping and regions overlays by intersecting polyhedra[END_REF] for order 1 reconstruction. The locality is inherent for efficiency and robustness. Once again for efficiency purposes, the proposed intersection algorithm is especially designed for simplicial meshes. The idea is to compute the intersection between two simplexes, mesh this intersection, and use a quadrature formulae to exactly compute the transferred mass. Moreover, the designed algorithm is highly scalable in parallel due to its locality.

The high-order accuracy is obtained by a solution gradient reconstruction from the discrete data and the use of Taylor formulae. This high-order interpolation can lead to loss of monotonicity. The maximum principle is then enforced by correcting the interpolated solution, thus the interpolated solution is free from any oscillations. Notice that much care has been taken while designing the localization algorithm as it is also critical for efficiency.

The proposed P 1 -conservative interpolation operator is suitable for solutions defined at elements or vertices. This paper is the extension of [START_REF] Alauzet | P1-conservative solution interpolation on unstructured triangular meshes[END_REF] to the three-dimensional case. As compared to [START_REF] Alauzet | P1-conservative solution interpolation on unstructured triangular meshes[END_REF], several novelties are presented in this work. They are required to ensure the robustness 1 , efficiency and accuracy of the method when dealing with 3D highly anisotropic meshes on complex geometries. These novelties concern:

• A specific tetra-tetra intersection procedure where floating point arithmetic is treated with a particular care, as it must be extremely robust and able to deal with highly anisotropic tetrahedra (anisotropic ratio up to 10 5 )

• A new method to mesh the tetra-tetra intersection ensuring extra numerical accuracy

• Solutions to deal with the case of non-matching boundaries between meshes (this is always the case for complex geometries), this is crucial to not introduce artifacts in the solution during the interpolation stage and preserve accuracy

• The parallelization of the interpolation stage for efficiency

• A convergence analysis pointing out the loss of convergence using P 1 interpolation which is recovered using the P 1 -conservative interpolation, thus advocating the use of the latter one.

As regards the differences and similarities between this work and [START_REF] Farrell | Conservative interpolation between volume meshes by local galerkin projection[END_REF]. In [START_REF] Farrell | Conservative interpolation between volume meshes by local galerkin projection[END_REF], they use local approach and an advancing front method for the intersection identification similar to the ones proposed in [START_REF] Alauzet | P1-conservative solution interpolation on unstructured triangular meshes[END_REF]. The first difference lies in the computation and the meshing of the intersection between two elements. In [START_REF] Farrell | Conservative interpolation between volume meshes by local galerkin projection[END_REF], they follow the Eberly clipping algorithm [START_REF] Eberly | 3D game engine design: a practical approach to real-time computer graphics[END_REF] (?which is not described in 3D while being the core of the algorithm?) while we propose specific algorithm to achieve greater numerical accuracy and to have a relative mass variation of the order of the round-off (≈ 10 -14 ). This algorithm is still accurate when dealing with highly anisotropic elements (ratio up to 10 5 ). The second difference is the projection operator. A local Galerkin projection is used in [START_REF] Farrell | Conservative interpolation between volume meshes by local galerkin projection[END_REF] while this work considers quadrature formulae and high-order solution reconstruction. As last, we observe that the proposed algorithm achieve better efficiency in serial, i.e., one order of magnitude faster, and it achieves good scaling in parallel.

The paper is outlined as follows. Section 2 introduces the main definitions and Section 3 recalls the standard linear interpolation operator. Then in Section 4, the proposed P 1conservative interpolation operator is described. First, the mesh intersection algorithm is presented and at a second stage, P 1 -conservative reconstruction is discussed. Finally, the accuracy of the proposed approach is emphasized on analytical examples in Section 5, and parallel efficiency is discussed in Section 6. The approach is successfully applied to adaptive numerical simulations in Section 7. Some concluding remarks close the paper.

Definitions and notations

This section provides notations, definitions and conventions used in this paper. Let us consider a bounded domain Ω ⊂ R 3 and denote by ∂Ω its boundary. Domain Ω is discretized by a tetrahedral mesh H = K i . A tetrahedron K i is defined by the list of its vertices which are locally numbered in a convenient way. This list, enriched with some conventions, provides the complete definition of the related element, including the definition of its faces, edges and neighbors, together with an orientation. In particular, the oriented local numbering of tetrahedron's vertices enables us to compute its volume while giving a sense to its sign, and to evaluate directional normals for each face.

The local numbering of vertices, edges, faces and neighboring tetrahedra is pre-defined in such a way that some properties are implicitly induced. In the case of a tetrahedron with vertices [p 0 , p 1 , p 2 , p 3 ], the vertices are ordered such that a positive signed volume is obtained while using the scalar triple product Formula (see Figure 1):

V K = |K| = 1 6 e 0 • (e 1 × e 2 ) . (1)
Faces topology is defined as 1, left). This numeration is such that the index of the face is the index of the viewing vertex, i.e., the opposite vertex. And, the four faces have the same inward orientation, in other words, the face normal n Fi = e i,1 × e i,2 is inward, where e i,1 is the edge composed by the first and the second face vertices and e i,2 is the edge composed by the first and the third face vertices. To improve the face normal computation accuracy, we choose the two edges forming an angle the closest as possible to 90 • , i.e., the minimal dot product. The surface area of the face is simply given by:

F 0 = [p 3 , p 2 , p 1 ], F 1 = [p 2 , p 3 , p 0 ], F 2 = [p 1 , p 0 , p 3 ], F 3 = [p 0 , p 1 , p 2 ] (Figure
A Fi = |F i | = 1 2 n Fi = 1 2 e i,1 × e i,2 .
Regarding the neighboring tetrahedra, they follow the face convention. We denote by K i the neighbor viewing vertex p i through face F i . For the edges topology, the following choice has been made:

e 0 = [p 0 , p 1 ], e 1 = [p 0 , p 2 ], e 2 = [p 0 , p 3 ], e 3 = [p 1 , p 2 ], e 4 = [p 1 , p 3 ], e 5 = [p 2 , p 3 ].
With the above notations, we now give some definitions used in the sequel. Let p be a point, we denote by K i the virtual tetrahedron where vertex p i is substituted by p. The F. Alauzet signed volumes V K i , for i = 0 . . . 3, are called the barycentrics of p. The four associated barycentric coordinates are given by:

β i = V K i V K for i = 0 . . . 3 .
The sign of the four barycentrics defines explicitly 15 regions of space where point p can be located with respect to element K. The possible combinations are given in Figure 1 (right). Now, we recall some definitions relative to the interpolation schemes. Let u be a solution defined on a mesh H 1 of a domain Ω. The mass of the solution over the mesh is simply m = H 1 u. We deduce the notion of mass on an element K given by m K = K u. An interpolation scheme is said to be conservative if it preserves the mass when transferring the solution field u from a mesh H 1 to another H 2 . Formally speaking, if we denote by Πu the interpolated field on H 2 , then such scheme verifies

H 1 u = H 2 Πu .
A scheme is said to be P k -exact if it is exact for polynomial solutions of degree lower than or equal to k. A P k -conservative interpolation scheme is a scheme satisfying both properties. Finally, an interpolation scheme verifies locally the maximum principle on a subset

Q ⊂ H 1 if for p ∈ H 2 , min q∈Q u(q) ≤ Πu(p) ≤ max q∈Q u(q) .

It verifies globally that property if

Q = H 1 .

Linear interpolation operator

Solution interpolation is a twofold process. First, vertices of the new mesh are located in the background mesh. Second, an interpolation scheme is applied for each vertex.

Localization algorithm

The localization problem or research of a point location consists in identifying the element of a simplicial mesh containing a given point. Here, we consider the simplified problem where the background and the new meshes are discretizations of the same domain Ω. This problem has to be dealt with great care in the case of simplicial meshes to handle difficult configurations. Indeed, background and current meshes can be non-convex and can contain holes. It is also possible that the overlapping of the current mesh does not coincide with the background mesh since their boundary discretization can differ. Consequently, some vertices of the current mesh can be outside of the background mesh and conversely. Moreover, efficient localization algorithms have to be implemented to avoid a quadratic complexity obtained with a naive algorithm.

The localization can be solved efficiently by traversing the background mesh using its topology, i.e., the neighboring elements of each element, thanks to a barycentric coordinatesbased [START_REF] Frey | Mesh generation. Application to finite elements, 2nd Edition[END_REF][START_REF] Löhner | Robust, vectorized search algorithms for interpolation on unstructured grids[END_REF]. The considered algorithm is presented in details for the 2D case in [START_REF] Alauzet | P1-conservative solution interpolation on unstructured triangular meshes[END_REF] with a discussion on how to death with pathological configurations Its extension to 3D is straightforward.

Classical linear interpolation

Once the localization has been performed, we can apply an interpolation algorithm. The easiest interpolation scheme is the classical P 1 interpolation:

Π 1 u(p) = 3 i=0 β i (p) u(p i ) .
where p is a vertex of the new mesh that has been localized in tetrahedron K = [p 0 , p 1 , p 2 , p 3 ] of the background mesh. β i are the barycentric coordinates of p w.r.t. K. This scheme is P 1 -exact and it is thus order 2. This scheme is monotone and satisfies the maximum principle. However, this scheme does not conserve the mass.

Matrix-free P 1 -exact conservative interpolation

In this section, a P 1 -exact conservative interpolation scheme is presented. The provided solution is considered to be piecewise (continuous or discontinuous) linear by element. The idea of the conservative interpolation is to compute the mass of each element of the new mesh H new knowing the mass of each element of the background mesh H back . For efficiency purposes, a local mesh intersection algorithm is utilized coupled with a matrix-free interpolation. Then, in the case of vertex-centered solution, the solution is transferred accurately and conservatively from elements to vertices using the mass of the elements of its ball. This process is summarized in Algorithm 1. The mesh intersection procedure corresponding to steps 3a and 3b is exposed in Section 4.1 and the conservative reconstruction, steps 3c, 4 and 5, is described in Section 4.

Mesh intersection algorithm

The mesh intersection algorithm consists in intersecting each tetrahedron of the current mesh with all the background mesh tetrahedra that it overlaps and in meshing the intersec-F. Alauzet tion region. The 3D mesh intersection is a lot harder than in 2D, thus the method of [START_REF] Alauzet | P1-conservative solution interpolation on unstructured triangular meshes[END_REF] cannot be applied. Specific attention to numerical accuracy is required to deal accurately with degenerated cases (occurring frequently for highly anisotropic meshes). To this end, the tetrahedron-tetrahedron intersection topology consistency is checked to detect numerical accuracy issues. If an inconsistency is detected a perturbation method is applied. Moreover, a dedicated meshing method of the intersection is proposed where extra accuracy is obtained for the orientation predicate.

In the following, we first describe our generic intersection algorithm between any pair of tetrahedra and how we discretize the intersection polyhedron (Section 4.1.1). Secondly, the algorithm to locate all background tetrahedra that are overlapped by the current element is presented (Section 4.1.2).

Tetrahedron-tetrahedron intersection

The tetrahedron-tetrahedron intersection procedure computes the intersection of two tetrahedra and meshes the intersection region if it is not empty. Notice that if the intersection exists, the intersection region of two tetrahedra is always a convex polygon given by the convex hull of the intersection points. This is a three steps procedure:

1. check all the degenerated intersection cases 2. evaluate the 48 possible edge-face intersections 3. steps 1. and 2. result in a cloud of points which is triangulated resulting in a mesh of the intersection polyhedron.

In the following, these three steps are described. But, first we give a few words about the robustness of the proposed algorithm as in 3D numerical accuracy and floating point arithmetic is of utmost importance.

Accuracy and robustness of the algorithm Accuracy and robustness of algorithms when implemented with floating point arithmetic is a major topic of study in computational geometry [START_REF] Hoffmann | The problems of accuracy and robustness in geometric computation[END_REF][START_REF] Stewart | Local robustness and its application to polyhedral intersection[END_REF]. However, this work follows the choice of the authors of [START_REF] Frey | Mesh generation. Application to finite elements, 2nd Edition[END_REF] where extended accuracy, interval analysis, arithmetic filters and sign or orientation predicates are discussed.

In our context, the most important is consistency meaning that the algorithm always deliver the same answer on the same configuration and the topology of the final result is correct. The consistency is crucial to avoid cycling and to have a unique answer in the localization process. It is also essential to obtain the same intersection point for a given pair edge-edge or edge-face as such intersection generally occurs many times, i.e., for many tetrahedron-tetrahedron intersections, in the presented algorithm. Checking the topology of the intersection is an efficient way to detect and to correct inconsistency due to floating point arithmetic. If the algorithm is consistent (and the initial data base correct) then the orientation predicate is reliable. The design of a consistent algorithm is based on handling properly geometric degeneracy and verifying the correctness of the topology of the intersection.

First, an adequate local ε is chosen to carefully handle degenerated cases without loss of accuracy. It needs to be consistent for all intersections with the current tetrahedron. Therefore, it is only based on the current tetrahedron size.

Second, the topology of the tetrahedron-tetrahedron intersection is stored and checked throughout and at the end of the process. Topology is powerful because it is exact as it only involves boolean operation. Figure 2 presents three degenerated intersections where floating point arithmetic is involved. The two left intersections are not valid topologically while the intersection on the right is valid. For all edges, we store all the entities (vertex, edge, face) they intersect. For instance, if an edge of K new intersects an edge of K back , then we store that these edges intersect each other and that they also intersect the faces they share. Moreover, for each intersection point, the entities (vertex, edge, face) of the current and of the background tetrahedra on which it lies on are stored. This set of data provides a complete view of the intersection polyhedron. To verify the topology validity, the number of intersection points on each entity is also checked. In particular, an edge and a face cannot

Valid configuration

Edge is intersecting two edges Figure 2: Three degenerated intersections are presented. The two intersections on the left are inconsistent topologically while the one on the right is valid. Indeed, for the left one, red edge e has an intersection with the top face but it has no other intersection with tetrahedron K. In such case, the red edge must intersect K twice or have a degenerated intersection with an edge of K. For the middle one, red edge e has an intersection with one edge of the side face and an intersection inside the side face which is not valid. Such an edge must either be coplanar with the side face, resulting in the configuration presented on the right, or have a non-degenerated intersection with the top face.

have more than 2 and 6 intersection points, respectively. And, the intersection polyhedron contains at most 12 intersection points [START_REF] Seidel | Convex hull computations[END_REF]. Cases with 12 intersection points are obtained, for instance, for the stella octangula polyhedron or by taking two opposite tetrahedra from Cundy and Rollett's construction of the tetrahedron 5-compound, [START_REF] Cundy | Mathematical Models[END_REF].

During the tetrahedron-tetrahedron intersection process, the topology consistency is checked to correct pathological cases due to floating point arithmetic and to update missed degenerated cases. After completing the tetrahedron-tetrahedron intersection, the topology consistency is analyzed. If an inconsistency occurs, then (i) the intersection is reset, (ii) a tiny perturbation of the background vertices position is performed and (iii) the intersection is resumed. The perturbation helps to get rid of unresolved degenerated cases that create the inconsistency. As the perturbation is of the order of ε, it doesn't affect the accuracy.

Remark:

The proposed choice to cope with floating point arithmetic is very efficient as it is able to accurately handle the intersection of highly anisotropic tetrahedra with aspect ratio of 1:100 000. and to recover the tetrahedron volume at the order of the round-off (≈ 10 -14 ).

Dealing with degenerated cases We first introduce the definition of the signed distance, also called power, of point p with respect to face F i = [p i,0 , p i,1 , p i,2 ]:

P(p, F i ) = p i,k p • n Fi n Fi
where k is either 0, 1 or 2 .

(

) 2
Notice that the barycentrics and the powers are linked by the relation:

V K i = 1 6 n Fi P(p, F i ) .
The distance of point p with respect to edge e i = p i,0 p i,1 is:

P(p, e i ) = p i,0 p i,1 × p i,0 p p i,0 p i,1 = p i,0 p × p i,1 p p i,0 p i,1 .
For the given pair of tetrahedra, the 32 vertices powers with respect to the faces, P(p j , F i ), are computed according to Equation [START_REF] Jiao | Common-refinement-based data transfer between non-matching meshes in multiphysics simulations[END_REF]. Then, all possible vertex degenerated cases are checked according to P(p j , F i ) values:

• is a vertex inside a face ?

• is a vertex on a edge ?

• are two vertices coinciding ?

If such a case happen, the vertex is snapped on the appropriate entity and added to the intersection points list. The intersection topology table is updated accordingly.

F. Alauzet

Afterward, edge-edge intersections -which are degenerated cases -are tested. Let p 0 p 1 and q 0 q 1 the two considered lines. A necessary condition is that the two lines are coplanar. If it is the case, the intersection point x is evaluated following Hill's approach [START_REF] Hill | The pleasures of 'perp dot' products[END_REF]:

x = p 0 + (p 0 q 0 × q 0 q 1 ) • (p 0 p 1 × q 0 q 1 ) p 0 p 1 × q 0 q 1 2 p 0 p 1 = p 0 + s p 0 p 1 , or x = q 0 + (p 0 p 1 × q 0 p 0 ) • (p 0 p 1 × q 0 q 1 ) p 0 p 1 × q 0 q 1 2 q 0 q 1 = q 0 + t q 0 q 1 .
Now, to check if the intersection point is an intersection between the two segments, we have to check that:

0 ≤ s ≤ 1 and 0 ≤ t ≤ 1 .
If there is intersection, point x is added to the intersection points list and the intersection topology table is updated.

Edge-face intersection Dealing with degenerated cases first simplifies the following edgeface intersection procedure. Indeed, the computation of all the above geometric degeneracy treats in particular all the possible coplanar edge-face intersections as depicted in Figure 3. Now, only non-coplanar edge-face intersection remains where the intersection point lies inside the face. The algorithm checks the 48 possible edge-face intersections. Let us denote the considered edge by e = [p 0 , p 1 ] and face by F = [q 0 , q 1 , q 2 ]. As it is a non-coplanar case, the edge's vertex powers with respect to the face are not zero: P(p 0 , F ) = 0 and P(p 1 , F ) = 0. A necessary condition for the edge-face intersection is:

P(p 0 , F ) P(p 1 , F ) < 0 .
If this condition is satisfied, the point of intersection between the edge and the plane defined by the face is computed:

x = p 0 + P(p 0 , F ) P(p 0 , F ) -P(p 1 , F ) p 0 p 1 .
Finally to verify if the edge-face intersection effectively occurs, i.e., point x lies inside face F , the signs of barycentrics of point x with respect to face F are analyzed. If there is intersection, point x is added to the intersection points list. Meshing the intersection polyhedron The intersection procedure results in a cloud of points formed by the intersection points list. This list is analyzed to mesh the intersection polyhedron. If the intersection process returns strictly less than 4 intersection points, then we face a degenerated tetrahedron-tetrahedron intersection where there is no geometric intersection or the geometric intersection is either a vertex, a line or a face. In this context, the algorithm has to make the distinction between the case where one tetrahedron is included in the other one, meaning that the intersection is this tetrahedron, and the case where the intersection is empty. Tetrahedron K p = [p 0 , p 1 , p 2 , p 3 ] is included inside tetrahedron K q = [q 0 , q 1 , q 2 , q 3 ] if and only if ∀i = 0, .., 3 and ∀j = 0, .., 3, P(p i , F q j ) ≥ 0, where F q j are the faces of K q . As all the powers are evaluated at the beginning of the intersection process, the included tetrahedron case can be immediately checked. Therefore, if the number of intersection points is strictly less than 4, it implies that no intersection occurs.

Otherwise, if more than 4 intersection points are returned, the convex hull of this cloud of points forms a convex polyhedron representing the region of intersection of the tetrahedra pair.

At this point, many strategies can be employed to mesh the intersection volume. The more intuitive idea is to create the Delaunay triangulation of the cloud of points -similar to the one used in [START_REF] Alauzet | P1-conservative solution interpolation on unstructured triangular meshes[END_REF] -as it returns the convex hull of the set of points. However, a basic algorithm lacks of robustness due to the presence of many coplanar points inside faces and the frequent nearly-degenerated configurations when two highly anisotropic tetrahedra are involved. In such context, orientation predicate may fail leading to a wrong configuration. To solve this issue, one can consider the use of a very advanced Delaunay triangulation algorithm such as the one proposed by [START_REF] George | Delaunay's mesh of a convex polygon in dimension d. Application to arbitrary polyedra[END_REF]. But, this seems incongruous to call such an evolved algorithm at each intersection to triangulate a dozen of points and it may face the same issues.

Therefore, we propose a clever strategy using the information gathered during the intersection process. Indeed, as an intersection topology table has been created, we know which intersection point belong to which tetrahedron's face. A first step consists in meshing all the faces of the convex polyhedron. Each of these faces corresponds to a face of one of the two tetrahedra. Thus, for each created triangle on these faces the orientation predicate is based on the orientation of the normal which must be identical to the orientation of the corresponding tetrahedron's face normal. This gives extra accuracy because this predicate is based on a cross product whereas a volume predicate is based on a triple product. Once all the tetrahedra's faces have been meshed, an oriented surface mesh of the intersecting polyhedron is obtained. In order to get a simplicial volume mesh of the polyhedron, its gravity center is added and tetrahedra are created by joining the barycenter to each triangle. In other words, the polyhedron surface mesh is star-shaped with respect to the barycenter.

Two examples of tetrahedron-tetrahedron are shown in Figure 4.

Overlapped tetrahedra detection

The conservative interpolation method consists in computing for each tetrahedron K new of the current mesh H new its intersection with all tetrahedra K back j of the background mesh H back that it overlaps. We present how this list of background elements is determined.

First of all, all vertices of the new mesh H new are localized in the background mesh H back using the algorithm presented in Section 3.1. Then, for each tetrahedron K new of H new , the initial list of background tetrahedra that are overlapped is given by the elements containing the vertices of K new . For degenerated cases where a new vertex lies on a background face or on a background edge or on a background vertex, we add to the initial list the background tetrahedra that are sharing the face or are in the edge's shell or are in the vertex's ball, respectively.

Then, the intersections between K new and the tetrahedra of the initial list are computed. New tetrahedra are added to the list during the intersection procedure as follows:

• if face F j of K back is intersected by K new then neighbor K back j of K back sharing face F j is added to the list • if K new is
intersected by an edge of K back then background tetrahedra of the edge's shell are added to the list

• if a vertex of K back lies inside K new then background tetrahedra of the vertex's ball are added to the list.

With this simple procedure, all overlapped elements are automatically detected while computing intersections. Overlapped elements are detected without any additional cost as powers and intersections have already been computed in the tetrahedron-tetrahedron intersection process.

P 1 -conservative reconstruction

In this section, we describe the P 1 -conservative solution reconstruction process. It mainly follows [START_REF] Alauzet | P1-conservative solution interpolation on unstructured triangular meshes[END_REF], but choices made to cope with the case of non-matching boundaries between meshes (this is always the case for complex geometries) are detailed. We consider a bounded domain Ω of R 3 and two tetrahedral meshes of this domain

H back = K back i and H new = K new i .
For sake of simplicity, we first make the assumption that the discrete boundaries of both meshes are the same, i.e., both meshes are discretization of the same polyhedral domain

Ω h : |H back | = |H new | where |H| = Ω h dx.
The case of non-matching discrete boundaries is addressed in Section 4.2.4. For each mesh, a dual partition of the domain is defined by associating with each vertex a control volume or cell (which is defined by some rules):

H back = C back i and H new = C new i . A P 1 discrete solution field u is given on the background mesh H back .

Now, we have to define a projection operator Π c

1 from H back to H new with the following properties:

• Π c
1 is conservative:

H back u = H new Π c 1 u • Π c 1 is P 1 -exact: if u is affine then Π c 1 u = u.
The projection operator is presented for solutions defined at elements and solutions defined at vertices.

Solution defined at elements

In this case, the solution is piecewise linear by elements and can be discontinuous. We have for each background tetrahedron K back :

• mass m K back = K back u = |K back | u(g), where g is the barycenter of K back • constant gradient ∇u K back .

Inria

For each tetrahedron K new of the current mesh, we compute the intersection with all tetrahedra of the background mesh {K back j } j it overlaps as described in the previous section. Each pair of tetrahedra K new and K back j provides a simplicial mesh of their intersection region denoted T j = K new ∩ K back j . The integrals Tj u and Tj ∇u are computed exactly using Gauss quadrature formulae. Consequently, we obtain for each tetrahedron of the current mesh a mass and a gradient given by:

m K new = K new Π c 1 u = j Tj u and (∇ Π c 1 u) |K new = j Tj ∇u |K new | .
This reconstruction is conservative and P 1 -exact. It gives a P 1 by element discontinuous solution. A specific treatment of the reconstruction is carried out to verify the maximum principle.

Verifying the maximum principle

Let K be a tetrahedron of the new mesh. In the following, for the sake of clarity, we denote by u K the P 1 -conservative interpolated solution Π c 1 u on K. The value at the barycenter and the gradient of the interpolated solution on K are given by:

u K (g K ) = 1 |K| K Π c 1 u and ∇u K = (∇Π c 1 u) |K .
Consequently, for each vertex p i of the new mesh, a value of the solution is obtained using Taylor expansion for each element K of its ball:

u K (p i ) = u K (g K ) + ∇u K • g K p i . (3)
A correction is applied to the linear representation of the solution on each element in order to verify the maximum principle. The interpolated solution is thus free from any oscillations. To this end, let K be the set of elements of the background mesh that K overlaps and let Q be the set of vertices of K:

K = {K back j | K ∩ K back j = ∅} and Q = {q j | q j ∈ K back such that K back ∈ K} .
Then, for each vertex p i of each element K of the new mesh, the nodal value u K (p i ) verify the maximum principle if:

u min = min q∈Q u(q) ≤ u K (p i ) ≤ max q∈Q u(q) = u max .
Notice that u K (g K ) always satisfies the maximum principle. If a vertex does not verify the maximum principle on an element K then the gradient value of this element is corrected. The proposed approach results from a minimization problem2 . We first reorder the indices such that u

K (p 0 ) ≤ u K (p 1 ) ≤ u K (p 2 ) ≤ u K (p 3
). Then, we set:

u M K (p3) = min uK (p3) , umax u M K (p2) = min uK (p2) + 1 3 (uK (p3) -u M K (p3)) , umax u M K (p1) = min uK (p1) + 1 2 3 i=2 (uK (pi) -u M K (pi)) , umax u M K (p0) = uK (p0) + 3 i=1 (uK (pi) -u M K (pi))
and

ũK (p0) = max u M K (p0) , umin ũK (p1) = max u M K (p1) + 1 3 (u M K (p0) -ũK (p0)) , umin ũK (p2) = max u M K (p2) + 1 2 1 i=0 (u M K (pi) -ũK (pi)) , umin ũK (p3) = u M K (p3) + 2 i=0 (u M K (pi) -ũK (pi)) .
These new nodal values ũK (p i ) define the corrected linear representation of the solution on K. For any points x included in K, its solution value is then given by: ũK (x) =

F. Alauzet 3 i=0 β i (x)ũ K (p i )
, where β i (x) are the barycentric coordinates of x with respect to K. The final interpolated solution verifies all required properties: Proposition 1: The reconstruction ũK satisfies the maximum principle, is linear preserving and is conservative. Moreover, we have:

u K (p 0 ) ≤ u K (p 1 ) ≤ u K (p 2 ) ≤ u K (p 3 ) =⇒ ũK (p 0 ) ≤ ũK (p 1 ) ≤ ũK (p 2 ) ≤ ũK (p 3 ) and if we have u min ≤ u K (p i ) ≤ u max for i = 0..3 then ũK (p i ) = u K (p i ) for i = 0..3.
Notice that this reconstruction comes from a minimization problem. Indeed, we have

Proposition 2: Suppose that u K (p 0 ) ≤ u K (p 1 ) ≤ u K (p 2 ) ≤ u K (p 3
) and that u max < u K (p 3 ). Then, we have

3 i=0 |u K (p i ) -u M K (p i )| 2 ≤ 3 i=0 |u K (p i ) -v K (p i )| 2
for all the linear reconstructions v

K satisfying v K (p 3 ) = u max , v K (p i ) ≤ u max for i = 0, .., 2 and K v K = K u K .
The proofs of these propositions are identical to the 2D proofs given in [START_REF] Alauzet | P1-conservative solution interpolation on unstructured triangular meshes[END_REF].

Solution defined at vertices

When the solution is given at vertices of the background mesh, i.e., nodal values are provided, it implicitly defines a piecewise linear continuous representation of the solution at the elements. Therefore, the P 1 -conservative interpolation defined in the previous section can be applied. However, a piecewise linear solution by elements, which is generally discontinuous, is obtained on the new mesh. Therefore, one more stage is required to retrieve a solution at vertices of the new mesh which consists in transferring this solution from elements to vertices while preserving the properties of the interpolation scheme. The solution is simply re-distributed to each vertex p of the new mesh by averaging:

ũ(p) = K new i p |K new i | ũK new i (p) K new i p |K i | ,
where ũ is the interpolated solution on the new mesh. Notice that after re-distribution to vertices the interpolated solution still satisfies the maximum principle, is linear preserving and is conservative.

Remark: As the mass of the solution is linked to the topology of the mesh, the P 1 -conservative interpolation operator Π c 1 depends on the mesh topology on which it is applied. In consequence, it cannot be applied to interpolate solution at any points of a given domain.

Non-matching discrete boundaries

Let Ω be a bounded domain of R 3 and H back and H new two meshes of Ω. We consider the case where the discrete boundaries of H back and H new do not match. In other words, H back and H new are meshes of two different polyhedral domains Ω back For the conservative interpolation, the non-matching discrete boundaries are handled differently depending on the solution behavior and the geometric configuration. This is crucial to not introduce artifacts in the solution during the interpolation stage and preserve accuracy. Let's take the simple example of a sphere in a constant flow field delimited by a box domain. If the geometric approximation of the sphere is increased (the surface mesh size on the sphere is divided by a factor 2), then the volume of the sphere increases (because it is convex) leading to a decrease of the domain volume. In that case, we obviously want to preserve the constant flow field when the solution is transferred from one mesh to the other, hence the conservation property must be violated. Satisfying the P 1 -exactness is more important than being conservative. Three specific solution cases may occur: the solution is constant, linear or something else. And, for the geometry configuration, we may have either a background tetrahedron or a new tetrahedron which is not entirely overlapped.

When a tetrahedron of the background mesh is not completely overlapped, we choose to do nothing specific thus some of the global mass is not recovered due to the volume change. The algorithm is no longer conservative but it still preserves the P 1 -exactness and the maximum principle properties.

When a tetrahedron of the new mesh is not completely overlapped, we have :

|K new | = j |K new ∩ K back j | = j |T j | = |T ∩ | .
This case is treated based on the solution behavior as follows.

If the solution is constant, the operator must recover a constant solution to satisfy the P 1 -exactness and the maximum principle properties. Indeed, if a constant uniform field is given as the initial set, we have to return an constant uniform solution set to not introduce any artifacts in the solution. The conservation principle must be violated. Formally speaking, we have:

m ∩ = c |T ∩ | and (∇ Π c 1 u
) |K new = 0 , and, the value at barycenter is obtained by dividing the mass by the intersection volume and not the new tetrahedron volume:

u K (g K ) = m ∩ |T ∩ | = c ⇒ m K new = c |K new | .
If the solution is linear, the operator must recover a linear solution to satisfy the P 1exactness, once more to not introduce any artifacts in the solution. The above method cannot be applied because the recovered mass is not correct to compute the barycenter value. However, the exact gradient value can be recovered:

(∇ Π c 1 u) |K new = j Tj c |T ∩ | = c .
Therefore, to retrieve the linear solution, we pick one of the new tetrahedron vertex that is inside the domain and a linear extrapolation is applied to obtain the other vertices values.

In such case, the maximum principle and the conservation properties are no more verified. Finally, if the solution is not in the two previous cases, we choose to preserve the maximum principle and to not be conservative. Values and gradients at barycenter are computed as:

u K (g K ) = j Tj u |T ∩ | and (∇ Π c 1 u) |K new = j Tj ∇u |T ∩ | .
and values at vertices are computed with Relation (3).

Accuracy and convergence study on analytical functions

In this section, the behavior of the P 1 -conservative interpolation is analyzed on four analytical functions defined on a cubic domain [-0.5, 0.5] 3 . These functions are representative of several physical phenomena encountered in computational fluid dynamics (CFD). The function's solutions are considered at vertices. The P 1 -conservative interpolation is compared to the linear interpolation, in particular, the mass conservation and the convergence order of the schemes are studied.

To perform this analysis, two meshes H 1 1 and H 2 1 , composed respectively of 906 and 918 vertices, are considered. These meshes are completely different and unconnected. In order F. Alauzet to study the convergence order of each interpolation method, each of these meshes spans a series of embedded meshes denoted (H 1 i ) i=1...5 and (H 2 i ) i=1...5 . Mesh H j i+1 is deduced from H j i by splitting each tetrahedra into eight tetrahedra in a Lagrangian fashion, i.e., in an isoparametric way. These series of meshes are summarized in Table 1.

For each case, the analytical function is applied on H 1 i providing a solution field u 1 i . This solution field is transferred from H 1 i to H 2 i , we get Πu 2 i . This solution transfer is called transfer

H 1 i → H 2 i .
The error is computed by comparing the interpolated solution Πu 2 i with the analytical function applied on H 2 i , i.e., u 2 i , in L 1 -norm:

ε i = H 2 i |u 2 i -Πu 2 i | .
The series of errors enable a convergence study. We also analyze the error when the solution field is re-interpolated back to H 1 i . More precisely, the function is applied on H

H 1 i → H 2 i → H 1 i .
For completeness, we analyze the behavior of each interpolation operator for many interpolation steps. This will point out how error due to solution interpolation accumulates. This is crucial for anisotropic mesh adaptation application where a large number of interpolation steps are done. To this end, transfer

H 1 i → H 2 i → H 1 i is performed five times leading to a total of 10 interpolation steps. It is denoted 5 × H 1 i → H 2 i → H 1 i .
The four studied analytical functions are represented in Figure 5. For each analytical function, a figure is given providing:

• left, relative mass variation

|m H 1 i -m H 2 i | m H 1 i for solution transfer H 1 i → H 2 i • middle, error ε i for solution transfer H 1 i → H 2 i
• and right, error i for solution transfers

H 1 i → H 2 i → H 1 i and 5 × H 1 i → H 2 i → H 1
i for a total of 2 (#I = 2) and 10 (#I = 10) interpolations, respectively.

The linear interpolation scheme is represented by the red (and green) lines and the P 1conservative interpolation is represented by the blue (and purple) lines.

A Gaussian function The first analytical function is a gaussian given by: u 1 (x, y, z) = exp(-30 (x 2 + y 2 + z 2 )) . Step

Inria

This smooth function is representative of the vortices encountered in CFD, Figure 6. The relative mass variation with the classical linear interpolation reduces up to 40, 000 vertices and then stabilizes around 0.005% of variation. Conversely, the relative mass variation with the P 1 -conservative interpolation is of the order of the round-off (≈ 5. 10 -14 ) for all interpolation steps.

As regards the accuracy and the convergence order, both interpolation scheme are converging at order 2 for solution transfers

H 1 i → H 2 i , H 1 i → H 2 i → H 1 i and 5 × H 1 i → H 2 i → H 1 i .
This fits to the theory. We notice that the P 1 -conservative interpolation is more accurate than the linear one in both cases. The difference in accuracy is almost a factor 1.7 and 2.4 for solution transfer H 1 i → H 2 i and H 1 i → H 2 i → H 1 i , respectively. However, analyzing the 5×H 1 i → H 2 i → H 1 i case, we notice that the error level stays almost the same with the P 1 -conservative interpolation when the number of interpolations increases while the error increases with the linear interpolation. It means that error accumulates in the classical case and this adverse effect is drastically reduced with the new interpolation scheme. As a consequence, the error with the linear interpolation is 7.7 times higher than the P 1 -conservative interpolation one.

A continuous sinusoidal shock This analytical function represents a continuous model of a shock which can be assimilated to the numerical capture of a shock with a dissipative flow solver, i.e., the solver captures the shock on several mesh elements. This smooth function is given by: u2(x, y, z) = tanh 20 x + 0.3 sin(-10 y) -0.3 sin(-5(z -0.1)) .

It contains two quasi-constant regions that are separated by sinusoidal interfaces in which strong gradient variation occurs continuously.

As previously, the relative mass variation with the P 1 -conservative interpolation is of the order of the round-off (≈ 10 -14 ) for all interpolation steps. For the linear interpolation the relative mass variation decreases with the mesh size. It varies from 13% for the first meshes to 0.02% for the last meshes with two million vertices.

The P 1 -conservative interpolation achieves an order 2 of convergence for all solution transfers whereas the linear interpolation has a convergence order less than 2. The convergence order asymptotically reaches 1.75 for all cases. As regards the accuracy, the P 1 -conservative interpolation is more accurate than the linear one in all cases and the difference increases while meshes are refined. For solution transfer H 1 i → H 2 i , the error is 2 times smaller with the P 1 -conservative interpolation on the finest meshes. This gap increases for multiple solution transfers and rises to 5.8 when ten interpolations are done.

A multi-scales smooth function This function presents smooth sinusoidal variations but at different scales. There are two order of magnitudes between small and large scales variations. This function reads:

u 3 (x, y, z) =          0.01 sin(200 x y z) if x y z ≤ -π 200 sin(200 x y z) if -π 200 < x y z ≤ 2 π 200 0.01 sin(200 x y z) if 2 π 200 < x y z
.

The relative mass variation with the P 1 -conservative interpolation is of the order of the round-off (≈ 10 -14 ) for all interpolation steps. For the linear interpolation, the mass variation is large from 5% for step 1 to 0.01% for step 5 for only one solution transfer. However, the mass variation seems to converge toward zero while the mesh size goes toward zero.

For this smooth case, the P 1 -conservative approach reaches an order 2 of convergence for all solution transfers. With the linear interpolation an order 2 is reached for one transfer and it is lower than 2 for multiple transfers. Concerning the accuracy, the P 1 -conservative interpolation achieves better accuracy than the standard linear interpolation. The difference in accuracy increases with the number of solution transfers. For the finest meshes (step 4 and 5), the error is 6 times smaller for ten solution transfers while it was only a factor 1.7 for one solution transfer. Again, error accumulates with linear interpolation while the effect is diminished with the conservative one.

F. Alauzet

H 1 i → H 2 i . Middle, error ε i for the transfer H 1 i → H 2 i . Right, error i for the transfer H 1 i → H 2 i → H 1 i (#I=2) and 5 × H 1 i → H 2 i → H 1 i (#I=10).
H 1 i → H 2 i . Middle, error ε i for the transfer H 1 i → H 2 i . Right, error i for the transfer H 1 i → H 2 i → H 1 i and 5 × H 1 i → H 2 i → H 1 i .
Figure 8: Multi-scales smooth analytical function u 3 . Left, mass variation for the transfer A discontinuous function The last analytical function is discontinuous and represents eight steps:

H 1 i → H 2 i . Middle, error ε i for the transfer H 1 i → H 2 i . Right, error i for the transfer H 1 i → H 2 i → H 1 i (#I=2) and 5 × H 1 i → H 2 i → H 1 i (#I=10).
u4(x, y, z) =                       
1 if x ≥ 0 and y ≥ 0 and z ≥ 0 2 if x ≥ 0 and y < 0 and z ≥ 0 3 if x < 0 and y ≥ 0 and z ≥ 0 4 if x < 0 and y < 0 and z ≥ 0 5 if x ≥ 0 and y ≥ 0 and z < 0 6 if x ≥ 0 and y < 0 and z < 0 7 if x < 0 and y ≥ 0 and z < 0 8 if x < 0 and y < 0 and z < 0

.
The solution is constant in eight cubic regions and is discontinuous at the interface of each region.

Again, the relative mass variation with the P 1 -conservative interpolation is of the order of the round-off for all the interpolation steps. For the linear interpolation, the relative mass variation varies from 0.2% for step 1 to 0.0002% for step 5 for one solution transfer. It seems to converge toward zero while the size approaches zero.

Even if the mass is preserved, the same accuracy is obtained for both approaches while transferring the solution from one mesh to another one, Nevertheless, for the multiple solution transfers, the P 1 -conservative interpolation performs better than the classical linear approach as it accumulates less error.

For this purely discontinuous case, the two approaches reach only an order 1 of convergence for all solution transfers.

Conclusions For all those analytical cases, while preserving the mass, the P 1 -conservative interpolation obviously achieves better accuracy than the classical linear interpolation and for some cases it converges at a faster rate. We also notice that it gathers less error when multiple solution transfers are performed. This points out the superiority of the conservative interpolation over the standard one.

Parallelization of the conservative interpolation algorithm

The CPU time overhead is minor in 2D but a major issue in 3D. In 2D, the conservative interpolation is three times slower than the polynomial linear interpolation [START_REF] Alauzet | P1-conservative solution interpolation on unstructured triangular meshes[END_REF]. In the context of 2D anisotropic mesh adaptation, the CPU times devoted to the interpolation stage is very small (less than one percent) relative to the total adaptive CPU time. Therefore, the time overhead induced by the conservative interpolation as compared to the polynomial one can be considered as negligible.

In 3D, the conservative interpolation is 50 times slower on average than the polynomial linear interpolation, meaning that it has non-negligible cost in the adaptive process, but numerical results shown in Section 7.2 will point out that it is well worth it. Let us explain why the overhead is so important. For all the performed test cases, there are 16 computed intersections on average for each tetrahedron. For each intersection, 32 barycentrics are computed F. Alauzet to solve the 48 edge-face intersections. Therefore, a mean of 512 barycentrics evaluations is needed for each tetrahedron which is 128 times the cost of the linear interpolation. Then, we have to generate a 3D mesh of these intersections which is composed of 8 tetrahedra on average.

Fortunately, the procedure is easily parallelized and scales very well because most of the CPU time is spent in step 3 of Algorithm 1. For this step, the intersection of each tetrahedron of the new mesh with the background mesh can be done independently. There is no dependency nor communication which makes it efficient to parallelize. In this work, we follow the strategy proposed in [START_REF] Alauzet | On the use of space filling curves for parallel anisotropic mesh adaptation[END_REF] where space filling curve based renumbering is used to minimize cache misses and memory contention, and the parallelization uses the p-thread paradigm [START_REF]A parallelization framework for numerical simulation[END_REF][START_REF]Handling unstructured meshes in multithreaded environments with the help of Hilbert renumbering and dynamic scheduling[END_REF].

Parallel performance has been analyzed on two different multi-cores computers with different processors and memory access speeds:

• Computer 1:

-2 chips: Xeon E5-2670 10 cores 2.5 GHz both chips are connected by 2 QPI links with a speed of 16 GB/s • Computer 2:

-4 chips: Xeon E7-4850 10 cores 2 GHz all chips are connected to all by 1 QPI link with a speed of 16 GB/s

The selected test case is the spherical blast of Section 7.2.1 where the state (five solution fields) at a-dimensioned time 0.6945 is interpolated. The background mesh size is 12, 993, 399 tetrahedra and the new mesh size is 13, 037, 975 tetrahedra. For that case, 226 millions tetrahedron-tetrahedron intersections have been computed and 1.7 billions tetrahedra have been generated to mesh the intersections.

Parallel timings are compared to the CPU time in serial on the same mesh, thus strong speed-ups are analyzed. In the timings analysis, I/Os and initializations are not taken into account. For each parallel run, hyper-threading has been used by launching a number of threads equal to twice the number of cores. For instance, the 4 HT run means than only 4 cores have been used but 8 threads have been launched. Timings and speed-ups for computer 1 and 2 are summarized in Table 2 and3, respectively.

We notice that the speed-ups are excellent on both computers and even super-linear for a low number of cores thanks to the hyper-threading which reduces memory latency. However, for computer 2, we notice that when more chips are used, the speed-up degrades. This is mainly due to slower memory access between the chips (only one link between each).

We can evaluate the algorithm performance by computing the average time spent to treat one element. On Computer 1, it takes in serial 0.1 ms of CPU time per element and on 20 cores with hyper-threading 0.0055 ms of CPU time per element.

Application to mesh adaptation

One of the main application is time-accurate anisotropic mesh adaptation in CFD. Mesh adaptation provides a way to control the accuracy of the numerical solution by modifying the domain discretization according to size and directional constraints. It is well known that mesh adaptation captures accurately physical phenomena in the computational domain while reducing significantly the CPU time, see [START_REF] Alauzet | 3D transient fixed point mesh adaptation for time-dependent problems: Application to CFD simulations[END_REF][START_REF] Bottasso | Anisotropic mesh adaption by metric-driven optimization[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] Gruau | Coupez, 3D tetrahedral, unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric[END_REF][START_REF] Hecht | Mesh adaptation by metric control for multi-scale phenomena and turbulence[END_REF][START_REF] Pain | Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations[END_REF].

Unsteady mesh adaptation scheme

Our goal is to solve an unsteady PDE which is set in the computational space-time domain Q = Ω × [0, T ] where T is the (positive) maximal time and Ω ⊂ R 3 is the spatial domain. Let Π h be the usual P 1 projector. The considered problem of mesh adaptation consists in Inria finding the space-time mesh H of Q that minimizes the space-time linear interpolation error u -Π h u in L p norm. The problem is thus stated in an a priori way:

Find H opt having N st space-time vertices such that E L p (H opt ) = min H u-Π h u L p (Ω h ×[0,T ]) .
In the continuous mesh framework, we rewrite this problem under the continuous form [START_REF] Loseille | Continuous mesh framework. Part I: well-posed continuous interpolation error[END_REF]:

Find M L p = (M L p (x, t)) (x,t)∈Q such that E L p (M L p ) = min M u -π M u L p (Ω×[0,T ]) , (4)
under the space-time constraint:

C st (M) = T 0 τ (t) -1 Ω d M (x, t) dx dt = N st . (5)
where τ (t) is the time step used at time t. To find the optimal space-time continuous mesh, Problem (4-5) is solved in two steps. First, a spatial minimization is done for a fixed t. Second, a temporal minimization is performed. The expression of the optimal space-time metric M L p for a prescribed time step is [START_REF] Alauzet | Contributions aux méthodes numériques pour l'adaptation de maillage et le maillage mobile[END_REF]:

M L p (x, t) = N 2 3 ST T 0 τ (t) -2p 2p+3 K(t)dt -2 3 τ (t) 2 2p+3 (det |H u (x, t)|) -1 2p+3 |H u (x, t)| , (6)
where H u is the Hessian of sensor u.

The previous analysis provides the optimal size of the adapted meshes for each time level. Hence, this analysis requires the mesh to be adapted at each flow solver time step which is inconceivable. Now, we want to extend the previous analysis to the fixed-point mesh adaptation algorithm context [START_REF] Alauzet | 3D transient fixed point mesh adaptation for time-dependent problems: Application to CFD simulations[END_REF]. The idea consists in splitting the simulation time frame [0, T ] into n adap adaptation sub-intervals:

[0, T ] = [0 = t 0 , t 1 ] ∪ . . . ∪ [t i , t i+1 ] ∪ . . . ∪ [t n adap -1 , t n adap ] , (7)
and to keep the same adapted spatial mesh M i for each time sub-interval [t i-1 , t i ]. On each sub-interval, the mesh is adapted to control the solution accuracy from t i to t i+1 . Consequently, the time-dependent simulation is performed with n adap different adapted meshes. This drastically reduces the number of remeshing during the simulation, hence the number of solution transfers. This can been seen as a coarse adapted discretization of the time axis, the spatial mesh being kept constant for each sub-interval when the global space-time mesh is visualized, thus providing a first answer to the adaptation of the whole space-time mesh.

Algorithm 2 Mesh Adaptation Loop for Unsteady Flows

Initial mesh and solution(H0, S 0 0 ) and set targeted space-time complexity Nst # Fixed-point loop to converge the global space-time mesh adaptation problem For j = 1, n ptf x # Adaptive loop to advance the solution in time on time frame [0, T ]

1. For i = 1, n adap (a) S j 0,i = Interpolate conservatively next sub-interval initial sol. from (H j i-1 , S j i-1 , H j i ); (b) S j i =
Compute solution on sub-interval from pair (S j 0,i , H j i ); (c) |H| j i = Compute sub-interval Hessian-metric from sol. sample (H j i , {S j i (k)} k=1,nk ); EndFor 2. C j = Compute space-time complexity from all Hessian-metrics ({|H| j i }i=1,n adap ); 3. {M j i }i=1,n adap = Compute all sub-interval unsteady metrics (C j , {|H| j i }i=1,n adap ); 4. {H j+1 i }i=1,n adap = Generate all sub-interval adapted meshes ({H j i , M j i }i=1,n adap ); EndFor

F. Alauzet

To converge the non-linear mesh adaptation problem, i.e., converging the mesh-solution couple, a fixed-point mesh adaptation algorithm is used. This is also a way to predict the solution evolution and to adapt the mesh accordingly. The previous error analysis can be also carried out in this context [START_REF] Alauzet | Contributions aux méthodes numériques pour l'adaptation de maillage et le maillage mobile[END_REF]. The time-accurate fixed-point mesh adaptation algorithm is schematized in Algorithm 2. where H, S, H and M denote respectively meshes, solutions, Hessian-metrics and metrics.

Flow solver In all the examples, the flow is modeled by the conservative Euler equations. Assuming that the gas is perfect, inviscid and that there is no thermal diffusion, the Euler equations for mass, momentum and energy conservation read:

∂W ∂t + ∇ • F (W ) = 0 ,
where W = t (ρ, ρu, ρE) is the conservative variables vector and vector F represents the convective operator:

F (W ) = t (ρu, ρuu + pex, ρvu + pey, ρwu + pez, u(ρE + p)) .
We have noted ρ the density, u = (u, v, w) the velocity vector, E = T + u 2 2 the total energy and p = (γ -1)ρT the pressure with γ = 1.4 the ratio of specific heats and T the temperature.

The Euler system is solved by means of a Finite Volume technique on unstructured meshes composed of tetrahedra. The proposed scheme is vertex-centered, achieved a second order accuracy in space and a third order accuracy in time with an explicit Runge-Kutta scheme. More details can be found in [START_REF] Alauzet | Size gradation control of anisotropic meshes[END_REF].

Local remesher The generation of the adapted anisotropic meshes is done using a metricbased adaptive local remeshing strategy [START_REF] Loseille | Adaptive anisotropic simulations in aerodynamics[END_REF] where the surface mesh is adapted conjointly with the volume mesh using local mesh modifications. One main advantage of this method is to be extremely robust. Indeed, if an invalid operation occurs, it is simply rejected. The core of the algorithm uses a unique cavity-based operator [START_REF] Loseille | Serial and parallel mesh modification through a unique cavitybased primitive[END_REF]. Thus, each meshing operateor is equivalent to a node insertion or reinsertion. It is the cavity initialization (definition) which defines the scope of the underlying local mesh modification. The node insertion is based on the extension of the Delaunay kernel to the metric-based anisotropic context. This strategy is very efficient to generate high-quality adapted meshes with a very high level of anisotropy O(1 : 10 6 ) [START_REF] Loseille | Adaptive anisotropic simulations in aerodynamics[END_REF].

Remark: For polynomial interpolation, it has been observed that interpolating physical variable (ρ, u, p) is more accurate than dealing with conservative variable (ρ, ρu, ρE). As ρ and p are interpolated, the positivity is preserved. But, the conservative interpolation works on (ρ, ρu, ρE) and, thus, there is no guarantee that p = (γ -1) ρE -ρ u 2 2 remains positive after the process even if it should be. Indeed, for 3D problems involving strong physics (e.g. blast wave) and complex geometry (e.g. sharp ridges in the city) the positivity issue may occur. This is very rare and concerns a few points over million. In that case, the resulting interpolated field contains negative pressure which is not admissible for the flow solver. A correction, using clipping or averaging with neighboring vertices, is then applied to restore a positive pressure.

Numerical simulations

The goal of the presented examples is to emphasized the benefits in terms of accuracy and convergence of using the P 1 -conservative interpolation instead of classical P 1 interpolation for item 1 (a) in the unsteady mesh adaptation loop presented in Algorithm 2. For each simulation, identical parameters are used for the flow solver and the local remesher, and a control of the interpolation error of the density variable in L 2 -norm is done.

To analyze the accuracy of each simulation and to perform a convergence analysis, we compute the L 1 -norm of the space-time error with respect to the reference solution:

errST = T 0 Ω |u ref (x, t) -u(x, t)|dxdt (8)
≈

n ref adap i adap =1 ∆t    N ref tet (i adap ) i tet =1 |Ki tet | |u ref (Gi tet , ti adap ) -u(Gi tet , ti adap )|   
with the notations:

• The total space-time mesh complexity, i.e., total number of vertices, of a simulation ran on n adap meshes is:

N ST = n adap i adap =1 N ver (i adap ) = n adap × Ñver (9)
where N ver (i adap ) is the number of vertices of the i th i adap adapted mesh and Ñver is the average number of vertices per mesh.

Spherical blast

The first example is a spherical Riemann problem between two parallel walls simulating a blast. Initially, the gas is at rest with density ρ out = 1 and pressure p out = 1 everywhere except in a sphere centered at (0, 0, 0.4) with radius 0.2. Inside the sphere the parameters are ρ in = 1 and p in = 5. For both regions, we have γ = 1.4. The initial pressure jump results in a strong outward moving shock wave, an outward contact discontinuity and an inward moving rarefaction wave. The main feature of the solution are the interactions between these waves. Another significant feature is the development of a low density region in the center of the domain. The solution remains cylindrically symmetric throughout the simulation and is computed until a-dimensioned time T = 0.7.

The adaptive reference solution has been computed using n ref adap = 128 sub-intervals, i.e., the number of adapted meshes used to run the simulation, each mesh having an average size of Ñ ref ver = 1, 712, 282 vertices. The total space-time number of vertices of the reference solution is N ref ST = 219 millions vertices. The final adapted mesh and density solution field are shown in Figure 10.

Four series of adaptive simulations have been run to make the comparison between the two interpolation methods. Two series with the P 1 interpolation and two series with the P 1conservative interpolation by setting a theoretical target mesh size of 120, 000 and 240, 000 vertices on average per sub-interval mesh. To increase the space-time mesh complexity and to perform a convergence analysis, the number of sub-intervals is increased. Notice that when the simulation complexity increased, the number of interpolation stages increases as well. For each series, six simulations have been run with 4, 8, 16, 32, 64 and 128 sub-intervals. The total theoretical space-time mesh size being the number of sub-intervals n adap multiplied by the average complexity N Avg . The resulting mesh sizes for the simulations with 240, 000 vertices on average per sub-interval mesh are presented in Table 4 using P 1 interpolation and in Table 5 using P 1 -conservative interpolation. We notice that for the same set of parameters, we obtain almost the same space-time mesh complexity for the P 1 and the P 1 -conservative interpolation.

F. Alauzet

The improvement in accuracy and convergence of the conservative interpolation with respect to the classical one is clearly pointed out in Figure 14 where the space-time errorgiven by Equation (8) -vs. the space-time complexity -given by Equation ( 9) -is plotted for each simulation.

As larger size adapted meshes are generated with Ñ 240K ver than Ñ 120K ver , it results (as expected) in more accurate solutions, therefore the blue and the pink curves are below the red and the green curves, respectively. This accounts for the discretization error that has been reduced by increasing the mesh size.

Besides, the difference between the blue (resp. red) curve and the pink (resp. green) curve points out the accumulation of the error due to the interpolation stage. For the P 1 interpolation, we observe that the accuracy of the solution degrades more and more with the number of interpolation stages. This degradation becomes visible when 16 or more interpolation steps are done. This is due to solution transfer errors which accumulates at each interpolation and spoils the solution accuracy. As this is not the case with the P 1 -conservative interpolation even if 128 solution transfers are done, the gap in accuracy between the two approaches grows quickly with the number of adaptations (sub-intervals) and becomes quickly significant. As shown in Section 5, the error introduced by the interpolation stage is linked to the mesh size. The smaller h, the smaller the error. This is why the pink curve is shift to the right with respect to the green curve. Hence, the accumulation of interpolation error is much more dramatic with coarser meshes. This significant difference is evident when we observe the final density solutions and adapted meshes, shown in Figures 11 and12, obtained with each method for a theoretical target mesh size of 240, 000 vertices on average per sub-interval mesh and n adap = 128. The solution using the P 1 -conservative interpolation is obviously more detailed with more shock waves present in the flow, less dissipated shock waves and the low density region at the domain center separated by a contact discontinuity is a lot less diffused. In comparison, we notice that both solutions for n adap = 16 -given in Figure 13 -are more similar. Indeed, less error due to the solution transfer has been accumulated. Moreover, we observe a time shift of the shock waves positions with the P 1 interpolation which increases with the number of sub-intervals, while no time shift occurs with P 1 -conservative interpolation. These comments are highlighted in Figure 14 by density solutions extraction along a line for the reference and these four cases solutions. This significant difference is also emphasized by the final adapted mesh size of each simulation. For a given number of sub-intervals and average complexity, we see in Tables 4 and5 that the average number of vertices per mesh in both cases is similar. However, the larger the number of adapted meshes, the larger the final adapted mesh size using the P 1 -conservative interpolation, see Tables 4 and5. A similar behavior is observed for the simulations where 120, 000 vertices on average per sub-interval mesh has been prescribed. This means that information on the solution is lost during the simulation using the P 1 interpolation.

Tables 4 and 5 also report the CPU time distribution between the interpolation stage, the flow solver stage, and the metric computation and the mesh generation stage. All stages are run in parallel using the same number of processors. Table 4 -for the P 1 interpolation -shows that the 3%/63%/33% (Interpolation/Flow/Metric-Mesh) CPU distribution using 4 adapted meshes evolves to a 11%/31%/58% CPU distribution using 128 adapted meshes. Hence, the flow solver CPU part decreases with respect to the interpolation and meshing stages when the number of sub-intervals increases. This change in the CPU distribution is due to the lower number of flow solver time steps performed on each mesh when a larger number of meshes is used to discretize the simulation time frame3 . Using a larger number of meshes increases the number of interpolation steps and the number of generated meshes. In this example, the average number of flow solver time steps per sub-interval is 525 for n adap = 4 and drops to 60 for n adap = 128. For the P 1 -conservative interpolation, see Table 5, the same behavior is observed: the 15%/41%/44% CPU distribution using 4 adapted meshes evolves to a 45%/23%/32% CPU distribution using 128 adapted meshes. As the P 1 -conservative interpolation is more costly in term of CPU time, it takes a larger part of the CPU distribution. An almost balanced distribution between the three stages is obtained for n adap = 32. Similar CPU distributions are obtained for the simulations where 120, 000 vertices on average per sub-interval mesh has been prescribed.

In conclusion, for time-accurate compressible flow simulations with a large number of adapted meshes, the conservative interpolation is mandatory even if it has an important CPU cost in the unsteady mesh adaptation loop.

A blast in a town

The second example is the propagation of a blast in a geometry representing a city plaza. The physics of this problem has a lot more energy than the one proposed in [START_REF] Alauzet | 3D transient fixed point mesh adaptation for time-dependent problems: Application to CFD simulations[END_REF], it results in stronger and more complex physical phenomena. The main feature is related to the random character of the blast wave propagation due to a large number of waves reflexions on the geometry and the interactions between the numerous blast waves. The computational domain size is 85 × 85 × 70 m 3 . Initially, the gas representing the ambient air is at rest with a density ρ out = 1 and p out = 1. To simulate the blast, a high pressure and density region is introduced in a quarter-circle centered at (6.5, 0) with a radius 0.25. In this region, the relevant parameters are ρ in = 10, p in = 25 and u in = 0. For both regions, we have γ = 1.

F. Alauzet

Two series of adaptive simulations have been run to make the comparison between the two interpolation methods. A theoretical target mesh size of 250, 000 vertices on average per sub-interval mesh has been prescribed. To increase the space-time mesh complexity and to perform a convergence analysis, the number of sub-intervals is increased. Thus, the number of interpolation stages increases with the simulation complexity. For each series, six simulations have been run with 4, 8, 16, 32, 64 and 128 sub-intervals. The resulting mesh size for each simulation are presented in Tables 6 and7. We notice that for the same set of parameters, 8)) vs. space-time complexity (Equation ( 9)) for the four series of adaptive simulations from 4 to 128 sub-intervals. In green and pink, adaptive simulations with the P 1 interpolation and, in red and blue, with the P 1 -conservative interpolation. Right, final density solution extraction along a line for the reference simulation and simulations at average complexity 240, 000 for both interpolations with 16 and 128 sub-intervals.

we obtain almost the same space-time mesh complexity for the P 1 and the P 1 -conservative interpolation. Again, when the number of sub-intervals increases, we observe that the size of the final adapted mesh obtained with the P 1 -conservative interpolation is (a lot) larger than the one obtained with P 1 interpolation. It denotes a loss of accuracy during the simulation while using the P 1 interpolation.

The space-time error vs. the space-time complexity for each simulation is plotted in Figure 15. Final density solutions and corresponding adapted meshes are shown in Figures 17 and18. Comments on the results are identical to the previous example. Again, for a large number of adaptation, the difference in accuracy between the two methods is significant. Thus, the conservative interpolation is mandatory to avoid the accumulation of errors due to the solution transfer which spoils considerably the solution accuracy.

Tables 6 and7 report the CPU time distribution between the interpolation stage, the flow solver stage, and the metric computation and the mesh generation stage. As previously, all stages are run in parallel using the same number of processors. For the P 1 interpolation, see Table 6, the 3%/60%/37% (Interpolation/Flow/Metric-Mesh) CPU distribution using 4 adapted meshes evolves to a 10%/51%/38% CPU distribution using 128 adapted meshes. Therefore, the same observations as Section 7.2.1 can be made. In this example, the average number of flow solver time steps per sub-interval is 750 for n adap = 4 and drops to 200 for n adap = 128. For the P 1 -conservative interpolation, see Table 7, the same behavior is observed: the 12%/55%/33% CPU distribution using 4 adapted meshes evolves to a 34%/37%/29% CPU distribution using 128 adapted meshes. This time, an almost balanced distribution between the three stages is obtained for n adap = 128.

Conclusion

In this work, we have proposed a matrix-free P 1 -conservative interpolation operator that satisfies the maximum principle. This operator is based on local mesh intersections and local operations that make it memory efficient and easy to parallelize. The properties of this new operator have been verified numerically on analytical examples and adaptive simulations. These examples also point out a significative improvement in accuracy and numerical convergence obtained with the conservative interpolation as compared to the classical one when a large number of adapted meshes are considered. The accuracy of the solution is not spoiled by this stage. Consequently, for long-time simulations, the conservative interpolation is mandatory as it is only slightly sensitive to the increase of the number of interpolations.

The proposed conservative interpolation scheme can be extend to P k -representation of the Table 7: City blast. Statistics for the unsteady mesh adaptation algorithm for different number of sub-intervals n adap using the P 1 -conservative interpolation. From left to right, average number of vertices per sub-interval Ñver and final meshes number of vertices N ver (n adap ), percentage of the total CPU time for the interpolation (Itp), the flow solver (Flow) and the metric computation and mesh generation (Mesh) stages. The total space-time mesh number of vertices is n adap × Ñver . solution if solutions are defined at the elements. In that case, a high-order gauss quadrature is considered to recover the mass of the solution field and all the derivatives (up to order k). However, the extension to P k -representation of the solution defined at vertices requires more work.

Inria

Figure 1 :

 1 Figure 1: Left, numeration convention for a tetrahedron K. Right, the 15 regions defined by the signs of the four barycentric coordinates of a point p with respect to element K.

4 .

 4 Correct the gradient to enforce the maximum principle 5. Set the solution values to vertices by an averaging procedure

Figure 3 :

 3 Figure3: All possible intersection configurations for coplanar pair edge-face (all kinds of intersection are gathered by edge-face configurations). Four configurations can occur. One intersection with an edge's vertex in the face and the other vertex out. Two intersections with zero (resp. one or two) edge's vertex in the face and two (resp. one or zero) edge's vertices outside of the face. In the pictures, an intersection marked by "i" or "d" represents a pure or a degenerated intersection, respectively.

Figure 4 :

 4 Figure 4: Two tetrahedron-tetrahedron intersection configurations where the cyan tet intersects the maroon tet. Left, a non-degenerated intersection and, right, a degenerated case where both tets share a common vertex. For each configuration, the mesh of the intersection is depicted. Meshes are composed of 12 (left) and 8 (right) tets, respectively.

  h . Therefore, the volume of each mesh differs:|H back | = |H new |.For instance, when dealing with complex geometries, the volume of the domain may change due to the change in the geometric approximation of the surface between both meshes.

Figure 5 :

 5 Figure 5: Representation of the four analytical functions. From left to right, gaussian function u 1 , continuous sinusoidal shock function u 2 , multi-scales smooth function u 3 and discontinuous function u 4 .

Figure 6 :

 6 Figure 6: Gaussian analytical function u 1 . Left, mass variation for the transfer H1 i → H 2 i . Middle, error ε i for the transfer H 1 i → H 2 i . Right, error i for the transfer H 1 i → H 2 i → H 1

Figure 7 :

 7 Figure 7: Continuous sinusoidal shock analytical function u 2 . Left, mass variation for the transferH 1 i → H 2 i . Middle, error ε i for the transfer H 1 i → H 2 i . Right, error i for the transfer H 1 i → H 2 i → H 1 i and 5 × H 1 i → H 2 i → H 1 i .

Figure 9 :

 9 Figure 9: Discontinuous analytical function u 4 . Left, mass variation for the transfer H 1i → H 2 i . Middle, error ε i for the transfer H 1 i → H 2 i . Right, error i for the transferH 1 i → H 2 i → H 1 i (#I=2) and 5 × H 1 i → H 2 i → H 1 i (#I=10).

4 .

 4 The solution is computed until a-dimensioned time T = 15. The adaptive reference solution has been computed using n ref adap = 128 adapted meshes, each mesh having an average size of Ñ ref ver = 3, 327, 382 vertices. The total space-time number of vertices of the reference solution is N ref ST = 426 millions vertices. The final adapted mesh and density solution field are shown in Figure 16.

  n

Figure 10 :

 10 Figure 10: Spherical blast. Reference density solution (left) at a-dimensioned time T = 0.7 and final adapted reference mesh (right) containing 2, 173, 612 vertices and 13, 037, 975 tetrahedra.

Figure 11 :

 11 Figure 11: Spherical blast. Density solution (left) at a-dimensioned time T = 0.7 and final adapted mesh (right) containing 367, 101 vertices and 2, 177, 486 tetrahedra obtained with the P 1 -conservative interpolation for n adap = 128 and Ñver = 240K.

Figure 12 :

 12 Figure 12: Spherical blast. Density solution (left) at a-dimensioned time T = 0.7 and final adapted reference mesh (right) containing 259, 879 vertices and 1, 535, 198 tetrahedra obtained with the P 1 interpolation for n adap = 128 and Ñver = 240K.

Figure 13 :Figure 14 :

 1314 Figure 13: Spherical blast. Density solution at a-dimensioned time T = 0.7 obtained with the P 1 -conservative interpolation (left) and the P 1 interpolation (right) for n adap = 16 and Ñver = 240K.

Figure 15 :

 15 Figure 15: City blast. Space-time error (Equation (8)) vs. space-time complexity (Equation (9)) for the adaptive simulations with the P 1 interpolation (green) and the P 1 -conservative interpolation (red) from 4 to 128 sub-intervals.

Figure 16 :

 16 Figure 16: City blast. Reference density solution (left) at a-dimensioned time T = 15 and final adapted reference mesh (right) containing 4, 187, 548 vertices and 25, 249, 618 tetrahedra.

Figure 17 :

 17 Figure 17: City blast. Density solution (left) at a-dimensioned time T = 15 and final adapted mesh (right) containing 487, 440 vertices and 2, 802, 472 tetrahedra obtained with the P 1 -conservative interpolation for n adap = 128.

Figure 18 :

 18 Figure 18: City blast. Density solution (left) at a-dimensioned time T = 15 and final adapted mesh (right) containing 314, 963 vertices and 1, 788, 719 tetrahedra obtained with the P 1 interpolation for n adap = 128.

continuous or discontinuous) representation of the solution on H back 1. Localize all vertices of H new in H back 2. ∀K back ∈ H back , compute solution mass m K back and gradient ∇ K back 3

  2. . ∀K new ∈ H new , recover solution mass m K new and gradient ∇ K new : K new and ∇ K new using Gauss quadrature formulae =⇒ get a piecewise linear discontinuous representation of the mass on H new

	Algorithm 1 Conservative Interpolation Process
	Piecewise linear ((a) compute the intersection of K new with all K back i (b) mesh the intersection polyhedron of each pair (K new , K back ∈ H back it overlaps i )
	(c) compute m

  The error i is obtained by computing the gap in L 1 -norm between u 1 i and Πu 1

	1 i giving u 1 i , i to H 1 i is interpolated from H 2 i and finally Πu 2 i giving Πu 2 then it is interpolated on H 2 i and we obtain Πu 1 i .

i on H 1 i . This double solution transfer is called transfer

Table 1 :

 1 Mesh sizes of the series of embedded meshes (H 1

		1	2	3	4	5
	# vertices H 1 i # vertices H 2 i	906 6 287 42 656 322 656 2 330 670 918 6 432 46 695 360 543 2 739 492

i ) i=1...5 and (H 2 i ) i=1...5 .

Table 2 :

 2 Timings and speed-up on computer 1 up to 20 cores with hyper-threading.

	1 i →

i . Middle, error ε i for the transfer H 1 i → H 2 i . Right, error i for the transfer

H 1 i → H 2 i → H 1 i (#I=2) and 5 × H 1 i → H 2 i → H 1 i (#I=10).

Table 3 :

 3 Timings and speed-up on computer 2 up to 40 cores with hyper-threading.

  adap ) is the number of tetrahedra of the i th i adap adapted mesh used to compute the reference solution for sub-interval[t i adap -1 , t i adap ] • |K itet | is the volume of the i th itet tetrahedron • u ref (G itet , t i adap )and u(G itet , t i adap ) are the reference solution and the solution at i th itet tetrahedron barycenter at time t i adap .

	n ref adap is the number of reference adapted meshes used for the simulation
	• ∆t = T adap n ref	is the sub-intervals time length
	• N ref tet (i

Table 4 :

 4 Spherical blast. Statistics for the unsteady mesh adaptation algorithm for different number of sub-intervals n adap using the P 1 interpolation. From left to right, average number of vertices per sub-interval Ñver and final meshes number of vertices N ver (n adap ), percentage of the total CPU time for the interpolation (Itp), the flow solver (Flow) and the metric computation and mesh generation (Mesh) stages. The total space-time mesh number of vertices is n adap × Ñver .

	n adap			P1-conservative
		Ñ 240K ver	N 240K ver (n adap ) Itp CPU Flow CPU Mesh CPU
	4	182, 330	206, 760	14.98% 41.12%	43.90%
	8	199, 328	229, 657	23.62% 53.70%	22.68%
	16	217, 308	259, 689	30.44% 43.94%	25, 62%
	32	236, 141	291, 980	35.77% 35.49%	28.74%
	64	256, 758	327, 466	41.11% 28.14%	30.75%
	128	279, 709	367, 101	44.80% 23.05%	32.15%

Table 5 :

 5 Spherical blast. Statistics for the unsteady mesh adaptation algorithm for different number of sub-intervals n adap using the P

1 -conservative interpolation. From left to right, average number of vertices per sub-interval Ñver and final meshes number of vertices N ver (n adap ), percentage of the total CPU time for the interpolation (Itp), the flow solver (Flow) and the metric computation and mesh generation (Mesh) stages. The total space-time mesh number of vertices is n adap × Ñver .

Table 6 :

 6 City blast. Statistics for the unsteady mesh adaptation algorithm for different number of sub-intervals n adap using the P 1 interpolation. From left to right, average number of vertices per sub-interval Ñver and final meshes number of vertices N ver (n adap ), percentage of the total CPU time for the interpolation (Itp), the flow solver (Flow) and the metric computation and mesh generation (Mesh) stages. The total space-time mesh number of vertices is n adap × Ñver .

	n adap			P1-conservative
		Ñ 120K ver	N 120K ver (n adap ) Itp CPU Flow CPU Mesh CPU
	4	210, 309	241, 620	12.02% 54.57%	33.41%
	8	230, 148	287, 236	16.76% 56.11%	27.13%
	16	247, 354	331, 257	19.56% 58.00%	22.44%
	32	270, 084	384, 186	23.69% 51.87%	24.44%
	64	296, 473	438, 169	28.20% 45.04%	26.76%
	128	323, 835	487, 440	33.72% 36.82%	29.46%

RR n°8785

Ensuring the robustness of the local mesh intersection and the meshing of these intersections is considerably more difficult in 3D than

in 2D. Indeed, volume positivity is harder to satisfy than area positivity when dealing with degenerated cases.

Inria

A method based on the notion of slope limiter widely used in numerical schemes has also been implemented but it turns out to be less accurate than the current method on all validation cases.

Indeed, more meshes considered to reach the simulation final time (see Expression[START_REF] Farrell | Conservative interpolation between unstructured meshes via supermesh construction[END_REF]) means less time steps performed on each mesh.

InriaRESEARCH CENTRE PARIS -ROCQUENCOURTDomaine de Voluceau, -Rocquencourt B.P. 105 -78153 Le Chesnay Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399