{'fraction_non_alphanumeric': 0.03480636954678225, 'fraction_numerical': 0.013059455384905305, 'mean_word_length': 4.454975328947368, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Parametric Shape Optimization using the Support Function

The optimization of shape functionals under convexity, diameter or constant width constraints shows numerical challenges. The support function can be used in order to approximate solutions to such problems by finite dimensional optimization problems under various constraints. We propose a numerical framework in dimensions two and three and we present applications from the field of convex geometry. We consider the optimization of functionals depending on the volume, perimeter and Dirichlet Laplace eigenvalues under the aforementioned constraints. In particular we confirm numerically Meissner's conjecture, regarding three dimensional bodies of constant width with minimal volume.


Introduction
Shape optimization problems are a particular class of optimization problems where the variable is a shape. A typical example of such a problem has the Finite dimensional approximations of convex bodies can be built using a truncation of a spectral decomposition of the support function: Fourier series decomposition for d = 2 and spherical harmonic decomposition for d = 3. The shape optimization problem becomes a parametric optimization problem in terms of the coefficients of the spectral decomposition. The convexity constraint is modeled by a set of linear pointwise inequalities for d = 2 or quadratic pointwise inequalities for d = 3. The constant width constraint is characterized by the fact that coefficients of the even basis functions are zero. Diameter constraints can also be translated into a set of pointwise linear inequalities. It can be noted that in some particular cases, functionals like volume and perimeter have explicit formulas in terms of the coefficients in the above decompositions.
In [9] the authors study numerically optimization problems under constant width constraint in dimension two, with the aid of the support function and Fourier series decomposition. They work with a global parametrization of the convexity constraint, which requires the use of specific semidefiniteprogramming techniques and software. We choose to work in a simplified framework, inspired from [10], in which the convexity constraint is imposed on a finite, sufficiently large, number of points distributed on the unit circle (d = 2) or on the unit sphere (d = 3). This gives rise to simpler constrained optimization problems that can be handled using standard optimization software.
In Section 3 we provide existence results for the problems we consider. Moreover, we prove that solutions obtained when using a truncated spectral decomposition of the support function converge to the solutions of the original problems as the number of coefficients goes to +∞.
The paper is organized as follows. In Section 2 various properties of the support function parametrization in dimension two and three are recalled. In Section 3 theoretical aspects regarding the existence of solutions and the convergence of the discrete solutions are investigated. Section 4 deals with the parametric representation of shapes using the spectral decomposition of the support function. The handling of convexity, constant width, diameter and inclusion constraints is discussed. Section 5 recalls the method of fundamental solutions used for solving the Dirichlet-Laplace eigenvalue problems.
Section 6 contains applications of the numerical framework proposed for various problems in convex geometry. In particular, a confirmation of the Meissner conjecture regarding bodies of constant width with minimal volume in dimension three is provided. The two different Meissner bodies are obtained by directly minimizing the volume under constant width constraint, starting from general random initializations. Further applications presented in Section 6 concern the minimization of eigenvalues of the Dirichlet-Laplace operator under convexity and constant width constraints, approximation of rotors of minimal volume in dimension three, approximation of Cheeger sets and the minimization of the area under minimal width constraint.
The goal of this paper is to present a general method for performing shape optimization under various non-standard constraints: convexity, fixed width, diameter bounds by transforming them into algebraic constraints in terms of a spectral decomposition of the support function. In order to illustrate the method, various numerical results are presented, some of which are new and are listed below: • optimization of the Dirichlet-Laplace eigenvalues under convexity constraint: the case k = 2 in dimension two was extensively studied (see for example [11], [12]). The case k ≥ 3 in dimension two and the simulations in dimension three are new and are presented in Section 6.1. • minimization of the Dirichlet-Laplace eigenvalues under fixed constant width constraint in dimension three: in Section 6.2 cases where the ball is not optimal are presented. • numerical confirmation of Meissner's conjecture: in Section 6.3 the proposed numerical framework allows to obtain Meissner's bodies when minimizing the volume under fixed width and convexity constraints, starting from randomized initial spherical harmonics coefficients. • rotors of minimal volume in dimension three: the numerical framework allows to approximate optimal rotors in the regular tetrahedron and the regular octahedron. (Section 6.4) Other results deal with the computation of Cheeger sets for various two dimensional and three dimensional domains and the minimization of area of a three dimensional body under minimal width constraint.

Support function parametrization
This section recalls some of the main properties of the support function, as well as the properties which will be used in order to implement numerically the various constraints of interest in this work. The references [9], [13], [14] and [15] contain more details about this subject. Let B be a convex subset of R d . The support function of B is defined on the unit sphere S d−1 by where the dot represents the usual Euclidean dot product. Geometrically, p(θ) represents the distance from the origin to the tangent plane α to B such that α is orthogonal to θ, taking into account the orientation given by θ. Therefore the sum of the values of the support function for two antipodal points gives the width or diameter of B in the direction defined by these two points. This shows that bounds on the width of B can be expressed by inequalities of the type w ≤ p(θ) + p(−θ) ≤ W for every θ ∈ S d−1 and a constant width constraint can be expressed by w = p(θ) + p(−θ) for every θ ∈ S d−1 .
As already shown in [10], it is possible to impose inclusion constraints when dealing with support functions. Consider two convex bodies B 1 , B 2 with support functions given by p 1 , p 2 . Then B 1 is included in B 2 if and only if p 1 (θ) ≤ p 2 (θ) for every θ ∈ S d−1 . In the case where B 2 is an intersection of half-spaces the inequality p 1 (θ) ≤ p 2 (θ) only needs to be imposed for a finite number of directions θ ∈ S d−1 , corresponding to the normals to the hyperplanes determining the hyperspaces.
Each convex body in R d has its own support function, however the converse is not true. The necessary assumptions for a function p : S d−1 → R to be the support function of a convex body together with precise regularity properties of the support function are discussed in [13,Section 1.7]. In the following, we make the assumption that the support functions are smooth enough so that first and second partial derivatives can be computed. In Section 3 it is recalled that this assumption is not too restrictive, since the class of convex shapes with smooth support function is dense in the class of convex sets with respect to the Hausdorff metric. Moreover, support functions used in the numerical simulation verify this assumption.
Given a convex set B and its support function p, a parametrization of ∂B is given by where ∇ τ represents the tangential gradient with respect to the metric in S d−1 . Note that for this parametrization the normal to ∂B at the point x(θ) ∈ ∂B is given by θ. The convexity constraint can be expressed by the fact that the principal curvatures of the surface are everywhere non-negative. In the following, the aspects related to convexity are detailed and the presentation is divided with respect to the dimension.

Dimension 2
In R 2 , the unit circle S 1 is identified to the interval [0, 2π], therefore the parametrization of the boundary of the shape in terms of the support function becomes It is immediate to see that (x (θ), y (θ)) = p(θ) + p (θ) and, as already underlined in [9], the convexity constraint in terms of the support function is equivalent to p + p ≥ 0, in the distributional sense.

Dimension 3
In R 3 it is classical to consider the parametrization of S 2 given by n = n(φ, ψ) → (sin φ sin ψ, cos φ sin ψ, cos ψ), φ ∈ [−π, π), ψ ∈ (0, π). (2) As recalled in [15], if p = p(φ, ψ) is a C 1 support function then a parametrization of the boundary is given by Moreover, the differential dx p on the basis n φ , n ψ of the corresponding tangent space to S 2 is given by Note that {n φ / sin(ψ), n ψ } is an orthonormal basis of the tangent space when sin ψ = 0. The convexity constraint is characterized by the non-negativity of the principal curvatures. Another characterization is given by the positive definiteness of the matrix having coefficients given by the differential (4) of x p for every φ ∈ [−π, π) and ψ ∈ [0, π). However, in dimension 3 it is enough to impose a simpler condition. Indeed, if a surface has non-negative Gaussian curvature in a neighborhood of a point, then the surface is locally convex around that point. Tietze's theorem states that if a three dimensional set is locally convex around each point then it is globally convex [16, p. 51-53]. As a consequence, a closed surface in dimension 3 which has positive Gaussian curvature everywhere bounds a convex body. This is also known as Hadamard's problem [17, p. 108]. Therefore, in dimension three, the convexity constraint can be imposed by assuring that the Gaussian curvature is positive at every point. More explicitly, the determinant of the matrix containing the coefficients of the differential (4) is positive: for every φ ∈ [−π, π), ∀ψ ∈ (0, π). Note that formulas (3), (4) and (5) contain sin ψ in some of the denominators. In the numerical simulations the discretization of the sphere is always chosen avoiding the north and south poles of the sphere S 2 where sin ψ cancels and where singular behavior may occur.

Theoretical aspects 3.1 Existence of optimal shapes
When dealing with shape optimization problems the existence of optimal shapes is a fundamental question. All problems studied numerically in this article deal with convex domains that are contained in a bounded set, so it is useful to define the class K d of closed convex sets in R d which are contained in a closed, large enough, ball B. The question of existence of solutions is greatly simplified due to the following classical result [13,Theorem 1.8.7].
Theorem 1 (Blaschke selection theorem) Given a sequence {Kn} of closed convex sets contained in a bounded set, there exists a subsequence which converges to a closed convex set K in the Hausdorff metric.
For the sake of completeness, recall that the Hausdorff distance between two convex bodies K 1 , K 2 is defined by A sequence of closed convex sets {K n } converges to K in the Hausdorff metric if and only if d(K, K n ) → 0 as n → ∞.
More details regarding this result and proofs can be found in [18,Chapter 2], [19,Chapters 6,7]. Existence results for all problems studied in this paper are a consequence of the properties listed below. These properties are classical, but are recalled below from the sake of completeness, with sketches of proof when the proofs were not readily found in the literature. If K 1 , K 2 have support functions p K1 and p K2 then the Hausdorff distance is simply Property 1 Convexity is preserved by the Hausdorff convergence.
Property 2 If {Kn} is a sequence of non-empty closed convex sets contained in a bounded set then the Hausdorff convergence of Kn to K is equivalent to the uniform convergence of the support functions p Kn to p K on S d−1 .
For a proof see [20,Theorem 6]. In the following, χ K denotes the characteristic function of the set K.
Property 3 Suppose that the sequence of convex sets {Kn} converges to the convex set K in the Hausdorff topology and that K has non-void interior. Then χ Kn converges to χ K in L 1 . As a consequence, |Kn| → |K| and H d−1 (∂Kn) → H d−1 (∂K) as n → ∞.
A proof of this fact can be found in [21,Proposition 2.4.3].
The Dirichlet-Laplace eigenvalues are solutions of It is classical that for Lipschitz domains, the spectrum of the Dirichlet-Laplace operator consists of a sequence of eigenvalues (counted with multiplicity) In particular, convex sets with non-void interior enter into this framework. Two basic properties of the Dirichlet-Laplace eigenvalues are the monotonicity with respect to inclusion and the scaling property: The following property deals with the continuity of these eigenvalues with respect to the Hausdorff metric.   For a proof see [18, p. 33].

Property 6
The diameter and width constraints are continuous with respect to the Hausdorff convergence of closed convex sets. In particular if the sequence of closed convex sets {Kn} converges to K in the Hausdorff metric and each Kn is of constant width w then K is also of constant width w.
Proof: Property 2 recalled above shows that the Hausdorff convergence implies the uniform convergence of support functions on S d−1 . Therefore, diameter and width constraints that can be expressed in pointwise form starting from the support functions are preserved, in particular, the constant width property.
In the following, {φ i } ∞ i=0 denotes an orthogonal basis of L 2 (S d−1 ) made of eigenvalues of the Laplace-Beltrami operator on S d−1 (the Fourier basis in 2D and the spherical harmonics in 3D). Denote by λ i ≥ 0, i ≥ 0, the corresponding eigenvalues: −∆ τ φ i = λ i φ i on S d−1 . In particular λ 0 corresponds to the constant eigenfunction φ 0 . When studying rotors, only some particular coefficients in the spectral decomposition are non-zero. Also, in numerical approximations a truncation of the spectral decomposition is used. Therefore, it is relevant to see if such a property is preserved by the Hausdorff convergence of convex sets. Let J ⊂ N be a non-empty, possibly infinite subset of indices and denote by i.e. convex shapes for which the coefficients of the support function in the basis {φ i } ∞ i=0 having indices that are not in J are zero. The following result holds: Property 7 For a fixed set J ⊂ N, let {Kn} be a sequence of closed convex sets contained in a bounded set B with support functions (p Kn ) contained in F J such that Kn converges to K in the Hausdorff metric. Then the support function of K also belongs to F J .
Proof: From Property 2 it follows that the support functions of K n converge uniformly to the support function of K, i.e. p Kn − p K ∞ → 0 as n → ∞. This obviously implies the convergence in L 2 (S d−1 ) of the support functions.
Since the left hand side converges to zero as n → ∞ for all i, it follows that α n i → α i for all i. In particular, if i / ∈ J then all α n i = 0. As a consequence α i = 0 for all i / ∈ J, which means that p K ∈ F J . It is possible to characterize constant width bodies by imposing that certain coefficients of the support function are zero. Proof: The support function of a body of constant width satisfies p(θ) + p(−θ) = w for every θ ∈ S 2 . The conclusion follows from the fact that φ i form an orthogonal basis of L 2 (S d−1 ) and S d−1 φ i = 0 for all i ≥ 1. Starting from these properties, the existence of solutions for all problems studied numerically in the following section can be proved. Recall that λ k (Ω) denotes the k-th eigenvalue of the Dirichlet-Laplace operator defined by (6).
The existence of solutions for this problem is proved in [22,Theorem 2.4.1]. It is a direct consequence of Theorem 1 and Properties 3, 4 above. For proving the existence in this case, choose a minimizing sequence {K n } which, by Theorem 1, up to a subsequence, converges to K. In view of the Properties 3, 4 and 6 above K is indeed a solution. The existence follows from Properties 3, 6 when working with a converging minimizing sequence, which exists by Theorem 1.
The next problem considered concerns rotors of minimal volume. A rotor is a convex shape that can be rotated inside a polygon (or polyhedron) while always touching every side (or face). A survey on rotors in dimension two and three can be found in [23]. In particular, the article [23] describes which coefficients are non-zero in the spectral decomposition of the support function of rotors, using Fourier series in 2D or spherical harmonics in 3D. It turns out that the earliest complete development on the subject was published in 1909 by Meissner [24]. More details and proofs of the claims in the papers described above can be found in [25].
In dimension two, every regular n-gon admits non-circular rotors and they are characterized by the fact that only the coefficients for which the index has the form nq ± 1 are non-zero, where q is a positive integer. In dimension three, there are only three regular polyhedra which admit rotors: the regular tetrahedron, the cube and the regular octahedron. The rotors in a cube are bodies of constant width. For rotors in a tetrahedron the only non-zero coefficients correspond to the spherical harmonics with indices 0, 1, 2 and 5, while in the case of the octahedron the non-zero coefficients have indices 0, 1 and 5. The constant term in the spectral decomposition of the support function of a rotor corresponds to the inradius of the domain.
Assuming P admits rotors, the existence of rotors of minimal volume is guaranteed by Theorem 1 and properties 5, 7. Moreover, the fact that the Hausdorff limit of rotors is still a rotor comes from Property 2. It is enough to choose the normal directions orthogonal to the sides of the polygon (polyhedron) and observe the limit of the corresponding support functions evaluated at these directions.
The following problem was considered in [7] and consists in minimizing the area under minimal width constraint.
Problem 5 Minimize the area of a convex three dimensional shape ω which has minimal width equal to 1.

In [21, Proposition 2.4.3], it is proved that of two convex bodies
. Every body with minimal width equal to 1 contains three mutually orthogonal segments of length ≥ 1. By convexity and the property above it is immediate to see that the volume of any body in K is at least 1 8 . Using the isoperimetric inequality, it can be seen that minimizing sequences exist and the existence of a solution comes from Theorem 1 and Property 6.
As an application for the inclusion constraint, the Cheeger set associated to some convex domains in dimension two and three is considered.
Problem 6 (Cheeger sets.) The Cheeger set associated to a convex domain Ω ⊂ R d is the solution of the problem where the minimum is taken over all convex sets X contained in Ω.
The Cheeger sets are extensively studied and it is not the objective to present the subject in detail here. In dimension two there is an efficient characterization which allows the analytical computation of Cheeger sets for a large class of domains [26]. Computational approaches based on various methods were introduced in [3], [27], [28] and [29]. Note that in dimension two, the convexity of Ω implies the convexity of the optimal Cheeger set. In dimension three this is no longer the case. However, one can prove that there exists at least one convex optimal Cheeger set [3]. Existence of Cheeger sets is a consequence of Theorem 1 and Properties 3, 5.

Convergence results
The numerical approach described in Section 4 uses a truncation of the spectral decomposition of the support function. Therefore, as underlined in [9], it is needed to prove that increasing the number of non-zero coefficients N in the parametrization gives optimal shapes which converge to the solution of the original problem.
In the following, denote by This corresponds to the notation F J used previously, with J = {0, 1, 2, ..., N }.
Denote by K d N the class of convex sets in R d whose support functions belong to F N . In [13,Appendix] it is proved that for N large enough, F N is not trivial and therefore K d N is also non-trivial. Property 7 proved in the previous section shows that K d N is closed in the Hausdorff metric. Therefore the existence of solutions can be shown for all problems recalled in the previous section by replacing the class of convex sets with F N .
An important question, which is not obvious at first sight, is whether a general convex body K can be approximated in the Hausdorff metric by convex bodies K n with support functions in some F Nn . This is proved in [13,Section 3.4]. Namely, the following property holds: Property 9 Let K be a convex body and ε > 0. Then there exists a positive integer Nε > 0 and a convex set Kε with support function in F Nε such that d H (K, Kε) < ε.
Remark 1 Following the remarks in [13, p. 185], starting from a body of constant width K, the smoothing procedure preserves the constant width. Moreover, the approximation of K in K d Nε is obtained by truncating the spectral decomposition of the regularized support function. This also preserves the constant width, which by Property 8 is equivalent to the fact that the coefficients of the even basis functions are zero. Therefore, if K is of constant width in the previous proposition, its approximation Kε ∈ F Nε can also be chosen of the same constant width.
In practice, however, it is often necessary to impose some other constraints, like fixed volume, area, minimal width, etc. Below we give another variant of this property for constraints of the form {C(K) ≥ c} where C is a continuous function with respect to the Hausdorff metric which is homogeneous of degree α > 0: C(ηK) = η α C(K) for η > 0. This includes many constraints of interest, like area, perimeter, minimal width, diameter, etc.
Property 10 Let C : K n → R + be a continuous functional, positively homogeneous of degree α > 0. Let K be a convex body which satisfies C(K) ≥ c, for some fixed c > 0 and ε > 0. Then there exists a positive integer Nε > 0 and a convex set Kε with support function in F Nε such that d H (K, Kε) < ε and C(Kε) ≥ c.
Proof: Property 9 implies the existence of a sequence K n converging to K in the Hausdorff metric such that K n ∈ K d Nn for some N n > 0. Since C is continuous it follows that C(K n ) → C(K). Define the new sets K n by The relation between the Hausdorff distance and the support functions implies that Therefore it is possible to approximate K in the Hausdorff metric with the sequence K n which also verifies the constraint C(K n ) ≥ c. Now we are ready to prove the following approximation result.
Theorem 2 Let G be a continuous functional defined on the class of closed convex sets. Consider a constraint function C which is continuous for the Hausdorff metric and let c > 0. In the following M denotes one of the following: Then any converging subsequence of K N converges in the Hausdorff metric to a solution K of min Proof: First, let us note that the existence of solutions to problems (7) and (8) is immediate using Theorem 1, Property 7 the fact that M is closed and the continuity of G.
Denote by K a solution of (8). By the results shown in Properties 9, 10 and Remark 1 there exists a sequence of convex sets L Nn ∈ K d Nn ∩ M such that L Nn → K in the Hausdorff metric.
In the following denote by K N a solution of (7) for N ≥ 1. It is obvious Therefore, the sequence {G(K N )} N ≥1 is nonincreasing and bounded from below, which implies the existence of the limit lim N →∞ G(K N ) = ≥ G(K). By the optimality of K N we have G(L Nn ) ≥ G(K Nn ) and since d H (K, L Kn ) → 0, by the continuity of G we obtain The inequality above implies that = G(K). Furthermore, if a subsequence of {K N } N ≥1 converges to K in the Hausdorff metric it follows that K ∈ M and G(K ) = = min K∈K d ∩M G(K). Therefore, every limit point for {K N } N ≥1 in the Hausdorff metric is a minimizer of (8).
Theorem 2 motivates our numerical approach. In order to obtain an approximation of solutions of the shape optimization problems considered, a truncation of the spectral decomposition is used. The theoretical result states that the solutions of the finite dimensional minimization problems obtained converge to the solution of the original problem.

Numerical framework
In numerical shape optimization, shapes are represented using a finite number of parameters. Theorem 2 shows that considering as parameters the coefficients of the truncation of a spectral decomposition is appropriate, since solutions of the resulting finite dimensional optimization problems converge to solutions of the convex shape optimization problem. This type of numerical approaches was already used in other contexts in [9], [10], [30], [31], [32]. Using such systems of orthogonal basis representations has further advantages which will be underlined below. Again, for the clarity of exposition, we divide the presentation following the dimension.

Dimension 2
We approximate the support function by a truncated Fourier series As stated in Section 2, in order for p to be the support function of a convex set in R 2 we need to have p (θ) + p(θ) ≥ 0 for every θ ∈ [0, 2π). In [9] the authors provide an exact characterization of this condition in terms of the Fourier coefficients, involving concepts from semidefinite programming. In [10] the author provides a discrete alternative of the convexity inequality which has the advantage of being linear in terms of the Fourier coefficients. We choose θ m = 2πm/M c , m = 1, 2, ..., M c for some positive integer M c and we impose the inequalities p(θ m ) + p (θ m ) ≥ 0 for m = 1, ..., M c . As already shown in [10] we obtain the following system of linear inequalities where α m,n = (1 − n 2 ) cos(nθ m ) and β m,n = (1 − n 2 ) sin(nθ m ).
An upper bound W on diameter can be introduced as a constraint for the support function as follows In the computations we consider a discrete version of the above inequality. Pick θ m = 2πm/M d , m = 1, 2, ..., M d for some positive integer M d and impose the following linear inequalities In order to impose a lower bound on the diameter it is enough to pick one direction θ and use the constraint It is also possible to consider variable lower and upper bounds on the width of the body which depend on θ.
Let us now recall the formulas for the area and perimeter of a two dimensional shape in terms of the Fourier coefficients of the support function. The perimeter is simply equal to P (p) = 2πa 0 , which is linear in terms of the Fourier coefficients. As already stated in [9] the area of a convex shape having support function p with the Fourier decomposition (9) is given by Note that a 1 and b 1 do not contribute to the area computations as modifying a 1 , b 1 only leads to translations of the shape defined by p.

Dimension 3
In [32] the authors parametrized three dimensional domains by their radial function using spherical harmonics. In our case we parametrize the support function using a finite number of spherical harmonics for a given positive integer N . The spherical harmonics are defined by where P m l are the associated Legendre polynomials and are normalization constants. The convexity constraint is imposed by considering a discrete version of (5). We construct a family of M c evenly distributed points on the unit sphere, for example like described in [33,Section 3]. We denote by (φ i , ψ i ) i = 1, ..., M c the corresponding pairs of angles given by the parametrization (2). We impose that the convexity condition (5) is satisfied at points given by (φ i , ψ i ), i = 1, ..., M c . As in the two dimensional case, width inequality constraints can be handled in a similar way, by imposing inequalities of the type (2)).
The constant width condition is p(θ) + p(−θ) = w for every θ ∈ S 2 . As shown in Property 8 this amounts to considering only odd spherical harmonics in the decomposition (12), except for the constant term. This corresponds to spherical harmonics Y m l for which the index l is odd. In the following, denote by h the part of the support function containing the non-constant terms: h = p − 1 4π S 2 pdσ. The area and volume of a convex body of constant width w in dimension three can be computed explicitly in terms of the spherical harmonics coefficients. Indeed, in [14, Theorem 2], the following formulas are provided: where The formulas in [14] are for a body of constant width 2w, which we transform so that they correspond to a body of width w. The spherical harmonics Y m l form an orthonormal family, therefore it follows that E(p) can be computed explicitly in terms of the coefficients a l,m and the eigenvalues λ l,m corresponding to the spherical harmonics Y l,m : As a consequence, when dealing with bodies of constant width, the volume and the area have explicit formulas in terms of the coefficients a l,m of the decomposition (12). We note that it is also possible to compute explicitly the area of a general convex body, using the coefficients of the support function. Indeed, Lemma 1 from [14, Section 5] is valid for general support functions h, not only those corresponding to a constant width body. Therefore the area of a convex body B is also given by (14), where w = 2a 0,0 Y 0 0 . Also following the results stated in [14] it should also be possible to compute the volume explicitly using the Gaunt coefficients involving integrals on the sphere of products of three spherical harmonics. In our computations, for general bodies parametrized using their support function, we use the divergence theorem. The volume of a convex body ω is computed as the integral on ∂ω of a vector field V with divergence equal to one. For simplicity choose V = 1 3 x = 1 3 (x, y, z) and we integrate V · n on ∂ω. It is straightforward from the parametrization (3) that x(θ) · n(θ) = p(θ) for θ ∈ S 2 .

Visualization of results
We briefly present how the results are visualized. The variables in the optimization algorithm, and therefore, the output obtained are coefficients of a truncated spectral decomposition of the support function. Given such a family of Fourier coefficients (spherical harmonics coefficients) it is possible to evaluate the support function and its derivatives at any point in the unit circle (unit sphere in dimension three).
Once the values of the support function and its derivatives are known at a family of discretization points it is possible to use the formula (1) in dimension two (formula (3) in dimension three) in order to find the associated points on the boundary of the domain. When such a family of points is known a simple contour plot is made in dimension two. In dimension three, the Matlab command convhull is used to generate a triangulation of the surface bounding the convex body. Then the command patch is used to plot the surface of the convex body.

Computation of the Dirichlet-Laplace eigenvalues
The Dirichlet-Laplace eigenvalue problem (6) is solved numerically using the Method of Fundamental Solutions (MFS) [34,35]. This method does not need the construction of a mesh, is precise and has low computational time in dimensions two and three. We consider a fundamental solution of the Helmholtz equation, and respectively for 2D and 3D cases, where H where the source points y j are placed on an admissible source set, for instance the boundary of a bounded open setω ⊃ω, with ∂ω surrounding ∂ω. By construction, the MFS approximation satisfies the PDE of the eigenvalue problem (6) and we can just focus on the approximation of the boundary conditions, which can be justified by density results (see, for example [35]). Next, we give a brief description of the numerical procedure for calculating the Dirichlet-Laplace eigenvalues. We define two sets of points W = {w i , i = 1, ..., n} and X = {x i , i = 1, ..., m}, almost uniformly distributed on the boundary ∂ω, with n < m and the set of source points, Y = {w i + αn i , , i = 1, ..., n} where α is a positive parameter and n i is the unitary outward normal vector at the point w i . We consider also some interior points z i , i = 1, ..., p with (p < m) randomly chosen in ω and used the Betcke-Trefethen subspace angle [36]. After defining the matrices (20) and A(λ) = M 1 (λ) M 2 (λ) we compute the QR factorization and calculate the smallest singular value of the block Q 1 (λ), which will be denoted by σ 1 (λ). The approximations for the Dirichlet-Laplace eigenvalues are the local minima λ, for which σ 1 (λ) ≈ 0.

Shape derivatives for the Dirichlet Laplace eigenvalues
Functionals like volume or area have explicit formulas in terms of the coefficients in the Fourier or spherical harmonics decomposition. This gives straightforward formulas for gradients and Hessians which can be used in optimization algorithms. When the shape functional is more complex, direct formulas are not available. Below, we present how the Hadamard shape derivatives can be used to obtain partial derivatives in terms of coefficients of the parametrization.
The Hadamard shape derivative formula shows how a shape functional F(ω) varies when considering some perturbation of the boundary given by a vector field V . One way to define this to consider the derivative of the functional t → F((Id + tV )(ω)) at t = 0. Under mild regularity assumptions it can be proved that the shape derivative may be written as a linear functional depending on the normal component of V . For more details one could consult [18,Chapter 5] or [19,Chapter 9]. In particular, in [18,Theorem 5.7.4], for the case of the eigenvalues of the Dirichlet-Laplace operator (6) the shape derivative is as soon as eigenvalue λ k (ω) is simple and the eigenfunction u k is in H 2 (ω). This is true in the particular case of convex sets. In the following, we suppose that the functional F(ω) has a Hadamard shape derivative which can be written in the form

Dimension 2
As already recalled in Section 4 a parametrization of the boundary of the convex shape defined by the support function p is given by and a straightforward computation gives x (θ) = −(p (θ) + p(θ)) sin θ y (θ) = (p (θ) + p(θ)) cos θ.
Therefore the norm of the velocity vector is given by (x (θ), y (θ)) = p (θ) + p(θ), which helps compute the Jacobian when changing variables. Moreover, the normal to the point corresponding to parameter θ is simply (cos θ, sin θ). In the following, (a k ) k≥0 and (b k ) k≥1 denote the Fourier coefficients of the support function: p(θ) = a 0 + k≥1 (a k cos(kθ) + b k sin(kθ)). In order to compute the partial derivatives of the functional with respect to the Fourier coefficients it is enough to transform the perturbation of the support function into a perturbation of the boundary and use the Hadamard formula. We summarize the derivative formulas below. For simplicity, we use the abuse of notation f (θ) = f (x(θ), y(θ)).

1.
Derivative with respect to a 0 . The corresponding boundary perturbation is V = (cos θ, sin θ) and the normal component is V · n = 1. Therefore the derivative is

Dimension 3
We differentiate now a functional F(ω) for 3D shapes parametrized using the coefficients of the spherical harmonic decomposition (12) of the support function. Given a general perturbation of the support function p → p + Y , in view of (3), we find that the boundary perturbation has the form V = Y n + P(θ, φ)n φ + Q(θ, φ)n ψ .
Since the vectors n, n φ and n ψ are orthogonal, the normal component is given by V · n = V · n = Y . Then, using the general Hadamard derivative formula (21) we find the partial derivatives of the objective function with respect to the coefficients in (12): where Jac(ψ, φ) is the Jacobian function given by (5). Indeed, the Jacobian this surface integral is computed by Jac(φ, ψ) = ∂ φ x p × ∂ ψ x p . Note that the vectors ∂ φ x p , ∂ ψ x p are orthogonal to the normal n to the surface. Therefore the Jacobian reduces to Jac(φ, ψ) = n · (∂ φ x × ∂ ψ x) and using the expressions of the differential of x in the tangent plane given by (4) we can conclude that Jac(φ, ψ) is indeed given by (5).

Applications
This section shows how the proposed numerical framework applies to the various problems presented in Section 3. For each problem the shapes are discretized using a truncated spectral decomposition of the support function, as shown in Section 4. The corresponding shape optimization problem is then approximated by a finite dimensional constrained optimization problem. The Matlab fmincon function with the interior-point algorithm is used in each of the various problems shown below. As shown in previous sections, all constraints on the shapes are transformed into algebraic constraints on the coefficients of a spectral decomposition. The full functionality of fmincon is used in order to handle: linear equality or inequality constraints and non-linear constraints. The gradient of the functional and the gradients of the constraints are computed and are used in the algorithm. When possible, the Hessian matrix is also computed, and in all the other computations a LBFGS approximation is used.
The optimization toolbox described previously is efficient when the gradients of the objective function and of the constraints are provided. The computation time mostly depends on the size of the problem (number of variables and constraints) and on the cost of the evaluation of the objective function and its gradient. The cost of one objective function evaluation is as follows: • geometric functions (area, volume, surface area): explicit formulas, fast evaluation of the objective function. • Dirichlet Laplace eigenvalues using MFS: roughly 30 seconds in dimension two and one minute in dimension three. The cost of computing higher eigenvalues is larger than for lower ones.
The cost of the optimization algorithm on a personal computer (Intel i7 processor, 4.2Ghz, 32GB RAM) ranges from under 10 minutes for functionals involving only geometric quantities to a couple of hours for functionals involving the Dirichlet Laplace eigenvalues in dimension three.

Minimize the Dirichlet-Laplace eigenvalues under volume and convexity constraints
This section presents numerical approximations of solutions to Problem 1. The theoretical and numerical study of minimization problems of the form gained a lot of interest in the recent years. Various problems were considered, like the optimization of eigenvalues under volume constraint [37], [38], the optimization under perimeter constraint [39] and recently, the minimization under diameter constraint [40]. For many of the problems considered, explicit solutions are not known, therefore various works, like [12], [30], [32] deal with the optimization of the eigenvalues for volume and perimeter constraints. Such constraints can naturally be incorporated in the functional, in view of the behaviour of the eigenvalue with respect to scaling, and therefore unconstrained optimization algorithms based on information given by the shape derivative are successfully used in practice. The challenges encountered when adding the convexity constraint are underlined in the study of the second eigenvalue: This problem is studied in [41] where it is shown that the optimizer is not the convex-hull of two tangent disks, as conjectured before. Moreover, the boundary of the optimal set cannot contain arcs of circles. An algorithm for finding numerically the minimizer of (22) was proposed, using a penalization of the difference between the volume of the shape and the volume of its convex hull. A more precise, parametric search for the minimum of (22) was done in [11], giving an optimal numerical value of λ 2 (ω) = 37.987. The numerical algorithm proposed in this article allows us to study problem (22) in dimensions two and three. The computation of the eigenvalues is done using the method of fundamental solutions described in Section 5. The partial derivatives with respect to the Fourier coefficients in the parametrization are computed using results shown in Section 5.1 and the convexity constraint is imposed using (10). Two dimensional results are summarized in Figure 1 and it can be noted that segments are sometimes present in the boundaries of the numerical minimizers. For k = 3 we find that the minimizer is the disk, which is in accord with the simulations performed using only a volume constraint in [12], [30]. The values presented in Figure 1 are obtained by rounding up the numerical optimal values and are thus upper bounds for the optimal values. The numerical simulations are made for N = 300 (601 Fourier coefficients) and M c = 5000 points where convexity constraints are imposed.
The numerical discretization proposed in Section 4.2 makes it possible to perform the same simulations in dimension three, with no additional difficulty. These results are shown in Figure 2

Minimize the Dirichlet-Laplace eigenvalues under constant width constraint
This section presents numerical approximations of solutions to Problem 2 in dimension three. The two dimensional case was studied in [40], where it was noted that the disk appears more often as a local minimizer. In fact, the precise list of indices k for which the disk is a local minimizer for λ k (ω) under constant width constraint in dimension two is given in [40]. The convexity and constant width constraints are imposed as indicated in Section 4.2. We use gradient information in order to perform the optimization as indicated in Section 5.1. In our computations we observe that the ball appears often as a minimizer, but as observed in the two dimensional case in [40], we expect that this only happens for finitely many indices k. Notable exceptions are the indices corresponding to a simple eigenvalue for the ball. Figure 3 shows the non-trivial shapes of constant width obtained with our algorithm for k ∈ {10, 46, 99}, the three smallest indices for which the corresponding eigenvalue of the ball is simple. For comparison the corresponding eigenvalue of the ball with the same width is shown in Figure 3. As in the two dimensional study in [40], one could investigate the local minimality of the eigenvalues of the Dirichlet-Laplace operator on the ball in the class of constant width bodies. The numerical simulations use 900 spherical harmonics and 3000 points on where convexity constraints are imposed.

Minimizing the volume under constant width constraint -the Meissner conjecture
This section presents numerical approximations of solutions to Problem 3. The convexity and constant width constraints are imposed as indicated in Section 4. The area is explicit in terms of the Fourier coefficients in dimension two (see (11)). Moreover, the volume of constant width bodies is explicit in terms of the spherical harmonics coefficients in dimension three (see (13) and (15)). Therefore, in this case the functional, its gradient and the corresponding Hessian matrix can be computed explicitly, leading to quickly converging numerical algorithms.
The two dimensional result, the Reuleaux triangle, is shown in Figure 4. The minimal value for the area obtained with our algorithm, for N = 250, corresponding to 501 Fourier coefficients and width w = 2 is 2.8196. This is slightly larger than but very close to the explicit area of the Reuleaux triangle of width 2 which is equal to 2(π − √ 3) = 2.8191. The minimization of the volume under constant width constraint in dimension three is a famous open problem. The conjectured optimizer is a Reuleaux tetrahedron with three rounded edges. There are two configurations, with the same volume, the difference being in the position of the rounded edges: all starting from one vertex or forming a triangle. These shapes are called Meissner's bodies. Various works deal with the analysis of 3D shapes of constant width which minimize the volume. Among these we cite [42], which presents many aspects related to the Meissner bodies and why they are conjectured to be optimal. It is mentioned that in [43] the author generates a million random three dimensional bodies of constant width, using techniques from [6]. Among these many bodies of constant width, none had a smaller volume than the ones of Meissner. In [7] the local optimality of the Meissner's body was verified using an optimization procedure with a different parametrization of constant width shapes.
The approach we present below allows us to obtain the Meissner bodies as results of a direct optimization procedure, starting from random initializations. The formulas (13) and (15) allow us to write the volume as a quadratic expression of the coefficients of the spherical harmonics decomposition (12). As indicated by the Property 8, the constant width condition is imposed by fixing the first coefficient and considering only odd spherical harmonics in the decomposition. The convexity condition is achieved by using a discrete version of (5). We note that even though the convexity condition is non-linear, it is explicit enough such that we may compute its gradient. In this way, the minimization of the volume under constant width condition in dimension three becomes a constrained optimization problem of a quadratic functional with non-linear quadratic constraints. The optimization algorithm uses the explicit expression of the gradient and the Hessian matrix. For unit width, in view of [42], the Meissner bodies have volume equal to In our computation, using 402 spherical harmonics, i.e. using Legendre polynomials up to degree 28, we obtain the shape represented in Figure 4, with volume 0.4224. The shape obtained strongly resembles Meissner's body and its volume is about 0.6% larger than V M presented above. This may be due to the fact that singularities in the surface of the Meissner bodies are not sufficiently well approximated using the number of spherical harmonics above. Starting from different random initial coefficients we always arrive at shapes which are close to one of the two Meissner bodies [42].

Rotors of minimal volume
As underlined in Section 3, for some convex domains P there exist convex shapes ω which can rotate inside P while touching all its sides (or faces in dimension three). Therefore, when rotors exist, it makes sense to consider the problem of finding rotors of minimal volume. Theoretical aspects are recalled in the definition of Problem 4 and numerical computations for the two dimensional case are presented in [9]. In dimension three there exist rotors for the cube (constant width bodies), the regular tetrahedron and the regular octahedron. A characterization of rotors in terms of the coefficients of the decomposition of the support function can be found in [23]. More precisely we have the following: • In dimension two every regular polygon admits non-circular rotors. If the regular polygon has n sides, n ≥ 3, then only the Fourier coefficients of the support function for which the index has the form nq ±1 are non-zero, where q is a positive integer. • The rotors in a cube are bodies of constant width.
• The support function of a rotor in a regular tetrahedron has non-zero coefficients for the spherical harmonics of indices 0, 1, 2 and 5 • The support function of a rotor in a regular octahedron has non-zero coefficients for the spherical harmonics of indices 0, 1 and 5.
The constant term in the spectral decomposition of the support function of a rotor in P corresponds to the inradius of P . We note that when taking the midpoints of the edges of the regular tetrahedron we obtain a regular octahedron with the same inradius. Therefore rotors in the octahedron are also rotors for the tetrahedron which was already apparent from the characterization using the spherical harmonic coefficients.
Computations of optimal rotors in dimension two were also made in [9], while the computations in dimension three are new. Note that rotors of maximal area and volume are the inscribed disc and the inscribed ball, respectively. Some two dimensional rotors of minimal area are shown in Figure 5. Minimal volume rotors obtained numerically for the regular tetrahedron and the regular octahedron are shown in Figure 6. In each case we consider an optimization problem depending only on the non-zero coefficients describing the rotors and we impose discrete convexity constraints like in Section 4.  The computations presented in Figure 6 are made for solids with inradius equal to 0.5, corresponding to an inscribed ball of diameter 1. Compared to the volume of the ball B with unit diameter which is equal to π/6 = 0.5236 the minimal volume found numerically of a rotor in the tetrahedron and the octahedron circumscribed to the same ball B are 0.3936 and 0.5041, respectively. Numerical minimizers for the tetrahedron and octahedron seem to be symmetric under a rotation of angle 2π/5. This is due to the fact that the only coefficients which may change the geometry of the rotors correspond to the spherical harmonics of order 2 or of order 5. We recall that changing coefficients for spherical harmonics of order 1 corresponds to translations. In particular, when searching for minimal rotors in the tetrahedron using only spherical harmonics of order 1 and 2 we get a radially symmetric minimizer of volume 0.4024 which is slightly larger than the result including spherical harmonics of order 5.

Minimal width constraint
Problem 5 was first considered in [7]. The area of a three-dimensional convex body is minimized under a minimal width constraint. The minimal width constraint is characterized by the inequality p(θ)+p(−θ) ≥ 1, for every θ ∈ S 2 .
As underlined in Section 4.2, diameter bounds can be imposed in an approximate way using a finite number of points uniformly distributed on the sphere. This is done in the same way as the discrete convexity condition. These bounds on the diameter correspond to a set of linear inequality constraints. In Figure  7 we present the result given by the algorithm. The shape resembles the optimizer given in [7] and the value of the functional is slightly improved. In the computations M c = 2000 points are used for the discrete convexity condition and M d = 1000 pairs of opposite diametral points for computing the discrete diameter inequalities. We used 250 spherical harmonics in the decomposition of the support function. The computation of the area is explicit in terms of the spherical harmonics coefficients, as shown in (14). The minimal area obtained with our algorithm is 2.9154 which is slightly smaller than 2.9249, the value of the minimal area in the result presented in [7].

Inclusion constraint
In this Section we show how to impose inclusion constraints for shape optimization problems. As recalled in Section 2, two convex bodies B 1 , B 2 in R n , with support functions p B1 , p B2 , respectively satisfy the inclusion constraint B 1 ⊂ B 2 if and only if the support functions verify p B1 (θ) ≤ p B2 (θ) for every θ ∈ S n−1 .
As in the case of the convexity and diameter constraints, we impose (23) on a sufficiently dense discrete subset of S n−1 . In dimension three, when dealing with Cheeger sets for polyhedra, it is enough to impose the inclusion constraints only for directions which are normals to the faces of the polyhedron. This simplifies the optimization algorithm by decreasing the number of constraints. Cheeger sets provide a classical example of shape optimization problems involving convexity and inclusion constraints. The theoretical formulation is given in Problem 6. The Cheeger sets of Ω minimize the ratio perimeter/area (or surface area/volume) for convex subsets of Ω. This objective function can be computed and optimized using the proposed algorithm. As shown in Section 4, convexity and inclusion constraints are discretized as linear inequalities in terms of the coefficients of the spectral decomposition of the support function. Some examples of computation of Cheeger sets for the square in the plane and for the regular tetrahedron, the cube and the regular dodecahedron in dimension three are shown in Figure 8.

Conclusions
In this work, the properties of the support function are used to deal numerically with various constraints in shape optimization problems. The spectral decomposition of the support function using Fourier series in dimension two and spherical harmonics in dimension three are particularly well suited in order to discretize convexity, constant-width, diameter and inclusion constraints. The numerical tests use standard tools readily available in optimization software like quasi-Newton or Newton methods with linear or non-linear constraints and cover a wide variety of shape optimization problems with various constraints.