{'fraction_non_alphanumeric': 0.0762273434390941, 'fraction_numerical': 0.03657610071994501, 'mean_word_length': 3.404397705544933, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Topological optimization via cost penalization

We consider general shape optimization problems governed by Dirichlet boundary value problems. The proposed approach may be extended to other boundary conditions as well. It is based on a recent representation result for implicitly defined manifolds, due to the authors, and it is formulated as an optimal control problem. The discretized approximating problem is introduced and we give an explicit construction of the associated discrete gradient. Some numerical examples are also indicated.


Introduction
Shape optimization is a relatively young branch of mathematics, with important modern applications in engineering and design. Certain optimization problems in mechanics, thickness optimization for plate or rods, geometric optimization of shells, curved rods, drag minimization in fluid mechanics, etc are some examples. Many appear naturally in the form of control by coefficients problems, due to the formulation of the mechanical models, with the geometric characteristics entering the coefficients of the differential operators. See [15], Ch. 6, where such questions are discussed in details.
It is the aim of this article to develop an optimal control approach, using penalization methods, to general shape optimization problems as investigated in [20], [23], [5], [10], [8], etc. We underline that our methodology allows simultaneous topological and boundary variations.
Here, we fix our attention on the case of Dirichlet boundary conditions and we study the typical problem (denoted by (P)): min Ω E j (x, y Ω (x)) dx, (1.1) −∆y Ω = f in Ω, (1.2) y Ω = 0 on ∂Ω, (1.3) where E ⊂⊂ D ⊂ R 2 are given bounded domains, D is of class C 1,1 and the minimization parameter, the unknown domain Ω, satisfies E ⊂ Ω ⊂ D and other possible conditions defining a class of admissible domains. Notice that the case of dimension two is of interest in shape optimization. Moreover, f ∈ L 2 (D), j : D × R → R is some Carathéodory mapping. More assumptions or constraints will be imposed later. Other boundary conditions or differential operators may be handled as well via this control approach and we shall examine such questions in a subsequent paper. For fundamental properties and methods in optimal control theory, we quote [12], [7], [2], [16]. The problem (1.1)-(1.3) and its approximation are strongly non convex and challenging both from the numerical and theoretical points of view. The investigation from this paper continues the one in [26] and is essentially based on the recent implicit parametrization method as developed in [25], [18], [24], that provides an efficient analytic representation of the unknown domains.
The Hamiltonian approach to implicitly defined manifolds will be briefly recalled together with other preliminaries in Section 2. The precise formulation of the problem and its approximation is analyzed in Section 3 together with its differentiability properties. In Section 4, we study the discretized version and find the general form of the discrete gradient. The last section is devoted to some numerical experiments, using this paper approach.
The method studied in this paper has a certain complexity due to the use of Hamiltonian systems and its main advantage is the possibility to extend it to other boundary conditions or boundary observation problems. This will be performed in a subsequent article.
Remark 2.1 In higher dimension, iterated Hamiltonian systems were introduced in [25] and uniqueness and regularity properties are proved. Some relevant examples in dimension three are discussed in [18]. In the critical case, generalized solutions can be obtained [24], [25].
We define now the family F of admissible functions g ∈ C 2 (D) that satisfy the conditions: (2.8) Condition (2.6) says that G ∩ ∂D = ∅ and condition (2.7) is an extension of (2.5). In fact, it is related to the hypothesis on the non existence of equilibrium points in the Poincare-Bendixson theorem, [11], Ch. 10, and the same is valid for the next proposition. The family F defined by (2.6)-(2.8) is obviously very rich, but it is not "closed" (we have strict inequalities). Our approach here, gives a descent algorithm for the shape optimization problem (P) and existence of optimal shapes is not discussed.
Following [26], we have the following two propositions: (2.9) We also introduce the neighborhood V , > 0 where d[(x 1 , x 2 ), G] is the distance from a point to G.

Proposition 2.3
If > 0 is small enough, there is λ( ) > 0 such that, for |λ| < λ( ) we have G λ ⊂ V and G λ is a finite union of class C 2 closed curves.

Remark 2.2
The inclusion G λ ⊂ V shows that G λ → G for λ → 0, in the Hausdorff-Pompeiu metric [15]. In a "small" neighborhood of each component of G there is exactly one component of G λ if |λ| < λ( ), due to this convergence property and the implicit functions theorem applied in the initial condition of the perturbed Hamiltonian system derived from (2.1)-(2.3) .
Assume that (x 1λ , x 2λ ) has the period T g+λr > jT g . Since (x 1 , x 2 ) is periodic with period T g and relations (2.12)-(2.13) show that (x 1λ , x 2λ ) is very close to (x 1 , x 2 ) in every t ∈ [0, jT g ] it yields that (x 1λ , x 2λ ) is, as well, surrounding E at least once. As it may have no self intersections, it yields that (x 1λ , x 2λ ) is as a limit cycle around E. Such arguments appear in the proof of the Poincaré-Bendixson theorem, [11], Ch. 10. That is (x 1λ , x 2λ ) cannot be periodic -and this is a false conclusion due to Proposition 2.3.
Consequently, we get T g+λr ≤ jT g , |λ| < λ( ). (2.14) On a subsequence, by (2.14), we obtain T g+λr → T * ≤ jT g . We assume that T * = T g . It is clear that (x 1 (T * ), x 2 (T * )) = (x 1 (T g ), x 2 (T g )) by the definition of the period. However, relation (2.12) and the related convergence properties give the opposite conclusion. This contradiction shows that T * = T g and the convergence is valid on the whole sequence. 2 Remark 2.3 Usually, the perturbation of a periodic solution may not be periodic and just asymptotic convergence properties are valid, under certain assumptions, Sideris [22]. A natural question, taking into account (2.12)-(2.13), is whether |T g+λr − T g | ≤ c|λ|, for c > 0 independent of λ, |λ| < λ( ).

The optimization problem and its approximation
that contains E. Clearly E ⊂ Ω g by (2.8). Notice as well that the domain Ω g defined by (3.1) (we use this notation for the domain as well) is not simply connected, in general. This is the reason why the approach to (1.1)-(1.3) that we discuss here is related to topological optimization in optimal design problems. But, it also combines topological and boundary variations. The penalized problem, > 0, is given by: 3) y = 0, on ∂D, (3.4) where z g = (z 1 g , z 2 g ) satisfies the Hamiltonian system (2.1)-(2.3) in I g with some (x 0 and I g = [0, T g ] is the period interval for (3.5)-(3.7), due to Proposition 2.2. The problem (3.2)-(3.7) is an optimal control problem with controls g ∈ F and u ∈ L 2 (D) distributed in D. The state is given by [y , z 1 g , z 2 g ] ∈ H 2 (D) × (C 2 (I g )) 2 . We also have y ∈ H 1 0 (D). In case the corresponding domain Ω g is not simply connected, in (3.7) one has to choose initial conditions on each component of ∂Ω g and the penalization term becomes a finite sum due to Proposition 2.2. The method enters the class of fixed domain methods in shape optimization and can be compared with [13], [14], [17]. It is essentially different from the level set method of Osher and Sethian [19], Allaire [1] or the SIMP approach of Bendsoe and Sigmund [3]. From the computational point of view, it is easy to find initial condition (3.7) on each component of G and the corresponding period intervals I g associated to (3.5)-(3.7). See the last section as well.
We have the following general suboptimality property: Proof. Let [y gm , g m ] ∈ H 2 (Ω gm )×F be a minimizing sequence in the problem (1.1)-(1.3), (3.1). By the trace theorem, since ∂Ω gm and D are at least C 1,1 under our assumptions, there is y gm ∈ H 2 (D) ∩ H 1 0 (D), not unique, such that y gm = y gm in Ω gm . We define the control u gm ∈ L 2 (D) as following: where Ω gm is the open set defined in (3.1). Notice that on the second line in the above formula, we have no singularity. It is clear that the triple [ y gm , g m , u gm ] is admissible for the problem (3.2)-(3.7) with the same cost as in the original problem (1.1)-(1.3) since the penalization term in (3.2) is null due to the boundary condition (1.3) satisfied by y gm . Consequently, there is n(m) sufficient big, such that E j x, y n(m) (x) dx + 1 Since j is bounded from below, we get from (3.8): with C a constant independent of > 0. Then, (3.9) shows that (1.3) is fulfilled with a perturbation of order 1/2 . Moreover, again by (3.8), we see the minimizing property of {y n(m) } in the original problem (P).
We notice that in the state equation (3.3), the right-hand side coincides with f in the set {(x 1 , x 2 ) ∈ D; g(x 1 , x 2 ) < − } , which is an approximation of Ω g n(m) . Namely, we notice that for any g ∈ F, the open sets ∈ Ω n = Ω g+ n , for n > 0 sufficiently small. Consequently, we have the desired convergence property by [15], p. 461. This ends the proof. 2 Remark 3.1 A detailed study of the approximation properties in the penalized problem is performed in [26], in a slightly different case.

We consider now variations
Notice that u + λv ∈ L 2 (D) and g + λr ∈ F for |λ| sufficiently small . The conditions (2.6), (2.7), (2.8) from the definition of F are satisfied for |λ| sufficiently small (depending on g) due to the Weierstrass theorem and the fact that E, ∂D and G are compacts. Here, we also use Proposition 2.3. Consequently, we assume |λ| "small". We study first the differentiability properties of the state system (3.3)-(3.7): where q = lim λ→0 being the solution of (3.3)-(3.4) corresponding to g + λr, u + λv, and z g+λr ∈ C 1 (I g ) 2 is the solution of (3.5)-(3.7) corresponding to g + λr. The limits exist in the above spaces. We denote by "·" the scalar product on R 2 .
Proof. We subtract the equations corresponding to y λ and y and divide by λ = 0, small: with 0 boundary conditions on ∂D. The regularity conditions on F and u, v ∈ L 2 (D) give the convergence of the right-hand side in (3.15) to the right-hand side in (3.10) (strongly in L 2 (D)) via some calculations. Then, by elliptic regularity, we have y λ −y λ → q strongly in H 2 (D) ∩ H 1 0 (D) and (3.10), (3.11) follows. For (3.12)-(3.14), the argument is the same as in Proposition 6, [24]. The convergence of the ratio z g+λr −zg λ is in C 1 (I g ) 2 on the whole sequence λ → 0, due to the uniqueness property of the linear system (3.12)-(3.14). Here, we also use Remark 2.2, on the convergence G λ → G and the continuity with respect to the perturbations of g in the Hamiltonian system (2.1)-(2.3), according to [24]. 2 Remark 3. 2 We have as well imposed the condition where (x 0 1 , x 0 2 ) ∈ D \ E is some given point. Similarly, constraints like (3.16) may be imposed on a finite number of points or on some curves in D \ E and their geometric meaning is that the boundary ∂Ω g of the admissible unknown domains should contain these points, curves, etc.
Proof. In the notations of Proposition 3.2, we compute In (3.18), λ > 0 is "small" and Proposition 2.3 ensures that g + λr ∈ F (see [25] as well). By Proposition 2.2 we know that the trajectories associated to g + λh are periodic, that is the functions in the second integrals are defined on I g . Moreover, since f, u ∈ L p (D), then y λ , y defined as in (3.3) due to the hypothesis on y (z g (t)). We can study term by term: By Proposition 3.2, each of the above three integrands are uniformly bounded and their limits can be easily computed, for instance on [0, 2T g ] due to Proposition 2.4. Notice, in the last term, that |z g (t)| = |∇g (z g (t)) | = 0 due to (3.5), (3.6) and (3.7), that is we can differentiate here as well.
Again by Proposition 2.4 and the above uniform boundedness, we infer that each of the above three terms has null limit as λ → 0, i.e. L = 0. Consequently, it is enough to study the limit (3.18).
We also have y λ → y in C 1 (D), for λ → 0, by elliptic regularity. Then, under the assumptions on j(·, ·), we get For the second integral in (3.18), we intercalate certain intermediary terms and we compute their limits for λ → 0: due to the convergence z g+λr → z g in C 1 (I g ) 2 by g, r ∈ C 2 (D) and the continuity properties in (2.1)-(2.3); where w = (w 1 , w 2 ) satisfies (3.12)-(3.14) and again we use the regularity and the convergence properties in C 1 (D), respectively where we recall that |z g (t)| = (z 1 g ) (t) 2 + (z 2 g ) (t) 2 is non zero by (2.7) and the Hamiltonian system, and standard derivation rules may be applied under our regularity conditions.

Remark 3.3
In the case that Ω g is not simply connected, the penalization integral in (3.2) is in fact a finite sum and each of these terms can be handled separately, in the same way as above, due to Proposition 2.3 and Remark 2.2. The significance of the hypothesis y (z g (t)) = 0, is that one should first minimize the penalization term with respect to the control u (this is possible due to the arguments in the proof of Proposition 3.1). Then, the obtained control should be fixed and the minimization with respect to g ∈ F is to be performed. In case T g+λr can be evaluated as in Remark 2.3, then the hypothesis can be relaxed to y (x 0 1 , x 0 2 ) = 0 via a variant of the above arguments.

Now, we denote by
the linear continuous operator given by r, v → q , defined in (3.10), (3.11). We also denote by B : C 2 (D) → C 1 (I g ) 2 the linear continuous operator given by (3.12)-(3.14), Br = [w 1 , w 2 ]. In these definitions, g ∈ C 2 (D) and u ∈ L p (D) are fixed. We have: Corollary 3.1 The relation (3.17) can be rewritten as: where the vector C(t) is explained below.
Proof. In the last integral in (3.17), we replace w (t) by the right-hand side in (3.12), (3.13). We compute: We denote by C(t) the (known) vector and together with (3.24), we get (3.23). This ends the proof. 2

Finite element discretization
We assume that D and E are polygonal. Let T h be a triangulation of D with vertices A i , i ∈ I = {1, . . . , n}. We consider that T h is compatible with E, i.e.
where T designs a triangle of T h and h is the size of T h . We consider a triangle as a closed set. For simplicity, we employ piecewise linear finite element and we denote We use a standard basis of W h , {φ i } i∈I , where φ i is the hat function associated to the vertex A i , see for example [6], [21]. The finite element approximations of g and u We set the vectors G = (G i ) T i∈I ∈ R n , U = (U i ) T i∈I ∈ R n and g h can be identified by G, etc. The function u is in L p (D), as in Proposition 3.3. Alternatively, for u h , we can use discontinuous piecewise constant finite element P 0 . In order to approach g ∈ C 2 (D), we can use high order finite elements.

Discretization of the optimization problem
We introduce V h = {ϕ h ∈ W h ; ϕ h = 0 on ∂D}, I 0 = {i ∈ I; A i / ∈ ∂D} and n 0 = card(I 0 ). The finite element weak formulation of (3.3)-(3.4) is: find In order to obtain the linear system, we take the basis functions ϕ h = φ i in (4.1) for i ∈ I 0 . Let us consider the vector The matrix K is symmetric, positive definite. The finite element approximation of the state system (3.3)-(3.4) is the linear system: Now, we shall discretize the objective function (3.2). We denote I E = {i ∈ I; A i ∈ E} and n E = card(I E ). For the first term of (3.2), we introduce We shall study the second term of (3.2). In order to solve numerically the ODE system (3.5)-(3.7), we use a partition [t 0 , . . . , t k , . . . , t m ] of [0, T g ], with t 0 = 0 and t m = T g . We can use the forward Euler scheme: for k = 0, . . . , m − 2. We set Z k = (Z 1 k , Z 2 k ) and we impose Z m = Z 0 . In fact, Z k is an approximation of z g (t k ). We do not need to stock Z 0 and we set Z = (Z 1 , Z 2 ) ∈ R m ×R m , In the applications, one can use more performant numerical methods for the ODE's, like explicit Runge-Kutta or backward Euler, but here we want to avoid a too tedious exposition.
Without risk of confusion, we introduce the function Z : [0, T g ] → R 2 defined by for k = 0, 1, . . . , m − 1. We have Z(t k ) = Z k and we can identify the function Z(·) by the vector Z ∈ R m × R m . We remark that Z(·) is derivable on each interval (t k , t k+1 ) and subject to (4.2). We point out that Y depends on G and U by (4.2) and Z depends on G by (4.3)-(4.5).

Discretization of the derivative of the objective function
Let r h , v h be in W h and R, V ∈ R n be the associated vectors. The finite element weak formulation of (3.10)-(3.11) is: find Let Q ∈ R n 0 be the associated vector to q h and we construct the n 0 ×n matrix C 1 (G, , U ) defined by The linear system of (4.7) is In order to approximate ∂ 2 j(x, y (x)), y given by (3.3)-(3.4), we consider the nonlinear application We define the n E × n 0 matrix M ED defined by The first term of (3.17) is approached by and the second term of (3.17) is approached by where the matrix N (Z) was introduced in the previous subsection. Next, we introduce the partial derivative for a piecewise linear function. Let g h ∈ W h and G ∈ R n its associated vector, i.e. g h ( where J i is the set of index j such that the triangle T j has the vertex A i . Since g h is a linear function in each triangle T j , then ∂ 1 g h |T j is constant. Similarly, we construct Π 2 h G ∈ R n for ∂ 2 . In fact, Π 1 h and Π 2 h are two n × n matrices depending on T h . Then, we set and similarly for ∂ h 2 g h . Finally, we put  1), A 3 = (0, 1) and the triangulation of two triangles . We shall present the discrete derivative of the hat function We have J 1 = {1, 2} and In order to solve the ODE system (3.12)-(3.14), we use the forward Euler scheme on the same partition as for (4.3)-(4.5): for k = 0, . . . , m − 1. We set W k = (W 1 k , W 2 k ) and now we have W m = W 0 generally. In fact, W k is an approximation of w(t k ). We do not need to stock W 0 and we set W = (W 1 , W 2 ) ∈ R m × R m , with W 1 = (W 1 k ) T 1≤k≤m and W 2 = (W 2 k ) T 1≤k≤m . As mentioned before, we can use more performant numerical methods for the ODE, like explicit Runge-Kutta or backward Euler.
We construct W : [0, T g ] → R 2 in the same way as for Z(t) If ψ k is the one-dimensional piecewise linear hat function associated to the point t k of the partition [t 0 , . . . , t k , . . . , t m ], we can write equivalently The third term of (3.17) is approached by where |Z k Z k+1 | is the length of the segment in R 2 with ends Z k and Z k+1 . We have We introduce the n 0 × n matrices N then (4.15) can be rewritten as and finally the third term of (3.17) is approached by We can introduce the linear operators T 1 (Z) and T 2 (Z) by (Z) (4.17) then (4.16) can be rewritten as The fourth term of term of (3.17) is approached by But Z (t) and W (t) are constants for t k ≤ t < t k+1 , then where |Z k Z k+1 | is the length of the segment in R 2 with ends Z k and Z k+1 . We introduce the n 0 × n 0 matrix R k (Z) defined by The (4.19) can be rewritten as The study of this subsection can be resumed as following: The discret version of (3.17) is Proof. We get (4.22) just by assembling (4.9), (4.10), (4.18) and (4.21). 2

Discretization of the formula (3.23)
From (4.8), we get and the discrete version of the operator A in the Corollary 3.1 is Replacing Q in the first two terms of (4.22), we get We denote The third term of (3.23) is approached by and using the trapezoidal quadrature formula on each sub-interval [t k , t k+1 ], we get Similarly, for the 4th and 5th terms of (3.23), we get In order to write (4.24)-(4.26) shorter, we introduce the vectors: . Also, we introduce the vectors in R n : Proposition 4.2 The discrete version of the (3.23) is Proof. We obtain (4.27) by summing (4.23)-(4.26). 2 Next, we give more details about the relationship between W and R. Let us introduce the 2 × 2 matrices , I 2 = 1 0 0 1 and the 2 × n matrice We remark that M 2 depends on G and Z and N 2 on Z. The system (4.11)-(4.12) can be written as

Proposition 4.3
We have the following equality where at the right-hand side, M 2m is a 2m × 2m matrix defined by and the size of the second matrix, which contains N 2 , is 2m × n.
Proof. From (4.13) and the recurrent relation, we have Since W depends on R by (4.28), we can introduce the linear operator approximation of B in the Corollary 3.1 If we denote by i the i-th line of the matrix M 2m at the right-hand side of (4.28), for 1 ≤ i ≤ 2m, then . . .
and B 2 (G, Z), B 3 (G, Z) are m × n matrices. The size of the matrix containing N 2 is 2m × n.

Gradient type algorithm
We start by presenting the algorithm.
Step 1 Start with k = 0, > 0 some given "small" parameter and select some initial (G k , U k ).
Step 2 Compute Y k the solution of (4.2) and Z k solution of (4.3)-(4.5). Step Step Step 5 If |J(G k+1 , U k+1 )−J(G k , U k )| is below some prescribed tolerance parameter, then Stop. If not, update k := k + 1 and go to Step 3.

In the
Step 3, we have to provide a descent direction. We present in the following a partial result. Let us introduce a simplified adjoint system: find for all ϕ h ∈ V h and with Z(t) given by (4.3)-(4.5). We have p h (x) = i∈I 0 P i φ i (x) and P = (P i ) T i∈I 0 ∈ R n 0 . The linear system associated to (4.29) is We recall that K and N (Z) are symmetric matrices.
Proposition 4.4 Given g h , u h ∈ W h , let y h ∈ V h the solution of (4.1). If r h = −p h u h and v h = −p h , where p h ∈ V h is the solution of (4.29), then where q h ∈ V h is the solution of (4.7) depending on r h and v h .
Proof. Putting ϕ h = p h in (4.7) and ϕ h = q h in (4.29), we get For v h = −p h , we have This ends the proof. 2

Remark 4.1
The left-hand side of (4.30) represents the first two terms of (4.22). We can obtain a similar result as in Proposition 4.4, without using the adjoint system, by taking in place of V and R in (4.23). In this case, (4.23) is different from the direction given by Proposition 4.4. Now, we present a descent direction, obtained from the complete gradient of the discrete cost (4.22). Proposition 4.5 For (R * * , V * ) ∈ R n × R n given by we obtain a descent direction for J at (G, U ).
Proof. In (4.27), we replace W 1 by B 2 (G, Z)R and W 2 by B 3 (G, Z)R, we obtain that

Numerical tests
Shape optimization problems and their penalization are strongly nonconvex. The computed optimal domain depends on the starting domain, but also on the penalization or other numerical parameters. It may be just a local optimal solution.
Moreover, the final computed value of the penalization integral is small, but not null. This allows differences between the optimal computed domain Ω g and the zero level curves of the computed optimal state y . Consequently, we compare the obtained optimal cost in the penalized problem with the costs in the original problem (1.1) -(1.3) corresponding to the optimal computed domain Ω g and the zero level curves of y . This is a standard procedure, to inject the approximating optimal solution in the original problem. Notice that in all the experiments, the cost corresponding to Ω g is the best one, but the differences with respect to the other computed cost values are small. This shows that the rather complex approximation/penalization that we use is reasonable. Its advantage is that it may be used as well in the case of boundary observation or for Neumann boundary conditions and this will be performed in a subsequent paper.
In the examples, we have employed the software FreeFem++, [9]. The mesh of D has 73786 triangles and 37254 vertices. The penalization parameter is = 10 −3 and the tolerance parameter for the stopping test at the Step 5 of the algorithm is tol = 10 −6 . The initial domain is the disk of center (0, 0) and radius 2.5 with a circular hole of center (−1, −1) and radius 0.5.
At the Step 3 of the Algorithm, we use (R k , V k ) given by Proposition 4.4. At the Step 4, in order to have E ⊂ Ω k , we use a projection P at the line search λ k ∈ arg min λ>0 J P(G k + λR k ), U k + λV k and G k+1 = P(G k + λ k R k ). If the value of g k h + λr k h at a vertex from E is positive, then we set this value to −0.1. We recall that the left-hand side of (4.30) represents only the first two terms of (4.22), not the whole gradient. If r h , v h are given by Proposition 4.4 and γ > 0 is a scaling parameter, then γr h and v h verify (4.30), that is they also give a descent direction. We take the scaling parameter for r h given by γ = 1 max(r h ) , that is a normalization of r h . In this way we avoid the appearance of very high values of the objective function, that may stop the algorithm even in the first iteration. For the line search at the Step 4, we use λ = ρ i λ 0 , with λ 0 = 1, ρ = 0.5 for i = 0, 1, . . . , 30. At the final iteration, the first term of the optimal objective function is 1.03796 and ∂Ωg y 2 (s)ds = 1.39471 × 10 −2 . We point out that the optimal Ω g has a hole and the penalization term is a sum of two integrals ∂Ωg y 2 (s)ds = 2 j=1 I j (y (z g (t))) 2 |z g (t)|dt where the integral over I 1 corresponding to the exterior boundary of Ω g and I 2 to the boundary of the hole. In Figure 1 in the bottom, we can see the computed optimal state y in iteration 94. We also compute the costs E (y 1 − y d ) 2 dx = 0.998189 where y 1 is the solution of the initial elliptic problem (1.2)-(1.3) in the domain Ω g with g obtained in iteration 94 and E (y 2 − y d ) 2 dx = 1.04032 where y 2 is the solution of the elliptic problem (1.2)-(1.3) in the domain bounded by the zero level sets of y in iteration 94, see Figure 1. 1 − x 2 2 + 1, we have the exact optimal state y = y d defined in the disk of center (0, 0) and radius 1, that gives an optimal domain of the problem (1.1)-(1.3).
We have used = 10 −1 and the starting configuration: the disk of center (0, 0) and radius 2.5 with the circular hole of center (−1, −1) and radius 0.5. We use (R k , V k ) given by Proposition 4.4. The parameters for the line search and γ are the same as in the precedent example.
The stopping test is obtained for k = 64. The initial and the final computed values of the objective function are 5368.84 and 11.2311. We obtain a local minimum that is different from the above global solution. The first term of the final computed objective function is 0.472856. The term ∂Ωg y 2 (s)ds is 1.07583 and it was computed over the exterior boundary as well as over the boundaries of two holes. The length of the total boundary of the optimal domain is 23.9714 and of the initial domain is 2π(2.5 + 0.5) = 18.8495. The domain changes its topology. The computed optimal state y is presented in Figure 2 in the bottom. At the left, we show y 1 the solution of the elliptic problem (1.2)-(1.3) in the domain Ω g which gives E (y 1 − y d ) 2 dx = 0.295178, at the right we show y 2 the solution of the elliptic problem (1.2)-(1.3) in the domain bounded by the zero level sets of y , which gives E (y 2 − y d ) 2 dx = 0.471788.

Example 3.
We have also used the descent direction given by Proposition 4.5, for the starting configuration the disk of center (0, 0) and radius 1.5, = 10 −1 , γ = 1 r h ∞ and a mesh of  Finally, we notice that the hypothesis of Proposition 3.3 is obviously fulfilled by the null level sets of y with a corresponding parametrization. In an approximate sense, it is also fulfilled by the computed optimal domain Ω g since the penalization integral is "small" in all the examples.