A Least-Squares/Relaxation Method for the Numerical Solution of the Three-Dimensional Elliptic Monge–Ampère Equation

In this article, we address the numerical solution of the Dirichlet problem for the three-dimensional elliptic Monge–Ampère equation using a least-squares/relaxation approach. The relaxation algorithm allows the decoupling of the differential operators from the nonlinearities. Dedicated numerical solvers are derived for the efficient solution of the local optimization problems with cubicly nonlinear equality constraints. The approximation relies on mixed low order finite element methods with regularization techniques. The results of numerical experiments show the convergence of our relaxation method to a convex classical solution if such a solution exists; otherwise they show convergence to a generalized solution in a least-squares sense. These results show also the robustness of our methodology and its ability at handling curved boundaries and non-convex domains.


Introduction
The Monge-Ampère equation can be considered as the prototypical example of fully nonlinear elliptic equations [11,21,29]. Theoretical investigations of fully nonlinear equations started many years ago [6,39] and have received a lot of attention lately [4,5,8,14,16,19,[35][36][37][38], due to the many applications involving this type of equations, for instance in finance [41], in seismic wave propagation [13], in geostrophic flows [15], in differential geometry [2,18], and in mechanics and physics. However, the Monge-Ampère equation in three space dimensions is a complicated problem, which is still lacking a full theoretical understanding, particularly when the domain of interest is not strictly convex.
From a computational point of view, various approaches have been identified, relying in particular on finite difference [5,13] or finite element [3,9,17,27] approximations.
Since the Monge-Ampère equation may not have smooth classical solutions, even for smooth data (see [11]), notions of generalized solutions have been introduced, such as Aleksandrov solutions [1], and viscosity solutions [31,38]. Actually, another notion of generalized solution for fully nonlinear elliptic equations has been introduced relatively recently, namely generalized solutions in a least-squares sense. These least-squares generalized solutions have been obtained via augmented Lagrangian or least-squares/relaxation approaches [9,12,24]. The least-squares approach will be the one used in this article.
We consider in this article three-dimensional bounded domains Ω with a Lipschitz continuous boundary ∂Ω. For data f and g sufficiently smooth, it makes sense from the existence, uniqueness and regularity results reported in, e.g., [10], to look for ψ belonging to H 2 (Ω) and convex, solution to the following Dirichlet problem for the Monge-Ampère equation Indeed, note that, if Ω is a bounded strictly convex domain of R 3 with a C ∞ boundary ∂Ω, this problem has a unique convex solution belonging to C ∞ (Ω) (see, e.g., [10]). On the other hand, if Ω is a bounded strictly convex domain of R 3 , the Dirichlet problem for the Monge-Ampère equation has a unique convex generalized solution belonging to C 0 (Ω) ∩ W 2,∞ loc (Ω) (see, e.g., [29]). Therefore, H 2 (Ω) is a reasonable choice to look for a solution to the Monge-Ampère equation in three space dimensions.
The chosen least-squares formulation we used here consists in minimizing the L 2 -distance between D 2 ψ and a matrix-valued function p ∈ (L 2 (Ω)) 3×3 , where ψ satisfies the boundary conditions of the problem (but not the Monge-Ampère equation) and p satisfies det p = f . Using a relaxation algorithm to minimize such a distance, we obtained a solution method where one solves alternatively, until convergence, a sequence of linear variational problems (to be approximated by mixed finite element methods) and a sequence of cubicly constrained algebraic optimization problems. Using a similar approach, we have been able to compute generalized solutions of the Monge-Ampère equation when this problem has no classical solutions in two dimensions of space [9]. In this article, the methods discussed in [9] have been generalized to three-dimensional problems.
In this article, the linear variational problems we mentioned above are solved by a preconditioned conjugate gradient algorithm, whose computer implementation relies on low order (P 1 or Q 1 ) mixed finite element approximations. The local cubicly constrained optimization problems are solved by Newton-like methods [42] or time-stepping methods associated with a dynamical flow (see, e.g., [30]). The main difference between the two-dimensional case discussed in [9] , and the present article is the cubic nature (vs quadratic) of the equality constraints.
The structure of this article is as follows. In Sect. 2, we describe the proposed methodology, while the relaxation algorithm is described in Sect. 3. Sections 4 and 5 detail the algebraic and differential solvers respectively. The mixed finite element discretization is discussed in Sect. 6. The method is applied in Sect. 7 to the solution of several numerical examples, including examples without a classical exact solution. Finally, in Sect. 8, some numerical results are presented where P 1 finite elements have been replaced with Q 1 ones.
The numerical solution of the 3D elliptic Monge-Ampère equation has been discussed in [8] using piecewise P 3 continuous finite element approximations. A fast multigrid scheme has been presented in [34], while smooth cases on structured meshes have been considered in [3] (with numerical results that are consistent with [9]). However, to the best of our knowledge, the method discussed in the present article is one of the very few able to solve the 3D elliptic Monge-Ampère equation on domains with curved boundaries, using piecewise P 1 continuous finite element approximations associated with unstructured meshes, while preserving optimal, or nearly optimal, orders of convergence for the approximation errors, including situations where the solution does not have the C 2 (Ω) regularity.

Mathematical Formulation and Least-Squares Approach
Let Ω be a bounded convex domain of R 3 ; we denote by Γ the boundary of Ω. The Dirichlet problem for the elliptic Monge-Ampère equation reads as follows: where is the Hessian of the unknown function ψ.
Among the various methods available for the solution of (1) discussed in the introduction, we advocate a nonlinear least-squares method that relies on the introduction of an additional auxiliary variable. Namely, we look for: where: using the Fröbenius norm and inner product defined by |T| = √ T : T, S : T = 3 i, j=1 s i j t i j , for all S = (s i j ), T = (t i j ) ∈ R 3×3 . The functional spaces in (2) are respectively defined by: We assume that f ∈ L 2/3 (Ω) and g ∈ H 3/2 (Γ ), so that V g and Q f are both non-empty. Note that the space Q is a Hilbert space for the inner product (q, q ) → Ω q : q dx, and the associated norm.

Relaxation Algorithm
In order to decouple differential operators from the nonlinearity, we suggest using a relaxation algorithm of the Gauss-Seidel-type. Furthermore, we wish to compute a convex solution to problem (2) (or at least to force the convexity of the solution), and thus advocate the following approach. First, the initialization is performed by solving: Then, for n ≥ 0 and assuming that ψ n is known, one computes p n , ψ n+1/2 and ψ n+1 as follows: with 1 ≤ ω ≤ ω max < 2. For the numerical experiments presented in Sect. 7, we used ω ≡ 1 (unless otherwise specified). The rationale behind the initialization procedure (4) is based on the above isotropic assumption. If we denote the eigenvalues of D 2 ψ by λ i , i = 1, 2, 3, the Monge-Ampère equation reads λ 1 λ 2 λ 3 = f . If λ 1 , λ 2 and λ 3 are "close" from each other (and thus all equal to, let's say, λ), we have λ 3 = f , and thus λ = 3 √ f . Therefore: Since the initialization procedure is based on an isotropic assumption, the case when the eigenvalues are very different from each other has to be put under scrutiny in the numerical experiments. In the sequel, the numerical algorithms used for the solution of problems (5) and (6) will be discussed in details. Despite several investigations, note that the question of the uniqueness of the solution to the local problem (5) still remains an open question.
In the case of the 2D Monge-Ampère equation, this step led to a class of quadratically constrained minimization problems, whose solution was addressed in [9,28] using a dedicated algorithm. This dedicated approach does not apply when the constraint is a cubic equation as here, implying that other approaches have to be considered. The two approaches below rely on an appropriate re-parameterization of the problem in order to transform the constrained minimization problem into an unconstrained one; both ending up with some kind of Newtonrelated methods.

A Reduced Newton Method
For a. e. x ∈ Ω, (9) is an algebraic optimization problem. Using a Cholesky decomposition of q, we write q = LDL t , where This re-parameterization is arbitrary but serves two purposes: first, it guarantees that all eigenvalues are strictly positive (convexity of the local solution). Second, the constraint detq = f (x) is automatically satisfied. It thus allows one to replace (9) by an unconstrained minimization problem in the variable X := (a, b, c, ρ 1 , ρ 2 ). For the sake of simplicity, we do not write the dependency on x ∈ Ω anymore. The problem becomes: The first order optimality conditions corresponding to (11) can formally be written as ∇ X G(X) = 0.
This nonlinear system can be solved with a safeguarded Newton method for the variable X. Namely, given X 0 ∈ R 5 , solve, for k ≥ 0: followed by where λ k ∈ R + is a step-length to be adapted according to some Armijo rule (see, e.g., [7]). Typically, we update the step-length if ∇G(X k+1 ) > (1 − αλ k ) ∇G(X k ) , where α = 10 −4 and ||·|| denotes the canonical Euclidian norm of R 5 , and set in that case λ k+1 = 1 2 λ k . The stopping criterion is based on the residual value ∇G(X k ) , and the iterations are stopped if ∇G(X k ) < ε Newton , where ε Newton is a given tolerance.

A Runge-Kutta Method for the Dynamical Flow Problem
An alternative re-parameterization of the nonlinear problem (9) can be considered, based on a eigenvalues-eigenvectors decomposition, in the spirit of the approach in [28]. Namely, consider q = Q Q t , where Q ∈ O(3) ⊂ R 3×3 is the orthogonal matrix whose columns represent the eigenvectors of q (O(3) being the group of the 3 × 3 orthogonal matrices), and ∈ R 3×3 is the diagonal matrix whose diagonal elements are the corresponding eigenvalues. We can denote Again, this parameterization is not unique; it ensures that the relation detq = f (x) is automatically satisfied and that the eigenvalues are positive in order to enforce convexity of the local solutions. The property Q ∈ O(3) implies the following constraints for its column vectors T i , i = 1, 2, 3 (where T i · T j denotes the dot product of vectors T i and T j ): Let us define the variables Y = (ρ 1 , ρ 2 , T 1 , T 2 , T 3 ) ∈ R 11 . Problem (9) can be rewritten as Below, for i = 1, 2, 3, we will denote by |T i | the quantity (T i · T i ) 1/2 . We penalize the equality constraints in order to obtain an unconstrained problem that can be solved by a Newton approach. Let ε 1 , ε 2 > 0 be two given (small) parameters. The constraints are taken into account by penalization, leading to the following unconstrained minimization problem: where Similarly to the solution of (11), the first order optimality conditions associated with (12) can be written as In order to smoothen the transition to critical point(s), we favor an evolutive formulation of the first order optimality conditions (in the sense of a flow problem in the dynamical systems terminology), which read as follows: find Y : (0, +∞) → R 11 such that The steady state solution of (13) (14) corresponds to the desired critical point. In order to increase the stability of the numerical scheme and allow larger time steps and therefore a faster convergence to the steady state solution, it is customary to modify (13) into a modified flow problem [32], namely: find Y : (0, +∞) → R 11 such that with the same initial condition. The stability of the scheme is important here since we are aiming at solving such a flow problem for a.e. x ∈ Ω, which requires an efficient numerical algorithm. The additional computational cost induced by the introduction of the term System (15) is solved by a two-stage (second order explicit) Runge-Kutta method (see, e.g., [30]) in order to capture steady state solutions. Let Δt be a given time step, t n = nΔt and Y n Y(t n ), n = 0, 1, . . .. Let us define Y 0 = Y 0 ; then, at each time step, solve An adaptive time stepping strategy for Runge-Kutta methods is incorporated to the numerical algorithm; numerical experiments will show that the adaptive time step is particularly useful at the beginning of the outer iterations loop, when the initial solution is not close to the final steady state solution.
It is worth noticing that, if we treat the modified flow problem with a first order Euler explicit scheme, it leads to solving at each time step this problem corresponds actually to a classical safeguarded Newton method, reminiscent to the one we presented in Sect. 4.2, with Δt playing the role of the step-length λ. With this remark, the adaptive time stepping algorithm for Runge-Kutta schemes can be seen as an adaptive Armijo-like rule, with Δt = λ. Furthermore, one can see that the Runge-Kutta approach is slower than the reduced Newton strategy (since it corresponds to solving two Newton-type systems at each time step), but it is more accurate since the two-step Runge-Kutta scheme is a higher order method than the Euler scheme. Finally, a study of the stability of Runge-Kutta schemes [30] shows that its stability properties are better than those of the Euler scheme.

Numerical Solution of the Linear Variational Problems
Written in variational form the Euler-Lagrange equation of the sub-problem (6) reads as follows: find ψ n+1/2 ∈ V g satisfying: The linear variational problem (16) is well-posed and belongs to the following family of linear variational problems: with the functional L(·) linear and continuous over H 2 (Ω); problem (17) is a bi-harmonic type problem, which can be solved by a conjugate gradient algorithm operating in wellchosen Hilbert spaces (see, e.g., [22,Chapter 3]). Here, our conjugate gradient algorithm operates in the spaces V 0 and V g , both spaces being equipped with the inner product defined by (v, w) → Ω ΔvΔwdx, and the corresponding norm. It reads as follows: Step 1 u 0 ∈ V g given. (18) Step 2 Solve: find g 0 ∈ V g satisfying and set Then, for k ≥ 0, u k , g k and w k being known, the last two different from zero, we compute u k+1 , g k+1 and, if necessary, w k+1 as follows.
and compute Step 4 Compute If δ k < ε, take u = u k+1 ; otherwise, compute: and Step 5 Do k + 1 → k and return to Step 3.

Mixed Finite Element Approximation
Considering the highly variational flavor of the methodology discussed in the preceding sections, it makes sense to look for finite element methods for the solution of (1). We will use a mixed finite element approximation (closely related to those discussed in, e.g., [25] for the solution of linear and nonlinear bi-harmonic problems) with low order (piecewise linear and globally continuous) finite elements on a partition of Ω made of tetrahedra. The modification of the numerical approximation method obtained when replacing the P 1 based finite element spaces on tetrahedra by the Q 1 based finite element spaces associated with partitions of Ω made of hexahedra is discussed in Sect. 8 with some numerical experiments.

Finite Element Spaces
For simplicity, let us assume that Ω is a bounded polyhedral domain of R 3 , and define T h as a finite element partition of Ω made out of tetrahedra (see, e.g., [23, with P 1 the space of the three-variable polynomials of degree ≤ 1. We define also V 0h as In the sequel, V 0h will be used to approximate both H 1 0 (Ω) and H 2 (Ω) ∩ H 1 0 (Ω).

Finite Element Approximation of the Monge-Ampère Equation
When solving (17) by the conjugate gradient algorithm (18)-(27), one has to (i) compute the discrete analogues of the second order derivatives, e.g., D 2 w k and D 2 u 0 , and (ii) solve biharmonic problems such as (19) and (21). Concerning (i) we will approximate the second order derivatives by functions belonging to V 0h , but this has to be handled carefully. For a function ϕ belonging to H 2 (Ω), it follows from Green's formula that, for i, j = 1, 2, 3: Consider now ϕ ∈ V h . We define the discrete analogue D 2 hi j ϕ ∈ V 0h of the second derivative The functions D 2 hi j ϕ are thus uniquely defined; in order to simplify the computation of the above discrete second order partial derivatives, we could use the trapezoidal rule to evaluate the integrals in the left hand sides (mass lumping). Since the Hessian matrix D 2 ψ is in (L 2 (Ω)) 3×3 , and since H 1 0 (Ω) is dense in L 2 (Ω), the approximation D 2 h ψ is considered to be in V 0h . Indeed, the underlying method being a collocation method, the Monge-Ampère equation is approximated and imposed at the internal vertices of Ω only.
As emphasized in [9,40], when using piecewise linear mixed finite elements, the approximation of the error on the second derivatives of the solution ψ is, in general, O(1) in the L 2 -norm. Therefore, the convergence properties of the global algorithm strongly depends on the type of partition of Ω one employs, and could be completely jeopardized in some situations. One way to improve the approximation properties of the discrete second order derivatives D 2 hi j ϕ is to use, as in [9], a Tychonoff-like regularization [43]. Let us introduce a stabilization constant C (to be calibrated in the numerical experiments), and replace the previous variational problem by: for all i, j, The rationale behind (30) is to correct the approximation errors of (29), which deteriorate when h → 0. This behavior being associated with the high modes of the approximate solutions. Similar ideas have been used for approximations of the Stokes problem in incompressible fluid mechanics when using essentially the same finite element spaces for the velocity components and the pressure (see, e.g., [23,Chapter 5]). Among the possible cures of this unwanted behavior, let us mention the use of higher degree polynomials to define the finite element spaces. However since piecewise affine approximations are optimal for handling domains of (almost) arbitrary shape and with curved boundaries, we will try to rescue them following the approach in [9]. The first step is to approximate (28) with ε > 0 a "small" positive number (of the order of h 2 in practice). Problem (31) is well-posed and one can easily show that In this article, we approximate both (28) and (31) by (30). Numerical experiments show that this regularization procedure provides approximations of optimal or nearly optimal orders.
, and assuming that the boundary function g is continuous over Γ , the affine space V g can be approximated by Concerning issue (ii), that is the solution of the bi-harmonic problems encountered in the conjugate gradient algorithm (18)- (27), after space discretization the resulting discrete bi-harmonic problems are all particular cases of with It follows from this definition that the discrete bi-harmonic problem (33) is equivalent to the following system of easy to solve discrete Poisson-Dirichlet problems:

Discrete Formulation of the Least-Squares Method
We define the discrete analogues of spaces Q and Q f as follows: We associate with V h (or V 0h and V gh ) and Q h , the discrete inner products: where A k is the volume of the polyhedral domain which is the union of those tetrahedra of T h which have P k as a common vertex.
The solution of (32) is then addressed via a nonlinear least-squares method, namely find where:

A Discrete Relaxation Algorithm
The discrete relaxation algorithm we employ reads as follows: First find For n ≥ 0, assuming that ψ n h is known, compute p n h , ψ n+1/2 h and ψ n+1 h as follows:

Finite Element Approximation of the Local Nonlinear Problems
The finite dimensional minimization problems, discrete analogues of (9), are approximated, at each grid point P k ∈ Σ 0h , by: The methods discussed in Sect. 4 still apply.

Finite Element Approximation of the Linear Variational Problems
The variational problems arising in the discrete version of the relaxation algorithm can be solved similarly as in the continuous case using a conjugate gradient algorithm. Let us point out however a particularity that arises in the discrete case. The discrete version of (16) reads as follows: find ψ n+1/2 h ∈ V gh satisfying: The linear problem (36) leads to excessive computer resource requirements, which could be acceptable for two-dimensional problems, but become prohibitive for three dimensional ones.
(Indeed, to derive the linear system equivalent to (36), we need to compute-via the solution of (30)-the matrix-valued functions D 2 h w j , where the functions w j form a basis of V 0h .) To avoid this difficulty, we are going to employ, as previously discussed in [9], an adjoint equation approach (see, e,g., [26]) to derive an equivalent formulation of (36), well-suited to a solution by a conjugate gradient algorithm. This adjoint approach reads as where ∂ J h ∂ϕ (ϕ, q), θ denotes the action of the partial derivative ∂ J h ∂ϕ (ϕ, q) on the test function θ . In order to solve (37), we first determine D 2 hi j ϕ via (30). Then, we find λ hi j ∈ V 0h , 1 ≤ i, j ≤ 3, by solving the (adjoint) systems: and we can show (see, e.g., [26]) that, for all (ϕ h , p h ) ∈ V gh × Q h : This last relation can be used directly in the conjugate gradient algorithm (18)- (27), to solve (19) and (21). For instance, the discrete equivalent of (21) consists in findingḡ k h ∈ V 0h such that (with obvious notation):

Numerical Results
The first numerical results we are going to report, in order to validate our methodology, are associated with the unit cube Ω = (0, 1) 3 . Two types of partitions of the unit cube are considered, in order to study the mesh-dependence of our methods. These partitions have been constructed by using either advancing front 3D procedures, or successive extrusions [20], and are visualized in Fig. 1. All experiments were performed on a desktop computer with Intel ® Xeon(R) CPU E5-1650 v3 @ 3.50 GHz × 12. In the numerical examples presented hereafter, we consider C = 0 for the structured mesh (i.e. no Tychonoff regularization) and C = 1 for the isotropic mesh. The local nonlinear problems are solved with a stopping criterion of ε Newton = 10 −9 on the residual for the Newton method, with a maximal number of iterations equal to 1000. When using the Runge-Kutta method for the dynamical flow approach, the time step is set to Δt = 0.1, and is reduced only if needed (only for the first 2-3 times steps usually); the maximal number of iterations is 20,000 and the stopping criterion is ε = 10 −7 on two successive iterates.
Unless otherwise specified, the relaxation parameter is set to ω = 1 at the beginning of the outer iterations, and gradually increased to 2 to speed up convergence. The conjugate gradient algorithm for the solution of the variational problems has a stopping criterion of ε = 10 −8 on successive iterates, with a maximal number of iterations equal to 100. Actually, numerical experiments show that the number of conjugate gradient iterations is never larger than 35. The outer relaxation algorithm has a stopping criterion of ε = 5×10 −4 on the residual D 2 h ψ n h − p n h 0h , or on successive iterates if the problem does not admit a classical solution (see Sect. 7.3), with a maximal number of iterations equal to 5000.

Polynomial Examples
Let us consider a first example involving a smooth, "reasonably" isotropic, exact solution. More precisely, let us consider as exact solution By a direct calculation, one obtains λ 1 = 1, λ 2 = 5, λ 3 = 15, and, therefore the data for the Monge-Ampère problem correspond to f (x, y) = 75 and g(x, y, z) = 1 2 (x 2 + 5y 2 + 15z 2 ).    the same results. For the structured mesh, the method is globally second-order, resp. firstorder, convergent for the L 2 (resp. H 1 ) norm of the approximation error. Table 1 confirms those convergence results, for both structured and isotropic meshes and for the approach based on the Runge-Kutta approximation for the dynamical flow problem in R 11 . We observe that, for the structured meshes, we have text-book second and first order convergence, while the orders of convergence deteriorate for the isotropic unstructured meshes. Table 2 provides CPU times versus the number of degrees of freedom involved in the numerical approximation. Comparing mainly with [8,34] (who also performed CPU times investigations), this test case is more stringent since the eigenvalues of the Hessian are not close from each other, which makes the problem less isotropic than examples used in the literature (typically the example presented in Sect. 7.2). The reduced Newton method is used for the solution of the local nonlinear problems  The second test problem that we consider still has a polynomial exact solution, but this solution is much more anisotropic than the one in (38), since it is given by

Isotropic Structured
Here λ 1 = 1, λ 2 = 10, λ 3 = 100, and the data for the Monge-Ampère problem are given by f (x, y) = 1000 and g(x, y, z) = 1 2 (x 2 + 10y 2 + 100z 2 ). This time, the initialisation (4) of the relaxation algorithm is not close to the solution, implying, as expected, that more iterations are needed to achieve convergence. Figure 3 illustrates the convergence orders for the computed approximation error for both types of meshes. Despite the anisotropy of the solution of this second test problem, the approximation errors are similar to those associated with the first test problem, that is perfect second and first orders with the structured meshes and slightly lower convergence orders for the anisotropic unstructured meshes. Note that here we have to choose ω = 0.5 initially (under-relaxation), and increase it gradually to 2, to ensure convergence of the relaxation algorithm, and C = 2.5 for the isotropic mesh.
The least-squares approach never enforces directly the solution ψ to be convex. However, it enforces explicitly the additional matrix-valued variable p to be symmetric positive definite. It is thus remarkable that the positive definiteness property of p translates automatically to the Hessian of the main variable ψ. More precisely, when the Monge-Ampère problem has a smooth classical solution, the converged iterate satisfy D 2 ψ = p, and thus the convexity of ψ is automatically verified. For this test problem, and for both types of meshes, D 2 ψ is indeed symmetric positive definite for all grid points.

A Smooth Exponential Example
The third test problem we consider has a smooth exponential exact solution, namely the radial function ψ defined by ψ(x, y, z) = e 1 2 (x 2 +y 2 +z 2 ) , (x, y, z) ∈ Ω. (40) This test problem generalizes to three dimensions a two-dimensional one commonly used in the community for Monge-Ampère solver benchmarking (see, e.g., [9,16]). Let us denote x 2 + y 2 + z 2 by r . Taking advantage of the fact that, if φ is a radial function, one has (with obvious notation) det D 2 φ = φ (φ /r ) 2 , the data for the Monge-Ampère-Dirichlet problem (1) associated with the above function ψ are f (x, y, z) = (1 + r 2 )e 3r 2 /2 , and g(x, y, z) = e r 2 /2 . The stopping criterion for the relaxation algorithm is D 2 h ψ n h − p n h 0h < 5 × 10 −4 , and C = 1 for the isotropic unstructured mesh (as for the first test problem). Figure 4 visualizes the L 2 (Ω) and H 1 (Ω) computed approximation errors for both approaches for the solution of the local nonlinear problems. The conclusions are similar: both Newton and Runge-Kutta methods provide exactly the same results, and the method is globally second-order convergent for the L 2 norm. Table 3 confirms these convergence results, showing in particular no loss of convergence orders for the unstructured isotropic mesh. Table 4 provides CPU times versus the number of degrees of freedom involved in the numerical approximation for both structured and unstructured discretizations of the unit cube. When using a structured mesh of the unit cube, the performance of the algorithm is comparable to the other algorithms from the literature, albeit slightly less efficient. Using an unstructured isotropic mesh degrades the performance of the algorithm; computational performance for the algebraic part is identical, the difference solely coming from the increased number of conjugate gradient iterations. Note that, for this test case, the Hessian matrix D 2 ψ The local optimization problems are solved using the Runge-Kutta based method described in Sect. 4.3 (Ω = (0, 1) 3 ) The reduced Newton method is used for the solution of the local nonlinear problems admits the eigenvalues λ 1 = e r 2 /2 , λ 2 = e r 2 /2 and λ 3 = (1 + r 2 )e r 2 /2 . This example is thus rather isotropic (the eigenvalues of the Hessian are close to each other).

Non-smooth Test Problems
Some of the test problems we are going to consider in this section do not have exact solution with the H 2 (Ω)-regularity or may have no solution at all (but may have generalized solutions). These non-smooth problems are therefore ideally suited to test the robustness of our methodology, and its ability at capturing generalized solutions when no exact solution does exist. The local optimization problems are solved using the Runge-Kutta based method described in Sect With Ω still being the unit cube (0, 1) 3 , the first problem that we consider is the particular case of problem (1) which has the convex function ψ defined, for R ≥ √ 3, by When R > √ 3, this function ψ belongs to C ∞ (Ω), while ψ ∈ C 0 (Ω) ∩ W 1,s (Ω), with 1 ≤ s < 2, if R = √ 3. It is therefore interesting to see how our methodology can handle the possible non-smoothness of the particular problem (1) associated with f and g defined by f (x, y, z) = R 2 (R 2 − r 2 ) 5/2 and g(x, y, z) = − R 2 − (x 2 + y 2 + z 2 ), with, as earlier, r = x 2 + y 2 + z 2 . (a simple way to compute f is to use, as before, the relation det D 2 φ = φ (φ /r ) 2 if φ is a radial function). Of particular interest will be the behavior of our methodology when R → √ 3 from above (or even when R = √ 3). On Tables 5 and 6, we have reported, for R = √ 6 and R = √ 3, computed approximation errors and orders of convergence as h varies, together with the number of iterations necessary to achieve convergence of the relaxation algorithm. For R = √ 6, the convergence orders of the approximation errors are the ones we expect, namely second order (resp., first order) for the L 2 -norm (resp., H 1 -norm), the number of iterations being pretty low. The case R = √ 3 is more interesting; indeed despite the solution singularity at point (1, 1, 1), the L 2 (Ω)-norm of the computed approximation error ψ h − ψ still decreases super-linearly with respect to h (for both mesh families), while the related H 1 (Ω)-norm stays stable around 0.62 for the same values of h.
To conclude this section, we will consider the particular problem (1) associated with Ω = (0, 1) 3 , f = 1 and g = 0. For these particular data, problem (1) has no smooth solution (the arguments developed in [9,23] for the related two-dimensional problem still apply here). Figure 5 shows different features of the approximated solution inside the unit cube. The stopping criterion for this particular case without a classical solution is ψ n+1 h − ψ n h 0,h < 10 −5 . When studying the number of outer iterations of the relaxation algorithm, we observe that the number of iterations is larger for structured meshes than isotropic ones, and that it The local optimization problems are solved using the Runge-Kutta based method described in Sect    increases as expected when h → 0. Figure 5 (bottom row) visualizes graphs of the computed solutions restricted to the lines y = z = 1/2 and x = y, z = 1/2 for x ∈ (0, 1), and shows little influence of the type of partition on the solution. We can also observe that D 2 ψ is symmetric positive definite for 100% of the grid points, independently of the nature of the discretization when h 0.04, even though the Monge-Ampère equations does not have a classical solution, that is D 2 ψ = p. The (necessary) loss of convexity of the solution is thus located (near the corners) in a region smaller than the mesh size. When arbitrarily refining the mesh in a corner of the domain, we observe that the Hessian D 2 ψ is not symmetric positive definite when evaluated in some grid points in a neighborhood of size 10 −3 around that corner. This effect is highlighted when calculating D 2 ψ h − p h L 2 , using a structured mesh of the unit cube, both on Ω, but also on Ω ⊂ Ω, as illustrated in Table 7 for Ω = (0.2, 0.8) 3 . These results show that the error inside the domain Ω = (0.2, 0.8) 3 is significantly smaller than the error on Ω, implying that the error is mainly committed near the boundary.

Curved Boundaries and Non Convex Domains
In order to further validate the robustness and flexibility of our methodology, we are going to consider test problems where Ω has a curved boundary and/or is nonconvex. The first domain with a curved boundary we consider is the unit ball B 1 = (x, y, z) ∈ R 3 , x 2 + y 2 + z 2 < 1 . Assuming that Ω = B 1 , f = 1 3 √ 3 and g = 0, the unique convex solution ψ of the related Monge-Ampère-Dirichlet problem (1) is given by On Fig. 6 (left) we have visualized a typical finite element mesh used for computation and some cuts of the computed solution. On Fig. 6 (right) and Table 8 we have provided information on the L 2 (Ω) and H 1 (Ω)-norms of the approximation error ψ h − ψ and of the related rates of convergence, and on the number of relaxation iterations necessary to achieve convergence. Albeit the L 2 (Ω)-approximation error is O(h 1.8 ), approximately, these numerical results show that our methodology can handle rather accurately domains Ω with curved boundaries. Table 9 provides CPU times versus the number of degrees of freedom involved in the numerical approximation of the solution on the unit sphere. Results are comparable to those obtained when using unstructured meshes on the unit cube, and thus show that the curved boundaries are handled appropriately. The non-convexity of Ω may prevent problem (1) to have solutions (see, e.g., [11]). However, it makes sense to assess the capabilities of our methodology at handling problems  The local optimization problems are solved using the reduced Newton method described in Sect. 4.2 (Ω = B 1 ) The reduced Newton method is used for the solution of the local nonlinear problems having smooth solutions despite the non-convexity of Ω. To do so, we consider the particular problem (1) where: (i) Ω is the subset of B 1 obtained by removing from this ball a part of angular size θ , symmetric about Ox and oriented along the Oz axis (as shown on Fig. 7 for θ = π/2 and θ = π/9), (ii) f = 1/(3 √ 3), g being the restriction to ∂Ω of the function ψ defined by (41). The function ψ defined by (41) is clearly a convex solution of the above problem (1). The solution methodology discussed in Sects. 3-6 still applies for this case where an exact smooth solution does exist, some of the numerical results we obtained being reported  (41), Ω being the truncated unit ball (left: θ = π/2, right: θ = π/9). Second row: visualization of the truncated balls (left: θ = π/2, right: θ = π/9). Third row: visualization of the restrictions of the computed solutions to the plane z = 0 (left: θ = π/2, right: θ = π/9) in Fig. 7. We observe in particular that the convergence orders are essentially independent of the value of the re-entrant angle θ .

An Alternative Discretization Method Based on Q 1 Finite Elements
We finally report some numerical results that were obtained using Q 1 finite elements for the space discretization instead of the P 1 finite elements used earlier. The finite element library libmesh [33] has been used for implementation. The discretization of the unit cube Ω = (0, 1) 3 is based on a structured mesh of elementary cubes, as visualized in Fig. 8. The least-squares/relaxation methodology is still applicable. The nonlinear problems (9) are solved for each vertex of the hexahedral mesh with the Newton and Runge-Kutta methods discussed in Sects. 4.2 and 4.3, respectively. The variational problem (17) is solved by a conjugate gradient algorithm, using Gauss quadrature rules (of order up to 4) for the numerical computation of integrals; all other techniques and approaches remain the same.  The space approximation relies on Q 1 based finite element spaces while the local optimization problems are solved using the Newton method described in Sect. 4.2 (Ω = (0, 1) 3 ) All the numerical results reported below are related to Ω = (0, 1) 3 . On Table 10 we have reported the variations with respect to h of the L 2 (Ω) and H 1 (Ω) norms of the computed approximation error ψ h − ψ, for ψ defined by ψ(x, y, z) = e 1 2 (x 2 +y 2 +z 2 ) and ψ(x, y, z) = 1 2 x 2 + 5y 2 + 15z 2 , the related convergence orders, and the number of relaxation iterations necessary to achieve convergence. The local optimization problems are solved using the Newton method described in Sect. 4.2. As expected, nearly optimal orders of convergence are obtained for both the L 2 (Ω) and H 1 (Ω) norms of the computed approximation error. Both solutions exhibit comparable orders of convergence, however, the larger anisotropy of The space approximation relies on Q 1 based finite element spaces while the local optimization problems are solved using the Newton method described in Sect. 4.2 (Ω = (0, 1) 3 ) The space approximation relies on Q 1 based finite element spaces while the local optimization problems are solved using the Newton method described in Sect. 4.2 (Ω = (0, 1) 3 ) the second one implies a larger number of iterations for the relaxation algorithm to achieve its convergence. Approximation errors and iteration numbers are consistent with those reported in Sect. 7.2 for the same test problems. The next test problem we consider, is the one, already investigated in Sect. 7.3, whose exact solution ψ is given by ψ(x, y, z) = − R 2 − (x 2 + y 2 + z 2 ), with R ≥ √ 3, Ω still being the unit cube (0, 1) 3 . Tables 11 and 12 show that, for R = √ 6 and R = √ 3, the L 2 (Ω) and H 1 (Ω)-norms of the approximation error ψ h − ψ are nearly of optimal order; moreover, the above tables show that 3. The numerical results we have just reported show that, as long as accuracy and number of iterations are concerned, Q 1 based finite element approximations of problem (1) compared well with P 1 based ones if Ω is a cube and uniform structured partitions of Ω are used to define the finite element spaces. However the P 1 based methods can easily handle domains Ω of arbitrary shapes and unstructured finite element partitions, properties that the Q 1 based methods do not share.

Further Comments and Conclusions
In this article, we have discussed a least-squares/relaxation/mixed finite element methodology for the numerical solution of the Dirichlet problem for the three-dimensional elliptic Monge-Ampère equation det D 2 ψ = f (> 0) in Ω. The results reported in Sects. 7 and 8 show the robustness and flexibility of this methodology and its ability at approximating smooth convex solutions (if such solutions do exist) with nearly optimal orders of accuracy for the L 2 (Ω) and H 1 (Ω) norms of the approximation error ψ h − ψ. To the best of our knowledge, the above methodology is one of the very few which can solve, with nearly optimal orders of accuracy, the three-dimensional Monge-Ampère equation on domains Ω with curved boundaries using P 1 based finite element approximations.