Partitions of Minimal Length on Manifolds
aa r X i v : . [ m a t h . O C ] J un PARTITIONS OF MINIMAL LENGTH ON MANIFOLDS
BENIAMIN BOGOSEL AND ´EDOUARD OUDET
Abstract.
We study partitions on three dimensional manifolds which minimize the total geodesicperimeter. We propose a relaxed framework based on a Γ-convergence result and we show somenumerical results. We compare our results to those already present in the literature in the case of thesphere. For general surfaces we provide an optimization algorithm on meshes which can give a goodapproximation of the optimal cost, starting from the results obtained using the relaxed formulation. Introduction
In this article we propose a theoretical and numerical framework for the study of the partitions( ω i ) ni =1 of a surface S ⊂ R which minimize the total geodesic perimeter while keeping a prescribedarea for each cell. Thus, we are interested in minimizing H ( ∪ ni =1 ∂ S ω i ) or equivalentlyPer( ω ) + ... + Per( ω n )in the class of partitions ( ω i ) of the surface S such that | ω i | = c i , with the compatibility constraint c + ... + c n = | S | . Here ∂ S ω denotes the boundary of a set ω as a subset of the surface S , Per( ω )denotes the geodesic perimeter of ω , i.e. the perimeter of ω regarded as a subset of the surface S and | ω | is the area of the subset ω . General theoretical results concerning these minimal partitioningproblems are presented by Morgan in [16]. This theoretical result states that boundaries of aminimal-perimeter partition are arcs of constant geodesic curvature and the boundaries of the setsmeet in threes with angles of measure 2 π/ n = 2the solution is the partition into two half-spheres. This was proved by Bernstein in 1905 [5]. Inthe case n = 3 the optimal candidate is the partition of the sphere into three slices correspondingto an angle of 2 π/
3. This was proved by Masters in [15]. The case n = 12 was solved by Hales in[13] using methods similar to the ones involved in the proof of the honeycomb conjecture [12]. Thecase n = 4 was treated by Engelstein in [11] and the corresponding optimal partition is the oneassociated to the regular tetrahedron.The case of the sphere has been studied numerically by Cox and Flikkema [9] using the SurfaceEvolver software [7]. They perform computations for n ∈ J , K and they confirm the naturalconjecture for n = 6: the optimal partition in this case is probably the one associated to the cube.Their algorithm performs the perimeter optimization after choosing a topological structure for thepartition. Thus, the optimization algorithm has to know a priori the topological structure in orderto find the corresponding local minimum. In the end we keep the configuration which gives the bestoptimal cost among the admissible combinatorial possibilities.The algorithm we propose is a generalization of the ideas in [17] to the case of surfaces. First,there is a theoretical result, similar to the theorem of Modica and Mortola, which we present inSection 2. This theoretical result justifies the use of the functional J ε ( u ) = ε Z S |∇ τ u | + 1 ε Z S u (1 − u ) as an approximation of the perimeter as ε →
0. The direct consequence of the Γ-convergence resultis that a sequence of minimizers u ε for J ε under the constraint R S u ε = c converges to a minimizer of the geodesic perimeter under area constraint. For the partitioning case we prove that functionalsof the type n X i =1 J ε ( u i )approximate the perimeter as ε →
0, where u i are functions associated to the sets ω i which satisfysome integral and non-overlapping constraints. We implement an optimization algorithm which isable to solve the above problem on a large class of surfaces. This is an advantage over the methodsused in [9] which can be used only in the case of the sphere.Working with the relaxed formulation does not provide an exact representation of the contours.Thus, we cannot directly provide the associated cost once we have the relaxed optimal partitions.The particular case of the sphere can be solved directly by noting that boundaries between twocells have constant geodesic curvature [16] and are, thus, arcs of circles. We recover all the resultspresented in [9] in the case of the sphere. On more complex surfaces it is complicated to explicitlywork with curves of constant geodesic curvature. Nevertheless, we can extract the contours fromthe density representation in order to compute the total perimeter. Since the extracted contoursare not smooth, we perform a constrained optimization stage on the triangulated surface preservingthe topology to obtain reliable approximations of the optimal costs.2. Theoretical result
As in [17] we would like to have a rigorous theoretical framework which justifies our numericalmethod. In the euclidean case it was an adapted version of the Modica-Mortola theorem to thecase of partitions which provided the needed result. In the case of surfaces we did not find anequivalent result in the literature. We did find the results in [4] which suggest that the relaxationwe consider is the right one on general manifolds. In the above reference a the authors do not provea Γ-convergence result, but only the convergence of minimisers. We are concerned here only withsmooth manifolds of codimension one and in this particular case it is possible to adapt classicalmethods in order to prove a Γ-convergence result.We start by defining the space of functions of bounded variations on a d − R d . Let S be a smooth d − R d . In the following weconsider the tangential gradient of a function u defined on S to be ∇ τ u = ∇ ˜ u − ( ∇ ˜ u.n ) n, where ˜ u is a regular extension of u in a neighbourhood of S n denotes the normal vector to thesurface. In the same way we define the tangential divergence of a vector field w ∈ C ( S ; R d ) bydiv τ w = tr( D τ w )where the matrix D τ w contains on line i the tangential gradient of the i -th component of w , i.e. ∇ τ w i . See [14, Section 5.4] for further details.We consider the space of functions with bounded variation on SBV ( S ) = { u ∈ L ( S ) : T V ( u ) < ∞} where T V ( u ) = sup { Z S u div τ g : | g | ∞ ≤ } . Using the divergence theorem on manifolds (see [14, Section 5.4]), we obtain that if u is C ( S ) then T V ( u ) = Z S |∇ τ u | . If ω is a subset of S we define its generalized perimeter as Per( ω ) = T V ( χ ω ), where χ ω representsthe characteristic function of ω . By mimicking the proof in the euclidean case we can prove that thetotal variation is lower semi-continuous for the L ( S ) convergence. We refer to [6] for more details. ARTITIONS OF MINIMAL LENGTH ON MANIFOLDS 3
Let ( C i ) is a set of local charts which cover S such that each C i is diffeomorphic to a connectedand bounded open subset D i of R d − . We denote by θ i : D i → C i these diffeomorphisms. Then itis possible to transfer a function u from C i to D i using the transformation ˜ u i = u ◦ θ i . These newfunctions ˜ u i , which lie now in Euclidean spaces, are functions of bounded variation. Therefore, it ispossible to transfer some of the theory of BV functions from Euclidean spaces to manifolds of co-dimension 1 by using local charts and partitions of unity. In particular, it is possible to approximatefinite perimeter sets ω ⊂ S with smooth sets ω n ⊂ S such that ω n → ω in the L ( S ) topology andPer( ω n ) → Per( ω ).We are now ready to state the relaxation result in the case of a single phase, which will begeneralized later to the case of a partition. To derive the theorem below we follow the approachprovided by Buttazzo in [8] and Alberti in [1]. Theorem 2.1.
Define F ε , F : L ( S ) → [0 , + ∞ ] as follows: F ε ( u ) = Z S (cid:18) ε |∇ τ u | + 1 ε u (1 − u ) (cid:19) dσ if u ∈ H ( S ) , R S u = c + ∞ otherwise. F ( u ) = ( Per( { u = 1 } ) if u ∈ BV ( S, { , } ) , R S u = c + ∞ otherwise.Then F ε Γ −→ F in the L ( S ) topology.Proof: We define φ ( t ) = R t | s (1 − s ) | ds . We consider a sequence ( u ε ) → u in L ( S ) such thatlim inf ε → F ε ( u ε ) < + ∞ . Since F ε ( u ε ) ≥ ε R S u ε (1 − u ε ) , if we take a subsequence of u ε whichconverges almost everywhere to u we obtain that Z S u (1 − u ) = 0 , and thus u ∈ { , } almost everywhere in S . Note that truncating u ε between 0 and 1 decreasesthe value of F ε ( u ε ) while preserving the fact that u ε → u in L ( S ). Also note that φ is Lipschitzon [0 ,
1] so we can conclude that φ ◦ u ε → φ ◦ u in L ( S ). By applying the classical inequality a + b ≥ ab we get that F ε ( u ε ) ≥ Z S |∇ τ u | φ ′ ( u ε ) = 2 Z S |∇ τ ( φ ◦ u ε ) | . Taking lim inf in the above inequality and using the semi-continuity of the total variation withrespect to the L ( S ) convergence we obtain thatlim inf ε → F ε ( u ε ) ≥ T V ( φ ◦ u ) = 2 φ (1) T V ( u ) . Since u is a characteristic function, it follows that the perimeter of { u = 1 } is bounded and therefore u ∈ BV ( S, { , } ). Note that φ (1) = 1 / − lim inf part of the theorem.For the Γ − lim sup part we need to exhibit a recovery sequence for each u such that F ( u ) < + ∞ .By a classical argument it is enough to find a recovery sequence only for functions u which arecharacteristic functions of smooth sets in S . See [6] for more details concerning the reduction toregular sets and [3, Theorem 3.42] for the BV approximation of finite perimeter sets with smoothsets.Let’s consider now u = χ ω where ω ⊂ S is a set with smooth boundary relative to S . We considerthe signed distance function d ω : S → R defined by d ω ( x ) = d τ ( x, S \ ω ) − d τ ( x, ω ) , B. BOGOSEL AND ´E. OUDET where d τ is the geodesic distance on S . Note that d ω is positive outside ω and negative inside.Consider the optimal profile problem c = min { Z R ( W ( v ) + | v ′ | : v ( −∞ ) = 0 , v (+ ∞ ) = 1 } . Any solution of this minimizing problem satisfies v ′ = p W ( v ) and we can impose the initialcondition v (0) = 1 / c = 2 R p W ( s ) ds . In our problem we have chosen W ( s ) = s (1 − s ) . In order to have a functionwhich goes from 0 to 1 in finite time we may choose v η = min { max { , (1 + 2 η ) v − η } , } . We see that as η → c η = Z R ( W ( v η ) + | ( v η ) ′ | ) → c as η → . All these considerations are inspired from [6]. We can define u ε ( x ) = v η ( d ω ( x ) /ε ) . We can see that F ε ( u ε ) = Z S (cid:18) ε |∇ τ u | + 1 ε W ( u ) (cid:19) = Z T/ε − T ε Z d ω ( x )= t (cid:18) ε | ( v η ) ′ ( d ω ( x ) /ε ) | |∇ τ d ω ( x ) | ε + 1 ε W ( v η ( d ω /ε )) (cid:19) d H d − ( x ) dt = Z T/ε − T/ε Z d ω ( x )= t ε ( | ( v η ) ′ ( t/ε ) | + W ( v η ( t/ε ))) d H d − ( x ) dt = Z T/ε − T/ε
Per( d ω ( x ) = t ) 1 ε ( | ( v η ) ′ ( t/ε ) | + W ( v η ( t/ε ))) dt = Z T − T Per( d ω ( x ) = tε )( | ( v η ) ′ ( t ) | + W ( v η ( t ))) dt where we have applied the co-area formula and T is chosen such that the support of v η is inside[ − T, T ]. Since lim s → Per( { d ω ( x ) = s } ) = Per( ω ) we see that for ε small enough there exists δ suchthat Per( d ω ( x ) = s ) < Per( ω ) + δ when | s | < T ε . Thereforelim sup ε → F ε ( u ε ) ≤ (Per( ω ) + δ ) Z T − T ( | ( v η ) ′ ( t ) | + W ( v η ( t ))) dt = (Per( ω ) + δ ) c η . Since this is true for any δ, η small enough, by letting δ, η → R S χ ω = c it is enough to consider a shift in the definitionof u ε : u ε ( x ) = v η (( d ω ( x ) + s ε ) /ε ) , where s ε ∈ [ − T ε, T ε ] . We can see that for s ε = T ε we have u ε = 1 on ω and thus R S u ε > c whilefor s ε = − T ε the support of u ε is included in ω and we have the opposite inequality. Thus, for each ε small enough we can change the definition of u ε so that R S u ε = c . The estimates presented aboveare carried with no difficulty in this setting. (cid:3) We can now state the result in the partitioning case. We denote by u an element in ( L ( S )) n . Inorder to simplify the notations we introduce the space X = { u ∈ ( L ( S )) n : Z S u i = c i , n X i =1 u i = 1 } ARTITIONS OF MINIMAL LENGTH ON MANIFOLDS 5 where c i satisfy the compatibility condition P ni =1 c i = H d − ( S ). It is easy to see that X is closedunder the convergence in ( L ( S )) n . Theorem 2.2.
Define F ε , F : ( L ( S ) n → [0 , + ∞ ] as follows: F ε ( u ) = n X i =1 Z S (cid:18) ε |∇ τ u i | + 1 ε u i (1 − u i ) (cid:19) dσ if u ∈ ( H ( S )) n ∩ X + ∞ otherwise F ( u ) = ( P ni =1 Per( { u i = 1 } ) if u ∈ ( BV ( S, { , } )) n ∩ X + ∞ otherwiseThen F ε Γ −→ F in the ( L ( S )) n topology.Proof: It is easy to see that the Γ − lim inf part follows at once from Theorem 2.1 and from thefact that X is closed under the topology of ( L ( S )) n .In order to construct the recovery sequence we reduce the problem to the case where the limit u is consists of piecewise smooth parts in S . In this case we define u i = v η ( d ω i ( x ) /ε ) as in the onephase case. Thus on each ω i we have u i ≥ / P ni =1 u i ≥ /
2. There are twopoints which need to be addressed:(1) The sum equal to 1 condition. Due to the symmetry of the optimal profile we deduce thatthere is only one zone where the sum condition is not satisfied and that is in the neighborhoodof singular points. Since an ε -neighborhood of the singular set is of order ε d − . Replacingeach u i by u i / ( P ni =1 u i ) in these problematic regions we preserve the regularity of each u i and we note that the functions have bounded gradient of order O (1 /ε ). We immediatelyfind that the corresponding energy Z N ε (cid:18) ε |∇ τ u i | + 1 ε u i (1 − u i ) (cid:19) vanishes as ε → u i so that they have the same integral over S . In orderto do this we apply a procedure found in [2] where we consider a family of balls in regionswhere u i ∈ { , } . On each such ball we can consider modifications of u i such that the sumis preserved and the integrals have the right value. As above, the sum of energies on theseballs will be negligible in the limit.Once these points are addressed, the lim sup estimates follows just like in the one dimensional caseand the proof of the theorem is completed. (cid:3) Finite Element framework
We wish to use this relaxation by Γ-convergence to perform numerical computations so we needa framework which allows us to compute the quantity ε Z S |∇ τ u | + 1 ε Z S u (1 − u ) , in fast, efficient way. In order to do this we triangulate the surface S and we compute the massmatrix M and the stiffness matrix K associated to the P finite elements on this triangulation. Then,if for the sake of simplicity, we use the same notation u for the P finite element approximation of u , we have Z S |∇ τ u | = u T Ku B. BOGOSEL AND ´E. OUDET and Z S u (1 − u ) = v T M v, where v = u. . × (1 − u ) . . We have used the Matlab convention that adding a point beforean operation means that we are doing component-wise vector computations. Note that once thematrices K, M are computed, we only have to perform matrix-vector multiplications, which is reallyfast. In this setting we use the discrete gradients of the above expressions given by: ∇ u u T Ku = 2 Ku, ∇ u v T M v = 2
M v. × (1 − u ) . The partition condition and the equal areas constraint are imposed by making an orthogonalprojection on the linear constraints as follows. We write the discrete vectors representing P dis-cretization of the density functions in the following matrix form M = ( ϕ ϕ ... ϕ n ) . The partition constraint implies that the sum of the elements on every line of M is equal to 1 andthe equal area constraint implies that for every column of the matrix M we have the relation h v, ϕ i i = A/n, where v = × N · M. Here the constant A is the total area of the surface, N is the total number of points in the trian-gulation and the notation p × q represents the p × q matrix whose entries are all equal to 1. Theseconditions are discretizations in the finite element setting of the conditions that the integrals of thedensity functions u i are all equal to A/n . Indeed, given a triangulation T of S and its associatedmass matrix M , we have Z S · u i = × N · M · ϕ i , where ϕ i is the vector containing the values of u i at the vertices of the triangulation. The projection routine can be found in Algorithm 1. Algorithm 1
Orthogonal projection on the partition and area constraints
Require: A = ( a ij ) ∈ R N × n , c ∈ R × n , d ∈ R N × , v ( e i ) = P j a ij − c i (line sum error; N × ( f i ) = P i v i a ij − d j (column scalar product error; n × Define the matrix C of size n × n by ( c kl = k v k /n if k = lc kk = k v k − k v k /n ( q j ) = ( f j ) − h v, e i /n ( n × Compute ( λ j ) ∈ R n × with λ n = 0 such that C | ( n − × ( n − ( λ j ) | n − = ( q j ) | n − . The indicesindicate a sub-matrix with the first n − n − S = P j λ j η i = ( e i − S · v i ) /n ( N × A orth = ( η i ) · × n + v · ( λ j ) T , where p × q is the p × q matrix with all entries equal to 1 A = A − A orth return A Once we have this discrete formulation we use an optimized LBFGS gradient descent procedure[19] to compute the numerical minimizers. In order to avoid local minima where one of the phases
ARTITIONS OF MINIMAL LENGTH ON MANIFOLDS 7 ϕ l is constant, which arise often when the number of phases is greater than 5, we add a Lagrangemultiplier which penalizes the constant functions. In this way, we optimize n X i =1 ε Z S |∇ τ ϕ i | + 1 ε Z S ( ϕ i ) (1 − ϕ i ) + λ (std( ϕ i ) − starget) , where std( ϕ l ) is the standard deviation of ϕ l and starget is the standard deviation of a characteristicfunction of area Area( S ) /n .In order to have a good approximation of the optimal partition, we want do decrease ε so thatthe width of the interface is small. We notice that if we chose ε of the same order as the sides ofthe mesh triangles the algorithm converges. Furthermore, we cannot make ε smaller, since then thegradient term will not contain any real information, as the width of the interface is of size ε . Inorder to avoid this problem, we consider refined meshes associated to each ε . At each step wherewe decrease ε we interpolate the values of the previous optimizer on a refined mesh and we considerthese interpolated densities as starting point for the descent algorithm on the new mesh. In thecase of the sphere we make four refinements ranging from 10000 to 160000 points. Some optimalconfigurations, in the case of the sphere, are presented in Figure 1. A detailed study of the case ofthe sphere along with a comparison with the known results of Cox and Flikkema [9] are presentedin the next section.As underlined before, our approach allows a direct treatment of any surface, as long as a qual-itative triangulation is found. We perform some numerical computations on various shapes likea torus, a double torus, and a more complex surface called Banchoff-Chmutov of order 4. A fewdetails about the definitions of these surfaces are provided below: • We consider a torus of outer radius R = 1 and inner radius 0 . f ( x, y, z ) = ( x + y + z + R − r ) − R ( x + y ) . • The double torus used in the computation (see Figure 3 is given by the zero level set of thefunction f ( x, y, z ) = ( x ( x − ( x −
2) + y ) + z − . . • The complex Banchoff-Chmutov surface (see Figure 4) is given by the zero level set of thefunction f ( x, y, z ) = T ( x ) + T ( y ) + T ( z ) , where T ( X ) = 8 X − X + 1 is the Tchebychev polynomial of order 4.4. Refined optimization in the case of the sphere
The costs associated to the relaxed functional do not provide a good enough approximation ofthe total length of the boundaries. In this section we propose a method to approximate the optimalcost in the case of the sphere. The results of [16] state that boundaries of the cells of the optimalpartitions have constant geodesic curvature. In the case of the sphere the only such curves are thearcs of circle. See for example [18, Exercise 2.4.9] for a proof. The results of Cox and Flikkema[9] show that optimal configurations are not made of geodesic polygons. In order to perform anoptimization procedure which captures this effect they chose to make an initial optimization inthe class of geodesic polygons and then divide each geodesic arc into 16 smaller arcs and restartthe procedure with more variable points. They manage to approximate well enough the generaloptimal structure but they still work in the class of geodesic polygons with additional vertices.Our approach presented below is different in the sense that we consider general circle arcs (notnecessarily geodesics) which connect the points.The first step is to extract the topology of the partition from the previous density results, i.e.locate the triple points, the edge connections and construct the faces. In order to perform the
B. BOGOSEL AND ´E. OUDET
Figure 1.
Minimal perimeter partitions on the sphere into n equal area cells for n ∈ { , , ..., , } . ARTITIONS OF MINIMAL LENGTH ON MANIFOLDS 9
Figure 2.
Minimal perimeter partitions on the torus with outer radius R = 1 andinner radius r = 0 . n ∈ [2 , Figure 3.
Minimal perimeter partitions on a double torus for n ∈ { , , } . Figure 4.
Minimal perimeter partitions on a Banchoff-Chmutov surface for n ∈ { , , , } . refined optimization procedure we need to be able to compute the areas of portions of the spheredetermined by arcs of circles. This is possible using the Gauss-Bonnet formula. If M is a smoothsubset of a surface then Z M KdA + Z ∂M k g = 2 πχ ( M ) , (4.1)where K is the curvature of the surface, k g is the geodesic curvature and χ ( M ) is the Euler char-acteristic of M . This result extends to piecewise smooth curves and in this case we have Z M KdA + Z ∂M k g + X θ i = 2 πχ ( M ) , (4.2)where θ i are the turning angles between two consecutive smooth parts of the boundary. In the caseof a polygon the turning angles are the external angles of the polygon. The formula (4.2) allows thecomputation of the area of a piece of the sphere bounded by arcs of circle. In this case the Eulercharacteristic is equal to 1, the curvature of the unit sphere is K = 1 and the geodesic curvature ispiecewise constant. For more details we refer to [10, Chapter 4].A first consequence of the Gauss-Bonnet theorem in connection to our problem is noting the factthat, apart from cases where we have a certain symmetry like n ∈ { , , , } the optimal cellsare not geodesic polygons. This is made clear in cases where we have a hexagonal cell. If the arcsforming the boundary of such a hexagonal cell would be geodesic polygons then its area would beequal to 6 · π/ − π = 0. Thus a spherical shape bounded by six arcs of circle can never be ageodesic polygon without being degenerate.In order to perform the optimization we take the vertices as variables and we add one supple-mentary vertex for each edge. This is enough to contain all the necessary information since an arcof circle is well defined by three distinct points on the sphere. In the sequel we denote P n the set ofpartitions of the sphere into n cells and with A n the partitions in P n having equal areas. In orderto have a simpler numerical treatment of the problem we can incorporate the area constraints inthe functional by defining for every partition ( ω i ) ∈ P n the quantity defined for every ε > G ε (( ω i )) = n X i =1 Per( ω i ) + 1 ε n − X i =1 n X j = i +1 (Area( ω i ) − Area( ω j )) . If we denote G (( ω i )) = (P ni =1 Per( ω i ) if ( ω i ) ∈ A n ∞ if ( ω i ) ∈ P n \ A n . then we have the following Γ-convergence result. Theorem 4.1.
We have that G ε Γ −→ G for the L ( S ) convergence of sets.Proof: For the (LI) property consider a sequence ( ω εi ) ⊂ P n which convergence in L ( S ) to ( ω i ).It is clear that we have Area( ω εi ) → Area( ω i ) and the perimeter is lower semicontinuous for the L convergence. Thus we have two situations. If ( ω i ) ∈ P n \A n then lim ε → G ε (( u εi )) = ∞ . If ( ω i ) ∈ A n then the lower semicontinuity of the perimeter implies that lim inf ε → G ε (( ω εi )) ≥ G (( ω i )).The (LS) property is immediate in this case. Choose ( ω i ) ∈ A n , or else there is nothing to prove.We may choose the recovery sequence equal to ( ω i ) for every ε >
0. Thus the property is verifiedimmediately. (cid:3)
Remark 4.2.
We note that in the above proof the simplicity of the proof of the (LS) property isdue to the fact that the functionals G ε are well defined on the space { G < ∞} , which makes possiblethe choice of constant recovery sequences. This is not the case in the results proved in Section 2.This Γ-convergence result proves that minimizers of G ε converge to minimizers of G . As aconsequence, in the numerical computations, we minimize G ε for ε smaller and smaller in order toapproach the minimizers of G , which are in fact the desired solutions to our problem. ARTITIONS OF MINIMAL LENGTH ON MANIFOLDS 11
Since the parameters are of two types: triple points and edge points, we prefer to use an opti-mization algorithm which is not based on the gradient. The algorithm is described below. • For each point P consider a family of m tangential directions ( v i ) mi =1 chosen as follows:the first direction is chosen randomly and the rest are chosen so that the angles betweenconsecutive directions are 2 π/m . • Evaluate the cost function for the new partition obtained by perturbing the point P in eachof the directions v i according to a parameter ε . • Choose the direction which has the largest decrease and update the partition accordingly. • Do the same procedure for each edge point by performing the two possible orthogonalperturbations of the point with respect to the edge. • If there is no decrease for each of the points of the partition, then decrease ε .This algorithm converges in each of the test cases and the results are presented in Table 1. Inthe optimization procedure we start with ε = 1 and we reiterate the optimization decreasing ε by afactor of 10 at each step until we reach the desired precision on the area constraints. We are ableto recover the same results as Cox and Flikkema for n ∈ [4 , π/
3. In Figure 5 you can see the results for n = 9 and n = 20. The red arcs are geodesic connectingthe points and are drawn to visually see that not all the boundaries of the optimal structure aregeodesic arcs. our results Cox-Flikkema N non-geo. area tol. non-geo.4 11 . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − .
465 our results Cox-Flikkema N non-geo. area tol. non-geo.19 28 . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . . × − . Table 1.
Comparison between our results and the results of Cox and Flikkema inthe case of the sphere.Thus we can conclude that the relaxed formulation presented in the previous section is ableto match the best known configurations in the literature. Furthermore for n ∈ [5 , ∪ { } thealgorithm finds the good configuration without much effort, while for n ∈ [26 ,
31] multiple tries withdifferent initial conditions were needed in order to find the best configuration. The fact that thestructure of the partition is not fixed is a great advantage offered by our method.5.
Computing the optimal cost - general surfaces
The approach used in the previous section cannot be applied to other surfaces than the sphere.Indeed, the general expression of curves of constant curvature is not known explicitly for other types
Figure 5.
The difference between optimal configuration (black) and the geodesicsconnecting the points (red).of surfaces. One way to approximate the total perimeter of the partition would be to extract thecontours of the optimal densities and evaluate the length of each discrete contour. A natural way toextract a contour corresponding to a density function would be taking a level set, for example thelevel 0 .
5. It is possible to extract such level sets by looking at which triangles contain values whichare both above and below the level set. On each triangle which is cut by the contour we make alinear interpolation which determines a segment in the contour of the level set.Once we have an idea on how to extract the contours, the first question arises: how to make surethat the level sets extracted form a partition of S ? We denote by T a triangulation of S . If wethink of extracting the 0 . . T by φ i ( x ) = ( u i ( x ) ≥ max i = j u j ( x )0 otherwise , (5.1)where u i are the optimal densities obtained numerically. These contour levels of the functions φ i almost realize a partition of S with the following issues:(1) There is a small void space around each triple point, but this void is included in one of thetriangles of the mesh, and can be dealt with.(2) Since we extract the level sets of a function which is either 0 or 1 on the vertices of thetriangulation, the contour lines will pass through the middle of the edges of the trianglessituated at the border between two phases. This creates some contours which are quitezigzagged and whose length is significantly larger than the optimal total perimeter.We illustrate these two issues in Figure 6.Nevertheless, once we have extracted these contours it is possible to make a direct optimizationof the total length of the boundaries with the constraint of fixed area of the cells. This optimizationis made directly on the triangulated surface. We describe the optimization algorithm below. Variables and representation of the partitions.
We denote ( x i ) hi =1 a generic family ofvariable points situated each on an edge of the triangulation T such that each edge contains exactlyone variable point. To these points we associate a family of parameters ( λ i ) hi =1 which gives theposition of each point x i on the corresponding edges. We take this global parametric approach sinceeach of these points belongs to at least two cells and we’ll need to evaluate its contribution in thegradient of the area and the for all the cells that contain it. Having a global sets of points avoidshaving to match points between different contours. ARTITIONS OF MINIMAL LENGTH ON MANIFOLDS 13 v Figure 6.
A small space left around triple points (left) and the non-regular initialextracted contours (right).Each cell of the partitions is represented by a structure of pairs of edges of triangles of T whichdetermine, along with the parameters ( λ i ), the segments which form the discrete contour of the cell.The pairs of edges is ordered so that the contour is continuous. Contours may have one or moreconnected components. Computation of the perimeters of the cells.
The perimeter of a cell is computed by followingthe segments forming the contour and incrementally adding their lengths to the total length. If thevertices of the segment are given by x i = λ i v + (1 − λ i ) v and x j = λ j v + (1 − λ j ) v then thelength of the segment [ x i , x j ] is ℓ ([ x i , x j ]) = k λ i v + (1 − λ i ) v − λ j v − (1 − λ j ) v k , expression which is differentiable if the length is not zero. The derivatives with respect to λ i and λ j are then added to the gradient vector. Note that for the points which are not vertices of somecontour the gradient is zero. Computation of the areas of the cells.
In order to compute the area of a cell we use theinformation given by the functions φ i defined in (5.1). The function φ i shows, among other things,what is the position of each triangle in T with respect to the cell i . Indeed, denoting by T a trianglein T , we have the following cases:(1) All the vertices v of the triangle T satisfy φ i ( v ) = 1. Then T is completely inside the cell i and we add its area to the total area of the cell.(2) Two vertices v , v of T satisfy φ i ( v , ) = 1 and the third satisfies φ i ( v ) = 0. Thus we onlyadd a portion of the area of T to the total area of cell i . Note that this value of the areadepends linearly of one parameter λ k and of another parameter λ l . The derivatives of thesecontributions are added to the vectors containing the gradient of the area of the cell i .(3) Two vertices v , v of T satisfy φ i ( v , ) = 0 and the third satisfies φ i ( v ) = 1. Again, weonly add a portion of the area of T to the total area of cell i which again depends linearly ofone parameter λ k and of another parameter λ l . The derivatives of these contributions areadded to the vectors containing the gradient of the area of the cell i .(4) If all the vertices of T satisfy φ i ( v ) = 0 then the triangle is outside the cell and we move on. The empty spaces around triple points.
As we have noted above and seen in Figure 6,around triple points we have some empty spaces determined by three points which belong to thethree sides of some of the triangles in T . In each configuration of this type we add a Steiner treecorresponding to the three variable points. Each of the three area regions which are formed are addedto the corresponding cell while the perimeter is modified with the length of two adjacent segmentsin the Steiner tree. See Figure 7 for further details. In order to find the gradient corresponding tothe lengths and area changes due to the addition of these Steiner points we use a finite differencesapproximation. Figure 7.
Treatment of empty space around triple points. We consider theFermat point X of the empty triangle ABC and we add corresponding area andperimeters to the corresponding cells. For example the area of
ABX is added toCell 3 and the quantity AX + BX − AB is added to the perimeter of Cell 3. Figure 8.
Contours after the constrained optimization algorithm. You can alsosee a zoom around the triple points: the segments which join the Fermat pointsalign themselves with the rest of the contour.
Constrained optimization algorithm.
We have the expressions and the gradients of theperimeters and areas of the cells as functions of the parameters ( λ i ) hi =1 . This allows us to use thealgorithm fmincon from the Matlab Optimization Toolbox in order to implement the constrainedoptimization algorithm. We use the interior-point algorithm with a low-memory hessian approxi-mation given by an LBFGS algorithm. The initial values of the parameters ( λ i ) hi =1 are all set to0 .
5. The algorithm manages to satisfy the constraints at machine precision while minimizing theperimeter and thus smoothing the zigzagged initial contours (like the ones in Figure 6). An exampleof result may be seen in Figure 8.It may be the case that some vertices of the contour would ”like” to switch to another side. Thiscan be the case if at the end of the optimization one of the parameters λ i is close to 0 or 1 or atriple point in one of the constructed Steiner trees is on the boundary of the corresponding meshtriangle. In this cases we modify the initial contours taking into the account these results and werestart the optimization procedure. The modification is done in the following way. ARTITIONS OF MINIMAL LENGTH ON MANIFOLDS 15 n Minimal length2 15 .
073 22 .
614 30 .
155 37 .
256 41 . n Minimal length7 47 .
128 50 .
779 53 . . Table 2.
Approximation of the optimal costs for minimal partitions of a torusinto equal area cells. These partitions are represented in Figure 2(1) If one of the λ i is equal to 0 or 1 then we add the corresponding point to the adjacent celland restart the algorithm.(2) If one of the triple points arrives on the edge of its corresponding mesh triangle then weallow it to move to the adjacent triangle.After a finite number of switches the configuration stabilizes and a local minimum is found.We test the presented algorithm on the results obtained in previous sections. In the case of thesphere we obtain the same values found in Table 1. The approximations of the optimal costs forpartitions presented in Figure 2 for a torus of radii R = 1 , r = 0 . Conclusions
We propose an algorithm for finding numerically the partitions which divide a surface into cells ofprescribed areas and minimize the sum of the corresponding perimeters. This algorithm is rigorouslyjustified by a Γ-convergence result which is a generalization of the Modica-Mortola theorem in thecase of smooth ( d − a priori . The cells emerge from random density configurations and place themselves inthe best positions.The Γ-convergence method is not limited to the case of the sphere. Once we have triangulated asurface the same algorithm applies. We present a few test cases of more complex surfaces. While therelaxed optimal partitions can easily be obtained, computing the optimal costs is not straightforwardsince the relaxed costs are not precise enough. In order to be able to compute an approximation ofthese optimal costs we extract the contours of the optimal densities and we perform a constrainedoptimization on the triangulated surface. References [1] Giovanni Alberti. Variational models for phase transitions, an approach via gamma-convergence. 1998.[2] Luigi Ambrosio and Andrea Braides. Functionals defined on partitions in sets of finite perimeter. II. Semiconti-nuity, relaxation and homogenization.
J. Math. Pures Appl. (9) , 69(3):307–333, 1990.[3] Luigi Ambrosio, Nicola Fusco, and Diego Pallara.
Functions of bounded variation and free discontinuity problems .Oxford Mathematical Monographs. The Clarendon Press, Oxford University Press, New York, 2000.[4] Sisto Baldo and Giandomenico Orlandi. Cycles of least mass in a Riemannian manifold, described through the“phase transition” energy of the sections of a line bundle.
Math. Z. , 225(4):639–655, 1997.[5] Felix Bernstein. ¨Uber die isoperimetrische Eigenschaft des Kreises auf der Kugeloberfl¨ache und in der Ebene.
Math. Ann. , 60(1):117–136, 1905.[6] Andrea Braides.
Approximation of Free-Discontinuity Problems . Springer, 1998. [7] Kenneth A. Brakke. The surface evolver.
Experiment. Math. , 1(2):141–165, 1992.[8] Giuseppe Buttazzo. Gamma-convergence and its Applications to Some Problems in the Calculus of Variations.
School on Homogenization ICTP, Trieste, September 6-17 , 1993.[9] 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.[10] Manfredo P. do Carmo.
Differential geometry of curves and surfaces . Prentice-Hall, Inc., Englewood Cliffs, N.J.,1976. Translated from the Portuguese.[11] Max Engelstein. The least-perimeter partition of a sphere into four equal areas.
Discrete Comput. Geom. ,44(3):645–653, 2010.[12] Thomas C. Hales. The honeycomb conjecture.
Discrete & Computational Geometry , 25(1):1–22, 2001.[13] Thomas C. Hales. The honeycomb problem on the sphere, 2002.[14] Antoine Henrot and Michel Pierre.
Variation et optimisation de formes , volume 48 of
Math´ematiques & Appli-cations (Berlin) [Mathematics & Applications] . Springer, Berlin, 2005. Une analyse g´eom´etrique. [A geometricanalysis].[15] Joseph D. Masters. The perimeter-minimizing enclosure of two areas in S . Real Anal. Exchange , 22(2):645–654,1996/97.[16] Frank Morgan. Soap bubbles in R and in surfaces. Pacific J. Math. , 165(2):347–361, 1994.[17] ´Edouard Oudet. Approximation of partitions of least perimeter by Γ-convergence: around Kelvin’s conjecture.
Exp. Math.
Laboratoire Jean Kuntzmann, Universit´e Grenoble Alpes, BˆatimentIMAG, 700 avenue centrale, 38400 Saint Martin d’H`eres France
E-mail address , Beniamin Bogosel: [email protected]
E-mail address , ´Edouard Oudet:, ´Edouard Oudet: