Mustapha Bahari

Abderrahmane Habbal

AND Ahmed Ratnani

Eric Sonnendr Ücker

ADAPTIVE ISOGEOMETRIC ANALYSIS USING OPTIMAL TRANSPORT

Keywords: Mathematics Subject Classification. Primary B-spline discretization, Optimal Transport, Isogeometric Analysis, multigrid, mesh generation, r-refinement. ©XXXX American Mathematical Society

Introduction Motivations.

Equidistribution mesh generation using Monge-Kantorovich Problem. Problem formulation Our contributions 2. B-Splines and geometrical representation 3. Equidistribution mesh generation based on Monge-Ampère equation 3.1. The Monge-Kantorovich Problem 3.2. New Mesh Equidistribution Technique 4. Benamou-Froese-Oberman method 4.1. Standard Monge Ampère Solver 4.2. Variational Formulation and A priori error estimate 4.3. Fast Diagonalization method 5. Numerical Results 5.1. Monge-Ampère solver using BFO method 5.2. Adaptive mesh generation 6. Mixed Finite Elements approach 6.1. Mixed formulation of Benamou-Froese-Oberman system 6.2. Fast Diagonalization method for Mixed Variational Formulation 7. Numerical Results 7.1. Monge-Ampère solver using mixed variational formulation 7.2. Adaptive mesh generation 8. Conclusions Current limitations and future work References

Introduction

The use of adaptive mesh methods is fundamental to the numerical solution of systems of partial differential equations that involve large solution variations or with different scales. Adaptive mesh methods are also a major tool for solving problems with high anisotropy, such those encountered in Computational Plasmas Physics. Grid Adaptation and moving meshes are used in different fields and the litterature is quite rich [START_REF] Cao | Approaches for generating moving adaptive meshes: location versus velocity[END_REF][START_REF] Huang | An stability analysis for the finite-difference solution of one-dimensional linear convection-diffusion equations on moving meshes[END_REF][START_REF] Huang | A study of moving mesh pde methods for numerical simulation of blowup in reaction diffusion equations[END_REF][START_REF] Huang | Variational mesh adaptation ii: error estimates and monitor functions[END_REF][START_REF]Adaptive mesh movement -the mmpde approach and its applications, Partial Differential Equations[END_REF][START_REF]Variational mesh adaptation: Isotropy and equidistribution[END_REF][START_REF] Huang | Practical aspects of formulation and solution of moving mesh partial differential equations[END_REF][START_REF]A high dimensional moving mesh strategy[END_REF][START_REF] Huang | A moving collocation method for solving time dependent partial differential equations[END_REF][START_REF] Huang | Moving mesh methods based on moving mesh partial differential equations[END_REF][START_REF] Lang | A two-dimensional moving finite element method with local refinement based on a posteriori error estimates[END_REF][START_REF] Budd | Adaptivity with moving grids[END_REF]. They are also used to generate Anisotropic meshes [START_REF] Li | An anisotropic mesh adaptation method for the finite element solution of heterogeneous anisotropic diffusion problems[END_REF][START_REF] Lu | The cutoff method for the numerical computation of nonnegative solutions of parabolic pdes with application to anisotropic diffusion and lubrication-type equations[END_REF][START_REF] Huang | A new anisotropic mesh adaptation method based upon hierarchical a posteriori error estimates[END_REF][START_REF] Huang | An anisotropic mesh adaptation method for the finite element solution of variational problems[END_REF][START_REF]Metric tensors for anisotropic mesh generation[END_REF]. The proposed methods in this work aim to construct a one-to-one mapping F that maps a logical domain (patch, computational domain) with the physical domain as shown in Fig. 1. The function F is constructed using B-splines or NURBS surfaces which are widely used in the Computer Aided Design (CAD) community. The use of these tools in numerical simulations was made popular thanks to the introdution of the IsoGeometric Analysis paradigm by Hughes [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF]. Because of the geometric interpretation of the control points, B-spline curves and surfaces have become very popular in CAD. We are interested in these features as they allow us to construct a set of mappings, where each of them will map the unit square onto a sub-domain of the physical domain. Using the geometric properties of B-splines curves and surfaces, it is easy to stick these mappings together, in order to have a global C 1 or even C 2 mapping. Local regularity of each mapping is ensured by construction.

Motivations. Many physical phenomena can be numerically simulated by the resolution of the associated partial differential equations. However, many numerical methods suffer from the local singularity of the solution of these equations. Such as shock waves, boundary layer in compressible flow problems. Thus, it is very interesting and advantageous to use a mesh that complies to the behavior of each phenomenon, and of course, because the refinement of a uniform mesh directly consumes a lot of energy the local refinement is mandatory.

Recently, using moving mesh methods, some authors are generally addressed the application to the field of weather prediction and climate simulation, which represent small-scale problems [START_REF] Andrew Tt Mcrae | Optimal-transport-based mesh adaptivity on the plane and sphere using finite elements[END_REF]; thus, based on the optimal transport problem the process requires the resolution of the Monge-Ampere equation on the sphere domain (earth shape), the investigation in this direction is given in the papers [START_REF] Weller | Mesh adaptation on the sphere using optimal transport and the numerical solution of a monge-ampère type equation[END_REF][START_REF] Andrew Tt Mcrae | Optimal-transport-based mesh adaptivity on the plane and sphere using finite elements[END_REF][START_REF] Froese | A convergence framework for optimal transport on the sphere[END_REF]. Moreover, the global adaptive meshes and conservation of the connection between the grid points are the most attractive properties of optimal transport problem [START_REF] Shakti | A moving mesh refinement based optimal accurate uniformly convergent computational method for a parabolic system of boundary layer originated reaction-diffusion problems with arbitrary small diffusion terms[END_REF][START_REF] Khosravi | Three-dimensional features of the lateral thermal plume discharge in the deep cross-flow using dynamic adaptive mesh refinement[END_REF].

Equidistribution mesh generation using Monge-Kantorovich Problem. The idea behind the equidistribution is very simple. It aims to generate a grid that equidistributes a given quantity (a density or a monitor function) along a surface (2D) or a volume (3D). More details about the equidistribution and minimization of total error can be found in [START_REF] Lapenta | Variational grid adaptation based on the minimization of local truncation error: time-independent problems[END_REF], [START_REF] Baines | Least squares and approximate equidistribution in multidimensions[END_REF]. An interesting introduction to equidistribution grid adaptation based on the Monge-Kantorovich problem (MKP) can be found in [START_REF] Delzanno | An optimal robust equidistribution method for two-dimensional grid adaptation based on monge-kantorovich optimization[END_REF] and the references therein. The connection between MKP and the resolution of the Monge-Ampère equation (MA) is now quite well-known [START_REF] Delzanno | An optimal robust equidistribution method for two-dimensional grid adaptation based on monge-kantorovich optimization[END_REF][START_REF] Huang | Adaptive moving mesh methods[END_REF]. Recently, many authors proposed equidistribution grid generation methods using the MA [START_REF] Delzanno | An optimal robust equidistribution method for two-dimensional grid adaptation based on monge-kantorovich optimization[END_REF][START_REF] Luca Delzanno | The fluid dynamic approach to equidistribution methods for grid adaptation[END_REF][START_REF] Chacòn | Robust, multidimensional mesh-motion based on monge-kantorovich equidistribution[END_REF][START_REF] Sulman | An efficient approach for the numerical solution of the monge-ampère equation[END_REF][START_REF] Sulman | Optimal mass transport for higher dimensional adaptive grid generation[END_REF][START_REF] Cossette | A monge-ampère enhancement for semilagrangian methods[END_REF][START_REF] Budd | Monge-ampère based moving mesh methods for numerical weather prediction, with applications to the eady problem[END_REF][START_REF]Numerical simulation of four-field extended magnetohydrodynamics in dynamically adaptive curvilinear coordinates via newton-krylov-schwarz[END_REF][START_REF] Yuan | Moving grids for magnetic reconnection via newton-krylov methods[END_REF] The numerical solution of the elliptic Monge-Ampère equation has been a subject of increasing interest. Many methods have been proposed, but only few of them converge for singular solutions. In this work, we use the Benamou-Froese-Oberman (BFO) method that we have implemented using the IGA approach. More details can be found in [START_REF] Benamou | Two numerical methods for the elliptic monge-ampère equation[END_REF]. In [START_REF] Edward | Numerical solution of the two-dimensional elliptic monge-ampère equation with dirichlet boundary conditions: an augmented lagrangian approach[END_REF], the authors solve the MA using an augmented Lagrangian approach. In [START_REF] Dean | Numerical solution of the two-dimensional elliptic monge-ampère equation with dirichlet boundary conditions: a least-squares approach[END_REF], they proposed a least square approach. Much more details can be found in [START_REF] Dean | An optimal robust equidistribution method for two-dimensional grid adaptation based on mongekantorovich optimization[END_REF]. In [START_REF] Loeper | Numerical solution of the monge-ampère equation by a newton's algorithm[END_REF], the authors proposed an iterative solver based on the divergence form and Newton method. While, in [START_REF] Delzanno | An optimal robust equidistribution method for two-dimensional grid adaptation based on monge-kantorovich optimization[END_REF] the authors used an iterative Newton-Krylov solver with preconditioning coupled with a Finite Difference method. In [START_REF]Spline element method for the monge-ampère equation[END_REF] [1], the author solves the Monge-Ampère equation as the limit of the solution of a singular perturbation problem, using triangular B-splines, which by the way can treat complex geometries. Another way of solving the Monge-Ampère equation, is to consider it as the steady state of a parabolic equation namely the Parabolic Monge-Ampère equation. This has been studied in [START_REF] Budd | Parabolic monge¿ampère methods for blow-up problems in several spatial dimensions[END_REF][START_REF] Sulman | An efficient approach for the numerical solution of the monge-ampère equation[END_REF][START_REF] Sulman | Optimal mass transport for higher dimensional adaptive grid generation[END_REF]. In [START_REF] Zheligovsky | The monge-ampère equation: Various forms and numerical solution[END_REF], the authors proposed three new methods based on Fourier integral, secondorder divergence and a convolution forms.

Problem formulation. Now, we consider a fixed computational domain Ω c , and Ω p the physical domain where the underlying PDE is posed and the initial mapping F : Ω c ξ -→ x ∈ Ω p . The main idea behind an optimum mesh is that it should be closest to a uniform mesh in a suitable norm [START_REF] Huang | Adaptive moving mesh methods[END_REF]. The equidistributed grids are constructed by solving the Monge-Ampere equation which is equivalent to the L 2 Monge-Kantorovich problem [18,[START_REF] Huang | Adaptive moving mesh methods[END_REF].

Hence, the adaptive meshes are obtained using a new mapping F : Ω c ξ -→ x ∈ Ω p , that maps the parametric domain Ω c into the physical domain Ω p , defined as the composition of initial mapping F : Ω c ξ -→ x ∈ Ω p and B-spline geometrical map Ψ : Ω c ξ -→ ξ ∈ Ω c derived from a resolution of Monge-Ampere equation. Indeed, for our application, we have Ω c = (0, 1) d , d = 2 or 3. We give the illustration of these ideas in Fig. 1.

Our contributions. Our contributions are the following:

(1) A new technique of mesh adaptation illustrated in Fig. 1, (2) A fast solver for Monge-Ampère equation, using fast diagonalization method, (3) Demonstrate the effectiveness of our method on complex geometry, (4) Derive a priori estimates.

The remainder of this paper is as follows. In Section 2, we describe the basic theory behind B-splines curves and surfaces. In Section 3, we describe the basic theory behind adaptive meshes methods and specially the optimal transport problem which is linked to the Monge-Ampere equation. In Section 4-5, we use the Benamou-Froese-Oberman (BFO) method to generate equidistributed meshes. In order to enhance the performance of the Picard algorithm, a multi grids method is used. In Section 6-7, we introduce a mixed Aligned-Equidistributed strategy to handle the boundary conservation problem and limitation of BFO method. With high accuracy and fast convergence using the multi grids method, two algorithms of the fast diagonalization method in two and three dimensions are developed for the mixed variational formulation. Concluding remarks are given in section 8.

B-Splines and geometrical representation

IsoGeometric analysis (IGA) has been introduced in [START_REF] Thomas | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF], the main idea was to improve the interoperability between Computer-Aided Design (CAD) and partial differential equation (PDE) solvers. Let's start by recalling the definition and some aspects of isogeometric analysis. B-splines. Given two positive integers p and n, we introduce the (ordered) knot vector

(2.1) Ξ = {0 = ξ 1 , ξ 2 , ..., ξ n+p+1 = 1}
where n is the number of basis functions necessary to describe it.

(2.3) B 0 i (ξ) = 1 if ξ ∈ [ξ i , ξ i+1 ) 0 otherwise B p i (ξ) = ξ -ξ i ξ i+p -ξ i B p-1 i (ξ) + ξ i+p+1 -ξ ξ i+p+1 -ξ i+1 B p-1 i+1 (ξ) (2.4)
where in 2.4, we adopt the convention 0/0 = 0. These basis functions are piecewise polynomials of degree p on the subdivision {ζ 1 , .., ζ m }, where at ζ j they have µ j := p -r j continuous derivative. Therefore, the vector µ = {µ 1 , ..., µ m } collects the regularity at the internal knots, with µ 1 = µ m = -1 for the boundary knots associated with the open knot vector.

An example of cubic B-splines is presented in Fig. 2. In this case µ = {-1, 2, 1, 0, -1}.

Respectively, by means of tensor products, a multi-dimensional B-spline can be constructed as B p1,p2,...,p d j1,j2,...,j d = ⊗ d i=1 B pi ji [START_REF] Thomas | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF][START_REF] Cottrell | Isogeometric analysis: toward integration of cad and fea[END_REF]. Geometrical representation. A single patch domain Ω p is a B-spline region associated with the control points C j1,j2,...,j d . We introduce the B-spline geometrical map F : Ω c -→ Ω p by (2.5) F(ξ) = j1,j2,...,j d B p1,p2,...,p d j1,j2,...,j d (ξ)C j1,j2,...,j d .

For our purpose we assume that the geometry mapping is continuous and bijective which is a natural assumption for CAD applications. We introduce Q h a family of meshes on Ω c by (2.6)

Q h = {Q = ⊗ d i=1 (ζ ji,i , ζ ji+1,i ), 1 ≤ j i ≤ m i -1, d = 2, 3} each element Q ∈ Q h is mapped into an element (2.7) k = F(Q) = {F(ξ), ξ ∈ Q},
and analogously Q, the support extension of Q, is mapped into

(2.8) k = F( Q).
We then introduce the mesh T h in the physical domain Ω p , defined by:

(2.9)

T h := {k = F(Q), Q ∈ Q h }.
Finally, we define B-spline spaces needed in the sequel as follows:

(2.10)

S p µ = span{B p i , for i = 1, ..n} and (2.11) dv dξ : v ∈ S p µ = S p-1 µ-1
where we adopt the notation µ -

1 = {µ 1 , µ 2 -1, ..., µ m-1 -1, µ m }.
Then, we consider the discrete subspace of H s (Ω c ), given by (2.12)

V h = ⊗ d i=1 S pi µi .
Let us introduce the discrete subspace of H(div, Ω) and L 2 (Ω) denoted by V h (div, Ω) and V h (L 2 , Ω) successively by -For 2D case:

         V h (div, Ω) = S p1 µ1 ⊗ S p2-1 µ2-1 S p1-1 µ1-1 ⊗ S p2 µ2 V h (L 2 , Ω) = S p1-1 µ1-1 ⊗ S p2-1 µ2-1 (2.13)
-For 3D case:

             V h (div, Ω) =    S p1 µ1 ⊗ S p2-1 µ2-1 ⊗ S p3-1 µ3-1 S p1-1 µ1-1 ⊗ S p2 µ2 ⊗ S p3-1 µ3-1 S p1-1 µ1-1 ⊗ S p2-1 µ2-1 ⊗ S p3 µ3    V h (L 2 , Ω) = S p1-1 µ1-1 ⊗ S p2-1 µ2-1 ⊗ S p3-1 µ3-1 (2.14)
The functional spaces that we choose are because they form an exact sequence (called DeRham sequence) needed for the mixed variational formulation. So, we recall that in 2d, we have

V h (grad, Ω) ∇× ---→ V h (div, Ω) ∇• ---→ V h (L 2 , Ω)
While in 3d, we have the following sequence

V h (grad, Ω) ∇ --→ V h (curl, Ω) ∇× ---→ V h (div, Ω) ∇• ---→ V h (L 2 , Ω)
where

V h (curl, Ω) =    S p1-1 µ1-1 ⊗ S p2 µ2 ⊗ S p3 µ3 S p1 µ1 ⊗ S p2-1 µ2-1 ⊗ S p3 µ3 S p1 µ1 ⊗ S p2 µ2 ⊗ S p3-1 µ3-1    (2.15)
is the discrete subspace of H(curl, Ω).

Equidistribution mesh generation based on Monge-Ampère equation

3.1. The Monge-Kantorovich Problem. Definition 3.1 (L 2 Monge-Kantorovich problem). Let ρ 0 and ρ 1 be two given densities of equal masses, defined in R d . Find a mapping x = T(x), x, x ∈ R d , that transfers the density ρ 0 to ρ 1 and minimizes the transport cost

J [T] = R d |T(x) -x | 2 ρ 0 (x) dx (3.1)
The density transfer means that, for all Borelian set B in R d

T -1 (B) ρ 0 (x) dx = B ρ 1 (x) dx.
Then, it follows from a change of variable that

ρ 1 (T) det(∇T) = ρ 0 .
Brenier's theorem [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF], shows that T = ∇φ for given convex potential function φ. Hence, we get the Monge-Ampere equation:

(3.2) ρ 1 (∇φ) det(H(φ)) = ρ 0 ,
where H(φ) is the Hessian matrix of φ.

Remark 3.1. For arbitrary ρ 0 and ρ 1 defined in the physical domain Ω p , let

(3.3) σ = Ωp ρ 1 (x) dx/ Ωp ρ 0 (x) dx,
we rewrite the Monge-Ampere equation as follows

(3.4) ρ 1 (∇φ) det(H(φ)) = σ.
3.2. New Mesh Equidistribution Technique. Since both density functions are defined in the physical domain, we explain in the following the mathematical background of the technique presented in Fig. 1.

Proposition 3.1. In terms of mesh generation, we assume that the initial mapping F is bijective, so we take the two densities defined in Ω c as follows (1) σ det(∇F): corresponding to a uniform mesh in Ω c , (2) ρ 1 • F det(∇F): corresponding to what we call the inverse image by the initial mapping F of the adapted mesh in Ω p , where σ is given by 3.3, the two densities are of equal masses. Then, according to Brenier's theorem [START_REF]Polar factorization and monotone rearrangement of vector-valued functions[END_REF] there exists an optimal mapping Ψ defined in the parametric domain Ω c which is the gradient of a given convex potential function φ, such that

(3.5) ρ 1 (F • ∇ φ) det(∇F) • ∇ φ det(H( φ)) = σ det(∇F).
Moreover, our mesh adaptation technique and the one existing in the literature [18] are equivalent.

Proof. Using ρ 0 = 1 correspond to uniform mesh in Ω p and assuming that Ω p is a convex domain, then there exists T : Ω p -→ Ω p satisfying the Monge-Ampere equation 3.4, which corresponds to the adaptation in the physical domain, hence we have

(3.6) ρ 1 (T) det(∇T) = σ.
using the fact that F is bijective, let T 0 : Ω c -→ Ω c such that

(3.7) T = F • T 0 • F -1 ,
as a consequence, by a change of variable in 3.6, we have

(3.8) ρ 1 (F • T 0 • F -1 ) det ∇(F • T 0 • F -1 ) = σ using the fact that ∇ F • T 0 • F -1 = ∇F -1 T ∇ F • T 0 • F -1 = ∇F -1 T ∇T 0 T • F -1 ∇F • T 0 • F -1 ,
and det(AB) = det(A) det(B), we get

(3.9) ρ 1 (F • T 0 • F -1 ) det (∇F -1 ) T det (∇T 0 ) T • F -1 det ∇F • T 0 • F -1 = σ
and therefore, by applying F we have

(3.10) ρ 1 (F • T 0 ) det (∇F -1 ) • F det ∇T 0 det ∇F • T 0 = σ.
From differential calculus, we know that (∇F -1 )

• F = ∇F -1 , then (3.11) ρ 1 (F • T 0 ) det ∇F • T 0 det ∇T 0 = σ det ∇F .
Finally, from the uniqueness of the optimal mapping we conclude that T 0 = Ψ.

Corollary 3.1. Assuming the initial mapping is uniform, we suppose that

(3.12) det ∇F • ∇ φ ≈ det ∇F ,
then, we simplify the equation 3.5 by considering the new form:

(3.13) ρ 1 (F • ∇ φ) det(H( φ)) = σ.
The boundary conditions are obtained by requiring that Ψ maps ∂Ω c to ∂Ω c , i.e ∇ φ(∂Ω c ) = ∂Ω c which writes:

(3.14) ∇ φ • n = x • n.
Finally, we denote

(3.15) f = σ ρ 1 (F) ,
the desired convex potential function is the solution of Monge-Ampere system:

(3.16)    det(H( φ)) = f (∇ φ) in Ω c , ∇ φ • n = x • n on ∂Ω c , φ is convex.
A solution of 3.16 exists and is unique up to an additive constant [START_REF] Lions | The neumann problem for equations of monge-ampère type[END_REF]. From now on, the notation φ stands for φ solution of the latest system 3.16 and ρ stands for ρ 1 .

4. Benamou-Froese-Oberman method 4.1. Standard Monge Ampère Solver. We adopt the second semi-implicit method proposed by Benamou Froese and Oberman [START_REF] Benamou | Two numerical methods for the elliptic monge-ampere equation[END_REF] which guarantees convergence when a second member is strictly positive and regular. The authors noticed that a solution of the Monge-Ampere equation is a fixed point of the operator T .

T : H 2 (Ω c ) -→ H 2 (Ω c ) v -→ T (v) = (∆) -1 G(v) (4.1)
where

G(v) = (∆v) d + d!(f (∇v)) -det(H(v))
   -∆φ = -G(φ) in Ω c , ∇φ • n = x • n on ∂Ω c , φ is convex.
Thus, the solution is given by Picard iteration as follow:

(1) Given an initial value φ 0 .

(2) Compute φ n+1 as the solution of the problem:

(4.3)    -∆φ n+1 + εφ n+1 = -G(φ n ) in Ω c , ∇φ n+1 • n = x • n on ∂Ω c , φ n is convex,
where ε = 10 -8 is a penalization parameter. (3) Repeat the iteration number (2) until L 2 residual is smaller than a given tolerance.

Remark 4.1.

-Firstly, we note that the initial guess for the Picard algorithm can be computed using a multi-grids (MG) method.

-We have tested different methods that exist in the literature such as the Lagrange multiplier, the Newton method, and the projection method, but we obtain the same convergence results as for the penalization method. Then, for reasons of simplicity and fast convergence, the penalization method is used for a robust numerical resolution of the Poisson equation with pure Neumann boundary conditions.

4.2.

Variational Formulation and A priori error estimate. The goal of this subsection is to validate our solver by two examples in both a two and threedimensional cases.

The variational formulation of 4.3, is:

(4.4) Ωc ∇φ n+1 ∇v dx+ε Ωc φ n+1 v dx = - Ωc G(φ n )v dx+ ∂Ωc x• n v dσ, for all v ∈ H 2 (Ω c )
and a Galerkin approximation of 4.3 is given by (4.5)

Ωc ∇φ n+1 h ∇v h dx+ε Ωc φ n+1 h v h dx = - Ωc G(φ n h )v h dx+ ∂Ωc x• n v h dσ, for all v ∈ V h ∩H 2 (Ω c )
where V h is given by 2.12. Without losing the generalization, we neglect the penalization part of the variational form in the following estimation. Proposition 4.1. Let {φ n } n be a solution sequence of Picard algorithm in H 2 (Ω c ), then it holds that:

(4.6) |e n,h | H 1 (Ωc) ≤ Ch l-1 ,
where e n,h = φ n -φ n h is the approximation error at iteration n and l is a Sobolev regularity exponent of

φ n φ n ∈ H l (Ω c ), l ≥ 2 .
Proof. By recurrence, we have for n = 0

|φ 0 -φ 0 h | H 1 (Ωc) ≤ Ch l-1
, holds by choice. So, we first prove that we have the following implication (4.7)

|φ n -φ n h | H 1 (Ωc) ≤ Ch l-1 =⇒ ||G(x, φ n ) -G(x, φ n h )|| L 2 (Ωc) ≤ Ch l-2 , for n ∈ N, recalling that ||H(φ n )|| 2 = (∆φ n ) d -d! det(H(φ n )) ≥ 0,
and, the triangular inequality leads to

|G(x, φ n ) -G(x, φ n h )| 2 = ||H(φ n )|| 2 + d! f (∇φ n ) -||H(φ n h )|| 2 + d! f (∇φ n h ) 2 ≤ 2 ||H(φ n )|| 2 + d! f (∇φ n ) -||H(φ n )|| 2 + d! f (∇φ n h ) 2 + 2 ||H(φ n )|| 2 + d! f (∇φ n h ) -||H(φ n h )|| 2 + d! f (∇φ n h ) 2 := 2A 2 + 2B 2 .
Using the inequality

√ a + b - √ c + b ≤ | √ a - √ c| for a, b, c ≥ 0, we have (4.8) B ≤ ||H(φ n )|| -||H(φ n h )|| and A ≤ √ d! f (∇φ n ) -f (∇φ n h ) ≤ √ d! √ σ ρ(∇φ n ) - √ σ ρ(∇φ n h ) ≤ √ d!σ ρ(∇φ n ) -ρ(∇φ n h ) ρ(∇φ n )ρ(∇φ n h ) since that ρ ≥ γ > 0, and t -→ √ t is 1-Lipschitz function in [1, +∞) then A ≤ d!σ γ ρ(∇φ n ) -ρ(∇φ n h ) if we assume that ρ is L ρ -Lipschitz (4.9) A ≤ L ρ d!σ γ ∇φ n -∇φ n h .
As a consequence of 4.8 and 4.9

||G(x, φ n ) -G(x, φ n h )|| L 2 (Ωc) ≤ 2|φ n -φ n h | H 2 (Ωc) + 2L ρ d!σ γ |φ n -φ n h | H 1 (Ωc) .

Now, taking the Projector Π

V h : L 2 (Ω c ) -→ V h , we get |φ n -φ n h | H 2 (Ωc) ≤ |(φ n -φ n h ) -Π V h (φ n -φ n h )| H 2 (Ωc) + |Π V h (φ n -φ n h )| H 2 (Ωc) = I + Π.
With C 2 regularity in B-spline function, from Proposition 3.2 [START_REF] Tagliabue | Isogeometric analysis and error estimates for high order partial differential equations in fluid dynamics[END_REF]

I = |(φ n -φ n h ) -Π V h (φ n ) + φ n h | H 2 (Ωc) = |φ n -Π V h (φ n )| H 2 (Ωc) ≤ Ch l-2 ||φ n || H l (Ωc) .
By the usual inverse inequality for polynomials, we have

Π 2 = k∈T h |Π V h (φ n -φ n h )| 2 H 2 (k) ≤ C shap k∈T h h -2 K |Π V h (φ n -φ n h )| 2 H 1 (k) ≤ C shap h -2 |φ n -φ n h | 2 H 1 (Ωc) , (4.10)
where h = min(h k ), so if we assume that for a positive constant θ, h ≤ θh.

According to the hypothesis of 4.7, we have

I + II ≤ Ch l-2 ,
where C depend only on |φ n | H l (Ωc) and θ.

Returning now to the proof of our proposition, let φ n+1 be the solution at n + 1 iteration in Picard algorithm and v h ∈ V h , we have for some constant c n+1

|φ n+1 -φ n+1 h | 2 H 1 (Ωc) = Ωc ∇(φ n+1 -φ n+1 h )∇(φ n+1 -v h + c n+1 ) dx + Ωc ∇(φ n+1 -φ n+1 h )∇(v h -c n+1 -φ n+1 h ) dx = Ωc ∇(φ n+1 -φ n+1 h )∇(φ n+1 -v h ) dx + Ωc (G(x, φ n h ) -G(x, φ n ))(v h -c n+1 -φ n+1 h
) dx, by Holder's inequality, we have

|φ n+1 -φ n+1 h | 2 H 1 (Ωc) ≤ |φ n+1 -φ n+1 h | H 1 (Ωc) |φ n+1 -v h | H 1 (Ωc) + ||G(x, φ n ) -G(x, φ n h )|| L 2 (Ωc) ||v h -c n+1 -φ n+1 h || L 2 (Ωc) ,
taking to account 4.7, we get by the triangular inequality

|φ n+1 -φ n+1 h | 2 H 1 (Ωc) ≤ |φ n+1 -φ n+1 h | H 1 (Ωc) |φ n+1 -v h | H 1 (Ωc) + Ch l-2 ||v h -c n+1 -φ n+1 h || L 2 (Ωc) ≤ |φ n+1 -φ n+1 h | H 1 (Ωc) |φ n+1 -v h | H 1 (Ωc) + Ch l-2 ||φ n+1 -v h || L 2 (Ωc) + Ch l-2 ||φ n+1 -φ n+1 h -c n+1 || L 2 (Ωc)
where v h ∈ V h , according to Proposition 3.2 [START_REF] Tagliabue | Isogeometric analysis and error estimates for high order partial differential equations in fluid dynamics[END_REF], for v h = Π V h (φ n+1 ) we have

|φ n+1 -Π V h (φ n+1 )| H 1 (Ωc) ≤ Ch l-1 ,
and

||φ n+1 -Π V h (φ n+1 )|| L 2 (Ωc) ≤ Ch l .
Then,

|φ n+1 -φ n+1 h | 2 H 1 (Ωc) ≤ Ch l-1 |φ n+1 -φ n+1 h | H 1 (Ωc) + Ch 2l-2 + Ch l-2 ||φ n+1 -φ n+1 h -c n+1 || L 2 (Ωc) . (4.11)
Let's now introduce the following problem:

(4.12) -∆ψ = φ n+1 -φ n+1 h -c n+1 in Ω c , ∇ψ • n = 0 on ∂Ω c .
The compatibility condition holds for c n+1 = Ωc φ n+1 -φ n+1 h dx, and the solution exists and unique up to an additive constant. Moreover, the variational formulation of 4.12 is written as follows :

(4.13) Ωc ∇ψ∇v dx = Ωc (φ n+1 -φ n+1 h -c n+1 )v dx for all v ∈ H 1 (Ω c ) by assuming that (4.14) Ωc (G(x, φ n ) -G(x, φ n h )) dx = 0 we introduce the constant c such that (4.15) Ωc (G(x, φ n ) -G(x, φ n h ))(Π V h (ψ) -c) dx = 0.
Hence, from the variational formulation 4.5, for

v h = Π V h (ψ) -c, we have Ωc ∇(φ n+1 -φ n+1 h -c n+1 )∇(Π V h (ψ) -c) dx = Ωc (G(x, φ n ) -G(x, φ n h ))(Π V h (ψ) -c) dx = 0. (4.16)
As a consequence, by replacing the results in 4.13, we get for

v = φ n+1 -φ n+1 h -c n+1 Ωc (φ n+1 -φ n+1 h -c n+1 ) 2 dx = Ωc ∇ψ∇(φ n+1 -φ n+1 h ) dx (4.17) = Ωc ∇(ψ -(Π V h (ψ) -c))∇(φ n+1 -φ n+1 h ) dx (4.18) ≤ |ψ -(Π V h (ψ) -c)| H 1 (Ωc) |φ n+1 -φ n+1 h | H 1 (Ωc) . (4.19) Furthermore, because φ n+1 -φ n+1 h ∈ H l (Ω c ) at least ψ ∈ H l+2 (Ω c
), then from the Proposition 3.2 [START_REF] Tagliabue | Isogeometric analysis and error estimates for high order partial differential equations in fluid dynamics[END_REF] we have the following inequality :

|ψ -(Π V h (ψ) -c)| H 1 (Ωc) = |ψ -Π V h (ψ)| H 1 (Ωc) ≤ Ch l+1 ||ψ|| H l+2 (Ωc)
which is equivalent to saying that we have :

(4.20) Ωc (φ n+1 -φ n+1 h -c n+1 ) 2 dx ≤ Ch l+1 ||ψ|| H l+2 (Ωc) |φ n+1 -φ n+1 h | H 1 (Ωc)
as a consequence, replacing the results in inequality 4.11, we obtain

|φ n+1 -φ n+1 h | 2 H 1 (Ωc) ≤ Ch l-1 |φ n+1 -φ n+1 h | H 1 (Ωc) + Ch 2l-2 + Ch 3(l-1) 2 |φ n+1 -φ n+1 h | 1 2 H 1 (Ωc) , (4.21)
with Young's inequality (ab ≤ a 2 + b 2 4 ) in the first and (ab ≤ 3a

4 3 4 + b 4 4
) in the last term of 4.21, we have

(1 - 1 2 )|φ n+1 -φ n+1 h | 2 H 1 (Ωc) ≤ Ch 2l-2 + Ch 2l-2 + Ch 2l-2 , (4.22)
finally, the estimation is given by 1) , where C constant depend on ||φ n+1 || H l (Ωc) .

|φ n+1 -φ n+1 h | 2 H 1 (Ωc) ≤ Ch 2(l-
Corollary 4.1. The error estimation for Monge-Ampere equation in H 2 -norm fulfills

(4.23) |φ n+1 -φ n+1 h | H 2 (Ωc) ≤ Ch l-2 .
Proof. By taking the Projector Π

V h : L 2 (Ω c ) -→ V h , we get |φ n+1 -φ n+1 h | H 2 (Ωc) ≤ |φ n+1 -φ n+1 h -Π V h (φ n+1 -φ n+1 h )| H 2 (Ωc) + |Π V h (φ n+1 -φ n+1 h )| H 2 (Ωc)
by the proposition 3.2 [52]

|φ n+1 -φ n+1 h -Π V h (φ n+1 -φ n+1 h )| H 2 (Ωc) = |φ n+1 -Π V h (φ n+1 )| H 2 (Ωc) ≤ Ch l-2 ||φ n+1 || H l (Ωc) .
By the usual inverse inequality as in 4.10, we have

|Π V h (φ n+1 -φ n+1 h )| H 2 (Ωc) ≤ C θ h -1 |Π V h (φ n+1 -φ n+1 h )| H 1 (Ωc) ≤ C θ h -1 |φ n+1 -φ n+1 h | H 1 (Ωc) .
Therefore, according to the last proposition our result is achieved.

Remark 4.2. If we denote by K T the constant of the contraction T , such that K T < 1, we have:

(4.24) |φ -φ n h | H 2 (Ωc) ≤ C( K T n + h l-2 ).
4.3. Fast Diagonalization method. In order to devise a fast solver the Poisson and Laplace problems, we choose to follow the work of Sangalli and Tani [START_REF] Sangalli | Isogeometric preconditioners based on fast solvers for the sylvester equation[END_REF], we describe in the sequel the fast diagonalization method in the case of Isogeometric Analysis. This method was first introduced in [START_REF] Lynch | Direct solution of partial difference equations by tensor product methods[END_REF].

For the seek of simplicity, we shall consider the following Laplace problem,

(4.25) -∇ 2 u + τ u = f, Ω u = 0, ∂Ω
The Poisson problem and its solver shall be retrieved with τ = 0. After discretizing the Laplace problem, we get the following linear system (4.26)

L τ x := (K 1 ⊗ M 2 ⊗ M 3 + M 1 ⊗ K 2 ⊗ M 3 + M 1 ⊗ M 2 ⊗ K 3 + τ M 1 ⊗ M 2 ⊗ M 3 ) x = b
where M and K are unidimensional parametric mass and stiffness matrices respectively, ⊗ is a Kronecker product.

We first consider the generalized eigendecompositions problems (4.27)

K 1 U 1 = M 1 U 1 D 1 , K 2 U 2 = M 2 U 2 D 2 , K 3 U 3 = M 3 U 3 D 3 ,
where D 1 , D 2 and D 3 are diagonal matrices such that U 1 , U 2 and U 3 fulfills (4.28)

U T 1 M 1 U 1 = I 1 , U T 2 M 2 U 2 = I 2 , U T 3 M 3 U 3 = I 3 Therefor, (4.26) can be written as (4.29) (U 1 ⊗ U 2 ⊗ U 3 ) -1 (D 1 ⊗ I 2 ⊗ I 3 + I 1 ⊗ D 2 ⊗ I 3 + I 1 ⊗ I 2 ⊗ D 3 + τ I 1 ⊗ I 2 ⊗ I 3 ) (U 1 ⊗ U 2 ⊗ U 3 ) -T x = b
The direct solver for the Laplace problem (4.26), is then given by the following algorithm, where we omit the initialization step achieved by solving the generalized eigendecompositions in (4.27):

Algorithm 1: fast diag: Fast diagonalization method for Laplace problem

Input : L τ , b Output: x 1 b ← (U 1 ⊗ U 2 ⊗ U 3 ) b 2 x ← (D 1 ⊗ I 2 ⊗ I 3 + I 1 ⊗ D 2 ⊗ I 3 + I 1 ⊗ I 2 ⊗ D 3 + τ I 1 ⊗ I 2 ⊗ I 3 ) -1 b 3 x ← (U 1 ⊗ U 2 ⊗ U 3 ) T x 5. Numerical Results
In this section, we show the performance and a limit of standard method by • presenting numerical convergence and CPU timing of our Monge-Ampère solver in 2D and 3D,

• showing examples of two-dimensional adaptive mesh generation,

• showing the limitation of using the monitor function with a variation near the boundary. Our numerical code is written in Python, where intensive computations parts were accelerated using Pyccel [START_REF] Ratnani | python static compiler for high performance computing[END_REF], a static compiler for Python, based on Fortran and C. A laptop with a quad 2.6 GHz Intel Core i7, on ubuntu OS is used throughout the numerical tests. 5.1. Monge-Ampère solver using BFO method. The goal of this subsection is to validate our solver 4.1 by a manufactured convex solutions of Monge-Ampere equation. The L 2 residual error tolerance in the Picard iteration is fixed at 10 -10 .

Example 5.1 (Two-dimensional case). We test our solver using a convex analytical solution in the unit square, defined as follows.

(5.1) φ(ξ, η) = exp ξ 2 + η 2 2 for which the source term is

(5.2) f (ξ, η) = (1 + ξ 2 + η 2 )exp(ξ 2 + η 2 ).
Table 1. CPU-time needed to solve the Monge-Ampere equation and error evolution for various polynomial degrees p along each direction. Example 5.2 (Three-dimensional case). We test our solver by a convex manufactured solution on the cube case, defined as :

(a) Degree p = 3 #cells H 1 -norm H 1 -
(5.3) φ(ξ, η, γ) = exp ξ 2 + η 2 + γ 2 3
where the source term is Remark 5.1. For numerical purpose d! is replaced by δ such that

(5.4) f (ξ, η, γ) = 8 81 (3 + 2(ξ 2 + η 2 + γ 2 )) exp(ξ 2 + η 2 + γ 2 ).
(5.5) G(φ) = (∆φ) d + δ(f (∇φ)) -det(H(φ)) 1 d
Then, we see that a choice of δ in three-dimension affects the number of iterations of the Picard algorithm to achieve the same precision Fig. 3. 5.2. Adaptive mesh generation. In the following examples, we show in different figures two-dimensional adaptive meshes obtained using bi-cubic B-splines with #cells = 32. Here, a tolerance of L 2 residual error in Picard's iteration is fixed at 10 -8 .

As in [18], we define the total error in the scheme and the quality of the adapted grid by measuring the global displacement of the grid points.

(5.6)

Err 2 = Ωc | det(H(φ))-f (∇(φ))| 2 dx Qual 2 = Ωc |∇(φ(x))-x| 2 / f (∇(φ(x))) dx.
Examples on the unit square.

Example 5.3. In this example, the monitor function represents a circle centered in the middle of the square Fig. 4. We display CPU-time, the min-max value of the Jacobian function, and the accuracy in the Table 3.

ρ(ξ, η) = 1 + 5 exp(-50 (ξ -0.5) 2 + (η -0.5) 2 -0.09 ), Table 3. Grid convergence analysis for various polynomial degrees p along each direction.

( Examples on more general geometries.

Example 5.5 (Circle). In this test we show the Mickey Mouse example, the mesh quality using a new technique of mesh adaptation can be seen clearly in the Figs.

7-8. However, to clear the idea about a limit using the standard method, we add the last colon in the Table 5 showing the error at the boundary. Results of adaptive meshes shown in Fig. 9, there is no mesh-distortion but there are some errors at the boundary Table 6 that vanish with refinement.

Table 6. Grid convergence analysis for various polynomial degrees p along each direction.

( 6. Mixed Finite Elements approach 6.1. Mixed formulation of Benamou-Froese-Oberman system. We can obtain a more accurate mapping with an exact normal boundary condition by using a mixed formulation that consists of taking into a count the boundary condition in strong form Ψ(∂Ω c ) = ∂Ω c . Hence, returning to the system 4.2, by introducing Ψ = ∇φ, we have

   ∇ • Ψ = G(Ψ) in Ω c , Ψ = ∇φ in Ω c , Ψ • n = x • n on ∂Ω c , (6.1)
where

(6.2) G(Ψ) = (∇.Ψ) d + d!(f (Ψ) -det ∇(Ψ) 1 d .
Hence, let's introduce the decomposition

Ψ = Ψ 0 + Ψ N such that Ψ 0 • n = 0 and Ψ N • n = x • n.
Therefor, the problem (6.1) writes

   ∇ • Ψ 0 = G(Ψ) -∇ • Ψ N in Ω c , Ψ 0 = ∇φ -Ψ N in Ω c , Ψ 0 • n = 0 on ∂Ω c . (6.3)
Then, the variational formulation of 6.3, is :

Find (Ψ, φ) ∈ H 0 (div, Ω) × L 2 0 (Ω) such that ∀(v, ϕ) ∈ H 0 (div, Ω) × L 2 0 (Ω) Ω Ψ 0 • v + Ω φ ∇ • v = -Ω Ψ N • v, Ω (∇ • Ψ 0 ) ϕ = Ω G(Ψ)ϕ -Ω (∇ • Ψ N ) ϕ, (6.4
) 6.2. Fast Diagonalization method for Mixed Variational Formulation. Following the same steps as in section 4.3, we derive a new algorithm of fast diagonalization in two and three dimensions.

We consider for simplicity Poisson equation with homogeneous Neumann boundary conditions.

-∇ 2 φ = f in Ω, ∇φ • n = 0 on ∂Ω, (6.5) As well known, the associated variational formulation can be written as a minimization problem over H 1 (Ω). Let us now, introduce an auxiliary variable u = ∇φ. The Poisson problem is then equivalent to the following system

   -∇ • u = f in Ω, u = ∇φ in Ω, u • n = 0 on ∂Ω. (6.6)
The mixed weak formulation for the Poisson equation is

Find (u, φ) ∈ H 0 (div, Ω) × L 2 0 (Ω) such that ∀(Ψ, ϕ) ∈ H 0 (div, Ω) × L 2 0 (Ω) Ω u • Ψ + Ω φ ∇ • Ψ = 0, Ω (∇ • u) ϕ = -Ω f ϕ. (6.7)
Using the discrete spaces defined in 2.13 and 2.14, we get the following variational problem

Find (u h , φ h ) ∈ V h (div, Ω) × V h (L 2 , Ω) such that ∀(Ψ, ϕ) ∈ V h (div, Ω) × V h (L 2 , Ω) Ω u h • Ψ + Ω φ h ∇ • Ψ = 0, Ω (∇ • u h ) ϕ = -Ω f ϕ. (6.8)
We then recover a saddle point problem of the form

A u h φ h = A B B T 0 u h φ h = 0 F (6.9)
with

F i = - Ω f φ i dx, 1 ≤ i ≤ N V h
where N V h is a number of B-spline basis functions.

Notations. Before expliciting the matrix forms in 2D and 3D, we start by introducing some 1D matrices that will be needed.

                 (M s ) is,js = 1 0 B ps µs, is B ps µs, js dx s (D s ) is,js = 1 0 B ps-1 µs-1, is B ps-1 µs-1, js dx s (R s ) is,js = 1 0 B ps-1 µs-1, is dB ps µs, js dξ dx s (6.10)
where

B ps µs, is = B ps is , such that B ps is ∈ S ps µs , s = 1, 2.., d.
Matrix form in 2D. For the sake of simplicity we shall introduce the auxiliary scalar functions

               Ψ 1 j = B p1 µ1, j1 B p2-1 µ2-1, j2 Ψ 2 j = B p1-1 µ1-1, j1 B p2 µ2, j2 φ j = B p1-1 µ1-1, j1 B p2-1 µ2-1, j2 (6.11)
where j = (j 1 , j 2 ), we also define the vectors e 1 = 1 0 and e 2 = 0 1 . Therefor, the expression of u h ∈ V h (div, Ω) becomes

u h = j u 1 j Ψ 1 j e 1 + u 2 j Ψ 2 j e 2
Because e i • e j = δ ij , we find that A is a block diagonal matrix of the form

A = A 11 0 0 A 22 (6.12)
where

A 11i,j = Ω Ψ 1 i Ψ 1 j dx = (M 1 ⊗ D 2 ) ij A 22i,j = Ω Ψ 2 i Ψ 2 j dx = (D 1 ⊗ M 2 ) ij
While for B we have

B = B 1 B 2 (6.13)
where

B 1i,j = Ω φ j ∂ x1 Ψ 1 i dx = (R 1 ⊗ D 2 ) ij B 2i,j = Ω φ j ∂ x2 Ψ 2 i dx = (D 1 ⊗ R 2 ) ij
The final matrix form is

A =   M 1 ⊗ D 2 0 R 1 ⊗ D 2 0 D 1 ⊗ M 2 D 1 ⊗ R 2 (R 1 ⊗ D 2 ) T (D 1 ⊗ R 2 ) T 0   Matrix form in 3D.
It is possible to extend the notation of the auxiliary scalar functions to the 3D case as follows:

Ψ 1 j = B p1 µ1, j1 B p2-1 µ2-1, j2 B p3-1 µ3-1, j3 Ψ 2 j = B p1-1 µ1-1, j1 B p2 µ2, j2 B p3-1 µ3-1, j3 Ψ 3 j = B p1-1 µ1-1, j1 B p2-1 µ2-1, j2 B p3 µ3, j3
where j = (j 1 , j 2 , j 3 ), and

φ j = B p1-1 µ1-1, j1 B p2-1 µ2-1, j2 B p3-1 µ3-1, j3
we also define the vectors

e 1 =   1 0 0   , e 2 =   0 1 0   and e 3 =   0 0 1   . The expression of u h ∈ V h (div, Ω) equivalent to u h = j u 1 j Ψ 1 j e 1 + u 2 j Ψ 2 j e 2 + u 3 j Ψ 3 j e 3
and as a consequence of e i • e j = δ ij , we have

A =   A 11 0 0 0 A 22 0 0 0 A 33   (6.14)
where

A 11i,j = Ω Ψ 1 i Ψ 1 j dx = (M 1 ⊗ D 2 ⊗ D 3 ) ij A 22i,j = Ω Ψ 2 i Ψ 2 j dx = (D 1 ⊗ M 2 ⊗ D 3 ) ij A 33i,j = Ω Ψ 3 i Ψ 3 j dx = (D 1 ⊗ D 2 ⊗ M 3 ) ij
While for B we have

B =   B 1 B 2 B 3   (6.15)
where

B 1i,j = Ω φ j ∂ x1 Ψ 1 i dx = (R 1 ⊗ D 2 ⊗ D 3 ) ij B 2i,j = Ω φ j ∂ x2 Ψ 2 i dx = (D 1 ⊗ R 2 ⊗ D 3 ) ij B 3i,j = Ω φ j ∂ x3 Ψ 3 i dx = (D 1 ⊗ D 2 ⊗ R 3 ) ij Finally, A =     M 1 ⊗ D 2 ⊗ D 3 0 0 R 1 ⊗ D 2 ⊗ D 3 0 D 1 ⊗ M 2 ⊗ D 3 0 D 1 ⊗ R 2 ⊗ D 3 0 0 D 1 ⊗ D 2 ⊗ M 3 D 1 ⊗ D 2 ⊗ R 3 (R 1 ⊗ D 2 ⊗ D 3 ) T (D 1 ⊗ R 2 ⊗ D 3 ) T (D 1 ⊗ D 2 ⊗ R 3 ) T 0     6.
2.1. Linear Solver. Based on the problem 6.9, we consider the following saddle point problem

A x 1 x 2 = A B B T 0 x 1 x 2 = b 1 b 2 (6.16)
Let us define à = B T A -1 B as the Schur complement of (6.16). Using the Uzawa method leads to

Ãx 2 = B T A -1 b 1 -b 2 Ax 1 = b 1 -Bx 2 (6.17)
Fast Diagonalization method (2D). The matrix à can be expressed as the sum of Kronecker matrices

à = D 1 ⊗ R T 2 M -1 2 R 2 + R T 1 M -1 1 R 1 ⊗ D 2 (6.18)
We have

B T A -1 b 1 = (R 1 ⊗ D 2 ) T (D 1 ⊗ R 2 ) T (M 1 ⊗ D 2 ) -1 0 0 (D 1 ⊗ M 2 ) -1 b 11 b 12 = R T 1 ⊗ D T 2 D T 1 ⊗ R T 2 M -1 1 ⊗ D -1 2 0 0 D -1 1 ⊗ M -1 2 b 11 b 12 = R T 1 ⊗ D T 2 M -1 1 ⊗ D -1 2 D T 1 ⊗ R T 2 D -1 1 ⊗ M -1 2 b 11 b 12 = R T 1 M -1 1 ⊗ I 2 I 1 ⊗ R T 2 M -1 2 b 11 b 12 = R T 1 M -1 1 ⊗ I 2 b 11 + I 1 ⊗ R T 2 M -1 2 b 12 Therefor Ãx 2 = b (6.19) where b = R T 1 M -1 1 ⊗ I 2 b 11 + I 1 ⊗ R T 2 M -1 2 b 12 -b 2
On the other hand, we have

Bx 2 = (R 1 ⊗ D 2 ) x 2 (D 1 ⊗ R 2 ) x 2 (6.20)
By introducing x 1 := (x 11 , x 12 )

T , the equation Ax 1 = b 1 -Bx 2 is equivalent to

x 11 = M -1 1 ⊗ D -1 2 b 11 -M -1 1 ⊗ D -1 2 (R 1 ⊗ D 2 ) x 2 x 12 = D -1 1 ⊗ M -1 2 b 12 -D -1 1 ⊗ M -1 2 (D 1 ⊗ R 2 ) x 2
which can be simplified as

x 11 = M -1 1 ⊗ D -1 2 b 11 -M -1 1 R 1 ⊗ I 2 x 2 x 12 = D -1 1 ⊗ M -1 2 b 12 -I 1 ⊗ M -1 2 R 2 x 2
In order to devise a fast solver (6.17), we choose to follow the approach of Sangalli and Tani [START_REF] Sangalli | Isogeometric preconditioners based on fast solvers for the sylvester equation[END_REF]. We describe in the sequel the fast diagonalization method in the case of Isogeometric Analysis. This method was first introduced in [START_REF] Lynch | Direct solution of partial difference equations by tensor product methods[END_REF]. We consider the generalized eigendecompositions problems (6.21)

R T 1 M -1 1 R 1 U 1 = D 1 U 1 d 1 , R T 2 M -1 2 R 2 U 2 = D 2 U 2 d 2
where d 1 and d 2 are diagonal matrices such that U 1 and U 2 fulfills (6.22)

U T 1 D 1 U 1 = I 1 , U T 2 D 2 U 2 = I 2 Therefor, (6.23) (U 1 ⊗ U 2 ) -T (d 1 ⊗ I 2 + I 1 ⊗ d 2 ) (U 1 ⊗ U 2 ) -1 x = b
The direct solver for (6.17), is then given by the following algorithm, where we omit the initialization step achieved by solving the generalized eigendecompositions in (6.21):

Algorithm 2: fast diag: Fast diagonalization method for (6.17) in the 2D case Input : b Output:

x 1 (b 11 , b 12 ) , b 2 ← unfold(b) 2 r 2 ← R T 1 M -1 1 ⊗ I 2 b 11 + I 1 ⊗ R T 2 M -1 2 b 12 -b 2 3 r2 ← (U 1 ⊗ U 2 ) T r 2 4 x2 ← (d 1 ⊗ I 2 + I 1 ⊗ d 2 ) -1 r2 5 x 2 ← (U 1 ⊗ U 2 ) x2 6 x 11 ← M -1 1 ⊗ D -1 2 b 11 -M -1 1 R 1 ⊗ I 2 x 2 7 x 12 ← D -1 1 ⊗ M -1 2 b 12 -I 1 ⊗ M -1 2 R 2 x 2 8 x ← fold(x 11 , x 12 , x 2 )
Fast Diagonalization method (3D). The matrix à in (3D) can be expressed as the sum of Kronecker matrices

à = D 1 ⊗ D 2 ⊗ R T 3 M -1 3 R 3 + D 1 ⊗ R T 2 M -1 2 R 2 ⊗ D 3 + R T 1 M -1 1 R 1 ⊗ D 2 ⊗ D 3 (6.24)
And as in (2D) case we can show that we have

B T A -1 b 1 = R T 1 M -1 1 ⊗ I 2 ⊗ I 3 b 11 + I 1 ⊗ R T 2 M -1 2 ⊗ I 3 b 12 + I 1 ⊗ I 2 ⊗ R T 3 M -1 3 b 13 .
Therefor

Ãx 2 = R T 1 M -1 1 ⊗ I 2 ⊗ I 3 b 11 + I 1 ⊗ R T 2 M -1 2 ⊗ I 3 b 12 + I 1 ⊗ I 2 ⊗ R T 3 M -1 3 b 13 -b 2 . (6.25)
On the other hand, we have

Bx 2 =   (R 1 ⊗ D 2 ⊗ D 3 ) x 2 (D 1 ⊗ R 2 ⊗ D 3 ) x 2 (D 1 ⊗ D 2 ⊗ R 3 ) x 2   (6.26) By introducing x 1 := (x 11 , x 12 , x 13 ) T , the equation Ax 1 = b 1 -Bx 2 is equivalent to      x 11 = M -1 1 ⊗ D -1 2 ⊗ D -1 3 b 11 -M -1 1 ⊗ D -1 2 ⊗ D -1 3 (R 1 ⊗ D 2 ⊗ D 3 ) x 2 x 12 = D -1 1 ⊗ M -1 2 ⊗ D -1 3 b 12 -D -1 1 ⊗ M -1 2 ⊗ D -1 3 (D 1 ⊗ R 2 ⊗ D 3 ) x 2 x 13 = D -1 1 ⊗ D -1 2 ⊗ M -1 3 b 13 -D -1 1 ⊗ D -1 2 ⊗ M -1 3 (D 1 ⊗ D 2 ⊗ R 3 ) x 2
which can be simplified to

     x 11 = M -1 1 ⊗ D -1 2 ⊗ D -1 3 b 11 -M -1 1 R 1 ⊗ I 2 ⊗ I 3 x 2 x 12 = D -1 1 ⊗ M -1 2 ⊗ D -1 3 b 12 -I 1 ⊗ M -1 2 R 2 ⊗ I 3 x 2 x 13 = D -1 1 ⊗ D -1 2 ⊗ M -1
Following the same steps, we consider the generalized eigendecompositions problems (6.27)

R T 1 M -1 1 R 1 U 1 = D 1 U 1 d 1 , R T 2 M -1 2 R 2 U 2 = D 2 U 2 d 2 , R T 3 M -1 3 R 3 U 3 = D 3 U 3 d 3 where d 1 , d 2 ,
and d 3 are diagonal matrices such that U 1 , U 2 and U 2 fulfills (6.28)

U T 1 D 1 U 1 = I 1 , U T 2 D 2 U 2 = I 2 , U T 3 D 3 U 3 = I 3 Therefor (6.29) (U 1 ⊗ U 2 ⊗ U 3 ) -T (d 1 ⊗ I 2 ⊗ I 3 + I 1 ⊗ d 2 ⊗ I 3 + I 1 ⊗ I 2 ⊗ d 3 ) (U 1 ⊗ U 2 ⊗ U 3 ) -1 x = b
The direct solver for (6.17), is then given by the following algorithm:

Algorithm 3: fast diag: Fast diagonalization method for (6.17) in the 3D case

Input : b Output: x 1 (b 11 , b 12 , b 13 ) , b 2 ← unfold(b) 2 r 2 ← R T 1 M -1 1 ⊗ I 2 ⊗ I 3 b 11 + I 1 ⊗ R T 2 M -1 2 ⊗ I 3 b 12 + I 1 ⊗ I 2 ⊗ R T 3 M -1 3 b 13 -b 2 3 r2 ← (U 1 ⊗ U 2 ⊗ U 3 ) T r 2 4 x2 ← (d 1 ⊗ I 2 ⊗ I 3 + I 1 ⊗ d 2 ⊗ I 3 + I 1 ⊗ I 2 ⊗ d 3 ) -1 r2 5 x 2 ← (U 1 ⊗ U 2 ⊗ U 3 ) x2 6 x 11 ← M -1 1 ⊗ D -1 2 ⊗ D -1 3 b 11 -M -1 1 R 1 ⊗ I 2 ⊗ I 3 x 2 7 x 12 ← D -1 1 ⊗ M -1 2 ⊗ D -1 3 b 12 -I 1 ⊗ M -1 2 R 2 ⊗ I 3 x 2 8 x 13 ← D -1 1 ⊗ D -1 2 ⊗ M -1 3 b 13 -I 1 ⊗ I 2 ⊗ M -1 3 R 3 x 2 9
x ← fold(x 11 , x 12 , x 13 , x 2 ) Remark 6.1. For the steps 6-7 in Algorithm 2 and steps 6-7-8 in Algorithm 3, we use Algorithm 1 instead of inverse matrices.

Numerical Results

To assess the efficiency of our method, the present section

• shows numerical convergence and CPU timing of our Monge-Ampère solver in 2D and 3D, • compares CPU timing in 3D for the standard solver 4.1 and Mixed variational formulation solver 6.2, • shows that we have solved the problem of errors at the boundary,

• shows examples of adaptive mesh generation using Mixed variational formulation in 2D, 2.5D and 3D.

7.1. Monge-Ampère solver using mixed variational formulation. We adopt the same test in both two and three dimension as in section 5.1. We also keep the L 2 residual error tolerance in the Picard iteration which is fixed at 10 -10 .

Example 7.1 (Two-dimensional case). In the following Table 7 we remark that we obtain more accuracy in the calculation using the mixed formulation compared to the standard method Table 1.

Table 7. CPU-time needed to solve the Monge-Ampere equation and error evolution for various polynomial degrees p along each direction. Example 7.2 (Three-dimensional case). The table 8 shows high accuracy and fast convergence in terms of the number of iterations performed using a multi-grid initial estimate. However, we show in the figure 10 CPU-time behavior for each method where we fix a tolerance of H 1 -norm at 10 -5 . Table 8. CPU-time needed to solve the Monge-Ampere equation and error evolution for various polynomial degrees p along each direction.. Figure 10. CPU time profiles between standard method (H 2 ) and mixed variational formulation (Hdiv).

(a) Degree p = 3 #cells H 1 -norm H 1 -
(a) Degree p = 3 #cells H 1 -norm H 1 -

Adaptive mesh generation.

Examples on the unit square. As it is for the standard method, in order to compare the result, we keep the same tolerance and #cell = 32 for all figures in the following examples.

Example 7.3. We test with the same monitor function as in [18] defined as :

(7.1) ρ(ξ, η) = 1 2 + cos(8π (ξ -0.5) 2 + (η -0.5) 2 ) .
Fig. 11 shows the adaptive meshes results and Table 9 shows a high quality of mesh adaptation.

Table 9. Grid convergence analysis for various polynomial degrees p along each direction.

(a) Degree p = 3 Example 7.4. In order to show that we have solved the problem of boundary errors and negative values, we adopt the same example 5.4 used in the standard method. As expected, the errors at the boundary are eliminated Fig. 12 with high accuracy validated by Table 10.

#
Table 10. Grid convergence analysis for various polynomial degrees p along each direction.

( 12 shows fast convergence, high accuracy and positive Jacobian function. Examples of two and half-dimensional adaptive meshes. Solving the Monge-Ampere equation directly in the torus, Möbius, and other 2.5-dimensional surfaces is complicated because of convexity and other assumptions necessary to ensure convergence. However, in what follows, we give a starting point for the moving mesh method into this complicated geometry.

Example 7.8 (Empty-Torus subdomain). In this example, we present a mesh adaptation in an extracted part of the torus domain. The geometry is represented in a three-dimensional space as shown in Fig. 18, but thinks to our new technique 1, the mesh adaptation is done in the square with a high quality of mesh adaption Table 14. Where the monitor function is defined in the torus by : (7.4) ρ(ξ, η, γ) = 1 2 + cos(8π (x -0.5) 2 + (y -0.5) 2 + (z -0.5) 2 ) . Example 7.9 (Non-Orientable Möbius Surface). From the application viewpoint Möbius strip is well-known as a tool for building more efficient mechanical belts. Hence, we defend a mesh adaptation in a non-orientable Möbius surface represented in a three-dimensional space Fig. 7.9. Table 15 shows that accurate and bijective mapping is always achieved using various regularities.. Therefore, here we use the exact form of non-axisymmetric function used in [START_REF] Andrew Tt Mcrae | Optimal-transport-based mesh adaptivity on the plane and sphere using finite elements[END_REF] into the Möbius strip as follows.

ρ(x, y) = 1 + 5/sech 5 (x -√

2 ) 2 + (y -0.5) 2 + (z -0.5) 2 - π 2 4 (7.5) + 5/sech 5 (x + √ 3 2 ) 2 + (y -0.5) 2 + (z -0.5) 2 - π 2 4 . 3
Table 15. Grid convergence analysis for various polynomial degrees p along each direction.

( Examples of three-dimensional adaptive meshes. In this subsection, we shall give some examples of three-dimensional adaptive meshes using different analytical density functions defined in the physical domain. However, To show the quality of meshes, we use the last formula 5.6. The tolerance is set to 10 -8 throughout the coming tests.

Example 7.10. We use the generalized form of the monitor function in 7.3 defined as :

(7.6) ρ(ξ, η, γ) = 1 2 + cos(8π (ξ -0.5) 2 + (η -0.5) 2 + (γ -0.5) 2 ) .
We can see a fast convergence and the accuracy of B-spline mapping in Table 16. Fig. 20 shows high quality of mesh adaptation. Example 7.11 (Examples of adapted mesh on more general geometries: sphere).

We assess the quality of mesh adaptation in the sphere case, see Fig. 21 and no mesh-distortion with high quality as it is demonstrated in Tab. 17 using a density function as cos-waves defined by ρ(x, y) = 3/ 2 + cos 8π (x -0.5) 2 + (y -0.5) 2 + (z -0.5) 2 . (7.7) Example 7.12 (Examples of adapted mesh on more general geometries: cylinder).

Using the same density function as in Example 7.7, we show numerical results for the case of the cylinder that has four singular points on each side of the boundary. The adaptive meshes in Fig. 22 (left) and a close-up Fig. 22 (right) illustrates the quality of mesh adaptation. However, even if the error does not change for different degrees, the regularity may be mandatory requirement by the end user, so we can guarantee that bijectivity is preserved with a high-quality of mesh adaptation, see Tab. 18

Table 18. Grid convergence analysis for various polynomial degrees p along each direction. Example 7.13 (Examples of adapted mesh on more general geometries: halftorus). A torus domain is one of the popular geometries used in many papers for the three-dimensional simulation of plasma physic problems. Therefore, we study this geometry by adding mesh adaptation. However, we use the same density function 7.7 presented in Fig. 23

Conclusions

We have presented the original technique of the adaptive meshes, and it works efficiently. The L 2 Monge-Kantorovich problem gave us a good strategy of mesh adaptation in the Isogeometric Analysis; solving the Monge-Ampere equation by the BFO method that is robust and reliable as shown in section 6, our method constructs a global and regular one-to-one mapping that maps the unit square into the physical domain, while ensuring the equidistribution property. Therefore, the bijectivity from the unit square to the physical domain is assured when we use an initial bijective mapping.

We started our study with a standard method which solves the BFO system directly by means of the Galerkin variational formulation, a fast convergence is always obtained using the fast diagonalization method Algorithm 1. However, the accuracy is provided for some examples, but unfortunately, it generates boundary errors due to the weak integration of the Neumann boundary conditions, which may lead to the loss of the geometry. We fixed this problem by introducing a mixed variational formulation, and the boundary errors are fully resolved. Moreover, for fast convergence, two new fast Alg. 2 and Alg. 3 based on the fast diagonalization method are proposed, and we obtain a more accurate bijective mapping, as shown in the tables for various B-spline regularities and geometric complexities in twoand three-dimensional cases.

The performance and stability of the technique of adaptation are demonstrated by several examples. The algorithm is simple to implement, accelerated by the fast diagonalization method, and has no free parameters. Moreover, it works for all geometry that we have tested.

Current limitations and future work.

Figure 1 .

 1 Figure 1. Illustration of the Adaptive Mapping idea. The old mapping F : ξ → x and the new mapping F : ξ → x map the logical domain (patch) onto the same physical domain Ωp. The mapping Ψ : ξ → ξ , maps the old (uniform) mesh onto the new one in parametric domain.

Figure 2 .

 2 Figure 2. (left) A cubic B-Spline curve and its control points using the knot vector Ξ = {0000 1 4 1 2 1 2 3 4 3 4 3 4 1111}, (right) the corresponding B-Splines.

Figure 3 .

 3 Figure 3. L 2 and H 1 -residual and H 1 -error profiles with different choices of δ = 6 (left), δ = 9 (middle) and δ = 27 (right) using bi-cubic B-splines.

Figure 4 .Example 5 . 4 .

 454 Figure 4. (left) Adaptive meshes and (right) the monitor function.

Figure 5 .

 5 Figure 5. (left) Adaptive meshes and (right) the monitor function.

Figure 6 .

 6 Figure 6. A close-up of the region near the boundary.

ρ

  (x, y) = 1. + 5 exp(-100|(x -0.45) 2 + (y -0.4) 2 -0.1|) + 5 exp(-100|x 2 + y 2 -0.2|) (5.8) + 5 exp(-100|(x + 0.45) 2 + (y -0.4) 2 -0.1|).

Figure 7 .

 7 Figure 7. (left) Adaptive meshes and (right) the monitor function.

Example 5 . 6 ( 2 +

 562 Circle). In this example we use a non-axisymmetric monitor function defined by: ρ(x, y) = 1 + 5/sech 5 (x -(y -0.5) 2 -π 2 4 .

Figure 8 .

 8 Figure 8. A close-up of the region where the two circles intersect (left) and on one of the singular points (right).

Figure 9 .

 9 Figure 9. (left) Adaptive meshes and (right) the monitor function.

Figure 11 .

 11 Figure 11. (left) Adaptive meshes and (right) the monitor function.

Figure 13 .

 13 Figure 13. (left) Adaptive meshes and (right) the monitor function.

Figure 14 .

 14 Figure 14. A close-up of the region where the two circles intersect (left) and on one of the singular points (right).

Example 7 . 6 (

 76 Circle). With no-errors at the boundary we adopt the example of non-axisymmetric monitor function used in 5.6 Fig. 15. On other hand, Table

Figure 15 .

 15 Figure 15. (left) Adaptive meshes and (right) the monitor function.

Example 7 . 7 (

 77 Examples of adapted mesh on more general geometries (quart-annulus)). In this test, we provide the quality of mesh adaptation in the quart-annulus domain as shown in the Tab. 13. However, there is no distortion in the mesh using a density function as two sin-waves, illustrated in Figs.16-17. ρ(x, y) = 1 + 5 exp -50|y -0.5 -0.25 sin(2πx) sin( -50|y -0.25 -0.25 sin(2πx)| .

Figure 16 .

 16 Figure 16. (left) Adaptive meshes and (right) the monitor function.

Figure 17 .

 17 Figure 17. A close-up of the region where the mesh is very fine.

Figure 18 .

 18 Figure 18. Adaptive meshes and the monitor function.

Figure 19 .

 19 Figure 19. Adaptive meshes and the monitor function.

Figure 20 .

 20 Figure 20. Location of the planes in the cube corresponding to x = 0.1, 0.6, 0.9.

Figure 21 .

 21 Figure 21. (left) A close-up of one extract part inside the sphere (right) Adaptive meshes and the monitor function.

Figure 22 .

 22 Figure 22. (left) A close-up of one extract part inside the sphere (right) Adaptive meshes and the monitor function.

Figure 23 .

 23 Figure 23. (left) A close-up of one extract part inside the sphere (right) Adaptive meshes and the monitor function.

Table 2 .

 2 CPU-time needed to solve the Monge-Ampere equation and error evolution for various polynomial degrees p along each direction.

			(a) Degree p = 3
	#cells H 1 -norm H 1 -norm MG CPU-time (s) CPU-time (s) MG Nbr-iter Nbr-iter MG
	8 5.885e-05	5.885e-05	0.214	0.187	14	10
	16 7.672e-06	7.672e-06	1.290	0.987	14	10
	32 9.793e-07	9.793e-07	9.671	6.683	14	9
	64 1.235e-07	1.236e-07	76.92	40.88	14	7
			(b) Degree p = 4
	#cells H 1 -norm H 1 -norm MG CPU-time (s) CPU-time (s) MG Nbr-iter Nbr-iter MG
	8 1.758e-06	1.758e-06	0.528	0.472	14	10
	16 1.184e-07	1.184e-07	3.547	2.680	14	10
	32 7.658e-09	7.658e-09	27.35	16.42	14	8
	64 9.609e-10	9.617e-10	216.1	98.91	14	6
	8 9.064e-08	9.064e-08	1.206	1.181	14	10
	16 3.266e-09	3.266e-09	9.118	6.787	14	10
	32 8.756e-10	8.757e-10	71.63	42.58	14	8
	64 8.561e-10	8.563e-10	582.9	270.6	14	6

(c) Degree p = 5 #cells H 1 -norm H 1 -norm MG CPU-time (s) CPU-time (s) MG Nbr-iter Nbr-iter MG

Table 4 .

 4 Grid convergence analysis for various polynomial degrees p along each direction.

			(a) Degree p = 3
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 4.816e-01	0.127	1.512e-01	-2.286e-02	3.587e+00
	32 1.892e-01	0.303	1.201e-01	-7.120e-02	3.644e+00
	64 4.322e-02	0.784	1.083e-01	1.936e-01	3.382e+00
	128 9.048e-03	2.259	1.081e-01	3.464e-01	3.128e+00
			(b) Degree p = 4
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 4.702e-01	0.190	1.414e-01	5.377e-02	4.142e+00
	32 1.288e-01	0.461	1.103e-01	4.119e-02	3.384e+00
	64 2.965e-02	0.978	1.081e-01	2.215e-01	3.300e+00
	128 4.821e-03	3.194	1.081e-01	3.511e-01	3.136e+00
			(c) Degree p = 5
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 4.640e-01	0.320	1.459e-01	-5.108e-02	5.528e+00
	32 1.267e-01	0.549	1.099e-01	6.682e-02	3.461e+00
	64 2.787e-02	1.355	1.081e-01	2.348e-01	3.302e+00
	128 4.112e-03	4.430	1.081e-01	3.510e-01	3.140e+00

Table 5 .

 5 Grid convergence analysis for various polynomial degrees p along each direction.

				(a) Degree p = 3
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ) max ∂Ωc	|(Ψ -x) • n|
	16 2.720e-01	0.070	5.179e-02	9.727e-02	1.719e+00	2.426e-03
	32 7.321e-02	0.146	5.786e-02	9.899e-02	1.412e+00	3.969e-04
	64 1.765e-02	0.414	5.807e-02	1.338e-01	1.314e+00	4.189e-05
	128 5.066e-03	1.516	5.811e-02	1.365e-01	1.280e+00	2.881e-07
				(b) Degree p = 4
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ) max ∂Ωc	|(Ψ -x) • n|
	16 2.353e-01	0.090	5.217e-02	1.835e-01	1.613e+00	1.848e-03
	32 5.832e-02	0.222	5.818e-02	1.554e-01	1.429e+00	4.023e-04
	64 1.375e-02	0.617	5.807e-02	1.373e-01	1.307e+00	4.301e-05
	128 3.558e-03	2.158	5.811e-02	1.353e-01	1.280e+00	3.245e-07
	16 2.205e-01	0.123	5.482e-02	1.490e-01	1.621e+00	2.606e-03
	32 4.607e-02	0.292	5.812e-02	1.487e-01	1.367e+00	2.842e-04
	64 1.273e-02	0.858	5.807e-02	1.414e-01	1.306e+00	3.999e-05
	128 3.496e-03	3.093	5.811e-02	1.354e-01	1.280e+00	5.209e-07

(c) Degree p = 5 #cells Err CPU-time (s) Qual min Jac(Ψ) max Jac(Ψ) max ∂Ωc |(Ψ -x) • n|

Table 11 .

 11 Figure 12. (left) Adaptive meshes and (right) the monitor function.Example 7.5 (Examples of adapted mesh on more general geometries (circle)). We add to the example 5.5 another circle touching the boundary like an unhappy mouth for Mickey Mouse. A good performance in the computation is shown in table11, with no self-intersection Figs.13-14. Grid convergence analysis for various polynomial degrees

			a) Degree p = 3
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 4.165e-01	0.196	1.428e-01	1.077e-01	2.950e+00
	32 1.231e-01	0.600	1.111e-01	3.123e-02	3.413e+00
	64 3.020e-02	1.279	1.082e-01	2.209e-01	3.323e+00
	128 4.589e-03	3.782	1.081e-01	3.509e-01	3.132e+00
			(b) Degree p = 4
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 4.249e-01	0.216	1.435e-01	-3.402e-03	2.963e+00
	32 1.182e-01	0.620	1.101e-01	8.919e-02	3.567e+00
	64 2.747e-02	1.539	1.081e-01	2.483e-01	3.324e+00
	128 3.411e-03	4.250	1.081e-01	3.513e-01	3.134e+00
			(c) Degree p = 5
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 3.764e-01	0.344	1.407e-01	-6.914e-02	3.361e+00
	32 1.065e-01	0.764	1.092e-01	8.918e-02	3.269e+00
	64 2.213e-02	1.832	1.081e-01	2.296e-01	3.308e+00
	128 3.138e-03	5.358	1.081e-01	3.511e-01	3.137e+00

Table 12 .

 12 Grid convergence analysis for various polynomial degrees

	p along each direction.
			(a) Degree p = 3
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 5.900e-02	0.114	9.608e-02	1.665e-01	1.982e+00
	32 1.096e-02	0.226	9.611e-02	1.702e-01	1.893e+00
	64 2.048e-03	0.671	9.611e-02	1.699e-01	1.867e+00
	128 1.819e-03	2.490	9.611e-02	1.699e-01	1.867e+00
			(b) Degree p = 4
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 4.695e-02	0.138	9.629e-02	1.676e-01	1.955e+00
	32 8.832e-03	0.265	9.611e-02	1.703e-01	1.882e+00
	64 1.922e-03	0.774	9.611e-02	1.699e-01	1.867e+00
	128 1.816e-03	2.819	9.611e-02	1.699e-01	1.867e+00
			(c) Degree p = 5
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 4.790e-02	0.204	9.617e-02	1.662e-01	1.942e+00
	32 9.055e-03	0.366	9.611e-02	1.696e-01	1.882e+00
	64 1.867e-03	1.057	9.611e-02	1.699e-01	1.867e+00
	128 1.815e-03	3.612	9.611e-02	1.699e-01	1.867e+00

Table 13 .

 13 Grid convergence analysis for various polynomial degrees

	p along each direction.
			(a) Degree p = 3
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 6.810e-02	0.117	7.786e-02	1.880e-01	1.795e+00
	32 1.765e-02	0.230	7.752e-02	2.531e-01	1.678e+00
	64 4.754e-03	0.682	7.754e-02	2.740e-01	1.614e+00
	128 1.763e-03	2.397	7.754e-02	2.715e-01	1.614e+00
			(b) Degree p = 4
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 5.320e-02	0.145	7.794e-02	1.696e-01	1.702e+00
	32 1.721e-02	0.278	7.751e-02	2.570e-01	1.660e+00
	64 4.538e-03	0.821	7.754e-02	2.739e-01	1.616e+00
	128 1.707e-03	2.806	7.754e-02	2.718e-01	1.614e+00
			(c) Degree p = 5
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 5.087e-02	0.174	7.779e-02	1.737e-01	1.784e+00
	32 1.612e-02	0.367	7.751e-02	2.660e-01	1.670e+00
	64 4.386e-03	1.027	7.754e-02	2.739e-01	1.614e+00
	128 1.678e-03	3.737	7.754e-02	2.716e-01	1.614e+00

Table 14 .

 14 Grid convergence analysis for various polynomial degrees p along each direction.

			(a) Degree p = 3
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 3.626e-01	0.104	7.822e-03	3.264e-01	1.652e+00
	32 2.698e-01	0.275	8.649e-03	3.806e-01	1.702e+00
	64 7.026e-02	0.845	6.952e-03	4.616e-01	1.780e+00
	128 8.173e-03	2.870	6.952e-03	5.763e-01	1.759e+00
			(b) Degree p = 4
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 3.560e-01	0.133	8.141e-03	3.137e-01	1.975e+00
	32 2.700e-01	0.423	8.663e-03	3.650e-01	1.786e+00
	64 7.102e-02	1.160	6.953e-03	4.550e-01	1.755e+00
	128 6.492e-03	3.619	6.952e-03	5.712e-01	1.749e+00
			(c) Degree p = 5
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 3.599e-01	0.159	7.748e-03	2.696e-01	1.733e+00
	32 2.713e-01	0.534	8.511e-03	3.011e-01	1.791e+00
	64 7.167e-02	1.480	6.953e-03	4.484e-01	1.754e+00
	128 5.981e-03	4.582	6.952e-03	5.696e-01	1.742e+00

Table 16 .

 16 Grid convergence analysis for various polynomial degrees

	p along each direction.
			(a) Degree p = 3
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 5.303e-03	1.503	4.925e-03	4.922e-01	1.791e+00
	32 3.662e-04	8.520	3.329e-03	5.718e-01	1.746e+00
	64 2.451e-05	55.52	2.322e-03	5.733e-01	1.725e+00
			(b) Degree p = 4
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 3.763e-03	3.769	4.922e-03	5.166e-01	1.779e+00
	32 1.934e-04	21.40	3.330e-03	5.720e-01	1.732e+00
	64 4.271e-06	129.6	2.322e-03	5.736e-01	1.721e+00
			(c) Degree p = 5
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 2.780e-03	9.072	4.917e-03	5.229e-01	1.763e+00
	32 1.252e-04	45.97	3.330e-03	5.720e-01	1.725e+00
	64 3.368e-06	270.2	2.322e-03	5.736e-01	1.721e+00

Table 17 .

 17 Grid convergence analysis for various polynomial degrees p along each direction.

			(a) Degree p = 3
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 2.148e-02	1.716	3.108e-03	3.825e-01	1.920e+00
	32 2.005e-03	11.26	2.198e-03	5.141e-01	1.796e+00
	64 1.195e-04	74.83	1.533e-03	5.725e-01	1.748e+00
			(b) Degree p = 4
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 2.030e-02	4.208	3.116e-03	3.668e-01	2.139e+00
	32 1.580e-03	29.18	2.198e-03	5.167e-01	1.786e+00
	64 4.936e-05	166.7	1.533e-03	5.731e-01	1.732e+00

Table 19 .

 19 (right) and to see what happens inside the domain, see Fg.23 (left), and for the different degrees we validate the bijectivity and quality with fast convergence results in Table19. Grid convergence analysis for various polynomial degrees p along each direction.

			(a) Degree p = 3
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 4.801e-02	2.637	2.029e-03	3.284e-01	1.893e+00
	32 7.936e-03	12.84	1.338e-03	4.300e-01	1.829e+00
	64 4.406e-04	76.00	9.252e-04	5.721e-01	1.774e+00
			(b) Degree p = 4
	#cells	Err	CPU-time (s)	Qual	min Jac(Ψ) max Jac(Ψ)
	16 4.680e-02	7.165	2.095e-03	2.872e-01	2.030e+00
	32 7.310e-03	34.75	1.339e-03	4.481e-01	1.841e+00
	64 2.774e-04	214.9	9.252e-04	5.705e-01	1.755e+00
	16 4.321e-02	15.52	2.097e-03	2.005e-01	2.032e+00
	32 6.685e-03	76.55	1.338e-03	4.481e-01	1.799e+00
	64 2.335e-04	454.7	9.252e-04	5.699e-01	1.743e+00

(c) Degree p = 5 #cells Err CPU-time (s) Qual min Jac(Ψ) max Jac(Ψ)

b 13 -I 1 ⊗ I 2 ⊗ M -1

Acknowledgments. We warmly thank G. Latu, G. Dif-Pradalier, V. Grandgirard (CEA-Cadarache), B. Nkonga (University of Nice), G. Huysmans (Iter Organization), B. Mourain (INRIA Sophi-Antipolis) and A. Buffa (EPFL), for useful discussions and advices. This project started during the ANEMOS-ANR and EFDAfellowship and IPAM long program Computational Methods in High Energy Density Plasmas in 2012.