Longest minimal length partitions
LLongest minimal length partitions
Beniamin Bogosel and Edouard OudetFebruary 8, 2021
Abstract
This article provides numerical evidence that under volume constraint the ball is the set whichmaximizes the perimeter of the least-perimeter partition into cells with prescribed areas. We introducea numerical maximization algorithm which performs multiple optimizations steps at each iterationto approximate minimal partitions. Using these partitions we compute perturbations of the domainwhich increase the minimal perimeter. The initialization of the optimal partitioning algorithm usescapacity-constrained Voronoi diagrams. A new algorithm is proposed to identify such diagrams, bycomputing the gradients of areas and perimeters for the Voronoi cells with respect to the Voronoipoints.
In [15] the authors prove that among planar convex sets of given area the disk maximizes the length ofthe shortest area-bisecting curve. In order to write the result in a more formal way, let us introduce thefollowing notations. Denote by Ω ⊂ R d an open, connected region with Lipschitz boundary. Moreover,consider ω ⊂ Ω an arbitrary Lipschitz subset. Consider c ∈ (0 ,
1) and denote with | · | the usual Lebesguemeasure (area in 2D, volume in 3D). Given Ω and c , define the shortest fence set to be SF (Ω , c ) = argmin { Per Ω ( ω ) : | ω | = c | Ω |} . (1)In other words, SF (Ω , c ) is one subset ω ⊂ Ω which minimizes the relative perimeter Per Ω ( ω ) when themeasure | ω | is fixed to c | Ω | . For simplicity, denote with L (Ω , c ) = Per Ω ( SF (Ω , c )). The paper [15] citedabove solves the problem of maximizing L (Ω , c ) with respect to Ω,max | Ω | = v d L (Ω , c ) , (2)in dimension two for c = 1 /
2. In the following v d denotes the volume of the unit ball in R d . The choiceof the volume constraint | Ω | = v d does not reduce the generality of the problem, since changing thisconstant only rescales the solution via a homothety. Classical details regarding the existence of SF (Ω , c )and the definition of the relative perimeter Per Ω ( · ) are recalled in the next section.This work was initiated by the note [24] published on the French CNRS site Images des math´ematiques,where the author asks what happens to the solution of (2) when the parameter c varies in (0 , ⊂ R d and n > ω , ..., ω h ) to be a partition of Ω, in the sense that the unionof ω i , i = 1 , ..., n is ω and ω i ∩ ω j = ∅ . Given a vector c = ( c , ..., c n ) ∈ R n with (cid:80) ni =1 c i = 1 consider the shortest partition of | Ω | with area constraints c to be SP (Ω , c ) = argmin { n (cid:88) i =1 Per Ω ( ω i ) : ( ω i ) partition of Ω , | ω i | = c i | Ω |} . (3)Figure 1: Examples of minimizers of problem (1) for various shapes Ω and various constraints.1 a r X i v : . [ m a t h . O C ] F e b s before, denote P L (Ω , c ) to be the minimal total relative perimeter of SL (Ω , c ), i.e. the sum of relativeperimeters of members of SP (Ω , c ). Now it is possible to formulate the problemmax | Ω | = v d P L (Ω , c ) , (4)where the total relative perimeter of a partitions with constraints c is maximized when Ω has fixedvolume. It is immediate to see that (2) is a particular case of (4) by considering n = 2 and c = ( c, − c ).In this paper, problems (2) and (4) are investigated from both numerical and theoretical points ofview. In order to approximate solutions of these problems multiple issues need to be addressed: • Compute reliably a numerical approximation of the shortest partition SP (Ω , c ) once the domainΩ and the constraints vector c are given. It is important to avoid local minimizers at this stage,since the objective is to maximize the shortest perimeter. Any local minimizer may give a falsecandidate for the solutions of (2) or (4). There are many works in the literature which deal withthe investigation of minimal length partitions. In [11] the authors use the surface Evolver softwareto approximate minimal partitions. In this work the approach presented in [30] is used, where theperimeter is approximated using a Γ-convergence result of Modica and Mortola [28]. This allowsus to work with density functions rather than sets of finite perimeter and simplifies the handling ofthe partition condition. Moreover, working with densities directly allows changes in the topologyof the partitions.Given a domain Ω, a mesh is constructed and finite elements are used in FreeFEM [20] in orderto approximate SF (Ω , c ) or SP (Ω , c ). When dealing with partitions, in order to accelerate theconvergence, an initialization based on Voronoi diagrams of prescribed areas is used. • Once the shortest partition SP (Ω , c ) is identified, the bounding set Ω needs to be modified in orderto increase the objective function P L (Ω , c ). In order to find a suitable ascent direction classicalresults related to the shape derivative will be used [13, 21]. • The family of star-shaped domains is parametrized using radial functions, and radial functions arediscretized by considering truncation of the associated Fourier series. Using the shape derivativeit is possible to compute the gradient of the objective function with respect to the discretizationparameters. Once the gradient is known, an optimization algorithm is used in order to searchfor solutions of (2) and (4). The choice of the optimization algorithm is also an important factorsince the computation of SP (Ω , c ) is highly sensitive to local minima. Moreover, when changingΩ following a perturbation field found using a shape derivative argument, the configuration of theoptimal partition might change. The chosen algorithm is a gradient flow with variable step size.Minimal length partitioning algorithms presented in [30] or [5] use random initializations. While thisillustrates the flexibility of Modica-Mortola type algorithms and the ability of the algorithm to avoidmany local minima, choosing random initializations leads to longer computation times required for theoptimization algorithm. A classical idea is to use Voronoi diagrams as initializations. However, theseVoronoi diagrams should consist of cells which verify the area constraints | ω i | = c i . In the literature, thenotion of capacity-constrained Voronoi diagrams is employed and results in this direction can be found in[4], [3] and [34]. In this work we propose a new way of computing capacity-constrained Voronoi diagramsby explicitly computing the gradients of the areas of the Voronoi cells with respect to variations in theVoronoi points. The gradient of the perimeters of the Voronoi cells is also computed, which allows thesearch of capacity-constrained Voronoi diagrams with minimal length.The numerical simulations give rise to the following conjectures: • The result of [15] seems to generalize to every c ∈ (0 ,
1) in dimensions two and three. • For n > c = ( c i ) ni =1 ∈ R n with (cid:80) ni =1 c i = 1 arbitrary, the solution of (4) is the disk in 2Dand the ball in 3D. What is surprising is that this result holds even when the area constraints ofthe cells of the partition are not the same. Outline and summary of results.
Section 2 presents classical theoretical results regarding ap-proximations of minimal perimeter partitions by Γ-convergence.Section 3.1 recalls basic aspects regarding the numerical computation of minimal length partitions.Section 3.3 presents the computation of the gradients of the areas and perimeters of Voronoi cells andshows how to use prescribed-area Voronoi cells in order to construct initializations for our optimiza-tion algorithm. Section 3.4 presents the computation of a descent direction for the shape optimizationalgorithm using the notion of shape derivative. The choice of the discretization and the optimizationalgorithm for approximating solutions of problems (2) and (4) are presented in Section 3.5. We under-line that the maximization algorithm approximates solutions to a max-min problem, and the optimalpartitioning algorithm presented in Section 3.1 is run at every iteration.2inally, results of the optimization algorithm in dimensions two and three are presented in Section 4.The numerical results suggest that the solution of problems (2) and (4) is the disk in dimension two andthe ball in dimension three. A brief discussion of the optimality conditions is presented in Section 5.
The appropriate framework to work with sets of finite relative perimeter in Ω is to consider the spaceof functions with bounded variation on Ω BV (Ω) = { u ∈ L (Ω) : T V ( u ) < ∞} where T V ( u ) = sup { (cid:90) Ω u div g : g ∈ C c (Ω) , (cid:107) g (cid:107) ∞ ≤ } . As usual C c (Ω) represents the space of C ∞ functions defined on Ω with compact support in Ω. Usingthe divergence theorem it is easy to observe that if u is C ( S ) then T V ( u ) = (cid:90) Ω |∇ u | . If ω is a subset of Ω its generalized perimeter is defined by Per( ω ) = T V ( χ ω ), where χ ω represents thecharacteristic function of ω . All these aspects are classical and can be found, for example, in [2, 6].The fact that problems (1) and (3) have solutions is classical and is a consequence of the fact thatthe generalized perimeter defined above is lower-semicontinuous for the L convergence of characteristicfunctions. For more aspects related to solutions of these problems see [26, Chapter 17]. The bookpreviously referenced also presents aspects related to the regularity of optimal partitions in Part Four.Aspects about optimal partitions in the smooth case are presented in [29] where qualitative properties ofminimal partitions in the plane or on surfaces are presented.Proving existence of solutions for problems (2) and (4) is more difficult since these are maximumproblems and the perimeter is lower semicontinuous. We recall that problem (2) was solved in [15] in thecase d = 2, c = 1 /
2. In particular, existence was proved exploiting results in [10] which show that in thiscase the minimal relative perimeter sets are convex. In the following we prove that solutions exist in theclass of convex domains for arbitrary area constraints.
Theorem 2.1.
Problem (2) has solutions in the class of convex sets.Proof:
We divide the proof into steps which allow us to apply classical methods in calculus of varia-tions.
Step 1: Upper bounds.
In the following denote by w (Ω) the minimal H d − measure of theprojection of Ω on a hyperplane (in dimension two this corresponds to the minimal width). For convexbodies the following reverse Loomis-Whitney inequality holds true:min { e ,...,e d } d (cid:89) i =1 H d − ( K | e ⊥ i ) ≤ Λ d | K | d − , where the minimum is taken over all orthonormal bases of R d and K | e ⊥ i represents the projection of K onto a hyperplane orthogonal to e i . In [25] it is shown that there exists a constant c such thatΛ d ≤ ( c √ d ) d . In particular, this shows that the minimal projection w (Ω) verifies w (Ω) d ≤ Λ d | Ω | d − .As a direct consequence w (Ω) is bounded above in the class of convex sets Ω which satisfy | Ω | = v d .It is immediate to see that the quantity w (Ω) gives an upper bound for SF (Ω , c ). To justify thischoose e the direction for which H d − (Ω | e ⊥ ) is minimal and slice Ω with a hyperplane orthogonal to e which divides Ω into two regions ω and Ω \ ω with volume | ω | = c . The relative perimeter of the set ω in Ω is at most equal to w (Ω), the H d − measure of the projection. Therefore, we may conclude thatin the class of convex sets with measure | Ω | = v n the quantity L (Ω , c ) is bounded from above, and theupper bound only depends on d and v d . This implies the existence of a maximizing sequence (Ω h ) h ≥ which verifies L (Ω h , c ) ≤ L (Ω h +1 , c ) and L (Ω h , c ) → sup | Ω | = v d L (Ω , c ), where the supremum is taken inthe class of convex sets. 3 tep 2: Compactness. When dealing with a sequence of convex sets we may extract a subsequenceconverging in the Hausdorff distance provided the sets are uniformly bounded. For classical aspectsrelated to the Hausdorff distance we refer to [21, Chapter 2]. Therefore, in the following we show thatthe diameters diam(Ω h ) of convex sets Ω h forming the maximizing sequence are uniformly bounded.First, let us note that since (Ω h ) is a maximizing sequence for L (Ω , c ) there exists a positive constant p > L (Ω , c ) > p . Since w (Ω) ≥ L (Ω , c ) we also have w (Ω h ) ≥ p > n ≥
1. The resultsin [16] show that the minimal perimeter projection, the diameter and the volume of a convex set Ω satisfy w (Ω) diam(Ω) ≤ | Ω | /d. It is now immediate to see that diam(Ω h ) ≤ | Ω h | / ( dw (Ω)) ≤ v d / ( dp ), and therefore the diameters of(Ω h ) are bounded. Without loss of generality we may assume that (Ω h ) are contained in a large enoughball. Applying the classical Blaschke selection theorem we find that there exists a maximizing sequence,denoted for simplicity by (Ω h ), such that Ω h converges, with respect to the Hausdorff distance, to theconvex set Ω . Moreover, the volume is continuous for the Hausdorff distance among bounded convexsets, so Ω also satisfies the volume constraint | Ω | = v d . Step 3. Continuity.
The last step is to prove that L (Ω , c ) is indeed equal to lim sup n →∞ L (Ω h , c ).This is a direct consequence of [32, Theorem 4.1], which states that if (Ω h ) is a sequence of convex bodiesin R d and Ω h → Ω in the Hausdorff distance then L (Ω h , c ) → L (Ω , c ) for every c ∈ [0 , (cid:3) Remark 2.2.
Removing the convexity assumption is not straightforward. Nevertheless, using the reg-ularity results regarding solutions of (1) it is possible that this result could be partially extended in thegeneral case. There are multiple difficulties which follow the structure of the proof above: • Proving there exists an upper bound in (2) . • Proving that a maximizing sequence is bounded: long tails may not intersect the minimizing set in (1) therefore cutting them may increase L (Ω , c ) . • Obtaining compactness results of a maximizing sequence: classically this should be possible whenworking in the class of sets of finite perimeter. • Proving that the maximizing sequence converges to an actual maximizer. This would involve obtain-ing some continuity properties regarding the perimeter of a sequence of sets. This is not straight-forward, as the perimeter is only lower-semicontinuous for the L convergence of characteristicfunctions. Nevertheless, using the regularity of minimal relative perimeter sets might help obtainthe desired results. The case of partitions can be handled using a similar strategy in the class of convex sets. The missingingredient is the convergence of the minimal perimeters of partitions, analogue to the results in [32].
Theorem 2.3.
Problem (4) has solutions in the class of convex sets.Proof:
As in the proof of 2.1 it is straightforward to give upper bounds for
P L (Ω , c ) in terms of w (Ω)(the minimal H d − measure of the projection on a hyperplane). A maximizing sequence (Ω h ) wouldhave a positive lower bound 0 < p ≤ w (Ω h ) for the sequence of minimal projections on hyperplanes.Therefore the diameters of (Ω h ) are bounded from above and we may assume that the convex sets Ω h converge to a convex set Ω (with respect to the Hausdorff distance). The set Ω also verifies the volumeconstraint | Ω | = v d .It only remains to prove the continuity of the minimal perimeters P L (Ω h , c ) for the convergence withrespect to the Hausdorff distance. In order to do this, the same tools as in the proof of Theorem 4.1 in[32] can be used.1. Lower-semicontinuity.
Theorem 3.4 in [32] shows that there exist bilipschitz maps f h : Ω h → Ω with Lipschitz constants Lip( f h ) converging to 1, the Lipschitz constants of the inverse mapsLip( f − h ) also converging to 1. The volumes and perimeters of the images of finite perimeter sets E h ⊂ Ω h have upper and lower bounds as follows:1Lip( f − h ) d | E h | ≤ | f h ( E h ) | ≤ Lip( f h ) d | E h | f − h ) d − Per Ω h ( E h ) ≤ Per Ω ( f h ( E h )) ≤ Lip( f h ) d − Per Ω h ( E h )Let ( ω ih ) ni =1 be a minimal perimeter partition for Ω h with constraint c ∈ R n . Then ( f h ( ω ih )) is apartition of Ω with lim h →∞ | f h ( ω ih ) | = c i | Ω | . Extracting a diagonal sequence, we may assume that4 ω ih ) ni =1 converges with respect to the Hausdorff distance to a partition ( ω i ) ni =1 of Ω as h → ∞ .Using the estimates above and the fact that the perimeter is lower semi-continuous with respect tothe convergence of finite perimeter sets we have P L (Ω , c ) ≤ n (cid:88) i =1 Per Ω ( ω i ) ≤ lim inf h →∞ n (cid:88) i =1 Per Ω ( f i ( ω ih ))= lim inf h →∞ n (cid:88) i =1 Per Ω h ( ω ih ) = lim inf h →∞ P L (Ω h , c ) . Upper-semicontinuity.
It remains to prove that
P L (Ω , c ) ≥ lim sup h →∞ P L (Ω h , c ). Reasoningby contradiction, suppose that P L (Ω , c ) < lim sup P L (Ω h , c ). Up to a subsequence we may assumethat P L (Ω h , c ) converges. Choose ( ω i ) ni =1 a minimal partition in Ω with constraints | ω i | = c i | Ω | .As in [32] using these sets it is possible to construct better competitors on some Ω h for large h thanthe corresponding optimal partition. This leads to a contradiction.Indeed, ( f − h ( ω ih )) ni =1 forms a partition of Ω h , which may fail to satisfy the volume constraints.Optimality conditions imply that common boundaries of the sets in the partition are regular hy-persurfaces. Therefore, it is possible to perturb these boundaries around regular points in orderto attain the desired volume constraints. Moreover, for h large enough this will produce partitionswhich verify n (cid:88) i =1 Per Ω h ( f − h ( ω h )) < P L (Ω h , c ) , contradicting the optimality of P L (Ω h , c ).This concludes the proof of the existence of solutions for the given problem. (cid:3) Remark 2.4.
Existence results obtained in this section may also be generalized to the case where Ω isthe boundary of a convex set in R d . In particular, there exist sets Ω which are surfaces of co-dimension that are boundaries of some convex set in R d and have fixed H d − measure which maximize the minimalrelative geodesic perimeter of a subset or partition with given H d − measure constraints. A key point in our approach is to approximate minimal length partitions SP (Ω , c ). In order to avoiddifficulties related to the treatment of the partition constraint it is convenient to represent each set inthe partition ω i as a density u i : Ω → [0 , (cid:80) ni =1 u i = 1 on Ω. The next aspect is to approximate the perimeter of a setrepresented via its density function. A well known technique is to use a Γ-convergence relaxation for theperimeter inspired by a result of Modica and Mortola [28]. The main idea is to replace the perimeter witha functional that, when minimized, yields minimizers converging to those that minimize the perimeter.Let us briefly recall the concept of Γ-convergence and the property that motivates its use when dealingwith numerical optimization. Remark 2.5.
Let X be a metric space. For ε > consider the functionals F ε , F : X → [0 , + ∞ ] . Wesay that F ε Γ -converges to F and we denote F ε Γ −→ F if the following two properties hold:(LI) For every x ∈ X and every ( x ε ) ⊂ X with ( x ε ) → x we have F ( x ) ≤ lim inf ε → F ε ( x ε ) (5) (LS) For every x ∈ X there exists ( x ε ) ⊂ X such that ( x ε ) → x and F ( x ) ≥ lim sup ε → F ε ( x ε ) . (6)An important consequence is the following result concerning the convergence of minimizers of a se-quence of functionals that Γ converge. Proposition 2.6.
Suppose that F ε Γ −→ F and x ε minimizes F ε on X . Then every limit point of ( x ε ) isa minimizer for F on X . F it is possible to search for mini-mizers of F ε , for ε small enough.Let us now state the two theoretical results that are used in this work concerning the Γ-convergencerelaxation of the perimeter and of the total perimeter of a partition, with integral constraints on thedensities. The first result is the classical Modica-Mortola theorem [28]. Various proofs can be found in[1, 6, 9]. In the following Ω is a bounded, Lipschitz open set and let W : R → [0 , ∞ ) is a continuousfunction such that W ( z ) = 0 if and only if z ∈ { , } . For a given double well potential W describedpreviously, denote γ = 2 (cid:82) (cid:112) W ( s ) ds . In the following c ∈ [0 ,
1] represents the fraction used for thevolume constraint.
Theorem 2.7 (Modica-Mortola) . Define F ε , F : L (Ω) → [0 , + ∞ ] by F ε ( u ) = (cid:90) Ω (cid:18) ε |∇ u | + 1 ε W ( u ) (cid:19) u ∈ H (Ω) , (cid:82) Ω u = c | Ω | + ∞ otherwiseand F ( u ) = (cid:40) γ Per Ω ( u − (1)) u ∈ BV (Ω; { , } ) , (cid:82) Ω u = c | Ω | + ∞ otherwise . Then F ε Γ −→ F in the L (Ω) topology. In [30] this result was generalized to the case of partitions and was used to compute approximationsfor SP (Ω , c ). For c ∈ R n with (cid:80) ni =1 c i = 1, in order to simplify notations, let us denote by X (Ω , c )the space of admissible densities which verify the integral constraints and the algebraic non-overlappingconstraint X (Ω , c ) = { u = ( u i ) ni =1 ∈ L (Ω) n : (cid:90) Ω u i = c i | Ω | , n (cid:88) i =1 u i = 1 } . The Γ-convergence result in the case of partitions is recalled below.
Theorem 2.8.
Define P ε , G : L (Ω) → [0 , + ∞ ] by G ε ( u ) = n (cid:88) i =1 (cid:90) Ω (cid:18) ε |∇ u i | + 1 ε W ( u i ) (cid:19) if u ∈ ( H (Ω)) n ∩ X (Ω , c )+ ∞ otherwise G ( u ) = (cid:40) γ (cid:80) ni =1 Per Ω ( { u i = 1 } ) if u ∈ ( BV (Ω , { , } )) n ∩ X (Ω , c )+ ∞ otherwiseThen G ε Γ −→ G in the ( L (Ω)) n topology. A proof of this result can be found in [30]. In the numerical simulations the double well potential ischosen to be W ( s ) = s (1 − s ) which gives the factor γ = 1 / Remark 2.9.
It can be seen that SF (Ω , c ) is a minimizer of F in Theorem and SP (Ω , c ) is aminimizer of G in Theorem . Using the result recalled in Proposition it is possible to approximatethese minimizers by those of F ε and G ε , respectively, for ε small enough. From a numerical point ofview, dealing with the minimization of F ε and G ε is easier since the variable densities are H regular. Inthe following we denote by L ε (Ω , c ) and P L ε (Ω , c ) the optimal costs obtained when minimizing F ε and G ε respectively. Remark 2.10.
It can immediately be seen that, assuming W is at least of class C , minimizers u of F ε verify an optimality condition of the form (cid:90) Ω (cid:18) ε ∇ u · ∇ ϕ + 1 ε W (cid:48) ( u ) ϕ + µϕ (cid:19) = 0 , for every ϕ ∈ H (Ω) where µ ∈ R is a Lagrange multiplier for the volume constraint. Classical regularity theory results thatcan be found in [17] allow us to employ a bootstrap argument and conclude that u is of class C ∞ in theinterior of Ω and u has the regularity of Ω up to the boundary. For example, for smooth domains Ω the ptimizer u is also smooth up to the boundary of Ω . Moreover, it can be proved that the minimizer u takesvalues in [0 , . We refer to [18] for more details. Moreover, when Ω is convex it is possible to deducethat the minimizer u is Lipschitz continuous.The same type of result holds true for minimizers of G ε in the partition case, with eventual singularitiesat junction points between three phase or more in the partition. Nevertheless, the contact between theoptimal partition and the boundary ∂ Ω has the desired regularity. Since in the numerical section we deal with the minimization of F ε , G ε for fixed Ω and with themaximization of L ε (Ω , c ) , P L ε (Ω , c ), we briefly recall existence results related to these problems. In thefollowing we suppose that the double well potential W is Lipschitz continuous on R . This is not restrictivesince minimizers of F ε , P ε are densities which take values in [0 , W far awayfrom this interval do not matter in the analysis. Theorem 2.11. (i) Problems min u ∈ L (Ω) F ε ( u ) and min u ∈ L (Ω) n P ε ( u ) admit solutions for Ω a Lipschitz domain with finite volume.(ii) Given c ∈ (0 , and c = ( c i ) ∈ R n , (cid:80) ni =1 c i = 1 problems max | Ω | = v d L ε (Ω , c ) and max | Ω | = v d P L ε (Ω , c ) admit solutions in the class of convex sets.Proof: The proof of (i) is classical. Note that the constraints on the density functions are embedded inthe definition of the functionals F ε , P ε to be minimized. We give the ideas for P ε as F ε is just a particularcase. The existence proof goes as follows: • The functional P ε is obviously bounded from below by zero. Moreover, truncating the densityfunctions ( u i ) to take values in [0 ,
1] does not increase the value of P ε . This allows us from now onto assume that the densities have values in this interval. • Minimizing sequences exist and they are bounded in H (Ω) n , which allows us to extract a subse-quence weakly converging in H . The constraints are stable under the L convergence. Moreover,the lower-semicontinuity of the H norm and Fatou’s lemma allow us to see that any weak H -limitpoint of the minimizing sequence is a minimizer.The proof of (ii) follows the same lines as the proofs of Theorems 2.1, 2.3. As in the proof ofthese theorems, we start by noticing that the minimal H d − measure w (Ω) of the projection of Ω on ahyperplane is bounded from above. We detail the proof for L ε , while the proof in the case of partitionsfollows the same path. In order to emphasize the dependence of F ε on Ω we write F ε ( u ) = F ε (Ω , u ). Upper bound.
Choose e the direction for which H d − (Ω | e ⊥ ) is minimal and equal to w (Ω). Givena hyperplane ζ orthogonal to e consider the function u ε = ϕ ε ( d ( x )), where d ( x ) is the signed distanceto the hyperplane ζ (choosing an orientation) and ψ ε , ϕ ε are given by ψ ε ( t ) = (cid:90) t ε (cid:112) ε + W ( s ) ds, ϕ ε ( t ) = t ≤ ψ − ε ( t ) 0 ≤ t ≤ ψ ε (1)1 t ≥ ψ ε (1) . This type of construction is standard when proving the limsup part of the Γ-convergence proof for theModica-Mortola type results in Theorems 2.7, 2.8 (see for example [27]). The coarea formula and thefact that |∇ d ( x ) | = 1 allows us to write (cid:90) Ω u ε = (cid:90) R (cid:90) { d ( x )= t } ϕ ε ( t ) d H d − dt. The definition of ϕ ε and a continuity argument allow us to deduce that there is a position of the hyperplanefor which the constraint (cid:82) Ω u ε = c | Ω | is verified.Using the coarea formula to evaluate F ε (Ω , u ε ) we obtain F ε (Ω , u ε ) = (cid:90) Ω (cid:18) ε | ϕ (cid:48) ε ( d ( x )) | + 1 ε W ( ϕ ε ( d ( x ))) (cid:19) = (cid:90) R (cid:90) { d ( x )= t } (cid:18) ε | ϕ (cid:48) ε ( t ) | + 1 ε W ( ϕ ε ( t )) (cid:19) d H d − dt ≤ w (Ω) (cid:90) ψ ε (1)0 (cid:18) ε | ϕ (cid:48) ε ( t ) | + 1 ε W ( ϕ ε ( t )) (cid:19) dt. { d ( x ) = t } is a slice of Ω orthogonal to e and its H d − measure is at most w (Ω). Moreover, we can restrict the bounds in the one dimensional integral to 0 and ψ ε (1) since for t not in this interval the integrand is zero. A simple computation gives ϕ (cid:48) ε ( t ) = 1 ψ (cid:48) ε ( ψ − ε ( t )) = 1 ε (cid:112) ε + W ( ϕ ε ) . Thus we obtain F ε (Ω , u ε ) ≤ w (Ω) ε (cid:90) ψ ε (1)0 ( ε + W ( ϕ ε )) dt = 2 w (Ω) (cid:90) (cid:112) ε + W ( s ) ds, where the last equality comes from the change of variables s = ϕ ε ( t ). This quantity depends only on W and w (Ω) and is bounded from above independently of Ω. Compactness.
The same argument used in the proof of Theorem 2.1 can be applied in order toconclude that there exists a maximizing sequence (Ω h ) converging in the Hausdorff distance to a convexset Ω with non-empty interior and volume | Ω | = v d . Moreover, we may assume that there exists abounded open set D such that (Ω h ) h ≥ , Ω ⊂ D .Following the ideas in [21, Chapter 2] we may assume that (Ω h ) and Ω satisfy an ε -cone condition,or equivalently that they are Lipschitz regular with a uniform Lipschitz constant. In this case, theconvergence with respect to the Hausdorff distance implies that | Ω h \ Ω | + | Ω \ Ω h | → Continuity.
It now remains to prove that L ε (Ω h , c ) → L ε (Ω , c ) as h → ∞ . Let us note first that sinceΩ h is a maximizing sequence we have L ε (Ω , c ) ≤ lim h →∞ L ε (Ω h , c ). Consider a minimizer u ∈ H (Ω)such that F ε (Ω , u ) = L ε (Ω , c ).Since (Ω h ) and Ω have a uniform Lipschitz constant L (as recalled above), using the extension theoremsrecalled in [8, Theorem 3.4], there exists an extension ˜ u ∈ W ,p ( D ) of u which verifies (cid:107) ˜ u (cid:107) W ,p ( D ) ≤ Const( L ) (cid:107) u (cid:107) W ,p (Ω h ) . Together with the results recalled in Remark 2.10 we find that ε |∇ ˜ u | + ε W (˜ u ) ∈ L ∞ ( D ). Combining this with the fact that | Ω \ Ω h | + | Ω h \ Ω | → F ε (Ω h , ˜ u ) → F ε (Ω , u ) . We cannot conclude yet, since ˜ u may not satisfy the integral constraints on Ω h .In order to fix this, let x be a point in the interior of Ω. For h large enough there exists a ball B δ ofradius δ > B δ ⊂ Ω ∩ Ω h . Denote by d δ the function which is equal to the distance to ∂B δ inside B δ and zero outside. We use this function to construct functions u h = ˜ u + x h d δ , for x h ∈ R , whichverify the integral constraints (cid:82) Ω h u h = c | Ω h | . Since (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω ˜ u − (cid:90) Ω h u (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω \ Ω h ˜ u (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω h \ Ω u (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = O ( | Ω \ Ω h | + | Ω h \ Ω | ) → , we necessarily have x h →
0. This immediately shows that | F ε (Ω h , u h ) − F ε (Ω h , ˜ u ) | → h → ∞ .Since L ε (Ω h , c ) ≤ F ε (Ω , u h ) we find that lim sup h →∞ L ε (Ω h , c ) ≤ F ε (Ω , u ) = L ε (Ω , c ). This concludesthe proof of the existence result. The case of partitions can be handled in a similar manner with theadditional difficulty that the area constraints and sum constraints need to be handled simultaneously.This can be achieved by modifying the candidate densities in a finite family of balls. (cid:3) In this section the numerical minimization of F ε and P ε is discussed. Since Ω is a general domain inthis work, we choose to work with finite element discretizations. Given T h a triangulation of Ω, denoteby ( x j ) Nj =1 the set of the nodes. Working with P Lagrange finite elements, a piecewise affine function u defined on the mesh T h is written (cid:80) Nj =1 u j φ j . As usual, φ j are the piece-wise linear functions on eachtriangle, characterized by φ j ( x k ) = δ jk . For a P finite element function, the values u j are given by u ( x j )8nd we denote u = ( u j ) = ( u ( x j )) ∈ R N . With these notations, it is classical to introduce the massmatrix M and the rigidity matrix K defined by M = (cid:18)(cid:90) T h φ i φ j (cid:19) ≤ i,j ≤ N and K = (cid:18)(cid:90) T h ∇ φ i · ∇ φ j (cid:19) ≤ i,j ≤ N As an immediate consequence of the linearity of the decompositions u = (cid:80) Nj =1 u j φ j , v = (cid:80) Nj =1 v j φ j wehave that (cid:90) T h uv = u T M v and (cid:90) T h ∇ u · ∇ v = u T K v . This immediately shows that the functionals F ε and P ε can be expressed in terms of the mass and rigiditymatrices M and K using the expression (cid:90) T h (cid:18) ε |∇ u | + 1 ε u (1 − u ) (cid:19) = ε u T K u + 1 ε v T M v =: F ( u ) (7)where v = ( u j (1 − u j )) Nj =1 . The gradient of this expression w.r.t. u can be computed and is given by ∇F ( u ) = 2 εK u + 2 ε M v (cid:12) (1 − u ) , (8)where (cid:12) denotes pointwise multiplication of two vectors: u (cid:12) v = ( u j v j ) Nj =1 .It is obvious that with (7) and (8) it is possible to implement a gradient-based optimization algorithmin order to minimize F ε and P ε . The software FreeFEM [20] is used for constructing the finite elementframework and the algorithm LBFGS from the package Nlopt [23] is used for the minimization of (7).We address the question of handling the constraints in the next section. The area or volume constraint can be expressed with the aid of the vector m = M e with e =(1 , , ..., ∈ R N . Indeed, with this notation, for a finite element function u we have (cid:82) T h u = m · u . Projection for one phase.
Let us start with the projection of one function onto the integralconstraint. Given a P finite element function u and its values u at the nodes we search a function u = u + α m with values at nodes u verifying the constraint m · u = c by solving( u + α m ) · m = c, which leads to α = ( c − u · m ) / ( m · m ).An alternative way of handling the constraint during the optimization process is to project the initialvector on the constraint and project the gradient onto the hyperplane x · m = 0 at each iteration. Thiscan simply be done by using c = 0 in the relation above. Such a modification of the gradient allows usto use efficient black-box optimization toolboxes, since quasi-Newton algorithms like LBFGS will performupdates based on a number of gradients stored in memory. If these gradients verify x · m = 0, the integralconstraint will be preserved throughout the optimization process. Projection for multiple phases.
In the case of partitions projections on the integral constraintswere already proposed in [30] (when using finite differences) and in [5] (when using finite elements). Adrawback of using orthogonal projections parallel to the vector m is the fact that the vector u is modifiedalmost everywhere in the domain Ω, which also includes the regions where it is 0 or 1. As observedin [7], this can cause resulting optimal densities to be non-zero at interfaces between two cells and attriple points. The solution proposed was to use instead projections parallel to (cid:112) W ( u i ), which basicallyconsists in modifying the function u i only where it is different from 0 or 1.Let us now describe the construction of the projection algorithm on the constraints (cid:90) Ω u i = c i | Ω | , n (cid:88) i =1 u i = 1 . with the compatibility condition (cid:80) c i = 1. Consider λ ∈ H (Ω) and ( µ i ) ∈ R n and perform thetransformation u i + λ (cid:112) W ( u i ) + µ i (cid:112) W ( u i )9n order to satisfy the constraints (cid:90) Ω u i + (cid:90) Ω λ (cid:112) W ( u i ) + µ i (cid:90) Ω (cid:112) W ( u i ) = c i | Ω | , i = 1 , ..., n (9)and n (cid:88) i =1 u i + λ n (cid:88) i =1 (cid:112) W ( u i ) + n (cid:88) i =1 µ i (cid:112) W ( u i ) = 1 (10)It is easy to note that: • in view of (10), given µ i we can find λ : λ = 1 − (cid:80) ni =1 u i − (cid:80) ni =1 µ i (cid:112) W ( u i ) (cid:80) ni =1 (cid:112) W ( u i ) . • in view of (9), given λ we can find µ i using the relations above.In the following we introduce the quantities ¯ λ i = (cid:82) Ω λ (cid:112) W ( u i ). Again, in view of (9), if ¯ λ i are known,then µ i are known and so is λ . In order to obtain a system for ¯ λ i , let us note that µ i = c i | Ω | − (cid:82) D u i − ¯ λ i (cid:82) D (cid:112) W ( u i ) . With this in mind we get¯ λ i = (cid:90) Ω λ (cid:112) W ( u i ) = (cid:90) Ω (cid:32) − (cid:80) nj =1 u j − (cid:80) nj =1 µ j (cid:112) W ( u j ) (cid:80) nj =1 (cid:112) W ( u j ) (cid:33) (cid:112) W ( u i )= (cid:90) Ω − (cid:80) nj =1 u j − (cid:80) nj =1 (cid:18) c j | Ω |− (cid:82) Ω u j − ¯ λ j (cid:82) D √ W ( u j ) (cid:19) (cid:112) W ( u j ) (cid:80) nj =1 (cid:112) W ( u j ) (cid:112) W ( u i )In order to further simplify the above expression, let’s make the following notations • (cid:112) W ( u i ) = w i • − (cid:80) nj =1 u i = E • c j | Ω | − (cid:82) D u j = F j , j = 1 , ..., n This gives ¯ λ i = (cid:90) Ω E − (cid:80) nj =1 (cid:16) F j − ¯ λ j (cid:82) Ω w j (cid:17) w j (cid:80) nj =1 w j w i = (cid:90) Ω E − (cid:80) nj =1 F j w j / (cid:82) Ω w j (cid:80) nj =1 w j w i + n (cid:88) j =1 (cid:90) Ω w i w j / (cid:82) Ω w j (cid:80) nj =1 w j ¯ λ j which can be written in the form ( I − A )¯ λ = b with ¯ λ = (¯ λ , ..., ¯ λ n ) and A = (cid:32)(cid:90) Ω w i w j / (cid:82) Ω w j (cid:80) nj =1 w j (cid:33) i,j =1 ,...,n , b = (cid:32)(cid:90) Ω E − (cid:80) nj =1 F j w j / (cid:82) Ω w j (cid:80) nj =1 w j w i (cid:33) j =1 ,...,n . One may note that the above system ( I − A )¯ λ = b is singular since the sum on the columns of A is equalto 1 and therefore ( I − A T ) e = 0 , where e = (1 , ..., ∈ R n . This is due to the fact that one of the constraints is redundant, in view of thecompatibility condition. In practice we simply discard one unknown and set it to zero.As noted previously, the same procedure can be applied to the gradients g i associated to each u i inorder to satisfy at every iteration (cid:90) Ω g i = 0 , n (cid:88) i =1 g i = 0 . This allows us to preserve the constraints when using a black-box LBFGS optimizer when initial param-eters satisfy the integral and sum constraints. 10 .3 Initializations for 2D partitions - Voronoi diagrams
The optimization algorithm for approximating L ε (Ω , c ) and P L ε (Ω , c ) is ready to be implemented,following the ideas shown in the previous sections. There is, however, the choice of the initialization whichis non-trivial and which has an impact on the performance of the optimization algorithm. It was alreadynoted in [30] and [5] that starting from random initializations is possible, but some additional work needsto be done in order to avoid constant phases, which are encountered at some local minimizers. Keepingin mind that the optimal partition problem needs to be solved multiple times during the optimizationalgorithm, we propose below a different initialization strategy, based on Voronoi diagrams. The use ofVoronoi diagrams for generating initializations is a rather natural idea when dealing with partitions andwas already mentioned in [11]. In this section Ω is assumed to be a polygon.Using random Voronoi diagrams is not very helpful, since area constraints are not verified in general.This led us to consider Voronoi diagrams which verify the area constraints, which in the literature arecalled capacity-constrained Voronoi diagrams. Algorithms for computing such diagrams were proposed in[4] for the discrete case and in [3] for the continuous case. In the continuous case the method employedin [3] was to optimize the position of one Voronoi point at a time using the gradient-free Nelder-Meadmethod. In [34] the authors propose efficient ways of generating such diagram, but for weighted Voronoidiagrams only. In the following we propose an alternative method for constructing capacity-constrainedVoronoi diagrams by computing the sensitivity of the areas of the Voronoi cells with respect to the positionof the points generating the respective Voronoi diagram. Since we are also interested in minimizing theperimeter, the computation of the sensitivity of the perimeter of Voronoi cells is also described.
Terminology related to Voronoi diagrams.
Given a set of points p , ..., p n ∈ R (called Voronoipoints) the associated Voronoi diagram consists of n Voronoi cells V , ..., V n defined for i = 1 , ..., n by V i = { x ∈ R : | x − p i | ≤ | x − p j | , j = 1 , ..., n, j (cid:54) = i } . Each Voronoi cell V i is a polygonal region (possibly infinite). The vertices of V i are simply called vertices in the following (please observe the difference between the Voronoi points and the
Voronoi vertices ). Theedges of V i are called ridges , some of which can be unbounded. Each ridge connects two Voronoi vertices(possibly at infinity, for unbounded ridges) called ridge vertices . Moreover, each ridge separates two ofthe initial points, called ridge points . All structure information of a Voronoi diagram associated to a setof points can be recovered as an output to some freely available software like scipy.spatial.Voronoi .The Voronoi diagrams are not restricted to a bounded domain. It is possible, however, to considerrestrictions of a Voronoi diagram to a bounded set Ω by simply intersectiong the regions V i with the setΩ. In our implementation the intersection of polygons is handled using the Shapely
Python package forcomputational geometry.In the following, given the points p i , i = 1 , ..., n , we consider the Voronoi regions restricted to a finitedomain Ω re-defined by V i = V i ∩ Ω. Note that in some cases, some V i may be void if Ω does not containthe associated point p i . We explain below how to compute the gradients of the areas and perimeters of V i with respect to positions of the points p i . Gradient of the areas of the Voronoi cells.
The derivative of a functional that can be representedas an integral over the Voronoi cell V i with respect to the Voronoi points can be computed if the normaldisplacement of the cell is known. This fact was recalled in [14] and [22] and is classical in the shapederivative theory. However, since the functionals considered there were sums over all Voronoi cells V i , thecontributions coming from the variations of the boundary cancelled themselves and only the variation ofthe integrand mattered.This is no longer the case in our situation. The area of the voronoi cell V i is A i = (cid:82) V i dx and itsdirectional derivative when perturbing a point p j in direction d is given by the integral on the boundaryof the normal variation of V i : A (cid:48) i ( d ) = (cid:82) ∂V i θ.n , where θ is the infintesimal displacement of the boundaryof V i when moving p j in the direction d i . More explicitly, if V i ( t ) is the Voronoi cell for p j + td then θ = lim t → v ( t ) t , where the vector field v ( t ) is defined by V i + v ( t ) = V i ( t ) on the boundary of V i .For a given ridge v k v l with associated ridge points p i , p j , we perturb the point p i (cid:55)→ p i + δ andinvestigate the derivative of the normal perturbation of v k v l as δ →
0. The two main perturbations arethe following: • δ is collinear with p i p j : in this case the perturbation induced on the ridge is just δ/
2. (see Figure2 (left)). The associated infinitesimal normal perturbation is constant equal to 1 / • δ is orthogonal to p i p j : in this case the infinitesimal perturbation induced on the ridge is a rotationaround the intersection m ij of p i p j and v k v l . The associated infinitesimal normal perturbation varieslinearly on v k v l from −| v l − m ij | / | p i − p j | to | v k − m ij | / | p i − p j | (the signs vary with respect to11he orientation of the orthogonal perturbation). (see Figure 2 (middle)). In order to prove this itis enough to consider the normal perturbation v ( t ) of the ridge v k v l illustrated in Figure 2 (right)and take the limit v ( t ) /t as t → δ of p i we denote by n the normal vector to v k v l pointing outwards to V i and by t the unit vector collinear with −−→ v l v k . Furthermore, consider the notations for the normal andtangential contributions (computed as one dimensional integrals on v k v l of the infinitesimal perturbationsdescribed above): z n = ( δ · n ) / | v k − v l | , z t = ( δ · t ) 12 | p i − p j | ( | v k − m ij | − | v l − m ij | ) . (11)By symmetry, these contributions will be similar, but with changed signs when perturbing p j with δ . Thecontributions to the gradients of the areas of the cells V i and V j when perturbing p i or p j are describedbelow V i V j p i z n + z t − z n − z t p j z n − z t − z n + z t (12)The algorithm for computing the gradient of the areas of the cells simply iterates over all the Voronoiridges that intersect Ω and for each ridge adds the contributions described in (12).Algorithm 1 describes the computation of the gradient of the areas of the cells. The coordinates ofthe n input points are given in the vector x ∈ R n and the output is the real matrix M of size 2 n × n containing as columns the gradients of the areas of the n cells with respect to the 2 n coordinates. Algorithm 1
Compute gradients of areas of Voronoi cells
Require: x = ( x , y , ..., x n , y n ), coordinates of points p , ..., p n , bounding polygon Ω Initialize M = 0 (of size 2 n × n ) Compute the Voronoi diagram associated to the points ( p i ) ni =1 and the intersections of the Voronoicells with the polygon Ω. Set
Voronoi ridges as the set of Voronoi ridges that intersect the bounding polygon Ω. for r in Voronoi ridges do For the Voronoi ridge r Find the associated Voronoi points p i and p j and the Voronoi vertices v k , v l . Set δ = (1 ,
0) and compute the contributions z n , z t as above. Perform the updates using (12): M i − ,i ← M i − ,i + z n + z t , M i − ,j ← M i − ,j − z n − z t ,M j − ,i ← M j − ,i + z n − z t , M j − ,j ← M j − ,j − z n + z t . Set δ = (0 ,
1) and compute the contributions z n , z t as above. Perform the updates using (12): M i,i ← M i,i + z n + z t , M i,j ← M i,j − z n − z t ,M j,i ← M j,i + z n − z t , M j,j ← M j,j − z n + z t . end forreturn M The explicit formulas for the gradients of the areas allow us to easily find capacity-constrained Voronoidiagrams as results of an optimization algorithm. For given constraints | V i | = c i | Ω | with (cid:80) ni =1 c i = 1, itis enough to minimize the functional( p , ..., p n ) (cid:55)→ n (cid:88) i =1 (Area( V i ) − c i ) . (13)In order to obtain more regular structures it is also possible to minimize the energy( p , ..., p n ) (cid:55)→ n (cid:88) i =1 (cid:90) V i | x − p i | (14)12 k v l p j p i v k v l p j p i v k v l p j p i Figure 2: Normal perturbation of the Voronoi ridge when moving one of the point in the normal andtangent directions.Figure 3: (top) Voronoi diagrams with more than 100 cells with equal areas obtained when minimizing(13). (bottom) Voronoi diagram obtained when minimizing (14) under capacity constraints.under the capacity constraints | V i | = c i . The energy (14) is employed for characterizing centroidal Voronoidiagrams where each Voronoi point p i coincides with the centroid of the cell V i . In particular, CentroidalVoronoi diagrams are critical points for (14). See [34] for more details regarding this functional. Examplesin this sense are shown in Figure 3. The constrained minimization is done using the MMA algorithm[33] from the NLOPT library [23]. Note that all constraints are coded as inequality constraints in thisalgorithm: | V i | ≤ c i . Since V i form a partition of Ω it is immediate to see that if the sets satisfy theinequality constraints, they, in fact, also satisfy the equality constraints | V i | = c i . Remark 3.1.
It is also possible to generalize the gradient formulas when a density is involved, whendealing with quantities of the type (cid:82) V i ρ , where ρ ∈ L (Ω) is a given density. The shape derivative of (cid:82) V i ρ is (cid:82) ∂V i ρθ.n , where θ is the perturbation of the boundary ∂V i . The boundary perturbations are obviously thesame, but the computations in (11) are no longer explicit, and a one-dimensional numerical integrationneeds to be performed for each Voronoi ridge. Gradient of the perimeter of the Voronoi cells.
We saw that in order to compute the gradientof the areas of the Voronoi cells, the normal displacement of the Voronoi ridges needed to be understood,when moving the Voronoi points. On the other hand, the variation of the perimeter of a Voronoi regiondepends on the tangential perturbation of the Voronoi ridges. In order to understand this perturbationone needs to see how the Voronoi vertices move when perturbing the Voronoi points. Moreover, it can13e observed that when two Voronoi vertices merge, i.e. a Voronoi ridge collapses, the total perimeter ofthe cells is not smooth. This behavior is illustrated by an example shown in Figure 4.
Figure 4: Variation of the perimeter corresponding to a four point singularity in a square. The Voronoipoints are ( − t, , (0 , − , ( t, , (0 ,
2) for t ∈ [1 . , . A x , A y ) , ( B x , B y ) , ( C x , C y ): O x = 1 D [( A x + A y )( B y − C y ) + ( B x + B y )( C y − A y ) + ( C x + C y )( A y − B y )] O y = 1 D [( A x + A y )( C x − B x ) + ( B x + B y )( A x − C x ) + ( C x + C y )( B x − A x )] (15)where D = 2[ A x ( B y − C y ) + B x ( C y − A y ) + C x ( A y − B y )]. The formulas above are well defined as long asthe three points A, B, C are not colinear. Moreover, it is immediate to see that in this case the circum-center varies smoothly with respect to the coordinates of the vertices of the triangle. The infinitesimalperturbation of the circumcenter when moving ( A x , A y ) can be computed by simply differentiating theabove formulas w.r.t. A x and A y . Once the derivative of the circumcenter is known, in order to find thegradient of the prerimeter it is enough to project this derivative on all the Voronoi ridges going throughthe respective circumcenter and add the contribution to the gradient of the perimeter of each cell withrespect to the corresponding coordinates. See Figure 5 for more details. p i p j p k δ V p i p j p (cid:48) j δ V(cid:96) Figure 5: (left) Perturbation of the circumcenter when moving one point and projections on the Voronoiridges. (right) Computing the perturbation of a boundary point by transforming it into a circumcenter.Variations induced by the Voronoi nodes are enough to compute the gradient of perimeters of Voronoicells that do not intersect the boundary of the bounding polygon. For the boundary cells, it is necessary todescribe the perturbation of intersections between Voronoi ridges and the bounding polygon. Fortunately,this can also be described using variations of circumcenters for some particular triangles.Indeed, let v k v l be a Voronoi ridge intersecting a side (cid:96) of the bounding polygon Ω at the point q andlet p i , p j be the associated Voronoi points. Consider now p (cid:48) j the reflection of the point p j with respect14o the line supporting (cid:96) . Then obviously q is the circumcenter of the triangle p i p j p (cid:48) j and the variationof q with respect to perturbations of p i can be found using the same procedure as above. See Figure 5for more details. The algorithm for computing the gradients for the perimeters of the Voronoi cells ispresented in Algorithm 2, assuming that every Voronoi vertex is a circumcenter of exactly one triangledetermined by the Voronoi points. Algorithm 2
Compute gradients of perimeters of Voronoi cells
Require: x = ( x , y , ..., x n , y n ), coordinates of points p , ..., p n , bounding polygon Ω Initialize M = 0 (of size 2 n × n ) Compute the Voronoi diagram associated to the points ( p i ) ni =1 and the intersections of the Voronoicells with the polygon Ω. Set
Voronoi vertices as the set of Voronoi ridges that intersect the bounding polygon Ω. for v in Voronoi vertices do Let p i , p j , p k be the three Voronoi points which are associated to ridges going through v . Compute the derivative (cid:126)d of the circumcenter of p i p j p k when moving p i in the direction δ = (1 , For all ridges r going through v project (cid:126)d on r and add this to the gradient w.r.t. the x coordi-nate of the perimeter of the cells { V , V } neighbors to the ridge r (determined by the ridge pointsassociated to the ridge r ): these are elements M i − ,V , M i − ,V in matrix M . Repeat the above with δ = (0 ,
1) in order to get the gradients with respect to the y -coordinates. Do the same instructions for p j and p k . end for Set
Voronoi ridges as the set of Voronoi ridges that intersect the boundary polygon Ω. for r in Voronoi ridges do Denote by p i , p j the associated ridge points and by (cid:96) the edge of the boudnary polygon Ω cut by r Let p (cid:48) j be the reflection of p j with respect to (cid:96) . For δ = (1 ,
0) compute the derivative (cid:126)d of the circumcenter of p i p j p (cid:48) j when moving p i in thedirection δ . See Figure 5. Project (cid:126)d on the ridge r and add this projection to the gradient of the cells i and j w.r.t. the x coordinate: M i − ,i and M i − ,j in matrix M . Project (cid:126)d on (cid:96) and add this to the gradient of the cells i and j (with the proper sign). Repeat the above with δ = (0 ,
1) in order to get the gradients with respect to the y coordinates. Do the same instructions for p j . end forreturn M Using the gradients for the area and perimeters of Voronoi cells it is possible to perform a constrainedminimization of the perimeter under area constraint starting from random Voronoi initializations. Theoptimization is performed with Nlopt [23] optimization toolbox in Python using the MMA [33] algorithm.Some examples of initializations obtained are shown in Figure 6. The initial Voronoi points are chosenrandomly inside the polygon Ω. In order to accelerate the convergence of the optimization algorithm a fewiterations of Lloyd’s algorithm are performed before starting the optimization process. Recall that Lloyd’salgorithm consists in replacing the Voronoi points by the centroids of the respective cells iteratively (seefor example [34] for more details). In order to deal with local minima multiple optimizations (typically 10)are performed for every polygon Ω and the one with the partition having the least perimeter is retainedas a valid initialization. Note that the algorithm gives similar topologies with the best known ones shownin [11] for the case of equal areas and in [19] for the case of cells with two different areas.
Initialization of a partition.
Having at our disposal the gradients of areas and perimeters ofVoronoi cells, we are now ready to propose initialization algorithms for optimal partitioning algorithm.In practice we use one of the options below:1. Compute minimizers of (13) starting from random Voronoi points ( p i ). Repeat the procedure anumber of times and keep the configuration having the smallest total perimeter. This works wellwhen the areas of the cells are the same.2. Optimize the total perimeter of the Voronoi cells under capacity constraints starting from randomVoronoi points ( p i ). Repeat the procedure a number of times and keep the configuration having thesmallest total perimeter. This approach gives good results when the areas of the cells are differentand the optimization process is more difficult, since more local minima are present.15. When n ≤ n ≤ n ≥ In order find perturbations of the domain Ω that decrease the value of a functional J (Ω) the conceptof shape derivative is useful. For a Lipschitz domain Ω, the functional J is said to be shape differentiableat Ω if there exists a linear form θ (cid:55)→ J (cid:48) (Ω)( θ ) such that for every vector field θ ∈ W , ∞ ( R d , R d ) we have J (( I + θ )(Ω)) = J (Ω) + J (cid:48) (Ω)( θ ) + o ( (cid:107) θ (cid:107) W , ∞ ) , where I denotes the identity mapping. Classical results from [13, 21] show that for a function f ∈ H (Ω)the functional J (Ω) = (cid:82) Ω f is shape differentiable with J (cid:48) (Ω)( θ ) = (cid:90) ∂ Ω f θ.n. (16)In the following, ( · ) (cid:48) Ω denotes the derivation with respect to the domain Ω. The case of one phase.
Consider u Ω which minimizes F ε ( u ) from Theorem 2.7 and suppose that u Ω is unique. In this case we may assume that u Ω varies smoothly with respect to perturbations of theboundary of Ω. Remark 2.10 underlines the fact that u Ω is C ∞ in the interior of Ω and has the regularityof Ω up to the boundary. In the following we suppose that Ω is at least of class C , which implies that u Ω is indeed in H (Ω) and the gradient ∇ u Ω has a well defined trace on ∂ Ω. Differentiating F ε ( u Ω )with respect to Ω gives two terms, one involving just the derivation of a volumic integral and the otherinvolving the shape derivative of the optimal density u Ω with respect to θ . Fortunately, the minimalityof u Ω eliminates the latter term, which is more delicate and depends implicitly on the perturbation θ .Indeed, L ε (Ω , c ) (cid:48) Ω ( θ ) = ( F ε ( u Ω )) (cid:48) Ω ( θ ) = (cid:90) ∂ Ω (cid:18) ε |∇ u Ω | + 1 ε W ( u Ω ) (cid:19) θ.n + (cid:90) Ω (cid:18) ε |∇ u Ω | + 1 ε W ( u Ω ) (cid:19) ( u (cid:48) Ω ( θ )) (equal to 0)= (cid:90) ∂ Ω (cid:18) ε |∇ u Ω | + 1 ε W ( u Ω ) (cid:19) θ.n (17)where u (cid:48) Ω ( θ ) is the shape derivative of u Ω with respect to θ . For more details one may consult [13].In order to understand why the term above involving u (cid:48) Ω ( θ ) is zero let us recall that minimizing F ε ( u )under the constraint (cid:82) Ω u = c implies the existence of a Lagrange multiplier µ ∈ R such that (cid:90) Ω (cid:18) ε ∇ u Ω · ∇ φ + 1 ε W (cid:48) ( u Ω )( φ ) + µφ (cid:19) = 0 , for every φ smooth enough. In the numerical simulation we consider W ( s ) = s (1 − s ) which, in particularverifies W (cid:48) (0) = W (cid:48) (1) = 0. For ε small enough, the set { u Ω ∈ { , }} has non-empty interior. Moreover, ∇ u Ω = 0 and W (cid:48) ( u Ω ) = 0 on { u Ω ∈ { , }} , so choosing φ with compact support in { u Ω ∈ { , }} gives µ = 0. As a direct consequence, given a vector field θ we have (cid:90) Ω (cid:18) ε |∇ u Ω | + 1 ε W ( u Ω ) (cid:19) (cid:48) ( u (cid:48) Ω ( θ )) = 016 .51.01.52.02.5 Figure 7: (left) Minimization of the Modica-Mortola functional F ε with integral constraint 0 . | Ω | . (right)Representation of the shape gradient and the normal perturbation producing an ascent direction for SF (Ω , . (cid:0) ε |∇ u Ω | + ε W ( u Ω ) (cid:1) is non-zero (and strictly positive) only in the neighborhood of the contact pointsof the minimal relative perimeter set with the boundary ∂ Ω. Moving the boundary outwards at thesepoints with a small enough step size will increase the minimal perimeter L ε (Ω , c ).In the case the solution u Ω is not unique, we cannot assume that u Ω varies smoothly with Ω. Indeed,perturbing the contact points of (cid:0) ε |∇ u Ω | + ε W ( u Ω ) (cid:1) may drastically change the topology of the minimalset. Nevertheless, perturbing the boundary of Ω with this normal velocity will either increase the valueof L ε (Ω , c ) or keep it constant (in case there are multiple solutions). Moreover, the perturbed domainwill have a reduced multiplicity for the family of optimal solutions u Ω . Assuming there are only finitelymany solutions for a given Ω, repeating this process will eventually decrease the objective function. InFigure 7 the numerical approximation of u Ω is shown together with the value of (cid:0) ε |∇ u Ω | + ε W ( u Ω ) (cid:1) .Perturbing Ω in the normal direction as shown will increase the minimal value, provided the solution u Ω is unique. Remark 3.2.
The hypothesis regarding the uniqueness of u Ω for the shape derivative to exist is similarto the hypothesis needed when differentiating the eigenvalue of an operator with respect to the shape.When dealing with multiple eigenvalues, the shape derivative does not exist, but directional derivativesare available. For more details see [21, Chapter 5]. It is possible that such theoretical results could beobtained in our case, but this goes outside the scope of this article. The case of partitions.
Consider u Ω = ( u i Ω ) ni =1 which minimizes P ε in Theorem 2.8. As in theprevious paragraphs, supposing that u Ω is unique and varies smoothly with respect to perturbations inΩ, using again the minimality of u Ω for P ε we find that P L ε (Ω , c ) (cid:48) Ω ( θ ) = ( P ε ( u Ω )) (cid:48) Ω ( θ ) = (cid:90) ∂ Ω n (cid:88) i =1 (cid:18) ε |∇ u i Ω | + 1 ε W ( u i Ω ) (cid:19) θ.n + (cid:90) Ω (cid:34) n (cid:88) i =1 (cid:18) ε |∇ u i Ω | + 1 ε W ( u i Ω ) (cid:19)(cid:35) (cid:48) ( u (cid:48) Ω ( θ )) (equal to 0)= (cid:90) ∂ Ω n (cid:88) i =1 (cid:18) ε |∇ u i Ω | + 1 ε W ( u i Ω ) (cid:19) θ.n. (18)Again, in order to justify the fact that the term underlined above is indeed zero, let us recall that whenminimizing P ε ( u ) under the constraints (cid:82) Ω u i = c i | Ω | and (cid:80) ni =1 u i = 1 we find that there exist Lagrangemultipliers µ i ∈ R , i = 1 , ..., n and λ ∈ L (Ω) such that (cid:90) Ω (cid:32) n (cid:88) i =1 ε ∇ u i Ω · ∇ φ i + 1 ε W (cid:48) ( u i Ω )( φ i ) + µ i φ i + λ ( φ + ... + φ n ) (cid:33) = 0 . Following similar arguments as in the case of one phase, we find that the Lagrange multipliers associated17
Figure 8: (left) Minimization of the minimal partition functional P ε for three cells with equal areas.(right) Representation of the shape gradient and the normal perturbation producing an ascent directionfor SF (Ω , c ).Figure 9: Symmetric domain with two minimizers L ε (Ω , . (cid:90) Ω (cid:34) n (cid:88) i =1 (cid:18) ε |∇ u i Ω | + 1 ε W ( u i Ω ) (cid:19)(cid:35) (cid:48) ( u (cid:48) Ω ( θ )) = − (cid:90) Ω λ (cid:32) n (cid:88) i =1 ( u i Ω ) (cid:48) ( θ ) (cid:33) . On the other hand, differentiating the constraint (cid:80) ni =1 u i Ω = 1 gives (cid:80) ni =1 ( u i Ω ) (cid:48) ( θ ) = 0, which proves thedesired result.As discussed above, in the case u Ω is not unique, perturbing the boundary of Ω with normal velocitygiven by (cid:80) ni =1 (cid:0) ε |∇ u i Ω | + ε W ( u i Ω ) (cid:1) will not decrease the value of P L ε (Ω , c ) and will reduce the multi-plicity of the family of optimal partitions u Ω . Repeating the process will eventually increase the length ofthe minimal partition. In Figure 8 the numerical approximation of u Ω is shown together with the valueof (cid:80) ni =1 (cid:0) ε |∇ u i Ω | + ε W ( u i Ω ) (cid:1) . Perturbing Ω in the normal direction as shown will increase the length ofthe minimal partition, provided the solution u Ω is unique. An example of non differentiability - multiplicity greater than one.
It is not difficult toimagine domains Ω for which there are multiple minimizers L ε (Ω , c ), given c >
0. It is enough to consider c small enough and a symmetric domain like in Figure 9. Let us show that in such a case the functional L ε (Ω , c ) does not admit a shape derivative. Indeed, suppose that Ω is symmetric like in Figure 9. Denote u , u the two solutions and θ , θ two vector fields such that θ i .n = ε |∇ u i | + ε u i (1 − u i ) (also illustratedin Figure 9).Suppose that J (Ω) := L ε (Ω , .
2) is differentiable at Ω. Then for { i, j } = { , } it is clear that F ε (( I + tθ i )(Ω) , u j ) = F ε (Ω , u j ), i.e. the minimal value of F ε does not change when modifying Ω withonly one of the vector fields θ i . This would imply that J (cid:48) (Ω)( θ ) = J (cid:48) (Ω)( θ ) = 0 and by linearity J (cid:48) (Ω)( θ + θ ) = 0. However, this last equality is clearly false, since if we modify Ω with the combinedvector field θ + θ we clearly have J (( I + ( θ + θ )(Ω)) = J (Ω) + (cid:90) ∂ Ω (cid:18) ε |∇ u | + 1 ε u (1 − u ) (cid:19) θ .n + o ( (cid:107) θ (cid:107) W , ∞ ) .
18n order to deduce the above equality it is enough to work with one half Ω i of the domain Ω. We use thefact that this half can be extended to a non-symmetric domain for which u i is a unique minimizer andapply the formula for the shape derivative found previously.Therefore, we arrive at a contradiction, showing that when multiple minimizers of F ε (Ω , u ) exist, thefunctional L ε (Ω , c ) is not shape differentiable. The same kind of argument can be applied for P L ε (Ω , c ). The results of [15] and existence results obtained in Section 2.1 are restricted to convex domainsΩ. We therefore choose to search for domains maximizing SF (Ω , c ) and SP (Ω , c ) in the larger class ofstar-shaped domains which includes the class of convex sets. These domains can be parametrized using aradial function in dimension two and three. Furthermore, a spectral decomposition of the radial functionwith a Finite number of Fourier coefficients is used in order to work with a finite, but sufficiently largenumber of parameters in the computations. Planar domains.
In dimension two, the radial function ρ : [0 , π ] → R + is discretized using 2 N + 1Fourier coefficients ρ ( t ) = a + N (cid:88) k =1 ( a k cos( kt ) + b k sin( kt )) . Consider a shape functional J (Ω) whose shape derivative is expressed by J (cid:48) (Ω)( θ ) = (cid:82) ∂ Ω G θ.n . Usingthe discretization above, given v = ( a , a , ..., a N , b , ..., b N ) that defines Ω via the radial function ρ , afinite dimensional function is obtained j ( v ) = J (Ω). It is classical to compute the gradient of j using theshape derivative, by choosing the appropriate boundary perturbation for each Fourier coefficient. Usingthe notation r = x/ | x | we have r .n = ρ/ (cid:112) ρ + ( ρ (cid:48) ) . Therefore, denoting v n = ρ/ (cid:112) ρ + ( ρ (cid:48) ) , we obtain ∂j∂a k = (cid:90) ∂ Ω G cos( kt ) v n and ∂j∂b k = (cid:90) ∂ Ω G sin( kt ) v n . (19) Domains in R . In dimension three we choose to parametrize the unit sphere using ( φ, ψ ) ∈ [ − π, π ] × [0 , π ] (cid:55)→ (cos ψ cos φ, sin ψ cos φ, sin φ ). Next, we are interested in parametrizing radial functions ρ :[ − π, π ] × [0 , π ] which are constant for φ ∈ {− π, π } . This is needed in order to be able to create 3Dmeshes in FreeFEM [20] by deforming two dimensional meshes. One way of attaining this objective isto use two dimensional Fourier parametrizations which contain only sines for the φ coordinate, togetherwith an affine function in φ in order to allow different values at the extremities φ ∈ {− π, π } : ρ ( φ, ψ ) = aφ + b + N (cid:88) k =1 M (cid:88) l =1 ( c k,l sin(2 kφ ) cos( lψ ) + d k,l sin(2 kφ ) sin( lψ )) . As in dimension two, it is straightforward to infer the gradient of the discretized functional with respectto each one of the parameters. A simple computation yields v n = r .n = ρ/ (cid:113) ρ + ( ρ (cid:48) θ ) / cos φ + ( ρ (cid:48) φ ) : ∂j∂a = (cid:90) ∂ Ω G φv n , ∂j∂b = (cid:90) ∂ Ω G v n ,∂j∂c k,l = (cid:90) ∂ Ω G sin(2 kφ ) cos( lψ ) v n , ∂j∂d k,l = (cid:90) ∂ Ω G sin(2 kφ ) sin( lψ ) v n (20) Optimization algorithm.
Given the discretization and the gradients expressed above it is straight-forward to implement a gradient descent algorithm. The delicate issue is the fact that at each iteration,the objective function and its gradient are computed as a result of a minimization algorithm. If L ε (Ω , c )or P L ε (Ω , c ) are replaced by a local minimum instead of the global one, then this may give a false candi-date for the maximization algorithm. Therefore, we choose to work with a gradient flow type algorithm,which consists in advancing at each iteration in the direction given by the gradient of the functionalwith a prescribed step, regardless of the fact that the objective function increases or decreases. In thisway, even the optimization algorithms solved at one of the iterations yields a local minimum, the globaloptimization algorithm may still correct itself at subsequent iterations. The area constraint is imposedby a projection algorithm: the next iterate is rescaled to have the desired area via a homothety. Theprecise description is given in Algorithm 3.An example of a result obtained with this algorithm for the maximization of L ε (Ω , c ) is shown inFigure 10 together with the graph of the objective function. It can be seen that the objective function19 lgorithm 3 Global maximization algorithm
Require:
Initial Fourier coefficients, area constraints ( c for one phase, a vector c for the partitions), thenumber of iterations Niter , ε , initial step α , the number of iterations Nmod after which the step ishalved for i in { Niter } do Construct the mesh of Ω from the Fourier coefficients v : the size of triangles/tetrahedra shouldbe at most ε/ Approximate L ε (Ω , c ) (or P L ε (Ω , c ) in the case of partitions) Compute the gradient ∇ j ( v ): use (19) or (20) with G given by (17) (or (18) for the partitionscase) Advance in the direction of the gradient in order to increase the value of j ( v ): v ← v + α ∇ j ( v ) . Project on the area/volume constraint of Ω using a homothety If i mod Nmod ≡ α ← α/ end forreturn the final set of Fourier coefficients v C o s t Figure 10: Maximization of L (Ω , .
3) in dimension two together with the evolution of the cost function.increases and stabilizes as the size of the step decreases. Oscillations in the curve describing the cost havetwo main causes: first, the optimization algorithm at the current iteration might yield a local minimuminstead of the global one and secondly, the size of the step may be too big. An example for the caseof partitions is shown in Figure 11. Multiple instances of the gradient flow maximization algorithm arerepresented in Figure 12 for n = 6 and in Figure 13 for n = 10. Numerical aspects.
When minimizing F ε and P ε it is classical to consider meshes with elementsthat have size smaller than ε . This is due to the fact that the phase transition from 0 to 1 typically takesplace in a region of width proportional to ε and the mesh needs to be fine enough to capture this. Indimension two we consider ε = 0 .
05 giving rise to meshes having around 23 k nodes.In dimension three using ε = 0 . k nodes. When dealing with more cells indimension three we start with ε = 0 .
07 and we interpolate and re-optimize the result on a finer meshcorresponding to ε = 0 .
04. This gives meshes with around 35 k nodes. For postprocessing and plottingpurposes, the final mesh is further refined using MMG3D [12] such that more tetrahedra are presentwhere phases change quickly. The final partition is interpolated and re-optimized (with ε = 0 . k nodes) before plotting. Code.
The finite element software used for the optimization algorithm described in Section 3.1 isFreeFEM [20], which provides an interface to the LBFGS optimizer from Nlopt [23].The partition initialization via Voronoi diagrams is coded in Python, where optimization algorithmsfrom
Scipy.optimize and
Nlopt are used for unconstrained and, respectively, constrained optimizations.Codes and examples are provided in the following Github repository: https://github.com/bbogo/LongestShortestPartitions/tree/main/GradientVoronoi .The visualization is done with Python using Matplotlib in dimension two and Mayavi [31] in dimensionthree. The graphical representation of partitions is done by extracting surface meshes of an iso-level foreach cell in the optimal partition using FreeFEM [20] and MMG3D [12]. These surface meshes are then20
20 40 60 80 100 120Iterations1.8001.8251.8501.8751.9001.9251.950 C o s t Figure 11: Maximization of
P L (Ω , (1 / , / , / .
314 Iter 6: 3 .
381 Iter 13: 3 .
391 Iter 20: 3 .
436 Iter 70: 3 .
446 Iter 150: 3 . n = 6: the numerical optimalpartition and its associated cost are represented for a couple of iterations.plotted with Mayavi [31].Some codes used for obtaining the results illustrated in the paper can be found in the Github reposi-tory: https://github.com/bbogo/LongestShortestPartitions/tree/main/FreeFEMcodes . In this section we use the algorithm described previously in order to study problems (2) and (4).Results from [15] show that problem (2) is solved by the disk in dimension two for c = 1 /
2. We performsimulations for various values of c < / c or 1 − c for the constraint gives the sameresult) and the numerical result is the disk in dimension two. In dimension three the same phenomenonoccurs: for various values of the volume fraction c the shape which maximized the relative minimalperimeter of a subset with volume c | Ω | is the ball. Some examples are shown in Figure 14.Surprisingly, the case of partitions shows similar results. When considering equal area constraints theset with fixed area maximizing the length of the minimal partition is still the disk (see Figure 15 for someexamples). In dimension three for n ∈ { , , , } we obtain similar results: the ball maximizes the totalsurface area of the smallest total perimeter partition. These results are illustrated in Figures 17 and 18.Note that simulations made in the one phase case already shows that when partitioning the domaininto two regions with two non-equal areas, the maximizer of the minimal length partition is still thedisk. This suggests that even in the case where the cells do not have the same prescribed area, the setΩ which maximizes the minimal perimeter of a partition is still the disk in 2D (the ball in 3D). Indeed,Iter 1: 4 .
663 Iter 6: 4 .
795 Iter 13: 4 .
794 Iter 20: 4 .
863 Iter 70: 4 .
892 Iter 150: 4 . n = 10: the numerical optimalpartition and its associated cost are represented for a couple of iterations.21igure 14: Maximization of the minimal relative perimeter in 2D and 3D with volume constraints c ∈{ . , . , . } . The optimal set Ω (the disk/ball) together with the set obtained numerically whenminimizing the relative perimeter for the given volume fraction.Figure 15: Maximization of the length of the minimal perimeter partition into equal areas for n ∈ { , , } Figure 16: Maximization of the length of the minimal perimeter partition into different areas: n = 3,ratios 1 : 2 : 3, n = 4, ratios 1 : 1 : 1 . . n = 4 : 1 : 1 : 1 . . n ∈ { , } . (right) Results obtained when the area constraints are not the same: n = 3: ratios 1 : 2 : 2, n = 4: ratios 1 : 2 : 2 : 2.Figure 18: Maximization of the length of the minimal perimeter partition into equal areas for n ∈ { , } .An expanded view of the optimal partition is also illustrated for each case.when considering more cells with different areas, the numerical result is the same: the disk seems to bethe maximizer (see Figure 16 for some examples). As already underlined in [19], the study of partitionsin cells with prescribed but different areas is more complex, since in this case there are even more localminima.The numerical simulations give rise to the following conjecture, which generalizes the results of [15]. Conjecture 4.1.
1. Given c ∈ (0 , , the set Ω maximizing L (Ω , c ) under the constraint | Ω | = v d (i.e.solving (2) ) is the ball.2. Given n > and c = ( c i ) ni =1 ∈ R n + with (cid:80) ni =1 c i = 1 , the set Ω maximizing P L (Ω , c ) under theconstraint | Ω | = v d (i.e. solving (4) ) is the ball. Remark 4.2.
The same type of results seem to hold when maximizing the minimal geodesic perimeterfor closed surfaces in Ω in R which are boundaries of convex sets with a constraint on the H ( ∂ Ω) . Thetechniques used in this case are those presented in [5] and the theoretical and numerical framework issimilar to what was done in dimension three. In this case the sphere seems to be the maximizing set,which is in accord with the conjecture stated above. As discussed in Section 3.4, existence of shape derivatives for L ε and P L ε depends on the uniquenessof the minimizers for these functionals. Therefore, it is not straightforward to obtain classical optimalityconditions. It is possible, however, to obtain some qualitative information about sets maximizing theminimal values of L ε and P L ε under volume constraint. Recall that the optimizer of shape differentiablefunctional J under volume constraint will verify an optimality condition of the form J (cid:48) (Ω)( θ ) + (cid:96) | Ω | (cid:48) ( θ ) = 0 , (21)where (cid:96) ∈ R is a Lagrange multiplier associated to the volume constraint. Recall that the shape derivativeof the volume functional is | Ω | (cid:48) ( θ ) = (cid:82) ∂ Ω θ.n . Non-uniqueness of the minimal relative perimeter set/partition at the optimum.
Resultsin Section 3.4 show that the shape derivatives of L ε and P L ε exist when they correspond to uniqueminimizers of the Modica-Mortola type functionals. In this case, the corresponding shape derivatives areboundary integrals of non-constant functions multiplied by the normal perturbation θ.n . Therefore, it is23traightforward to see that a relation of the type (21) cannot hold. This allows us to conclude that of Ω ∗ is a minimizer of L ε (Ω , c ) the optimal minimal relative perimeter set is not unique. The same happens inthe case of partitions: if Ω ∗ minimizes P L ε (Ω , c ) and the minimal length partition of Ω with constraints c is not unique.Suppose now that Ω is a domain with fixed volume | Ω | = v d such that SF (Ω , c ) is unique. Thenfor ε > L ε (Ω , c ) will also be unique. Thus L ε (Ω , c ) admits a shape derivative. Also thecorresponding optimal density u Ω is not constant on the boundary, and therefore the optimality relation(21) cannot hold. This shows that such a domain Ω is not a solution of problem (2). The same argumentcan be applied for problem (4).As a conclusion, a domain Ω that solves (2) must have multiple minimal relative perimeter sets (ormultiple minimal length partitions for problem (4)). Infinitely many minimal relative perimeter sets/partitions at the optimum.
Suppose thatΩ is a domain with volume | Ω | = v d such that there are a finite number of minimizers SF (Ω , c ). Thenthere are relatively open subsets Γ ⊂ ∂ Ω such that Γ does not intersect the boundary of any of theseminimizers. It is possible to perturb ∂ Ω inwards on Γ without creating additional solutions SF (Ω , c ).This process diminishes the volume of | Ω | . Rescaling the perturbed domain to have the initial volumeincreases the length of minimal relative perimeter set SF (Ω , c ).The same argument can be repeated in the case where there are subsets of ∂ Ω which do not intersectany of the minimizers SF (Ω , c ). The same type of arguments applies also the case of partitions SP (Ω , c ).Therefore any set Ω which solves (2) or (4) should verify the following observations: • there are infinitely many minimal relative perimeter sets (partitions) SF (Ω , c ) ( SP (Ω , c )) • for any relatively open set Γ ⊂ ∂ Ω there exist a minimal relative perimeter set (partition) SF (Ω , c )( SP (Ω , c )) which intersects ΓThese observations show that sets with radial symmetry are natural candidates to be solutions of problems(2), (4). This fact is also confirmed by the numerical simulations shown in Section 4. The theoretical considerations and numerical simulations presented in this paper suggest that theresults of [15] are valid in more general settings: under volume constraint the ball is the set Ω whomaximizes • the minimal relative perimeter of a subset ω ⊂ Ω with volume constraint | ω | = c | Ω | for all c ∈ (0 , • the minimal relative perimeter of a partition of Ω into sets ( ω i ) ni =1 with volume constraints | ω i | = c i | Ω | given c i ∈ (0 ,
1) with (cid:80) ni =1 c i = 1. The result seems to hold even in the case where the sets | ω i | do not have the same volume constraints.The numerical maximization algorithm consisted in solving at each iteration an optimization problemwhich approximates the least perimeter set or partition under the given constraint. Then a perturbationof the set which does not decrease the minimal perimeter is found and the set is modified. In all cases,the numerical result was close to the disk/ball.The initialization phase for the computation of the optimal partitions is made using Voronoi diagramswith prescribed capacity. We provide a new way of generating such Voronoi diagrams using the gradientsof the areas with respect to the Voronoi points. The gradient of the perimeter of the Voronoi cells is alsocomputed. Acknowledgments
The authors were partially supported by the project ANR-18-CE40-0013 SHAPO financed by theFrench Agence Nationale de la Recherche (ANR).
References [1] G. Alberti. Variational models for phase transitions, an approach via Γ-convergence. In
Calculus ofvariations and partial differential equations (Pisa, 1996) , pages 95–114. Springer, Berlin, 2000.[2] L. Ambrosio, N. Fusco, and D. Pallara.
Functions of bounded variation and free discontinuity prob-lems . Oxford Mathematical Monographs. The Clarendon Press, Oxford University Press, New York,2000. 243] M. Balzer. Capacity-constrained voronoi diagrams in continuous spaces. In , pages 79–88, 2009.[4] M. Balzer, T. Schl¨omer, and O. Deussen. Capacity-constrained point distributions: A variant oflloyd’s method. In
ACM SIGGRAPH 2009 Papers , SIGGRAPH ’09, New York, NY, USA, 2009.Association for Computing Machinery.[5] B. Bogosel and E. Oudet. Partitions of minimal length on manifolds.
Exp. Math. , 26(4):496–508,2017.[6] A. Braides.
Approximation of Free-Discontinuity Problems . Springer, 1998.[7] E. Bretin, R. Denis, J.-O. Lachaud, and E. Oudet. Phase-field modelling and computing for a largenumber of phases.
ESAIM Math. Model. Numer. Anal. , 53(3):805–832, 2019.[8] V. Burenkov. Extension theorems for Sobolev spaces. In
The Maz’ya anniversary collection, Vol. 1(Rostock, 1998) , volume 109 of
Oper. Theory Adv. Appl. , pages 187–200. Birkh¨auser, Basel, 1999.[9] G. Buttazzo. Gamma-convergence and its Applications to Some Problems in the Calculus of Varia-tions.
School on Homogenization ICTP, Trieste, September 6-17 , 1993.[10] A. Cianchi. On relative isoperimetric inequalities in the plane.
Boll. Un. Mat. Ital. B (7) , 3(2):289–325, 1989.[11] S. J. Cox and E. Flikkema. The minimal perimeter for N confined deformable bubbles of equal area. Electron. J. Combin. , 17(1):Research Paper 45, 23, 2010.[12] C. Dapogny, C. Dobrzynski, and P. Frey. Three-dimensional adaptive domain remeshing, implicit do-main meshing, and applications to free and moving boundary problems.
J. Comput. Phys. , 262:358–378, 2014.[13] M. C. Delfour and J.-P. Zol´esio.
Shapes and geometries , volume 22 of
Advances in Design andControl . Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second edition,2011. Metrics, analysis, differential calculus, and optimization.[14] Q. Du, V. Faber, and M. Gunzburger. Centroidal voronoi tessellations: Applications and algorithms.
SIAM Review , 41(4):637–676, Jan. 1999.[15] L. Esposito, V. Ferone, B. Kawohl, C. Nitsch, and C. Trombetti. The longest shortest fence andsharp Poincar´e-Sobolev inequalities.
Arch. Ration. Mech. Anal. , 206(3):821–851, 2012.[16] W. J. Firey. Lower bounds for volumes of convex bodies.
Archiv der Mathematik , 16(1):69–74, Dec.1965.[17] D. Gilbarg and N. S. Trudinger.
Elliptic partial differential equations of second order . Classics inMathematics. Springer-Verlag, Berlin, 2001. Reprint of the 1998 edition.[18] M. E. Gurtin and H. Matano. On the structure of equilibrium phase transitions within the gradienttheory of fluids.
Quart. Appl. Math. , 46(2):301–317, 1988.[19] F. Headley and S. Cox. Least perimeter partition of the disc into n bubbles of two different areas.
The European Physical Journal E , 42(7), July 2019.[20] F. Hecht. New development in FreeFem++.
J. Numer. Math. , 20(3-4):251–265, 2012.[21] A. Henrot and M. Pierre.
Shape variation and optimization , volume 28 of
EMS Tracts in Mathemat-ics . European Mathematical Society (EMS), Z¨urich, 2018. A geometrical analysis, English versionof the French publication [ MR2512810] with additions and updates.[22] M. Iri, K. Murota, and T. Ohya. A fast voronoi-diagram algorithm with applications to geographicaloptimization problems. In
System Modelling and Optimization , pages 273–288. Springer-Verlag.[23] S. G. Johnson. The nlopt nonlinear-optimization package. http://github.com/stevengj/nlopt .2524] B. Kloeckner. Dans quelle forme la plus petite paroi enfermant un volume donn´e est-elle la plus grande ? - Images des math´ematiques : http://images.math.cnrs.fr/Dans-quelle-forme-la-plus-petite-paroi-enfermant-un-volume-donne-est-elle-la.html ,12 2019.[25] A. Koldobsky, C. Saroglou, and A. Zvavitch. Estimating volume and surface area of a convex bodyvia its projections or sections.
Studia Math. , 244(3):245–264, 2019.[26] F. Maggi.
Sets of finite perimeter and geometric variational problems , volume 135 of
CambridgeStudies in Advanced Mathematics . Cambridge University Press, Cambridge, 2012. An introductionto geometric measure theory.[27] L. Modica. The gradient theory of phase transitions and the minimal interface criterion.
Arch.Rational Mech. Anal. , 98(2):123–142, 1987.[28] L. Modica and S. Mortola. Un esempio di Γ − -convergenza. Boll. Un. Mat. Ital. B (5) , 14(1):285–299,1977.[29] F. Morgan. Soap bubbles in R and in surfaces. Pacific J. Math. , 165(2):347–361, 1994.[30] ´E. Oudet. Approximation of partitions of least perimeter by Γ-convergence: around Kelvin’s con-jecture.
Exp. Math. , 20(3):260–270, 2011.[31] P. Ramachandran and G. Varoquaux. Mayavi: 3d visualization of scientific data.
Computing inScience Engineering , 13(2):40–51, 2011.[32] M. Ritor´e and E. Vernadakis. Isoperimetric inequalities in Euclidean convex bodies.
Trans. Amer.Math. Soc. , 367(7):4983–5014, 2015.[33] K. Svanberg. A class of globally convergent optimization methods based on conservative convexseparable approximations.
SIAM Journal on Optimization , 12(2):555–573, Jan. 2002.[34] S.-Q. Xin, B. L´evy, Z. Chen, L. Chu, Y. Yu, C. Tu, and W. Wang. Centroidal power diagrams withcapacity constraints: Computation, applications, and extension.