# Adaptive isogeometric analysis using optimal transport and their fast solvers

Mustapha Bahari, Abderrahmane Habbal, Ahmed Ratnani, Eric

Sonnendrücker

► To cite this version:

Mustapha Bahari, Abderrahmane Habbal, Ahmed Ratnani, Eric Sonnendrücker. Adaptive isogeometric analysis using optimal transport and their fast solvers. Computer Methods in Applied Mechanics and Engineering, 2024, 418, pp.116570. 10.1016/j.cma.2023.116570. hal-03936627

HAL Id: hal-03936627

https://inria.hal.science/hal-03936627

Submitted on 12 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ADAPTIVE ISOGEOMETRIC ANALYSIS USING OPTIMAL TRANSPORT

MUSTAPHA BAHARI, ABDERRAHMANE HABBAL, AHMED RATNANI, AND ERIC SONNENDRÜCKER

## ABSTRACT.

## CONTENTS

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

Received by the editor 26 October 2021.

1991 Mathematics Subject Classification. Primary.

Key words and phrases. B-spline discretization, Optimal Transport, Isogeometric Analysis, multigrid, mesh generation, r-refinement.

©XXXX American Mathematical Society
# 1. 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 as those encountered in Computational Plasmas Physics. Grid Adaptation and moving meshes are used in different fields and the literature is quite rich [10, 26, 30, 32]. They are also used to generate Anisotropic meshes [4, 41, 25, 26, 24].

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 introduction of the IsoGeometric Analysis paradigm by Hughes [37].

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 [46]; thus, based on the optimal transport problem the process requires the resolution of the Monge-Ampère equation on the sphere domain (earth shape), the investigation in this direction is given in the papers [53, 46, 21]. Moreover, the global adaptive meshes and conservation of the connection between the grid points are the most attractive properties of optimal transport problem [49, 38].

**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 [40], [3]. An interesting introduction to equidistribution grid adaptation based on the Monge-Kantorovich problem (MKP) can be found in [20] and the references therein. The connection between MKP and the resolution of the Monge-Ampère equation (MA) is now quite well-known [20, 33]. Recently, many authors proposed equidistribution grid generation methods using the MA [20, 19, 12, 51, 50, 13, 10, 55, 54].
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 [4].

In [15], the authors solve the MA using an augmented Lagrangian approach. In [16], they proposed a least square approach. Much more details can be found in [17]. In [43], the authors proposed an iterative solver based on the divergence form and Newton method. While, in [20] the authors used an iterative Newton-Krylov solver with preconditioning coupled with a Finite Difference method. In [2] [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 [8, 51, 50].

In [56], the authors proposed three new methods based on Fourier integral, second-order divergence and a convolution forms.

**Problem formulation.** Now, we consider a fixed computational domain $\Omega_c$, and $\Omega_p$ the physical domain where the underlying PDE is posed and the initial mapping $\mathbf{F}: \Omega_c \ni \xi \rightarrow x \in \Omega_p$. The main idea behind an optimum mesh is that it should be closest to a uniform mesh in a suitable norm [32]. The equidistributed grids are constructed by solving the Monge-Ampère equation which is equivalent to the $L^2$ Monge-Kantorovich problem [18, 32].

Hence, the adaptive meshes are obtained using a new mapping $\mathbf{F}': \Omega_c \ni \xi \rightarrow x' \in \Omega_p$, that maps the parametric domain $\Omega_c$ into the physical domain $\Omega_p$, defined as the composition of initial mapping $\mathbf{F}: \Omega_c \ni \xi \rightarrow x' \in \Omega_p$ and B-spline geometrical map $\mathbf{\Psi}: \Omega_c \ni \xi \rightarrow \xi' \in \Omega_c$ derived from a resolution of Monge-Ampère equation. Indeed, for our application, we have $\Omega_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-Ampère 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.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Illustration of the Adaptive Mapping idea. The old mapping $\mathbf{F}: \xi' \to x$ and the new mapping $\mathbf{F}': \xi \to x'$ map the logical domain (patch) onto the same physical domain $\Omega_p$. The mapping $\Psi: \xi \to \xi'$, maps the old (uniform) mesh onto the new one in parametric domain.</figcaption></figure>

## 2. B-SPLINES AND GEOMETRICAL REPRESENTATION

IsoGeometric analysis (IGA) has been introduced in [36], 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) \qquad \Xi = \{0 = \xi_1, \xi_2, \dots, \xi_{n+p+1} = 1\} $$

where $n$ is the number of basis functions necessary to describe it. Here we work only with open knot vectors, which means that first and last knots in $\Xi$ have multiplicity $p+1$, so that

$$ (2.2) \qquad \Xi = \underbrace{\{\zeta_1, \dots, \zeta_1, \zeta_2, \dots, \zeta_2, \dots, \zeta_m, \dots, \zeta_m\}}_{p+1} $$

where $1 \le r_j \le p+1$ for $j = 2, \dots, m-1$ and $n = p+1 + \sum_{j=2}^{m-1} r_j$. Univariate B-spline basis functions $B_i^p$, $i = 1, \dots, n$ are defined recursively by the well known Cox-de Boor recursion formula:

$$ (2.3) \qquad B_i^0(\xi) = \begin{cases} 1 & \text{if } \xi \in [\xi_i, \xi_{i+1}) \\ 0 & \text{otherwise} \end{cases} $$

$$ (2.4) \qquad B_i^p(\xi) = \frac{\xi - \xi_i}{\xi_{i+p} - \xi_i} B_i^{p-1}(\xi) + \frac{\xi_{i+p+1} - \xi}{\xi_{i+p+1} - \xi_{i+1}} B_{i+1}^{p-1}(\xi) $$
where in 2.4, we adopt the convention $0/0 = 0$. These basis functions are piecewise polynomials of degree $p$ on the subdivision $\{\zeta_1, ..., \zeta_m\}$, where at $\zeta_j$ they have $\mu_j := p - r_j$ continuous derivative. Therefore, the vector $\mu = \{\mu_1, ..., \mu_m\}$ collects the regularity at the internal knots, with $\mu_1 = \mu_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 $\mu = \{-1, 2, 1, 0, -1\}$.

Respectively, by means of tensor products, a multi-dimensional B-spline can be constructed as $B_{j_1, j_2, ..., j_d}^{p_1, p_2, ..., p_d} = \otimes_{i=1}^d B_{j_i}^{p_i} [36, 14]$.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. (left) A cubic B-Spline curve and its control points using the knot vector $\Xi = \{0000 \quad \frac{1}{4} \quad \frac{1}{2} \quad \frac{1}{2} \quad \frac{3}{4} \quad \frac{3}{4} \quad 1111\}$, (right) the corresponding B-Splines.</figcaption></figure>

*Geometrical representation.* A single patch domain $\Omega_p$ is a B-spline region associated with the control points $C_{j_1, j_2, ..., j_d}$. We introduce the B-spline geometrical map $\mathbf{F}: \Omega_c \to \Omega_p$ by

$$ (2.5) \qquad \mathbf{F}(\xi) = \sum_{j_1, j_2, \dots, j_d} B_{j_1, j_2, \dots, j_d}^{p_1, p_2, \dots, p_d}(\xi) C_{j_1, j_2, \dots, 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 $\mathcal{Q}_h$ a family of meshes on $\Omega_c$ by

$$ (2.6) \qquad \mathcal{Q}_h = \{Q = \otimes_{i=1}^d (\zeta_{j_i, i}, \zeta_{j_i+1, i}), 1 \le j_i \le m_i - 1, d = 2, 3\} $$

each element $Q \in \mathcal{Q}_h$ is mapped into an element

$$ (2.7) \qquad k = \mathbf{F}(Q) = \{\mathbf{F}(\xi), \xi \in Q\}, $$

and analogously $\tilde{Q}$, the support extension of $Q$, is mapped into

$$ (2.8) \qquad \tilde{k} = \mathbf{F}(\tilde{Q}). $$

We then introduce the mesh $\mathcal{T}_h$ in the physical domain $\Omega_p$, defined by:

$$ (2.9) \qquad \mathcal{T}_h := \{k = \mathbf{F}(Q), Q \in \mathcal{Q}_h\}. $$

Finally, we define B-spline spaces needed in the sequel as follows:
$$ (2.10) \qquad S_{\mu}^{p} = \operatorname{span}\{B_{i}^{p}, \text{ for } i = 1, \dots, n\} $$

and

$$ (2.11) \qquad \left\{ \frac{dv}{d\xi} : v \in S_{\mu}^{p} \right\} = S_{\mu-1}^{p-1} $$

where we adopt the notation $\mu - 1 = \{\mu_1, \mu_2 - 1, \dots, \mu_{m-1} - 1, \mu_m\}$.

Then, we consider the discrete subspace of $H^s(\Omega_c)$, given by

$$ (2.12) \qquad V_h = \bigotimes_{i=1}^{d} S_{\mu_i}^{p_i}. $$

Let us introduce the discrete subspace of $H(\operatorname{div}, \Omega)$ and $L^2(\Omega)$ denoted by $V_h(\operatorname{div}, \Omega)$ and $V_h(L^2, \Omega)$ successively by

-For 2D case:

$$ (2.13) \qquad \left\{ \begin{array}{l} V_h(\operatorname{div}, \Omega) = \begin{pmatrix} S_{\mu_1}^{p_1} \otimes S_{\mu_2}^{p_2-1} \\ S_{\mu_1-1}^{p_1-1} \otimes S_{\mu_2}^{p_2} \end{pmatrix} \\ V_h(L^2, \Omega) = S_{\mu_1-1}^{p_1-1} \otimes S_{\mu_2-1}^{p_2-1} \end{array} \right. $$

-For 3D case:

$$ (2.14) \qquad \left\{ \begin{array}{l} V_h(\operatorname{div}, \Omega) = \begin{pmatrix} S_{\mu_1}^{p_1} \otimes S_{\mu_2-1}^{p_2-1} \otimes S_{\mu_3-1}^{p_3-1} \\ S_{\mu_1-1}^{p_1-1} \otimes S_{\mu_2}^{p_2} \otimes S_{\mu_3-1}^{p_3-1} \\ S_{\mu_1-1}^{p_1-1} \otimes S_{\mu_2-1}^{p_2-1} \otimes S_{\mu_3}^{p_3} \end{pmatrix} \\ V_h(L^2, \Omega) = S_{\mu_1-1}^{p_1-1} \otimes S_{\mu_2-1}^{p_2-1} \otimes S_{\mu_3-1}^{p_3-1} \end{array} \right. $$

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(\operatorname{grad}, \Omega) \xrightarrow{\nabla \times} V_h(\operatorname{div}, \Omega) \xrightarrow{\nabla \cdot} V_h(L^2, \Omega) $$

While in 3d, we have the following sequence

$$ V_h(\operatorname{grad}, \Omega) \xrightarrow{\nabla} V_h(\operatorname{curl}, \Omega) \xrightarrow{\nabla \times} V_h(\operatorname{div}, \Omega) \xrightarrow{\nabla \cdot} V_h(L^2, \Omega) $$

where

$$ (2.15) \qquad V_h(\operatorname{curl}, \Omega) = \begin{pmatrix} S_{\mu_1-1}^{p_1-1} \otimes S_{\mu_2}^{p_2} \otimes S_{\mu_3}^{p_3} \\ S_{\mu_1}^{p_1} \otimes S_{\mu_2-1}^{p_2-1} \otimes S_{\mu_3}^{p_3} \\ S_{\mu_1}^{p_1} \otimes S_{\mu_2}^{p_2} \otimes S_{\mu_3-1}^{p_3-1} \end{pmatrix} $$

is the discrete subspace of $H(\operatorname{curl}, \Omega)$.

# 3. EQUIDISTRIBUTION MESH GENERATION BASED ON MONGE-AMPÈRE EQUATION

## 3.1. The Monge-Kantorovich Problem.
**Definition 3.1** ($L^2$ Monge-Kantorovich problem). Let $\rho_0$ and $\rho_1$ be two given densities of equal masses, defined in $\mathbb{R}^d$.

Find a mapping $\mathbf{x}' = \mathbf{T}(\mathbf{x})$, $\mathbf{x}, \mathbf{x}' \in \mathbb{R}^d$, that transfers the density $\rho_0$ to $\rho_1$ and minimizes the transport cost

$$ (3.1) \qquad \mathcal{J}[\mathbf{T}] = \int_{\mathbb{R}^d} |\mathbf{T}(\mathbf{x}) - \mathbf{x}'|^2 \rho_0(\mathbf{x}) \, d\mathbf{x} $$

The density transfer means that, for all Borelian set $B$ in $\mathbb{R}^d$

$$ \int_{\mathbf{T}^{-1}(B)} \rho_0(x) \, dx = \int_B \rho_1(x) \, dx. $$

Then, it follows from a change of variable that

$$ \rho_1(\mathbf{T}) \det(\nabla \mathbf{T}) = \rho_0. $$

Brenier's theorem [6], shows that $\mathbf{T} = \nabla\phi$ for given convex potential function $\phi$. Hence, we get the Monge-Ampere equation:

$$ (3.2) \qquad \rho_1(\nabla\phi) \det(H(\phi)) = \rho_0, $$

where $H(\phi)$ is the Hessian matrix of $\phi$.

**Remark 3.1.** For arbitrary $\rho_0$ and $\rho_1$ defined in the physical domain $\Omega_p$, let

$$ (3.3) \qquad \sigma = \int_{\Omega_p} \rho_1(x) \, dx / \int_{\Omega_p} \rho_0(x) \, dx, $$

we rewrite the Monge-Ampere equation as follows

$$ (3.4) \qquad \rho_1(\nabla\phi) \det(H(\phi)) = \sigma. $$

**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 $\mathbf{F}$ is bijective, so we take the two densities defined in $\Omega_c$ as follows

(1) $\sigma \det(\nabla \mathbf{F}):$ corresponding to a uniform mesh in $\Omega_c$,

(2) $\rho_1 \circ \mathbf{F} \det(\nabla \mathbf{F}):$ corresponding to what we call the inverse image by the initial mapping $\mathbf{F}$ of the adapted mesh in $\Omega_p$,

where $\sigma$ is given by 3.3, the two densities are of equal masses. Then, according to Brenier's theorem [7] there exists an optimal mapping $\mathbf{\Psi}$ defined in the parametric domain $\Omega_c$ which is the gradient of a given convex potential function $\tilde{\phi}$, such that

$$ (3.5) \qquad \rho_1(\mathbf{F} \circ \nabla \tilde{\phi}) \det(\nabla \mathbf{F}) \circ \nabla \tilde{\phi} \det(H(\tilde{\phi})) = \sigma \det(\nabla \mathbf{F}). $$

Moreover, our mesh adaptation technique and the one existing in the literature [18] are equivalent.

*Proof.* Using $\rho_0 = 1$ correspond to uniform mesh in $\Omega_p$ and assuming that $\Omega_p$ is a convex domain, then there exists $\mathbf{T}: \Omega_p \to \Omega_p$ satisfying the Monge-Ampere equation 3.4, which corresponds to the adaptation in the physical domain, hence we have

$$ (3.6) \qquad \rho_1(\mathbf{T}) \det(\nabla \mathbf{T}) = \sigma. $$
using the fact that $\mathbf{F}$ is bijective, let $\mathbf{T}_0 : \Omega_c \to \Omega_c$ such that

$$ (3.7) \qquad \mathbf{T} = \mathbf{F} \circ \mathbf{T}_0 \circ \mathbf{F}^{-1}, $$

as a consequence, by a change of variable in 3.6, we have

$$ (3.8) \qquad \rho_1(\mathbf{F} \circ \mathbf{T}_0 \circ \mathbf{F}^{-1}) \det (\nabla(\mathbf{F} \circ \mathbf{T}_0 \circ \mathbf{F}^{-1})) = \sigma $$

using the fact that

$$ \begin{aligned} \nabla(\mathbf{F} \circ \mathbf{T}_0 \circ \mathbf{F}^{-1}) &= (\nabla \mathbf{F}^{-1})^T \nabla(\mathbf{F} \circ \mathbf{T}_0) \circ \mathbf{F}^{-1} \\ &= (\nabla \mathbf{F}^{-1})^T (\nabla \mathbf{T}_0)^T \circ \mathbf{F}^{-1} (\nabla \mathbf{F}) \circ \mathbf{T}_0 \circ \mathbf{F}^{-1}, \end{aligned} $$

and $\det(AB) = \det(A)\det(B)$, we get

$$ (3.9) \qquad \rho_1(\mathbf{F} \circ \mathbf{T}_0 \circ \mathbf{F}^{-1}) \det ((\nabla \mathbf{F}^{-1})^T) \det ((\nabla \mathbf{T}_0)^T) \circ \mathbf{F}^{-1} \det (\nabla \mathbf{F}) \circ \mathbf{T}_0 \circ \mathbf{F}^{-1} = \sigma $$

and therefore, by applying $\mathbf{F}$ we have

$$ (3.10) \qquad \rho_1(\mathbf{F} \circ \mathbf{T}_0) \det ((\nabla \mathbf{F}^{-1}) \circ \mathbf{F}) \det (\nabla \mathbf{T}_0) \det (\nabla \mathbf{F}) \circ \mathbf{T}_0 = \sigma. $$

From differential calculus, we know that $(\nabla \mathbf{F}^{-1}) \circ \mathbf{F} = (\nabla \mathbf{F})^{-1}$, then

$$ (3.11) \qquad \rho_1(\mathbf{F} \circ \mathbf{T}_0) \det (\nabla \mathbf{F}) \circ \mathbf{T}_0 \det (\nabla \mathbf{T}_0) = \sigma \det (\nabla \mathbf{F}). $$

Finally, from the uniqueness of the optimal mapping we conclude that $\mathbf{T}_0 = \mathbf{\Psi}$. $\square$

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

$$ (3.12) \qquad \det (\nabla \mathbf{F}) \circ \nabla \tilde{\phi} \approx \det (\nabla \mathbf{F}), $$

then, we simplify the equation 3.5 by considering the new form:

$$ (3.13) \qquad \rho_1(\mathbf{F} \circ \nabla \tilde{\phi}) \det(H(\tilde{\phi})) = \sigma. $$

The boundary conditions are obtained by requiring that $\mathbf{\Psi}$ maps $\partial\Omega_c$ to $\partial\Omega_c$, i.e

$$ \nabla \tilde{\phi}(\partial \Omega_c) = \partial \Omega_c $$

which writes:

$$ (3.14) \qquad \nabla \tilde{\phi} \cdot \vec{n} = x \cdot \vec{n}. $$

Finally, we denote

$$ (3.15) \qquad f = \frac{\sigma}{\rho_1(\mathbf{F})}, $$

the desired convex potential function is the solution of Monge-Ampere system:

$$ (3.16) \qquad \left\{ \begin{array}{ll} \det(H(\tilde{\phi})) & = f(\nabla \tilde{\phi}) \quad \text{in } \Omega_c, \\ \nabla \tilde{\phi} \cdot \vec{n} & = x \cdot \vec{n} \quad \text{on } \partial \Omega_c, \\ \tilde{\phi} \text{ is convex.} & \end{array} \right. $$

A solution of 3.16 exists and is unique up to an additive constant [42]. From now on, the notation $\phi$ stands for $\tilde{\phi}$ solution of the latest system 3.16 and $\rho$ stands for $\rho_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 [5] which guarantees convergence when a second member is strictly positive and regular. The authors noticed that a solution of the Monge-Ampère equation is a fixed point of the operator $T$.

$$ (4.1) \qquad \begin{aligned} T: H^2(\Omega_c) &\to H^2(\Omega_c) \\ v &\to T(v) = (\Delta)^{-1}\mathcal{G}(v) \end{aligned} $$

where $\mathcal{G}(v) = ((\Delta v)^d + d! (f(\nabla v)) - \det(H(v)))^{\frac{1}{d}}$, $d$ is spatial dimension and $f$ is given by 3.15, that rewrite as

$$ (4.2) \qquad \left\{ \begin{array}{ll} -\Delta\phi & = -\mathcal{G}(\phi) & \text{in } \Omega_c, \\ \nabla\phi \cdot \vec{n} & = x \cdot \vec{n} & \text{on } \partial\Omega_c, \\ \phi & \text{is convex.} \end{array} \right. $$

Thus, the solution is given by Picard iteration as follow:

(1) Given an initial value $\phi^0$.

(2) Compute $\phi^{n+1}$ as the solution of the problem:

$$ (4.3) \qquad \left\{ \begin{array}{ll} -\Delta\phi^{n+1} + \varepsilon\phi^{n+1} & = -\mathcal{G}(\phi^n) & \text{in } \Omega_c, \\ \nabla\phi^{n+1} \cdot \vec{n} & = x \cdot \vec{n} & \text{on } \partial\Omega_c, \\ \phi^n & \text{is convex,} \end{array} \right. $$

where $\varepsilon = 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 three-dimensional cases.

The variational formulation of 4.3, is:

$$ (4.4) \qquad \int_{\Omega_c} \nabla\phi^{n+1}\nabla v \, dx + \varepsilon \int_{\Omega_c} \phi^{n+1}v \, dx = - \int_{\Omega_c} \mathcal{G}(\phi^n)v \, dx + \int_{\partial\Omega_c} x \cdot \vec{n} v \, d\sigma, \text{ for all } v \in H^2(\Omega_c) $$

and a Galerkin approximation of 4.3 is given by

$$ (4.5) \qquad \int_{\Omega_c} \nabla\phi_h^{n+1}\nabla v^h \, dx + \varepsilon \int_{\Omega_c} \phi_h^{n+1}v^h \, dx = - \int_{\Omega_c} \mathcal{G}(\phi_h^n)v^h \, dx + \int_{\partial\Omega_c} x \cdot \vec{n} v^h \, d\sigma, \text{ for all } v \in V_h \cap H^2(\Omega_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 {$\phi^n$}$_n$ be a solution sequence of Picard algorithm in $H^2(\Omega_c)$, then it holds that:

$$
(4.6) \qquad |e^{n,h}|_{H^1(\Omega_c)} \leq Ch^{l-1},
$$

where $e^{n,h} = \phi^n - \phi_h^n$ is the approximation error at iteration $n$ and $l$ is a Sobolev regularity exponent of $\phi^n$ ($\phi^n \in H^l(\Omega_c)$, $l \ge 2$).

*Proof.* By recurrence, we have for $n = 0$

$$
|\phi^0 - \phi_h^0|_{H^1(\Omega_c)} \leq Ch^{l-1},
$$

holds by choice. So, we first prove that we have the following implication

$$
(4.7) \qquad |\phi^n - \phi_h^n|_{H^1(\Omega_c)} \le Ch^{l-1} \implies \|\mathcal{G}(x, \phi^n) - \mathcal{G}(x, \phi_h^n)\|_{L^2(\Omega_c)} \le Ch^{l-2}, \text{ for } n \in \mathbb{N},
$$

recalling that

$$
\|H(\phi^n)\|^2 = (\Delta\phi^n)^d - d! \det(H(\phi^n)) \ge 0,
$$

and, the triangular inequality leads to

$$
\begin{align*}
|\mathcal{G}(x, \phi^n) - \mathcal{G}(x, \phi_h^n)|^2 &= \left| \sqrt{||H(\phi^n)||^2 + d! \hat{f}(\nabla \phi^n)} - \sqrt{||H(\phi_h^n)||^2 + d! \hat{f}(\nabla \phi_h^n)} \right|^2 \\
&\le 2 \left| \sqrt{||H(\phi^n)||^2 + d! \hat{f}(\nabla \phi^n)} - \sqrt{||H(\phi^n)||^2 + d! \hat{f}(\nabla \phi_h^n)} \right|^2 \\
&\quad + 2 \left| \sqrt{||H(\phi^n)||^2 + d! \hat{f}(\nabla \phi_h^n)} - \sqrt{||H(\phi_h^n)||^2 + d! \hat{f}(\nabla \phi_h^n)} \right|^2 \\
&:= 2A^2 + 2B^2.
\end{align*}
$$

Using the inequality $|\sqrt{a+b} - \sqrt{c+b}| \le |\sqrt{a} - \sqrt{c}|$ for $a, b, c \ge 0$, we have

$$
(4.8) \qquad B \le \frac{\|H(\phi^n)\| - \|H(\phi_h^n)\|}{\|H(\phi^n)\|}
$$

and

$$
\begin{align*}
A &\le \sqrt{d!} \left| \sqrt{\hat{f}(\nabla \phi^n)} - \sqrt{\hat{f}(\nabla \phi_h^n)} \right| \\
&\le \sqrt{d!} \left| \frac{\sqrt{\sigma}}{\sqrt{\rho(\nabla \phi^n)}} - \frac{\sqrt{\sigma}}{\sqrt{\rho(\nabla \phi_h^n)}} \right| \\
&\le \sqrt{d! \sigma} \frac{\left| \sqrt{\rho(\nabla \phi^n)} - \sqrt{\rho(\nabla \phi_h^n)} \right|}{\sqrt{\rho(\nabla \phi^n) \rho(\nabla \phi_h^n)}}
\end{align*}
$$

since that $\rho \ge \gamma > 0$, and $t \to \sqrt{t}$ is 1-Lipschitz function in $[1, +\infty)$ then

$$
A \le \sqrt{\frac{d! \sigma}{\gamma}} |\rho(\nabla \phi^n) - \rho(\nabla \phi_h^n)|
$$

if we assume that $\rho$ is $L_\rho$-Lipschitz

$$
(4.9) \qquad A \le L_{\rho} \sqrt{\frac{d! \sigma}{\gamma}} |\nabla \phi^n - \nabla \phi_h^n|.
$$

As a consequence of 4.8 and 4.9

$$
\|\mathcal{G}(x, \phi^n) - \mathcal{G}(x, \phi_h^n)\|_{L^2(\Omega_c)} \le 2|\phi^n - \phi_h^n|_{H^2(\Omega_c)} + 2L_{\rho}\sqrt{\frac{d!\sigma}{\gamma}}|\phi^n - \phi_h^n|_{H^1(\Omega_c)}.
$$
Now, taking the Projector $\Pi_{\mathcal{V}_h} : L^2(\Omega_c) \to \mathcal{V}_h$, we get

$$ \begin{aligned} |\phi^n - \phi_h^n|_{H^2(\Omega_c)} &\le |(\phi^n - \phi_h^n) - \Pi_{\mathcal{V}_h}(\phi^n - \phi_h^n)|_{H^2(\Omega_c)} + |\Pi_{\mathcal{V}_h}(\phi^n - \phi_h^n)|_{H^2(\Omega_c)} \\ &= I + \Pi. \end{aligned} $$

With $C^2$ regularity in B-spline function, from Proposition 3.2 [52]

$$ \begin{aligned} I &= |(\phi^n - \phi_h^n) - \Pi_{\mathcal{V}_h}(\phi^n) + \phi_h^n|_{H^2(\Omega_c)} \\ &= |\phi^n - \Pi_{\mathcal{V}_h}(\phi^n)|_{H^2(\Omega_c)} \\ &\le Ch^{l-2}||\phi^n||_{H^l(\Omega_c)}. \end{aligned} $$

By the usual inverse inequality for polynomials, we have

$$ \begin{aligned} \Pi^2 &= \sum_{k \in \mathcal{T}_h} |\Pi_{\mathcal{V}_h}(\phi^n - \phi_h^n)|_{H^2(k)}^2 \\ &\le C_{shap} \sum_{k \in \mathcal{T}_h} h_K^{-2} |\Pi_{\mathcal{V}_h}(\phi^n - \phi_h^n)|_{H^1(k)}^2 \\ (4.10) \qquad \qquad \qquad \qquad \qquad \qquad \qquad \le C_{shap} \underline{h}^{-2} |\phi^n - \phi_h^n|_{H^1(\Omega_c)}^2, \end{aligned} $$

where $\underline{h} = \min(h_k)$, so if we assume that for a positive constant $\theta$,

$$ h \le \theta \underline{h}. $$

According to the hypothesis of 4.7, we have

$$ I + II \le Ch^{l-2}, $$

where $C$ depend only on $|\phi^n|_{H^l(\Omega_c)}$ and $\theta$.

Returning now to the proof of our proposition, let $\phi^{n+1}$ be the solution at $n+1$ iteration in Picard algorithm and $v_h \in \mathcal{V}_h$, we have for some constant $c_{n+1}$

$$ \begin{aligned} |\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)}^2 &= \int_{\Omega_c} \nabla(\phi^{n+1} - \phi_h^{n+1})\nabla(\phi^{n+1} - v_h + c_{n+1}) \,dx \\ &\quad + \int_{\Omega_c} \nabla(\phi^{n+1} - \phi_h^{n+1})\nabla(v_h - c_{n+1} - \phi_h^{n+1}) \,dx \\ &= \int_{\Omega_c} \nabla(\phi^{n+1} - \phi_h^{n+1})\nabla(\phi^{n+1} - v_h) \,dx \\ &\quad + \int_{\Omega_c} (\mathcal{G}(x, \phi_h^n) - \mathcal{G}(x, \phi^n))(v_h - c_{n+1} - \phi_h^{n+1}) \,dx, \end{aligned} $$

by Holder's inequality, we have

$$ \begin{aligned} |\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)}^2 &\le |\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)} |\phi^{n+1} - v_h|_{H^1(\Omega_c)} \\ &\quad + ||\mathcal{G}(x, \phi^n) - \mathcal{G}(x, \phi_h^n)||_{L^2(\Omega_c)} ||v_h - c_{n+1} - \phi_h^{n+1}||_{L^2(\Omega_c)}, \end{aligned} $$

taking to account 4.7, we get by the triangular inequality

$$ \begin{aligned} |\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)}^2 &\le |\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)} |\phi^{n+1} - v_h|_{H^1(\Omega_c)} + Ch^{l-2}||v_h - c_{n+1} - \phi_h^{n+1}||_{L^2(\Omega_c)} \\ &\le |\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)} |\phi^{n+1} - v_h|_{H^1(\Omega_c)} + Ch^{l-2}||\phi^{n+1} - v_h||_{L^2(\Omega_c)} \\ &\quad + Ch^{l-2}||\phi^{n+1} - \phi_h^{n+1} - c_{n+1}||_{L^2(\Omega_c)} \end{aligned} $$

where $v_h \in \mathcal{V}_h$, according to Proposition 3.2 [52], for $v_h = \Pi_{\mathcal{V}_h}(\phi^{n+1})$ we have

$$ |\phi^{n+1} - \Pi_{\mathcal{V}_h}(\phi^{n+1})|_{H^1(\Omega_c)} \le Ch^{l-1}, $$
and

$$||\phi^{n+1} - \Pi_{\nu_h}(\phi^{n+1})||_{L^2(\Omega_c)} \leq Ch^l.$$

Then,

$$(4.11) \qquad |\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)}^2 \leq Ch^{l-1}|\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)} + Ch^{2l-2} + Ch^{l-2}||\phi^{n+1} - \phi_h^{n+1} - c_{n+1}||_{L^2(\Omega_c)}.$$

Let's now introduce the following problem:

$$(4.12) \qquad \begin{cases} -\Delta\psi = \phi^{n+1} - \phi_h^{n+1} - c_{n+1} & \text{in } \Omega_c, \\ \nabla\psi \cdot \vec{n} = 0 & \text{on } \partial\Omega_c. \end{cases}$$

The compatibility condition holds for $c_{n+1} = \int_{\Omega_c} (\phi^{n+1} - \phi_h^{n+1}) 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) \qquad \int_{\Omega_c} \nabla\psi\nabla v \,dx = \int_{\Omega_c} (\phi^{n+1} - \phi_h^{n+1} - c_{n+1})v \,dx \quad \text{for all } v \in H^1(\Omega_c)$$

by assuming that

$$(4.14) \qquad \int_{\Omega_c} (\mathcal{G}(x, \phi^n) - \mathcal{G}(x, \phi_h^n)) \,dx \neq 0$$

we introduce the constant $c$ such that

$$(4.15) \qquad \int_{\Omega_c} (\mathcal{G}(x, \phi^n) - \mathcal{G}(x, \phi_h^n))(\Pi_{\nu_h}(\psi) - c) \,dx = 0.$$

Hence, from the variational formulation 4.5, for $v^h = \Pi_{\nu_h}(\psi) - c$, we have

$$(4.16) \qquad \int_{\Omega_c} \nabla(\phi^{n+1} - \phi_h^{n+1} - c_{n+1})\nabla(\Pi_{\nu_h}(\psi) - c) \,dx = \int_{\Omega_c} (\mathcal{G}(x, \phi^n) - \mathcal{G}(x, \phi_h^n))(\Pi_{\nu_h}(\psi) - c) \,dx \\ = 0.$$

As a consequence, by replacing the results in 4.13, we get for $v = \phi^{n+1} - \phi_h^{n+1} - c_{n+1}$

$$(4.17) \qquad \int_{\Omega_c} (\phi^{n+1} - \phi_h^{n+1} - c_{n+1})^2 \,dx = \int_{\Omega_c} \nabla\psi\nabla(\phi^{n+1} - \phi_h^{n+1}) \,dx \\ (4.18) \qquad = \int_{\Omega_c} \nabla(\psi - (\Pi_{\nu_h}(\psi) - c))\nabla(\phi^{n+1} - \phi_h^{n+1}) \,dx \\ (4.19) \qquad \leq |\psi - (\Pi_{\nu_h}(\psi) - c)|_{H^1(\Omega_c)} |\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)}.$$

Furthermore, because $\phi^{n+1} - \phi_h^{n+1} \in H^l(\Omega_c)$ at least $\psi \in H^{l+2}(\Omega_c)$, then from the Proposition 3.2 [52] we have the following inequality :

$$|\psi - (\Pi_{\nu_h}(\psi) - c)|_{H^1(\Omega_c)} = |\psi - \Pi_{\nu_h}(\psi)|_{H^1(\Omega_c)} \\ \leq Ch^{l+1}||\psi||_{H^{l+2}(\Omega_c)}$$

which is equivalent to saying that we have :

$$(4.20) \qquad \int_{\Omega_c} (\phi^{n+1} - \phi_h^{n+1} - c_{n+1})^2 \,dx \leq Ch^{l+1}||\psi||_{H^{l+2}(\Omega_c)} |\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)}$$
as a consequence, replacing the results in inequality 4.11, we obtain
(4.21)

$$
|\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)}^2 \leq C h^{l-1} |\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)} + C h^{2l-2} + C h^{\frac{3(l-1)}{2}} |\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)}^{\frac{1}{2}}
$$

with Young's inequality (ab ≤ a² + b²/4) in the first and (ab ≤ 3a³/4 + b⁴/4) in the last
term of 4.21, we have

$$
(4.22) \quad \left(1 - \frac{1}{2}\right) \left|\phi^{n+1} - \phi_h^{n+1}\right|^2_{H^1(\Omega_c)} \le C h^{2l-2} + C h^{2l-2} + C h^{2l-2},
$$

finally, the estimation is given by

$$
|\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)}^2 \leq C h^{2(l-1)},
$$

where C constant depend on ||φⁿ+¹||ₙ⁽¹⁾(Ωₑ). □

**Corollary 4.1.** The error estimation for Monge-Ampere equation in H²-norm
fulfills

$$
(4.23) \qquad |\phi^{n+1} - \phi_h^{n+1}|_{H^2(\Omega_c)} \le C h^{l-2}.
$$

Proof. By taking the Projector Π<sub>V<sub>h</sub></sub> : L<sup>2</sup>(Ω<sub>c</sub>) —> V<sub>h</sub>, we get

$$
|\phi^{n+1} - \phi_h^{n+1}|_{H^2(\Omega_c)} \leq |\phi^{n+1} - \phi_h^{n+1} - \Pi_{V_h}(\phi^{n+1} - \phi_h^{n+1})|_{H^2(\Omega_c)} + |\Pi_{V_h}(\phi^{n+1} - \phi_h^{n+1})|_{H^2(\Omega_c)}
$$

by the proposition 3.2 [52]

$$
\begin{align*}
|\phi^{n+1} - \phi_h^{n+1} - \Pi_{V_h}(\phi^{n+1} - \phi_h^{n+1})|_{H^2(\Omega_c)} &= |\phi^{n+1} - \Pi_{V_h}(\phi^{n+1})|_{H^2(\Omega_c)} \\
&\le C h^{l-2} \|\phi^{n+1}\|_{H^1(\Omega_c)}.
\end{align*}
$$

By the usual inverse inequality as in 4.10, we have

$$
\begin{align*}
|\Pi_{V_h}(\phi^{n+1} - \phi_h^{n+1})|_{H^2(\Omega_c)} &\le C_\theta h^{-1} |\Pi_{V_h}(\phi^{n+1} - \phi_h^{n+1})|_{H^1(\Omega_c)} \\
&\le C_\theta h^{-1} |\phi^{n+1} - \phi_h^{n+1}|_{H^1(\Omega_c)}.
\end{align*}
$$

Therefore, according to the last proposition our result is achieved. □

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

$$
(4.24) \qquad |\phi - \phi_h^n|_{H^2(\Omega_c)} \le 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 [48], we describe in the sequel the fast diagonalization method in the case of Isogeometric Analysis. This method was first introduced in [45]. For the seek of simplicity, we shall consider the following Laplace problem,

$$
(4.25) \qquad \begin{cases} -\nabla^2 u + \tau u = f, & \Omega \\ u = 0, & \partial\Omega \end{cases}
$$

The Poisson problem and its solver shall be retrieved with τ = 0. After discretizing
the Laplace problem, we get the following linear system

$$
(4.26) \quad \mathcal{L}_{\tau} x := (K_1 \otimes M_2 \otimes M_3 + M_1 \otimes K_2 \otimes M_3 + M_1 \otimes M_2 \otimes K_3 + \tau M_1 \otimes M_2 \otimes 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) \qquad K_1 U_1 = M_1 U_1 D_1, \quad K_2 U_2 = M_2 U_2 D_2, \quad 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) \qquad U_1^T M_1 U_1 = I_1, \quad U_2^T M_2 U_2 = I_2, \quad U_3^T M_3 U_3 = I_3 $$

Therefor, (4.26) can be written as

$$ (4.29) \qquad (U_1 \otimes U_2 \otimes U_3)^{-1} (D_1 \otimes I_2 \otimes I_3 + I_1 \otimes D_2 \otimes I_3 + I_1 \otimes I_2 \otimes D_3 + \tau I_1 \otimes I_2 \otimes I_3) (U_1 \otimes U_2 \otimes 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):

<table><caption>**Algorithm 1: fast_diag**: Fast diagonalization method for Laplace problem</caption><thead><tr><th>Input</th><th>: $\mathcal{L}_{\tau}, b$</th></tr></thead><tbody><tr><td>Output:</td><td>: $x$</td></tr><tr><td>1</td><td>$\tilde{b} \leftarrow (U_1 \otimes U_2 \otimes U_3) b$</td></tr><tr><td>2</td><td>$\tilde{x} \leftarrow (D_1 \otimes I_2 \otimes I_3 + I_1 \otimes D_2 \otimes I_3 + I_1 \otimes I_2 \otimes D_3 + \tau I_1 \otimes I_2 \otimes I_3)^{-1} \tilde{b}$</td></tr><tr><td>3</td><td>$x \leftarrow (U_1 \otimes U_2 \otimes U_3)^T \tilde{x}$</td></tr></tbody></table>

## 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 [47], 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) \qquad \phi(\xi, \eta) = \exp\left(\frac{\xi^2 + \eta^2}{2}\right) $$

for which the source term is

$$ (5.2) \qquad f(\xi, \eta) = (1 + \xi^2 + \eta^2)\exp(\xi^2 + \eta^2). $$

<table><caption>TABLE 1. CPU-time needed to solve the Monge-Ampere equation and error evolution for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>H<sup>1</sup>-norm</th><th>H<sup>1</sup>-norm MG</th><th>CPU-time (s)</th><th>CPU-time (s) MG</th><th>Nbr-iter</th><th>Nbr-iter MG</th></tr></thead><tbody><tr><td>8</td><td>1.783e-04</td><td>1.783e-04</td><td>0.067</td><td>0.017</td><td>35</td><td>26</td></tr><tr><td>16</td><td>2.308e-05</td><td>2.308e-05</td><td>0.032</td><td>0.033</td><td>35</td><td>28</td></tr><tr><td>32</td><td>2.934e-06</td><td>2.934e-06</td><td>0.106</td><td>0.094</td><td>35</td><td>27</td></tr><tr><td>64</td><td>3.692e-07</td><td>3.694e-07</td><td>0.399</td><td>0.310</td><td>35</td><td>25</td></tr><tr><td>128</td><td>4.579e-08</td><td>4.595e-08</td><td>1.671</td><td>1.186</td><td>35</td><td>23</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>H<sup>1</sup>-norm</th><th>H<sup>1</sup>-norm MG</th><th>CPU-time (s)</th><th>CPU-time (s) MG</th><th>Nbr-iter</th><th>Nbr-iter MG</th></tr></thead><tbody><tr><td>8</td><td>6.037e-06</td><td>6.037e-06</td><td>0.075</td><td>0.020</td><td>35</td><td>26</td></tr><tr><td>16</td><td>4.020e-07</td><td>4.020e-07</td><td>0.046</td><td>0.041</td><td>35</td><td>26</td></tr><tr><td>32</td><td>2.571e-08</td><td>2.572e-08</td><td>0.148</td><td>0.121</td><td>35</td><td>24</td></tr><tr><td>64</td><td>2.187e-09</td><td>2.183e-09</td><td>0.564</td><td>0.393</td><td>35</td><td>22</td></tr><tr><td>128</td><td>1.693e-09</td><td>1.684e-09</td><td>2.297</td><td>1.404</td><td>35</td><td>20</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>H<sup>1</sup>-norm</th><th>H<sup>1</sup>-norm MG</th><th>CPU-time (s)</th><th>CPU-time (s) MG</th><th>Nbr-iter</th><th>Nbr-iter MG</th></tr></thead><tbody><tr><td>8</td><td>3.584e-07</td><td>3.584e-07</td><td>0.042</td><td>0.026</td><td>35</td><td>26</td></tr><tr><td>16</td><td>1.325e-08</td><td>1.325e-08</td><td>0.064</td><td>0.058</td><td>35</td><td>26</td></tr><tr><td>32</td><td>1.874e-09</td><td>1.875e-09</td><td>0.220</td><td>0.175</td><td>35</td><td>24</td></tr><tr><td>64</td><td>1.711e-09</td><td>1.712e-09</td><td>0.831</td><td>0.571</td><td>35</td><td>22</td></tr><tr><td>128</td><td>1.710e-09</td><td>1.710e-09</td><td>3.336</td><td>2.040</td><td>35</td><td>20</td></tr></tbody></table>

**Example 5.2** (Three-dimensional case). We test our solver by a convex manufactured solution on the cube case, defined as :

$$ (5.3) \qquad \phi(\xi, \eta, \gamma) = \exp\left(\frac{\xi^2 + \eta^2 + \gamma^2}{3}\right) $$

where the source term is

$$ (5.4) \qquad f(\xi, \eta, \gamma) = \frac{8}{81}(3 + 2(\xi^2 + \eta^2 + \gamma^2)) \exp(\xi^2 + \eta^2 + \gamma^2). $$
<table><caption>TABLE 2. CPU-time needed to solve the Monge-Ampere equation and error evolution for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>H<sup>1</sup>-norm</th><th>H<sup>1</sup>-norm MG</th><th>CPU-time (s)</th><th>CPU-time (s) MG</th><th>Nbr-iter</th><th>Nbr-iter MG</th></tr></thead><tbody><tr><td>8</td><td>5.885e-05</td><td>5.885e-05</td><td>0.214</td><td>0.187</td><td>14</td><td>10</td></tr><tr><td>16</td><td>7.672e-06</td><td>7.672e-06</td><td>1.290</td><td>0.987</td><td>14</td><td>10</td></tr><tr><td>32</td><td>9.793e-07</td><td>9.793e-07</td><td>9.671</td><td>6.683</td><td>14</td><td>9</td></tr><tr><td>64</td><td>1.235e-07</td><td>1.236e-07</td><td>76.92</td><td>40.88</td><td>14</td><td>7</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>H<sup>1</sup>-norm</th><th>H<sup>1</sup>-norm MG</th><th>CPU-time (s)</th><th>CPU-time (s) MG</th><th>Nbr-iter</th><th>Nbr-iter MG</th></tr></thead><tbody><tr><td>8</td><td>1.758e-06</td><td>1.758e-06</td><td>0.528</td><td>0.472</td><td>14</td><td>10</td></tr><tr><td>16</td><td>1.184e-07</td><td>1.184e-07</td><td>3.547</td><td>2.680</td><td>14</td><td>10</td></tr><tr><td>32</td><td>7.658e-09</td><td>7.658e-09</td><td>27.35</td><td>16.42</td><td>14</td><td>8</td></tr><tr><td>64</td><td>9.609e-10</td><td>9.617e-10</td><td>216.1</td><td>98.91</td><td>14</td><td>6</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>H<sup>1</sup>-norm</th><th>H<sup>1</sup>-norm MG</th><th>CPU-time (s)</th><th>CPU-time (s) MG</th><th>Nbr-iter</th><th>Nbr-iter MG</th></tr></thead><tbody><tr><td>8</td><td>9.064e-08</td><td>9.064e-08</td><td>1.206</td><td>1.181</td><td>14</td><td>10</td></tr><tr><td>16</td><td>3.266e-09</td><td>3.266e-09</td><td>9.118</td><td>6.787</td><td>14</td><td>10</td></tr><tr><td>32</td><td>8.756e-10</td><td>8.757e-10</td><td>71.63</td><td>42.58</td><td>14</td><td>8</td></tr><tr><td>64</td><td>8.561e-10</td><td>8.563e-10</td><td>582.9</td><td>270.6</td><td>14</td><td>6</td></tr></tbody></table>

**Remark 5.1.** For numerical purpose $d!$ is replaced by $\delta$ such that

$$ (5.5) \qquad \mathcal{G}(\phi) = \left( (\Delta\phi)^d + \delta(f(\nabla\phi)) - \det(H(\phi)) \right)^{\frac{1}{d}} $$

Then, we see that a choice of $\delta$ in three-dimension affects the number of iterations of the Picard algorithm to achieve the same precision Fig. 3.

<figure><img src="image_4.png" /><figcaption>FIGURE 3. $L^2$ and $H^1$-residual and $H^1$-error profiles with different choices of $\delta = 6$ (left), $\delta = 9$ (middle) and $\delta = 27$ (right) using bi-cubic B-splines.</figcaption></figure>
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² residual error in Picard's iteration is fixed at 10⁻⁸.

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 = \int_{\Omega_c} |\det(H(\phi)) - \hat{f}(\nabla(\phi))|^2 dx \qquad Qual^2 = \int_{\Omega_c} |\nabla(\phi(x)) - x|^2 / \hat{f}(\nabla(\phi(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.

$$\rho(\xi, \eta) = 1 + 5 \exp(-50|(\xi - 0.5)^2 + (\eta - 0.5)^2 - 0.09|),$$

<table><caption>TABLE 3. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>4.195e-02</td><td>0.059</td><td>6.030e-02</td><td>2.784e-01</td><td>1.796e+00</td></tr><tr><td>32</td><td>1.124e-02</td><td>0.080</td><td>6.003e-02</td><td>2.759e-01</td><td>1.627e+00</td></tr><tr><td>64</td><td>3.023e-03</td><td>0.257</td><td>6.005e-02</td><td>2.707e-01</td><td>1.625e+00</td></tr><tr><td>128</td><td>1.090e-03</td><td>0.934</td><td>6.005e-02</td><td>2.711e-01</td><td>1.625e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>2.391e-02</td><td>0.079</td><td>5.993e-02</td><td>2.864e-01</td><td>1.674e+00</td></tr><tr><td>32</td><td>6.038e-03</td><td>0.123</td><td>6.005e-02</td><td>2.803e-01</td><td>1.626e+00</td></tr><tr><td>64</td><td>2.345e-03</td><td>0.393</td><td>6.005e-02</td><td>2.757e-01</td><td>1.625e+00</td></tr><tr><td>128</td><td>7.728e-04</td><td>1.419</td><td>6.005e-02</td><td>2.730e-01</td><td>1.625e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>2.494e-02</td><td>0.085</td><td>5.996e-02</td><td>2.945e-01</td><td>1.671e+00</td></tr><tr><td>32</td><td>6.633e-03</td><td>0.178</td><td>6.005e-02</td><td>2.795e-01</td><td>1.628e+00</td></tr><tr><td>64</td><td>2.119e-03</td><td>0.569</td><td>6.005e-02</td><td>2.755e-01</td><td>1.625e+00</td></tr><tr><td>128</td><td>8.364e-04</td><td>1.936</td><td>6.005e-02</td><td>2.733e-01</td><td>1.625e+00</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>FIGURE 4. (left) Adaptive meshes and (right) the monitor function.</figcaption></figure>

**Example 5.4.** This example in the unit square corresponds to a spiral centered at (0.7, 0.5) with very tight arms Fig. 5 (right).

$$ (5.7) \qquad \rho(\xi, \eta) = 1 + 9/(1 + (10r(\xi, \eta) \cos(\theta(\xi, \eta) - 20r(\xi, \eta)))^2), $$

where $\theta(\xi, \eta) = \arctan((\eta - 0.5)/(\xi - 0.7))$ and $r(\xi, \eta) = \sqrt{(\xi - 0.7)^2 + (\eta - 0.5)^2}$.

A monitor function close to the boundary can generate a negative value in the Jacobian function Table 4 and boundary errors Fig 6.

<table><caption>TABLE 4. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>4.816e-01</td><td>0.127</td><td>1.512e-01</td><td>-2.286e-02</td><td>3.587e+00</td></tr><tr><td>32</td><td>1.892e-01</td><td>0.303</td><td>1.201e-01</td><td>-7.120e-02</td><td>3.644e+00</td></tr><tr><td>64</td><td>4.322e-02</td><td>0.784</td><td>1.083e-01</td><td>1.936e-01</td><td>3.382e+00</td></tr><tr><td>128</td><td>9.048e-03</td><td>2.259</td><td>1.081e-01</td><td>3.464e-01</td><td>3.128e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>4.702e-01</td><td>0.190</td><td>1.414e-01</td><td>5.377e-02</td><td>4.142e+00</td></tr><tr><td>32</td><td>1.288e-01</td><td>0.461</td><td>1.103e-01</td><td>4.119e-02</td><td>3.384e+00</td></tr><tr><td>64</td><td>2.965e-02</td><td>0.978</td><td>1.081e-01</td><td>2.215e-01</td><td>3.300e+00</td></tr><tr><td>128</td><td>4.821e-03</td><td>3.194</td><td>1.081e-01</td><td>3.511e-01</td><td>3.136e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>4.640e-01</td><td>0.320</td><td>1.459e-01</td><td>-5.108e-02</td><td>5.528e+00</td></tr><tr><td>32</td><td>1.267e-01</td><td>0.549</td><td>1.099e-01</td><td>6.682e-02</td><td>3.461e+00</td></tr><tr><td>64</td><td>2.787e-02</td><td>1.355</td><td>1.081e-01</td><td>2.348e-01</td><td>3.302e+00</td></tr><tr><td>128</td><td>4.112e-03</td><td>4.430</td><td>1.081e-01</td><td>3.510e-01</td><td>3.140e+00</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>FIGURE 5. (left) Adaptive meshes and (right) the monitor function.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6. A close-up of the region near the boundary.</figcaption></figure>

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.

$$ (5.8) \qquad \rho(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 \exp(-100|(x + 0.45)^2 + (y - 0.4)^2 - 0.1|). $$
<table><caption>TABLE 5. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th><th>max<sub>∂Ω<sub>c</sub></sub>|(Ψ - x) ⋅ n̄|</th></tr></thead><tbody><tr><td>16</td><td>2.720e-01</td><td>0.070</td><td>5.179e-02</td><td>9.727e-02</td><td>1.719e+00</td><td>2.426e-03</td></tr><tr><td>32</td><td>7.321e-02</td><td>0.146</td><td>5.786e-02</td><td>9.899e-02</td><td>1.412e+00</td><td>3.969e-04</td></tr><tr><td>64</td><td>1.765e-02</td><td>0.414</td><td>5.807e-02</td><td>1.338e-01</td><td>1.314e+00</td><td>4.189e-05</td></tr><tr><td>128</td><td>5.066e-03</td><td>1.516</td><td>5.811e-02</td><td>1.365e-01</td><td>1.280e+00</td><td>2.881e-07</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th><th>max<sub>∂Ω<sub>c</sub></sub>|(Ψ - x) ⋅ n̄|</th></tr></thead><tbody><tr><td>16</td><td>2.353e-01</td><td>0.090</td><td>5.217e-02</td><td>1.835e-01</td><td>1.613e+00</td><td>1.848e-03</td></tr><tr><td>32</td><td>5.832e-02</td><td>0.222</td><td>5.818e-02</td><td>1.554e-01</td><td>1.429e+00</td><td>4.023e-04</td></tr><tr><td>64</td><td>1.375e-02</td><td>0.617</td><td>5.807e-02</td><td>1.373e-01</td><td>1.307e+00</td><td>4.301e-05</td></tr><tr><td>128</td><td>3.558e-03</td><td>2.158</td><td>5.811e-02</td><td>1.353e-01</td><td>1.280e+00</td><td>3.245e-07</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th><th>max<sub>∂Ω<sub>c</sub></sub>|(Ψ - x) ⋅ n̄|</th></tr></thead><tbody><tr><td>16</td><td>2.205e-01</td><td>0.123</td><td>5.482e-02</td><td>1.490e-01</td><td>1.621e+00</td><td>2.606e-03</td></tr><tr><td>32</td><td>4.607e-02</td><td>0.292</td><td>5.812e-02</td><td>1.487e-01</td><td>1.367e+00</td><td>2.842e-04</td></tr><tr><td>64</td><td>1.273e-02</td><td>0.858</td><td>5.807e-02</td><td>1.414e-01</td><td>1.306e+00</td><td>3.999e-05</td></tr><tr><td>128</td><td>3.496e-03</td><td>3.093</td><td>5.811e-02</td><td>1.354e-01</td><td>1.280e+00</td><td>5.209e-07</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>FIGURE 7. (left) Adaptive meshes and (right) the monitor function.</figcaption></figure>

**Example 5.6 (Circle).** In this example we use a non-axisymmetric monitor function defined by:

$$ (5.9) \qquad \rho(x, y) = 1 + 5/\operatorname{sech}\left(5\left(x - \frac{\sqrt{3}}{2}\right)^2 + \left(y - 0.5\right)^2 - \frac{\pi^2}{4}\right) \\ \qquad \qquad + 5/\operatorname{sech}\left(5\left(x + \frac{\sqrt{3}}{2}\right)^2 + \left(y - 0.5\right)^2 - \frac{\pi^2}{4}\right). $$
<figure><img src="image_9.png" /><figcaption>FIGURE 8. A close-up of the region where the two circles intersect (left) and on one of the singular points (right).</figcaption></figure>

*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.*

<caption>TABLE 6. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption>

(A) Degree $p=3$

<table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th><th>max|<span style="text-decoration: overline">(Ψ - x)</span> · <span style="text-decoration: overline>n</span>|</th></tr></thead><tbody><tr><td>16</td><td>7.015e-02</td><td>0.070</td><td>9.650e-02</td><td>1.742e-01</td><td>2.167e+00</td><td>1.453e-03</td></tr><tr><td>32</td><td>1.646e-02</td><td>0.137</td><td>9.627e-02</td><td>1.706e-01</td><td>1.926e+00</td><td>3.079e-04</td></tr><tr><td>64</td><td>3.468e-03</td><td>0.446</td><td>9.620e-02</td><td>1.699e-01</td><td>1.869e+00</td><td>4.823e-05</td></tr><tr><td>128</td><td>8.166e-04</td><td>1.695</td><td>9.619e-02</td><td>1.701e-01</td><td>1.869e+00</td><td>6.217e-06</td></tr></tbody></table>

(B) Degree $p=4$

<table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th><th>max|<span style="text-decoration: overline">(Ψ - x)</span> · <span style="text-decoration: overline>n</span>|</th></tr></thead><tbody><tr><td>16</td><td>6.839e-02</td><td>0.096</td><td>9.577e-02</td><td>1.671e-01</td><td>2.015e+00</td><td>1.665e-03</td></tr><tr><td>32</td><td>1.168e-02</td><td>0.214</td><td>9.619e-02</td><td>1.704e-01</td><td>1.899e+00</td><td>1.381e-04</td></tr><tr><td>64</td><td>9.126e-04</td><td>0.689</td><td>9.619e-02</td><td>1.701e-01</td><td>1.869e+00</td><td>1.590e-05</td></tr><tr><td>128</td><td>7.788e-05</td><td>2.609</td><td>9.619e-02</td><td>1.701e-01</td><td>1.868e+00</td><td>9.362e-07</td></tr></tbody></table>

(C) Degree $p=5$

<table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th><th>max|<span style="text-decoration: overline">(Ψ - x)</span> · <span style="text-decoration: overline>n</span>|</th></tr></thead><tbody><tr><td>16</td><td>4.822e-02</td><td>0.131</td><td>9.631e-02</td><td>1.660e-01</td><td>1.992e+00</td><td>9.295e-04</td></tr><tr><td>32</td><td>9.567e-03</td><td>0.284</td><td>9.620e-02</td><td>1.702e-01</td><td>1.888e+00</td><td>1.254e-04</td></tr><tr><td>64</td><td>6.496e-04</td><td>0.945</td><td>9.619e-02</td><td>1.701e-01</td><td>1.869e+00</td><td>4.599e-06</td></tr><tr><td>128</td><td>1.833e-05</td><td>3.512</td><td>9.619e-02</td><td>1.701e-01</td><td>1.868e+00</td><td>2.851e-07</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>FIGURE 9. (left) Adaptive meshes and (right) the monitor function.</figcaption></figure>
## 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 $\Psi(\partial\Omega_c) = \partial\Omega_c$. Hence, returning to the system 4.2, by introducing $\Psi = \nabla\phi$, we have

$$ (6.1) \qquad \left\{ \begin{array}{ll} \nabla \cdot \Psi = \mathcal{G}(\Psi) & \text{in } \Omega_c, \\ \Psi = \nabla\phi & \text{in } \Omega_c, \\ \Psi \cdot n = x \cdot n & \text{on } \partial\Omega_c, \end{array} \right. $$

where

$$ (6.2) \qquad \mathcal{G}(\Psi) = \left( (\nabla \cdot \Psi)^d + d! (f(\Psi) - \det (\nabla(\Psi))) \right)^{\frac{1}{d}}. $$

Hence, let's introduce the decomposition $\Psi = \Psi_0 + \Psi_N$ such that $\Psi_0 \cdot n = 0$ and $\Psi_N \cdot n = x \cdot n$. Therefore, the problem (6.1) writes

$$ (6.3) \qquad \left\{ \begin{array}{ll} \nabla \cdot \Psi_0 = \mathcal{G}(\Psi) - \nabla \cdot \Psi_N & \text{in } \Omega_c, \\ \Psi_0 = \nabla\phi - \Psi_N & \text{in } \Omega_c, \\ \Psi_0 \cdot n = 0 & \text{on } \partial\Omega_c. \end{array} \right. $$

Then, the variational formulation of 6.3, is:

Find $(\Psi, \phi) \in H_0(\text{div}, \Omega) \times L_0^2(\Omega)$ such that $\forall (v, \varphi) \in H_0(\text{div}, \Omega) \times L_0^2(\Omega)$

$$ (6.4) \qquad \left\{ \begin{array}{ll} \int_{\Omega} \Psi_0 \cdot v & + \int_{\Omega} \phi \nabla \cdot v = - \int_{\Omega} \Psi_N \cdot v, \\ \int_{\Omega} (\nabla \cdot \Psi_0) \varphi & = \int_{\Omega} \mathcal{G}(\Psi) \varphi - \int_{\Omega} (\nabla \cdot \Psi_N) \varphi, \end{array} \right. $$

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.

$$ (6.5) \qquad \left\{ \begin{array}{ll} -\nabla^2 \phi = f & \text{in } \Omega, \\ \nabla \phi \cdot n = 0 & \text{on } \partial\Omega, \end{array} \right. $$

As well known, the associated variational formulation can be written as a minimization problem over $H^1(\Omega)$. Let us now, introduce an auxiliary variable $u = \nabla\phi$. The Poisson problem is then equivalent to the following system

$$ (6.6) \qquad \left\{ \begin{array}{ll} -\nabla \cdot u = f & \text{in } \Omega, \\ u = \nabla\phi & \text{in } \Omega, \\ u \cdot n = 0 & \text{on } \partial\Omega. \end{array} \right. $$

The mixed weak formulation for the Poisson equation is

Find $(u, \phi) \in H_0(\text{div}, \Omega) \times L_0^2(\Omega)$ such that $\forall (\Psi, \varphi) \in H_0(\text{div}, \Omega) \times L_0^2(\Omega)$

$$ (6.7) \qquad \left\{ \begin{array}{ll} \int_{\Omega} u \cdot \Psi & + \int_{\Omega} \phi \nabla \cdot \Psi = 0, \\ \int_{\Omega} (\nabla \cdot u) \varphi & = - \int_{\Omega} f \varphi. \end{array} \right. $$

Using the discrete spaces defined in 2.13 and 2.14, we get the following variational problem
Find $(\mathbf{u}_h, \phi_h) \in \mathbf{V}_h(\mathrm{div}, \Omega) \times \mathbf{V}_h(L^2, \Omega)$ such that $\forall (\mathbf{\Psi}, \varphi) \in \mathbf{V}_h(\mathrm{div}, \Omega) \times \mathbf{V}_h(L^2, \Omega)$

$$ (6.8) \qquad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \mathbf{u}_h \cdot \mathbf{\Psi} \quad + \int_{\Omega} \phi_h \nabla \cdot \mathbf{\Psi} \quad = 0, \\ \displaystyle \int_{\Omega} (\nabla \cdot \mathbf{u}_h) \varphi \quad = - \int_{\Omega} f \varphi. \end{array} \right. $$

We then recover a saddle point problem of the form

$$ (6.9) \qquad \mathcal{A} \begin{pmatrix} \mathbf{u}_h \\ \phi_h \end{pmatrix} = \begin{pmatrix} A & B \\ B^T & 0 \end{pmatrix} \begin{pmatrix} \mathbf{u}_h \\ \phi_h \end{pmatrix} = \begin{pmatrix} 0 \\ F \end{pmatrix} $$

with

$$ F_i = - \int_{\Omega} f \phi_i \, dx, \quad 1 \le i \le 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.

$$ (6.10) \qquad \left\{ \begin{aligned} (M_s)_{i_s, j_s} &= \int_0^1 B_{\mu_s, i_s}^{p_s} B_{\mu_s, j_s}^{p_s} \, dx_s \\ (D_s)_{i_s, j_s} &= \int_0^1 B_{\mu_s-1, i_s}^{p_s-1} B_{\mu_s-1, j_s}^{p_s-1} \, dx_s \\ (R_s)_{i_s, j_s} &= \int_0^1 B_{\mu_s-1, i_s}^{p_s-1} \frac{d B_{\mu_s, j_s}^{p_s}}{d\xi} \, dx_s \end{aligned} \right. $$

where

$$ B_{\mu_s, i_s}^{p_s} = B_{i_s}^{p_s}, \text{ such that } B_{i_s}^{p_s} \in S_{\mu_s}^{p_s}, s = 1, 2, \dots, d. $$

**Matrix form in 2D.** For the sake of simplicity we shall introduce the auxiliary scalar functions

$$ (6.11) \qquad \left\{ \begin{aligned} \Psi_j^1 &= B_{\mu_1, j_1}^{p_1} B_{\mu_2-1, j_2}^{p_2-1} \\ \Psi_j^2 &= B_{\mu_1-1, j_1}^{p_1-1} B_{\mu_2, j_2}^{p_2} \\ \phi_j &= B_{\mu_1-1, j_1}^{p_1-1} B_{\mu_2-1, j_2}^{p_2-1} \end{aligned} \right. $$

where $j = (j_1, j_2)$, we also define the vectors $e_1 = \begin{pmatrix} 1 \\ 0 \end{pmatrix}$ and $e_2 = \begin{pmatrix} 0 \\ 1 \end{pmatrix}$. Therefore, the expression of $\mathbf{u}_h \in \mathbf{V}_h(\mathrm{div}, \Omega)$ becomes

$$ \mathbf{u}_h = \sum_j (u_j^1 \Psi_j^1 e_1 + u_j^2 \Psi_j^2 e_2) $$

Because $e_i \cdot e_j = \delta_{ij}$, we find that $A$ is a block diagonal matrix of the form

$$ (6.12) \qquad A = \begin{pmatrix} A_{11} & 0 \\ 0 & A_{22} \end{pmatrix} $$
where

$$A_{11i,j} = \int_{\Omega} \Psi_i^1 \Psi_j^1 \mathrm{d}\mathbf{x} = (M_1 \otimes D_2)_{ij}$$

$$A_{22i,j} = \int_{\Omega} \Psi_i^2 \Psi_j^2 \mathrm{d}\mathbf{x} = (D_1 \otimes M_2)_{ij}$$

While for $B$ we have

$$ (6.13) \qquad B = \begin{pmatrix} B_1 \\ B_2 \end{pmatrix} $$

where

$$B_{1i,j} = \int_{\Omega} \phi_j \partial_{x_1} \Psi_i^1 \mathrm{d}\mathbf{x} = (R_1 \otimes D_2)_{ij}$$

$$B_{2i,j} = \int_{\Omega} \phi_j \partial_{x_2} \Psi_i^2 \mathrm{d}\mathbf{x} = (D_1 \otimes R_2)_{ij}$$

The final matrix form is

$$ \mathcal{A} = \begin{pmatrix} M_1 \otimes D_2 & 0 & R_1 \otimes D_2 \\ 0 & D_1 \otimes M_2 & D_1 \otimes R_2 \\ (R_1 \otimes D_2)^T & (D_1 \otimes R_2)^T & 0 \end{pmatrix} $$

Matrix form in 3D. It is possible to extend the notation of the auxiliary scalar functions to the 3D case as follows:

$$\Psi_j^1 = B_{\mu_1, j_1}^{p_1} B_{\mu_2, j_2}^{p_2-1} B_{\mu_3, j_3}^{p_3-1}$$

$$\Psi_j^2 = B_{\mu_1-1, j_1}^{p_1-1} B_{\mu_2, j_2}^{p_2} B_{\mu_3-1, j_3}^{p_3-1}$$

$$\Psi_j^3 = B_{\mu_1-1, j_1}^{p_1-1} B_{\mu_2-1, j_2}^{p_2-1} B_{\mu_3, j_3}^{p_3}$$

where $j = (j_1, j_2, j_3)$, and

$$\phi_j = B_{\mu_1-1, j_1}^{p_1-1} B_{\mu_2-1, j_2}^{p_2-1} B_{\mu_3-1, j_3}^{p_3-1}$$

we also define the vectors $e_1 = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}$, $e_2 = \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}$ and $e_3 = \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}$. The expression of $u_h \in V_h(\mathrm{div}, \Omega)$ equivalent to

$$u_h = \sum_j (u_j^1 \Psi_j^1 e_1 + u_j^2 \Psi_j^2 e_2 + u_j^3 \Psi_j^3 e_3)$$

and as a consequence of $e_i \cdot e_j = \delta_{ij}$, we have

$$ (6.14) \qquad A = \begin{pmatrix} A_{11} & 0 & 0 \\ 0 & A_{22} & 0 \\ 0 & 0 & A_{33} \end{pmatrix} $$

where

$$A_{11i,j} = \int_{\Omega} \Psi_i^1 \Psi_j^1 \mathrm{d}\mathbf{x} = (M_1 \otimes D_2 \otimes D_3)_{ij}$$

$$A_{22i,j} = \int_{\Omega} \Psi_i^2 \Psi_j^2 \mathrm{d}\mathbf{x} = (D_1 \otimes M_2 \otimes D_3)_{ij}$$

$$A_{33i,j} = \int_{\Omega} \Psi_i^3 \Psi_j^3 \mathrm{d}\mathbf{x} = (D_1 \otimes D_2 \otimes M_3)_{ij}$$
While for $B$ we have

$$ (6.15) \qquad B = \begin{pmatrix} B_1 \\ B_2 \\ B_3 \end{pmatrix} $$

where

$$ B_{1i,j} = \int_{\Omega} \phi_j \partial_{x_1} \Psi_i^1 dx = (R_1 \otimes D_2 \otimes D_3)_{ij} $$

$$ B_{2i,j} = \int_{\Omega} \phi_j \partial_{x_2} \Psi_i^2 dx = (D_1 \otimes R_2 \otimes D_3)_{ij} $$

$$ B_{3i,j} = \int_{\Omega} \phi_j \partial_{x_3} \Psi_i^3 dx = (D_1 \otimes D_2 \otimes R_3)_{ij} $$

Finally,

$$ A = \begin{pmatrix} M_1 \otimes D_2 \otimes D_3 & 0 & 0 & R_1 \otimes D_2 \otimes D_3 \\ 0 & D_1 \otimes M_2 \otimes D_3 & 0 & D_1 \otimes R_2 \otimes D_3 \\ 0 & 0 & D_1 \otimes D_2 \otimes M_3 & D_1 \otimes D_2 \otimes R_3 \\ (R_1 \otimes D_2 \otimes D_3)^T & (D_1 \otimes R_2 \otimes D_3)^T & (D_1 \otimes D_2 \otimes R_3)^T & 0 \end{pmatrix} $$

6.2.1. Linear Solver. Based on the problem 6.9, we consider the following saddle point problem

$$ (6.16) \qquad A \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} A & B \\ B^T & 0 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \end{pmatrix} $$

Let us define $\tilde{A} = B^T A^{-1} B$ as the Schur complement of (6.16). Using the Uzawa method leads to

$$ (6.17) \qquad \begin{cases} \tilde{A}x_2 = B^T A^{-1}b_1 - b_2 \\ Ax_1 = b_1 - Bx_2 \end{cases} $$

Fast Diagonalization method (2D). The matrix $\tilde{A}$ can be expressed as the sum of Kronecker matrices

$$ (6.18) \qquad \tilde{A} = D_1 \otimes (R_2^T M_2^{-1} R_2) + (R_1^T M_1^{-1} R_1) \otimes D_2 $$
We have

$$
\begin{align*}
B^T A^{-1} b_1 &= ((R_1 \otimes D_2)^T \quad (D_1 \otimes R_2)^T) \begin{pmatrix} ((M_1 \otimes D_2)^{-1} & 0 \\ 0 & (D_1 \otimes M_2)^{-1} \end{pmatrix} \begin{pmatrix} b_{11} \\ b_{12} \end{pmatrix} \\
&= (R_1^T \otimes D_2^T \quad D_1^T \otimes R_2^T) \begin{pmatrix} M_1^{-1} \otimes D_2^{-1} & 0 \\ 0 & D_1^{-1} \otimes M_2^{-1} \end{pmatrix} \begin{pmatrix} b_{11} \\ b_{12} \end{pmatrix} \\
&= ((R_1^T \otimes D_2^T) (M_1^{-1} \otimes D_2^{-1}) \quad (D_1^T \otimes R_2^T) (D_1^{-1} \otimes M_2^{-1})) \begin{pmatrix} b_{11} \\ b_{12} \end{pmatrix} \\
&= ((R_1^T M_1^{-1}) \otimes I_2 \quad I_1 \otimes (R_2^T M_2^{-1})) \begin{pmatrix} b_{11} \\ b_{12} \end{pmatrix} \\
&= ((R_1^T M_1^{-1}) \otimes I_2) b_{11} + (I_1 \otimes (R_2^T M_2^{-1})) b_{12}
\end{align*}
$$

Therefor

$$
(6.19) \qquad \tilde{A}x_2 = b
$$

where

$$
b = ((R_1^T M_1^{-1}) \otimes I_2) b_{11} + (I_1 \otimes (R_2^T M_2^{-1})) b_{12} - b_2
$$

On the other hand, we have

$$
(6.20) \qquad Bx_2 = \begin{pmatrix} (R_1 \otimes D_2) x_2 \\ (D_1 \otimes R_2) x_2 \end{pmatrix}
$$

By introducing $x_1 := (x_{11}, x_{12})^T$, the equation $Ax_1 = b_1 - Bx_2$ is equivalent to

$$
\begin{cases}
x_{11} = (M_1^{-1} \otimes D_2^{-1}) b_{11} - (M_1^{-1} \otimes D_2^{-1}) (R_1 \otimes D_2) x_2 \\
x_{12} = (D_1^{-1} \otimes M_2^{-1}) b_{12} - (D_1^{-1} \otimes M_2^{-1}) (D_1 \otimes R_2) x_2
\end{cases}
$$

which can be simplified as

$$
\begin{cases}
x_{11} = (M_1^{-1} \otimes D_2^{-1}) b_{11} - (M_1^{-1} R_1 \otimes I_2) x_2 \\
x_{12} = (D_1^{-1} \otimes M_2^{-1}) b_{12} - (I_1 \otimes M_2^{-1} R_2) x_2
\end{cases}
$$

In order to devise a fast solver (6.17), we choose to follow the approach of Sangalli and Tani [48]. We describe in the sequel the fast diagonalization method in the case of Isogeometric Analysis. This method was first introduced in [45]. We consider the generalized eigendecompositions problems

$$
(6.21) \qquad (R_1^T M_1^{-1} R_1) U_1 = D_1 U_1 d_1, \quad (R_2^T M_2^{-1} 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) \qquad U_1^T D_1 U_1 = I_1, \quad U_2^T D_2 U_2 = I_2
$$

Therefor,

$$
(6.23) \qquad (U_1 \otimes U_2)^{-T} (d_1 \otimes I_2 + I_1 \otimes d_2) (U_1 \otimes 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):

<caption>**Algorithm 2: fast_diag:** Fast diagonalization method for (6.17) in the 2D case</caption>

Input : b
Output: x
1 (b<sub>11</sub>, b<sub>12</sub>), b<sub>2</sub> ← unfold(b)
2 r<sub>2</sub> ← ((R<sub>1</sub><sup>T</sup>M<sub>1</sub><sup>-1</sup>) ⊗ I<sub>2</sub>) b<sub>11</sub> + (I<sub>1</sub> ⊗ (R<sub>2</sub><sup>T</sup>M<sub>2</sub><sup>-1</sup>)) b<sub>12</sub> - b<sub>2</sub>
3 r̃<sub>2</sub> ← (U<sub>1</sub> ⊗ U<sub>2</sub>)<sup>T</sup> r<sub>2</sub>
4 x̃<sub>2</sub> ← (d<sub>1</sub> ⊗ I<sub>2</sub> + I<sub>1</sub> ⊗ d<sub>2</sub>)<sup>-1</sup> r̃<sub>2</sub>
5 x<sub>2</sub> ← (U<sub>1</sub> ⊗ U<sub>2</sub>) x̃<sub>2</sub>
6 x<sub>11</sub> ← (M<sub>1</sub><sup>-1</sup> ⊗ D<sub>2</sub><sup>-1</sup>) b<sub>11</sub> - (M<sub>1</sub><sup>-1</sup>R<sub>1</sub> ⊗ I<sub>2</sub>) x<sub>2</sub>
7 x<sub>12</sub> ← (D<sub>1</sub><sup>-1</sup> ⊗ M<sub>2</sub><sup>-1</sup>) b<sub>12</sub> - (I<sub>1</sub> ⊗ M<sub>2</sub><sup>-1</sup>R<sub>2</sub>) x<sub>2</sub>
8 x ← fold(x<sub>11</sub>, x<sub>12</sub>, x<sub>2</sub>)

Fast Diagonalization method (3D). The matrix $\tilde{A}$ in (3D) can be expressed as the sum of Kronecker matrices

$$
(6.24) \qquad \tilde{A} = D_1 \otimes D_2 \otimes (R_3^T M_3^{-1} R_3) + D_1 \otimes (R_2^T M_2^{-1} R_2) \otimes D_3 + (R_1^T M_1^{-1} R_1) \otimes D_2 \otimes D_3
$$

And as in (2D) case we can show that we have

$$
B^T A^{-1} b_1 = ((R_1^T M_1^{-1}) \otimes I_2 \otimes I_3) b_{11} + (I_1 \otimes (R_2^T M_2^{-1}) \otimes I_3) b_{12} + (I_1 \otimes I_2 \otimes (R_3^T M_3^{-1})) b_{13}.
$$

Therefor

$$
(6.25) \qquad \tilde{A}x_2 = ((R_1^T M_1^{-1}) \otimes I_2 \otimes I_3) b_{11} + (I_1 \otimes (R_2^T M_2^{-1}) \otimes I_3) b_{12} + (I_1 \otimes I_2 \otimes (R_3^T M_3^{-1})) b_{13} - b_2.
$$

On the other hand, we have

$$
(6.26) \qquad Bx_2 = \begin{pmatrix} (R_1 \otimes D_2 \otimes D_3) x_2 \\ (D_1 \otimes R_2 \otimes D_3) x_2 \\ (D_1 \otimes D_2 \otimes R_3) x_2 \end{pmatrix}
$$

By introducing $x_1 := (x_{11}, x_{12}, x_{13})^T$, the equation $Ax_1 = b_1 - Bx_2$ is equivalent to

$$
\begin{cases}
x_{11} = (M_1^{-1} \otimes D_2^{-1} \otimes D_3^{-1}) b_{11} - (M_1^{-1} \otimes D_2^{-1} \otimes D_3^{-1}) (R_1 \otimes D_2 \otimes D_3) x_2 \\
x_{12} = (D_1^{-1} \otimes M_2^{-1} \otimes D_3^{-1}) b_{12} - (D_1^{-1} \otimes M_2^{-1} \otimes D_3^{-1}) (D_1 \otimes R_2 \otimes D_3) x_2 \\
x_{13} = (D_1^{-1} \otimes D_2^{-1} \otimes M_3^{-1}) b_{13} - (D_1^{-1} \otimes D_2^{-1} \otimes M_3^{-1}) (D_1 \otimes D_2 \otimes R_3) x_2
\end{cases}
$$

which can be simplified to

$$
\left\{
\begin{aligned}
x_{11} &= (M_1^{-1} \otimes D_2^{-1} \otimes D_3^{-1}) b_{11} - (M_1^{-1} R_1 \otimes I_2 \otimes I_3) x_2 \\
x_{12} &= (D_1^{-1} \otimes M_2^{-1} \otimes D_3^{-1}) b_{12} - (I_1 \otimes M_2^{-1} R_2 \otimes I_3) x_2 \\
x_{13} &= (D_1^{-1} \otimes D_2^{-1} \otimes M_3^{-1}) b_{13} - (I_1 \otimes I_2 \otimes M_3^{-1} R_3) x_2
\end{aligned}
\right.
$$
Following the same steps, we consider the generalized eigendecompositions problems

(6.27)

$$ (R_1^T M_1^{-1} R_1) U_1 = D_1 U_1 d_1, \quad (R_2^T M_2^{-1} R_2) U_2 = D_2 U_2 d_2, \quad (R_3^T M_3^{-1} 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_3$ fulfills

$$ (6.28) \qquad U_1^T D_1 U_1 = I_1, \quad U_2^T D_2 U_2 = I_2, \quad U_3^T D_3 U_3 = I_3 $$

Therefor

(6.29)

$$ (U_1 \otimes U_2 \otimes U_3)^{-T} (d_1 \otimes I_2 \otimes I_3 + I_1 \otimes d_2 \otimes I_3 + I_1 \otimes I_2 \otimes d_3) (U_1 \otimes U_2 \otimes U_3)^{-1} = b $$

The direct solver for (6.17), is then given by the following algorithm:

<table><caption>**Algorithm 3: fast_diag**: Fast diagonalization method for (6.17) in the 3D case</caption><thead><tr><th></th><th>Input : b</th></tr></thead><tbody><tr><td>1</td><td>b<sub>11</sub>, b<sub>12</sub>, b<sub>13</sub>) , b<sub>2</sub> &larr; unfold(b)</td></tr><tr><td>2</td><td>r<sub>2</sub> &larr; ((R<sub>1</sub><sup>T</sup>M<sub>1</sub><sup>-1</sup>) &otimes; I<sub>2</sub> &otimes; I<sub>3</sub>) b<sub>11</sub> + (I<sub>1</sub> &otimes; (R<sub>2</sub><sup>T</sup>M<sub>2</sub><sup>-1</sup>) &otimes; I<sub>3</sub>) b<sub>12</sub> + (I<sub>1</sub> &otimes; I<sub>2</sub> &otimes; (R<sub>3</sub><sup>T</sup>M<sub>3</sub><sup>-1</sup>)) b<sub>13</sub> - b<sub>2</sub></td></tr><tr><td>3</td><td>r&#x0303;<sub>2</sub> &larr; (U<sub>1</sub> &otimes; U<sub>2</sub> &otimes; U<sub>3</sub>)<sup>T</sup> r<sub>2</sub></td></tr><tr><td>4</td><td>x&#x0303;<sub>2</sub> &larr; (d<sub>1</sub> &otimes; I<sub>2</sub> &otimes; I<sub>3</sub> + I<sub>1</sub> &otimes; d<sub>2</sub> &otimes; I<sub>3</sub> + I<sub>1</sub> &otimes; I<sub>2</sub> &otimes; d<sub>3</sub>)<sup>-1</sup> r&#x0303;<sub>2</sub></td></tr><tr><td>5</td><td>x<sub>2</sub> &larr; (U<sub>1</sub> &otimes; U<sub>2</sub> &otimes; U<sub>3</sub>) x&#x0303;<sub>2</sub></td></tr><tr><td>6</td><td>x<sub>11</sub> &larr; (M<sub>1</sub><sup>-1</sup> &otimes; D<sub>2</sub><sup>-1</sup> &otimes; D<sub>3</sub><sup>-1</sup>) b<sub>11</sub> - (M<sub>1</sub><sup>-1</sup>R<sub>1</sub> &otimes; I<sub>2</sub> &otimes; I<sub>3</sub>) x<sub>2</sub></td></tr><tr><td>7</td><td>x<sub>12</sub> &larr; (D<sub>1</sub><sup>-1</sup> &otimes; M<sub>2</sub><sup>-1</sup> &otimes; D<sub>3</sub><sup>-1</sup>) b<sub>12</sub> - (I<sub>1</sub> &otimes; M<sub>2</sub><sup>-1</sup>R<sub>2</sub> &otimes; I<sub>3</sub>) x<sub>2</sub></td></tr><tr><td>8</td><td>x<sub>13</sub> &larr; (D<sub>1</sub><sup>-1</sup> &otimes; D<sub>2</sub><sup>-1</sup> &otimes; M<sub>3</sub><sup>-1</sup>) b<sub>13</sub> - (I<sub>1</sub> &otimes; I<sub>2</sub> &otimes; M<sub>3</sub><sup>-1</sup>R<sub>3</sub>) x<sub>2</sub></td></tr><tr><td>9</td><td>x &larr; fold(x<sub>11</sub>, x<sub>12</sub>, x<sub>13</sub>, x<sub>2</sub>)</td></tr></tbody></table>
**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.
## 7. 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><caption>TABLE 7. CPU-time needed to solve the Monge-Ampère equation and error evolution for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>H<sup>1</sup>-norm</th><th>H<sup>1</sup>-norm MG</th><th>CPU-time (s)</th><th>CPU-time (s) MG</th><th>Nbr-iter</th><th>Nbr-iter MG</th></tr></thead><tbody><tr><td>8</td><td>5.232e-05</td><td>5.232e-05</td><td>0.029</td><td>0.036</td><td>41</td><td>32</td></tr><tr><td>16</td><td>6.508e-06</td><td>6.507e-06</td><td>0.061</td><td>0.067</td><td>41</td><td>35</td></tr><tr><td>32</td><td>8.155e-07</td><td>8.153e-07</td><td>0.178</td><td>0.200</td><td>35</td><td>32</td></tr><tr><td>64</td><td>1.018e-07</td><td>1.017e-07</td><td>0.602</td><td>0.625</td><td>35</td><td>31</td></tr><tr><td>128</td><td>1.279e-08</td><td>1.273e-08</td><td>2.434</td><td>2.343</td><td>35</td><td>30</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>H<sup>1</sup>-norm</th><th>H<sup>1</sup>-norm MG</th><th>CPU-time (s)</th><th>CPU-time (s) MG</th><th>Nbr-iter</th><th>Nbr-iter MG</th></tr></thead><tbody><tr><td>8</td><td>1.786e-06</td><td>1.786e-06</td><td>0.032</td><td>0.042</td><td>38</td><td>29</td></tr><tr><td>16</td><td>1.150e-07</td><td>1.150e-07</td><td>0.069</td><td>0.068</td><td>35</td><td>26</td></tr><tr><td>32</td><td>7.383e-09</td><td>7.369e-09</td><td>0.222</td><td>0.204</td><td>35</td><td>24</td></tr><tr><td>64</td><td>8.965e-10</td><td>9.467e-10</td><td>0.763</td><td>0.598</td><td>35</td><td>22</td></tr><tr><td>128</td><td>7.297e-10</td><td>8.337e-10</td><td>2.937</td><td>2.074</td><td>35</td><td>20</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>H<sup>1</sup>-norm</th><th>H<sup>1</sup>-norm MG</th><th>CPU-time (s)</th><th>CPU-time (s) MG</th><th>Nbr-iter</th><th>Nbr-iter MG</th></tr></thead><tbody><tr><td>8</td><td>8.527e-08</td><td>8.527e-08</td><td>0.049</td><td>0.047</td><td>35</td><td>26</td></tr><tr><td>16</td><td>3.057e-09</td><td>3.056e-09</td><td>0.085</td><td>0.085</td><td>35</td><td>26</td></tr><tr><td>32</td><td>7.407e-10</td><td>7.475e-10</td><td>0.281</td><td>0.254</td><td>35</td><td>24</td></tr><tr><td>64</td><td>7.293e-10</td><td>7.368e-10</td><td>0.983</td><td>0.757</td><td>35</td><td>22</td></tr><tr><td>128</td><td>7.292e-10</td><td>7.368e-10</td><td>3.846</td><td>2.641</td><td>35</td><td>20</td></tr></tbody></table>

**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><caption>TABLE 8. CPU-time needed to solve the Monge-Ampere equation and error evolution for various polynomial degrees *p* along each direction.</caption><thead><tr><th>#cells</th><th>H<sup>1</sup>-norm</th><th>H<sup>1</sup>-norm MG</th><th>CPU-time (s)</th><th>CPU-time (s) MG</th><th>Nbr-iter</th><th>Nbr-iter MG</th></tr></thead><tbody><tr><td>8</td><td>1.966e-05</td><td>1.966e-05</td><td>0.238</td><td>0.268</td><td>14</td><td>10</td></tr><tr><td>16</td><td>2.519e-06</td><td>2.519e-06</td><td>1.398</td><td>1.309</td><td>14</td><td>13</td></tr><tr><td>32</td><td>3.190e-07</td><td>3.190e-07</td><td>9.747</td><td>8.134</td><td>14</td><td>12</td></tr><tr><td>64</td><td>4.011e-08</td><td>4.016e-08</td><td>74.71</td><td>57.22</td><td>14</td><td>11</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>H<sup>1</sup>-norm</th><th>H<sup>1</sup>-norm MG</th><th>CPU-time (s)</th><th>CPU-time (s) MG</th><th>Nbr-iter</th><th>Nbr-iter MG</th></tr></thead><tbody><tr><td>8</td><td>6.622e-07</td><td>6.622e-07</td><td>0.521</td><td>0.583</td><td>14</td><td>10</td></tr><tr><td>16</td><td>4.322e-08</td><td>4.321e-08</td><td>3.686</td><td>2.702</td><td>14</td><td>10</td></tr><tr><td>32</td><td>2.750e-09</td><td>2.754e-09</td><td>27.55</td><td>15.50</td><td>14</td><td>8</td></tr><tr><td>64</td><td>1.965e-10</td><td>4.997e-10</td><td>214.5</td><td>90.37</td><td>14</td><td>6</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>H<sup>1</sup>-norm</th><th>H<sup>1</sup>-norm MG</th><th>CPU-time (s)</th><th>CPU-time (s) MG</th><th>Nbr-iter</th><th>Nbr-iter MG</th></tr></thead><tbody><tr><td>8</td><td>2.294e-08</td><td>2.294e-08</td><td>1.300</td><td>1.312</td><td>14</td><td>10</td></tr><tr><td>16</td><td>8.125e-10</td><td>8.118e-10</td><td>8.569</td><td>6.011</td><td>14</td><td>10</td></tr><tr><td>32</td><td>1.199e-10</td><td>1.202e-10</td><td>65.68</td><td>35.52</td><td>14</td><td>8</td></tr><tr><td>64</td><td>1.104e-10</td><td>1.112e-10</td><td>507.8</td><td>210.6</td><td>14</td><td>6</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>FIGURE 10. CPU time profiles between standard method (H<sup>2</sup>) and mixed variational formulation (Hdiv).</figcaption></figure>

## 7.2. 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) \qquad \rho(\xi, \eta) = \frac{1}{2 + \cos(8\pi\sqrt{(\xi - 0.5)^2 + (\eta - 0.5)^2})} $$

Fig. 11 shows the adaptive meshes results and Table 9 shows a high quality of mesh adaptation.

<table><caption>TABLE 9. Grid convergence analysis for various polynomial degrees *p* along each direction.</caption><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>4.949e-02</td><td>0.086</td><td>1.905e-02</td><td>5.017e-01</td><td>1.769e+00</td></tr><tr><td>32</td><td>4.665e-03</td><td>0.142</td><td>1.902e-02</td><td>5.681e-01</td><td>1.722e+00</td></tr><tr><td>64</td><td>5.005e-04</td><td>0.430</td><td>1.902e-02</td><td>5.697e-01</td><td>1.713e+00</td></tr><tr><td>128</td><td>1.128e-04</td><td>1.516</td><td>1.902e-02</td><td>5.699e-01</td><td>1.710e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>3.221e-02</td><td>0.101</td><td>1.903e-02</td><td>4.997e-01</td><td>1.748e+00</td></tr><tr><td>32</td><td>2.744e-03</td><td>0.174</td><td>1.902e-02</td><td>5.679e-01</td><td>1.718e+00</td></tr><tr><td>64</td><td>8.635e-05</td><td>0.534</td><td>1.902e-02</td><td>5.699e-01</td><td>1.710e+00</td></tr><tr><td>128</td><td>2.102e-05</td><td>1.879</td><td>1.902e-02</td><td>5.699e-01</td><td>1.710e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>4.077e-02</td><td>0.161</td><td>1.903e-02</td><td>5.005e-01</td><td>1.721e+00</td></tr><tr><td>32</td><td>2.289e-03</td><td>0.206</td><td>1.902e-02</td><td>5.673e-01</td><td>1.713e+00</td></tr><tr><td>64</td><td>3.299e-05</td><td>0.650</td><td>1.902e-02</td><td>5.699e-01</td><td>1.710e+00</td></tr><tr><td>128</td><td>3.737e-05</td><td>2.168</td><td>1.902e-02</td><td>5.699e-01</td><td>1.710e+00</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>FIGURE 11. (left) Adaptive meshes and (right) the monitor function.</figcaption></figure>
**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><caption>TABLE 10. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>4.165e-01</td><td>0.196</td><td>1.428e-01</td><td>1.077e-01</td><td>2.950e+00</td></tr><tr><td>32</td><td>1.231e-01</td><td>0.600</td><td>1.111e-01</td><td>3.123e-02</td><td>3.413e+00</td></tr><tr><td>64</td><td>3.020e-02</td><td>1.279</td><td>1.082e-01</td><td>2.209e-01</td><td>3.323e+00</td></tr><tr><td>128</td><td>4.589e-03</td><td>3.782</td><td>1.081e-01</td><td>3.509e-01</td><td>3.132e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>4.249e-01</td><td>0.216</td><td>1.435e-01</td><td>-3.402e-03</td><td>2.963e+00</td></tr><tr><td>32</td><td>1.182e-01</td><td>0.620</td><td>1.101e-01</td><td>8.919e-02</td><td>3.567e+00</td></tr><tr><td>64</td><td>2.747e-02</td><td>1.539</td><td>1.081e-01</td><td>2.483e-01</td><td>3.324e+00</td></tr><tr><td>128</td><td>3.411e-03</td><td>4.250</td><td>1.081e-01</td><td>3.513e-01</td><td>3.134e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>3.764e-01</td><td>0.344</td><td>1.407e-01</td><td>-6.914e-02</td><td>3.361e+00</td></tr><tr><td>32</td><td>1.065e-01</td><td>0.764</td><td>1.092e-01</td><td>8.918e-02</td><td>3.269e+00</td></tr><tr><td>64</td><td>2.213e-02</td><td>1.832</td><td>1.081e-01</td><td>2.296e-01</td><td>3.308e+00</td></tr><tr><td>128</td><td>3.138e-03</td><td>5.358</td><td>1.081e-01</td><td>3.511e-01</td><td>3.137e+00</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>FIGURE 12. (left) Adaptive meshes and (right) the monitor function.</figcaption></figure>
**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 table 11, with no self-intersection Figs. 13-14.

$$ \rho(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|) \\ (7.2) \\ + 5 \exp(-100|(x + 0.45)^2 + (y - 0.4)^2 - 0.1|) + 7 \exp(-100|x^2 + (y + 1.25)^2 - 0.4|). $$

<table><caption>TABLE 11. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>3.258e-01</td><td>0.119</td><td>4.369e-02</td><td>1.831e-01</td><td>1.574e+00</td></tr><tr><td>32</td><td>8.837e-02</td><td>0.373</td><td>4.479e-02</td><td>4.882e-02</td><td>1.508e+00</td></tr><tr><td>64</td><td>2.684e-02</td><td>0.785</td><td>4.392e-02</td><td>1.274e-01</td><td>1.408e+00</td></tr><tr><td>128</td><td>5.374e-03</td><td>2.342</td><td>4.389e-02</td><td>1.358e-01</td><td>1.333e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>3.325e-01</td><td>0.142</td><td>4.553e-02</td><td>1.556e-01</td><td>1.658e+00</td></tr><tr><td>32</td><td>7.359e-02</td><td>0.419</td><td>4.440e-02</td><td>9.047e-02</td><td>1.504e+00</td></tr><tr><td>64</td><td>2.246e-02</td><td>0.930</td><td>4.390e-02</td><td>1.420e-01</td><td>1.430e+00</td></tr><tr><td>128</td><td>5.135e-03</td><td>2.777</td><td>4.389e-02</td><td>1.362e-01</td><td>1.338e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>2.955e-01</td><td>0.176</td><td>4.407e-02</td><td>1.692e-01</td><td>1.613e+00</td></tr><tr><td>32</td><td>7.086e-02</td><td>0.610</td><td>4.437e-02</td><td>7.912e-02</td><td>1.492e+00</td></tr><tr><td>64</td><td>2.491e-02</td><td>1.108</td><td>4.392e-02</td><td>1.333e-01</td><td>1.405e+00</td></tr><tr><td>128</td><td>4.987e-03</td><td>3.461</td><td>4.389e-02</td><td>1.359e-01</td><td>1.341e+00</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>FIGURE 13. (left) Adaptive meshes and (right) the monitor function.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>FIGURE 14. A close-up of the region where the two circles intersect (left) and on one of the singular points (right).</figcaption></figure>

**Example 7.6 (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 12 shows fast convergence, high accuracy and positive Jacobian function.

<table><caption>TABLE 12. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>5.900e-02</td><td>0.114</td><td>9.608e-02</td><td>1.665e-01</td><td>1.982e+00</td></tr><tr><td>32</td><td>1.096e-02</td><td>0.226</td><td>9.611e-02</td><td>1.702e-01</td><td>1.893e+00</td></tr><tr><td>64</td><td>2.048e-03</td><td>0.671</td><td>9.611e-02</td><td>1.699e-01</td><td>1.867e+00</td></tr><tr><td>128</td><td>1.819e-03</td><td>2.490</td><td>9.611e-02</td><td>1.699e-01</td><td>1.867e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>4.695e-02</td><td>0.138</td><td>9.629e-02</td><td>1.676e-01</td><td>1.955e+00</td></tr><tr><td>32</td><td>8.832e-03</td><td>0.265</td><td>9.611e-02</td><td>1.703e-01</td><td>1.882e+00</td></tr><tr><td>64</td><td>1.922e-03</td><td>0.774</td><td>9.611e-02</td><td>1.699e-01</td><td>1.867e+00</td></tr><tr><td>128</td><td>1.816e-03</td><td>2.819</td><td>9.611e-02</td><td>1.699e-01</td><td>1.867e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>4.790e-02</td><td>0.204</td><td>9.617e-02</td><td>1.662e-01</td><td>1.942e+00</td></tr><tr><td>32</td><td>9.055e-03</td><td>0.366</td><td>9.611e-02</td><td>1.696e-01</td><td>1.882e+00</td></tr><tr><td>64</td><td>1.867e-03</td><td>1.057</td><td>9.611e-02</td><td>1.699e-01</td><td>1.867e+00</td></tr><tr><td>128</td><td>1.815e-03</td><td>3.612</td><td>9.611e-02</td><td>1.699e-01</td><td>1.867e+00</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>FIGURE 15. (left) Adaptive meshes and (right) the monitor function.</figcaption></figure>

**Example 7.7** (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.

$$ (7.3) \qquad \rho(x, y) = 1 + 5 \exp \left( -50|y - 0.5 - 0.25 \sin(2\pi x) \sin(\frac{2\pi}{3})| \right) + 5 \exp \left( -50|y - 0.25 - 0.25 \sin(2\pi x)| \right). $$

<table><caption>TABLE 13. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption><thead><tr><th colspan="7">(A) Degree p = 3</th></tr><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>6.810e-02</td><td>0.117</td><td>7.786e-02</td><td>1.880e-01</td><td>1.795e+00</td><td>1.795e+00</td></tr><tr><td>32</td><td>1.765e-02</td><td>0.230</td><td>7.752e-02</td><td>2.531e-01</td><td>1.678e+00</td><td>1.678e+00</td></tr><tr><td>64</td><td>4.754e-03</td><td>0.682</td><td>7.754e-02</td><td>2.740e-01</td><td>1.614e+00</td><td>1.614e+00</td></tr><tr><td>128</td><td>1.763e-03</td><td>2.397</td><td>7.754e-02</td><td>2.715e-01</td><td>1.614e+00</td><td>1.614e+00</td></tr></tbody><thead><tr><th colspan="7">(B) Degree p = 4</th></tr><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>5.320e-02</td><td>0.145</td><td>7.794e-02</td><td>1.696e-01</td><td>1.702e+00</td><td>1.702e+00</td></tr><tr><td>32</td><td>1.721e-02</td><td>0.278</td><td>7.751e-02</td><td>2.570e-01</td><td>1.660e+00</td><td>1.660e+00</td></tr><tr><td>64</td><td>4.538e-03</td><td>0.821</td><td>7.754e-02</td><td>2.739e-01</td><td>1.616e+00</td><td>1.616e+00</td></tr><tr><td>128</td><td>1.707e-03</td><td>2.806</td><td>7.754e-02</td><td>2.718e-01</td><td>1.614e+00</td><td>1.614e+00</td></tr></tbody><thead><tr><th colspan="7">(C) Degree p = 5</th></tr><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>5.087e-02</td><td>0.174</td><td>7.779e-02</td><td>1.737e-01</td><td>1.784e+00</td><td>1.784e+00</td></tr><tr><td>32</td><td>1.612e-02</td><td>0.367</td><td>7.751e-02</td><td>2.660e-01</td><td>1.670e+00</td><td>1.670e+00</td></tr><tr><td>64</td><td>4.386e-03</td><td>1.027</td><td>7.754e-02</td><td>2.739e-01</td><td>1.614e+00</td><td>1.614e+00</td></tr><tr><td>128</td><td>1.678e-03</td><td>3.737</td><td>7.754e-02</td><td>2.716e-01</td><td>1.614e+00</td><td>1.614e+00</td></tr></tbody></table>
<figure><img src="image_17.png" /><figcaption>FIGURE 16. (left) Adaptive meshes and (right) the monitor function.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 17. A close-up of the region where the mesh is very fine.</figcaption></figure>

*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) \qquad \rho(\xi, \eta, \gamma) = \frac{1}{2 + \cos(8\pi\sqrt{(x-0.5)^2 + (y-0.5)^2 + (z-0.5)^2})} $$
<table><caption>TABLE 14. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>3.626e-01</td><td>0.104</td><td>7.822e-03</td><td>3.264e-01</td><td>1.652e+00</td></tr><tr><td>32</td><td>2.698e-01</td><td>0.275</td><td>8.649e-03</td><td>3.806e-01</td><td>1.702e+00</td></tr><tr><td>64</td><td>7.026e-02</td><td>0.845</td><td>6.952e-03</td><td>4.616e-01</td><td>1.780e+00</td></tr><tr><td>128</td><td>8.173e-03</td><td>2.870</td><td>6.952e-03</td><td>5.763e-01</td><td>1.759e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>3.560e-01</td><td>0.133</td><td>8.141e-03</td><td>3.137e-01</td><td>1.975e+00</td></tr><tr><td>32</td><td>2.700e-01</td><td>0.423</td><td>8.663e-03</td><td>3.650e-01</td><td>1.786e+00</td></tr><tr><td>64</td><td>7.102e-02</td><td>1.160</td><td>6.953e-03</td><td>4.550e-01</td><td>1.755e+00</td></tr><tr><td>128</td><td>6.492e-03</td><td>3.619</td><td>6.952e-03</td><td>5.712e-01</td><td>1.749e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>3.599e-01</td><td>0.159</td><td>7.748e-03</td><td>2.696e-01</td><td>1.733e+00</td></tr><tr><td>32</td><td>2.713e-01</td><td>0.534</td><td>8.511e-03</td><td>3.011e-01</td><td>1.791e+00</td></tr><tr><td>64</td><td>7.167e-02</td><td>1.480</td><td>6.953e-03</td><td>4.484e-01</td><td>1.754e+00</td></tr><tr><td>128</td><td>5.981e-03</td><td>4.582</td><td>6.952e-03</td><td>5.696e-01</td><td>1.742e+00</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>FIGURE 18. Adaptive meshes and the monitor function.</figcaption></figure>

**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 [46] into the Möbius strip as*
follows.

$$ (7.5) \qquad \rho(x, y) = 1 + 5/\operatorname{sech}\left(5\left((x - \frac{\sqrt{3}}{2})^2 + (y - 0.5)^2 + (z - 0.5)^2 - \frac{\pi^2}{4}\right)\right) \\ \qquad \qquad + 5/\operatorname{sech}\left(5\left((x + \frac{\sqrt{3}}{2})^2 + (y - 0.5)^2 + (z - 0.5)^2 - \frac{\pi^2}{4}\right)\right). $$

<table><caption>TABLE 15. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>2.202e-01</td><td>0.236</td><td>9.284e-02</td><td>4.959e-02</td><td>2.023e+00</td></tr><tr><td>32</td><td>7.748e-02</td><td>0.332</td><td>8.802e-02</td><td>1.635e-01</td><td>1.944e+00</td></tr><tr><td>64</td><td>1.367e-02</td><td>0.915</td><td>8.745e-02</td><td>1.627e-01</td><td>1.853e+00</td></tr><tr><td>128</td><td>1.149e-03</td><td>3.073</td><td>8.743e-02</td><td>1.626e-01</td><td>1.790e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>2.173e-01</td><td>0.175</td><td>9.307e-02</td><td>7.994e-02</td><td>2.169e+00</td></tr><tr><td>32</td><td>8.067e-02</td><td>0.420</td><td>8.812e-02</td><td>1.576e-01</td><td>1.954e+00</td></tr><tr><td>64</td><td>1.486e-02</td><td>1.181</td><td>8.745e-02</td><td>1.627e-01</td><td>1.856e+00</td></tr><tr><td>128</td><td>1.057e-03</td><td>3.552</td><td>8.743e-02</td><td>1.626e-01</td><td>1.789e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>2.082e-01</td><td>0.268</td><td>9.200e-02</td><td>2.356e-02</td><td>2.060e+00</td></tr><tr><td>32</td><td>7.130e-02</td><td>0.498</td><td>8.794e-02</td><td>1.599e-01</td><td>1.949e+00</td></tr><tr><td>64</td><td>1.312e-02</td><td>1.517</td><td>8.745e-02</td><td>1.628e-01</td><td>1.841e+00</td></tr><tr><td>128</td><td>8.070e-04</td><td>4.622</td><td>8.743e-02</td><td>1.626e-01</td><td>1.792e+00</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>FIGURE 19. Adaptive meshes and the monitor function.</figcaption></figure>
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) \qquad \rho(\xi, \eta, \gamma) = \frac{1}{2 + \cos(8\pi\sqrt{(\xi - 0.5)^2 + (\eta - 0.5)^2 + (\gamma - 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.

<table><caption>TABLE 16. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>5.303e-03</td><td>1.503</td><td>4.925e-03</td><td>4.922e-01</td><td>1.791e+00</td></tr><tr><td>32</td><td>3.662e-04</td><td>8.520</td><td>3.329e-03</td><td>5.718e-01</td><td>1.746e+00</td></tr><tr><td>64</td><td>2.451e-05</td><td>55.52</td><td>2.322e-03</td><td>5.733e-01</td><td>1.725e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>3.763e-03</td><td>3.769</td><td>4.922e-03</td><td>5.166e-01</td><td>1.779e+00</td></tr><tr><td>32</td><td>1.934e-04</td><td>21.40</td><td>3.330e-03</td><td>5.720e-01</td><td>1.732e+00</td></tr><tr><td>64</td><td>4.271e-06</td><td>129.6</td><td>2.322e-03</td><td>5.736e-01</td><td>1.721e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>2.780e-03</td><td>9.072</td><td>4.917e-03</td><td>5.229e-01</td><td>1.763e+00</td></tr><tr><td>32</td><td>1.252e-04</td><td>45.97</td><td>3.330e-03</td><td>5.720e-01</td><td>1.725e+00</td></tr><tr><td>64</td><td>3.368e-06</td><td>270.2</td><td>2.322e-03</td><td>5.736e-01</td><td>1.721e+00</td></tr></tbody></table>

<figure><img src="image_21.png" /><figcaption>FIGURE 20. Location of the planes in the cube corresponding to $x = 0.1, 0.6, 0.9$.</figcaption></figure>
**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

$$ (7.7) \qquad \rho(x, y) = 3 / \left( 2 + \cos \left( 8\pi \sqrt{(x-0.5)^2 + (y-0.5)^2 + (z-0.5)^2} \right) \right). $$

<table><caption>TABLE 17. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>2.148e-02</td><td>1.716</td><td>3.108e-03</td><td>3.825e-01</td><td>1.920e+00</td></tr><tr><td>32</td><td>2.005e-03</td><td>11.26</td><td>2.198e-03</td><td>5.141e-01</td><td>1.796e+00</td></tr><tr><td>64</td><td>1.195e-04</td><td>74.83</td><td>1.533e-03</td><td>5.725e-01</td><td>1.748e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>2.030e-02</td><td>4.208</td><td>3.116e-03</td><td>3.668e-01</td><td>2.139e+00</td></tr><tr><td>32</td><td>1.580e-03</td><td>29.18</td><td>2.198e-03</td><td>5.167e-01</td><td>1.786e+00</td></tr><tr><td>64</td><td>4.936e-05</td><td>166.7</td><td>1.533e-03</td><td>5.731e-01</td><td>1.732e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>1.969e-02</td><td>9.441</td><td>3.115e-03</td><td>3.566e-01</td><td>1.919e+00</td></tr><tr><td>32</td><td>1.315e-03</td><td>58.24</td><td>2.198e-03</td><td>5.306e-01</td><td>1.760e+00</td></tr><tr><td>64</td><td>3.708e-05</td><td>349.7</td><td>1.533e-03</td><td>5.724e-01</td><td>1.729e+00</td></tr></tbody></table>

<figure><img src="image_22.png" /><figcaption>FIGURE 21. (left) A close-up of one extract part inside the sphere (right) Adaptive meshes and the monitor function.</figcaption></figure>
**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><caption>TABLE 18. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td colspan="6" style="text-align: left;">(A) Degree p = 3</td></tr><tr><td>16</td><td>3.196e-02</td><td>2.083</td><td>2.769e-03</td><td>3.252e-01</td><td>1.845e+00</td></tr><tr><td>32</td><td>5.415e-03</td><td>10.90</td><td>1.891e-03</td><td>4.868e-01</td><td>1.818e+00</td></tr><tr><td>64</td><td>3.313e-04</td><td>61.59</td><td>1.313e-03</td><td>5.751e-01</td><td>1.760e+00</td></tr><tr><td colspan="6" style="text-align: left;">(B) Degree p = 4</td></tr><tr><td>16</td><td>2.847e-02</td><td>4.734</td><td>2.753e-03</td><td>3.562e-01</td><td>1.863e+00</td></tr><tr><td>32</td><td>4.518e-03</td><td>28.79</td><td>1.890e-03</td><td>5.060e-01</td><td>1.785e+00</td></tr><tr><td>64</td><td>1.945e-04</td><td>176.1</td><td>1.313e-03</td><td>5.747e-01</td><td>1.741e+00</td></tr><tr><td colspan="6" style="text-align: left;">(C) Degree p = 5</td></tr><tr><td>16</td><td>2.599e-02</td><td>11.62</td><td>2.780e-03</td><td>3.034e-01</td><td>1.894e+00</td></tr><tr><td>32</td><td>4.037e-03</td><td>70.16</td><td>1.890e-03</td><td>4.848e-01</td><td>1.774e+00</td></tr><tr><td>64</td><td>1.518e-04</td><td>359.5</td><td>1.313e-03</td><td>5.738e-01</td><td>1.737e+00</td></tr></tbody></table>

<figure><img src="image_23.png" /><figcaption>FIGURE 22. (left) A close-up of one extract part inside the sphere (right) Adaptive meshes and the monitor function.</figcaption></figure>
**Example 7.13** (Examples of adapted mesh on more general geometries: half-torus). 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 (right) and to see what happens inside the domain, see Fig. 23 (left), and for the different degrees we validate the bijectivity and quality with fast convergence results in Table 19.

<table><caption>TABLE 19. Grid convergence analysis for various polynomial degrees $p$ along each direction.</caption><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>4.801e-02</td><td>2.637</td><td>2.029e-03</td><td>3.284e-01</td><td>1.893e+00</td></tr><tr><td>32</td><td>7.936e-03</td><td>12.84</td><td>1.338e-03</td><td>4.300e-01</td><td>1.829e+00</td></tr><tr><td>64</td><td>4.406e-04</td><td>76.00</td><td>9.252e-04</td><td>5.721e-01</td><td>1.774e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>4.680e-02</td><td>7.165</td><td>2.095e-03</td><td>2.872e-01</td><td>2.030e+00</td></tr><tr><td>32</td><td>7.310e-03</td><td>34.75</td><td>1.339e-03</td><td>4.481e-01</td><td>1.841e+00</td></tr><tr><td>64</td><td>2.774e-04</td><td>214.9</td><td>9.252e-04</td><td>5.705e-01</td><td>1.755e+00</td></tr></tbody></table><table><thead><tr><th>#cells</th><th>Err</th><th>CPU-time (s)</th><th>Qual</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>4.321e-02</td><td>15.52</td><td>2.097e-03</td><td>2.005e-01</td><td>2.032e+00</td></tr><tr><td>32</td><td>6.685e-03</td><td>76.55</td><td>1.338e-03</td><td>4.481e-01</td><td>1.799e+00</td></tr><tr><td>64</td><td>2.335e-04</td><td>454.7</td><td>9.252e-04</td><td>5.699e-01</td><td>1.743e+00</td></tr></tbody></table>

<figure><img src="image_24.png" /><figcaption>FIGURE 23. (left) A close-up of one extract part inside the sphere (right) Adaptive meshes and the monitor function.</figcaption></figure>
8. 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 two-
and 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.

*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 EFDA-fellowship and IPAM long program *Computational Methods in High Energy Density Plasmas* in 2012.
REFERENCES

1. G Awanou, *Pseudo transient continuation and time marching methods for monge-ampère type equations*, Submitted.

2. __________, *Spline element method for the monge-ampère equation*, Submitted.

3. M. J. Baines, *Least squares and approximate equidistribution in multidimensions*, Numerical Methods for Partial Differential Equations **15** (1999), no. 5, 605–615.

4. Jean-David Benamou, Brittany D. Froese, and Adam M. Oberman, *Two numerical methods for the elliptic monge-ampère equation*, ESAIM: Mathematical Modelling and Numerical Analysis **44** (2010), 737–758.

5. Jean-David Benamou, Brittany D Froese, and Adam M Oberman, *Two numerical methods for the elliptic monge-ampère equation*, ESAIM: Mathematical Modelling and Numerical Analysis **44** (2010), no. 4, 737–758.

6. Yann Brenier, *Polar factorization and monotone rearrangement of vector-valued functions*, Communications on pure and applied mathematics **44** (1991), no. 4, 375–417.

7. __________, *Polar factorization and monotone rearrangement of vector-valued functions*, Communications on Pure and Applied Mathematics **44** (1991), no. 4, 375–417.

8. C J Budd and J F Williams, *Parabolic monge-ampère methods for blow-up problems in several spatial dimensions*, Journal of Physics A: Mathematical and General **39** (2006), no. 19, 5425.

9. Chris J. Budd, Weizhang Huang, and Robert D. Russell, *Adaptivity with moving grids*, Acta Numerica **18** (2009), 111–241. MR 2506041 (2010c:65165)

10. C.J. Budd, M.J.P. Cullen, and E.J. Walsh, *Monge-ampère based moving mesh methods for numerical weather prediction, with applications to the eady problem*, Journal of Computational Physics **236** (2013), no. 0, 247 – 270.

11. Weiming Cao, Weizhang Huang, and Robert D. Russell, *Approaches for generating moving adaptive meshes: location versus velocity*, Applied Numerical Mathematics **47** (2003), no. 2, 121 – 138, Second International Workshop on Numerical Linear Algebra - Numerical Methods for Partial Differential Equations and Optimization.

12. L. Chacòn, G.L. Delzanno, and J.M. Finn, *Robust, multidimensional mesh-motion based on monge-kantorovich equidistribution*, Journal of Computational Physics **230** (2011), no. 1, 87 – 103.

13. Jean-François Cossette and Piotr K. Smolarkiewicz, *A monge-ampère enhancement for semi-lagrangian methods*, Computers & Fluids **46** (2011), no. 1, 180 – 185, 10th ICFD Conference Series on Numerical Methods for Fluid Dynamics (ICFD 2010).

14. J Austin Cottrell, Thomas JR Hughes, and Yuri Bazilevs, *Isogeometric analysis: toward integration of cad and fea*, John Wiley & Sons, 2009.

15. Edward J Dean and Roland Glowinski, *Numerical solution of the two-dimensional elliptic monge-ampère equation with dirichlet boundary conditions: an augmented lagrangian approach*, Comptes Rendus Mathematique **336** (2003), no. 9, 779 – 784.

16. Edward J. Dean and Roland Glowinski, *Numerical solution of the two-dimensional elliptic monge-ampère equation with dirichlet boundary conditions: a least-squares approach*, Comptes Rendus Mathematique **339** (2004), no. 12, 887 – 892.

17. E.J. Dean and R. Glowinski, *Numerical methods for fully nonlinear elliptic equations of the monge-ampère type*, Computer Methods in Applied Mechanics and Engineering **195** (2006), no. 13-16, 1344 – 1386, A Tribute to Thomas J.R. Hughes on the Occasion of his 60th Birthday.

18. Gian Luca Delzanno, Luis Chacón, John M Finn, Y Chung, and Giovanni Lapenta, *An optimal robust equidistribution method for two-dimensional grid adaptation based on monge-kantorovich optimization*, Journal of Computational Physics **227** (2008), no. 23, 9841–9864.

19. Gian Luca Delzanno and John M. Finn, *The fluid dynamic approach to equidistribution methods for grid adaptation*, Computer Physics Communications **182** (2011), no. 2, 330 – 346.

20. G.L. Delzanno, L. Chacòn, J.M. Finn, Y. Chung, and G. Lapenta, *An optimal robust equidistribution method for two-dimensional grid adaptation based on monge-kantorovich optimization*, Journal of Computational Physics **227** (2008), no. 23, 9841 – 9864.

21. Brittany Froese Hamfeldt and Axel GR Turnquist, *A convergence framework for optimal transport on the sphere*, Numerische Mathematik (2022), 1–31.

22. Weizhang Huang, *Practical aspects of formulation and solution of moving mesh partial differential equations*, Journal of Computational Physics **171** (2001), no. 2, 753 – 775.
23. __________, *Variational mesh adaptation: Isotropy and equidistribution*, Journal of Computational Physics **174** (2001), no. 2, 903–924.

24. __________, *Metric tensors for anisotropic mesh generation*, Journal of Computational Physics **204** (2005), no. 2, 633–665.

25. Weizhang Huang, Lennard Kamenski, and Jens Lang, *A new anisotropic mesh adaptation method based upon hierarchical a posteriori error estimates*, Journal of Computational Physics **229** (2010), no. 6, 2179–2198.

26. Weizhang Huang and Xianping Li, *An anisotropic mesh adaptation method for the finite element solution of variational problems*, Finite Elements in Analysis and Design **46** (2010), no. 1-2, 61–73, Mesh Generation - Applications and Adaptation.

27. Weizhang Huang, Jingtang Ma, and Robert D. Russell, *A study of moving mesh pde methods for numerical simulation of blowup in reaction diffusion equations*, Journal of Computational Physics **227** (2008), no. 13, 6532–6552.

28. Weizhang Huang, Yuhe Ren, and Robert D. Russell, *Moving mesh methods based on moving mesh partial differential equations*, Journal of Computational Physics **113** (1994), no. 2, 279–290.

29. Weizhang Huang and Robert D. Russell, *A moving collocation method for solving time dependent partial differential equations*, Applied Numerical Mathematics **20** (1996), no. 1-2, 101–116, Method of Lines for Time-Dependent Problems.

30. __________, *A high dimensional moving mesh strategy*, Applied Numerical Mathematics **26** (1998), no. 1-2, 63–76.

31. __________, *Adaptive mesh movement - the mmpde approach and its applications*, Partial Differential Equations (D. Sloan, S. Vandewalle, and E. Süli, eds.), Elsevier, Amsterdam, 2001, pp. 383–398.

32. Weizhang Huang and Robert D Russell, *Adaptive moving mesh methods*, vol. 174, Springer Science & Business Media, 2010.

33. Weizhang Huang and Robert D. Russell, *Adaptive moving mesh methods*, Applied mathematical sciences, Springer, New York, Heidelberg, London, 2011.

34. Weizhang Huang and Forrest Schaeffer, *An stability analysis for the finite-difference solution of one-dimensional linear convection-diffusion equations on moving meshes*, Journal of Computational and Applied Mathematics **236** (2012), no. 13, 3338–3348.

35. Weizhang Huang and Weiwei Sun, *Variational mesh adaptation ii: error estimates and monitor functions*, Journal of Computational Physics **184** (2003), no. 2, 619–648.

36. Thomas JR Hughes, John A Cottrell, and Yuri Bazilevs, *Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement*, Computer methods in applied mechanics and engineering **194** (2005), no. 39-41, 4135–4195.

37. T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs, *Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement*, Computer Methods in Applied Mechanics and Engineering **194** (2005), no. 39-41, 4135–4195.

38. Milad Khosravi and Mitra Javan, *Three-dimensional features of the lateral thermal plume discharge in the deep cross-flow using dynamic adaptive mesh refinement*, Theoretical and Computational Fluid Dynamics (2022), 1–18.

39. Jens Lang, Weiming Cao, Weizhang Huang, and Robert D. Russell, *A two-dimensional moving finite element method with local refinement based on a posteriori error estimates*, Applied Numerical Mathematics **46** (2003), no. 1, 75–94.

40. Giovanni Lapenta, *Variational grid adaptation based on the minimization of local truncation error: time-independent problems*, Journal of Computational Physics **193** (2004), no. 1, 159–179.

41. Xianping Li and Weizhang Huang, *An anisotropic mesh adaptation method for the finite element solution of heterogeneous anisotropic diffusion problems*, Journal of Computational Physics **229** (2010), no. 21, 8072–8094.

42. P-L Lions, Neil Sydney Trudinger, and John IE Urbas, *The neumann problem for equations of monge-ampère type*, Communications on pure and applied mathematics **39** (1986), no. 4, 539–563.

43. Grégoire Loeper and Francesca Rapetti, *Numerical solution of the monge-ampère equation by a newton's algorithm*, Comptes Rendus Mathematique **340** (2005), no. 4, 319–324.
44. Changna Lu, Weizhang Huang, and Erik S. Van Vleck, *The cutoff method for the numerical computation of nonnegative solutions of parabolic pdes with application to anisotropic diffusion and lubrication-type equations*, Journal of Computational Physics (2013), no. 0, —.

45. Robert E Lynch, John R Rice, and Donald H Thomas, *Direct solution of partial difference equations by tensor product methods*, Numerische Mathematik **6** (1964), no. 1, 185–199.

46. Andrew TT McRae, Colin J Cotter, and Chris J Budd, *Optimal-transport-based mesh adaptivity on the plane and sphere using finite elements*, SIAM Journal on Scientific Computing **40** (2018), no. 2, A1121–A1148.

47. Ahmed Ratnani, Said Hadjout, Yaman Guclu, and Emily Bourne, *Pyccel, python static compiler for high performance computing*, 2019.

48. Giancarlo Sangalli and Mattia Tani, *Isogeometric preconditioners based on fast solvers for the sylvester equation*, SIAM Journal on Scientific Computing **38** (2016), no. 6, A3644–A3671.

49. Deepti Shakti, Jugal Mohapatra, Pratibhamoy Das, and Jesus Vigo-Aguiar, *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*, Journal of Computational and Applied Mathematics **404** (2022), 113167.

50. Mohamed Sulman, J.F. Williams, and R.D. Russell, *Optimal mass transport for higher dimensional adaptive grid generation*, Journal of Computational Physics **230** (2011), no. 9, 3302–3330.

51. Mohamed M. Sulman, J.F. Williams, and Robert D. Russell, *An efficient approach for the numerical solution of the monge-ampère equation*, Applied Numerical Mathematics **61** (2011), no. 3, 298–307.

52. Anna Tagliabue, Luca Dede, and Alfio Quarteroni, *Isogeometric analysis and error estimates for high order partial differential equations in fluid dynamics*, Computers & Fluids **102** (2014), 277–303.

53. Hilary Weller, Philip Browne, Chris Budd, and Mike Cullen, *Mesh adaptation on the sphere using optimal transport and the numerical solution of a monge-ampère type equation*, Journal of Computational Physics **308** (2016), 102–123.

54. Xuefei Yuan, Stephen C. Jardin, and David E. Keyes, *Moving grids for magnetic reconnection via newton-krylov methods*, Computer Physics Communications **182** (2011), no. 1, 173–176, Computer Physics Communications Special Edition for Conference on Computational Physics Kaohsiung, Taiwan, Dec 15-19, 2009.

55. —, *Numerical simulation of four-field extended magnetohydrodynamics in dynamically adaptive curvilinear coordinates via newton-krylov-schwarz*, Journal of Computational Physics **231** (2012), no. 17, 5822–5853.

56. V. Zheligovsky, O. Podvigina, and U. Frisch, *The monge-ampère equation: Various forms and numerical solution*, Journal of Computational Physics **229** (2010), no. 13, 5043–5061.
MODELING, SIMULATION AND DATA ANALYSIS, MOHAMMED VI POLYTECHNIC UNIVERSITY, BENGUERIR 43150, LOT 660, HAY MOULAY RACHID MOROCCO

*Current address:* Modeling, Simulation and Data Analysis, Mohammed VI Polytechnic University, Benguerir 43150, Lot 660, Hay Moulay Rachid Morocco

*Current address:* Laboratoire J.A.Dieudonné, Université de Nice Sophia-Antipolis, France

*Email address:* mustapha.bahari@um6p.ma

LABORATOIRE J.A.DIEUDONNÉ, UNIVERSITÉ DE NICE SOPHIA-ANTIPOLIS, FRANCE

*Current address:* Laboratoire J.A.Dieudonné, Université de Nice Sophia-Antipolis, France

*Email address:* abderrahmane.habbal@unice.fr

MODELING, SIMULATION AND DATA ANALYSIS, MOHAMMED VI POLYTECHNIC UNIVERSITY, BENGUERIR 43150, LOT 660, HAY MOULAY RACHID MOROCCO

*Current address:* Modeling, Simulation and Data Analysis, Mohammed VI Polytechnic University, Benguerir 43150, Lot 660, Hay Moulay Rachid Morocco

*Email address:* ahmed.ratnani@um6p.ma

MAX PLANCK INSTITUTE FOR PLASMA PHYSICS, BOLTZMANNSTR. 2, 85748 GARCHING, GERMANY

*Current address:* Max Planck Institute for Plasma Physics, Boltzmannstr. 2, 85748 Garching, Germany

*Email address:* Eric.Sonnendruecker@ipp.mpg.de