# Two Algorithms for a Class of Elliptic Problems in Shape Optimization

Authors: Li Tian; Dai Xiaoxia; Zhang Chengwei
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

We propose two algorithms for elliptic boundary value problems in shape optimization. With the finite element method, the optimization problem is replaced by a discrete variational problem. We give rules and use them  to decide which elements are to be reserved. Those rules are determined by the optimization; as a result, we get the optimal design in shape. Numerical examples are provided to show the effectiveness of our algorithms.


---

## Body

## 1. Introduction

With numerous applications in structural engineering and fluid mechanics, shape optimization has received much more attention in recent years [1–5]. Typical problems are to find optimal shapes that minimize or maximize some design objectives with certain PDE and geometry constraints. The essential difficulty for such problems is that the optimal shape or topology is unknowna priori. Many different numerical methods are used to solve this type of problems. The classical method of shape sensitivity analysis has been studied in [6, 7]. One drawback of this method is its computational effort due to remeshing.The most recent topology approach is the homogenization method [8]. A relaxed formulation of the original shape optimization problem should be introduced. Then, penalization techniques are needed to remove composite designs in the final shape. The method is very promising in that it is nearly independent on the initial guess for design. Moreover, it has a rigorous mathematical result of the existence of the solution. However, this method can be applied to limited areas and particular objective functions.Level set methods are popular for solving shape optimization problems (see [5, 9–13]). Level set method first proposed by Osher and Sethian [14] is a versatile tool when dealing with problems requiring tracing interfaces separating a domain into subregions. For details of the method and its wide applications, we refer to the monograph [15] and the survey [16]. Level set methods can easily handle certain types of shape and topological changes, such as merging, splitting, and developing sharp corners, which can be done by the evolution of the level set function. Therefore, level set methods can naturally be applied to solve shape and topology optimization problems. Its cost is moderate since the shape is captured on a fixed Eulerian mesh. However, by its nature, the obtained level set function decreases the number of connected components of resulting geometrical domains during the optimization process. Besides, the convergence to the optimal geometrical domain is slow due to the necessary reinitialization process.In this paper, we propose a novel method for a class of linear elliptic problems in shape optimization. The main difference from other methods is that we perform optimization in each element, which means that we deal with this problem in a discrete formula. First, let us introduce the model problem. LetU be a bounded domain of R
2 with a boundary ∂
U
=
Γ
D. Let V also be a subset of R
2 and satisfy V
⊂
⊂
U. For any open set ω
⊂
R
2, we denote by #
ω
¯ the number of connected components of ω
¯ and we consider the set of admissible domains: (1)
O
k
=
Ω
=
U
∖
ω
¯
;
ω
⊂
V
,
#
ω
¯
<
k
.For any Ω
∈
O
k, k
⩾
1, the boundary of Ω can be denoted by ∂
Ω
=
Γ
D
∪
Γ
N with Γ
N
=
∂
ω. The boundary Γ
N receives Neumann boundary conditions, while we have Dirichlet conditions on Γ
D. We illustrate the details in Figure 1.Figure 1
An admissible domainΩ.Now, let us consider the following elliptic boundary value problem:(2)
-
Δ
u
+
u
=
f
in
Ω
,
u
=
0
on
⁡
Γ
D
,
∂
u
∂
n
=
0
in
Γ
N
,where f
∈
C
∞
(
U
¯
) and n is the unit outwards normal vector to ∂
Ω. Our goal is to optimize the shape of domain Ω which maximizes a given cost functional depending on the domain and the solution u. Here, we consider the energy functional: (3)
E
Ω
,
u
≔
1
2
∫
Ω
∇
u
2
+
u
2
d
x
-
∫
Ω
f
u
d
x
=
-
1
2
∫
Ω
∇
u
2
+
u
2
d
x
=
-
1
2
∫
Ω
f
u
d
x
,and let us introduce the regularized shape functional: (4)
J
Ω
,
u
≔
E
Ω
,
u
+
λ
A
Ω
,with A
(
Ω
) defined by A
(
Ω
)
≔
|
Ω
|, which is the Lebesgue measure of Ω in R
2. Constant λ are positive and allow us to take into account the area. Then, for any fixed integer k
⩾
1, the shape optimization problem is the following:(5)
max
⁡
J
Ω
:
Ω
∈
O
k
.Since the exterior boundary ofΩ is fixed once for all and equals ∂
U, the domains Ω satisfying (5) are actually determined by their moving boundaries Γ
N.Our algorithm comes from greedy algorithm. A greedy algorithm might also be called a “single-minded” algorithm or an algorithm that gobbles up all of its favorites first. It may not completely solve the problem or if it produces a solution, it may not be the very best one, but it is one way of approaching the problem and sometimes yields very good results. It is well known that greedy algorithms are usually used in discrete mathematics, for example, the salesman problem, graph coloring problem, and other NP-complete problems [17]. And a few papers have reported that the greedy algorithms can also work well for combinatorial optimization problems [18, 19]. In this paper, we propose numerical algorithms based on the greedy algorithm. In one step of our algorithm, we compute the shape functional in each element to find some elements according to the given rules. Then, we throw those elements away (they may come back at some next step). The sum of remain elements is just the domain Ω that is the approximation of our problem.We will propose two types of algorithms in the next section and present some numerical results in Section3.


## 2. Algorithms

In [20], the authors proved that there exists a domain Ω satisfying problem (5). They obtain this result when they apply the results from [21, 22] to this problem with a few changes. Here, the results are listed and the details can be found in [20].Theorem 1.
LetΩ
i
=
U
∖
ω
¯
i be a sequence of open sets such that the number of connected components of ω
¯
i is uniformly bounded. If ω
¯
i converges in the sense of the Hausdorff metric to ω
¯, then the solution u
i of the Neumann problem(6)
-
Δ
u
i
+
u
i
=
f
,
i
n
Ω
i
,
u
i
=
0
o
n
∂
U
,
∂
n
u
i
=
0
,
o
n
∂
ω
iconverges to the solution u of problem (6) on Ω if and only if |
Ω
i
|
→
|
Ω
|.Theorem 2.
For anyk
⩾
0, problem (5) admits at least one solution Ω
∈
O
k.For the sake of convenience, we chooseU
=
(
0,1
)
×
(
0,1
). Let T
h be a partition of the domain U into rectangular elements of mesh size h. Then, we define the bilinear finite element space V
h by (7)
V
h
=
∀
v
∈
H
0
1
U
:
v
k
∈
Q
1
k
,
∀
k
∈
T
h
.We take V
h as finite element space of problem (2); then, we can get the following variational formulation: find u
h
∈
V
h such that(8)
∫
Ω
∇
u
h
·
∇
v
h
+
u
h
v
h
d
x
=
∫
Ω
f
v
h
d
x
,
∀
v
h
∈
V
h
.If we change domain Ω into element k, there is no such type of equation set up since the function does not vanish on boundary of each element. In the element K, we obtain a measurement:(9)
J
K
=
-
1
2
∫
K
f
u
K
d
x
+
λ
K
.Now, we use the greedy strategy: we drop the element K, in which J
(
K
)
=
m
i
n
K
J
(
K
). Accordingly, we can rewrite the shape functional as the sum of each element; namely, (10)
J
Ω
≔
-
1
2
∫
Ω
f
u
d
x
+
λ
A
Ω
=
∑
k
-
1
2
∫
k
f
u
k
d
x
+
λ
K
,where |
K
| is the area of the element K. In our quadrilateral mesh case, |
K
|
=
h
2. In this way, we can propose an iterative algorithm; in each step, we seek and throw away the element whose value of shape functional is the minimum. Since each time we throw away only one element, which is very small compared to the whole domain, this algorithm can increase the shape functional in general. We list the algorithm as in Algorithm 1.Algorithm 1:
(1)i = 0.
(2) Compute the solutionu
i of  (8) with domain Ω
i
=
U∖ω
i, S
J
=
0.
Fork = 1 : Nelement
Find no marked elements, and compute their shape functional
J
i
(
K
)
=
-
(
1
/
2
)
∫
K
f
u
i
d
x
+
λ
|
K
|, and the sum of shape functional
S
J
=
S
J
+
J
i
(
K
);
end
Find the minimum ofJ
i
(
K
) and mark the element.
(3)i = i + 1, Go to Step  2.Here, “i” and “S
J” denote the iteration of the algorithm and the value of shape functional on domain Ω
i, respectively. We denote by “N
e
l
e
m
e
n
t” the number of finite elements of domain U. “ω
i” represents the sum of elements that are thrown away in the preceding i iterations. When i
=
0, we solve (8) with the domain U, because we have no marked elements at first. And this is our initial case. From results of the numerical examples (in Section 3), we can see that Algorithm 1 can obtain an optimal domain. Because we throw away one element in each iteration, the convergence is very slow. And the number of iterations gets growth when the mesh size h is refined. What is more, it is difficult to get the conditions for stopping the iterations. Therefore, we make some changes in Algorithm 1. At each iteration, we throw away several elements to accelerate convergence. But the optimal number of elements that should be thrown away is not knowna priori. Sometimes we may throw away some elements which should be left. Therefore, we should design a process to get such elements back. That is our algorithm (Algorithm 2).Algorithm 2:
(1)i = 0.
(2) Compute the solutionu
i of  (8) with domain Ω
i
=
U∖ω
i, S
J
=
0.
Fork = 1 : Nelement
ComputeJ
K
=
-
(
1
/
2
)
∫
K
f
u
i
d
x;
Find no marked elements and compute their sum of shape functional
S
J
=
S
J
+
J
K
+
λ
|
K
|;
ifJ
K
<
-
λ
|
K
|
θ
Yes: mark;
No: remove the mark.
end
end
(3)i = i + 1, Go to Step  2.The parameterθ controls how much elements will be thrown away in each step. The optimum of θ is different for different problems. According to this algorithm, in each step, we discard not only one element but a batch of elements in general. Furthermore, the elements which are thrown away maybe come back in some next step. In Figure 2, there are three elements and k
1 is the element and we keep it to Ω
i; the other two elements are thrown away. In the next step, after solving the equation u
1, u
2 will be changed and therefore the value of shape functional in element k
2 is changed. Thus, if this value does not satisfy the rule that we set up should be back. This procedure is necessary for our algorithm.Figure 2
Three elements.Both of Algorithms1 and 2 are simple and easy for coding. In the next section, we will see Algorithm 2 enjoys fast convergence for shape optimization problems.


## 3. Numerical Examples

In this section, we present some numerical tests for our algorithms. As mentioned above, we chooseU
=
(
0,1
)
×
(
0,1
) and use bilinear finite element space V
h as our space. We start with the entire domain U in all examples. In order to make numerical comparisons, we use numerical examples in [20]. The data of the first example are (11)
f
=
10
sin
2
⁡
4
π
x
,
for
x
,
y
∈
U
=
0,1
×
0,1
,
λ
=
0.5
.First, we apply Algorithm1 to this example and let h
=
1
/
64. We denote by Ω
1 with green the optimal domain obtained in this example. The value of shape functional J in Ω
1 is J
(
Ω
1
)
=
0.274468. The results are shown in Figure 3. We find that Algorithm 1 requires many iterations and converges very slowly. Besides, there are some oscillations in iterative process. These two flaws are unfavorable for a good algorithm. In the following, let us see the effects of Algorithm 2. Note that Algorithm 2 involves choice of a parameter θ. Considering the little dependence on h, we choose some θ for a relatively large h and adjust it for a small h. Numerical results are shown in Figure 4.Figure 3
Shape optimization with the dataf
=
10
sin
2
⁡
4
π
x, λ
=
0.5 using Algorithm 1. Optimal value: J
(
Ω
)
=
2.7446853.Figure 4
Shape optimization with the dataf
=
10
sin
2
⁡
(
4
π
x
), λ
=
0.5, θ
=
0.7 using Algorithm 2. Optimal value: J
(
Ω
)
=
0.278974.In [20], in order to make comparisons with the level set method with topological derivatives, conventional level set method was used to solve the model problem. In that case, the Hamilton-Jacobi equation is solved on a 101
×
101 grid. At last, the level set method converges to a local maximum J
(
Ω
)
=
0.2640695. The value J
(
Ω
2
)
=
0.278974 obtained by using Algorithm 2 is better than J
(
Ω
)
=
0.2640695. Our method is simple and enjoys fast convergence, although the final functional value 0.282319 by the level set method with topological derivatives is a little better than ours.The second example tests the choice ofλ with f
(
x
,
y
)
=
10
x
+
y. Thus, we get two groups of figures, Figures 5 and 6. When λ is large, from Figure 6, we see the oscillation exists. This phenomenon is caused by the increase of parameters, which also appears in [20]. Our visible result is similar as that in [20].Figure 5
Shape optimization with the dataf
=
10
x
+
y, λ
=
1, θ
=
0.64 using Algorithm 2. Optimal value: J
(
Ω
)
=
0.597596.Figure 6
Shape optimization with the dataf
=
10
x
+
y, λ
=
1.7, θ
=
0.55 using Algorithm 2. Optimal value: J
(
Ω
)
=
1.189091.


## 4. Conclusion

Based on the greedy idea, we have presented two interesting and promising algorithms for solving the linear elliptic problems in shape optimization. Our algorithms are simple but effective for the model problem. Algorithm2 is especially efficient compared with the level set method. We will try to apply our method to other shape and topology problems involving geometry and topology in the future.



---

*Source: 241837-2015-07-16.xml*