A consistent approximation of the total perimeter functional for topology optimization algorithms

This article revolves around the total perimeter functional, one particular version of the perimeter of a shape [[EQUATION]] contained in a fixed computational domain D
measuring the total area of its boundary [[EQUATION]], as opposed to its relative perimeter, which only takes into account the regions of [[EQUATION]] strictly inside [[EQUATION]].
We analyze approximate versions of the total perimeter which make sense for general ``density functions'' [[EQUATION]]. Their use in the context of density-based topology optimization is particularly convenient as they do not involve the gradient of the optimized function [[EQUATION]]. Two different constructions are proposed: while the first one involves the convolution of the function [[EQUATION]] with a smooth mollifier, the second one relies on the resolution of an elliptic boundary-value problem featuring Robin boundary conditions. The ``consistency'' of these approximations is appraised from various points of view. At first, we prove the pointwise convergence of our approximate functionals, then the convergence of their derivatives, as the level of smoothing tends to [[EQUATION]], when the considered density function [[EQUATION]] is the characteristic function of a shape [[EQUATION]]. Then, we focus on the [[EQUATION]]-convergence of the second type of approximate total perimeter functional. Several numerical examples are eventually presented.

. The relative perimeter Per R (Ω) of the set Ω ⊂ D is the length of the yellow curves, while its total perimeter Per T (Ω) is the sum of the lengths of the yellow and red curves.
Here, the optimized design Ω is a domain in R d (d = 2, 3); as often, it is contained in a fixed "hold-all" domain D ⊂ R d , which serves as computational domain in numerical applications. The objective criterion J(Ω) measures the physical performance of the shape Ω, and P (Ω) is a constraint -enforced by a fixed penalization of J(Ω) to set ideas -which is for instance related to its geometric features.
Among such geometric constraint functionals, the perimeter Per(Ω), measuring the length of the contour ∂Ω of Ω in 2d or the area of the surface ∂Ω in 3d, has received a particular attention. From the theoretical point of view, including a perimeter constraint into a shape optimization problem is a well-known means to overcome the homogenization phenomenon, and thereby to ensure the well-posedness of the problem [10,46]. In numerical practice, perimeter constraints tend to prevent the emergence of sharp features, thus promoting designs with "smooth" boundaries, see e.g. [64] about this point and for a more general discussion about common sources of numerical instabilities in shape and topology optimization and possible remedies.
Before proceeding, let us stress that when the considered shapes Ω are subsets of a fixed domain D, the terminology "perimeter" actually encompasses two subtly different notions, which are illustrated in Figure 1: -The relative perimeter Per R (Ω) of a shape Ω ⊂ D is defined by: (1.1) the term "relative" reflects the fact that Per R (Ω) only takes into account the part of the boundary ∂Ω lying strictly inside the "hold-all" domain D. -The total perimeter Per T (Ω) of Ω ⊂ D reads: Per T (Ω) = ∂Ω ds, (1.2) so that Per T (Ω) takes into account the way ∂Ω intersects ∂D.
Obviously, when D stands for the whole ambient space R d , the above two notions coincide with the more familiar perimeter of a subset Ω ⊂ R d : Despite their intuitive definition and their "simple" mathematical formulations, perimeter functionals raise interesting challenges depending on the elected numerical framework for optimal design. Let us briefly describe a few of them, which are frequently encountered in the literature.
-In shape optimization problems where, as in the above model example, the optimized design is a domain Ω ⊂ D, the numerical treatment of perimeter constraints naturally brings into play the shape derivatives of the functionals Ω → Per R (Ω) or Ω → Per T (Ω). These are expressed in terms of the mean curvature κ of ∂Ω, which is difficult to evaluate precisely, see e.g. Section 3.3.2 in [42] for a discussion about this issue. -In the same framework, the topological derivative of Per R (Ω) or Per T (Ω), measuring the effect of the nucleation of a tiny hole inside Ω, is also uneasy to handle. Indeed, the sensitivity of the perimeter functions is of lower order in terms of the size of the vanishing hole, when compared with that of most "usual" functionals of the domain; see for instance [44,58]. -In the relaxed setting of topology optimization, where the design is accounted for by a density function u : D → [0, 1] defined on the fixed domain D, a consistent counterpart to the perimeter functionals Per R (Ω) or Per T (Ω) has to be found, which makes sense for such density functions u (and not only for characteristic functions χ Ω of "true" shapes Ω ⊂ D), and enjoys "nice" algorithmic properties.
The latter task of devising suitable approximate perimeter functionals, defined for general density functions u : D → [0, 1] has been the focus of many investigations in the literature, in the particular instances where the bounding box D is the whole ambient space R d , or where D is a strict subset of R d and the version of the perimeter at hand is Per R (Ω); we sketch an overview of the existing literature about these subjects in Section 2.3 below.
On the contrary, very little attention has been paid to the total perimeter Per T (Ω), to the best of our knowledge. It turns out that this perimeter functional is rarely used in numerical practice; often, its relative counterpart is tacitly retained, as in the practice of the level set method for shape optimization [7]. This is actually quite surprising, since in a fair number of applications, the total perimeter appears more natural than its relative counterpart: the boundary of the "hold-all" domain D has to be taken into account in the evaluation of the perimeter of the shape Ω when it has a "physical meaning". This is the case, for instance, in situations involving contact mechanics, where at least one part of ∂D is in contact with an external device; see [50] about this physical setting and [51] for examples of related shape optimization problems.
The present article is a contribution to the mathematical and numerical approximation of the total perimeter functional Per T (Ω) in density-based topology optimization frameworks. More precisely, we aim to construct approximate total perimeter functionals P ε (u) with "nice" algorithmic properties, and which are "consistent" with Per T (Ω) as the small parameter ε vanishes. These are defined for a general density function u : D → [0, 1] and they involve a smoothing of the latter, whose magnitude is controlled by ε. Two such constructions are proposed: on the one hand, the smoothing of u is realized by convolution with a smooth mollifier, and on the other hand, it is based on the regularizing effect of an elliptic Partial Differential Equation (PDE) supplemented with Robin boundary conditions. The adequate and consistent characters of our approximate perimeter functionals will be given precise meanings in Section 2.3 below, where our main results are described in more technical details.
This article is organized as follows. In Section 2, we present the shape optimization problems of interest in this article, in an informal way. We then introduce density-based reformulations of such problems, and we overview the contributions from the literature devoted to the perimeter functional in this context, before describing ours. In Section 3, we recall some necessary background material about functions with bounded variations and sets with finite perimeter, and then introduce our approximate total perimeter functionals. The next Section 4 discusses the quality of these approximations in terms of their pointwise convergence and of the convergence of their derivatives in the particular case where the considered density u is the characteristic function χ Ω of a smooth shape Ω. Section 5 then deals with the Γ-convergence of the PDE-based approximate total perimeter functional. In Section 6, we slip into the context of numerical shape and topology optimization: several numerical examples are presented to appraise the main points of our study. In Section 7, some conclusions and perspectives arising from our study are outlined. The article ends with a technical appendix.

Shape and topology optimization under perimeter constraints
This section introduces the shape and topology optimization setting which serves as a motivation and guide throughout the paper. After introducing informally shape optimization problems and recalling basic related facts in Section 2.1, density-based reformulations of such problems are broached in Section 2.2. The available methods from the literature for dealing with perimeter functionals in the context of density-based topology optimization are then summarized in Section 2.3 and the main contributions of the present article are described.

Classical stakes about shape optimization
Let us consider a model shape optimization problem of the form: min Ω⊂D J(Ω) + P (Ω). (2.1) In the above formulation, -The shape Ω is sought as a subset of a fixed computational domain D in R d (d = 2 or 3).
-The objective function J(Ω) brings into play the physical or the mechanical behavior of the optimized shape Ω (as described by e.g. the heat equation, the linearized elasticity system, etc.). We are not too specific for the moment about the precise definition of J(Ω), and we shall present several instances of such in Section 6.4.1. -The functional P (Ω) stands for either of the two perimeter functionals Per R (Ω), Per T (Ω) introduced in the previous section.
Admittedly, there are multiple variants of (2.1), featuring for instance P (Ω) as a constraint (instead of a mere penalization of the objective function J(Ω)), or bringing into play additional constraints on Ω (e.g. on its volume). Shape optimization problems of this kind will be studied numerically in Section 6, but for the simplicity of the exposition, we stick to the simple formulation (2.1) for the moment.
Remark 2.1. The penalization of J(Ω) by the perimeter functional in (2.1) is a well-known stratagem to enforce the existence of an optimal shape, at least when particular instances of the objective function J(Ω) are considered, see e.g. [10] when J(Ω) is the compliance functional, and the work [37] dedicated to spectral functionals.
Most numerical methods for problems of the form (2.1) rely on the derivatives of the objective and constraint functions with respect to the domain. This notion may be defined in a variety of frameworks, and we retain Hadamard's boundary variation method, giving rise to the concept of shape derivative; see for instance [9,46,57,58] about the alternative notion of topological derivative. Briefly, and omitting for simplicity the requirement that the considered shapes Ω should be contained inside D, variations of a bounded, Lipschitz domain Ω ⊂ R d are considered under the form: Ω θ := (Id + θ)(Ω), θ ∈ W 1,∞ (R d , R d ), ||θ|| W 1,∞ (R d ,R d ) < 1. (2.2) One function of the domain F (Ω) is then said to be shape differentiable at Ω if the underlying mapping θ → F (Ω θ ), from W 1,∞ (R d , R d ) into R, is Fréchet differentiable at θ = 0. The corresponding derivative F (Ω)(θ) is the shape derivative of F (Ω) at Ω, and the following expansion holds: From the numerical vantage, shape and topology optimization algorithms can be based on shape derivatives, using for instance mesh deformation algorithms as in [8,59], or the level set method as in [7,61,66] for the numerical representation of the shape and its evolution; see [6] for a recent overview.
Let us conclude this section by recalling the following classical result about the shape differentiability of perimeter functionals (see e.g. [46], Sect. 5.4.2): Let Ω be a bounded domain of class C 3 . The mapping θ → Per(Ω θ ) from C 2 (R d , R d ) ∩ W 1,∞ (R d , R d ) into R is Fréchet differentiable at θ = 0, with shape derivative: where κ : ∂Ω → R is the mean curvature of ∂Ω (see Appendix A.1). Remark 2.3. Similar results hold for the relative and total perimeter functionals, except that the considered deformations θ in the practice of Hadamard's method (see (2.2)) should satisfy θ · n = 0 on the boundary ∂D of the computational domain.

Density-based reformulation of shape optimization problems
Optimal design problems of the form (2.1) pose difficulties of various natures. From the mathematical viewpoint, the existence of an optimal shape Ω is often difficult to guarantee; in addition, the expression of the optimality conditions for (2.1) requires the technical calculation of shape or topological derivatives, as the building blocks of optimality conditions, see e.g. [9,46] about these points. From the numerical viewpoint, the implementation of efficient shape optimization algorithms is far from straightforward [6,63].
These challenges legitimate the popularity of density-based reformulations of (2.1). In a nutshell, the considered set of designs is extended from true "black-and-white" shapes Ω ⊂ D (or equivalently their characteristic functions χ Ω , taking values 1 in Ω and 0 in D \ Ω) to more general density functions u : D → [0, 1], which are allowed to take intermediate "grayscale" values between 0 and 1.
The problem (2.1) is recast as: Here, J(u) and P (u) stand for "suitable" extensions of the objective and perimeter functions J(Ω) and P (Ω) in (2.1) to density functions u, that is, J(u) = J(Ω) and P (u) = P (Ω) when u = χ Ω is the characteristic function of a shape Ω ⊂ D. The formulation (2.3) certainly presents a number of advantages over (2.1); in particular, it is much simpler to handle numerically: for instance, algorithmic strategies based on the derivatives of the mappings u → J(u) and u → P (u) are easy to carry out. This question of how to conveniently extend an objective function J(Ω) accounting for the "physical performance" of Ω (that is, involving the solution to a partial differential equation posed on Ω such as the heat equation or the linear elasticity system) to accommodate density functions has been extensively studied in the literature. Without entering into details, let us solely mention that the homogenization theory is often invoked [4] to this end, or simplified, heuristic versions are used, such as the famous Solid Isotropic Material with Penalization (SIMP) method in structural mechanics [25], the so-called "porosity method" in fluid mechanics [26], etc. see also our previous works [5,16,19] in this direction, where approximate versions J(u) of typical functionals J(Ω) of the domain are constructed which are "consistent" with the shape and topological derivatives of the latter.
2.3. Perimeter functionals in density-based topology optimization: existing works and contributions of the present article A great deal of work has been devoted to the relative perimeter functional Per R (Ω) and its extensions to the context of density functions. One natural mathematical setting features functions u ∈ BV (D) with bounded variations inside the computational domain D, see Section 3.1 below for a brief reminder about this notion. The natural density-based counterpart of Per R (Ω) then reads: where |Du|(D) is the total variation of u ∈ BV (D). Unfortunately, this neat mathematical framework is difficult to exploit from the numerical vantage, since functions u ∈ BV (D) generally present jumps, which are difficult to capture. This drawback has aroused interest for appropriate, smoothed versions F ε (u) of the (extended) relative perimeter P R (u), the perhaps most famous of which being the Cahn-Hilliard functional, defined for u ∈ L 1 (D) by: Here, the function R s → W (s) ∈ R is a potential with two wells in s = 0 and s = 1, satisfying suitable growth conditions. The functional F ε (u) essentially assumes H 1 (D) functions and ε > 0 is a "small" parameter tuning the level of smoothing. Intuitively, a minimizer u of F ε (u) is expected to strike a good compromise between showing "not too steep" variations (which are penalized by the gradient term) while avoiding to take values "too different" from 0 and 1 (those being penalized by the potential W (u)). It was proved by Modica and Mortola [54,55] that F ε Γ-converges to the functional F (u) defined by where c W is a real-valued constant which explicitly depends on the potential W . We refer to Section 5.1 below for a short reminder of the notion of Γ-convergence; see also [2] for an enlightning discussion, and [30] for different instances of the Modica-Mortola approximation. The approximate relative perimeter functional (2.5) (or variants of it) has been a popular means to formulate constraints on the length of codimension 1 subsets in density-based optimal design problems. Let us mention a few such applications: -In image segmentation, the objects composing an image, supplied as a grayscale intensity function I : D → [0, 1], can be identified by minimizing the Mumford-Shah functional, defined over functions u ∈ BV (D, [0, 1]). The latter is a weighted sum of a proximity term between u and I and the measure of the jump set S u of the optimized function; see [24] for an overview. The so-called Ambrosio-Tortorelli functional approximates this jump part by a modified version of (2.5) [14,15]. -A quite similar approach is used in [28] in the context of an optimal repartition problem of a solid phase, made of an elastic material, a liquid phase and void. These three phases are described by a single "phasefield" function u ∈ BV (D), whose values are imposed to be close to −1, 0, or 1 by the presence of a three-well potential in the minimized energy. Furthermore, the length of the interface between each phase is penalized in a manner reminiscent to the aforementioned Modica-Mortola procedure.
-The functional F ε (u) in (2.5) also makes it possible to approximate the gradient flow of the perimeter functional -the mean curvature flow -by the gradient flow of the minimization of u → F ε (u) -the Allen-Cahn equation. A degenerate case of this last equation leads to the Bence-Merriman-Osher algorithm, also referred to as "threshold dynamics". In a nutshell, a shape Ω is updated by alternating convolutions of the characteristic function of Ω with a kernel (for instance, the fundamental solution to the heat equation) with thresholdings of the resulting function. We refer to [52] for the seminal article about threshold dynamics, to [39] for an extension to the multi-phase context and to e.g. [32] and references therein for more details on these approaches.
The use of the approximate relative perimeter functional (2.5) is certainly more convenient than that of the exact formula (2.4) from a practical viewpoint. Yet, a number of numerical difficulties persist, mainly because the defining formula (2.5) involves the gradient of the optimized function u. As a consequence, (2.5) only makes sense when u belongs to H 1 (D), and so its use in optimal design prevents u from being a characteristic function; instead, u will inevitably present intermediate "grayscale" values between 0 and 1. Moreover, the derivative of u → F ε (u), which is the key ingredient of most first-order optimization algorithm for (2.3) (where P (u) is replaced by F ε (u)), involves second-order derivatives of u. Hence, a high-order discretization has to be used for u to ensure the accuracy of the computation of this derivative, or a sophisticated semi-implicit scheme has to be used to update u, as in [28].
In order to circumvent these concerns, elaborating on ideas present in [3,53] (see also [39] in the multi-phase context), an alternative approximation of the relative perimeter P R (u) was proposed in [21]: Using the jargon of topology optimization, the optimized design u is filtered thanks to the regularizing effect of the elliptic operator (2.7); see [27] for a theoretical justification and [34] for a related numerical investigation. The analysis in [21] reveals that F ε Γ-converges in L ∞ (D, [0, 1]) (when the latter is equipped with the L 1 (D) distance) to the following multiple of the extended relative perimeter: The above contributions are dedicated to the relative perimeter Per R (Ω) of a domain Ω, or its extension (2.4) to density functions u. As we have mentioned, in sharp contrast with this abundant literature, very few contributions have been dealing with the total perimeter functional Per T (Ω), to the best of our knowledge.
The latter is readily extended to the density-based topology optimization setting via the following formula: where u stands for the extension of u by 0 outside D. Understandably enough, this extended function suffers from the same drawbacks as its relative counterpart in (2.4).
The main goal of the present work is precisely to construct, analyze and put into practical use suitable approximate versions of (2.8). More precisely, we introduce approximate versions of the total perimeter P T (u) whose common structure resembles very much that of (2.6): In the above definition, L ε is a regularizing operator, about which two versions are proposed: -One of them, L conv ε u = η ε * u is based on the convolution of u with a smooth mollifier η ε ; -The second one, L Rob ε , resembles much the operator in (2.6), except that the homogeneous Neumann boundary condition in (2.7) is replaced by a Robin boundary condition: L Rob ε u is now the unique solution v ε ∈ H 1 (D) to the boundary-value problem: We denote by P conv ε (u) (resp. P Rob ε (u)) the instance of (2.9) corresponding to the choice of the operator L conv ε (resp. L Rob ε ). These constructions echo to our previous contributions [5,16,19], where density-based approximate topology optimization problems were introduced, showing "consistent" optimality conditions with those of the original shape optimization problem.
We then analyze thoroughly the mathematical properties of the functions P ε (u) defined by (2.9), as well as the quality of the induced approximation of P T (u). Our main contributions can be summarized as follows: (i) We first address the pointwise convergence of the functions P ε (u): when u = χ Ω is the characteristic function of a "smooth enough" shape Ω ⊂ D, we prove that both quantities P conv ε (u) and P Rob ε (u) converge (up to an explicit multiplicative constant) to the total perimeter Per T (Ω) as ε → 0, see Sections 4.1 and 4.2. This result is generalized later in Section 5, where we show that, in the case of the functional P Rob ε (u), the pointwise convergence of P Rob ε (u) to (a multiple of) P T (u) actually holds when u is an arbitrary function in BV (D, {0, 1}). (ii) We then consider in Section 4.3 the convergence of the "derivative" g u,ε of the mapping u → P ε (u) to the shape derivative of Ω → Per(Ω) as ε → 0 in both cases where L ε = L conv ε and L ε = L Rob ε . When u = χ Ω is the characteristic function of a "smooth enough" shape Ω ⊂ D, we prove that g u,ε (x) converges (up to a multiplicative constant) to the mean curvature κ(x) of Ω at all points x ∈ ∂Ω ∩ D. This convergence result is important in practice since it accounts for the convergence of the optimality conditions attached to P conv ε (u) and P Rob ε (u) to those corresponding to Per T (Ω), and also for the consistency of the associated gradient flows. (iii) Finally, we show in Section 5 that the PDE-based approximate total perimeter P Rob ε (u) Γ-converges to (a multiple of) P T (u); this feature accounts for the behavior of the minimizers of P Rob ε (u) as ε → 0, or more generally for that of solutions to optimization problems involving P Rob ε (u) such as (2.3), see Remark 5.6 for further comments about this point.
Last but not least, these results are illustrated in Section 6 by several numerical examples, ranging from mere validation experiments to more concrete resolutions of shape and topology optimization problems.

The total and relative perimeter functionals and their approximations
In this section, we properly introduce our approximate total perimeter functionals. Since they are built in the density-based setting introduced in Section 2.2, we first recall in Section 3.1 some basic material about the chosen mathematical framework, namely that of functions with bounded variations and sets with finite perimeter. Meanwhile, we set some notations used throughout the article. The approximate total perimeter functionals are then formally presented in Section 3.2.
Thenceforth, D stands for a bounded "hold-all" domain of R d (where d = 2, 3 in applications), with at least Lipschitz regularity.

Preliminaries about geometric measure theory
Our main reference about the topics discussed in this section is [22]; see also [12].

Functions with bounded variations
Let us start with a definition. Definition 3.1. A function u : D → R with bounded variations on D is a function u ∈ L 1 (D) whose distributional gradient Du, defined by: We denote by BV (D) the space of functions with bounded variations in D.
Note that, as in Chapter 1 of [12], this definition implicitly assumes that the total variation of Du is bounded when u ∈ BV (D), that is: As a result of the last fact, the integral D f · Du can be defined for all continuous and bounded functions f ∈ C(D, R d ).
The quantity |Du|(D) (be it finite or infinite) is called the total variation of a function u ∈ L 1 (D); since it reads as the supremum of a family of linear functionals, the mapping L 1 (D) u → |Du|(D) ∈ R ∪ {∞} is lower semi-continuous.
The following proposition asserts that smooth functions are dense in the space BV (D), when an adapted notion of convergence is considered on the latter space, see Theorem 10.1.2 in [22]. In this context, an adapted version of the Green's formula holds for functions u ∈ BV (D), (see [22], Thm. 10.2.1): Proposition 3.4. Let u be a function in BV (D); then for any function ϕ ∈ C 1 (D, R d ), (3.1) Here and thenceforth, we denote by H t the t-dimensional Hausdorff measure on R d for any non negative real number t ≥ 0. We recall that H d−1 is simply the surface measure ds when it is restricted to the boundary ∂D, while H d coincides with the Lebesgue measures L d of Borel subsets of R d . The vector field n D in (3.1) (or simply n when the context is clear) is the unit normal vector to ∂D, pointing outward D.
Throughout the paper, we denote by v the extension by 0 to R d of a function v ∈ L 1 (D); the previous identity then rewrites equivalently: and so we have the following identity between R d -valued measures on R d :

Approximate limits
Let E ⊂ R d be a Borel subset, and let u : R d → R be a measurable function. if and only if for every ε > 0, the following relation holds: Here and throughout the text, {|f − | < ε} is a shortcut for the set of points x ∈ R d such that |f (x) − | < ε. Also, B ρ (x 0 ) is the open ball with center x 0 and radius ρ.
For x 0 ∈ R d and a given unit vector a in the unit sphere S d−1 ⊂ R d , we introduce the half-space oriented by a: Definition 3.6. One point x 0 ∈ R d is said to be a regular point of the measurable function u if there exists a direction a ∈ S d−1 such that the following two approximate limits exist:  The following proposition sheds some light about the nature of regular points; see ( [22], Thm. 10.3.1).
Definition-Proposition 3.7. Let x 0 be a regular point of the measurable function u and let the direction a ∈ S d−1 be as in (3.4); two situations may occur: -If u a (x 0 ) = u −a (x 0 ), then for any direction b ∈ S d−1 , the approximate limit u b (x 0 ) exists and equals u a (x 0 ); x 0 is said to be a point of approximate continuity of u; -If u a (x 0 ) = u −a (x 0 ), then a is the unique vector (up to a sign) in S d−1 such that u ±a (x 0 ) exist; x 0 is said to be a jump point of u.
We denote by S u the jump set of u, i.e. the set of its jump points.

Sets with finite perimeter
In the following, we shall be interested in those functions in BV (D) which only take the values 0 and 1, that is, the characteristic functions χ E of Borel subsets E ⊂ D. Those sets E ⊂ D such that χ E belongs to BV (D, {0, 1}) are called subsets of D with finite (relative) perimeter.
The topological boundary ∂E of a set E with finite perimeter in D may be a quite singular object; see for instance ([45], Sect. 1.10). Two other notions of boundary are actually better suited for the study of such sets: -The first one of them is the measure theoretical (or essential) boundary ∂ M E, defined for a general Borel set E ⊂ R d , not necessarily contained in D, by: -The second notion of importance is the jump set S χ E ∩ D of the characteristic function of E. In this case, it coincides with the set of points x 0 ∈ D where there exists a direction a ∈ S d−1 such that: Finally, we list below a few useful properties of sets E with finite perimeter in D, see ( [22], Thm. 10.3.2).
Proposition 3.8. Let E ⊂ D be a set with finite perimeter; then -There exists a vector field ν(x) with values in S d−1 , the inner measure theoretic normal vector to E, defined for H d−1 almost every point x ∈ ∂ M E, such that: In particular, It follows from (3.2) that ν(x) coincides with −n(x) at points x ∈ ∂ M E ∩ ∂D, where n is the unit normal vector to ∂D pointing outward D.

Perimeter functionals in the setting of functions with bounded variations
Let us now consider the functional defined for a general function with bounded variations in R d by: (3.7) Note that we have used (3.2) in passing from the first to the second line of the above equality. These functions are tailored in such a way that, if Ω is a bounded Lipschitz domain of R d , and if Ω is a bounded Lipschitz subset of D, Per R (Ω) = P R (χ Ω ), and Per T (Ω) = P T (χ Ω ).
Hence, with a small abuse of language, these functions will still be referred to as the perimeter (resp. the relative and total perimeter) of u.

Definition of the approximate perimeter functionals
We are now in position to introduce approximate versions P ε (u) of the total perimeter P T (u). Two different constructions are presented, with competing assets and drawbacks. They have the common structure and when possible, we shall deal with them in a unified fashion. In the above, L ε : is a regularizing operator of the integral form The kernel N ε (x, y) has the decomposition: where -the function µ ∈ L 1 (R d ) is radial, smooth on R d \ {0}, and it has unit integral: Besides, µ has (at most) a weak singularity near 0: and it has fast decay at infinity: -For any compact set K D, the remainder R ε (x, y) satisfies for some positive constants C and α depending on K.
We now discuss two possible definitions L conv ε and L Rob ε of the regularizing operator L ε , which induce via (3.8) two different approximate total perimeter functions P conv ε and P Rob ε .

Approximate perimeter functionals based on a mollifying kernel
Let η : R d → R be a smooth mollifier, that is: -η is of class C ∞ ; -η has compact support inside the unit ball B 1 (0); -η has unit integral R d η dx = 1; -η has radial symmetry: there exists a smooth function η : (3.14) For ε > 0, we set η ε (x) = 1 ε d η( x ε ) and we introduce the convolution operator Inspired from [17,21], we now define the approximation P conv ε (u) of the total perimeter P T (u) by: Remark 3.9.
-A similar approximation could be worked out for the absolute perimeter functional P (u), but not for the relative perimeter P R (u). Indeed, doing so would rely on a procedure to extend u from D to R d without creating "artificial" jumps across ∂D. -We shall see that, despite its intuitive character, the convolution-based functional P conv ε (u) is not so easy to handle as the PDE-based counterpart presented in the next section, from both theoretical and numerical viewpoints.

Approximate functionals using the regularizing effect of an elliptic partial differential equation
Our PDE-based approximation of the total perimeter P T (u) is based on the following operator: for a given function u ∈ L ∞ (D, [0, 1]), let where v ε ∈ H 1 (D) is the unique solution to the boundary value problem: The approximate total perimeter P Rob At the moment, it is not obvious that the functional P Rob ε has the desired structure (3.8). To see this more clearly, we rely on the integral representation of v ε : In the above formula, N ε (x, y) is the fundamental solution of the equation (3.17): This function Φ ε has the analytic expression (iv) (Estimates of the derivative near 0) For m > 0, there exists a constant C > 0 depending on m such that (v) (Estimates at ∞) For M > 0 and a multi-index α ∈ N d , there exist constants a > 0, C > 0 depending on M and p such that The function N ε (x, y) can now be decomposed as: where for given x ∈ D, the remainder y → R ε (x, y) satisfies: It then follows from the structure (3.19) of the fundamental solution Φ ε , the decay of modified Bessel functions at infinity (see (v) in Lem. 3.10) and classical elliptic estimates that for a given compact subset K D, for some positive constants a, C depending on K.
Remark 3.11. A similar approximate counterpart was proposed for the relative perimeter (3.7) in [21]: (3.20) in which w ε is the H 1 (D) solution to the problem Let us emphasize the flexibility of the use of PDE-based regularizing operators L ε in the device of our gradientfree perimeter functionals (3.8): an approximate version of the total perimeter is obtained from its relative counterpart by simply changing boundary conditions in the defining boundary-value problem.

Pointwise convergence properties of the approximate total perimeter functionals under a regularity assumption
We first discuss the convergence properties of the approximate functionals P conv ε (u) and P Rob ε (u) towards the total perimeter functional P T (u) in the particular case where u = χ Ω is the characteristic function of a domain  Roughly speaking, part (iii) of the above assumption means that the boundaries ∂O and ∂D are not exactly parallel at the points where they intersect. In particular, Ω is Lipschitz regular, see Figure 2 for an illustration. We start in Section 4.1 by proving the pointwise convergence P conv ε (u) → P T (u) as ε → 0 when u = χ Ω is the characteristic function of a domain Ω ⊂ D. Then, in Section 4.2, we deal with the pointwise convergence of P Rob ε (u). Note that the results of this section will be retrieved in a more general context in Section 5.5; yet, we believe that the calculations conducted in here are interesting per se, since they allow for a clearer intuition of the convergence process. We finally examine in Section 4.3 the behavior of the derivatives of the mappings u → P conv ε (u) and u → P Rob ε (u) as ε → 0.

4.1.
Pointwise convergence of the convolution-based functional P conv ε (u) towards P T (u) Let us start with a technical lemma: there exists a function µ : R → R such that η(x) = µ(|x|). The following identity holds: where y := (y 1 , . . . , y d−1 ) ∈ R d−1 is the collection of the (d − 1) first components of a point y = (y 1 , . . . , y d ) ∈ R d , and the function η : R → R is defined by: In the following, we denote by C(η) the common value of both sides of (4.2).
Proof. From the definition (4.3) of η, we obtain: where we have operated a switch to polar coordinates in R d−1 to pass from the second to the third line, denoting by ω d−1 the surface of the unit sphere S d−2 ⊂ R d−1 . Using the change of variables z = √ r 2 + t 2 in the last integral in the above right-hand side, this rewrites in turn: An application of the Fubini theorem now yields: and so, using another polar change of variables in R d−1 : which is the desired identity.
It will be useful in the following to calculate the value C(Φ) of the constant attached via (4.2) to the fundamental solution Φ of the operator v → −∆v + v defined in (3.19).
As a consequence of Lemma 3.10, the function Φ is in L 1 (R d ), and it easily follows that belongs to L 1 (R). In particular, the functions Φ and Φ are tempered distributions on R d and R, respectively. Moreover, the Fourier transform FΦ of Φ, defined by is a continuous function on R d , which is related to the one-dimensional Fourier transform FΦ of Φ via the relation: where, again, we have denoted by 0 the point (0, . .
Taking the (one-dimensional) inverse Fourier transform of the previous relation, it follows that the function Φ ∈ L 1 (R) satisfies the following ordinary differential equation in the sense of distributions on R: A direct resolution of the latter yields: With this result at hand, we are now in position to calculate the constant C(Φ): as desired.
We now prove the main result of this subsection.
Proposition 4.3. Let D ⊂ R d and Ω ⊂ D be two domains satisfying the assumptions (4.1); the following convergence holds: where C(η) is the constant defined by (4.2).
Proof. For clarity, we limit ourselves to proving the proposition in the case where Ω D is a smooth domain. The full statement of Proposition 4.3 is obtained by applying the argument below to the each individual smooth portion of ∂Ω. The definition (3.8), and (3.15) of the approximate perimeter P conv ε (χ Ω ) yields: Since the function η has compact support inside the unit ball vanishes unless y ∈ B 1 (0) and x belongs to the tubular neighborhood {x ∈ Ω, d(x, ∂Ω) < ε}. Using the coarea formula of Lemma A.3, the above expression then rewrites: where κ 1 , . . . , κ d−1 are the principal curvatures of ∂Ω, see Appendix A.1.
For given x ∈ ∂Ω and t ∈ (−1, 0), we are now led to investigate the behavior of the quantity: as ε → 0. Without loss of generality, we assume that where e d is the d th coordinate vector in R d . Since Ω is smooth around x, there exists an open neighborhood U of x in R d and a smooth function f : R d−1 → R such that Ω coincides with the subgraph of f around x, that is: see Appendix A.1. In order to calculate I ε , we approximate it by the following integral, defined for fixed x ∈ ∂Ω, t ∈ (−1, 0): whose definition mimicks that of I ε except that χ Ω is replaced by the characteristic function χ of the lower half-space {x = (x 1 , . . . , x d ), x d < 0}. The error between both quantities is estimated as: A simple calculation yields, for ε small enough and y ∈ B 1 (0), 0 otherwise, and Hence, it is easy to see that A rough estimate of (4.6) based on the Fubini theorem now implies: Noting that f (0) = 0 and ∇f (0) = 0 because of (4.5), there exists a constant C > 0 such that: and so: We are now in position to calculate the limit of I ε . Using again (4.7) and the Fubini theorem, I ε reads: where the function η is defined by (4.3). Therefore, Collecting (4.4), (4.8) and (4.9), we obtain: Eventually, an integration by parts yields the alternative expression for the constant involved in the above convergence result: which is exactly the constant C(η) defined in (4.2). This terminates the proof. In this section, we prove that, for a given domain Ω ⊂ D fulfilling (4.1), the approximate quantity P Rob ε (χ Ω ) converges to Per T (Ω), up to a factor 1 2 . This conclusion will be improved in Section 5 (see notably Thm. 5.4) by using completely different techniques. However, we believe that the following argument is interesting in itself as it sheds light about why Robin (and not e.g. Dirichlet) boundary conditions are somehow "natural" in the definition of our PDE-based approximate total perimeter functional.  Proof. Let us recall the expression of P Rob ε (χ Ω ) from Section 3.2.2: where v ε is the unique solution in H 1 (D) to the boundary-value problem The latter rewrites, under variational form: A formal integration by parts in the above definition of P Rob ε (χ Ω ) yields: A careful understanding of the behavior of v ε near the boundary ∂Ω is thus the key ingredient in passing to the limit in (4.12). As we shall see, grossly speaking, v ε behaves like the characteristic function χ Ω as ε → 0, near points x ∈ D lying "far" from ∂Ω as ε → 0. On the contrary, the description of v ε near ∂Ω will involve boundary layer correctors.
For the sake of simplicity and without loss of generality, we limit ourselves with the case of two space dimensions in the present proof. Let us recall that, according to our geometric assumption (4.1) whereby Ω arises as the intersection of D with a piecewise smooth domain O ⊂ R d , the boundary ∂Ω is decomposed as: the former region being the part of ∂Ω lying strictly inside D and the latter being the part of ∂Ω which is included in ∂D, see again Figure 2. Also, the intersection ∂O ∩ ∂D consists of finitely many points p 1 , . . . , p N , and there exist positive constants c, δ such that: (4.13) Our calculations proceed as in Chapter 2 of [49], and the proof is decomposed into four steps.
Step 1: We analyze the difference z ε := v ε − χ Ω from a formal point of view.
To achieve this, we first derive a variational characterization of z ε . This function belongs to H 1 (Ω) and H 1 (D \ Ω); it is discontinuous across Γ and its jump across this boundary reads, for x ∈ Γ: is the difference between the one-sided limits z ± ε (x) := lim t→0 ± z ε (x + tn(x)). Likewise, the jump of the normal derivative of z ε across ∂Ω reads: The function z ε additionnally satisfies: As we have mentioned, we expect z ε to be concentrated around the boundary ∂Ω and to show "fast" decay away from the latter as ε → 0. To verify this, let δ > 0 be the small thickness parameter as in the statement of Proposition A.2, and let us take arbitrary test functions w with compact support inside D ∩ ∂Ω δ in (4.16), where we recall the definition of the tubular neighborhood ∂Ω δ := x ∈ R d , |d Ω (x)| < δ . Applying the coarea formula from Lemma A.3, we obtain: In particular, we may select functions w of the form where ζ ∈ C ∞ c (Γ) and g ∈ C ∞ c (−1, 1) are arbitrary. Doing so and using the change of variables t = εu in (4.17), we obtain: Neglecting the terms weighted by ε and ε 2 in the above identity and using the fact that ζ ∈ C ∞ c (Γ) and g ∈ C ∞ c (−1, 1) are arbitrary, we see that for any given point x ∈ Γ, the function on the real line f x (u) := z ε (x + εun(x)) should satisfy: where we recall the jump relations (4.14) and (4.15). A direct resolution now yields: This supplies our candidate as for the behavior of z ε near the region Γ ⊂ ∂Ω: where d O is the signed distance function to O, see Section A.2 about this notion, and notably Proposition A.2.
Using similar considerations, we expect that: Step 2: Definition of the boundary layer correctors. Let us denote by ∂O δ := x ∈ R d , d(x, ∂O) < δ the tubular neighborhood of ∂O with width 2δ (and likewise for ∂D δ ), where δ > 0 is chosen as in the statement of Proposition A.2. Let ϕ, ψ : R → R be two smooth cut-off functions with the following properties: and Following the analysis in Step 1, we introduce profile functions r Ω and r D , accounting for the local behavior of the difference (v ε − χ Ω ) in the neighborhoods of Γ and Γ D , respectively: The corrector function r ε : D → R is eventually defined by: that is, r ε coincides with r Ω near Γ and with r D near Γ D , outside neighborhoods of the junction points of ∂O ∩ ∂D with diameter of order ε.
We aim to prove that v ε behaves likes χ Ω + r ε as ε → 0, i.e. we estimate the residual in a way yet to be specified. To achieve this, let us observe that s ε ∈ H 1 (D) since both functions v ε and r ε + χ Ω do, by construction of r ε . A simple calculation based on Lemma A.3 yields the following estimates about r Ω and r D : here and until the end of the proof C stands for a positive constant which does not depend on ε.
Moreover, it follows from assumption (4.13) that, for i = 1, . . . , N , the following estimate holds: Step 3: Estimate of the function s ε . For an arbitrary test function w ∈ H 1 (D), we calculate: where we have used the variational formulation (4.11) to pass from the first line to the second one, and we have defined: Grossly speaking, the contribution I 1 ε (w) appraises how well the corrector r Ω accounts for the harmonic behavior of (v ε − χ Ω ) inside Ω and the jump of this function across Γ. The contribution I 2 ε (w) measures how well r D represents the Robin boundary condition fulfilled by (v ε − χ Ω ) on Γ D .
We proceed to estimate each term separately. Throughout the proof, we denote by q ε (w) a remainder term, possibly changing from one line to the next, which satisfies the estimate: At first, the definition of I 1 ε (w) rewrites: Using (4.20), the first integral in the above right-hand side is estimated as: Likewise, relying on the estimate (4.21), we prove that: As far as the remaining two terms in (4.23) are concerned, an integration by parts yields: From the definition (4.18) of r Ω , the first integral in the above right-hand side vanishes. Moreover, using again (4.20) and (4.21), the second and third integrals are estimated as: We now consider the last integral in the right-hand side of (4.24). Recalling the expression (A.3) of the Laplace operator in the local basis (τ, n) attached to the tubular neighborhood ∂O δ a simple calculation based on the explicit expression (4.18) of r Ω and the properties of the signed distance function recalled in Appendix A.2 yields: (4.26) Hence, we obtain: We now estimate the term I 2 ε (w) for an arbitrary function w ∈ H 1 (D). This quantity rewrites: (4.29) The first integral in the above right-hand side is easily estimated thanks to (4.20): (4.30) As for the second integral in the right-hand side of (4.29), we get similarly: We now proceed to estimate where we have used the coarea formula of Lemma A.3. Taking advantage of (4.13), there exists a constant c > 0 such that: and so: This entails: and so, recalling (4.31): When it comes to the third term in the right-hand side of (4.29), an integration by parts, followed by estimates similar to those previously derived yields: Gathering (4.29), (4.30), (4.33), and (4.34), we obtain: Eventually, explicit calculations similar to (4.26), based on the definition (4.18) of r D yield: Using the estimates (4.28) and (4.35) in (4.22), then taking w = s ε as test function, it follows that: We now apply the Cauchy-Schwarz inequality with the following inner product over R 3 to obtain: 1/2 ε 2 + εε + ε 2 1/2 , and so: Step 4: Asymptotic analysis of the approximate perimeter P Rob ε (χ Ω ). Let us recall the alternative formula (4.12) for P Rob ε (χ Ω ): A straightforward calculation based on the explicit expression (4.18), and (4.19) of the corrector r ε reveals that: Therefore, the conclusion of the theorem follows from the convergence: which we now proceed to show. To this end, let σ ε := ε∇s ε ∈ L 2 (D) 2 ; the following estimate is a direct consequence of (4.37): We now prove that ||divσ ε || L 2 (Ω) ≤ C. To achieve this, let us calculate: in the sense of distributions in Ω. Therefore: The first term in the above right-hand side is bounded in L 2 (Ω) owing to (4.37), and the second one is also bounded in L 2 (Ω) as is immediately revealed by an explicit calculation based on the expression (4.18), and (4.19) of the corrector r ε . Hence, σ ε is a bounded sequence in H div (Ω), and so there exists a subsequence (still indexed by ε for simplicity) and σ 0 ∈ H div (Ω) such that σ ε ε→0 − −− → σ 0 weakly in L 2 (Ω), and div(σ ε ) ε→0 − −− → div(σ 0 ) weakly in L 2 (Ω).
Because of (4.39), the weak limit σ 0 is necessarily 0, so that along the same subsequence: By a classical argument relying on the uniqueness of the limit, the above convergence actually holds for all the sequence ε (and not only along a subsequence). This shows (4.38), and terminates the proof of the theorem.

4.3.
Convergence of the derivative of the approximate perimeter functionals P conv ε (u) and P Rob ε (u) As we have mentioned, we intend to use P conv ε (u) and P Rob ε (u) in topology optimization algorithms, as generalized counterparts of the usual total perimeter Per T (Ω) of shapes Ω. Recalling the common structure (3.8) of these functionals, where L ε stands for either L conv ε in (3.15) or L Rob ε in (3.16), we have just proved the consistency of the values of P ε (u) with those of Per T (Ω): when u = χ Ω is the characteristic function of a shape Ω ⊂ D (satisfying mild assumptions), the quantity P ε (χ Ω ) converges to C(µ)Per T (Ω) as ε → 0; see (3.10) and (4.2) for the definition of the constant C(µ).
As we shall see in more details in Section 6.4.1, a key quantity in the numerical treatment of density-based topology optimization problems involving P ε (u) is the derivative of this functional. Its expression follows from a standard calculation: where we have used the fact that the operators L conv ε and L Rob ε are self-adjoint with respect to the L 2 (R d ) and L 2 (D) inner product, respectively.
In this section, we examine the consistency of the derivative P ε (u) of our approximate perimeter functionals with the shape derivative Per T (Ω)(θ) (supplied by Prop. 2.2) of the exact total perimeter in the particular case where u = χ Ω is the characteristic function of a given domain Ω. Our main result reads as follows.
Proposition 4.8. Let D ⊂ R d and Ω ⊂ D be two bounded domains satisfying the assumption (4.1), and let x ∈ ∂Ω ∩ D be a point such that the boundary ∂Ω is smooth around x. Let where L ε is one of the regularizing operators (3.9) and (3.10). Then, the following convergence holds: where κ(x) is the mean curvature of ∂Ω at x and C(µ) is the constant associated with the radial principal part µ of L ε via (4.2).
Remark 4.9. This result expresses the quite intriguing fact that when u = χ Ω is the characteristic function of a smooth subdomain Ω ⊂ D satisfying (4.1), the derivative of the approximate perimeter P ε (u), which is a quantity able to measure "topological changes" in the shape accounted for by the density u, is consistent with the shape derivative of the exact perimeter functional, describing geometric variations of Ω.
Proof. Without loss of generality, we assume that the considered point x coincides with the origin 0, and that the normal vector to ∂Ω at x is n(x) = e d . A simple calculation, based on the structure (3.9) of L ε , yields: It follows from the decay estimate (3.13) that the second term in the above right-hand side vanishes as ε → 0. Then, a simple calculation based on the properties of the function µ in (3.10) yields: so that: Let us now introduce a local representation of Ω around x. As in Appendix A.1, there exists a real number ρ > 0 and a smooth function f ∈ C ∞ (R d−1 ) (which we assume to be defined on the whole space R d−1 for commodity) such that: Our convention whereby x = 0 and the normal vector n(0) equals e d implies that: where II 0 is the second fundamental form of ∂Ω at 0. Hence, for ε small enough and y ∈ B ρ ε (0), we obtain: Now taking advantage of the decay (3.12) of µ at infinity, it follows: and (4.41) then becomes: where we have introduced the smooth function b(ε) := − 1 ε f (−εy 1 , . . . , −εy d−1 ), whose dependence on y 1 , . . . , y d−1 is omitted for notational simplicity. In order to pass to the limit ε → 0 in the above expression, we first need to calculate the values b(0) and b (0); to achieve this, a Taylor expansion at ε = 0 yields: and using (4.42), we obtain: We then infer that b(0) = 0 and b (0) = 1 2 II 0 ( y, y). Now using the smoothness of the partial mapping y d → µ(y 1 , . . . , y d−1 , y d ) for fixed y = (y 1 , . . . , y d−1 ) = 0, we obtain for y = 0: (y 1 , . . . , y d−1 , b(0)).
The estimates (3.11), and (3.12) about µ then make it possible to apply the Lebesgue dominated convergence theorem in (4.43), which results in: where we have denoted by II 0,i,j the entries of the (d − 1) × (d − 1) matrix II 0 , and we have used the radial character of µ, which implies that and that all the integrals Remark 4.10. The above analysis stands only in the neighborhood of the points x ∈ ∂Ω ∩ D where the boundary ∂Ω is smooth, and it seems difficult to extend to intersection points of the boundaries ∂Ω and ∂D. The behavior of the approximate perimeter functional P Rob ε (u) near those intersection points as ε → 0 will be better understood in the framework of Γ-convergence, as we shall see in the next Section 5.

Convergence properties of the PDE-based approximate
perimeter functional P Rob ε (u) In this section, we look at the approximation of the total perimeter P T (u) in (3.7) from another perspective, namely, we investigate the Γ-convergence of the PDE-based functional P Rob ε (u) to its exact counterpart P T (u), which indicates how the minimum value and the minimizers of P Rob ε (u) behave as ε → 0. We first recall in the next Section 5.1 some material about the notion of Γ-convergence, and we state Theorem 5.4, which is the main result in this section. Technical preliminaries for the proof of this theorem are discussed in Sections 5.2, 5.3; Sections 5.4, 5.5 are then devoted to the proof of the Γ-limsup (resp. Γ-liminf) inequalities. Finally, we deal with the equi-coercivity of P Rob ε (u) in Section 5.6.

A brief reminder of Γ-convergence
For the convenience of the reader, we summarize in this section the basic ingredients of the Γ-convergence theory, referring to Chapter 12 of [22] and [29] for further details.
Let us start with a definition: Definition 5.1. Let (X, d) be a metric space, F ε : X → R ∪ {∞} be a sequence of functions on X, and F : X → R ∪ {∞} be another such function. The sequence F ε is said to Γ-converge to F if for every x ∈ X, the following two conditions hold at every point x ∈ X: -(Γ-limsup inequality) There exists a sequence x ε ∈ X with d(x ε , x) → 0 such that: The relevance of the notion of Γ-convergence comes from the next theorem, according to which if (quasi-) minimizers x ε of a Γ-convergent sequence of functions F ε accumulate at some point x ∈ X, then x is a minimizer of the Γ-limit F . Theorem 5.2. Let (X, d) be a metric space, and let F ε : X → R ∪ {∞} be a sequence of functions which Γ-converges to F : X → R ∪ {∞}. Then, (i) Let x ε ∈ X be a sequence of quasi-minimizers of F ε , i.e. there exists a sequence of positive real numbers λ ε converging to 0 such that: If the set {x ε } is relatively compact, then every of its cluster points is a minimizer of F over X, and: (ii) If G : X → R ∪ {∞} is a continuous function, then (G + F ε ) Γ-converges to (G + F ).
An attached concept with utmost interest is that of equi-coercivity, which is precisely about the compactness of sequences of (quasi-)minimizers of F ε assumed in the previous theorem. Definition 5.3. Let (X, d) be a metric space; a sequence of functions F ε : X → R ∪ {∞} is said to be equicoercive if from every sequence x ε ∈ X such that sup ε F ε (x ε ) < ∞, one can extract a subsequence ε j such that x εj converges to some element x * ∈ X.
With these notions at hand, we are in position to state the main result in this section:    (3.20) to the relative perimeter functional P R (u). As we shall see, the present analysis of the functional P Rob ε (u) involves much more technicalities.
Remark 5.6. Theorem 5.2 and 5.4 exemplify why the notion of Γ-convergence is well-adapted to our original purpose of approximating the total perimeter P T (u) by P Rob ε (u) in the context of a density-based optimal design problem. Let us indeed consider the problem introduced in Section 2. The Γ-convergence of P Rob ε (u) to 1 2 P T (u) then ensures that minimizers of (5.2) converge to minimizers of (5.1), under suitable assumptions on J(u); see [21] for more details and examples about this point, in the context of the relative perimeter functional.

ε (u)
Our first task is to give a variational characterization of the approximate total perimeter P Rob ε ; this simply relies on the variational formulation associated to (3.17) and on a duality argument. Throughout this section, for a given function u ∈ L ∞ (D, [0, 1]), we denote by v ε the unique solution in H 1 (D) to:

4)
where the infimum is uniquely attained at v = v ε . Alternatively, P Rob ε (u) rewrites as a supremum:
Proof. The variational formulation for (5.3) reads: As a consequence of the standard Lax-Milgram theory, the following relation holds: where v ε is the unique minimum point in the above right-hand side. Combining this expression with the definition (3.18) of P Rob ε (u), we immediately obtain (5.4). In order to obtain (5.5), we use a duality argument, which is more systematically described in [47]. The latter relies on a simple algebraic identity, taking place in a Hilbert space (H, ·, · H ), involving an invertible, positive operator A : H → H: where the supremum is uniquely attained at σ = Av. Thence, it easily follows from the previously established identity (5.4) that: where In addition, the min-max problem (5.6) has exactly one saddle point, namely (v, σ, w, z) = (v ε , ε∇v ε , 1 ε v ε , v ε ); hence, the min and max in (5.6) can be interchanged. We then arrive at: The condition (5.8) imposes in turn that −divσ = 1 ε u − w, and σ · n = −z.

Study of an auxiliary function
In this section, we consider a given function u ∈ L ∞ (D, [0, 1]), as well as its extension u by 0 to R d . In the forthcoming derivations, we shall encounter on several occasions the function w ε := u * Φ ε , where we recall that see (3.19). By definition, w ε is the unique solution in H 1 (R d ) to the following equation: In this section, we aim to analyze the precise asymptotic behavior of w ε as ε → 0, and we start with a few remarks. (i) For all ε > 0, the function w ε belongs to H 2 (R d ).
(ii) The sequence w ε converges to u strongly in L 1 (R d ).
(iii) The following estimates hold: Proof. Proof of (i): Since u ∈ L 2 (R d ), this point is a simple consequence of the classical elliptic regularity theory, see for instance ( [33], Sect. 9.6).
Proof of (ii): Since Φ has unit integral (see Lem. 3.10), it holds for almost every point x ∈ R d , and so For an arbitrary real number ρ ≥ 0, we may estimate: and it follows from the decay of Φ at infinity (see again Lem. 3.10) that, for given δ > 0, there exists ρ > 0 such that for every ε > 0: Therefore, using the Cauchy-Schwarz inequality and considering the behavior of Φ in the neighborhood of 0 recalled in Lemma 3.10, it follows: for a constant C > 0 which may depend on δ and ρ, but which is independent of ε. Since u takes values in [0, 1], this entails: It follows from the Lebesgue differentiation theorem and the Lebesgue dominated convergence theorem that the integral in the above right-hand side converges to 0 as ε → 0. In particular, for ε sufficiently small, it is smaller than δ, and then: Since δ is arbitrary, this reveals that w ε → u strongly in L 1 (R d ).
Proof of (iii): This fact follows from quite precise estimates about the fundamental solution Φ ε . Since the proof has already been tackled in the 2d case in [20] (see Prop. 2.7 in there), and since it is similar, albeit more involved in the 3d case, we only focus on the latter, which is new to the best of our knowledge. To this end, let us first recall the explicit expressions of the modified Bessel function of the second kind K 1 2 (r) and the fundamental solution Φ(x) in 3d (see [1], Sect. 9.6.23): A simple change of variables in the definition of w ε now yields: and the first part of the statement then follows from the simple calculation: where we have introduced the spherical coordinates (r, θ, ϕ) centered at 0 in passing from the second line to the third one. Now let x ∈ R 3 be given, and let ν ∈ R 3 be an arbitrary unit vector; we obtain from the definition of w ε and a change of variables that: Since the function Φ(z) is radial, the chain rule allows to derive the following relation, which is valid for all point y ∈ R 3 \ {0} and any unitary matrix R ∈ R 3×3 such that Rν coincides with the first coordinate vector e 1 of R 3 : Using another change of variables, the quantity ε∇w ε (x) · ν may now be estimated as follows: Hence, we have proved that: which leads to the desired conclusion.
It is quite natural to try and approximate the solution v ε to our boundary-value problem (5.3) by the function w ε . Indeed, the latter satisfies −ε 2 ∆w ε + w ε = u inside D. The purpose of the next lemma is to show that, as ε → 0, w ε is "close" to fulfilling the Robin boundary condition of v ε in (5.3): ε ∂wε ∂n + w ε ≈ 0 as ε → 0. This fact may be understood as a geometric measure theory counterpart of the derivation of the boundary layer profile r D in the proof of Theorem 4.5.
Let us start with a lemma.
Lemma 5.9. Let f ∈ L ∞ (R d ), and let x 0 ∈ R d be a regular point of f . Let a ∈ S d−1 be a unit vector such that both approximate limits f a (x 0 ) and f −a (x 0 ) exist, with f a (x 0 ) = 0. Then the function w ε = f * Φ ε ∈ H 1 (R d ) satisfies: Remark 5.10. Despite its quite technical proof, the above result is fairly intuitive: roughly speaking, if the function f jumps from some value α = f −a (x) to 0 across the hyperplane with normal vector a, the smoothed approximation w ε should take a value close to 1 2 α on this hyperplane, and it should have a normal derivative behaving as − α 2ε . Proof. The definition of w ε immediately implies that: We now fix ρ > 0 and δ > 0. Then, where the remainder r 1 is defined by:

Denoting by
α := f a (x 0 ) = 0, and α := f −a (x 0 ), we further decompose: where we have set (5.13) and Let us now introduce the functions r(y) = f (y) − α and r (y) = f (y) − α . Then, where (5.14) and Finally, since α = 0, we arrive at the following expression: We now proceed to estimate each term in the right-hand side of (5.15).
Then a change of variables and the Fubini theorem yield: Ha,t ∇Φ(y) · a + Φ(y) dy dt. Because Φ(z) depends only on the radial part |z| of z, the function ψ a is independent of the vector a. In particular, ψ a (t) = ψ e d (t), where e d is the d th coordinate vector.
On a different note, we may calculate the derivative of ψ a by considering the expression Changing variables in the first integral in the above right-hand side and using the Lebesgue dominated convergence theorem, we obtain: Based on these ingredients, (5.17) rewrites: As in the statement of Lemma 4.1, for a given point x = (x 1 , . . . , x d ), we use the shorthand x = (x 1 , . . . , x d−1 ) ∈ R d−1 . We also denote by F x (resp. F t ) the partial Fourier transform in R d−1 (resp. in R) with respect to the d − 1 first variables (resp. to the last one). Since ψ a (t) = ψ e d (t), it follows that Therefore, taking the partial Fourier transform F t in both sides, we obtain: On the other hand, since −∆Φ + Φ = δ 0 in the sense of distributions in R d , the Fourier transform of Φ reads FΦ( ξ, ω) = 1 1+| ξ| 2 +ω 2 , and so F t ψ a (ω) = 1 1 + ω 2 .
The Fourier inversion formula then yields: In particular, for t > 0, we have ψ a (t) + ψ a (t) = 0, which implies that: We now proceed to estimate the remainders r 1 , r 2 , r 2 , r 3 , r 3 , r 4 featured in the formula (5.15) in terms of δ. This task demands a careful choice of the radius ρ of the ball B ρ (x 0 ) in their definitions and we start by introducing two functions which prove useful in this perspective.
Estimate of the remainder r 1 . A change of variables in the definition (5.12) yields: Since f ∈ L ∞ (R d ), ζ(ε) ε → ∞ and because of the exponential decay of the fundamental solution Φ(r) as r → ∞ (see Lem. 3.10), it follows that r 1 → 0 as ε → 0.
Estimate of the remainder r 2 . Likewise, a change of variables in (5.13) yields: Considering the behavior of Φ and its derivative near 0 and ∞ (see Lem. 3.10) and the fact that f ∈ L ∞ (R d ), a simple application of Hölder's inequality yields: for some q > 2. Here and until the rest of the proof, C stands for a constant which depends on f but not on δ and ε. We now estimate this last integral as: where we have used the definition (5.18) of θ from the second to the third line, as well as (5.19) to achieve the last line. Using again (5.19), we see that r 2 → 0, as desired.
Estimate of the remainder r 2 . One obtains r 2 → 0 as ε → 0 by using exactly the same arguments as in the case of r 2 .
Estimate of the remainder r 3 . From the definition (5.14), we obviously have: and so, using a change of variables together with the fact that Φ and ∇Φ are integrable on R d (see Lem. 3.10): Estimate of the remainder r 3 . By the same token, there exists a constant C > 0 which does not depend on δ and ε such that: Estimate of the remainder r 4 . A change of variables in the definition (5.16) leads to: We then decompose r 4 as: and I 2 := R d ∇Φ(z) · a + Φ(z) (χ π−a(x0)∩Bρ(x0) − χ π−a(x0)∩Bρ(x0)∩{|r |<δ} )(x 0 − εz) dz. Now, on the one hand, using again Hölder's inequality, we see that there exists q > 2 such that: it follows from (5.19) that I 2 → 0 as ε → 0.
On the other hand, which converges to 0 as ε → 0 for a.e. z ∈ R d since ζ(ε) ε → ∞. This, together with the integrability of Φ and ∇Φ on R d and the Lebesgue dominated convergence theorem, imply that I 1 → 0. Hence, r 4 → 0 as ε → 0.
Gathering all these ingredients, we have finally proved that: since this holds fo arbitrary δ > 0, the desired result is proved.
We now provide the missing link in the previous proof. is lower semi-continuous.
Proof. Let ρ > 0 be arbitrary, and let ρ n be any sequence of elements in [0, ∞) such that lim n→∞ ρ n = ρ. For any given point r 0 < ρ and for n large enough, we have r 0 ≤ ρ n , and so Taking the lim inf, we thus obtain: and recalling that the sequence ρ n → ρ is arbitrary, this results in Since 0 ≤ r 0 ≤ ρ is arbitrary, and since L is continuous, we actually end up with: that is: whence the claim in the case where ρ > 0. The case ρ = 0 is handled in a similar way, which completes the proof.
Corollary 5.12. Let u ∈ BV (D, {0, 1}) and set w ε = u * Φ ε ∈ H 2 (D). Then, In particular, the following limit holds: Proof. Since the function u is in BV (R d , {0, 1}), there exists a set Ω ⊂ D with finite perimeter such that u = χ Ω . We know from Proposition 3.8 that H d−1 almost every point x 0 ∈ R d is a regular point of u. On a different note, since D is a Lipschitz domain, the unit normal vector n to ∂D exists H d−1 -almost everywhere on ∂D.
Let then x 0 ∈ ∂D be a point where the latter two conditions hold and let a ∈ S d−1 be a direction such that the approximate limits f a (x) and f −a (x) exist. Two cases occur about x 0 , up to an H d−1 negligible set: -If x 0 belongs to the measure theoretical boundary ∂ M Ω of Ω, then x 0 ∈ ∂ M Ω \ D, and so we may take a = n, the unit normal vector to ∂D, see again Proposition 3.8. Obviously, the approximate limit u a (x 0 ) equals 0, and Lemma 5.9 implies that ε ∂wε ∂n + w ε (x 0 ) → 0.
-If x 0 does not belong to ∂ M Ω, it does not belong to the jump set of χ Ω , and in particular, the approximate limit u a (x 0 ) exists for all a ∈ S d−1 . Obviously, for a = n, it equals 0, and so, again ε ∂wε ∂n + w ε (x 0 ) → 0. This proves (5.20).
The second statement of the lemma follows directly from this, the uniform bounds on w ε in Lemma 5.8 and the Lebesgue dominated convergence theorem.

Proof of the Γ-limsup inequality
The following proposition implies that the Γ-limsup inequality in Denition 5.1 is satisfied by the sequence P Rob ε (u), which is one part of the statement of Theorem 5.4. Proof. Let us first consider a function u ∈ H 1 (D, [0, 1]) ⊂ BV (D, [0, 1]). As in the previous developments, we let u be the extension of u by 0 outside D, and w ε = Φ ε * u. The expression (5.4) of P Rob ε (u) as an infimum immediately implies that: Performing integration by parts in the right-hand side of the above identity and using the fact that −ε 2 ∆w ε + w ε = u in D, we obtain: Using the fact that u ∈ H 1 (D), the identity −ε 2 ∆w ε + w ε = u in D and integration by parts in the first integral in the above right-hand side, it follows that: We now rely on the estimates of the function w ε and its gradient contained in Lemma 5.8 above to find: We now wish to prove that the above inequality actually holds when u ∈ BV (D, [0, 1]), and to this end, we rely on a density argument. Recalling the density of Hence, we may apply the inequality (5.21) with u replaced by v n , leading to: We now wish to pass to the limit n → 0 in the above identity. This is possible thanks to (5.22) and the following facts: -The mapping u → P Rob ε (u) is continuous from L 1 (D) into R; -Since v n and u take values in [0, 1], the following convergence holds: Hence, (5.21) indeed holds when u ∈ BV (D, [0, 1]). Eventually, when u ∈ BV (D, {0, 1}), this yields: Now using Corollary 5.12, we arrive at: which is the desired result.

Proof of the Γ-liminf inequality
We eventually turn to the proof of the Γ-liminf inequality for the functional P Rob ε (u). The result of interest reads as follows.
The proof hinges on a localization and slicing argument, and we carefully follow the trail exposed in Chapter 15 of [29]. To this end, we first introduce a "localized" version of the approximate perimeter functional P Rob ε (u) which emphasizes the spatial dependence of the defining integrals. For any open subset A ⊂ R d and u ∈ L ∞ (D, [0, 1]), let us define the quantity F ε (A, u) by: Our aim is to prove that: where we denote as usual by u the extension of u ∈ L ∞ (D, [0, 1]) by 0 outside D. In the above, F (A, u) is defined as the Γ-liminf of the sequence of functions u → F ε (A, u), that is, Note that it is easily verified from the definition that for a given u ∈ L ∞ (D, We then express G ε (A, u, v) in terms of one-dimensional integrals for an arbitrary open subset A ⊂ R d . To this end, for an arbitrary direction ξ ∈ S d−1 , we denote by π ξ the linear hyperplane orthogonal to ξ: π ξ := z ∈ R d , ξ · z = 0 , and for y ∈ π ξ , we introduce the two subsets of R: D ξ,y = {t ∈ R, y + tξ ∈ D} , A ξ,y = {t ∈ R, y + tξ ∈ A} , whose characteristic functions satisfy: χ D ξ,y (t) = χ D (y + tξ), and χ A ξ,y (t) = χ A (y + tξ).
For given u ∈ L ∞ (D, [0, 1]) and v ∈ H 1 (D), we decompose G ε (A, u, v) along the normal lines to π ξ : In order to deal with the first integral in the defintion (5.24) of G ε (A, u, v), we have used the Fubini theorem; the treatment of the second integral relies on the following technical lemma (with the remark that t ∈ ∂D ξ,y implies that (y + tξ) ∈ ∂D), whose proof is postponed to the end of the present section.
Recalling from the definition that u ξ,y = 0 on R \ (a ξ,y , b ξ,y ), the above calculations allow to infer the following lower bound of H ε (A, u, v, y): Grossly speaking, we have here added a "thin skin" near the points of ∂D belonging to A.
Introducing the C 1 function ψ : R → R defined by: it follows: where we have used the results in Section 3.1.2 of [46] to calculate the derivative of the composite function ψ • v ε,ξ,y in the context where ψ is of class C 1 and v ε,ξ,y is in H 1 ( A ξ,y ). Let us set w ε,ξ,y := ψ • v ε,ξ,y and w ξ,y := ψ • u ξ,y ; introducing the total variation |Dw ε,ξ,y | (see Sect. 3.1), this rewrites: The conclusion (5.32) of the previous step immediately implies that w ε,ξ,y → w ξ,y strongly in L 1 ( A ξ,y ) as ε → 0. The lower semi-continuity of the total variation of functions in L 1 ( A ξ,y ) implies that: Note that the right-hand side of the above inequality may be infinite; the meaning of this inequality is that the liminf in the left-hand side is also infinite in such a case. Since u ξ,y ∈ L ∞ ( A ξ,y , {0, 1}), we eventually observe that As a result: where we have used the fact that |Du ξ,y |(A ξ,y \ A ξ,y ) = 0 to obtain the last equality.
Step 4: We recover the multi-dimensional case from the one-dimensional result. To this end, let us recall from the definitions that: Using Fatou's lemma yields: Hence, using the result of Step 3: We now rely on the following theorem about the characterization of functions with bounded variations by their one-dimensional sections, see Remark 3.104 and Theorem 3.107 in [12]. . . , ξ d ∈ R d such that u ξi,y ∈ BV (A ξi,y ) for dH d−1 a.e. y ∈ π ξi , and ∀i = 1, . . . , d, Then, for any direction ξ ∈ S d−1 , Now, if u is not a function with bounded variations in A, the integrand in the right-hand side of (5.33) takes an infinite value for a set of y ∈ π ξ with positive H d−1 Hausdorff measure, and we conclude that F (A, u) is infinite.
We finally suppose that, on the contrary, u ∈ BV (A, {0, 1}), that is, the restriction of u to A is the characteristic function of a subset E ⊂ A with finite perimeter. We then obtain: and since ξ is arbitrary: in the sense that if the right-hand side in the above inequality is infinite, so is the left-hand side. It follows that the Γ liminf function F (A, u) satisfies: Using Proposition 3.8 for the expression D u = νH d−1 | ∂ M E , where ν(x) is the inner measure theoretic normal vector to E, taking values in the unit sphere S d−1 , it follows: Since the mapping A → F (A, u) is a super-additive function, it then follows from Lemma 15.2 in [29] that the supremum and the integral sign may be interchanged in the above identity, so that: Eventually, taking A = R d yields the desired result.
We now provide the proof of the technical Lemma 5.15.
Proof of Lemma 5.15. Since ∂D is a compact and Lipschitz hypersurface in R d , it is enough to prove that, for any non negative function g ∈ L 1 (S), the following relation holds In order to prove (5.34), let p ξ : R d → π ξ be the orthogonal projection onto π ξ ; we introduce the Lipschitz function φ = p ξ • f : U → π ξ as well as g = g • f ∈ L 1 (U ). Using the Lipschitz change of variables formula (see e.g. [40], Sect. 3.3.3), we obtain: where the Jacobian determinant J(φ) is defined by J(φ)(x) = det(∇φ(x) T ∇φ(x)), for a.e. x ∈ U.
Since φ = p ξ • f , it follows ∇φ(x) = ∇p ξ ∇f (x); now using that ∇p ξ is the matrix of an orthogonal projection, we obtain where we have used the following elementary identity, valid for any orthogonal projection matrix P ∈ R d×d , and any M ∈ R d×(d−1) : Hence, we obtain from the non negativity of g that and so, using again the change of variables formula to transform the left-hand side of the above inequality, we arrive at which is the desired result.
This result of Proposition 5.14 applied with the constant sequence u ε = u, together with Proposition 5.13 immediately imply the following pointwise convergence of the functional P Rob    these functions show a transition between the values 0 and 1 over a sharper and sharper interface as ε → 0, which represents a blurred version of ∂Ω ∩ D in the former case, and of ∂Ω in the latter.

Numerical validation of the derivative of the perimeter functional
We now turn to illustrate the consistency of the derivative P Rob ε (u) of the approximate perimeter functional P Rob ε (u) with respect to the shape derivative of the exact function Per T (Ω). More precisely, let us recall the formula for P Rob ε (u) obtained in Section 4.3: In this section, the computational domain D is the unit two-dimensional square D = (0, 1) 2 , and Ω ⊂ D is the disk Ω = x = (x 1 , x 2 ) ∈ D, x 2 1 + x 2 2 ≤ R 2 with R = 0.25, so that Ω does not intersect the boundary ∂D, see Figure 7.
Letting u := χ Ω , we aim to study the "derivative" g u,ε , and notably to verify the conclusion of Proposition 4.8 whereby is not explicitly discretized in the mesh of D, but it is represented as the negative subdomain of a "level set function".
A rather fine triangular mesh T of D is used in this experiment, with size h = 1/160. Contrary to the study conducted in Section 6.1, in this experiment as in all the forthcoming ones, the shape Ω under scrutiny is not explicitly discretized in T . Rather, Ω is defined as the negative subdomain of a so-called "level set function" ψ : D → R, that is, a function such that In practice, ψ is discretized at the vertices of the mesh T , and the boundary-value problem (6.2) for v Rob ε is solved by applying the finite element method with the mesh T of D, with a careful assembly of the right-hand side in its finite element resolution. We refer to e.g. [60] about the general idea of handling implicit geometries in numerical computations.
6.2.1. Behavior of the field g u,ε inside Ω and D \Ω The graph of the derivative g u,ε on D is represented in Figure 8 for different values of the ratio ε/h, and its values on the line L := {x = (x 1 , x 2 ) ∈ D, 0.5 ≤ x 1 ≤ 1, x 2 = 0.5} are depicted on Figure 9.
We observe that the derivative takes negative values in Ω, positive values in D \Ω, and that it modulus blows up as the ratio ε/h decreases.  6.2.2. Behavior of the field g u,ε on ∂Ω We now focus on points x ∈ ∂Ω and we appraise the convergence of g u,ε (x) to (half) the mean curvature κ(x) of ∂Ω as ε → 0, recalling that in the present geometric situation κ(x) = 1/R for all x ∈ ∂Ω.
We evaluate this behavior in terms of the mean value m ε and the variance σ ε of g u,ε over ∂Ω: m ε := 2 |∂Ω| ∂Ω g u,ε ds and σ ε : where, again, the factor 2 is introduced out of consistency with the convergence result of Proposition 4.8. In practice, g u,ε (x) is calculated at the vertices of the mesh T by a direct application of the formula (6.3). The calculation of integral quantities involving the values of g u,ε on ∂Ω, as those in (6.4), then relies on linear interpolation from these values. Both quantities m ε and σ ε are represented in Figure 10 for several values of the ratio ε/h. As expected, the mean value m ε is approaching values close to the constant mean curvature κ of ∂Ω as ε/h decreases with a small bias. Concurrently, the variance σ ε increases, which can be explained by the fact that the values of g u,ε tend to −∞ (resp. +∞) at vertices of the mesh lying inside Ω (resp. inside D \ Ω) in this regime, and so large numerical errors occur when calculating the values of g u,ε on ∂Ω by interpolation. Roughly speaking, ε/h ≈ 4 seems to strike a fair compromise between a fair accuracy of the approximation of κ by g u,ε and a reasonable level of oscillations of the latter quantity. An interesting conclusion of this experiment is that the mean curvature of a shape Ω, which is in general difficult to compute (see [41] for a discussion about this point), may be accurately approximated by the quantity 2g u,ε .

Minimization of the perimeter under a volume constraint
In this section and the next, we appraise the numerical behavior of the approximate total perimeter functional P Rob where P ε (u) stands for either the approximate total perimeter P Rob ε (u) or the approximate relative perimeter P Neu ε (u), and with a small abuse of notation, we have denoted by Vol(u) = D u dx the straightforward extension of the volume functional to density functions.
In principle, the numerical resolution of (6.6) could be conducted thanks to a standard algorithm, such as a projected gradient method. In this context, even though the optimized function u is a density, a priori taking values between 0 and 1, the Γ-convergence result Theorem 5.4 (which in itself expresses a penalization of intermediate densities) suggests that, as ε gets small, u will present less and less such intermediate regions.
This is indeed observed, and for brevity, we do not report on these simulations, referring to [21] for analogous calculations and conclusions.
Rather, we treat the problem (6.6) with a fixed-point algorithm from the previous work [18] of the first author, which only features "black-and-white" designs. This method is based on the first-order necessary optimality conditions associated to (6.6). Recalling the expression (6.3) of the derivative of P ε (u), these optimality conditions reads: if u ∈ L ∞ (D, [0, 1]) is a local solution to (6.6), Vol(u) = V T , and ∃λ ∈ R, s.t. (6.7) If we restrict ourselves to considering characteristic functions of shapes u = χ Ω ∈ L ∞ (D, {0, 1}), we are thus led to search for a shape Ω ⊂ D such that: To carry out this program, we represent the shape Ω as the negative subdomain of a "level set" function ψ : D → R: x ∈ D. (6.9) The first-order necessary conditions 6.8 are thus equivalent to the fact that the function (g χΩ,ε + λ) be one level set function for Ω. They can be rewritten as: ψ L 2 (D) = 1 and ψ = g χΩ,ε + λ g χΩ,ε + λ L 2 (D) with Ω = {x ∈ D, ψ(x) < 0}. (6.10) Note that, in the above, since ψ and cψ define the same shape Ω via (6.9) for any constant c > 0, there is no loss of generality in imposing that ψ have unit norm. We rely on a fixed point algorithm with relaxation to enforce the optimality conditions (6.10). The latter is summarized in Algorithm 1. Briefly, a sequence of level set functions ψ n , and associated shapes Ω n := {x ∈ D, ψ n (x) < 0} is produced, which are updated via a linear interpolation scheme on the unit sphere [62] (see the definition of the coefficients a n and b n in Algorithm 1). The value λ n of the Lagrange multiplier ensuring that the volume constraint is satisfied is found by bisection at each iteration of the process.
In practice, as was noted in the previous work [18], solving the problem (6.6) with Algorithm 1 for a small value of the ratio ε/h may precipitate the optimization path in a local minimum with bad performance. To circumvent this issue, we rely on a continuation strategy: starting from a fairly "large" value of ε/h, we solve a sequence of problems (6.6), associated to decreasing values of this ratio. Each sub-problem is deemed to reach convergence at the iteration n where the angle θ n between the current level set function ψ n and the derivative g n satisfies θ n < 1 • .
We solve the problem (6.6) in both situations where P ε (u) is the approximate relative perimeter P Neu ε (u), or the approximate total perimeter P Rob ε (u). The considered 2d physical setting is the following: the computational domain D is the unit square (0, 1) 2 , which is discretized with a Cartesian mesh Q composed of 10, 000 Q 1 elements. In both cases, the volume target is set to V T = 0.4|D|; the initial guess Ω 0 is the total computational domain D minus a square hole so that Ω 0 fulfills the volume constraint. Figure 12. Evolution of the quantities involved in the (upper row) relative and (lower row) total perimeter minimization problem under volume constraint of Section 6.3. The values of the perimeter (normalized with that of the initial design), volume and Lagrange multiplier are represented on the left column; the exact relative and total perimeters Per R (Ω) or Per T (Ω) of the shape Ω are compared with the approximate values P Neu ε (χ Ω ) and P Rob ε (χ Ω ) through the iterations in the central column, and the algorithm parameters θ and ε/h are represented on the right column. the (approximate) relative perimeter. On the contrary, when the total perimeter is minimized, the optimized shape is a disk which does not touch the boundary ∂D, in agreement with the classical isoperimetric inequality.
The histories of the perimeter, volume and Lagrange multiplier associated to both experiments are reported in Figure 12. In particular, we observe that the approximate perimeter functionals P Neu ε (χ Ω ) and P Rob ε (χ Ω ) converge to their exact counterparts as the ratio ε/h decreases. The images of the second column illustrate that minimizing the sequence of approximate perimeters P Neu ε (χ Ω ) and P Rob ε (χ Ω ) for a decreasing value of the ratio ε/h leads to a minimization of the exact perimeters 1 2 P R (Ω) and 1 2 P T (Ω).

Using the perimeter functional in structural optimization
We finally turn to the use of our perimeter approximations in the more realistic setting of shape and topology optimization of elastic structures.

Shape and topology optimization of elastic structures
The optimized shape Ω ⊂ D stands for a linearly elastic structure, whose boundary is decomposed into three disjoint regions Γ D , Γ N , Γ: The shape is attached on the part Γ D of its boundary and surface loads g ∈ H 1 (R d ) d are applied on Γ N . For commodity, we assume that both regions are imposed by the context; they are non optimizable subsets of the boundary ∂D of the computational domain. The only optimized part Γ of ∂Ω is traction-free; moreover, body forces are omitted. In this context, the displacement of Ω is the unique solution u Ω ∈ H 1 (Ω) d to the linearized elasticity system: The constitutive tensor A Ω is constructed according to the so-called ersatz material approximation: where A is an isotropic Hooke's tensor with unit Young's modulus, Poisson's ratio ν = 0.3, and η 1 is a very small parameter (in practice we take η = 10 −4 ), so that the phase D \ Ω is made of a very soft material.
We consider the minimization problem where Per(Ω) stands for either the relative or the total perimeter of Ω, V T is a volume target, and C(Ω) is the compliance of Ω, that is: A Ω e (u Ω ) : e (u Ω ) dx = Γ N g · u Ω ds.
The density-based counterpart of this problem reads: where we have introduced a suitable relaxation C(u) of the compliance functional, devised in our previous work [19] and further analyzed in [43]. Again, P ε (u) stands for either P Neu ε (u) or P Rob ε (u). The numerical resolution of this problem relies on the same continuation strategy as in the previous Section 6.3 for the ratio ε/h, and a straightforward adaptation of Algorithm 1.

Optimization of a two-dimensional cantilever
We first consider the benchmark 2d cantilever test case. The computational domain D is a rectangle with size 2 × 1 and the clamping region of shapes Γ D is the left-hand side of ∂D. A unit vertical load g = (0, −1) is applied on Γ N = {x = (x 1 , x 2 ) ∈ D, x 1 = 2, 0.4 ≤ x 2 ≤ 0.6}. The domain is discretized with a regular Figure 13. Several iterations of the minimization of a weighted sum of the compliance and the perimeter of a 2d cantilever under a volume constraint in Section 6.4.2. (Upper row) the relative perimeter functional is used; (lower row) the total perimeter functional is used. triangular mesh composed of 7200 P 1 elements. The volume target is set to V T = 0.4|D| and the initial design Ω 0 is the solution of a standard compliance minimization problem under the volume constraint Vol(Ω) = V T (corresponding to the version of (6.11) without perimeter penalization). When solving (6.11), we normalize the compliance by replacing C(u) with C(u) = C(u)/C(χ Ω0 ).
We solve both instances of (6.11) featuring the approximate relative and total perimeter functionals, with a value of the weight coefficient α = 0.1. The results are displayed in Figure 13. In both cases, the small features of the initial design are removed in the course of the process and admittedly, both optimized designs look similar, except for the fact that the shape Ω obtained by using the total perimeter functional presents more natural, "straight" members in the junction region with the boundary of the computational domain D. This expected effect shows that it is actually more natural to rely on the total version of the perimeter functional in such a context, so as to ensure that the boundary ∂D is "seen" as a part of the boundary of shapes, and that the region ∂Ω ∩ ∂D is optimized as such.

Optimization of a 3d cantilever
We eventually deal with the three-dimensional counterpart of the example of the previous section. The computational domain D is now a box with size 4 × 1 × 1; the region Γ D where shapes are clamped corresponds to the left side of ∂D, and a unit vertical load g = (0, 0, −1) is applied on Γ N := {x = (x 1 , x 2 , x 3 ) ∈ D, x 1 = 4, 0.4 ≤ x 2 , x 3 ≤ 0.6}.
We solve the minimization problem (6.11) with a volume target V T = 0.15|D| and α = 0.1. The initial design Ω 0 is that obtained for minimum compliance under the volume constraint Vol(Ω) = V T , i.e. we solve the instance of (6.11) where the perimeter functional is omitted. The resulting optimized shapes are represented in Figure 14.
Obviously, the minimization of the compliance of Ω without any perimeter constraint tends to favor plates and laminates as structural elements. In contrast, beams and bars naturally appear when a perimeter penalization is incorporated into the problem. As expected, when compared with that obtained when using the relative perimeter functional, the optimized design using the total perimeter functional shows a lesser contact area with the boundary ∂D of the computational domain.

Conclusions and perspectives
The object at the place of honor in this article is the total perimeter Per T (Ω). This functional evaluates the total measure of the boundary of a subset Ω contained in a fixed "hold-all" domain D ⊂ R d , as opposed Figure 14. Optimization of a weighted sum of the compliance and the perimeter of a 3d cantilever under volume constraint in Section 6.4.3; (upper row) several views of the initial guess and (middle row) corresponding views of the optimized shape obtained when using the relative perimeter functional; (lower row) corresponding views of the optimized shape obtained when using the total perimeter functional.
to the relative perimeter Per R (Ω), which only takes into account the region ∂Ω ∩ D lying strictly inside D. We have introduced and analyzed two approximate versions P conv ε (u) and P Rob ε (u) of Per T (Ω), which are in addition defined on general "density functions" u : D → [0, 1]. These prove particularly relevant in densitybased topology optimization frameworks, since their expressions do not involve the gradient of the function u, but rather a "smoothing" of u (the degree of smoothing being controlled by the "small" parameter ε). The "consistency" of the approximate formulas P conv ε (u) and P Rob ε (u) with the original notion of total perimeter as ε tends to 0 has been proved by considering their pointwise convergence, the convergence of their derivatives, and their Γ-convergence. These features of our approximate total perimeter functionals have been appraised by numerical simulations.
The present work opens the way to several natural directions for future research: -It would be natural to try and extend the present study to the case of anisotropic perimeter functionals, of the form ∂Ω ϕ(n Ω ) ds, where ϕ : R d → R is a given function. Grossly speaking, these functionals are used to account for growth models where some expansion directions are preferred, see e.g. [67]. Let us mention that variants of the Modica-Mortola approximation have been proposed to deal with anisotropic perimeter functionals [11,65]. -The result of Proposition 4.8 suggests approximate versions of shape functionals involving the mean curvature of ∂Ω, which would be of great numerical interest, since calculating curvatures, and even worse, their derivatives, is well-known to be a difficult issue (again, see [41] for a related discussion).
-As we have mentioned, the total perimeter functional, or an approximate counterpart as those devised in this article, finds a natural use in optimal design involving contact mechanics: the contact region is often imposed by the context as a fixed subset of the boundary of the computational domain D. However, this region should be counted as part of the boundary of the optimized structure Ω and its area has to be taken into account in the evaluation of the perimeter of Ω.

A.1 A wee bit of geometry
For the convenience of the reader, we recall in the present section a few basic facts from differential geometry, referring to e.g. [38] for a more exhaustive treatment. Throughout the article, we consider Lipschitz bounded domains (resp. domains of class C k , k ≥ 1) Ω ⊂ R d : for any point x ∈ ∂Ω, up to rotating the canonical orthonormal frame, there exist an open neighborhood U of x and a Lipschitz function (resp. a function of class C k ) f : R d−1 → R such that: roughly speaking, Ω can be described as the subgraph of a (d − 1)-dimensional Lipschitz function, locally around any point x ∈ ∂Ω. Likewise, letting the function φ : U → R be defined by φ(x) = x d − f (x 1 , . . . , x d−1 ), it holds: x ∈ Ω ∩ U ⇔ φ(x) < 0.
When Ω is a bounded Lipschitz domain, the unit normal vector n to ∂Ω pointing outward Ω exists for dH d−1 a.e. point x ∈ ∂Ω. If Ω is of class C k , n is defined as a vector field of class C k−1 on ∂Ω. In either case, it is conveniently expressed in terms of the function f (or alternatively φ) if a local representation of the form (A.1) is assumed: Let us now suppose Ω to be of class C 2 ; if x is a point on ∂Ω and a local representation of ∂Ω of the form (A.1) is assumed, the second fundamental form of ∂Ω at x is given by the matrix II x defined by: In the particular case where n(x) = e d , these formulas rewrite: ∇f (x) = 0, and II x = −∇ 2 f (x).
The eigenvalues κ 1 (x), . . . , κ d−1 (x) of the matrix II x are the principal curvatures of the surface ∂Ω at x, and the associated eigenvectors are tangent vectors to ∂Ω at x called the principal directions of ∂Ω at x. Eventually, the mean curvature κ(x) of ∂Ω at x is the sum κ(x) := κ 1 (x) + . . . + κ d−1 (x).

A.2 A few words about the signed distance function
Let us start with a definition.
Definition A.1. Let Ω ⊂ R d be a bounded Lipschitz domain.
-The signed distance function d Ω to Ω is defined by: is the usual Euclidean distance from x to ∂Ω. -For a given x ∈ R d , any point y ∈ ∂Ω realizing the minimum in the definition (A.2) is called a projection point of x onto ∂Ω. When there exists a unique such point, it is denoted by p ∂Ω (x) and called the projection of x onto ∂Ω.
The following proposition expresses the fact that, provided Ω is regular enough, there exists a tubular neighborhood of ∂Ω whose points can be parametrized by their distance to ∂Ω; see [13] for a proof.
When Ω is at least of class C 2 , the above considerations allow to define natural extensions of the unit normal vector n to ∂Ω pointing outward Ω, and of any tangential vector field τ : ∂Ω → R d to the whole tubular neighborhood ∂Ω δ via the formulas: ∀x ∈ ∂Ω δ , n(x) := n(p ∂Ω (x)), and τ (x) = τ (p ∂Ω (x)).
Eventually, let Ω ⊂ R 2 be a 2d bounded domain of class C 2 . Let us denote by τ the unit tangential vector field to ∂Ω, so that (τ (x), n(x)) defined a local orthonormal frame of the plane at any point x ∈ ∂Ω. For any function f of class C 2 on D, the following formula holds:

A.3 A version of the coarea formula
The following result was proved in [5] as an application of the coarea formula (see [36]). (1 + d Ω (z)κ i (y)) dH 1 (z) ds(y), (A.4) where d Ω is the signed distance function to Ω, z denotes a point in the ray ray ∂Ω (y) := {x ∈ D, p ∂Ω (x) = y} emerging from y ∈ ∂Ω and dH 1 (z) is the line integration along that ray.