Franck Boyer
email: fboyer@cmi.univ-mrs.fr

Celine Lapuerta
email: celine.lapuerta@irsn.fr

Sebastian Minjeaud
email: sebastian.minjeaud@irsn.fr

Bruno Piar
email: bruno.piar@irsn.fr

A LOCAL ADAPTIVE REFINEMENT METHOD WITH MULTIGRID PRECONDITIONNING ILLUSTRATED BY MULTIPHASE FLOWS SIMULATIONS

The aim of this paper is to describe some numerical aspects linked to incompressible three-phase flow simulations, thanks to Cahn-Hilliard type model. The numerical capture of transfer phenomenon in the neighborhood of the interface require a mesh thickness which become crippling in the case where it is applied to the whole computational domain. This suggests the use of a local refinement method which allows to dynamically focus on problematic areas. The notion of refinement pattern, introduced for Lagrange finite elements, allows to build a conceptual hierarchy of nested conformal approximation spaces which is then used to implement the so-called CHARMS local refinement methods. Properties of these methods are proved ensuring in particular the conformity of approximation spaces at every time of simulations. Furthermore, the multilevel structure obtained by this method, is used to construct multigrid preconditioners. Finally, after a validation on a model problem, the performance of the whole method is illustrated on an example of a liquid lens spreading between two stratified fluids.

Résumé. L'objectif de l'article est de décrire certains aspects numériques liés à la simulation d'écoulements incompressibles à trois phases non miscibles, à l'aide de modèles à interfaces diffuses de type Cahn-Hilliard. La capture numérique des phénomènes de transfert au voisinage des interfaces requiert une finesse de maillage qui devient rédhibitoire si elle est appliquée à l'ensemble du domaine de calcul. Ceci suggère l'utilisation de méthodes de raffinement local adaptatif qui permettent de se focaliser dynamiquement sur les zones sensibles. La notion de motif de raffinement, introduite pour des éléments finis de Lagrange, permet de construire une hiérarchie conceptuelle d'espaces d'approximation conformes emboîtés qui est alors utilisée pour mettre en oeuvre les techniques de raffinement local dites CHARMS. Les propriétés de la méthode sont prouvées assurant en particulier la conformité de l'espace d'approximation à tout instant des simulations. En outre, la structure multiniveaux obtenue par cette méthode est exploitée pour construire des préconditionneurs multigrilles. Enfin, après une validation sur un problème modèle, les performances de l'ensemble sont illustrées sur un exemple d'étalement d'une lentille piégée entre deux phases stratifiées.

Introduction

During a hypothetical major accident in a nuclear reactor, the degradation of the core may produce multiphase flows where interfaces undergo extreme topological changes, e.g. break-up and coalescence. Because of their ability to capture interfaces implicitly, phase field models are attractive for the numerical simulation of such phenomena [START_REF] Boyer | Study of a three component Cahn-Hilliard flow model[END_REF]. They replace sharp interfaces by thin but nonzero thickness transition regions where capillary effects and changes in fluid properties are smoothed. The implied burden on the numerical method lies in the required resolution of very thin moving internal layers. We present here a general adaptive local refinement strategy and an associated multigrid solver in order to tackle such issues. We illustrate the method on the ternary Cahn-Hilliard model taken from [START_REF] Boyer | Study of a three component Cahn-Hilliard flow model[END_REF].

The challenge of local refinement can be expressed in that way: increase the spatial resolution of some part of the domain which are dynamically selected at each time step. Of course, the choice of these parts is not obvious, nevertheless, for the description of the general method, we assume that a refinement criterion is available so that we know the area to refine. We will give, in Section 3.3.3, the precise criterion we used for the resolution of the Cahn-Hilliard model.

To increase the spatial resolution in these selected areas, a solution is to use elements with smaller diameters, i.e. to split some cells into smaller ones. The main difficulty is to preserve, at the same time, the geometric conformity of the mesh and then its "good" geometric quality. The geometric conformity prevents the so-called hanging nodes which for different reasons are undesired in many applications. For instance, hanging nodes do not represent degree of freedom and are somewhat difficult to handle, because the local correlation pattern of the stiffness matrix is disturbed. When they exist, their taking into account may be carried out by many ways, for instance, by direct elimination of these "false" unknowns, or by adding constraint and using penalty methods or Lagrange multiplier methods. In these cases, numerical methods and schemes are modified.

Another method consists in eliminating non-conformity of the mesh by splitting cells until there is no nonconformal edge left. In two dimensions, for triangular meshes, we may for instance cite the red-green refinement method [START_REF] Bank | Some refinement algorithms and data structures for regular local mesh refinement[END_REF]. This technique consists in first using a "regular" splitting (called red refinement) of triangles in four congruent ones by connecting the midpoints of its edges. This refinement preserves geometric properties of triangles but creates non-conformal edges when a refined triangle and an unrefined one are adjacent. Therefore, a second type of refinement (called green refinement) is first used connecting one edge midpoint and the opposite corner. This yields a simple bisection which is "irregular" but only used for solving conformity issues. Bey [START_REF] Bey | Tetrahedral grid refinement[END_REF] and Zhang [START_REF] Zhang | Multi-level iteratives techniques[END_REF] proposed a generalisation of this method to three dimensions. Other methods based only on bisection have been introduced by Rivara [START_REF] Rivara | Mesh refinement processes based on the generalized bisection of simplices[END_REF][START_REF] Rivara | The 4-triangles longest-side partition of triangles and linear refinement algorithms[END_REF] or Mitchell [START_REF] Mitchell | Adaptive refinement for arbitrary finite-element spaces with hierarchical bases[END_REF] in two dimension and Bänsch [START_REF] Bänsch | Local mesh refinement in 2 and 3 dimensions[END_REF] or Maubach [START_REF] Maubach | Local bisection refinement for n-simplicial grids generated by reflection[END_REF] in three dimensions. All these methods depend on the choice of the element and on the dimension. Moreover, their implementation tends to be quite complex, notably in three dimensions.

A possible alternative considered in this paper is to adopt the point of view of basis functions refinement instead of cells refinement. In this approach [START_REF] Krysl | Natural hierarchical refinement for finite element methods[END_REF], a stock of basis functions with increasing spatial resolution is assumed to be given through a nested sequence of approximation spaces X 0 ⊂ • • • ⊂ X J , J

1. Local refinement is then performed by using multilevel approximation spaces containing few basis functions suitably selected in each approximation spaces X j , j ∈ 0, J . In this framework, Krysl, Grinspun and Schröder in [START_REF] Krysl | Natural hierarchical refinement for finite element methods[END_REF] (see also [START_REF] Grinspun | CHARMS: A simple framework for adaptive simulation[END_REF][START_REF] Hammon | Implementation of a general mesh refinement technique[END_REF][START_REF] Krysl | Object-oriented hierarchical mesh refinement with CHARMS[END_REF]) proposed procedures called CHARMS (Conforming Hierarchical Adaptive Refinement MethodS) which enable to (un)refine multilevel approximation spaces. We give here a precise construction of the nested sequence X 0 ⊂ • • • ⊂ X J for suitable Lagrange finite element (e.g. P k , Q k , k 1) using a unique refinement pattern and we study a slightly modified version of the quasi-hierarchical CHARMS adaptation technique (the "one-level-difference refinement rule" is used as a separate criterion, and is not involved in the proof of the properties of the procedure). Finally, we incorporate multigrid preconditioners [START_REF] Bank | The hierarchical basis multigrid method[END_REF][START_REF] Bramble | Parallel multilevel preconditioners[END_REF][START_REF] Hackbusch | Multigrid Methods and Applications[END_REF][START_REF] Yserentant | Hierarchical bases[END_REF][START_REF] Zhang | Multi-level iteratives techniques[END_REF] in the method. Attractive features follow from this methodology:

• there is no modification of the discrete problem due to the mesh adaptation,

• the cells are divided into cells of the same type, uniformly applying the same subdivision pattern, • the possible geometric non-conformity of the adapted meshes are implicitly handled, • there is no specific treatment due to particular Lagrange finite elements, e.g.

P k , Q k , k 1,
• all the procedure is independent of the space dimension,

• transfer operators between the different grids are not needed for time evolution problems provided that a suitable definition of elementary integration domains is used when assembling the system.

The detailed outline of the paper is the following. In Section 1, we introduce the notion of refinement pattern for suitable Lagrange elements (e.g. P k and Q k ) and we show how a hierarchy of nested conforming meshes can be built recursively by applying a unique refinement pattern to each cell of a given initial mesh. We investigate in details the parent-child relationship between basis functions of successive refinement levels. In particular, we show that all the coefficients in this linear relationship can be simply computed by considering the refinement pattern on the reference element. Consequently, the meshes hierarchy is never explicitly created since every necessary information is available on the reference element.

In Section 2, we describe the local adaptation procedure in the above framework and we establish its main properties. In particular, we give precise sufficient conditions ensuring that:

• the selected basis functions on the various levels are always linearly independent.

• no information is lost when refining a basis function: the refinement algorithm produces increasing approximation spaces sequences. • the approximation space obtained by refining (resp. unrefining) a set of basis functions is independent of the order in which successive refinements (resp. unrefinements) are performed.

Some counterexamples illustrate the fact that these properties are not satisfied in the general case.

In the second part of Section 2, we show that the multilevel structure of the approximation spaces built following the above methodology can be exploited to derive efficient multigrid preconditioners. All the results in this section are illustrated by numerical results on a model problem.

Section 3 is devoted to illustrations to more complicated partial differential equations. We concentrate on the simulation of three-phase systems thanks to a diffuse-interface model of Cahn-Hilliard type. The model is a non-linear parabolic problem describing the evolution of thin interfaces. This is the reason why such a local adaptive approach is particularly fruitful. Interestingly enough, we focus on the fact that, even though the approximation spaces may vary from a time step to another, it is not needed to make use of transfer operators between the different grids when assembling the matrix of the system. In the conclusion of the paper, we give even more complex simulation results we obtained by solving the complete three-phase flows model constituted by a coupling between the Cahn-Hilliard system and the Navier-Stokes equations, which is the main motivation of this study.

Multilevel Finite Element Spaces

Preliminary notation and definitions

This subsection introduces some classical notation, definitions and properties, following [START_REF] Ern | Theory and Pratice of Finite Elements[END_REF] and [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF], about meshes and Lagrange finite elements. These would be useful in Sections 1.2 and 1.3.

Definition 1.1 (Lagrange Finite Element [20, §4.1]). A Lagrange finite element is a triple (K, Σ, P ) where • K is a compact, connected, Lipschitz subset of R d (d = 1, 2 or 3), • Σ = {a k ∈ K; 1 k N } is a set of N distinct points belonging to K, called the Lagrange nodes, • P is a vector space of functions p : K → R such that Σ is P -unisolvent [20, Def 4.1-1], i.e. ∀(α 1 , . . . , α N ) ∈ R N , ∃!p ∈ P, ∀k ∈ 1, N , p(a k ) = α k .
The element (K, Σ, P ) is called a polygonal Lagrange finite element iff K is a polygon.

• ω = Ne e=1
K e and,

• Ke ∩ Kf = ∅ for any pair of distinct cells (K e , K f ).

Definition 1.3 (Finite element mesh generated by a reference element [9, §1.3.2] ). Let T = {K e ; 1 e N e } be a mesh of a domain ω. Let ( K, Σ, P ) be a polygonal Lagrange finite element, hereafter called reference element. We say that the mesh T is generated using the reference element ( K, Σ, P ) iff:

For all e ∈ 1, N e , there exists a C 1 -diffeomorphism T e from K to K e .

When the transformations T e , 1 e N e , are affine, the mesh is said to be affine.

For all e ∈ 1, N e , we can define a Lagrange finite element (K e , Σ e , P e ) by setting:

• Σ e = T e ( Σ) and

• P e = {p • T -1 e ; p ∈ P }. Remark 1.4. We always assume that the reference element is polygonal. However, depending on the mapping T e , a cell K e may be non polygonal [START_REF] Ern | Theory and Pratice of Finite Elements[END_REF]Figure 1.13].

The following definitions are commonly used in order to build H 1 (ω)-conformal finite element approximation spaces X, i.e. such that X ⊂ H 1 (ω), see Proposition 1.10.

Definition 1.5 (C 0 -class element [20, Def 5.1-2]). The polygonal Lagrange finite element ( K, Σ, P ) is a C 0 -class element iff • P ⊂ C 0 ( K) and • for any face F of K, Σ ∩ F is P | b F -unisolvent where P | b F = { ϕ | b F ; ϕ ∈ P }.
Definition 1.6 (Compatibility requirements [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF]). The compatibility requirements for the Lagrange reference element ( K, Σ, P ) hold iff for all faces F 1 and F 2 of K, for all affine invertible function A such that F 2 = A( F 1 ), we have,

• Σ ∩ F 2 = A( Σ ∩ F 1 ) and • { ϕ | b F1 ; ϕ ∈ P 1 } = { ϕ • A | b F1 ; ϕ ∈ P 2 }.
Definition 1.7 (Geometrically conformal mesh [9, Def 1.55]). A mesh T = {K e ; 1 e N e } of a domain ω is said to be geometrically conformal iff for all cells K e and K f having a non-empty (d-1)-dimensional intersection, say F = K e ∩ K f , T -1 e (F ) and T -1 f (F ) are faces of K, and there exists a bijective affine transformation

A : T -1 e (F ) → T -1 f (F ) such that A • T -1 e = T -1 f on F .
Remark 1.8. Definition 1.7 implies, in particular, that for any pair of distinct cells (K e , K f ), the intersection K e ∩ K f is:

• either empty or a common vertex in dimension 1,

• either empty, or a common vertex, or a common face in dimension 2,

• either empty, or a common vertex, or a common edge, or a common face in dimension 3. An example of a geometrically non-conformal mesh is shown in Figure 1. Remark 1.9. In case that the mesh is geometrically conformal, Definition 1.6 implies, in particular, that the Lagrange nodes on a common face belong to each element sharing the face. An example of incompatible node positions is given in Figure 2.

T f K e T e Σ e = { } Σ f = { } K f K Figure 2. Example of incompatible node positions.
Let T = {K e ; 1 e N e } be a mesh generated using a Lagrange reference finite element ( K, Σ, P ). Owing to Definition 1.3, we can associate to this mesh:

• a set of geometric mappings {T e : K → K e ; 1 e N e }.

• a set of Lagrange finite elements {(K e , Σ e , P e ); 1 e N e }.

• a set of Lagrange nodes Σ = Furthermore, since Σ is P -unisolvent, we have:

∀k ∈ 1, N , ∃ ϕ k ∈ P , ∀ℓ ∈ 1, N , ϕ k ( a ℓ ) = δ kℓ .
We can then associate to each node a i , 1 i N dof , a basis function ϕ i defined elementwise by

ϕ i|K e = ϕ k • T -1 e if there exists k ∈ 1, N such that I(e, k) = i 0 otherwise.
Owing to this definition, we denote by supp[ϕ i ] the following subset of ω:

supp[ϕ i ] = e∈E K e where E = e ∈ 1, N e ; ∃k ∈ 1, N , i = I(e, k) . • a H 1 (ω)-conformal approximation space [9, Prop 1.74] X = v ∈ C 0 (ω); ∀e ∈ 1, N e , v |Ke ∈ P e .
The advantage of this construction is not only to produce a H 1 (ω)-conformal approximation space but most important to provide an explicit basis of this space. Indeed, we have the following result: Proposition 1.10. Let ( K, Σ, P ) be a C 0 -class polygonal Lagrange reference element satisfying the compatibility requirements 1.6. Let T = {K e ; 1 e N e } be a geometrically conformal mesh of ω generated using the reference element ( K, Σ, P ). Then, we have ϕ k (a ℓ ) = δ kℓ , and

ϕ 1 , . . . , ϕ N dof is a basis of X. Proof. See [9, Prop 1.78]

Refinement pattern

Let us now define the notion of refinement pattern and the associated compatibility requirements that will be useful to generate geometric conformal uniformly refined meshes in Section 1.3.1.

Definition 1.11 (Refinement pattern). Let

• ( K, Σ, P ) be a C 0 -class polygonal Lagrange reference finite element satisfying the compatibility requirements 1.6,

• T = { K [1]
e ; 1 e N [START_REF] Bank | The hierarchical basis multigrid method[END_REF] e } be a geometrically conformal affine mesh of the interior of K generated using the reference element ( K, Σ, P ) itself.

Then, we say that ( K, Σ, P , T ) is a refinement pattern.

Owing to the previous subsection let us denote by

• { T e : K → K [1]
e ; 1 e N [START_REF] Bank | The hierarchical basis multigrid method[END_REF] e } the set of geometric mappings,

• {( K [1] e , Σ [1]
e , P [START_REF] Bank | The hierarchical basis multigrid method[END_REF] e ); 1 e N [START_REF] Bank | The hierarchical basis multigrid method[END_REF] e } the set of Lagrange finite elements,

• Σ [1] = { a [START_REF] Bank | The hierarchical basis multigrid method[END_REF] k ; 1 k N [1] } the set of Lagrange nodes, • { ϕ [START_REF] Bank | The hierarchical basis multigrid method[END_REF] k ; 1 k N [1] } the set of basis functions, associated to the mesh T . Definition 1.12 (Compatibility requirements). We say that a refinement pattern ( K, Σ, P , T ) satisfies the compatility requirements iff

• for all e ∈ 1, N [1] e , {ϕ | b K [1] e ; ϕ ∈ P } ⊂ P [1] e ,
• for all faces F 1 , F 2 of K, for all affine invertible function A such that F 2 = A( F 1 ), for all faces

F [1] 1 ⊂ F 1 of an element K [1] e , A( F [1]
1 ) ⊂ F 2 is exactly a face of an other element K [START_REF] Bank | The hierarchical basis multigrid method[END_REF] f , • any node of Σ is also a node of Σ [1] , that is Σ ⊂ Σ [1] .

Remark 1.13. The Definition 1.12 is used to avoid refinement patterns which would lead in the sequel to non-conformal meshes. An example is given with Figure 3. The compatibility requirements of refinement patterns 1.12 are also used to etablish the so-called refinement equation on the reference element which is the cornerstone of CHARMS method [START_REF] Krysl | Natural hierarchical refinement for finite element methods[END_REF]. Proposition 1.14 (Refinement equation for the refinement pattern). Let ( K, Σ, P , T ) be a refinement pattern satisfying the compatibility requirements 1.12. We have the following relationship ∀k ∈ 1, N , ϕ k = b N [1] ℓ=1

K K [1] 4 K [1] 3 K [1] 1 K [1] 2 K e K f
β kℓ ϕ [1] ℓ where β kℓ = ϕ k ( a [1] ℓ ).
Proof. Owing to Proposition 1.10, { ϕ

[1] 1 , . . . , ϕ [1] b N [1] } is a basis of {v ∈ C 0 ( K); ∀e ∈ 1, N [1] e , v | b K [1] e ∈ P [1] e }.
However, for all k ∈ {1, . . . , N }, the basis function ϕ k is in C 0 ( K) and using compatibility requirements 1.12, we have ∀e ∈ 1, N

[1] e , ϕ k | b K [1] e ∈ P [1]
e . Hence, there exist coefficients β kℓ such that

ϕ k = b N [1] ℓ=1 β kℓ ϕ [1] ℓ .
Finally, the coefficients β kt can be obtained thanks to the relationship ϕ

[1] ℓ ( a [1] t ) = δ ℓt .
Notations of this section are illustrated by Figures 4 and5 which show the complete refinement patterns and all the refinement equations of .Q 1 -square and P 1 -triangle elements.

Since the geometrical configuration is more complicated for the Q 2 -square element, we only give in Figure 6 the non-zero coefficients in the refinement equation associated to three coarse basis functions. More precisely, on each picture, the coarse basis function represented by a black bullet is the linear combination of the finer basis functions represented by circles with the coefficients mentionned nearby. The other refinement equations can be readily obtained by symmetry.

a [0] 1 K a [1] 0 a [1] 1 a [1] 2 a [1] 3 a [1] 4 a [1] 5 a [1] 6 a [1] 7 a [1] 8 K [1] 2 K [1] 1 K [1] 3 K [1] 4 T 2 a [0] 3 a [0] 2 a [0] 0 ϕ 0 = ϕ [1] 0 + 1 2 ϕ [1] 1 + 1 2 ϕ [1] 3 + 1 4 ϕ [1] 4 ϕ 1 = ϕ [1] 2 + 1 2 ϕ [1] 1 + 1 2 ϕ [1] 5 + 1 4 ϕ [1] 4 ϕ 2 = ϕ [1] 8 + 1 2 ϕ [1] 5 + 1 2 ϕ [1] 7 + 1 4 ϕ [1] 4 ϕ 3 = ϕ [1] 6 + 1 2 ϕ [1] 3 + 1 2 ϕ [1] 7 + 1 4 ϕ [1] 4 Figure 4
. Q 1 -square refinement pattern and refinement equations.

a [1] 1 a [0] 0 a [0] 1 a [0] 2 a [1] 0 a [1] 2 a [1] 4 a [1] 5 K [1] 4 K T 4 K [1] 2 K [1] 3 a [1] 3 K [1] 1 ϕ 0 = ϕ [1] 0 + 1 2 ϕ [1] 1 + 1 2 ϕ [1] 3 ϕ 1 = ϕ [1] 2 + 1 2 ϕ [1] 1 + 1 2 ϕ [1] 4 ϕ 2 = ϕ [1] 5 + 1 2 ϕ [1]
3 + 1 2 ϕ

[1] 4

Figure 5. P 1 -triangle refinement pattern and refinement equations.

-

Multilevel finite element approximation spaces

Let Ω be a bounded domain of R d , d = 1, 2 or 3. The purpose of this section is to give an automatic way to construct H 1 (Ω)-conformal multilevel finite element approximation spaces from an initial geometrically conformal mesh T 0 and a given refinement pattern. Let J ∈ N * . First, we construct a hierarchy of nested meshes T 0 , T 1 , . . . , T J from T 0 by applying uniformly and recursively the refinement pattern. Then, a multilevel approximation space is obtained by selecting some basis functions associated to each mesh T j , 0 j J in a way that guarantees linear independence of the selected basis functions. 1.3.1. Hierarchy of nested H 1 (Ω)-conformal approximation spaces. Parent-child relationship Let ( K, Σ, P , T ) be a refinement pattern and j ∈ 0, J -1 . In this subsection, we assume that a geometrically conformal mesh T j = {K [j] e ; 1 e N [j] e } of Ω is given, that this mesh T j is generated using the reference element ( K, Σ, P ), and then we explain how we can build the mesh T j+1 .

In the sequel, all mathematical objects associated to the mesh T j will be marked with the sign j as follows:

• T

[j]

e is the geometric mapping used to generate K

[j]
e , i.e. K

[j] e = T

[j] e ( K), • {a

[j] 1 , . . . , a [j] N [j] dof
} is the set of the Lagrange nodes of the mesh T j , called level-[j] nodes,

• B j = {ϕ [j] 1 , . . . , ϕ [j] N [j] dof
} is the set of basis functions of the mesh T j , called level-[j] basis functions,

• X j = {v ∈ C 0 (Ω); ∀e ∈ 1, N [j] e , v |K [j] e ∈ P [j]
e } is the H 1 (Ω)-conformal approximation space associated to T j .

Owing to Proposition 1.10, the following result holds: B j is a basis of the H 1 (Ω)-conformal approximation space X j . Definition 1.15. We define the set T j+1 as follows:

T j+1 = K [j+1] ef = T [j] e ( K [1]
f ); 1 e N [j] e , 1 f N [1] e .

Proposition 1.16. The set T j+1 is a geometrically conformal mesh of Ω generated by the reference element K.

Sketch of the proof. For the sake of simplicity, we only give here a sketch of the proof. It is straightforward to see that T j+1 is a mesh of Ω. Hence, we only have to prove that it is geometrically conformal. Let K

[j+1] ef

and

K [j+1]
e ′ f ′ be two cells which have a non empty (d -1)-dimensional intersection, say

F = K [j+1] ef ∩ K [j+1]
e ′ f ′ . The proof is based on different arguments depending on whether e = e ′ or e = e ′ . In the first case, we use the geometric conformity of T whereas in the second case the result is deduced from the geometrical conformity of T j and compatibility requirements 1.12.

The last point of the compatibility requirements 1.12 obviously ensures that a level-[j] node is also a level-

[j + 1] node. That is ∀k ∈ 1, N [j] dof , ∃ℓ ∈ 1, N [j+1] dof , a [j] k = a [j+1] ℓ . (1)
We can now prove that our construction leads to embedded approximation spaces:

Proposition 1.17. It holds that X j ⊂ X j+1 . Proof. Let v ∈ X j , e ∈ {1, . . . , N [j] dof } and f ∈ {1, . . . , N [1] e }. By definition, v |K [j] e ∈ P [j]
e . This is equivalent to v•T

[j]

e ∈ P . Using the compatibility requirements 1.12, we have v•T

[j] e | b K [1] f ∈ P [1] f . Hence, we get v•T [j] e • T f ∈ P , which exactly means that v |K [j+1] ef ∈ P [j+1] ef
, and the claim is proved.

For the following result, we have to introduce relevant indexation of level-[j] and level-[j +1] nodes. Note that level-[j] and level-[j + 1] nodes belonging to K

[j]
e are, by definition, image of nodes of Σ and Σ [1] respectively under the mapping T

[j]
e . Hence, we denote by: • I [j] (e, k) the index of the level-[j] node which is the image of a k under the mapping T

[j] e . That is ∀(e, k) ∈ 1, N [j] e × 1, N , a [j] I [j] (e,k) = T [j]
e ( a k ).

• I [j,1] (e, ℓ) the index of the level-[j + 1] node which is the image of a [START_REF] Bank | The hierarchical basis multigrid method[END_REF] ℓ under the mapping T

[j]
e . That is

∀(e, ℓ) ∈ 1, N [j]
e × 1, N [1] , a

[j+1]

I [j,1] (e,ℓ) = T [j] e ( a [1] ℓ ).
Proposition 1.18 (Refinement equation). The following relationship holds:

∀i ∈ 1, N [j] dof , ϕ [j] i = N [j+1] dof t=1 β [j] it ϕ [j+1] t (RE)
where the coefficients β

[j]
it are given by: 1] (e, ℓ), 0 otherwise. Remark 1.19. Note that the coefficients β kℓ only depend on the refinement pattern. Hence, these coefficients can be computed beforehand. They are in small number and therefore, requiring a low memory cost, they can be stored. Thus, the above refinement equations can be deduced without computation of any coefficient. In practice, coefficients β

∀(i, t) ∈ 1, N [j] dof × 1, N [j+1] dof , β [j] it = β kℓ if ∃(e, k, ℓ) ∈ 1, N [j] e × 1, N × 1, N [1] s.t. i = I [j] (e, k) and t = I [j,
[j]
it are obtained thanks to a loop on level-[j] cells included in supp[ϕ

[j] i ] by setting: for all e ∈ 1, N [j] e such that K [j] e ⊂ supp[ϕ [j] i ], for all (k, ℓ) ∈ 1, N × 1, N [1] , β [j] I [j] (e,k)I [j,1] (e,ℓ) = β kℓ ,
and other coefficients are zero. Remark that, such a loop may lead to consider several times the same pair of indices (I [j] (e, k), I [j,1] (e, ℓ)) for distinct e, k, ℓ. Proposition 1.18 ensures that corresponding coefficients β kℓ are the same.

Proof of Proposition 1.18. Let i ∈ 1, N [j]
dof . The basis function ϕ

[j] i belongs to X j . Since X j ⊂ X j+1 and B j+1 is a basis of X j+1 , the existence of the coefficients β

[j] it is straightforward. Let (i, t) ∈ 1, N [j] dof × 1, N [j+1] dof
. We have:

ϕ [j] i = N [j+1] dof t=1 β [j] it ϕ [j+1] t . (2)
• Case 1 : there exists (e, k, ℓ) ∈ 1, N

[j] e × 1, N × 1, N [1] such that i = I [j] (e, k) and t = I [j,1] (e, ℓ). The restriction of (2) to K

[j] e yields:

ϕ k • T [j] e -1 = b N [1] l=1 β [j] iI [j,1] (e,l) ϕ [j+1] I [j,1] (e,l) .
Owing to Proposition 1.14, we have: b N [1] l=1

β kl ϕ [1] l • T [j] e -1 = b N [1] l=1 β [j] iI [j,1] (e,l) ϕ [1] l • T [j] e -1 .
Evaluating this equality at T

[j] e ( a [1] ℓ ), it follows that β [j] iI [j,1] (e,ℓ) = β kℓ . This is exactly β [j] it = β kℓ . • Case 2: ∀(e, k, ℓ) ∈ 1, N [j] e × 1, N × 1, N [1] , i = I [j] (e, k) or t = I [j,1] (e, ℓ). Let e ∈ 1, N [j] e
and ℓ ∈ 1, N [1] such that t = I [j,1] (e, ℓ). We necessarily have, by assumption, ∀k ∈ 1, N , i = I [j] (e, k). Hence, we get:

0 = N [j+1] dof s=1 β [j] is ϕ [j+1] s |K [j] e = b N [1] v=1 β [j] iI [j,1] (e,v) ϕ [1] v • (T [j] e ) -1 .
Evaluating this equality at T [j] ( a

[1] ℓ ) yields β [j] iI [j,1] (e,ℓ) = 0. This is exactly, β [j] it = 0.
The refinement equation (RE) introduces a relationship between a level-[j] basis function and some level-[j+1] basis functions which are called its children.

Definition 1.20 (Parent-child relationship for basis functions). In the case where β

[j] it = 0, we say that:

• the level-[j] basis function ϕ [j] i is a parent of the level-[j + 1] basis function ϕ [j+1] t , • the level-[j + 1] basis function ϕ [j+1] t is a child of the level-[j] basis function ϕ [j] i .
For this reason, the refinement equation (RE) is also called the parent-child relationship. Along the same lines, we can define a parent-child relationship for cells. • For all f ∈ {1, . . . , N [START_REF] Bank | The hierarchical basis multigrid method[END_REF] e }, we say that the level-

[j + 1] cell K [j+1] ef is a child cell of the level-[j] cell K [j] e .
• Conversely, we say that the level-

[j] cell K [j] e is the parent cell of each level-[j + 1] cell K [j+1] ef , for f ∈ {1, . . . , N [1]
e }. Remark 1.22. A cell has at most one parent cell whereas a basis function may have several parents. Nevertheless, let us identify some basis functions which have only one parent.

Proposition 1.23 (Private Child). Let (k, ℓ) ∈ 1, N [j] dof × 1, N [j+1] dof such that a [j] k = a [j+1] ℓ Then, ϕ [j] k is the unique parent of ϕ [j+1] ℓ . Proof. For 1 i N [j]
dof , the parent-child relationship yields:

ϕ [j] i (a [j] k ) = N [j+1] dof t=1 β [j] it ϕ [j+1] t (a [j+1] ℓ ).
This is exactly:

δ ik = β [j]
iℓ . Thus, the basis function ϕ

[j+1] ℓ has a unique parent which is ϕ [j] k .
Summary 1.24. Let ( K, Σ, P , T ) be a refinement pattern. Let T 0 be a geometrically conformal mesh of Ω generated using the reference element ( K, Σ, P ). By applying uniformly the refinement pattern as described in this section, we are able to construct:

• a hierarchy of nested meshes, T 0 , T 1 , . . . ,

T J (Figures 7, 8, 9), • a hierarchy of nested H 1 (Ω)-conformal finite element approximation spaces X 0 ⊂ X 1 ⊂ • • • ⊂ X J ,
• basis function sets B 0 , B 1 , . . . , B J , spanning the above approximation spaces such that two consecutive sets are linked by refinement equations,

Mesh Basis Function Set Approximation Space

Level 0

T 0 B 0 = {ϕ [0] k ; k = 1, . . . , N [0] dof } X 0 = span B 0 Level 1 T 1 B 1 = {ϕ [1] k ; k = 1, . . . , N [1] dof } X 1 = span B 1 . . . . . . . . . . . . . . . Level J T J B J = {ϕ [J] k ; k = 1, . . . , N [J] dof } X J = span B J Table 1. Conceptual hierarchy of nested conformal FE spaces. level 0 level 1 level 2 level 3 Figure 7. Square-Q 1 . Nested meshes T j . level 0 level 1 level 2 level 3 Figure 8. Tri/Quadr-angle-P 1 /Q 1 . Nested meshes T j .
Table 1 gives a summary of notations used in the sequel.

Note that the hierarchy of nested meshes is never explicitly created. This conceptual structure is introduced in order to explain the refinement method but in practice, the refinement pattern can be applied only locally where basis functions have to be effectively refined.

Multilevel basis and multilevel approximation spaces

We assume that a structure as presented in Summary 1.24 is given and we use the same notation (Table 1).

The aim of this subsection is to explain how we can select some basis functions in

  ∀(j, j ′ ) ∈ 0, J 2 , ∀(k, k ′ ) ∈ 1, N [j] dof × 1, N [j ′ ]
dof such that a

[j] k = a [j ′ ] k ′ , (ϕ [j] k ∈ B, ϕ [j ′ ] k ′ ∈ B) =⇒ (j = j ′ , k = k ′ ),   (P LI )
then B is linearly independent.

Remark 1.26. Note that (a

[j] k = a [j ′ ] k ′ and j = j ′ ) =⇒ k = k ′ . Proof. The property (P LI ) implies that   ∀(j, j ′ ) ∈ 0, J 2 , ∀(k, k ′ ) ∈ 1, N [j] dof × 1, N [j ′ ] dof such that j ′ > j, (ϕ [j]
k ∈ B, ϕ

[j ′ ] k ′ ∈ B) =⇒ ϕ [j ′ ] k ′ (a [j] k ) = 0,   (3)
because, owing to (1), a

[j]
k is also a level-[j ′ ] node and, by (P LI ), this node is certainly different from a

[j ′ ] k ′
(otherwise j ′ = j ). Consider a linear combination of basis functions ϕ belonging to B such that

ϕ∈B λ ϕ ϕ = 0, (4)
and assume that E = {ϕ ∈ B; λ ϕ = 0} is not empty. We can then define

j m = min{j ∈ 0, J ; ∃k ∈ 1, N [j] dof such that ϕ [j] k ∈ E},
and select a k m ∈ 1, N [jm] dof such that ϕ [jm] km ∈ E. Let j ∈ 0, J and k ∈ 1, N [j] dof such that ϕ [j]
k ∈ E and (j, k) = (j m , k m ). Owing to (3), we have ϕ Owing to Proposition 1.25, we give the following definition of a multilevel basis. Definition 1.27 (Multilevel basis and multilevel approximation space). We call a multilevel level basis a subset B of J j=0 B j satisfying the property (P LI ). By Proposition 1.25, this set is actually independent. A space spanned by a multilevel basis is called multilevel approximation space.

Remark 1.28. Let V = span B be a multilevel approximation space and u ∈ V. The coordinates of the expansion of u in the multilevel basis B are not necessarily the values of u at the nodes associated to the corresponding basis function, since two basis functions of different levels may have overlapping supports.

Adaptation procedure and multigrid preconditioners

The adaptation consists in adding or removing some basis functions of a given multilevel basis B ⋆ in order to produce a new multilevel basis B whose spatial resolution will be better suited to the problem. The main points are to ensure that the Un/Refinement algorithm will actually produce a linearly independent family of basis functions, and that no information is loss during the refinement process. Section 2.1 is devoted to the proofs of such properties. Then, in Section 2.2, we show how to also use the multilevel structure obtained by the adaptation algorithm in order to build a multigrid preconditioning algorithm which let us solve a multilevel linear system with a moderate computational cost. k is said to be B-refined iff:

∃j ′ ∈ j + 1, J , ∃k ′ ∈ 1, N [j ′ ] dof , such that ϕ [j ′ ] k ′ ∈ B and a [j] k = a [j ′ ] k ′ .
Moreover, if the above condition holds for j ′ = j + 1 then we say that the basis function ϕ

[j] k is B-refined only once.
Remark 2.2. Owing to Property (P LI ), notice that:

• the indices j ′ and k ′ are necessarily unique.

• a B-refined basis function does not belong to B.

Remark also that, by Proposition 1.23, if ϕ

[j]
k is only once B-refined then ϕ

[j]
k is the unique parent of ϕ

[j ′ ]
k ′ . Let us give the following lemma which will be useful in the following proofs.

Lemma 2.3. Let j ∈ 0, J and j ′ ∈ 0, J -1 such that j ′ j. Let (k, ℓ, ℓ ′ ) ∈ 1, N [j] dof × 1, N [j+1] dof × 1, N [j ′ ] dof . If a [j+1] ℓ = a [j ′ ] ℓ ′ and ϕ [j] k is a parent of ϕ [j+1] ℓ then the node a [j]
k is necessarily at the same position that the nodes a

[j ′ ] ℓ ′ and a [j+1] ℓ .
Proof. Since j ′ j, owing to (1) and a straightforward recurrence, there exists t ∈ 1, N

[j]
dof such that a

[j] t = a [j+1] ℓ
. We can apply Proposition 1.23 which proves that ϕ

[j] t is the unique parent of ϕ [j+1] ℓ . Hence, ϕ [j] k = ϕ [j] t and then, k = t.
We can now describe the refinement and unrefinement procedure. Algorithm 2.4 (Quasi-Hierarchical Un/Refinement). Let B ⋆ be a multilevel basis.

• Refinement: Let ϕ This can be written in a compact way as follows:

B = B ⋆ \{ϕ [j] k } ∪ {children of ϕ [j] k not B ⋆ -refined}. • Unrefinement: Let ϕ [j]
k be an only once B ⋆ -refined basis function without B ⋆ -refined children. Unrefining the given basis function ϕ k which have no other B ⋆ -refined parent . This can be written in a compact way as follows:

B = B ⋆ \{children of ϕ [j] k without other B ⋆ -refined parent} ∪ {ϕ [j] k }.
Proof. As claimed in the above algorithm, we have to prove that refinement and unrefinement procedures produce actually a multilevel basis. That is to say that B satisfies the property (P LI ). Indeed, let ϕ

[j] k , ϕ [j ′ ] k ′ two basis functions belonging to B such that a [j] k = a [j ′ ]
k ′ . We have to show that j = j ′ .

• Refinement: Assume that B is obtained from B ⋆ by the refinement of a basis function belonging to B ⋆ , say ϕ

[j0] k0 ∈ B ⋆ , so that we have

B = B ⋆ \{ϕ [j0] k0 } ∪ {children of ϕ [j0] k0 not B ⋆ -refined}. -If ϕ [j]
k and ϕ

[j ′ ]
k ′ belong both to B ⋆ then, since B ⋆ satisfies the property (P LI ), we readily find j = j ′ . -Otherwise, assume for instance that ϕ

[j] k does not belong to B ⋆ . By definition of B, ϕ [j] k is then a child of ϕ [j0] k0 which is not B ⋆ -refined, say ϕ [j0+1] ℓ0 , i.e. k = ℓ 0 , j = j 0 + 1. * If ϕ [j ′ ] k ′ ∈ B ⋆ then, since ϕ [j0+1] ℓ0 is not B ⋆ -refined, we have j 0 +1 j ′ . Assume that j 0 +1 > j ′ . Lemma 2.3 yields a [j0] k0 = a [j ′ ]
k ′ . Since B ⋆ satisfies the property (P LI ), we obtain j 0 = j ′ and

k 0 = k ′ , but ϕ [j0] k0 ∈ B and ϕ [j ′ ] k ′ ∈ B.
This is a contradiction and we get j = j 0 + 1 = j ′ . * Otherwise, ϕ

[j ′ ] k ′ is a child of ϕ [j0] k0 , say ϕ [j0+1] ℓ0 ′ , i.e. k ′ = ℓ 0
′ and j = j 0 + 1. Hence, we have j = j ′ . • Unrefinement: Assume that B is obtained from B ⋆ by the unrefinement of a only once B ⋆ -refined basis function, say ϕ

[j0]
k0 , so that we have

B = B ⋆ \{children of ϕ [j0] k0 without other B ⋆ -refined parent} ∪ {ϕ [j0] k0 }. -If ϕ [j]
k and ϕ

[j ′ ]
k ′ belong to B ⋆ then, since B ⋆ satisfies the property (P LI ), we get j = j ′ . -Otherwise, assume for instance that ϕ

[j]
k does not belong to B ⋆ . By definition of B, we have ϕ

[j] k = ϕ [j0]
k0 , i.e. k = k 0 , j = j 0 . Arguing by contradiction, we assume that j ′ = j. We have j ′ = j 0 and so, ϕ

[j ′ ] k ′ ∈ B ⋆ . However, ϕ [j0]
k0 is an only once B ⋆ -refined basis function, so there exists

ℓ 0 ∈ 1, N [j0] dof such that ϕ [j0+1] ℓ0 ∈ B ⋆ and a [j0+1] ℓ0 = a [j0] k0 . Owing to Proposition 1.23, a [j0+1] ℓ0 = a [j0] k0 implies that ϕ [j0] k0 is the unique parent of ϕ [j0+1] ℓ0 . Therefore, since ϕ [j0+1] ℓ0 is a child of ϕ [j0]
k0 without other B ⋆ -refined parent, we get ϕ

[j0+1] ℓ0 ∈ B. However, since B ⋆ satisfies the property (P LI ) and since

a [j0+1] ℓ0 = a [j ′ ]
k ′ , we have j ′ = j 0 + 1 and k ′ = ℓ 0 . Finally, we get ϕ

[j ′ ] k ′ ∈ B.
This is a contradiction and we obtain j = j ′ . Remark 2.5. This algorithm is consistent with Definition 2.1. Indeed, owing to Proposition 1.23,

• if B is obtained from B ⋆ by the refinement of the basis function ϕ

[j] k , then ϕ [j]
k is a B-refined basis function in sense of Definition 2.1 ;

• if B is obtained from B ⋆ by the unrefinement of the only once B ⋆ -refined basis function ϕ

[j] k , then ϕ [j] k
is no longer a B-refined basis function in the sense of Definition 2.1.

ϕ [1] 6 ϕ [0] 2 B ⋆ of ϕ [0] 1 ϕ [0] 2 ϕ [1] 3 B ϕ [1] 4 ϕ [1] 5 ϕ [1] 2 ϕ [1] 1 ϕ [1] 3 ϕ [0] 1 Level-[0] ϕ [1] 6 Level-[1] Refinement ϕ [0] 2 = ϕ [1] 3 + 1 2 ϕ [1] 2 + 1 2 ϕ [1]
6 + 1 4 ϕ

[1] 5

Figure 10. Refinement do not preserve the linear independence of multilevel set not satifying (P LI ).

Remark 2.6. Refinement and unrefinement procedure described in Algorithm 2.4 do not preserve in general the linear independence of multilevel basis function set

J j=0
B j (with B j ⊂ B j ) which do not satisfy the property (P LI ). An example is given in Figure 10. The family B ⋆ represented on the left hand-side is linearly independent (but do no satisfy (P LI )) whereas the family B, on the right hand-side, obtained from B ⋆ by refinement of ϕ

[0]
1 , is not a linear independent family.

Conservation of information

A desirable property of a refinement procedure is that it does not involve a loss of information. It means that, if B is obtained from B ⋆ by refinement of a basis function, then span B ⋆ ⊂ span B, i.e. the refined basis B allows for any function in the original basis B ⋆ to be reproduced exactly. However, the refinement procedure described in Algorithm 2.4 is not lossless. An example is given in Figure 11. Nevertheless, we can prove the following results. Proposition 2.7. Let B be a multilevel basis satisfying the following property:

A child of a B-refined basis function either belongs to B or is itself B-refined.

(P LO )

Then, all B-refined basis functions belong to span B.

Proof. By a recurrence on the level j of basis functions, we prove the following statement (H j ) for all j ∈ 0, J :

All level-[j] B-refined basis functions belong to span B. (H j )

A level-[J] basis function cannot be B-refined. Hence, the statement (H J ) is straightforward. Let j ∈ 0, J -1 . Assume that the statement (H j+1 ) holds and let ϕ

[j]
k be a level-[j] B-refined basis function. Owing to Proposition 1.18, we have ϕ

[j] k = ℓ|ϕ [j+1] ℓ is a child of ϕ [j] k β [j] kℓ ϕ [j+1] ℓ . Level-[2] Level-[0] ϕ [1] 4 ϕ [1] 5 ϕ [1] 1 ϕ [2] 3 ϕ [0] 1 Refinement ϕ [0] 1 of Level-[1] ϕ [2] 3 ϕ [0] 1 ∈ span {ϕ [1] 1 , ϕ [1] 4 , ϕ [1] 5 , ϕ [2]
3 }

Figure 11. Refinement of multilevel basis is not lossless.

Moreover, property (P LO ) implies that all basis functions ϕ

[j+1] ℓ involving in the above sum are either in B or B-refined. In the last case, the recurrence assumption (H j+1 ) yields ϕ Proof. Assume that B is obtained from B ⋆ by the refinement of the basis function ϕ

[j0] k0 ∈ B ⋆ . Let ϕ [j] k ∈ B ⋆ . • If ϕ [j] k = ϕ [j0] k0 then ϕ [j] k ∈ B, • else ϕ [j] k = ϕ [j0]
k0 which is B-refined. Proposition 2.7 ensure that ϕ

[j] k ∈ span B. Then span B ⋆ ⊂ span B.
Moreover, the refinement and unrefinement procedures preserve the property (P LO ). Let us begin with the following lemma. Lemma 2.9. Let B ⋆ and B be two multilevel basis. Let ϕ

[j ⋆ ]
k ⋆ be a B ⋆ -refined basis function and ϕ

[j] k be a B-refined basis function.

1) If B is obtained from B ⋆ by the refinement of a basis function ϕ

[j0] k0 ∈ B ⋆ then (i) ϕ [j ⋆ ] k ⋆ is also B-refined, (ii) ϕ [j] k is either B ⋆ -refined or equal to ϕ [j0] k0 .
2) If B is obtained from B ⋆ by the unrefinement of an only once B ⋆ -refined basis function ϕ

[j0] k0 then (i) ϕ [j ⋆ ] k ⋆ is either B-refined or is equal to ϕ [j0] k0 , (ii) ϕ [j] k is also B ⋆ -refined. (iii) Moreover, if we assume that ϕ [j ⋆ ] k ⋆ is only once B ⋆ -refined then ϕ [j ⋆ ] k ⋆ is either also only once B-refined or is equal to ϕ [j0] k0 . Proof. Since ϕ [j ⋆ ] k ⋆ is B ⋆ -refined, there exists j ⋆′ > j ⋆ and k ⋆′ ∈ 1, N [j ⋆ ′ ] dof such that ϕ [j ⋆ ′ ] k ⋆ ′ ∈ B ⋆ and a [j ⋆ ] k ⋆ = a [j ⋆ ′ ] k ⋆ ′ . Along the same lines, since ϕ [j] k is B-refined, there exists j ′ > j and k ′ ∈ 1, N [j ′ ] dof such that ϕ [j ′ ] k ′ ∈ B and a [j] k = a [j ′ ] k ′ . 1) Refinement of ϕ [j0] k0 . (i) If ϕ [j ⋆ ′ ] k ⋆ ′ ∈ B then ϕ [j ⋆ ] k ⋆ is B-refined; else ϕ [j ⋆ ′ ] k ⋆ ′ = ϕ [j0] k0 and since ϕ [j0] k0 is B-refined, a fortiori ϕ [j ⋆ ] k ⋆ is B-refined. (ii) If ϕ [j ′ ] k ′ ∈ B ⋆ then ϕ [j]
k is B ⋆ -refined and the claim is proved. Otherwise, ϕ

[j ′ ] k ′ is a child of ϕ [j0] k0 , say ϕ [j0+1] ℓ0 , i.e. j ′ = j 0 + 1, k ′ = ℓ 0 . We have a [j] k = a [j0+1] ℓ0
and j 0 + 1 > j, then owing to Lemma 2.3,

a [j] k = a [j0] k0 . If j 0 = j then ϕ [j] k = ϕ [j0] k0 ; else j 0 > j, since ϕ [j0] k0 ∈ B ⋆ , ϕ [j] k is B ⋆ -refined. 2) Unrefinement of ϕ [j0] k0 . (i) If ϕ [j ⋆ ′ ] k ⋆ ′ ∈ B then ϕ [j ⋆ ] k ⋆ is B-refined and the claim is proved. Otherwise, ϕ [j ⋆ ′ ] k ⋆ ′ is a child of ϕ [j0] k0 , say ϕ [j0+1] ℓ0
, i.e. j ⋆′ = j 0 + 1, k ⋆′ = ℓ 0 . We have a

[j ⋆ ] k ⋆ = a [j0+1] ℓ0 and j 0 + 1 > j ⋆ , then owing to Lemma 2.3, a [j ⋆ ] k ⋆ = a [j0] k0 . If j 0 = j ⋆ then ϕ [j ⋆ ] k ⋆ = ϕ [j0] k0 ; else j 0 > j ⋆ , since ϕ [j0] k0 ∈ B, ϕ [j ⋆ ] k ⋆ is B-refined. (ii) If ϕ [j ′ ] k ′ ∈ B ⋆ then ϕ [j] k is B ⋆ -refined; else ϕ [j ′ ] k ′ = ϕ [j0] k0 and since ϕ [j0] k0 is B ⋆ -refined, a fortiori ϕ [j] k is k ⋆ = ϕ [j0] k0 .
Proposition 2.10. The refinement and unrefinement procedures of multilevel basis (in the sense of Algorithm 2.4) preserves the property (P LO ).

Proof. Assume that the multilevel basis B ⋆ satisfies the property (P LO ). Let ϕ k0 ∈ B ⋆ . By Lemma 2.9 property 1) (ii), ϕ

k is either B ⋆ -refined or equal to ϕ

[j0] k0 . Consider the two cases:

-

If ϕ [j]
k is B ⋆ -refined. Since B ⋆ satisfies the property (P LO ), only the two following cases are possible: * ϕ

[j+1] ℓ is B ⋆ -refined. Owing to Lemma 2.9 property 1)(i), this implies that ϕ

[j+1] ℓ is B-refined. * ϕ [j+1] ℓ ∈ B ⋆ . And then, ϕ [j+1] ℓ ∈ B or ϕ [j+1] ℓ = ϕ [j0] k0 which is B-refined. -If ϕ [j] k = ϕ [j0]
k0 , all its children are either in B or are B ⋆ -refined. Owing to Lemma 2.9 property 1)(i), they are either in B or B-refined.

• Unrefinement: Assume that B is obtained from B ⋆ by the unrefinement of an only once refined basis function ϕ

[j0] k0 of B ⋆ . By Lemma 2.9 property 2) (ii), ϕ

[j]

k is B ⋆ -refined. Since B ⋆ satisfies the property (P LO ), only the two following cases are possible:

ϕ

[j+1] ℓ is B ⋆ -refined. Owing to Lemma 2.9 property 2)(i), this implies that ϕ

[j+1] ℓ is B-refined or ϕ [j+1] ℓ = ϕ [j0] k0 ∈ B. -ϕ [j+1] ℓ ∈ B ⋆ . And then, we have ϕ [j+1] ℓ ∈ B or ϕ [j+1] ℓ is a child of ϕ [j0]
k0 with no other B ⋆ -refined parent. The last case is impossible. Indeed, ϕ

[j]

k is B ⋆ -refined and then, by unicity, we would have ϕ

[j] k = ϕ [j0] k0 . However, ϕ [j0]
k0 is not a B-refined function. This is a contradiction and we get ϕ

[j+1] ℓ ∈ B.
Remark 2.11. Note that the properties (P LI ) or (P LO ) are not so restrictive since they are preserved by Un/Refinement procedures and since it is straightforward to see that they are satisfied by the coarse basis B 0 which is used, in practice, for starting the adaptation algorithm.

Combining ( 9) and ( 12), we get B = B ⋆ \{children of ϕ and ψ which have no B ⋆ -refined parents except possibly ϕ and ψ} ∪ {ϕ, ψ}.

This expression shows that B do not depend of the order in which the basis functions ϕ and ψ were refined.

With this definition at hand, we can give the refinement algorithm.

Algorithm 2.13 (Adaptation procedure). Let B ⋆ a multilevel basis. Assume that, thanks to a refinement criterion, we are given the set E ⋆ ⊂ B ⋆ of basis functions to refine and the set F ⋆ of only once B ⋆ -refined basis functions (without B ⋆ -refined children) to unrefine. The adaptation procedure consists in the two following steps:

1) Refine the set E ⋆ , thus producing a new multilevel basis B.

2) Unrefine the set of basis functions of F ⋆ which are still only once B-refined basis functions without B-refined children.

One-level-difference rule

In practice, the following criterion is used to ensure the common rule of "one-level-difference refinement" [START_REF] Krysl | Natural hierarchical refinement for finite element methods[END_REF].

Criterion 2.14 (One-level-difference refinement rule). Let B ⋆ a multilevel basis. A basis function of B ⋆ may be refined only if all its parents are B ⋆ -refined.

Typically, this criterion is used to avoid an important difference of refinement level between "neighboring" basis functions. An example of refinement sequence forbidden by this criterion is given in Figure 12.

Level-[0]

Level- [START_REF] Bank | Some refinement algorithms and data structures for regular local mesh refinement[END_REF] Figure 12. Refinement sequence forbidden by the One-level-difference rule. This criterion can be taken into account by adding a recursive step in the refinement procedure: to refine a basis function, we first refine all its parents which are not B ⋆ -refined and then we refine the basis function itself. This is illustrated, for instance, by the spreading of the level- [START_REF] Bank | The hierarchical basis multigrid method[END_REF] refined area between the first two pictures in Figure 17.

Multigrid preconditioner

In this section, we assume that we have the following genuine variational problem to solve: Find u ∈ V h such that ∀v ∈ V h , a(u, v) = l(v), where V h = span B is a multilevel approximation space, a : H 1 (Ω) × H 1 (Ω) → R is a bilinear continuous and coercive form and b : H 1 (Ω) → R is a linear continuous form. We denote by A J the stiffness matrix associated to this problem, that is A J = [a(ϕ, ψ)] ϕ,ψ∈B . The multilevel space V h is built to achieve a given level of accuracy without increasing too much the number of degrees of freedom in the discrete problem. In fact, we are going to describe how to naturally take advantage of the multilevel structure of the approximation space in order to build a sequence of nested multilevel grids finally leading to a multigrid preconditioner.

Coarsening

From a given "fine" multilevel basis B F , the following algorithm is used to construct a "coarser" multilevel basis B C . Proof. Note first that any step in Algorithm 2.15 consists in an unrefinement of a level-[j M -1] basis function, this implies that added basis functions are certainly on level-[j M -1] and that removed basis functions are certainly on level-[j M ]. Hence, the set of added and removed basis functions are disjoint. In Algorithm 2.15, a basis function which is removed (or added) by an unrefinement procedure can not be added (or respectively removed) by an other unrefinement procedure. Furthermore, remark that by definition of j M , there is no B Frefined function on level-[j M ]. Hence, since they have no B F -refined children, all B F -refined basis functions can actually be unrefined. Unrefinement of a basis functions involved that it is added and then the set of added basis functions in Algorithm 2.15 is exactly the set of B F -refined basis functions. It remains to show that all level-[j M ] basis functions of B F are removed. Arguing by contradiction, assume that a level-[j M ] basis function belongs to coarsen(B F ). Since unrefinement preserves the property (P HI ) (Proposition 2.17), this basis function has at least one coarsen(B F )-refined parent. Owing to Lemma 2.9 property 2) (ii), this is also a B F -refined basis function (of level-[j M -1]). This is a contradiction and the claim is proved.

The property (P HI ) ensures the desirable fact that all level-[j M ] basis functions of B F are removed. This additional property is not so restrictive since it is preserved by the refinement and unrefinement procedures.

Proposition 2.17. The refinement and unrefinement procedures described in Algorithm 2.4 preserve the property (P HI ).

Proof. Assume that B ⋆ is a multilevel basis satisfying the property (P HI ) and let j ∈ 1, J .

• Consider first a basis function ϕ [j]

k which belongs to B.

-Refinement: Assume that the multilevel basis B is obtained from B ⋆ by the refinement of a basis function ϕ

[j0] k0 ∈ B ⋆ . * Case 1: ϕ [j]
k ∈ B ⋆ . Since B ⋆ satisfies the property (P HI ), ϕ

[j]

k has at least one B ⋆ -refined parent. Owing Lemma 2.9 property 1), this parent is B-refined. * Case 2: ϕ

[j] k is a child of ϕ [j0] k0 and ϕ [j0] k0 is B-refined.
-Unrefinement: Assume that B is obtained from B ⋆ by the unrefinement of an only once refined basis function ϕ

[j0] k0 of B ⋆ . * Case 1: ϕ [j] k ∈ B ⋆ \{children of ϕ [j]
k without other B ⋆ -refined parent}. Since B ⋆ satisfies the property (P HI ), ϕ

[j]

k has at least one B ⋆ -refined parent. Owing Lemma 2.9 property 2), either this parent is B-refined or this is ϕ

[j0] k0 .
In the first case, the proof is finished, in the second case ϕ

[j] k is a child of ϕ [j0] k0 which belongs to B ⋆ \{children of ϕ [j] k without other B ⋆ -refined parent}. Therefore, ϕ [j]
k has an other B ⋆refined parent. This parent is then B-refined because it is not ϕ

[j0] k0 . * Case 2: ϕ [j] k = ϕ [j0] k0 , ϕ [j0]
k0 is B ⋆ -refined. Since B ⋆ satisfies the property (P HI ), ϕ

[j0]
k0 has at least one B ⋆ refined parent. Owing Lemma 2.9 property 2), this parent is B-refined because it can not be ϕ

[j0] k0 .

• Consider now a basis function ϕ

[j]
k which is B-refined.

-Refinement: Assume that the multilevel basis B is obtained from B ⋆ by the refinement of a basis function ϕ

[j0] k0 ∈ B ⋆ . By Lemma 2.9 property 1) (ii), ϕ

[j] k is either B ⋆ -refined or equal to ϕ [j0]
k0 ∈ B ⋆ . Since B ⋆ satisfies the property (P HI ), in the two cases, ϕ

[j]
k has at least one B ⋆ -refined parent. Owing to Lemma 2.9 property 1)(i), this parent is also B-refined.

-Unrefinement: Assume that B is obtained from B ⋆ by the unrefinement of an only once refined basis function ϕ

[j0] k0 of B ⋆ without B ⋆ -refined children. By Lemma 2.9 property 2) (ii), ϕ [j]
k is B ⋆ -refined. Since B ⋆ satisfies the property (P HI ), ϕ

[j]

k has at least one B ⋆ -refined parent. Owing to Lemma 2.9 property 2) (i), either this parent is B-refined or is equal to ϕ The last main property of the coarsening procedure, is that it produced nested vector spaces. • Additive version [START_REF] Bramble | Parallel multilevel preconditioners[END_REF]:

P a = J k=0 I k S k I t k
where S k , k ∈ 0, J , is the Jacobi smoother. • Multiplicative version [START_REF] Hackbusch | Multigrid Methods and Applications[END_REF]: This correspond to the classical V -cycle. In this section S k is the Gauss-Seidel smoother. We define recursively, for all k ∈ 0, J , the linear operator

M G k : R # B k → R # B k . We first set M G 0 = A -1 0 and for all k ∈ 1, J , we define M G k (f k ), f k ∈ R # B k by the following steps: (0) v k 0 Initialisation (1) v k v k + S k (f k -A k v k ), Pre-smoothing step (2) v k v k + I k-1 M G k-1 (I t k-1 (f k -A k v k )), Coarse Grid Correction (3) v k v k + S k (f k -A k v k ), Post-smoothing step Set M G k (f k ) = v k .
The multiplicative preconditioner is then

P m = M G J .

Validation on a stationnary model problem

We first validate the local refinement method and the multigrid preconditioner on the following stationnary model problem.

The practical implementation has been performed using the software object-oriented component library PELICANS [START_REF] Pelicans | Collaborative Development environment[END_REF], developed at the "Institut de Radioprotection et de Sûreté Nucléaire (IRSN)" and distributed under the CeCILL-C license agreement (an adaptation of LGPL to the French law).

Continuous model problem

Let Ω = [0, 1] d . Consider the Laplace problem with homogeneous Dirichlet boundary conditions:

-∆u = f in Ω, u = 0 on ∂Ω. ( 13
)
The source term f is chosen so that the exact solution u is:

∀x ∈ R d , u(x) = H ε R -|x -x C | , where R ∈ R, ε ∈ R * + , x C ∈ R d are parameters and H ε : R -→ R is defined by ∀x ∈ R, H ε (x) =          0 if x < -ε, 1 2 1 + x ε + 1 π sin π x ε if |x| ε, 1 if x > ε.
The function H ε is represented in Figure 14 and the interpretation of parameters R, ε, x C is explained in Figure 15 which represents the exact solution u when d = 1.

For numerical simulations given in the sequel, we have set:

Ω = [0, 1] d , x C = (

Initial meshes and refinement criterion

In order to show the capability of the local refinement method, we used several types of meshes: square, triangle, general quadrangle, in two dimension and cube in three dimension. Figure 16 shows the initial meshes, i.e. before refinement steps, which have been used for the validation in two dimension. In particular, note that the method allows to combine several kind of geometric elements (triangle P 1 and square Q 1 in the example here) provided that the corresponding reference elements, and refinement patterns are compatible. In three dimension, the initial mesh is a cubic regular mesh obtained by dividing all sides in 15 segments. In this section, the shape of the exact solution is known. This is the reason why a geometric refinement criterion is used. It consists in choosing a priori a number of refinement steps and the location of the nodes associated to the basis functions to refine. Here the following criterion is used:

Criterion 2.21 (Geometric criterion). ϕ [j] k will be refined iff R -ε < a [j] k -x C < R + ε 2.3.

Local refinement and multigrid preconditioner

For each type of initial mesh, we performed six computations increasing the number of refinement steps from one to six and following the geometric refinement criterion 2.21. In Figures 17,18, 19 and 20, for smaller number of refinement steps (one, two and three steps), we represent the refined meshes and the plot of the function u h (1u h ) where u h is the computed approximate solution. Hence, the colored zone in the figures represents the computed "interface", that is the areas where local refinement takes place. Note that meshes in our method are obtained as support of basis functions involved in the formulation of the problem (see Section 3.3.2). In particular, we remark the effect of the "one-level-difference" refinement rule. Indeed, for each type of mesh, the "one level refinement" areas spreads during the second refinement step. This is due to the fact that during this refinement step, few parents of level- [START_REF] Bank | The hierarchical basis multigrid method[END_REF] basis functions are refined to enable the refinement of their childs.

One Refinement step

Two refinement steps Three refinement steps Table 2 shows the convergence rates in the L 2 -norm, computed with respect to the mesh size of the finest refinement level used in the computation. We obtain a convergence rate in the L 2 -norm equal to 2 for first order finite element (P 1 and Q 1 ) and equal to 3 for second order finite element (P 2 and Q 2 ) as expected.

Dimension

Mesh-Finite Element Convergence rate

2D Square-Q 1 1.99 Quadrangle-Q 1 1.99 Tri/Quadr-angle-P 1 /Q 1 1.99 3D Cube-Q 1 1.99 2D Square-Q 2 2.99 Quadrangle-Q 2 2
.90 Tri/Quadr-angle-P 2 /Q 2 2.97 Table 2. Convergence rates in the L 2 -norm, computed with respect to the mesh size of the finest refinement level used in the computation.

In order to illustrate the capability of multigrid preconditioners, we represent in Tables 3, 4 and 5 the number of necessary iterations of the conjugate gradient solver as a function of the number of unknowns. In Table 3, we show the results for two dimensional first order elements, in Table 4, for two dimensional second order elements and in Table 5 for three dimensional first order elements. Various preconditioners are compared : the classical Incomplete LU factorization (ILU0) and the additive (P a ) and multiplicative (P m ) versions of the multigrid algorithm. Note that in all those computations we say that convergence in the conjugate gradient method is achieved as soon as the relative L ∞ -norm of the residual is less than 10 -10 . We limit the total number of iterations to 600, and we denote by "-" in Tables 3, 4 and 5 the case where the convergence of the conjugate gradient method is not achieved before 600 iterations.

Without preconditionning, the number of iterations rapidly exceeds the limit when increasing the number of unknowns. The use of preconditioners enable to reduce the number of iterations and we observe in particular that the number of necessary iterations using (P m ) and (P a ) is almost independent of the size of the problem and is significantly smaller than using, for example, (ILU0) preconditioner.

Refinement

Application to a ternary Cahn-Hilliard system

We are interested in the numerical computation of the solution of a three-phase Cahn-Hilliard model proposed in [START_REF] Boyer | Study of a three component Cahn-Hilliard flow model[END_REF]. After a brief introduction of this model, we describe the time and space discretization of the problem that we are using. Then, a numerical example is given using local adaptive refinement and multigrid preconditioner.

Introduction

The Cahn-Hilliard approach consists in assuming that the interfaces between phases in the system have a small but positive thickness ε. Each phase i is represented by a smooth function c i called the order parameter Refinement level 1 2 3

Cube-Q (which is taken to be the volumic fraction of the component in the mixture). Since the mixture is supposed to be perfect, the three unknowns c 1 , c 2 and c 3 are linked though the relationship

c 1 + c 2 + c 3 = 1. ( 14
)
The model we consider has been fully derived and studied in [START_REF] Boyer | Study of a three component Cahn-Hilliard flow model[END_REF] as a generalization of the two-phase Cahn-Hilliard model. In the diphasic case, the free energy of the mixture depends on two parameters: the interface width ε and the surface tension σ. It can be written as follows:

F diph σ,ε (c) = Ω 12 σ ε c 2 (1 -c) 2 + 3 4 σε |∇c| 2 dx.
Therefore, in [START_REF] Boyer | Study of a three component Cahn-Hilliard flow model[END_REF], the authors have postulated that the three-phase free energy can be written as follows:

F triph Σ,ε (c 1 , c 2 , c 3 ) = Ω 12 ε F (c 1 , c 2 , c 3 ) + 3 8 εΣ 1 |∇c 1 | 2 + 3 8 εΣ 2 |∇c 2 | 2 + 3 8 εΣ 3 |∇c 3 | 2 dx. ( 15
)
The triple of constant parameters Σ = (Σ 1 , Σ 2 , Σ 3 ) and the bulk energy F have been determined so that the model fits with the prescribed surface tension σ 12 , σ 13 and σ 23 and is "consistent" with the two-component situation (see below). The evolution of the system is then driven by the gradient of the total free energy F triph Σ,ε and the time evolution of c = (c 1 , c 2 , c 3 ) is governed by the following system of equations:

       ∂c i ∂t = ∇ • M 0 Σ i ∇µ i , for i = 1, 2, 3 µ i = f F i (c) - 3 4 εΣ i ∆c i , for i = 1, 2, 3 (16)
where M 0 is a constant (it could also depend on c) and

f F i (c) = 4Σ T ε j =i 1 Σ j (∂ i F (c) -∂ j F (c)) with Σ T defined by 3 Σ T = 1 Σ 1 + 1 Σ 2 + 1 Σ 3 .
This choice of f F i , obtained by the use of a Lagrange multipliers technique, enforces the condition [START_REF] Jobelin | A finite element penalty-projection method for incompressible flows[END_REF] all along the time. Thus, one of the unknowns can be arbitrarily eliminated from the system [START_REF] Krysl | Object-oriented hierarchical mesh refinement with CHARMS[END_REF] which is then equivalent to:

                     ∂c i ∂t = ∇ • M 0 Σ i ∇µ i , for i = 1, 2 c 3 = 1 -c 1 -c 2 µ i = f F i (c) - 3 4 εΣ i ∆c i , for i = 1, 2 µ 3 = - Σ 3 Σ 1 µ 1 + Σ 3 Σ 2 µ 2 . (17)
For simplicity, we will later omit to write the equations giving c 3 and µ 3 as functions of c 1 , c 2 and µ 1 , µ 2 respectively. The consistency (or algebrical consistency) of the model defined by ( 15) and ( 16) with the diphasic systems corresponding to the given surface tensions σ 12 , σ 13 , σ 23 respectively means that the following properties hold:

• When the component i is not present, that is c i ≡ 0, the total free energy F triph Σ,ε (c 1 , c 2 , c 3 ) of the system has to be exactly equal to the total free energy F diph σ jk ,ε (c j ) of the diphasic system corresponding to the two other phases.

• When the component i is not present in the mixture at the initial time, the component i must not appear during the time evolution of the system. Finally, it is shown in [START_REF] Boyer | Study of a three component Cahn-Hilliard flow model[END_REF] that the model defined by ( 15) and ( 16) is algebraically consistent with the diphasic systems of surface tensions σ 12 , σ 13 , σ 23 respectively if and only if we have

∀i ∈ {1, 2, 3}, Σ i = σ ij + σ ik -σ jk , (18)
and there exists a smooth function G such that

∀c ∈ R 3 such that c 1 + c 2 + c 3 = 1, F (c) = σ 12 c 2 1 c 2 2 + σ 13 c 2 1 c 2 3 + σ 23 c 2 2 c 2 3 + c 1 c 2 c 3 (Σ 1 c 1 + Σ 2 c 2 + Σ 3 c 3 ) + c 2 1 c 2 2 c 2 3 G(c).
In the physical literature, the coefficient S i = -Σ i defined by ( 18) is well known and called the spreading coefficient of the phase i at the interface between phases j and k. If S i is positive (that is Σ i < 0), the spreading is said to be total and if S i is negative, it is said to be partial.

Notice that, in the following study, the coefficients Σ i are not assumed to be positive, so that the model presented above lets us cope with some total spreading situations. However, in order to have a well-posed system, we assume that the following condition holds:

Σ 1 Σ 2 + Σ 1 Σ 3 + Σ 2 Σ 3 > 0, (19)
In this paper we will consider Cahn-Hilliard potentials with the following form

F (c) = σ 12 c 2 1 c 2 2 + σ 13 c 2 1 c 2 3 + σ 23 c 2 2 c 2 3 + c 1 c 2 c 3 (Σ 1 c 1 + Σ 2 c 2 + Σ 3 c 3 ) F0(c) + 3Λc 2 1 c 2 2 c 2 3 P (c)
Under these assumptions, there exists a unique weak solution for the problem ( 16) (in 2D and with a slightly modified potential in 3D).

Discretization procedure

In this section, we present the time and space discretizations of problem (17) that we used. We also state two results ensuring existence of the discrete solution and convergence of this solution towards the weak solution of the continuous problem.

Time discretization

Let N ∈ N * . The temporal interval [0, t f ] is uniformly discretized with a fixed time step ∆t = t f N . For n ∈ 0, N , we define t n = n∆t. We use a semi-implicit discretization in time with a special care for nonlinear terms. The scheme is written as follows,

       c n+1 i -c n i ∆t = ∇ • M 0 Σ i ∇µ n+1 i , for i = 1, 2, µ n+1 i = D F i (c n+1 , c n ) - 3 4 εΣ i α∆c n i + (1 -α)∆c n+1 i , for i = 1, 2,
where 0 α < 1 2 and

D F i (c n+1 , c n ) = 4Σ T ε j =i 1 Σ j d F i (c n+1 , c n ) -d F j (c n+1 , c n ) .
The functions d F i represent a semi-implicit discretization of ∂ ci F . Since F = F 0 +P , we choose a discretization of the form d F i = d F0 i + d P i where d F0 i and d P i represent the discretization of ∂ ci F 0 and ∂ ci P respectively. We choose a fully implicit discretization of ∂ ci F 0 . It corresponds to the following definition of d F0 : for i ∈ {1, 2, 3},

d F0 i (c n+1 , c n ) = ∂ i F 0 (c n+1 ). ( 20
)
For ∂ ci P , we choose a semi-implicit discretization: for i ∈ {1, 2, 3},

d P i (c n+1 , c n ) = 2Λc n+1 i (c n j ) 2 (c n k ) 2 + 1 2 (c n+1 j ) 2 (c n k ) 2 + 1 2 (c n j ) 2 (c n+1 k ) 2 + (c n+1 j ) 2 (c n+1 k ) 2 . ( 21
)
We add the following Neumann boundary conditions for each of the unknown which ensures, in particular, that the volume of each phase is conserved along the time evolution:

∇c n+1 i • n = M 0 ∇µ n+1 i • n = 0, on ∂Ω.
Other boundary conditions of interest will be presented in [START_REF] Boyer | Numerical schemes for a three component Cahn-Hilliard model[END_REF].

Space discretization of problem (17)

We first give the variational formulation of the semi-discrete problem. ) ∈ (H 1 (Ω)) 4 so that ∀ν ∈ H 1 (Ω), we have, for i = 1 and 2,

       Ω c n+1 i -c n i ∆t ν dx = - Ω M 0 Σ i ∇µ n+1 i • ∇ν dx, Ω µ n+1 i ν dx = Ω D F i (c n+1 , c n )ν dx + Ω 3 4 Σ i ε α∇c n i + (1 -α)∇c n+1 i • ∇ν dx.
where

c n = (c n 1 , c n 2 , 1 -c n 1 -c n 2 )
. For the Galerkin approximation in space we introduce V h a finite element (possibly multilevel) approximation subspace of H 1 (Ω).

Problem 3.2 (Galerkin formulation). Given (c n 1h , c n 2h ) ∈ (V h ) 2 , find (c n+1 1h , c n+1 2h , µ n+1 1h , µ n+1 2h ) ∈ (V h ) 4 such that ∀ν h ∈ V h , we have, for i = 1, 2,        Ω c n+1 ih -c n ih ∆t ν h dx = - Ω M 0 Σ i ∇µ n+1 ih • ∇ν h dx, Ω µ n+1 ih ν h dx = Ω D F i (c n+1 h , c n h )ν h dx + Ω 3 4 Σ i ε α∇c n ih + (1 -α)∇c n+1 ih • ∇ν h dx.
where

c n h = (c n 1h , c n 2h , 1 -c n 1h -c n 2h ).
Remark 3.3. Note that we use the same approximation space V h for each of the unknowns c n+1 1h , c n+1 2h , µ n+1 1h and µ n+1 2h . Moreover, remark that in the formulation of Problem 3.2 and in Theorems 3.4 and 3.5, we assume that the approximation space V h does not depend of n. Because of adaptivity, this is not the case in practice (see Section 3.3.1).

Existence -Convergence

Theorem 3.4. Let (c n 1h , c n 2h ) ∈ (V h ) 2 . We assume that, for all i ∈ {1, 2, 3}, Σ i > 0. There exists a solution

(c n+1 1h , c n+1 2h , µ n+1 1h , µ n+1 2h ) ∈ (V h ) 4 for Problem 3.2. We define    c N ih (t, •) = t n+1 -t ∆t c n ih (•) + t -t n ∆t c n+1 ih (•) for t ∈ [t n , t n+1 ], µ N ih (t, •) = µ n+1 ih (•) for t ∈ [t n , t n+1 ].
The following convergence result holds.

Theorem 3.5. Let (c 1 , c 2 , µ 1 , µ 2 ) be the weak solution of the continuous problem [START_REF] Maubach | Local bisection refinement for n-simplicial grids generated by reflection[END_REF]. We assume that, for all i ∈ {1, 2, 3}, Σ i > 0 and that the finite element approximation spaces are such that

• inf ν h ∈V h |ν -ν h | H 1 (Ω) → h→0 0, for any ν ∈ H 1 (Ω).
• There exists a constant C such that ∀ν ∈ H

1 (Ω), Π V h 0 ν H 1 (Ω) C|ν| H 1 (Ω) where Π V h 0 is the L 2 (Ω)- projection on V h . Then, for i = 1 and 2, c N ih → c i in C 0 ([0, T ], L 2 (Ω)) strong, c N ih ⇀ c i in L ∞ ([0, T ], H 1 (Ω))weak- * , µ N ih ⇀ µ i in L 2 loc ([0, T ], H 1 (Ω))
weak. when (∆t, h) → (0, 0). These two theorems are not proved here; complete proofs will be presented in [START_REF] Boyer | Numerical schemes for a three component Cahn-Hilliard model[END_REF]. The conerstone of these proofs is a suitable discrete energy estimate. In the case where all Σ i are not assumed to be positive, this estimate is not satisfied for the implicit time discretization (20) of F 0 . Thus, other semi-implicit discretizations of this term (of the same kind that the given discretization (21) of P ) are introduced in [START_REF] Boyer | Numerical schemes for a three component Cahn-Hilliard model[END_REF] to ensure the validity of a such estimate and then to prove the convergence theorem in the case where we only assume that the condition (19) is fulfilled.

Practical issues

Newton linearization method

In practice, for the numerical resolution of the problem we choose for V h a multilevel approximation space as described in Section 1.3.2. Because of adaptivity, the unknowns are discretized in two different multilevel spaces at time step n and at time step n + 1, say V n h = span B n and V n+1 h = span B n+1 . Hence, the problem we have to solve is in fact:

Given (c n 1h , c n 2h ) ∈ (V n h ) 2 , find (c n+1 1h , c n+1 2h , µ n+1 1h , µ n+1 2h ) ∈ (V n+1 h ) 4 such that ∀ν h ∈ V n+1 h , we have, for i = 1, 2,        Ω c n+1 ih -c n ih ∆t ν h dx = - Ω M 0 Σ i ∇µ n+1 ih • ∇ν h dx, Ω µ n+1 ih ν h dx = Ω D F i (c n+1 h , c n h )ν h dx + Ω 3 4 Σ i ε α∇c n ih + (1 -α)∇c n+1 ih • ∇ν h dx.
This problem is solved by the Newton-Raphson linearization method. The algorithm is initialized by a L 2 -projection of (c n 1h , c n 2h , µ n 1h , µ n 2h ) on V n+1 h and at each time step we have to solve the following linear system:

    J c1c1 J c1µ1 J c1c2 0 J µ1c1 J µ1µ1 0 0 J c2c1 0 J c2c2 J c2µ2 0 0 J µ2c2 J µ2µ2          c (k+1) 1h -c (k) 1h µ (k+1) 1h -µ (k) 1h c (k+1) 2h -c (k) 2h µ (k+1) 2h -µ (k) 2h      =     -R c1 -R µ1 -R c2 -R µ2     , (22)
where (c

(k) 1h , c (k) 2h , µ (k) 1h , µ (k)
2h ) is the k-th iterate in the Newton iterative method and the matrix blocks are defined as follows:

∀(I, J) ∈ 1, N n+1 dof 2 , [J cici ] IJ = Ω - 3 4 Σ i ε(1 -α)∇ν J • ∇ν I dx - Ω ∂D i ∂c i (c (k) h , c n h )ν J ν I dx, [J cicj ] IJ = Ω - ∂D i ∂c j (c (k) h , c n h )ν J ν I dx, [J ciµi ] IJ = Ω ν J ν I dx, [J µi ci ] IJ = Ω 1 ∆t ν J ν I dx, [J µi µi ] IJ = Ω M 0 Σ i ∇ν J • ∇ν I dx, [R ci ] I = Ω µ (k) ih ν I dx - Ω D i (c (k) h , c n h )ν I dx - Ω 3 4 Σ i ε α∇c n ih + (1 -α)∇c (k) ih • ∇ν I dx, [R µi ] I = Ω c (k) ih -c n ih ∆t ν I dx + Ω M 0 Σ i ∇µ (k) ih • ∇ν I dx. (23)
This linear system [START_REF] Rivara | The 4-triangles longest-side partition of triangles and linear refinement algorithms[END_REF] is solved by the preconditionned GMRES iterative method. We use a diagonal bloc preconditioner:

    P c1c1 0 0 0 0 P µ1µ1 0 0 0 0 P c2c2 0 0 0 0 P µ2µ2    
where P c1c1 , P µ1µ1 , P c2c2 and P µ2µ2 are preconditioners built as explain in Section 2.2 (either additive or multiplicative version) associated respectively to the matrices A J = J c1c1 , J µ1µ1 , J c2c2 and J µ2µ2 .

Assembly avoiding field tranfers

Among the integrals in ( 23), some involve discrete unknowns which belong to two distinct approximation spaces; for example the integral

Ω c n ih ν I dx from [R µi ] I involves c n ih ∈ V n h and ν I ∈ V n+1 h
. Such an integral can be exactly computed avoiding field tranfers thanks to the following notions. Note that until the end of this section, the time step n + 1 is fixed. Definition 3.6 (Active DOF). We say that k ∈ {1, . . . , N

[j] dof } is an active degree of freedom of level-[j], if and only if ϕ

[j] k ∈ B n ∪ B n+1 . We denote by A [j]
dof the set of active degrees of freedom of level-[j]. Definition 3.7. Let j ∈ 0, J . We say that a level-[j] cell K [j] satisfies the (A K

[j] ) property iff ∀j ′ ∈ j + 1, J , ∀k ∈ A [j ′ ] dof , • K [j] ∩ • supp[ϕ [j ′ ] k ]= ∅. (A K [j] )
Remark 3.8. This definition ensures that on an active cell, all basis functions involved in the discretization of the problem, has a "unique" expression (opposed to an elementwise definition). Thus, in the case of polynomial finite elements (e.g. P k , Q k ), suitable quadrature rules can be used in order to exactly compute all the integrals. Definition 3.9 (Active cells). A level-[j] cell K [j] of T j is called active iff:

• K [j] satisfies the property (A K [j] ), and • in the case of j > 0, its parent cell P (K [j] ) (of level [j -1]) does not satisfy the property A P (K [j] ) .

With this definition at hand, we can introduce the notion of a multilevel mesh.

Proposition 3.10 (Multilevel mesh). Let T j be the set of active cells of T j . The set T = J j=0

T j is a mesh of Ω called multilevel mesh.

Proof. Let K e ′ are two distinct cells of the mesh T j . Then,

• K [j] e ∩ • K [j ′ ] e ′ = ∅. • Case 2 : j > j ′ .
Arguing by contradiction, assume that

• K [j] e ∩ • K [j ′ ] e ′ = ∅. Since j > j ′ , we have K [j] e ⊂ K [j ′ ] e ′ . But, the cell K [j]
e is active and ¬ A P (K [j] e ) yields that there exits j 0 j and k 0 ∈ A

[j0] dof such that • P (K [j] e ) ∩ • supp[ϕ [j0] k0 ] = ∅.
Moreover, we have j 0 > j ′ and A K

[j ′ ] e ′ yields • K [j ′ ] e ′ ∩ • supp[ϕ [j0]
k0 ]= ∅. However, we have (j > j ′ and K

[j] e ⊂ K [j ′ ] e ′ ) =⇒ • P (K [j] e )⊂ • K [j ′ ] e ′ . Hence, ∅ =   • P (K [j] e ) ∩ • supp[ϕ [j0] k0 ]   ⊂   • K [j ′ ] e ′ ∩ • supp[ϕ [j0] k0 ]
  = ∅. This is a contradiction.

Hence, the interiors of two different active cells are disjoint. Let x ∈ Ω. Since T J is a mesh of Ω, there exists a level-[J] cell K

[J]

eJ which contains x. Then, for all j ∈ 0, J -1 , we define K ). Hence, for all j ∈ 0, J , x belongs to the cell K ) . We have J ∈ E, so E = ∅. Let j m = min Remark 3.11. Multilevel meshes are not geometrically conformal, but multilevel finite element spaces are H 1 (Ω)-conformal by construction since span B ⊂ X J ⊂ H 1 (Ω). The non-conformity of multilevel meshes is not a problem since they are only used as elementary integration domains.

Since T is a partition of Ω, the integral can be decomposed into a sum over all cells of the multilevel mesh T on which a quadrature rule is then applied leading to exact computations (see Remark 3.8).

Refinement criterion

For the Cahn-Hilliard model, we want to refine in the interface zone (that is to say where the order parameters have important variations). For this reason, we need to introduce an indicator measuring whether or not a given cell (resp. basis function) lies in the interface at a given time step. Definition 3.12 (Per-cell indicator). At a given time step n, for any active cell K, we define the following indicator:

η K = max 1 |K| K c n 1 , 1 |K| K c n 2 , 1 |K| K c n 3
Remark 3.13. This per-cell indicator can be interpreted as follows:

• η K = 1 means that the cell K is completely filled with one of the bulk phases.

• η K < 1 means that the cell K contains an interface.

We can deduce a criterion for deciding whether or not a given basis function may be (un-)refined.

Definition 3.14 (Basis function indicator). At a given time step n, for a basis function ϕ, we define the following indicator:

η ϕ = 1 |supp[ϕ]| K∈T , K∩supp[ϕ] =∅ |K|η K
Criterion 3.15 ((Un)Refinement criterion). Given an anticipated cell size h interface for the interface neighborhood, the two following criterion let us decide if a basis function ϕ has to be refined or unrefined.

• Refinement criterion: η ϕ < 0.90 and diam(K) > h interface for at least one cell K ⊂ supp[ϕ].

• Unrefinement criterion: η ϕ > 0.95.

Numerical experiments

Finally, we give an application of the local refinement procedure and the multigrid methods to the resolution of the above ternary Cahn-Hilliard model in the case of a liquid lens spreading between two stratified fluids. Figure 21 shows the initial position of the interface and the numbering of phases used to defined the surface tensions σ ij between the phases i and j. Table 6 gives the parameters we used for numerical simulations. The refinement criterion is given in Section 3.3.3; h interface is specified for each mesh. Figures 22,[START_REF] Yserentant | Hierarchical bases[END_REF][START_REF] Zhang | Multi-level iteratives techniques[END_REF] shows three steps of the time marching. In each case, the value of the parameter h interface is voluntarily chosen large enough for the ease of visualisation.

As in Section 2. Table 7 illustrates that the multigrid process can be used to develop convergent iterative algorithms for the solution of the Cahn-Hilliard equations on local refined meshes. The simulations have been done during twenty time steps. For each time step, we calculate an average of the number of iterations required to solve the linear system by the GMRES method, over all iterations of the Newton algorithm. In Table 7, the numbers which face the preconditioner name are the minimum, the maximum and the average of this numbers over all the twenty time steps. For these simulations, we used h interface = 0.2ε. Note that in these simulations, we say that iteration 0 iteration 10 iteration 20 t = 0 t = 1. 10 -3 t = 2. 10 -3

Figure 24. Tri/Quadr-angle-P 1 /Q 1 . Time marching. h interface = 0.4ε convergence in the GMRES method is achieved as soon as the relative L 2 -norm of the residual (and not the preconditionned residual) is less than 10 -10 . In particular, this criterion is independent of the preconditioners. We observe that the number of necessary iterations is the almost the same from a time step to another and that the number of necessary iterations using (P m ) and (P a ) is significantly smaller than using (ILU0) preconditioner.

Min

Conclusion

This paper is devoted to the description and the resolution of some numerical issues linked to local adaptive refinement methods. We show how a hierarchy of nested conforming meshes can be built recursively applying a same refinement pattern to each cell of an initial mesh possibly unstructured. This hierarchy is then used in a slightly modified version of the CHARMS method [START_REF] Krysl | Natural hierarchical refinement for finite element methods[END_REF] in order to perform local refinement. Furthermore, the multilevel structure of the approximation spaces built with this method is exploited to derive suitable multigrid preconditioners for solving the corresponding linear system. We illustrate possibilities of the whole method with simultations of a lens spreading between two stratified liquid phases thanks to a ternary Cahn-Hilliard model. We particularly point out the practical way to implement the method for a time dependent problem. More complex situation have been considered [START_REF] Boyer | Cahn-Hilliard / Navier-Stokes model for the simulation of three-phase flows[END_REF] with simulations of three phases flows by solving the coupling of the Cahn-Hilliard system and the Navier-Stokes system. As an example, we give in Figure 25 the simulation of the rise of a gaz bubble across a light/heavy liquids interface. In future works, we will incorporate multigrid preconditioners in finite element penalty projection method for solving Navier-Stokes system [START_REF] Guermond | On the approximation of the unsteady Navier-Stokes equations by finite element projection methods[END_REF][START_REF] Jobelin | A finite element penalty-projection method for incompressible flows[END_REF] in order to build efficient solver for the complete model.

Figure 1 .

 1 Figure 1. Example of a geometrically non-conformal mesh.

Σ

  e . Let us denote by N dof = # Σ the number of Lagrange nodes, so that we can write Σ = {a i ; 1 i N dof }. • a set of basis functions {ϕ i ; 1 i N dof } defined as follows. First, for e ∈ 1, N e and k ∈ 1, N where N = # Σ, we denote by I(e, k) ∈ {1, . . . , N dof } the corresponding index, in the global numerotation, of the k-th local Lagrange node in the e-th cell of T . That is ∀e ∈ 1, N e , ∀k ∈ 1, N , a I(e,k) = T e ( a k ).

Figure 3 .

 3 Figure 3. Example of an incompatible refinement pattern.

Figure 6 .

 6 Figure 6. Q 2 -square refinement coefficients.

Definition 1 .

 1 21 (Parent-child relationship for cells). Let e ∈ 1, N [j] e .

BFigure 9 .

 9 Figure 9. Quadrangle-Q 1 . Nested meshes T j .

  km ) = 0. Hence, evaluating the linear combination (4) at a[jm] km yields λ ϕ [jm ] km = 0. This is a contradiction and the claim is proved.

2. 1 .

 1 Adaptation 2.1.1. Refinement/Unrefinement procedures Given a multilevel basis B, let us first introduce the notion of B-refined basis functions. Definition 2.1 (B-refined basis functions). Let B a multilevel basis. Let j ∈ 0, J and k ∈ 1, N [j] dof . The basis function ϕ [j]

k

  be a basis function belonging to B ⋆ . Refining the given basis function ϕ [j] k ∈ B ⋆ consists in producing a new multilevel basis B by removing this basis function ϕ [j] k , and adding all its children ϕ [j+1] ℓ which are not B ⋆ -refined.

  [j] k ∈ B ⋆ consists in producing a new multilevel basis B by adding this basis function ϕ [j] k , and removing those children of ϕ [j]

  [j+1] ℓ ∈ span B. Hence, ϕ [j] k ∈ span B and the recurrence is established. Theorem 2.8. Let B be a multilevel basis satisfying property (P LO ) and obtained from the multilevel basis B ⋆ by the refinement procedure (Algorithm 2.4), then span B ⋆ ⊂ span B.

  have to show that ϕ [j+1] ℓ either belong to B or is B-refined. • Refinement: Assume that the multilevel basis B is obtained from B ⋆ by the refinement of a basis function ϕ [j0]

and B =

 = B\{children of ψ which have no other B-refined parent} ∪ {ψ}[START_REF] Grinspun | CHARMS: A simple framework for adaptive simulation[END_REF] Applying Lemma 2.9 property 2) (i) and (ii), we get {B ⋆ -refined basis functions} = {B-refined basis functions} ∪ {ϕ}[START_REF] Guermond | On the approximation of the unsteady Navier-Stokes equations by finite element projection methods[END_REF] Hence, combining[START_REF] Grinspun | CHARMS: A simple framework for adaptive simulation[END_REF] and[START_REF] Guermond | On the approximation of the unsteady Navier-Stokes equations by finite element projection methods[END_REF], we obtain B = B\{children of ψ which have no B ⋆ -refined parent except possibly ϕ} ∪ {ψ}

Algorithm 2 .

 2 15 (Coarsening). Let B F a multilevel basis. Let j M = max {j ∈ 0, J ; B F ∩ B j = ∅} be the highest refinement level in B F . A "coarser" multilevel basis B C , denoted byB C = coarsen(B F ), is obtained from B F by unrefinement (Algorithm 2.4) of the set of B F -refined basis functions of level [j M -1].The following proposition gives an equivalent formulation of the above algorithm. Proposition 2.16. Assume that B F is a multilevel basis satisfying the following property: Any level-[j] basis function, j 1, which either belongs to B or is B-refined, has at least one B-refined parent . (P HI ) Let j M = max {j ∈ 0, J ; B F ∩ B j = ∅}. The multilevel basis B C = coarsen(B F ) defined in Algorithm 2.15 can be obtained by the following equivalent algorithm: • remove all level-[j M ] basis functions of B F , • add all B F -refined basis functions of level [j M -1].

  k0 has no B ⋆ -refined children.

Figure 14 .Figure 15 .Figure 16 .

 141516 Figure 14. Function H ε .

Figure 17 .Figure 18 .Figure 19 .

 171819 Figure 17. Square-Q 1 . Refined meshes.

Figure 20 .

 20 Figure 20. Refined mesh in 3D. One refinement step.

Problem 3 . 1 ( 1 , c n+1 2 , µ n+1 1 , µ n+1 2

 311212 Variational formulation). Given (c n 1 , c n 2 ) ∈ (H 1 (Ω)) 2 , find (c n+1

  ]e ′ be two distinct active cells and let us show that• K ] e ′ = ∅.• Case 1 : j = j ′ . Necessarily, we have e = e ′ . In this case, K [j] e and K [j ′ ]

  the set E = j ∈ 0, J , ∀ℓ j, K [ℓ] e ℓ satisfies (A K [ℓ] e ℓ

  ej m is active and contains x.

Figure 21 . 3 Figure 22 .

 21322 Figure 21. Initial position of interface. Numbering of phases.

3 . 3 ,- 3 Figure 23 .

 33323 Figure 23. Quadrangle-Q 1 . Time marching. h interface = 0.8ε

Figure 25 .

 25 Figure 25. Entrainment of an heavy liquid during a bubble rising.

  Definition 1.2 (Mesh [9, Def 1.49]). Let ω ⊂ R d (d = 1, 2 or 3) be a domain [9, Def 1.46]. A mesh T of ω is a set {K e ⊂ ω; 1 e N e } of compact, connected, Lipschitz subsets of ω with non-empty interior (called cells), such that

  Proposition 2.18. Let B F be a multilevel basis satisfying the property (P LO ) and B C = coarsen(B F ). The following embedding holds: span B C ⊂ span B F . Remark 2.19. Let us notice that, thanks to the refinement equation, it is straightforward to build the matrix representation in the corresponding basis of the natural embedding from span B C to span B F . Indeed, if ϕ ∈ B C , we have • either ϕ ∈ B F• or all the children of ϕ belong to B F and then, the refinement equation yields an expression of ϕ as a linear combination of elements in B F .Remark 2.20. Note that properties (P LO ) and (P HI ) used in this section are not so restrictive since they are preserved by Un/Refinement procedures and since it is straightforward to see that they are satisfied by the coarse basis B 0 which is used, in practice, to start the adaptation algorithm. (see alsoRemark 2.11)

Proof. Let ϕ ∈ B C . If ϕ ∈ B F then ϕ has been unrefined and so this is a B F -refined basis function of level [j M -1]. Owing Proposition 2.7, ϕ ∈ span B F .

Table 3 .

 3 Two dimensional test cases. Number of iterations in the conjugate gradient solver as a function of the number of unknowns.

	level		1	2	3	4	5	6
	Number of unknowns 893 3053 11021 41757 161233 633629
	No preconditioner	53	99	197	372	-	-
	Square-Q 1	ILU0	31	61	114	221	457	-
		P a	18	23	27	29	33	34
		P m	8	8	9	9	10	10
	Number of unknowns 935 3217 11609 43593 168347 660710
	No preconditioner	176 357	-	-	-	-
	Quadrangle-Q 1	ILU0	45	84	183	388	-	-
		P a	46	67	83	97	111	120
		P m	15	16	20	21	23	24
	Number of unknowns 869 2821 9893 37577 144969 569757
	No preconditioner	68 131	279	579	-	-
	Tri/Quadr-angle-P 1 /Q 1	ILU0	36	67	131	272	543	-
		P a	21	29	35	40	43	47
		P m	9	10	10	11	12	12
	Refinement level		1	2	3	4
		Number of unknowns 3653 12405 44485 166253
		No preconditioner	202	404	-	-
	Square-Q 2	ILU0		79	135	259	508
		P a		36	56	69	77
		P m		16	18	20	20
		Number of unknowns 3827 13129 46881 175126
		No preconditioner	595	-	-	-
	Quadrangle-Q 2	ILU0		112	195	389	-
		P a		83	127	158	178
		P m		27	31	33	34
		Number of unknowns 3569 11441 40293 151013
		No preconditioner	180	361	-	-
	Tri/Quadr-angle-P 2 /Q 2	ILU0		80	152	309	-
		P a		38	54	73	82
		P m		17	18	20	21

Table 4 .

 4 Two dimensional test cases. Number of iterations in the conjugate gradient solver as a function of the number of unknowns.

Table 5 .

 5 Three dimensional test case. Number of iterations in the conjugate gradient solver as a function of the number of unknowns.

		Number of unknowns 9942 63329 459063
		No preconditioner	45	224	-
	1	ILU0	28	55	109
		P a	21	31	35
		P m	10	13	14

Table 7 .

 7 Number of iterations in the GMRES solver. h interface = 0.2ε

	Max Average

B ⋆ -refined. (iii) Here, we assume thatj ⋆′ = j ⋆ + 1. If ϕ [j ⋆ ′ ] k ⋆ ′ ∈ B then ϕ [j ⋆ ] k ⋆ is only once B-refined and the claim is proved. Otherwise, ϕ [j ⋆ ′ ] k ⋆ ′ is a child of ϕ [j0]k0 . However,owing to Remark 2.2, ϕ[j ⋆ ] k ⋆ is the unique parent of ϕ [j ⋆ ′ ] k ⋆ ′ . Hence, we get ϕ [j ⋆ ]

Acknowlegments. The authors would like to thank the referees for their careful reading of the paper and their valuable remarks.

Adaptation procedure

Proposition 2.12. Let B ⋆ be a multilevel basis.

1) Let E ⋆ ⊂ B ⋆ . It is possible to refine successively all basis functions belonging to E ⋆ , hence producing a multilevel basis B which is independent of the order in which the basis functions were refined. We say that B is obtained from B ⋆ by the refinement of the set of basis functions E ⋆ . 2) Let F ⋆ be a set of only once B ⋆ -refined basis functions which have no B ⋆ -refined child. It is possible to unrefine successively all basis functions belonging to F ⋆ , hence producing a multilevel basis B which is independent of the order in which the basis functions were unrefined. We say that B is obtained from B ⋆ by the unrefinement of the set of basis functions F ⋆ .

Proof. In the two cases, we first have to prove that successive (un)refinements are possible and then that the obtained multilevel basis B is independent of the order in which the basis functions are (un)refined. 1)

The set E ⋆ \{ϕ} is included in the multilevel basis produced by the refinement of ϕ since, in this procedure, only ϕ is removed from B ⋆ . Hence, all basis functions in E ⋆ \{ϕ} can be then refined. -It is sufficient to prove that B is independent of the order in which the basis functions are refined in the case where # E ⋆ = 2, say E ⋆ = {ϕ, ψ}. Denote by B the multilevel basis obtained by the refinement of ϕ in B ⋆ and then by B the multilevel basis obtained by the refinement of ψ in B. By definition, we have

and

Applying Lemma 2.9 property 1) (i) and (ii), we get

Hence, combining ( 6) and ( 7), we obtain

Combining ( 5) and ( 8), we get B = B ⋆ ∪ {children of ϕ and ψ which are not B ⋆ -refined} \{ϕ, ψ}.

This expression shows that B do not depend of the order in which the basis functions ϕ and ψ were refined.

2)

-Let ϕ ∈ F ⋆ . Denote by B the multilevel basis obtained by the unrefinement of ϕ in B ⋆ . We have to prove that all basis functions of F ⋆ \{ϕ} can be unrefined in B,i.e. that any basis function of F ⋆ \{ϕ} is only once B-refined and has no B-refined children. Let ψ ∈ F ⋆ \{ϕ}. Owing to Lemma 2.9 property 2) (iii), since ψ is only once B ⋆ -refined basis function different from ϕ, ψ is only once B-refined. Furthermore, if ψ has a B-refined child then applying Lemma 2.9 property 2) (i), this child is also B ⋆ -refined. This a contradiction and the claim is proved. -It is sufficient to prove that B is independent of the order in which the basis functions are unrefined in the case where # F ⋆ = 2, say F ⋆ = {ϕ, ψ}. Denote by B the multilevel basis obtained from B ⋆ by the unrefinement of ϕ and then by B the multilevel basis obtained from B by the unrefinement of ψ. By definition, we have

Multigrid framework

We recursively define a sequence {V 0 , . . . , V J } of nested spaces built upon V h as follows:

• we first take B J = B and V J = span B J = V h ,

• then, for k = J, . . . , 1, we define a coarser multilevel basis B k-1 from B k by:

and the corresponding multilevel approximation space:

Owing to Proposition 2.18, we have:

Note that the auxiliary sequence V 0 ⊂ • • • ⊂ V J introduced here usually do not reflect the dynamic refinement process, although a such sequence can always be a posteriori deduced from any multilevel approximation space V J . An example with four refinement level is given in Figure 13.

V 3

Level- [START_REF] Bank | The hierarchical basis multigrid method[END_REF] basis functions Level- [START_REF] Bänsch | Local mesh refinement in 2 and 3 dimensions[END_REF] basis functions

Level- [START_REF] Bank | Some refinement algorithms and data structures for regular local mesh refinement[END_REF] basis functions Level-[0] basis functions

Example of coarsening : from V 3 to V 0 .

We do not need to know explicitly any information about the spaces V k , except the intergrids operators. Owing to Remark 2.19, it is staightforward to construct the matrix representation, denoted by I k k-1 , of the natural embedding from V k-1 to V k in the basis B k-1 and B k . Then, intergrids operators are defined in the following way, for all k ∈ 0, J :

We can also define approximate operators on each space V k , for all k ∈ 0, J by:

At last, in the sequel, we used Jacobi and Gauss-Seidel smoothers defined for all k ∈ 0, J as follows:

• Jacobi: S k = D k where D k is the diagonal part of A k .

• Gauss-Seidel: S k = T k where T k is the superior triangular part of A k .

Multigrid Algorithm

We used the two following multigrid preconditioners: