Hervé V J Le Meur

Optimization under convexity constraints: are the finite element discretizations consistent ?

Keywords: convexity, finite elements, interpolation, conformal approximation, minimization subjclass:26B25, 52A41, 65K10, 65N30, 90C25, 91B24

 that the problem of minimizing a Dirichlet-like functional of the function u h discretized with P 1 Finite Elements, under the constraint that u h be convex, cannot converge. Here, we first improve this result by proving that non-convergence is due to the mesh refinment lack of richness, remains local and is true even for any mesh. Then, we investigate the consistency of various natural discretizations (P 1 and P 2 ) of second order constraints (subharmonicity and convexity) without discussing the convergence. We also numerically illustrate convergence of a method proposed in the literature that is simpler than existing methods.

INTRODUCTION

This paper is devoted to the numerical discretization of optimization with constraints on the second order derivative, namely problems of the form inf J(u),

D 2 u ∈ K, ( 1
)
where J is a functional and K is a subset of the set of symmetric matrices. Such problems appear in various contexts, in particular in physics and economics.

Since Newton, the shape of minimal resistance has been a topic of interest. It is called the Newton's problem and is of the type of [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF]. With some additional assumptions, it amounts to looking for a concave function that minimizes a nonlinear functional. See for instance the original book [2], the historical survey in [START_REF] Goldstine | A history of the calculus of variations from the 17th to the 19th century[END_REF], more recent theoretical results in [START_REF] Buttazzo | Minimum problems over sets of concave functions and related questions[END_REF][START_REF] Lachand-Robert | An Example of Non-convex Minimization and an Application to Newton's Problem of the Body of Least Resistance[END_REF][START_REF] Lions | Identification of the dual cone of convex functions and applications[END_REF] or numerical results [START_REF] Lachand-Robert | Minimizing within convex bodies using a convex hull method[END_REF][START_REF] Wachsmuth | The numerical solution of Newtons problem of least resistance[END_REF].

1

The associated problem of discretizing convex functions or bodies is indeed wider than expected. For instance, the Alexandrov's problem (see [START_REF] Carlier | On a theorem of Alexandrov[END_REF] and [START_REF] Lachand-Robert | Minimizing within convex bodies using a convex hull method[END_REF]), the Cheeger's constant [START_REF] Lachand-Robert | Minimizing within convex bodies using a convex hull method[END_REF][START_REF] Carlier | Approximation of maximal Cheeger sets by projection, M2AN Math[END_REF] and the Newton's problem [START_REF] Lachand-Robert | Minimizing within convex bodies using a convex hull method[END_REF][START_REF] Wachsmuth | The numerical solution of Newtons problem of least resistance[END_REF] can be numerically studied.

In economics, it suffices to remember that utility functions are concave to see the very wide applicability of these problems. For instance, in [START_REF] Rochet | Ironing, sweeping and multidimensionnal screening[END_REF], the authors are interested in finding the minimum of a convex and quadratic functional J over the set of convex functions min u∈K J(u) for J(u)

= Ω 1 2 ∇u T C ∇u -x • ∇u + (1 -α)u dx, (2)
where 0 ≤ α ≤ 1 and C is a (2,2) symmetric and positive-definite matrix and K is given by K = {u ∈ H 1 (Ω), u ≥ 0, u x ≥ 0, u y ≥ 0, u convex }.

(

) 3
So J is strictly convex and the set K is convex. It is then easy to check that there exists a unique minimizer of J over K. The regularity of the solutions is studied in [START_REF] Carlier | Regularity of solutions for some variational problems subject to a convexity constraint[END_REF]. Note that, when α = 1 and in a domain Ω with nonnegative coordinates, the problem degenerates: an exact convex solution exists, up to an additive constant

∇u(x) = C -1 x ⇒ u(x) = x C -1 x/2 -Cst, (4)
that can be fixed if we force u to vanish at a point (if α = 1).

In [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF], Choné and Le Meur exhibit an obstruction to the convergence of the mere discretization of this problem through conformal P 1 Finite Element (FE) method on regular meshes. They also extensively use the explicit solution (4) to illustrate the lack of convergence when a predicted condition is no more satisfied. Here "conformal" means that the discretized function is supposed to satisfy exactly the continuous constraint. As a consequence of this result, no optimization process that would use such conformal P 1 discretization of both the functional and the constraint may converge for any solution on these meshes. Yet approximation theory easily proves convergence of every separate discretization of these terms.

A possible approach to circumvent the mesh problem is, first, to test whether a sample of values on given points (and not mesh) may be associated to a convex function or body. Then one must construct the associated mesh (so, function dependent !) so as to interpolate. This was done in [START_REF] Leung | C 1 convexity-preserving interpolation of scattered data[END_REF] for C 1 FE by Leung and Renka. But such a regularity is too restrictive for us. This article [START_REF] Leung | C 1 convexity-preserving interpolation of scattered data[END_REF] reviews other papers, some of them being commented as false. It proves that the issue is not so simple.

The H 1 0 projection of a given function is addressed in [START_REF] Carlier | H 1 -projection into the set of convex functions: a saddle-point formulation[END_REF] through the saddle-point method. In this article, the authors even weight the convexity constraint to tune the convergence to the solution. The computations appear to be more robust than expected by theory. An attempt of explanation is given in [START_REF] Maury | Version continue de l'algorithme d'Uzawa. (French) [Continuous version of the Uzawa algorithm[END_REF].

In [START_REF] Carlier | A numerical approach to variational problems subject to convexity constraint[END_REF], the authors describe and implement an algorithm that optimizes not in the set of discretized and convex (i.e. conformal) functions but in the set of the convex functions after discretization (so that they may be non-convex once discretized). More precisely they characterize, for a specific structured set of cartesian points, the image through the P 1 discretization of a (continuous) convex function. This yields such a huge number of constraints on the function that it may not be recommended.

In [START_REF] Lachand-Robert | Minimizing within convex bodies using a convex hull method[END_REF], Lachand-Robert and Oudet address the problem of discretizing a convex body. They use the parametric generalization in order to discretize convex functions' graphs. In the functional, they isolate the dependence on the point x, the unit normal ν at x and the signed distance φ = ν.x. Although they notice that the three variables are "somehow redundant", they implement a gradient like method, based on the variation of only φ = x.ν. They address both the Alexandrov's problem, the Cheeger's sets and the Newton's problem. This was revisited since in [START_REF] Carlier | Approximation of maximal Cheeger sets by projection, M2AN Math[END_REF] and solved in a very elegant way in [START_REF] Wachsmuth | The numerical solution of Newtons problem of least resistance[END_REF].

Separately, Aguilera and Morin [START_REF] Aguilera | Approximating optimization problems over convex functions[END_REF] prove convergence of a Finite Difference (FD) "approximation using positive semidefinite programs and discrete Hessians". In [START_REF] Aguilera | On convex functions and the finite element method[END_REF], the same authors prove even convergence of a Finite Element (FE) discretization of the weak Hessian. Since FD are included in FE for convenient meshes, these two articles seem to contradict both [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] and the present article. Indeed, they do not, but we discuss them below (section 5).

In [START_REF] Ekeland | An algorithm for computing solutions of variational problems with global convexity constraints[END_REF], Ekeland and Moreno propose a non-local discretization that relies on the representation of any convex function as the supremum of affine functions (its minorants). So their "discrete" representation of continuous functions is conformal but non-local. The major drawback is that the complexity increases drastically with dimension, due to the nonlocality, but it works.

More recently, two very different means of resolution were proposed. Both are nonconformal in the sense that the discretized solution is not convex. In the first one, Mérigot and Oudet [START_REF] Mérigot | Handling convexity-like constraints in variational problems[END_REF] propose to discretize the convexity constraint by under sampling it. More precisely, they prove convergence of their algorithm if the constraint is forced only on a subset of the sampling points. In the second one, Mirebeau [START_REF] Mirebeau | Adaptive, Anisotropic and Hierarchical cones of Discrete Convex functions[END_REF] proves that for a given grid, he may provide a sequence of sets of (less than four) points on which it suffices to force convexity to ensure convexity at convergence on this given grid. This makes only O(N log 2 N ) linear constraints.

There seems to be a convergent opinion through [START_REF] Mirebeau | Adaptive, Anisotropic and Hierarchical cones of Discrete Convex functions[END_REF], [START_REF] Mérigot | Handling convexity-like constraints in variational problems[END_REF], [START_REF] Aguilera | On convex functions and the finite element method[END_REF] that indeed the number of constraints must be not too large. A deeper discussion of this crucial point is postponed to a forthcoming article.

In [START_REF] Wachsmuth | Conforming Approximation of Convex Functions with the Finite Element Method[END_REF], the author goes further than [START_REF] Aguilera | On convex functions and the finite element method[END_REF]. He uses approximation theory, where there is no reason why a function (should it be convex or not) should not be approximable. So as to remain inside this theory, he uses two auxilliary functions that make his solution strictly convex, so that the constraint of convexity may not be saturated. Nevertheless, we are surprised that these functions are not needed in his numerical experiment. It could be due to the fact that the hidden constraints (that will appear below) are filled by the exact solution and so do not trigger any discrepancy. This explains why our results are not contradictory with Oberman's or Wachsmuth's.

A totally different tool has emerged since. It consists in using an evolution PDE with a u 0 as an initial condition. In [START_REF] Oberman | Computing the convex envelope using a nonlinear partial differential equation[END_REF] (and other articles), the author computes the convex envelope of a graph using such a nonlinear PDE. Later [START_REF] Carlier | Exponential convergence for a convexifying equation[END_REF] exhibits an evolution PDE, close to the previous, which, starting from a given u 0 , is proved to converge to the convexified of u 0 . They use stochastic control representation to prove an exponential convergence. This set of evolution PDE is a new idea and very important in optimal transport theory.

The present article deals with two main types of non-convergence results. On the one hand, the non-convergence in L 2 for conformal P 1 Finite Element (FE) is revisited after [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] from a theoretical point of view. On the other hand, the consistency of discretized linear optimization over second order constraints is investigated for various discretizations (P 1 and P 2 ) and various constraints: linear (subharmonicity) or nonlinear (convexity). The motivation for convexity has been stressed. Subharmonicity is only studied as an intermediate step toward convexity because it is linear.

In Section 2, first we restate some already known results, then we prove that nonconvergence is purely local and so is true for any mesh. In Section 3, we investigate the consistency of various P 1 FE discretizations of our model problem. Section 4 is devoted to the consistency of P 2 FE discretizations (strong or weak convexity). We discuss the articles [START_REF] Aguilera | Approximating optimization problems over convex functions[END_REF][START_REF] Aguilera | On convex functions and the finite element method[END_REF] that could seem to contradict [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] in Section 5 and we conclude in Section 6.

THEORETICAL RESULTS ON THE P 1 FEM

In the present section, we first recall some already known results and make them more explicit.

Already known results

In [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF], various results are proved. Since the goal in the present subsection is to extend this study, we need to remind the reader of these results. We consider an open convex and bounded domain Ω of R 2 , but we use only Ω = (a, b) 2 (and a and b positive) in applications. The Lemma 3 of [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] states: Lemma 1 A function u h , P 1 in every triangle of Ω's mesh, is convex if and only if, for any pair of adjacent triangles

(q 2 -q 1 ) • n ≥ 0, (5)
where q 1 (resp. q 2 ) is the (constant) gradient of u h inside triangle 1 (resp. 2) and n is the unit normal pointing from triangle 1 to 2.

As reminded in [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF], convexity of u can be defined dually as -∇u ⊗ ∇v 0, where A 0 means the matrix A is positive semi-definite, v is a nonnegative test function in

C ∞ 0 and -∇u ⊗ ∇v = - u x , v x u x , v y u y , v x u y , v y . ( 6
)
This definition is weak (in the sense of distributions) but it is also equivalent to a strong one when u is sufficiently regular. We will use a weaker definition where the trial and test functions will be in a (finite dimensionnal) subspace of C ∞ 0 : the set of some FE functions.

Last we will report on [START_REF] Aguilera | On convex functions and the finite element method[END_REF] which uses an even weaker definition where the test functions are in a much smaller subspace than u h 's. One can easily check that the weak definition of convexity (with C ∞ 0 test functions) implies the weaker definition of convexity (with only the basis functions in P 1 associated to interior points as test functions) but they are not equivalent.

The proof of Lemma 1 is based on the following formula written for the P 1 function

u h : ∂ 2 u h ∂a∂b , ϕ = e (q 2 -q 1 ) • n (n • a)(n • b) e ϕ(s)ds, (7)
where the summation is taken over all interior edges e of the mesh, a and b are two unit vectors, ϕ is a C ∞ function with compact support in Ω, q 1 and q 2 designate the two (constant) gradients of u h in the two triangles that share the edge e and n is the unit normal from triangle 1 to triangle 2. This formula is right thanks to the property that the gradient of a continuous function across the edge has its tangential derivative along the edge continuous. Namely (q 2 -q 1 ) • t = 0 where t is the unit tangent vector. The proof also relies on the fact that for all unit vector n, the bilinear form (a, b)

→ (n • a)(n • b) is positive semi-definite, whatever n.
Through a proof very similar to the one of Theorem 4 in [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF], one may prove a wider Theorem that even applies locally:

then any function u h convex and P 1 will satisfy the following equation in the sense of distributions on Ω :

∂ 2 u h ∂a∂b ≥ 0. (9)
The proof relies on the fact that if a P 1 function u h is convex, then by Lemma 1, the gradient's jumps are non-negative. As a consequence, in Formula (7), if u h is convex, the scalar coefficients (q 2 -q 1 ) • n are all non-negative. But separately, if condition (P M ) is satisfied, (n•a)(n•b) ≥ 0 and so even the non-diagonal terms of the hessian (∂ 2 u h /(∂a∂b)) are sign-constrained ! Notice first that property (P M ) only depends on the geometry of the mesh and not on any function. Notice then that the property (P M ) is at given h but, if it remains for a sequence of h → 0 in Ω ⊂ Ω, then the associated functions u h given by Theorem 1 satisfy (9) even at the limit in Ω ⊂ Ω. Yet such a property ( 9) is (sometimes) contradictory with the property of convexity. For instance (x, y) → (x 2 -2ρxy + µy 2 )/(µ -ρ 2 ), where µ -ρ 2 > 0 and 1 + µ > 0, is strictly convex. Yet ∂ 2 u/(∂x∂y) = -ρ/(µ -ρ 2 ) has no sign prescribed if u is only strictly convex. This is the key of the obstruction to the convergence stated in [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF]. The next subsection is devoted to discussing the generality of (P M ).

Is property (P M ) frequent ?

In this subsection, we start from some elementary observations on particular meshes, then we argue on whether such meshes are frequent and state a Theorem. Three types of meshes may be considered that are depicted in Figure 1.

Concerning mesh 1, three different unit normals (up to a multiplicative factor -1) exist in all the domain:

{(0, 1); (1, 0); (-1/ √ 2, 1/ √ 2)}. If we choose a = (-1, 0) and b = (0, 1), the values taken by (n • a)(n • b) (with the possible n) are {0; 0; (1/ √ 2)(1/ √ 2)}. They are all nonnegative.
Concerning mesh 2, four types of unit normals can be found in all the domain:

{(0, 1); (1, 0); (-1/ √ 2, 1/ √ 2); (1/ √ 2, 1/ √ 2)}. If we choose a = (1, 0) and b = (1/ √ 2, -1/ √ 2), the values taken by (n • a)(n • b) (with the possible n) are {1/ √ 2; 0; (-1/ √ 2)(-1); 0}. They are all nonnegative.
Any structured mesh like mesh 1 or 2, if refined while keeping the same structure, will satisfy (P M ) with the very same a and b for any h. But what can be stated about a more general mesh like mesh 3 ?

We are going to prove the following Theorem.

Theorem 2 Let Ω be an open and convex set of R 2 , Ω ⊂ Ω an open and convex subset and T h a triangulation of Ω.

For any given h and Ω , there exists (a, b) such that (P M ) is true in Ω and for T h . Moreover, if the refinement process does not enrich the edges' directions of T hn in Ω , then the property (P M ) is true for any T hn in Ω .

Proof of Theorem 2

First, one must notice that the condition (n • a)(n • b) ≥ 0 is invariant through changing n into -n. So, it does not depend on the choice of the hyperplane's unit normal.

Let us take a very general mesh like mesh 3 for which we may isolate a single triangle. Then, for each of the three unit normals (up to the multiplicative coefficient -1), having (n • a)(n • b) ≥ 0 is equivalent to having a and b in the same closed half plane whose normal is n. Equivalently, we can take either a and b on the same side as n or on the opposite side. Gathering the conditions associated to each of the three edges, we are led to choosing both a and b in one of the six half cone that intersects none of the three hyperplanes. Such a configuration is depicted on Figure (2). The three unit normals are drawn and called n1, n2, n3 and their associated hyperplanes are called P1, P2, P3. So, one may easily find such a couple (a, b) for any given triangle.

Even better, whatever might be the finite number of edges for a larger Ω and a given T h , one sees obviously that such a choice of a and b is easy since there always exists a finite number of hyperplanes. As a consequence, a and b can be choosen to be both in any of the cones which partition the whole unit circle. So the property (P M ) is true even on any given and potentially unstructured mesh of a subdomain of Ω.

Could a refinement process preserve the property (P M ) for a given couple (a, b), a given Ω and any h → 0 ?

In case of any general given mesh T h we have just proved that one may find a and b, should they suit the property only in Ω . Then, for instance, if the refinement process is such that every triangle in Ω is refined into four homothetic subtriangles, it is obvious that no more edges' direction will be provided. As a consequence, the very same a and b will then suit property (P M ) in Ω for any refined mesh and any h. The proof is complete.

In a sense to be defined, the set of refinement processes that enable (P M ) is open and non-empty. The main question is "Does the refinement process enrich in normals n ?".

More information on the lack of convergence

In [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF], the authors state that the "conformal method may not converge for some limit function" because the second derivative of the limit is forced to satisfy an unnatural condition. In this subsection, we give a more precise result in L 2 and still use our Lemma 1 that enables to identify the convexity of u h and the non-negativity of its gradients' jumps (for P 1 FE). This enables us to state our main Theorem: Theorem 3 Let Ω be an open and convex domain in R 2 , Ω an open and convex subdomain of Ω and a family of meshes (T h ) h→0 + . If the property (P M ) is satisfied in Ω and for the meshes (T h ) h→0 , then there exists ε > 0 and a C ∞ convex function u exact such that

min u h ∈P 1 and u h convex | u h -u exact | L 2 ≥ ε.
This Theorem uses that if Property (P M ) is satisfied in a domain Ω for a family of meshes (T h ) h→0 + , should it be only locally in Ω , there will be two unit vectors a, b such that (9) holds. Given those a, b, there exists a function u exact that may not be the limit in L 2 of any sequence of functions both P 1 and convex (in a strong definition) on this family of meshes. The improvement of this Theorem with respect to [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] is that the obstruction is local. As a consequence we may apply our obstruction on any mesh and not only regular ones.

We need some more preliminary results and definitions before proving this Theorem 3 concerning P 1 FE.

Definition 1

Let Ω be an open subset of R 2 , and u ∈ C 0 (Ω). For any x = (x, y) ∈ Ω and {a, b} two independent unit vectors, there exist two positive numbers {α 0 , β 0 } such that

∀α, 0 ≤ α ≤ α 0 , ∀β, 0 ≤ β ≤ β 0 , x + αa + βb = (x + αa 1 + βb 1 , y + αa 2 + βb 2 ) ∈ Ω.
For such (α 0 , β 0 , a, b), one defines:

φ (α 0 ,β 0 ,a,b) (u) = (u(x + α 0 a + β 0 b) -u(x + α 0 a) -u(x + β 0 b) + u(x)) /(α 0 β 0 ). If α 0 = β 0 → 0 + , then φ (α 0 ,β 0 ,a,b) (u) → ∂ 2 u/(∂a∂b).
We are going to prove that such a quantity φ (α 0 ,β 0 ,a,b) (u) will be overconstrained by the mere discretization and the limit u will not satisfy the right sign of the second derivative.

We want now to define an explicit solution depending on some parameters. A wellchosen combination of these parameters will provide a non-approximable function. For instance, let Ω = (1, 2) 2 and the problem (2, 3) for a positive definite real matrix C be such that:

C = µ 2 ρ ρ µ 1 . (10)
Then there exists an exact solution thanks to (4):

u exact = 1 µ 1 µ 2 -ρ 2 (µ 1 (x 2 -1)/2 + µ 2 (y 2 -1)/2 -ρ(xy -1)) in Ω. ( 11
)
The function u exact is such that it is zero at the corner (1, 1) of the domain Ω chosen but it can easily be generalized to other domains. Simple computations prove the following formula:

φ (α 0 ,β 0 ,a,b) (u exact ) = 1 µ 1 µ 2 -ρ 2 (µ 1 a 1 b 1 + µ 2 a 2 b 2 -ρ(a 1 b 2 + a 2 b 1 )) = a C -1 b, (12)
where a 1 , a 2 , b 1 , b 2 are the components of a, b. It is then useful to state the following Lemma.

Lemma 2 Let a, b be two given independent unit vectors in R 2 . Let Ω ⊂ R 2 and η > 0 given. Then, there exists

(µ 1 , µ 2 , ρ) such that µ 1 ≥ 0, µ 2 ≥ 0, µ 1 µ 2 -ρ 2 ≥ 0 and φ (α 0 ,β 0 ,a,b) (u exact ) = a T C -1 b ≤ -η in Ω, (13)
provided (α 0 , β 0 ) ∈ R + * 2 are such that there exists φ (α 0 ,β 0 ,a,b) (u exact ).

Roughly speaking, Lemma 2 states that for any independent a, b, one may find a positive definite matrix C -1 such that the associated u exact satisfies the following property on its associated bilinear form:

φ (α 0 ,β 0 ,a,b) (u exact ) = a T C -1 b < 0.

Proof of Lemma 2

From the formula [START_REF] Carlier | Regularity of solutions for some variational problems subject to a convexity constraint[END_REF], one sees that it is sufficient to find a positive definite symmetric bilinear form (C -1 ) such that, for independent a, b given, a T C -1 b < 0. We will build up C -1 from its eigenvalues and eigenvectors.

Let e 1 be a unit vector between a and b (normalized mean for instance). Then let e 2 be a unit vector normal to e 1 . With the same notations for a and b in this new basis, a 1 b 1 > 0 and a 2 b 2 < 0.

Let now C -1 be a matrix in the canonical basis with the eigenvectors e 1 , e 2 and the associated positive eigenvalues λ 1 , λ 2 . Obviously, C is positive semi-definite. Then

a T C -1 b = λ 1 a 1 b 1 + λ 2 a 2 b 2
which may be less than -η given, for an appropriate choice of λ 1 , λ 2 . Then the (µ 1 , µ 2 , ρ) are the coefficients of the matrix C in the canonical basis as written in [START_REF] Carlier | Approximation of maximal Cheeger sets by projection, M2AN Math[END_REF].

We have now obtained all what is needed to start the following proof.

Proof of Theorem 3. Since we assume (P M ) is satisfied in Ω for (T h ) h→0 , this provides a couple of vectors (a, b). Theorem 1 enables us to claim that any P 1 function u h satisfies (9) in the sense of distributions in Ω :

0 ≤ ∂ 2 u h ∂a∂b , Ψ Ω = u h , ∂ 2 Ψ ∂a∂b Ω , (14)
for any nonnegative Ψ ∈ C ∞ 0 (Ω ). This property obviously remains for open and convex subdomains of Ω .

Should it be needed, one could decrease α 0 , β 0 > 0 and find a subdomain Ω ⊂ Ω such that for any nonnegative ϕ ∈ C ∞ 0 (Ω ), the function

Ψ : x → 1 0 1 0 ϕ(x -tα 0 a -t β 0 b) dt dt , (15)
is nonnegative and in C ∞ 0 (Ω ). So the function Ψ is eligible for equation ( 14) whatever φ. The second derivative of Ψ may be computed:

∂ 2 Ψ ∂a∂b = (ϕ(x -α 0 a -β 0 b) -ϕ(x -β 0 b) -ϕ(x -α 0 a) + ϕ(x)) /(α 0 β 0 ) = φ (α 0 ,β 0 ,-a,-b) (ϕ).
So we have, for any nonnegative ϕ ∈ C ∞ 0 (Ω ):

0 ≤ ∂ 2 u h ∂a∂b , Ψ Ω = u h , ∂ 2 Ψ ∂a∂b Ω = u h , φ (α 0 ,β 0 ,-a,-b) (ϕ) Ω = φ (α 0 ,β 0 ,a,b) (u h ), ϕ Ω . (16)
Separately, given a, b, η > 0, Lemma 2 enables us to claim there exists a convex quadratic function u exact such that (13) is true for η > 0 given. Moreover, because of [START_REF] Carlier | A numerical approach to variational problems subject to convexity constraint[END_REF], since the open set Ω ⊂ Ω is of non-zero measure, and φ (α 0 ,β 0 ,a,b) (u exact ) is a constant known by [START_REF] Carlier | Regularity of solutions for some variational problems subject to a convexity constraint[END_REF], we have

φ (α 0 ,β 0 ,a,b) (u h -u exact ), ϕ Ω ≥ -φ (α 0 ,β 0 ,a,b) (u exact ) Ω ϕ.
As a consequence, for convenient Ω , we are led to

φ (α 0 ,β 0 ,a,b) (u h -u exact ), ϕ Ω ≥ η Ω ϕ, (17)
for any nonnegative ϕ ∈ C ∞ 0 (Ω ) and u h convex and P 1 . This last inequality contradicts any possible convergence in L 2 (Ω ) of a sequence u h of convex functions P 1 in (T h ) h to u exact given by Lemma 2.

The previous Theorem adds one more argument to the need of convenient numerical methods to discretize the constraint of convexity.

3 Consistency of the P 1 FEM

The non-convergence of the convexity problem (2, 3) discretized by conformal P 1 Finite Elements is proved and numerically illustrated in [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF]. It is even proved to apply to nonregular meshes in the subsection 2.3. We give here a very different argument by studying the consistency of various discretizations of convexity.

One could wonder why study consistency. The reason is that when one discretizes convexity (or even its linear form of subharmonicity), there exist meshes and FE on which the discretization is not consistent. As a consequence, any proof of convergence must exclude at least these cases. The fact that the FE-discretized Dirichlet functional and the FE-discretized Hessian matrix converge separately is not sufficient. Already [START_REF] Kawohl | Convergent finite elements for a class of nonconvex variationnal problems[END_REF] claimed convergence, but was later proved to be wrong in [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF].

Use of the gradients' jumps for convexity

Given a sufficiently smooth function u, one may interpolate it in P 1 FE as u h = N i=1 u hi φ i (x). If the mesh is structured as in Figure 3, should it be local, one may compute the gradients' jumps accross the edges as functions of the values at the various involved nodes. Then one may compute the series expansion of the sampled values u i from the exact initial function u. The jumps between triangles 1 and 2, 2 and 3 and last 3 and 4 are respectively:

Jump(1/2) = +u 6 -u 7 -u 1 +u 2 h = (u xx + u xy ) x=x 1 ,y=y 1 h + O(h 2 ), Jump(2/3) = -u 1 +u 7 -u 2 +u 3 h = (u xy + u yy ) x=x 1 ,y=y 1 h + O(h 2 ), Jump(3/4) = √ 2 -u 1 +u 2 -u 3 +u 4 h = - √ 2(u xy ) x=x 1 ,y=y 1 h + O(h 2 ), (18)
where we denote u x the partial derivative with respect to x and u xx the second order derivative with respect to x. When h → 0, instead of forcing the solution of a problem to be convex, we force it to some mesh-dependent combination of its second order derivatives to have a sign or another. In addition, this combination is meaningless since for the convex limit function, it may have whatever sign. We will say that such a discretization is not consistent.

We may repeat here that our non-consistency result is not contradictory with the convergence results of [START_REF] Wachsmuth | Conforming Approximation of Convex Functions with the Finite Element Method[END_REF] (for instance) since these results rely on a regularization that makes the algorithm look for strictly convex functions. Looking in the interior of the set of convex functions, the constraint may not be saturated and so disappears.

Below, we investigate the consistency of various other discretizations.

Use of a weak P 1 version

We use here a weak P 1 definition of convexity which is similar to the one of [START_REF] Aguilera | On convex functions and the finite element method[END_REF], except that their trial basis and test basis functions are different. We use the same trial and test P 1 basis. Their method is fully discussed and illustrated in Section 5.

The subharmonicity constraint

Subharmonicity is an interesting property simpler than convexity since it is only linear. A weak definition of subharmonicity (∆u ≥ 0) is, for any test function φ i in the discrete basis (⊗ is defined in ( 6)):

∆u h , φ i = Tr D 2 u h , φ i = -Tr ∇u h ⊗ ∇φ i ≥ 0. ( 19
)
One may then prove the following Proposition:

Proposition 3.1 The weak P 1 discretization of the subharmonicity constraint on a mesh like in Figure 3 is consistent:

∆u h , φ i = -4u 1 + u 2 + u 5 + u 7 + u 4 = (u xx + u yy ) x=x 1 ,y=y 1 h 2 + O(h 3 ), ( 20
)
where i is the node at the center of the cell's group (x 1 , y 1 ).

The proof is very easy and left to the reader. Indeed, it is only the discretization of the Laplacian which is known to be consistent and even convergent ! In order to test this discretization, we used the Matlab package optim to minimize the functional Ω |∇u| 2 /2 + ∇u exact • ∇u, where f = ∆u exact , over the set of subharmonic functions. The initial condition is x(x -1)y(y -1). The exact solution of this H 1 0 projection is u exact = x 2 y 2 -(x 4 + y 4 )/6 and its Laplacian vanishes. The convergence with the procedure quadprog of quadratic programming can be seen on Figure 4. It is quite satisfactory.

The convexity constraint

One may give a weak definition of convexity:

Tr D 2 u h , φ i ≥ 0 and det D 2 u h , φ i ≥ 0, (21)
for any φ i , basis function of the FEM at the node i. One may then prove the following Proposition:

Proposition 3.2 The weak P 1 discretization of the convexity constraint (21) on a mesh as in Figure 3 is consistent as can be seen from ( 20) and (⊗ is defined in ( 6)):

det D 2 u h , φ i = det ∇u h ⊗ ∇φ i = = (-2u 1 + u 2 + u 5 )(u 7 + u 4 -2u 1 ) -(-u 7 + u 3 -u 4 + u 6 -u 5 -u 2 + 2u 1 ) 2 /4 = (u xx u yy -u xy u yx ) x=x 1 ,y=y 1 h 4 + O(h 5 ). ( 22
)
where i is the node at the center of the cell's group (point of coordinates (x 1 , y 1 ). This discretization [START_REF] Mérigot | Handling convexity-like constraints in variational problems[END_REF] for the (nonlinear) convexity contraint is consistent. Nevertheless, the numerical treatment of the linear and nonlinear constraints should raise inaccuracies since they are of very different orders of magnitude (O(h 2 ) and O(h 4 )). A very natural question would be to investigate its convergence. Such a study is postponed to an other article.

The consistency of the P 1 FEM is tested for various meshes. The weak P 1 Hessian is reported in Table 1 where one may claim that two meshes are non-consistent with the determinant but all are consistent with the trace (laplacian).

Consistency of the P 2 FEM

In the present section, we investigate the consistency of various discretizations through P 2 FEM of two second order derivative constraints: subharmonicity and convexity.

First, we interpolate a continuous function u to a P 2 function u h = i u i φ i (x) in a domain Ω meshed with triangles as in Figure 5. Here, the index i denotes both vertices and edge midpoints. Then we compute the discretized version of the second order term constrained to be nonnegative (various versions are treated) and compute its series expansion. If forcing it to be nonnegative amounts to forcing the continuous limit function to the correct constraint, then we claim the discretization is consistent. Otherwise it is inconsistent.

For the whole section, we assume the mesh is (at least locally) structured around the point of coordinates (x 1 , y 1 ). The local numbering of triangles is depicted in Figure 5. The node (x 1 , y 1 ) is localy numbered 1 in every triangle and the local numbering of nodes is depicted in triangle 3 of Figure 5.

Gradients' jumps for the convexity constraint

In a way similar to the P 1 case, one may prove for P 2 functions u h :

∂ 2 u h ∂a∂b , ϕ = e e (q 2 (s) -q 1 (s)) • n ϕ(s)ds (n • a)(n • b) + K ∂ 2 u h ∂a∂b | K K ϕ,
for any e interior edge of the mesh, a and b are two unit normal vectors, ϕ is a C ∞ 0 function with compact support in Ω. By taking ϕ localized along the edge e, one may state that this definition of convexity (with C ∞ 0 test functions) implies the non-negativity of the gradients' jumps.

In order to test the consistency of such a discretization, we compute the gradients' jumps accross the edges common to triangles 1 and 2, 2 and 3 and last 3 and 4. Of course, they are not constant since they are P 1 FE. After an exact computation and a series expansion (details left to the reader), one may state the following Proposition. Proposition 4.1 The discretization of the convexity constraint with the jump of the gradients between triangles 1 and 2, 2 and 3, 3 and 4 of a P 2 function on a mesh like in Figure 5 gives terms

Jump(1/2) = (u xxy + u xyy ) x=x 1 ,y=y 1 (y -y 1 + h/2)h/2 + O(h 3 ); y ∈ [y 1 -h, y 1 ] Jump(2/3) = (u xxy + u xyy ) x=x 1 ,y=y 1 (x -x 1 -h/2)h/2 + O(h 3 ); x ∈ [x 1 , x 1 + h] Jump(3/4) = -(u xxy + u xyy ) x=x 1 ,y=y 1 (x -x 1 -h/2) √ 2h/2 + O(h 3 ); x ∈ [x 1 , x 1 + h]. (23) Such a discretization is non-consistent.
Since for instance between triangles 1 and 2, (y -y 1 + h/2) changes sign (but remains O(h)), one deduces from ( 23) that forcing the non-negativity of the gradients' jump along the whole edge forces the limit function to satisfy the non-natural equality condition u xxy + u xyy = 0 ! So the weak C ∞ 0 (or strong !) definition of convexity implies the non-negativity of gradient's jumps which is non-consistent and so must be rejected. We check below that the weak P 2 definition may (subsection 4.3) or maynot be consistent (subsection 4.2).

Weak P 2 version of the second derivative at a vertex

We define the weak P 2 -P 2 (i.e. with P 2 trial and test functions) version of the Hessian matrix at vertex i as:

D 2 u h , φ i = -∇u h ⊗ ∇φ i , ( 24
)
where ⊗ is defined in [START_REF] Lions | Identification of the dual cone of convex functions and applications[END_REF]. Strictly speaking, this may not provide a correct weak version of non-negativity since φ i changes sign. Anyway, one may state the following Proposition which proof is left to the reader:

Proposition 4.2
The discretization of the linear part of the convexity constraint according to (24) of a P 2 function on a mesh as in Figure 5 (h = ∆x = ∆y) and a P 2 test function centered at the vertex (x 1 , y 1 ) gives:

Tr ( D 2 u h , φ i ) = -(u xxxx + u yyyy ) x=x 1 ,y=y 1 * h 4 /48 + O(h 5 ). ( 25
)
Such a discretization of an inequality constraint is non-consistent.

It appears that such a discretization is not even consistent for the linear part of the constraint. More precisely, while one could believe one forces the solution to be subharmonic, indeed, one forces it to be such that u xxxx + u yyyy ≤ 0 ! The full nonlinear convexity constraint on the same vertices may only work worse.

One must notice that the fourth order of the expansion in ( 25) is meaningful. Indeed, on the one hand the three midpoints quadrature is exact for P 2 functions in a triangle. On the other hand the basis function for a vertex vanishes on these midpoints. As a result, the order two term is identically zero. So the first non-zero term is the fourth one and the constraint is of fourth order too.

Weak P 2 version of the second derivative at an edge midpoint

We define the weak P 2 version of the second derivative at an edge midpoint (denoted by index j) in a way similar to [START_REF] Carlier | Exponential convergence for a convexifying equation[END_REF]. Indeed, the function φ i in ( 24) is replaced by the P 2 basis function φ j associated to an edge midpoint indexed by j.

As in subsection 3.2, the basis functions (of the edge midpoints) are non-negative. So it is a priori an admissible weak formulation of semi-definite positiveness.

The subharmonicity constraint

Let us assume we discretize the constraint ∆u ≥ 0 by

Ω ∇u h • ∇φ j ≤ 0, (26)
for all j index of an edge midpoint interior to Ω. One may then state a Proposition (which proof is left to the reader):

Proposition 4.3 The discretization of the linear part of the convexity constraint according to (26) on a mesh as in Figure 5 (h = ∆x = ∆y) gives the same series expansion, whether the edge is vertical, horizontal or diagonal:

Tr ( D 2 u h , φ j ) = (∆u) x=x 1 ,y=y 1 h 2 /3 + O(h 3 ), (27)
for any j index of an edge midpoint interior to Ω. Such a discretization of an inequality constraint is consistent.

Various questions remain. If we discretize the subharmonicity (or the convexity) constraint only at edge midpoints, is it enough constraints or not ? More generaly, the amount of constraints compared to the amount of degrees of freedom (dof) should be discussed. The relevance of this question seems to be the conclusion of most recent articles like [START_REF] Mirebeau | Adaptive, Anisotropic and Hierarchical cones of Discrete Convex functions[END_REF][START_REF] Mérigot | Handling convexity-like constraints in variational problems[END_REF][START_REF] Aguilera | On convex functions and the finite element method[END_REF].

The convexity constraint

Like in the subharmonic case, we take a weak P 2 version of the continuous nonlinear constraint det D 2 u ≥ 0 with (nonnegative) test functions associated to every edge midpoint j in the interior. One may then easily prove the following Proposition: Proposition 4. [START_REF] Buttazzo | Minimum problems over sets of concave functions and related questions[END_REF] The discretization of the nonlinear part of the convexity constraint on a mesh as in Figure 5 (h = ∆x = ∆y) gives:

det( D 2 u h , φ j ) = (u xx u yy -u xy u yx ) x=x 1 ,y=y 1 h 4 /9 + O(h 5 ), (28)
for any j index of an edge midpoint interior to Ω. Such a discretization of an inequality constraint is consistent.

So the convergence order of the linear constraint is two while the nonlinear one's is four. Such a discrepancy between the orders of convergence of those two constraints should be managed in numerical simulation using both constraints. At least this discretization is consistent.

Discussion on the Aguilera and Morin's articles

Let us notice that [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] seems not to have been known of the authors of [START_REF] Aguilera | Approximating optimization problems over convex functions[END_REF][START_REF] Aguilera | On convex functions and the finite element method[END_REF].

Roughly speaking, the first article [START_REF] Aguilera | Approximating optimization problems over convex functions[END_REF] proves convergence of the FD discretization of problems like ours. The second [START_REF] Aguilera | On convex functions and the finite element method[END_REF] proves convergence of the FE discretization of the same problems but the authors imagine the clever trick of using different basis for the trial and the test functions.

Finite Differences

In [START_REF] Aguilera | Approximating optimization problems over convex functions[END_REF], Aguilera and Morin deal with FD and "prove convergence under very general conditions, even when the continuous solution is not smooth". Their proof relies on approximation theory and on a weak definition of convexity, namely the FD approximation of the Hessian matrix is forced to be positive. They stress that their FD-convexity is not equivalent to a continuous convexity. As a consequence, their approximation is not conformal. So there is no contradiction with [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF]. In addition, one must notice that the number of constraints is roughly twice the number of interior points. It is large, but grows only linearly with the degrees of freedom. They also provide numerical experiments.

Unfortunately the numerical experiments are not very convincing. Concerning the monopolist problem in 2D, they claim that "the error in the L ∞ norm is smaller or approximately equal to h" (p. 27). Yet, their Table 1 provides errors that we draw in a loglog scale in the left part of Figure 6. The order of convergence is not clear.

For the monopolist problem in 3D, the authors notice that "the L ∞ error is not converging to zero with order O(h)" (p. 29). Since they only prove convergence and do not forecast any order, there is no contradiction. But when their Table 2 is drawn in a loglog scale (see Figure 6 right), convergence seems not to be reached yet. Moreover the execution time grows faster than polynomialy both for 2D and 3D monopolist.

As reminded by the authors, FD discretization gives the very same matrices as the P 1 FE (both for trial and test functions), on a mesh like our mesh 1 on Figure 1 (except boundary conditions irrelevant here). In this sense, one may claim that their FD discretization is somehow equivalent to a P 1 (for trial) -P 1 (for test) discretization which was proved to be non-consistent. Indeed, the next article [START_REF] Aguilera | On convex functions and the finite element method[END_REF] gives one example (Example 3.7 p. 3150) of P 1 discretization and a weak definition of convexity (with P 1 trial and test functions, so alike FD discretization), where " (u h ) does not converge to u as h → 0" (p. 3151 and Figure 3.1).

Discussing the way these results on FD are coherent with the results on FE is postponed to a forthcoming article.

Finite Elements

In [START_REF] Aguilera | On convex functions and the finite element method[END_REF], Aguilera and Morin prove convergence of the discretization of the full problem: not only the approximation of the Hessian, but also of the functional together.

Especially, they define u h , interpolation of u on a given mesh with respect to the trial functions basis {φ h r } r as u h (x) = r∈I h trial u hr φ h r (x). They also use a (possibly) different FE basis for the test functions {ϕ h s } s∈I h test . Then, the discrete Hessian matrix is defined weakly by H h s u h = -∂ i u h , ∂ j ϕ h s i,j and FE convexity is defined as:

H h s u h = -   r∈I h trial u hr ∂ i φ h r , ∂ j ϕ h s   i,j 0,
for all ϕ h s in the test functions basis, where they denote A 0 if A is positive semidefinite. Notice that the matrix ( ∂ i φ h r , ∂ j ϕ h s ) r,s is a rectangle if the trial and test basis do not have the same number of functions ! Then their main result states convergence of the discretization once one assumes at least (p. 3147): C.2 There exists a linear operator I h with values in the discretization space V h (the interpolant), an integer m ≥ 2, and a constant C independent of u and h such that

u -I h u H 1 (Ω) ≤ Ch m-1 u H m (Ω) ,
C.3 The basis test-functions are such that ϕ h s (x) ≥ 0 for all x ∈ Ω, all h > 0, and all s.

The condition C.3 merely states that the test functions, that are supposed to approximate positive C ∞ 0 functions, remain nonnegative. So the basis functions are at most P 1 or a subset of P k .

The condition C.2 is very classical and states that the trial space is P m-1 . In addition, the proof requires m > 2. In other words, the trial functions must at least be P 2 . This is seen by the authors who say they "need to elaborate" on the condition m > 2. On their p. 3150, they notice the FE discretized Hessian (P 1 -P 1 and so m = 2) expands into the FD scheme (for one specific mesh !) at first order and so they refer to their FD article to conclude that in the case of some specific meshes, their Theorem "also holds for m = 2" (or P 1 -P 1 FE).

Yet, the very next example they provide is P 1 -P 1 and they "report numerical evidence supporting the necessity of assuming m > 2 in [their] Theorem 3.6" (p. 3150). They stress that "since u is convex, the projection u h should converge to u as h → 0, but this is not the case in this example" (p. 3150) and "(u h ) does not converge to u as h → 0" (p. 3151). They summarize this result by saying "although there is some sort of super-convergence for some meshes, for general meshes [...] FE-convex piecewise linear [m=2!] function may not suffice". Discussing how this numerical example articulates with [START_REF] Aguilera | Approximating optimization problems over convex functions[END_REF] is postponed to a forthcoming article.

All this is even complexified when the authors cope with numerical experiments of P 2 -P 2 FE discretization. As they notice, the FE test functions in the basis are not all nonnegative. So they "considered the usual piecewise linear nodal basis for the vertices" and "the usual quadratic bubbles" for the edges' midpoints (p. 3152).

Their table of convergence (Table 5.1) is transformed into loglog graphic in Figure 7. The last simulation with adaptive refinement gives a good point, but the uniformrefinement simulations give non-decreasing errors. Again convergence is not obvious. So as to make this convergence obvious, we illustrate their method in Section 5.4. [START_REF] Aguilera | Approximating optimization problems over convex functions[END_REF][START_REF] Aguilera | On convex functions and the finite element method[END_REF] contradictory with [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] ?

Are

The answer is no. But it deserves to be explained.

Both [START_REF] Aguilera | On convex functions and the finite element method[END_REF] and [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] use a dual definition of convexity and have conclusions that could seem to be contradictory. The only difference is that [START_REF] Aguilera | On convex functions and the finite element method[END_REF] uses different basis for trial and test functions:

-∇u h ⊗ ∇ϕ h s 0, while [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] uses the same basis for trial and test functions, both in P 1 :

-∇u h ⊗ ∇φ h 0.

Both discretizations are weak in a sense, but the one of [START_REF] Aguilera | On convex functions and the finite element method[END_REF] has less test functions than trial functions and so less contraints than the number of Degrees Of Freedom (DOF). Indeed m > 2 means that the trial basis is rather large (at least P 2 ) and condition C.3 that the test basis is at most P 1 . The matrix of constraints is rectangle and there are less constraints than the unknowns whereas the discretization of [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] uses the same basis for trial and test and so involves square matrices. Then, it is not surprising that the problem discretized with the weak P 1 -P 1 definition in [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] is overconstrained : there are as many constraints as DOF ! Indeed it appears also from the literature that the amount of constraints must be not too large. The article [START_REF] Mirebeau | Adaptive, Anisotropic and Hierarchical cones of Discrete Convex functions[END_REF] defines a stencil of constraints that is coarser than the grid and [START_REF] Mérigot | Handling convexity-like constraints in variational problems[END_REF] defines a subgrid of points on which the constraint is applied. Then it seems coherent that the FE proof of [START_REF] Aguilera | On convex functions and the finite element method[END_REF] requires more DOF than the number of constraints. The methods of [START_REF] Wachsmuth | Conforming Approximation of Convex Functions with the Finite Element Method[END_REF] and similar works look for functions in the interior of the convex functions set where the boundary are irrelevant.

Our Proposition 4.2 proved in subsection 3.2 states that if we discretize u h with P 2 FE and use a P 2 set of test functions (but then some test functions change sign), then the discretization of the linear part of the constraint is non-consistent. It proves that the proof of convergence of FEM by [START_REF] Aguilera | On convex functions and the finite element method[END_REF], which fails if the test functions change sign (as in P 2 and higher), maynot be improved on that point. So we do believe their theorem is optimal.

If we use the consistency test for some discretizations, we can compute the first nonzero term in the expansion of -∇u h ⊗ ∇ϕ h s , for P 1 -P 1 discretizations and report these computations in Table 1 (lower part) on various meshes depicted in Table 1 (upper part). It appears that some discretizations are not consistent with convexity (two on the left) while others are (two on the right).

Does the P 2 -P 1 FEM converge ?

The convergence is proved in [START_REF] Aguilera | On convex functions and the finite element method[END_REF] for P k (k ≥ 2) trial functions and P 1 test functions with additionnal assumptions. In order to ensure the numerical convergence, we wrote a Python code to mesh a square Ω = [1, 2] 2 , discretize the functional (2), and discretize the nonlinear constraints of inequality (convexity). Then we use the minimize function in scipy to minimize the functional. The execution time are meaningless since they rely on the level of precision we ask and we only want to illustrate convergence.

In a very first step, we look for the solution to (2,[START_REF] Goldstine | A history of the calculus of variations from the 17th to the 19th century[END_REF][START_REF] Carlier | Approximation of maximal Cheeger sets by projection, M2AN Math[END_REF], which exact solution is given in [START_REF] Rochet | Ironing, sweeping and multidimensionnal screening[END_REF]. It is not surprising that, whatever the mesh, we can get a relative error of 10 -12 . It only justifies that a quadratic function may be approximated by a P 2 numerical code without error.

In order to justify convergence, one needs a more complex convex function to find.

Convergence to an exact solution

We choose u exact (x) = exp((x 1 + x 2 )/2) -e, and take the functional

Ω | ∇u | 2 +∆u exact u dx - ∂Ω ∂u exact ∂n u,
which is also, up to a constant, Ω | ∇(u -u exact ) | 2 . The function is constrained to be convex and be such that u(1, 1) = 0, which makes sense since the function is convex and so regular enough for evaluating at (1, 1). The solution is convex, but we are going to prove that its P 2 interpolate is not convex. The geometrical argument is that we P 2 -interpolate a function of x 1 + x 2 on a mesh like mesh 1 in Figure 1. Let us take the P 2 interpolate of u exact on the two triangles depicted in Figure 8, close to the point (x 1 , y 1 ) at the center as before. Let us prove now that the P 2 interpolate u P 2 h = Π P 2 u ∈ P 2 is strictly convex inside T 2 . Along the segment between the point numbered 1 in T 2 and the point 5 in T 2 , one has:

u P 2 h (x 1 + αh, y 1 -αh) = 2u T 2 1 (-2α + 1)(-2α + 1/2) + 2u T 2 2 α(-α + 1/2) +2u T 2
3 α(α -1/2) + 4u T 2 4 α(-2α + 1) +4u T 2 5 α 2 + 4u T 2 6 α(-2α + 1) = e (x 1 +y 1 )/2 α 2 (12 + 2e -h/2 + 2e h/2 -8e -h/4 -8e h/4 ) +α(-6 -e -h/2 -e h/2 + 4e -h/4 + 4e h/4 ) + 1 , where u T 2 i is the value of the interpolate at the DOF locally numbered i in triangle T 2 . These values are the same as those of u exact by definition of the P 2 interpolate and enable to write the last equality. It is then easy to check that this function is strictly convex for α ∈ (0, 1/2). The same computations in the triangle T 1 gives for α ∈ (0, 1/2) u P 2 h (x 1 -αh, y 1 + αh) = e (x 1 +y 1 )/2 α 2 (12 + 2e -h/2 + 2e h/2 -8e -h/4 -8e h/4 ) +α(-6 -e -h/2 -e h/2 + 4e -h/4 + 4e h/4 ) + 1 .

Since both u P 2 ,T 2 h and u P 2 ,T 1 h are strictly convex in their domains and have the same value on three aligned points, the function u P 2 h maynot be convex at least on the segment between the points numbered 5 in T 2 and the point numbered 5 in T 1 (see Figure 8).

Despite the non-convexity of the P 2 interpolate of the exact solution (which is convex), the convergence of the FEM method proposed and proved in [START_REF] Aguilera | On convex functions and the finite element method[END_REF] is clearly illustrated in Figure 9 where structured and unstructured meshes are used.

The order of convergence seems to be the same on structured and unstructured meshes. It is roughly 1.6.

We also test the gradients' jumps of P 2 functions (strong discretization with C ∞ 0 test functions in subsection 4.1) and various weak P 2 discretizations for u h and for some P 2 test functions (subsections 4.2 and 4.3). Some of them are consistent and some are not.

We also discuss the literature, but let a deeper study of how [START_REF] Aguilera | Approximating optimization problems over convex functions[END_REF][START_REF] Aguilera | On convex functions and the finite element method[END_REF] interact with [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] to a forthcoming article. In a last subsection, we numerically illustrate the convergence of [START_REF] Aguilera | On convex functions and the finite element method[END_REF]'s FEM where the trial and test basis are different. Even if other methods exist, such as the recent methods using an evolution PDE converging to the convexified function [START_REF] Oberman | Computing the convex envelope using a nonlinear partial differential equation[END_REF][START_REF] Carlier | Exponential convergence for a convexifying equation[END_REF], this FEM is maybe one of the simplest. Comparisons of all these methods would be a thrilling challenge. Relative error

Theorem 1

 1 Let Ω be an open and convex subset of R 2 and T h a triangulation of Ω. If there exists an open and convex Ω ⊂ Ω such that the following property is satisfied in Ω and for the mesh T h : (P M ) ∃(a, b) two independent unit vectors such that (n • a)(n • b) ≥ 0 for all n unit normal to an edge of T h Ω ,

Figure 1 :Figure 2 :

 12 Figure 1: Three meshes.

1 Figure 3 :

 13 Figure 3: Local shape of the mesh close to (x 1 , y 1 ).

Figure 4 :

 4 Figure 4: Convergence of P 1 FE in case of a subharmonicity constraint.

Figure 5 :

 5 Figure 5: Local shape of the mesh close to (x 1 , y 1 ). Local numbering in triangle 3.

23

 23

Figure 6 :

 6 Figure 6: The monopolist numerical solution error versus DOF.

  ref. L 2 uniform ref. L 2 adaptive ref. L ∞ uniform ref. L ∞

Figure 7 :

 7 Figure 7: Error according to the refinement either adaptive (solid lines) or uniform (dashed lines) and the norm (L 2 or L ∞ ) versus Degrees of Freedom (DOF).

  h

  Figure 8: Local numbering.

Figure 9 :

 9 Figure 9: Convergence of the FEM P 2 -P 1 method.

Figure 10 :

 10 Figure 10: The monopolist's graph and its gradient (P 2 -P 2 ) on a 25 × 25 regular mesh.

Table 1 :

 1 2 u xx 2u xy 2u xy u yy h 2 u xx 0 0 u yy 4h 2 /3 u xx u xy u xy u yy 2h 2 /3 u xx u xy u xy u yy First order of the Hessian for various meshes (P 1 -P 1 ).

  Figure 11: The eigenvalues.

						0.010	1.8
						0.008
						0.004 0.006	1.6
						0.002
						0.000	1.4
	1.0 1.2 1.4 1.6 1.8 2.0 1.0	1.2	1.4	1.6	1.8	0.002 2.0	1.2	1.2	1.4	1.6	1.8

Acknowledgments

The author thanks G. Carlier for drawing his attention to the problem of subharmonicity approximation.

Can this method be used ?

In order to use this method justified theoretically (in [START_REF] Aguilera | On convex functions and the finite element method[END_REF]) and numerically (above), we look for a solution to the problem of the monopolist in Ω = [1, 2] 2 :

There exists no explicit solution to this problem but some properties of the solution are known (see [START_REF] Rochet | Ironing, sweeping and multidimensionnal screening[END_REF]).

The computed P 2 function on a structured mesh 25 × 25, its gradient and its second eigenvalue look similar to the ones computed in [START_REF] Mirebeau | Adaptive, Anisotropic and Hierarchical cones of Discrete Convex functions[END_REF]. The graph of the full P 2 function would be difficult to understand because of the number of DOF. That is why we only use the P 1 information from u h in Figure 10 (left).

So as to display the gradient, one might compute the gradient Ω ∇u h φ i for i associated to an interior point, with a P 2 function u h , and φ i ∈ P 2 . But since the quadrature uses the three mid-points, it is right only up to P 2 functions, and so the integral being over a P 3 function is not exact. Yet, since the P 2 -P 1 discretization, which is exact, gives the same graphical result with less points, we consider the Figure with P 2 test functions to be acceptable. The result may be seen in Figure 10 (right).

In Figure 11, one may see the second eigenvalue of the Hessian, computed as a P 2 function tested on P 1 functions and its level sets on the right. The first eigenvalue is not significant because of numerical artifacts. They are computed as the two roots per interior vertex of:

-Ω ∂u h /∂x ∂φ i /∂x ∂u h /∂x ∂φ i /∂y Ω ∂u h /∂y ∂φ i /∂x ∂u h /∂y ∂φ i /∂y .

CONCLUSION

In this article we prove (Theorem 3) that the P 1 discretization of a function that satisfies a strong definition of convexity (with C ∞ 0 test functions), which is equivalent to the gradients' jumps positivity (for P 1 FE), and so is conformal, leads to an additional constraint on the limit function. While [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] dealt mainly with regular meshes, our result extends explicitely to any mesh. The error of such a discretization of the constraint does not vanish with the space step. We justify that it is localized where the additional constraint on the limit function is not satisfied. The condition for existence of such a counter-example requires information both from the mesh and its refinement but is general.

In addition, among the P 1 discretization of u with a strong (section 3.1) or weak (section 3.2) definition of convexity, not all are consistent. The definition of consistency is very similar to the one of partial differential equations and it is used above to discriminate likely discretizations and unlikely ones. But this does not guarantee good numerical results even when the discretization is consistent. We also discuss some structured meshes, the P 1 discretization of u and the use of P 1 test functions for convexity (say P 1 -P 1 convexity). we prove they may or maynot be consistent (Table 1).