Journal of Computational Physics 194 (2004) 363–393

[Link]/locate/jcp

Structural optimization using sensitivity analysis

and a level-set method q
a,*
Gr
egoire Allaire , Francßois Jouve a, Anca-Maria Toader b

a
Centre de Math
ematiques Appliqu
ees (UMR 7641), Ecole Polytechnique, 91128 Palaiseau, France
b
CMAF, Faculdade de Ci^encias da Universidade de Lisboa, Av. Prof. Gama Pinto 2, 1699 Lisboa, Portugal
Received 9 September 2002; received in revised form 2 July 2003; accepted 12 September 2003

Abstract

In the context of structural optimization we propose a new numerical method based on a combination of the
classical shape derivative and of the level-set method for front propagation. We implement this method in two and three
space dimensions for a model of linear or nonlinear elasticity. We consider various objective functions with weight and
perimeter constraints. The shape derivative is computed by an adjoint method. The cost of our numerical algorithm is
moderate since the shape is captured on a fixed Eulerian mesh. Although this method is not specifically designed
for topology optimization, it can easily handle topology changes. However, the resulting optimal shape is strongly
dependent on the initial guess.
Ó 2003 Elsevier B.V. All rights reserved.

Keywords: Shape optimization; Topology optimization; Sensitivity analysis; Shape derivative; Level-set

1. Introduction

Shape optimization of elastic structures is a very important and popular field. The classical method of
shape sensitivity (or boundary variation) has been much studied [24,27,31,33]. It is a very general method
which can handle any type of objective functions and structural models, but it has two main drawbacks: its
computational cost (because of remeshing) and its tendency to fall into local minima far away from global
ones. The homogenization method (and its variants, such as power-law materials or SIMP method, see e.g.
[1,2,5,9–11,19]) is an adequate remedy to these drawbacks but it is mainly restricted to linear elasticity and

q
This work has been supported by the grant CNRS/GRICES no. 2002-12163 of the Centre National de la Recherche Scientifique
(France) and the Gabinete de Relacßo ~ es Internacionais da Ci^encia e do Ensino Superio (Portugal).
*
Corresponding author. Tel.: +33-1-69-33-46-11; fax: +33-1-69-33-30-11.
E-mail addresses: [Link]@[Link] (G. Allaire), [Link]@[Link] (F. Jouve), amtan@[Link].
[Link] (A.-M. Toader).

0021-9991/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved.
doi:10.1016/[Link].2003.09.032
364 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

particular objective functions (compliance, eigenfrequency, or compliant mechanism). Recently yet another
method appeared in [25,29] based on the level-set method which has been devised by Osher and Sethian
[26,28] for numerically tracking fronts and free boundaries. The level-set method is versatile and compu-
tationally very efficient: it is by now a classical tool in many fields such as motion by mean curvature, fluid
mechanics, image processing, etc.
The work [25] studied a two-phase optimization of a membrane (modeled by a linear scalar partial
differential equation), i.e. the free boundary was the interface between two constituents occupying a given
domain. It combined the level-set method with the shape sensitivity analysis framework. On the other
hand, the work [29] focused on structural optimization within the context of two-dimensional linear
elasticity. The shape of the structure was the free boundary which was captured on a fixed mesh using the
immersed interface method. However, [29] did not rely on shape sensitivity analysis: rather the structural
rigidity was improved by using an ad hoc criteria based on the Von Mises equivalent stress. Another
relevant reference is [12] which proposed a phase-field model (close to the level-set method) for structural
optimization.
In this paper we generalize these two previous works in many aspects. We propose a systematic im-
plementation of the level-set method where the front velocity is derived from a shape sensitivity analysis.
We investigate different objective functions in two and three space dimensions: compliance (rigidity), least
square deviation from a target (compliant mechanism), design dependent loads (pressure loads). Other
objective functions like eigenfrequencies, multiple loads and minimal stress will be addressed in another
paper. We focus on shape optimization rather than two-phase optimization, and we replace the immersed
interface method by the simpler ‘‘ersatz material’’ approach which amounts to fill the holes by a weak
phase. This is a well-known approach in topology optimization which can be rigorously justified in some
cases [1]. We also consider the case of a nonlinear elasticity model. For all such problems we compute a
shape derivative by using an adjoint problem. Then, the shape derivative is used as the normal velocity of
the free boundary which is moved during the optimization process. Front propagation is performed by
solving a Hamilton–Jacobi equation for a level-set function. We study the effect of first-order or second-
order discretization and of reinitialization in the numerical convergence toward an optimal shape. We
also investigate the strong dependence of the computed optimal shape on the initial design. Our results
were announced in [4]. Related results have been independently obtained in [34].

2. Setting of the problem

We start by describing a model problem in linearized elasticity. There is no conceptual difficulty in

choosing another model, and in particular Section 8 deals with a nonlinear elasticity problem. Let X
Rd
(d ¼ 2 or 3) be a bounded open set occupied by a linear isotropic elastic material with HookeÕs law A. Recall
that, for any symmetric matrix n, A is defined by

An ¼ 2ln þ kðTr nÞId;

where l and k are the Lam
e moduli of the material. The boundary of X is made of two disjoint parts
oX ¼ CN [ CD ; ð1Þ
with Dirichlet boundary conditions on CD , and Neumann boundary conditions on CN . The two boundary
parts CD and CN are allowed to vary in the optimization process, although it is possible to fix some portion
of it (see the numerical examples below).
We denote by f the vector-valued function of the volume forces and by g that of the surface loads. The
displacement field u in X is the solution of the linearized elasticity system
 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 365
8
< divðAeðuÞÞ ¼ f in X;
u¼0 on CD ; ð2Þ
:
ðAeðuÞÞn ¼ g on CN :

Since X is varying during the optimization process, f and g must be known for all possible configurations of
X. We therefore introduce a working domain D (a bounded open set of Rd ) which contains all admissible
shapes X.
d d
To give a precise mathematical meaning to (2), we choose f 2 L2 ðDÞ and g 2 H 1 ðDÞ and we assume
that CD 6¼ ; (otherwise we should impose an equilibrium condition on f and g). In such a case it is well
d
known that (2) admits a unique solution in H 1 ðXÞ .
The objective function is denoted by J ðXÞ. In this paper, we shall mostly focus on two possible choices of
J (these are merely examples, and much more freedom is allowed). A first classical choice is the compliance
(the work done by the load)
Z Z Z
J1 ðXÞ ¼ f  u dx þ g  u ds ¼ AeðuÞ  eðuÞdx; ð3Þ
X CN X

which is very common in rigidity maximization. A second choice is a least square error compared to a target
displacement

Z 1=a
a
J2 ðXÞ ¼ kðxÞju  u0 j dx ; ð4Þ
X

which is a useful criterion for the design of compliant mechanisms [3,30]. We assume a P 2, u0 2 La ðDÞ and
k 2 L1 ðDÞ, a non-negative given weighting factor. In both formulas (3) and (4), u ¼ uðXÞ is the solution of
(2). We define a set of admissible shapes that must be open sets contained in the working domain D and of
fixed volume V
Uad ¼ fX
D such that jXj ¼ V g: ð5Þ
Our model problem of shape optimization is
inf J ðXÞ: ð6Þ
X2Uad

It is well known that the minimization problem (6) is usually not well posed on the set of admissible shapes
defined by (5) (i.e. it has no solution). In order to obtain existence of optimal shapes some smoothness or
geometrical or topological constraints are required. For example, a variant of (6) with a perimeter con-
straint turns out to be a well-posed problem [7]. The perimeter P ðXÞ of an open set X is defined as the
ðd  1Þ-dimensional
R Hausdorff measure of its boundary oX, i.e. P ðXÞ ¼ Hd1 ðoXÞ, which reduces to
P ðXÞ ¼ oX ds for smooth domains. Then, if ‘ > 0 is a positive Lagrange multiplier, the minimization
problem
inf ðJ ðXÞ þ ‘P ðXÞÞ ð7Þ
X2Uad

admits at least one optimal solution. There are other regularized variants of (6) which are well-posed and
we refer to [15,18] for such existence theories. Note that, even if existence is not an issue of the present
paper, we shall work with a smoother subset of (5) in order to define properly a notion of shape derivative.

Remark 1. We described our shape optimization problem with a single state equation, i.e. the shape is
optimized for a single set of loads. Our approach can easily be generalized to the more meaningful case of
shape optimization for multiple loads.
366 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

3. Shape derivative

In order to apply a gradient method to the minimization of (6) we recall a classical notion of shape
derivative. This notion goes back, at least, to Hadamard, and many have contributed to its development
(see e.g. the reference books [27,33]). Here, we follow the approach of Murat and Simon [24,31]. Starting
from a smooth reference open set X, we consider domains of the type
Xh ¼ ðId þ hÞðXÞ;
with h 2 W 1;1 ðRd ; Rd Þ. It is well known that, for sufficiently small h, ðId þ hÞ is a diffeomorphism in Rd .

Definition 2. The shape derivative of J ðXÞ at X is defined as the Fr
echet derivative in W 1;1 ðRd ; Rd Þ at 0 of
the application h ! J ððId þ hÞðXÞÞ, i.e.
joðhÞj
J ððId þ hÞðXÞÞ ¼ J ðXÞ þ J 0 ðXÞðhÞ þ oðhÞ with limh!0 ¼ 0;
khk
where J 0 ðXÞ is a continuous linear form on W 1;1 ðRd ; Rd Þ.

A classical result states that the directional derivative J 0 ðXÞðhÞ depends only on the normal trace h  n on
the boundary oX.

Lemma 3. Let X be a smooth bounded open set and J ðXÞ a differentiable function at X. Its derivative satisfies
J 0 ðXÞðh1 Þ ¼ J 0 ðXÞðh2 Þ
if h1 ; h2 2 W 1;1 ðRd ; Rd Þ are such that h2  h1 2 C1 ðRd ; Rd Þ and
h1  n ¼ h2  n on oX:

We give two examples of shape derivative that will be useful in the sequel.

Lemma 4. Let X be a smooth bounded open set and /ðxÞ 2 W 1;1 ðRd Þ. Define
Z
J ðXÞ ¼ /ðxÞdx:
X

Then J is differentiable at X and

Z Z
0
J ðXÞðhÞ ¼ divðhðxÞ/ðxÞÞdx ¼ hðxÞ  nðxÞ/ðxÞds
X oX

for any h 2 W 1;1 ðRd ; Rd Þ.

Lemma 5. Let X be a smooth bounded open set and /ðxÞ 2 W 2;1 ðRd Þ. Define
Z
J ðXÞ ¼ /ðxÞds:
oX

Then J is differentiable at X and

Z

0 o/
J ðXÞðhÞ ¼ hn þ H / ds;
oX on
 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 367

for any h 2 W 1;1 ðRd ; Rd Þ, where H is the mean curvature of oX defined by H ¼ div n. Furthermore, this result
still holds true if one replaces oX by C, a smooth open subset of oX, and assumes that / ¼ 0 on the surface
boundary oC.

Remark 6. In particular Lemma 4 is useful in order to compute the shape derivative of a volume constraint
V ðXÞ ¼ C. Indeed, we have
Z Z
V ðXÞ ¼ dx and V 0 ðXÞðhÞ ¼ hðxÞ  nðxÞds:
X oX

Similarly, Lemma 5 is useful in order to compute the shape derivative of a perimeter constraint P ðXÞ ¼ C.
Indeed, we have
Z Z
0
P ðXÞ ¼ ds and P ðXÞðhÞ ¼ hðxÞ  nðxÞH ds:
oX oX

Theorem 7. Let X be a smooth bounded open set and h 2 W 1;1 ðRd ; Rd Þ. Assume that the data f and g as well
d d d
as the solution u of (2) are smooth, say f 2 H 1 ðXÞ , g 2 H 2 ðXÞ , u 2 H 2 ðXÞ . The shape derivative of (3) is
Z
   Z
oðg  uÞ
J10 ðXÞðhÞ ¼ 2 þ Hg  u þ f  u  AeðuÞ  eðuÞ h  n ds þ AeðuÞ  eðuÞh  n ds: ð8Þ
CN on CD

The shape derivative of (4) is

Z

C0 a oðg  pÞ
J20 ðXÞðhÞ ¼ kju  u0 j þ AeðpÞ  eðuÞ  f  p   Hg  p h  n ds
CN a on
Z

C0 a
þ kju  u0 j  AeðuÞ  eðpÞ h  n ds; ð9Þ
CD a
d
where p is the adjoint state, assumed to be smooth, i.e. p 2 H 2 ðXÞ , defined as the solution of
8
< divðAeðpÞÞ ¼ C0 kðxÞju  u0 ja2 ðu  u0 Þ in X;
p¼0 on CD ; ð10Þ
:
ðAeðpÞÞn ¼ 0 on CN ;
where C0 is a constant given by

Z 1=a1
a
C0 ¼ kðxÞjuðxÞ  u0 ðxÞj dx :
X

Remark 8. There is no adjoint state involved in (8). Indeed the minimization of (3) is a self-adjoint problem
which turns out to be easier to solve than (4).

Proof. Although Theorem 7 is a classical result [24,27,31,33] we briefly sketch its proof for the sake of
completeness. To simplify we give a short, albeit formal, proof due to C
ea [17]. We consider a general
objective function
Z Z
J ðXÞ ¼ jðx; uðxÞÞdx þ lðx; uðxÞÞds;
X oX
368 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

2
for which we introduce the Lagrangian defined for ðv; qÞ 2 ðH 1 ðRd ; Rd ÞÞ by
Z Z Z Z
LðX; v; qÞ ¼ jðvÞdx þ lðvÞds þ AeðvÞ  eðqÞdx  q  f dx
X oX X X
Z Z
 q  g ds  ðq  AeðvÞn þ v  AeðqÞnÞds: ð11Þ
CN CD

In (11) q is a Lagrange multiplier for the state equation and its boundary conditions. It is worth no-
ticing that v and q belong to a functional space that does not depend on X, so we can apply the usual
differentiation rule to the Lagrangian L. The stationarity of the Lagrangian is going to give the op-
timality conditions of the minimization problem. For a given X, we denote by ðu; pÞ such a stationary
point. The partial derivative of L with respect to q, in the direction / 2 H 1 ðRd ; Rd Þ, after integration by
parts leads to
  Z Z
oL
ðX; u; pÞ; / ¼ 0 ¼  /  ðdivðAeðuÞÞ þ f Þdx þ /  ððAeðuÞÞn  gÞds
oq X CN
Z
 u  Aeð/Þn ds: ð12Þ
CD

Taking first / with compact support in X gives the state equation. Then, varying the trace function / on CN
gives the Neumann boundary condition for u, while varying the corresponding normal stress ðAeð/ÞÞn on
CD gives the Dirichlet boundary condition for u. On the other hand, in order to find the adjoint equation,
we differentiate L with respect to v in the direction / 2 H 1 ðRd ; Rd Þ. This yields
  Z Z Z Z
oL
ðX; u; pÞ; / ¼ 0 ¼ j0 ðuÞ  / dx þ l0 ðuÞ  / ds þ Aeð/Þ  eðpÞdx  ðp  Aeð/Þn
ov X oX X CD

þ /  AeðpÞnÞds:

Integrating by parts we obtain

  Z Z
oL
ðX; u; pÞ; / ¼ ðj0 ðuÞ  divðAeðpÞÞÞ  / dx þ /  ðAeðpÞn þ l0 ðuÞÞds
ov X CN
Z
þ ð/  l0 ðuÞ  p  Aeð/ÞnÞds:
CD

Taking first / with compact support in X gives the adjoint state equation
divðAeðpÞÞ ¼ j0 ðuÞ in X:
Then, varying the trace of / on CN yields the Neumann boundary condition

ðAeðpÞÞn ¼ l0 ðuÞ on CN :

Finally, varying the normal stress ðAeð/ÞÞn on CD gives
p¼0 on CD :
We have therefore found a well-posed boundary value problem for the adjoint state p.
The shape derivative of the objective function is obtained by differentiating

J ðXÞ ¼ LðX; uðXÞ; pðXÞÞ;

 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 369

which, by the chain rule theorem, reduces to the partial derivative of L with respect to X in the direction h
oL
J 0 ðXÞðhÞ ¼ ðX; u; pÞðhÞ:
oX
Applying Lemma 4 and 5 we obtain
Z Z

oL olðuÞ
ðX; u; pÞðhÞ ¼ h  nðjðuÞ þ AeðuÞ  eðpÞ  p  f Þds þ hn þ HlðuÞ ds
oX oX oX on
Z
 Z

oðg  pÞ oh
 hn þ Hg  p ds  hn þ Hh ds; ð13Þ
CN on CD on

with h ¼ u  AeðpÞn þ p  AeðuÞn. Taking into account the boundary condition u ¼ p ¼ 0 on CD which also
implies



ou op ou op
AeðuÞ  eðpÞ ¼ l  þ ðl þ kÞ n n on CD ;
on on on on
we deduce
Z

oL oðg  pÞ
ðX; u; pÞðhÞ ¼ h  n jðuÞ þ AeðuÞ  eðpÞ  p  f   Hg  p ds
oX CN on
Z Z

olðuÞ
þ h  nðjðuÞ  AeðuÞ  eðpÞÞds þ hn þ HlðuÞ ds:
CD oX on

This proof is merely a formal computation (in particular it assumes that u and p are differentiable with
respect to the shape X) but it can be rigorously justified (see the references quoted above). Of course, if the
objective function is the compliance, i.e. jðuÞ ¼ f  u, lðuÞ ¼ g  u on CN and lðuÞ ¼ 0 on CD , we find that
p ¼ u and the problem is self-adjoint.

We now give a variant of Theorem 7 when the surface loading is a pressure load which is oriented in the
direction of the normal vector. In other words, we replace Eq. (2) by
8
< divðAeðuÞÞ ¼ f in X;
u¼0 on CD ; ð14Þ
:
ðAeðuÞÞn ¼ p0 n on CN ;

where n is the unit normal vector and p0 ðxÞ is a given scalar function in H 2 ðDÞ.

Corollary 9. Let X be a smooth bounded open set and h 2 W 1;1 ðRd ; Rd Þ. Assume that the solution u of (14) is
d
smooth, say u 2 H 2 ðXÞ . The shape derivative of the compliance,
Z Z
J3 ðXÞ ¼ f  u dx þ p0 n  u ds
X CN

is
Z Z
J30 ðXÞðhÞ ¼ h  nð2f  u þ 2divðp0 uÞ  AeðuÞ  eðuÞÞds þ AeðuÞ  eðuÞh  n ds: ð15Þ
CN CD
370 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

Proof. We rewrite the objective function as

Z Z
J3 ðXÞ ¼ f  u dx þ divðp0 uÞdx:
X X
2
The Lagrangian of the problem is defined for ðv; qÞ 2 ðH 1 ðRd ; Rd ÞÞ by
Z Z Z Z Z
LðX; v; qÞ ¼ f  v dx þ divðp0 vÞdx þ AeðvÞ  eðqÞdx  q  f dx  divðp0 qÞdx
X X X X X
Z
 ðq  AeðvÞn þ v  AeðqÞnÞds: ð16Þ
CD

One can check that the adjoint state of the problem is p ¼ u (self-adjoint problem). A computation similar
to that of Theorem 7 shows that
oL
J30 ðXÞðhÞ ¼ ðX; u; uÞðhÞ;
oX
which yields the desired result (15).

Remark 10. We can generalize Theorem 7 to more general objective functions, including functions of the
strain or stress. It is also possible to consider non homogeneous Dirichlet boundary conditions in the state
equation. The case of a nonlinear model is discussed in Section 8.

Remark 11. It is possible to further restrict the class of domains by asking that some parts of the boundary
Cfixed do not move. In such a case, the vector field h must satisfy the constraint (or boundary condition)
h  n ¼ 0 on Cfixed .

We now have all the necessary theoretical ingredients to describe a gradient method for the minimization
of an objective function J ðXÞ. As we have just seen, the general form of its shape derivative is
Z
J 0 ðXÞðhÞ ¼ vh  n ds;
oX

where the function v is given by a result like Theorem 7. Ignoring smoothness issues, a descent direction is
found by defining a vector field
h ¼ vn; ð17Þ
and then we update the shape X as
Xt ¼ ðId þ thÞX;
where t > 0 is a small descent step. Formally, we obtain
Z
J ðXt Þ ¼ J ðXÞ  t v2 ds þ Oðt2 Þ;
oX

which guarantees the decrease of the objective function.

There are other possible choices for the definition of the descent direction. Let us first remark that, from
a mathematical point of view, formula (17) makes sense only if the resulting vector field h belongs to
W 1;1 ðRd ; Rd Þ. In view of typical definitions of v, (8) or (9), this is the case only if the state u and adjoint state
p, as well as the boundary of X, are smooth enough. It is not clear that the optimal shapes, if they exist, are
smooth. (However, a recent result of Chambolle and Larsen [16] proves that, in two space dimensions and
for a scalar problem, the optimal shape under a perimeter constraint has indeed C 1 regularity.)
 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 371

If either v or the normal n is not smooth, then it may be desirable to smooth the velocity field vn (this is a
classical issue in shape optimization; see e.g. [6,14], or Chapter 5 in [23]). The main idea is to change the
scalar product with respect to which we evaluate a descent direction. For example, working with the H01
scalar product instead of the L2 one, we need to solve

Dh ¼ 0 in X;
oh
on
¼ vn on oX:

In other words, we apply the Neumann-to-Dirichlet map to vn which has the effect of increasing of one
order the regularity of h on C (with respect to that of vn). Integrating by parts, we find
Z Z
jrhj2 dx ¼  vh  n ds;
X oX

which shows that h is a descent direction which guarantees again the decrease of J .

Remark 12. Another possibility is to use the Laplace–Beltrami operator DS on oX. We first compute a
1
regularization ~v ¼ ðDS Þ v and then take the descent direction h  n ¼ ~v. By integration by parts, we find
Z
0 2
J ðXÞðhÞ ¼  jrS ~vj ds;
oX

which clearly shows that h is a smoother descent direction.

4. Shape representation by the level-set method

As described above, the method of shape sensitivity can be (and has been) implemented in a Lagrangian
framework. It suffices to mesh X and to advect the mesh according to the descent direction h. However, this
implementation suffers at least from two drawbacks. First, if the shape is deformed too much, then it is
necessary to remesh which can be very costly (especially in 3-d). Second, different parts of the boundary of
the shape may want to merge or split, but as is well known topology changes are very difficult to handle
with such Lagrangian or front-tracking methods. Therefore, we favor an Eulerian approach and, following
[25,29], we use a level-set method to capture the shape X on a fixed mesh.
Let a bounded domain D
Rd be the working domain in which all admissible shapes X are included, i.e.
X
D. In numerical practice, the domain D will be uniformly meshed once and for all. We parameterize the
boundary of X by means of a level-set function, following the idea of Osher and Sethian [26]. We define this
level-set function w in D such that
8
< wðxÞ ¼ 0 () x 2 oX \ D;
wðxÞ < 0 () x 2 X;
:
wðxÞ > 0 () x 2 ðD n XÞ:

The normal n to the shape X is recovered as rw=jrwj and the curvature H is given by the divergence of the
normal div n (these quantities are evaluated by finite differences since our mesh is uniformly rectangular).
Remark that, although n and H are theoretically defined only on oX, the level-set method allows to define
easily their extension in the whole domain D (this will be useful in the sequel).
The elasticity equations for the state u (as well as the adjoint state p) are extended to the whole domain D
by using the so-called ‘‘ersatz material’’ approach. It amounts to fill the holes D n X by a weak phase
mimicking void but avoiding the singularity of the rigidity matrix. This is a well-known procedure in to-
pology optimization which can be rigorously justified in some cases [1]. Our method is thus simpler than the
372 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

immersed interface method proposed in [29]. More precisely, we define an elasticity tensor A ðxÞ which is a
mixture of A in X and of the weak material mimicking holes in D n X

 1 in X;
A ðxÞ ¼ qðxÞA with q ¼ ð18Þ
103 in D n X:

In numerical practice, q is piecewise constant in each cell and is adequately interpolated in those cells cut by
the zero level-set w ¼ 0 (the shape boundary). Note that, unlike the homogenization method (or any other
generalized material method), the ‘‘material density’’ q in (18) is almost always equal to its extreme values
(1 or 103 ) and the zone around the shape boundary where it takes intermediate values does not increase in
size during the optimization process.
To be more specific, let us consider a simple example for which there is no body force, f ¼ 0, and the
boundary oD of the working domain is decomposed in three parts
oD ¼ oDD [ oDN [ oD0 ;
such that oDD corresponds to Dirichlet boundary conditions, oDN to non-homogeneous Neumann
boundary conditions (surface loads g 6¼ 0), and oD0 to homogeneous Neumann boundary conditions
(traction-free), respectively. Recall the decomposition (1) of the shape boundary, oX ¼ CD [ CN . Admis-
sible shapes X are further constrained to satisfy
CD
oDD ; CN ¼ oDN [ C0 ; ð19Þ
where C0 supports an homogeneous Neumann boundary conditions (traction-free). In other words, the
surface loads g are applied only on a fixed subset of the boundary CN , while the boundary CD , with zero
displacement, must be a subset of the fixed boundary oDD . Consequently, the only optimized part of the
shape boundary is C0 which is traction free. These conditions are precisely those assumed in all numerical
examples of Section 6. Then, the displacement u is computed as the solution of
8
>
> divðA eðuÞÞ ¼ 0 in D;
>
<u ¼ 0 on oDD ;

ð20Þ
>
> ðA eðuÞÞn ¼ g on oDN ;
>
: 
ðA eðuÞÞn ¼ 0 on oD0 :

A similar boundary value problem holds for the adjoint p. The homogeneous Neumann boundary con-
dition on C0 is automatically taken into account in the weak formulation of (20), at least in the limit when
the ersatz material goes to zero. The case of body forces and surface loads on the optimized boundary is
treated in Section 7.
Following the optimization process, the shape is going to evolve according to a fictitious time which
corresponds to descent stepping (we shall come back to this issue in the next section). As is well-known, if
the shape is evolving in time, then the evolution of the level-set function is governed by a simple Hamilton–
Jacobi equation. To be precise, assume that the shape XðtÞ evolves in time t 2 Rþ with a normal velocity
V ðt; xÞ. Then

wðt; xðtÞÞ ¼ 0 for any xðtÞ 2 oXðtÞ:

Differentiating in t yields

ow ow
þ x_ ðtÞ  rw ¼ þ Vn  rw ¼ 0:
ot ot
 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 373

Since n ¼ rw=jrwj we obtain

ow
þ V jrwj ¼ 0:
ot
This Hamilton–Jacobi equation is posed in the whole box D, and not only on the boundary oX, if the
velocity V is known everywhere (as will be the case in the sequel). As is well-known [26], Hamilton–Jacobi
equations do not usually admit smooth solutions. Existence and uniqueness can be obtained instead in the
framework of viscosity solutions which allows for a convenient definition of generalized shape motion. In
numerical practice we actually compute such viscosity solutions by using an upwind scheme (see Section 5).

5. Optimization algorithm

For the minimization problem

inf J ðXÞ;
X2Uad

we computed a shape derivative

Z
0
J ðXÞðhÞ ¼ vh  n ds;
oX

where the function vðu; p; n; H Þ is given by a result like Theorem 7. Since n and H , as well as the state u and
the adjoint state p, are computed everywhere in D, the integrand v in the shape derivative is defined
throughout the domain D and not only on the free boundary oX. Therefore, we can define a descent di-
rection in the whole domain D by
h ¼ vn:
(It is also possible to regularize h but that does not change the sequel.) The normal component h  n ¼ v is
therefore the advection velocity in the Hamilton–Jacobi equation
ow
 vjrwj ¼ 0: ð21Þ
ot
Transporting w by (21) is equivalent to move the boundary of X (the zero level-set of w) along the descent
gradient direction J 0 ðXÞ. Our proposed algorithm is an iterative method, structured as follows:
1. Initialization of the level-set function w0 corresponding to an initial guess X0 .
2. Iteration until convergence, for k P 0:
(a) Computation of the state uk and adjoint state pk through two problems of linear elasticity posed in
Xk , approximated by (20).
(b) Deformation of the shape by solving the transport Hamilton–Jacobi equation (21). The new shape
Xkþ1 is characterized by the level-set function wkþ1 solution of (21) after a time step Dtk starting from
the initial condition wk ðxÞ with velocity vk computed in terms of uk and pk . The time step Dtk is cho-
sen such that J ðXkþ1 Þ 6 J ðXk Þ.
3. From time to time, for stability reasons, we also reinitialize the level-set function w by solving (22).
The Hamilton–Jacobi equation (21) is solved by an explicit first order upwind scheme [28] on a Cartesian
grid. In one space dimension, the scheme reads

wnþ1  wni n  n þ n  n
i
þ minðVi n ; 0Þg ðDþ n þ
x wi ; Dx wi Þ þ maxðVi ; 0Þg ðDx wi ; Dx wi Þ ¼ 0;
Dt
374 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

n wniþ1 wni n wn wn

with Dþ
x wi ¼ , D
Dx x wi ¼
i
Dx
i1
, and
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2 2
gþ ðd þ ; d  Þ ¼ minðd þ ; 0Þ þ maxðd  ; 0Þ ;

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
g ðd þ ; d  Þ ¼ maxðd þ ; 0Þ2 þ minðd  ; 0Þ2 :
We also implemented a second order scheme in order to improve accuracy (see Fig. 12 for some com-
parisons). The boundary conditions for w are of Neumann type. Since this scheme is explicit in time, its
time stepping must satisfy a CFL condition. Remark that the time step issued from this CFL condition is
usually much smaller than the time step Dtk which plays the role of the descent step in the minimization of
J ðXÞ. Remark also that one explicit time step for (21) is much cheaper, in terms of CPU time and memory
requirement, than the solution of the state equation (2) or adjoint state equation (10). Therefore, for each
iteration k in the above algorithm (corresponding to a single evaluation of uk and pk ), we perform several
explicit time steps of the Hamilton–Jacobi equation (21). The number of such time steps per iteration k is
monitored by the decrease of J ðXk Þ.
In practice, our algorithm never creates new holes or boundaries in 2-d if the Hamilton–Jacobi equation
(21) is solved under a strict CFL condition. Indeed, there are only two possibilities for creating a hole in 2-
d: either an initial hole splits in two parts (which is very unlikely when minimizing a stiffness criterion), or a
new hole nucleates inside a material region. This last mechanism is impossible because of the maximum
principle satisfied by the Hamilton–Jacobi equation (21). In other words, there is no nucleation mechanism
for new holes in our method. However the level-set method is well known to handle easily topology
changes, i.e. merging or cancellation of holes. Therefore, our algorithm is able to perform topology opti-
mization. In 2-d, the best results are obtained if the number of holes of the initial design is sufficiently large
(see Fig. 1). The situation is different in 3-d where new holes easily appear by pinching a thin wall (see
Fig. 25), and then the initial topology is less important.

Remark 13. One of the main advantage of the level-set method is that we never have to know where
precisely is the boundary oX. In particular, the same numerical scheme for solving the Hamilton–Jacobi
equation (21) is applied everywhere in the working domain D. Another advantage of the level-set method
comes from the simple formula to compute the normal n and the mean curvature H ¼ div n (which plays an
important role in perimeter penalization).

In order to regularize the level-set function (which may become too flat or too steep), we reinitialize it
periodically by solving
 ow
þ signðw0 Þðjrwj  1Þ ¼ 0 in D  Rþ ;
ot ð22Þ
wðt ¼ 0; xÞ ¼ w0 ðxÞ in D;

which admits as a stationary solution the signed distance to the initial interface fw0 ðxÞ ¼ 0g (see for ex-
ample Fig. 7). In numerical practice, reinitialization is very important because the level-set function often
becomes too steep which implies a bad approximation of the normal n or of the curvature H .

Fig. 1. Boundary conditions and two initializations of a 2-d cantilever.

 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 375

6. Numerical examples

In all computations we use a quadrangular mesh for both the level-set function and the elastic dis-
placement. We use Q1 finite elements for the elasticity analysis. All test cases have the following data, unless
otherwise specified. The Young modulus E of material A is normalized to 1 and the Poisson ratio m is fixed
to 0.3. The void or holes are mimicked by an ersatz material with the same Poisson ratio and Young
modulus 103 . For each elasticity analysis (that we shall call iteration in the sequel) we perform 20 explicit
time steps of the first-order scheme for the Hamilton–Jacobi transport equation. This number is auto-
matically reduced if the objective function is not decreasing. We also reinitialize the level-set function every
5 time steps of transport by performing five explicit time steps of Eq. (22).

6.1. 2-d cantilever

In the two-dimensional setting d ¼ 2 we first study a medium cantilever problem. The working domain is
a rectangle of size 2  1 discretized with a rectangular 80  40 mesh, with zero displacement boundary
condition on the left side and a unit vertical point load at the middle of the right side (see Fig. 1). Ad-
missible shapes must satisfy the constraint (19), i.e.
oX ¼ CD [ oDN [ C0 ;
where the shape Dirichlet boundary CD is a subset of the fixed boundary oDD , the surface loads g are
applied only the fixed boundary oDN , and C0 , which is traction free, is varying during the optimization
process. There are no body forces, i.e. f  0 in (2). The objective function is a combination of the com-
pliance and of the weight of the structure
Z Z
J ðXÞ ¼ g  u ds þ ‘ dx; ð23Þ
oDN X

where ‘ ¼ 100 is a fixed Lagrange multiplier for the weight constraint, and g ¼ 0 on C0 . The boundary
conditions and two initial configurations with different number of holes are displayed in Fig. 1. Under these
assumptions, the shape derivative of (23) is a special case of Theorem 7, namely
Z
0
J ðXÞðhÞ ¼ ð‘  AeðuÞ  eðuÞÞh  n ds; ð24Þ
C0

since h  n ¼ 0 on CD and on oDN where g 6¼ 0.

The algorithm converges smoothly to a (local) minimum which strongly depends, of course, on the initial
topology as can be checked in Figs. 2 and 3, displaying the optimal shapes as well as an intermediate result.
We run 50 iterations in order to show the good convergence and stability properties of our algorithm, but it
is clear that it has converged in a much smaller number of iterations (see Fig. 4). One can also check in
Fig. 5 that the L2 norm of the gradient of the objective function is decreasing, although it does not converge

Fig. 2. Iterations 10 and 50 of the two-dimensional cantilever initialized as in Fig. 1 (middle).

376 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

Fig. 3. Iterations 10 and 50 of the two-dimensional cantilever initialized as in Fig. 1 (right).

Figure 2
Figure 3

250
Objective function

200

150
5 10 15 20 25 30 35 40 45 50
Iterations

Fig. 4. Convergence of the objective function for the two-dimensional cantilever of Figs. 2 and 3.

0.035

0.03
L2 norm of the gradient

0.025

0.02

0.015

0.01

5 10 15 20 25 30 35 40 45 50
Iterations

Fig. 5. Convergence of the L2 norm of the gradient of the objective function for the two-dimensional cantilever of Fig. 2.
 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 377

to zero, due to numerical approximations. Our algorithm is just a steepest descent gradient method; of
course, its convergence can be speed up by using, for example, a quasi-Newton algorithm.
For a given initial guess our results are mesh-independent in 2-d. This is most probably due to the fact that
our method cannot create new holes in 2-d. For example, we performed the same cantilever problem on a
finer 160  80 mesh. Since we did not change any other parameters and because of the CFL condition for the
Hamilton–Jacobi transport equation, the time step is two times smaller and we need to run 100 iterations
instead of 50. One can check that the final shape in Fig. 6 is almost identical to the final one in Fig. 2.
The effect of reinitialization is best seen on fine meshes and is less important on coarse meshes. For the
previous cantilever example on a 160  80 mesh, we can see in Fig. 8 that no reinitialization yields a poorer
convergence than reinitializing every 5 time steps of transport (i.e. 4 times every elasticity iteration). Re-
mark in Fig. 7 the effect of reinitialization on the final level-set function w. On the other hand, Fig. 8 shows
that there is no clear difference on the objective function between a first-order or second-order scheme for
the Hamilton–Jacobi equations (both the transport equation (21) and the reinitialization equation (22)).
However, second-order accuracy improves greatly the efficiency of the reinitialization process (it is less
important on the transport equation). The optimal shapes for first or second order schemes are slightly
different, mainly near the boundary of the working domain D (compare Figs. 6 and 11).
As is well known, a necessary condition of optimality for a shape minimizing the objective function (23)
is that the ‘‘velocity’’ ð‘  AeðuÞ  eðuÞÞ is zero on the boundary. Fig. 9 displays the isocontours of this scalar
velocity. One can check that the velocity is approximately zero on the boundary.

6.2. Comparison with the homogenization method

We compare the level-set method with the homogenization method as described in [1,2]. For ‘‘good’’
initializations, the numerical result of our level-set method are similar to those obtained by the homoge-
nization method (see Chapter 5 in [1]). To make a precise comparison we study the medium cantilever, as in

Fig. 6. Iterations 50 and 100 of the two-dimensional cantilever on a finer 160  80 mesh initialized as in Fig. 1 (middle).

Fig. 7. Level-set function for the 2-d cantilever without reinitialization (left) and with reinitialization (right).
378 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

300
1st order
2nd order
no reinitialization

250
Objective function

200

150
0 10 20 30 40 50 60 70 80 90 100
Iterations

Fig. 8. Convergence of the objective function for the two-dimensional cantilever of Fig. 6 with various options.

Fig. 9. Scalar velocity ð‘  AeðuÞ  eðuÞÞ for the 2-d cantilever of Fig. 6.

the previous subsection, with a fine mesh (160  80). The Lagrange multiplier is now fixed to the value
‘ ¼ 150 and we use a second order scheme. Fig. 10 shows the optimal cantilever obtained by the ho-
mogenization method, while Fig. 11 displays two results of the level-set method corresponding to two
different initializations. The shapes are slightly different but the best result is obtained by homogenization
 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 379

Fig. 10. Cantilever computed with the homogenization method. Composite (left) and penalized solution (right).

Fig. 11. Level-set method: initializations 1 and 2 (left) and resulting optimal designs (right).

205
Homogenization
Level-set: initialization 1
Level-set: initialization 2
200
Objective function

195

190

185

180
0 50 100 150
Iterations

Fig. 12. Convergence of the objective function for the two-dimensional cantilevers of Figs. 10 and 11.

as can be checked in Fig. 12. Remark that the bump on the homogenization convergence curve in Fig. 12 is
due to the penalization process, while the small bump in the convergence curve of the level-set method for
initialization 2 is due to a change of topology (bar elimination).
380 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

Fig. 13. Boundary conditions and initialization of the bridge problem.

Fig. 14. Iterations 30 and 50 of the two-dimensional bridge initialized as in Fig. 13.

6.3. 2-d bridge

The next example in dimension d ¼ 2 is a bridge problem. The working domain is a rectangle of size
2  1.2 discretized with a rectangular 80  48 mesh, at the two lower corners the vertical displacement is
zero, and a unit vertical force is applied at the middle of the bottom side (see Fig. 13). We again impose that
the shape Dirichlet boundary CD is a subset of the fixed boundary oDD , that the surface loads g are applied
only the fixed boundary oDN , and that only the traction free boundary C0 is varying during the optimi-
zation process. We neglect body forces, i.e. f  0. The objective function is (23) and its shape derivative is
(24). The Lagrange multiplier for the weight constraint is ‘ ¼ 30. The optimal design is displayed in Fig. 14.
We run the same problem with a perimeter constraint, namely we minimize
Z Z Z
J ðXÞ ¼ g  u ds þ ‘ dx þ ‘0 ds; ð25Þ
oDN X oX

with ‘0 ¼ 1. Recall that the derivative of the perimeter is the curvature as stated in Remark 6. We clearly see
in Fig. 15 that the resulting optimal shape has fewer holes although the initialization was the same. The
performances of the shapes in Figs. 14 and 15 are close: the objective function (23) is equal to 72.9 for the
former and 73.1 for the latter (their weights are different however: 0.77 and 0.84, respectively).

6.4. 3-d examples

One advantage of the level-set method is its easy extension to three space dimensions. We therefore turn
to a first three-dimensional test case: the 3-d cantilever. The domain is a parallelepiped of dimensions
5  2.4  3 discretized with 50  24  30 elements. The right wall is fixed and a unit force is applied
downward on the middle of the left wall. Due to symmetry, the computation is performed on half of the
 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 381

Fig. 15. Iterations 50 and 200 of the two-dimensional bridge, initialized as in Fig. 13, with a perimeter constraint.

Fig. 16. Initial guess and optimal shape for the three-dimensional cantilever.

domain. We still minimize the weighted sum (23) of the compliance and the weight for a Lagrange mul-
tiplier ‘ ¼ 15. Fig. 16 shows the initial and optimal designs.
A second example is that proposed in [2]: the optimal electrical mast. The workspace is a T-like box. Two
symmetric vertical loads are applied in the middle of the lower edges of the horizontal part of the T and
represent the force exerted by the wires on the mast. Simply supported boundary conditions are imposed at
the corners of the base of the T. We still minimize the weighted sum (23) of the compliance and the weight
for a Lagrange multiplier ‘ ¼ 2. Only a quarter of the object is computed, by virtue of the symmetries. The
mesh of the T-box is made of a 24  12  28 bar and a 12  12  48 foot, for a total number of 14,976

Fig. 17. Initial guess and optimal shape for the three-dimensional electrical mast.
382 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

hexaedral elements. Starting with an initial design with regularly distributed holes, the resulting optimal
shape, displayed in Fig. 17, evokes the shape of actual electric masts.

6.5. Another objective function

We now give a numerical example for the minimization of the least square objective function (4). This is
a classical gripping mechanism test case which is described, e.g., in [1,30]. The working domain is a rect-
angle of size 5  4 with a rectangular hole of size 0.2  1.4 on the middle of the left side. By symmetry, only
the upper half of the domain is meshed with a rectangular 50  20 mesh.
In the objective function (4) the localizing factor kðxÞ is zero except on the black zone on the left side (the
jaws of the mechanism) where it is equal to 1, whereas the target displacement u0 ðxÞ is set to (100; 0) in the
upper left black zone (see Fig. 18). This objective function has been cooked up in order to obtain a gripping
mechanism for which the jaws close. From a mathematical point of view the main difference between this

input force

output force

Fig. 18. Boundary conditions for a plane gripping mechanism.

Fig. 19. Initialization, and deformed optimal shape of a plane gripping mechanism.
 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 383

test case and the previous ones is that it needs an adjoint system in order to evaluate the derivative. Fig. 19
shows the resulting optimal shape.
We perform the same optimization in 3-d. A uniform pressure load is applied on the left side of the box,
while the upper and lower sides are fixed. The initialization and the optimal deformed shape are displayed
in Fig. 20. One can see in Fig. 21 that the algorithm is stable and convergence occurs quite early during the
process.

Fig. 20. Initialization and deformed optimal configuration of a 3-d gripping mechanism.

1600

1500

1400
Objective function

1300

1200

1100

1000

900

800

5 10 15 20 25 30
Iterations

Fig. 21. Convergence of the objective function for the three-dimensional gripping mechanism of Fig. 20.
384 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

7. Application to design-dependent loads

Until now and in the whole Section 6 we neglected the body forces f and the surface loads g were applied
on a fixed part of the shape boundary oDN . In such a case the elasticity system was simply (20). We now
explain how to generalize it for more complicated loadings, including design-dependent loads as discussed
for example in [12].
In the case of body forces f , it suffices to extend them by 0 in D n X. This is easily done by multiplying f
by the characteristic function v defined in terms of the level-set function w by

0 if wðxÞ > 0;
vðxÞ ¼
1 if wðxÞ < 0:

Since our code discretizes the body forces as piecewise constants in each cell, the characteristic function v is
also constant in each cell. In the cells cut by the zero level-set we use the same interpolation procedure as
that used for the density q defined in (18).
Surface loads g which are applied on the moving shape boundary C0 (see (19) for this notation) are called
design-dependent since they depend on the position of the shape (recall that g is defined a priori everywhere
in D). Numerically, we shall replace these surface loads by equivalent volume forces using the fact that, in
the variational formulation of (2), they appear as
Z Z
g  u ds ¼ dC0 g  u dx;
C0 D

where dC0 ðxÞ is the Dirac mass function concentrated on C0 . Since

rðsignðwÞÞ ¼ 2dC0 n; ð26Þ

we have
1
2
jrðsignðwÞÞj ¼ dC 0 :

Introducing the following approximation of the sign function

wðxÞ
sðxÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ;
2
wðxÞ þ
2

where
> 0 is a small parameter chosen to regularize the jump over a few mesh elements. In all our com-
putations,
is equal to h=20 where h is the typical element size. This value spreads the jump in the sign
function over 2 cells in average. It allows classically to define an approximation dðxÞ of the Dirac function dC0
dðxÞ ¼ 12jrsðxÞj:

Finally, recalling definition (18) of the mixture A of the true and ersatz materials, the elasticity equations in
the working domain D are given by
8
>
> divðA eðuÞÞ ¼ vf þ dg in D;
<
u¼0 on oDD ;
 ð27Þ
>
> ðA eðuÞÞn ¼ g on oDN ;
: 
ðA eðuÞÞn ¼ 0 on oD0 :
 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 385

As an example we computed the optimal shape of a vertical mast submitted to a uniform horizontal surface
load g ¼ ð1; 0; 0Þ. This is a reinforcement problem for a column which is not subject to optimization. The
objective function is the compliance
Z Z
J ðXÞ ¼ g  u ds þ ‘ dx;
C0 X

with ‘ ¼ 100. The size of the working domain D is 4  1  4. The bottom of the box is fixed, i.e. equal to
oDD , the other sides are free, i.e. equal to oD0 , and oDN ¼ ;. In truth we applied g only where its scalar
product with the exterior normal n is negative: this is a very crude model of the effect of the wind. The
resulting optimal shape can be seen in Fig. 22.
Another type of design-dependent loads is the case of a pressure load p0 n where p0 is a given pressure
distribution. In other words, we now consider the elasticity system (14). Note that this load depends on
the normal to the boundary and therefore, if two different configurations of the free boundary C0 pass
through the same point, the loads at this point may be different for each configuration according to the
orientation of C0 . Once again, the surface loads give the following contribution in the variational for-
mulation
Z Z
p0 n  u ds ¼ p0 dC0 n  u dx;
C0 D

and from (26) we define an approximation of the function ðdC0 nÞ by

1
2
rs:

Fig. 22. Vertical mast reinforced to support a uniform horizontal load (from the left) modeling the effects of the wind.
386 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

Fig. 23. Optimal shape under a uniform pressure load with two anchor points (left). Applied forces on a zoomed area of the same
structure (right).

Fig. 24. Optimal shape under a uniform pressure load with five anchor points and a vertical force at center.

We therefore replace the surface pressure load by the approximate volume force p20 rs which is zero in the
whole domain D except in a neighborhood of C0 . One can see on the 2-d example of Fig. 23 that the
approximate volume force is really close to a surface pressure load.
The objective function is again the compliance
Z Z
J ðXÞ ¼ p0 n  u ds þ ‘ dx;
C0 X

with ‘ ¼ 100. Its shape derivative is computed by Corollary 9. The size of the working domain D is
4  4  1. In 3-d we impose five fixed (anchor) points on the bottom of the domain and a uniform pressure
load p0 ¼ 1 on the free boundary C0 . If there is an additional unit vertical load in the middle of the
 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 387

Fig. 25. Two views of the optimal shape under a uniform pressure load with five anchor points.

bottom side we obtain a nice starfish as can be seen of Fig. 24. Without this vertical load, we obtain a more
complex topology as can be checked in Fig. 25. Remark that the initialization in both cases was a convex
domain and that the level-set algorithm was able to create holes in the case of Fig. 25.

8. Generalization to nonlinear elasticity

8.1. Model and shape derivative

In this section we consider a nonlinear hyperelastic model [20]. If u is the displacement field,
F ¼ ðI þ ruÞ denotes the deformation gradient and the boundary value problem writes in the reference
configuration X (i.e. the undeformed configuration):
8
< divðT ðF ÞÞ ¼ f
> in X;
u¼0 on CD ; ð28Þ
>
:
T ðF Þn ¼ g on CN :

For the sake of simplicity, the volume forces f and the surface loads g are supposed to be dead loads, i.e.
independent of the displacement u.
If the material is hyperelastic, the first Piola–Kirchhoff stress tensor T derives from a potential W ðF Þ
(supposed smooth enough) called stored energy function

oW ðF Þ
Tij ¼ ; i; j 2 f1; . . . ; dg;
oFij

and for any deformation tensor F , the elasticity tensor AðuÞ (which is the tangential operator around F ) is
defined by

oTij ðF Þ o2 W ðF Þ
Aijkl ðuÞ ¼ ¼ :
oFkl oFij oFkl
388 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

The variational formulation of (28) is

Z Z Z
d
T ðF Þ  rv dx  f  v dx  g  v ds ¼ 0 8v 2 H 1 ðXÞ ; v ¼ 0 on CD ;
X X CN

and it is formally equivalent to a stationarity condition on the energy

Z Z Z
IðuÞ ¼ W ðI þ ruÞdx  f  u dx  g  u ds:
X X CN

Remark 14. The minimization problem

min IðvÞ
v

has a solution if W satisfies some convexity, growth and regularity conditions [8], but the question of
existence of solutions to the boundary value problem (28) is still open.

We give a variant of Theorem 7 for nonlinear elasticity and for the following objective function
Z Z
J ðXÞ ¼ jðx; uðxÞÞdx þ lðx; uðxÞÞds; ð29Þ
X oX

where j and l are smooth functions satisfying suitable growth conditions (so that J ðXÞ, as well as the
adjoint problem (30), make sense). In order to avoid unnecessary technical complications, we assume en-
ough smoothness of the data and existence and uniqueness of solutions in simple functional spaces.

Theorem 15. Let X be a smooth bounded open set and h 2 W 1;1 ðRd ; Rd Þ. Assume that the data f and g are
smooth ðf 2 H 1 ðXÞd , g 2 H 2 ðXÞd Þ and that there exists a unique smooth solution u 2 H 2 ðXÞd of (28). Assume
also that there exists a unique smooth solution p 2 H 2 ðXÞd of the following (linear) adjoint problem
8
< divðAðuÞrpÞ ¼ j0 ðuÞ in X;
p¼0 on CD ; ð30Þ
:
ðAðuÞrpÞn ¼ l0 ðuÞ on CN :

Assuming that AðuÞ and AðpÞ are coercive, the shape derivative of (29) is
Z Z

0 olðuÞ
J ðXÞðhÞ ¼ h  nð jðuÞ þ T ðI þ ruÞ  rp  p  f Þds þ hn þ HlðuÞ ds
oX on
Z
 Z
oX 
oðg  pÞ oh
 hn þ Hg  p ds  hn ds; ð31Þ
CN on CD on

where h ¼ u  T ðI þ rpÞn þ p  T ðI þ ruÞn.

Remark 16. In a nonlinear setting, the compliance is obtained by taking jðx; uÞ ¼ f  u and lðx; uÞ ¼ g  u in
(29). In such a case, the adjoint state p is different from ðuÞ.

Remark 17. The adjoint problem (30) is linear. The elasticity tensor AðuÞ involved is the tangential
operator around the solution u. If the numerical method used to solve the problem (28) is Newton-like,
this elasticity tensor has to be computed during the resolution of (28). Furthermore, if the linear systems
are solved by a direct solver (Cholesky factorization), then the additional cost of the adjoint problem is
very small since the rigidity matrix has already being factorized during the resolution of the direct
problem.
 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 389

Remark 18. Usually the stored elastic energy W ðF Þ is strongly convex near the identity, so that AðuÞ is
coercive if u is close to 0 (and thus, automatically, there exists a unique solution p to the adjoint
problem). Such a condition is met if, for example, the applied loads f and g are sufficiently small.
Furthermore, if j0 ð0Þ ¼ 0 and l0 ð0Þ ¼ 0, then u small implies that p is also small and thus AðpÞ is
coercive.

Proof. The (formal) proof is very similar to that of Theorem 7. We introduce the Lagrangian defined for
2
ðv; qÞ 2 ðH 1 ðRd ; Rd ÞÞ by
Z Z Z Z Z
LðX; v; qÞ ¼ jðvÞdx þ lðvÞds þ T ðI þ rvÞ  rq dx  q  f dx  q  g ds
X oX X X CN
Z
 ðq  T ðI þ rvÞn þ v  T ðI þ rqÞnÞds;
CD

and we investigate its stationarity conditions.

The partial derivative of L with respect to q, in the direction / 2 H 1 ðRd ; Rd Þ, after integration by parts,
leads to
  Z Z
oL
ðX; u; pÞ; / ¼ 0 ¼  /  ðdivðT ðI þ ruÞÞ þ f Þdx þ /  ððT ðI þ ruÞÞn  gÞds
oq X CN
Z
 u  ðAðpÞr/Þn ds:
CD

Taking first / with compact support in X gives the state equation of (28). Then, varying the trace function /
on CN gives the Neumann boundary condition for u, while varying the corresponding normal stress
ðAðpÞr/Þn on CD gives the Dirichlet boundary condition for u. (At this point we use the coercivity of AðpÞ
which is a sufficient condition for ðAðpÞr/Þn to span a dense subset of L2 ðCD Þd .)
Then we differentiate L with respect to v in the direction / 2 H 1 ðRd ; Rd Þ. This yields after integration by
parts
 
oL
ðX; u; pÞ; / ¼ 0
ov
Z Z
¼ ðj0 ðuÞ  divðAðuÞrpÞÞ  / dx þ /  ðAðuÞrpÞn þ l0 ðuÞ ds
X CN
Z
þ /  l0 ðuÞ þ AðuÞrpn  T ðI þ rpÞn  p  ðAðuÞr/Þn ds:
CD

Taking first / with compact support in X gives the adjoint state equation

divðAðuÞrpÞ ¼ j0 ðuÞ in X;

varying the trace of / on CN yields the Neumann boundary condition

ðAðuÞrpÞn ¼ l0 ðuÞ on CN ;
then varying the normal stress ðAðuÞr/Þn with / ¼ 0 on CD gives

p¼0 on CD :
390 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

Finally, just like in the linear case, we have

oL
J 0 ðXÞðhÞ ¼ ðX; u; pÞðhÞ;
oX
and the claimed expression of J 0 ðXÞðhÞ is obtained by applying Lemma 4 and 5 and remarking that h ¼ 0
on CD .

8.2. A numerical example

We come back to the medium cantilever problem described at the beginning of Section 6. The working
domain is a rectangle of size 2  1 discretized with a rectangular 160  80 mesh, with zero displacement
boundary condition on the left side and a downward vertical point load at the middle of the right side. The
initialization is that of Fig. 1 (middle). We use a Saint Venant–Kirchhoff material [20] with a Young
modulus equal to 1000 and a PoissonÕs coefficient equal to 0.3. This constitutive law is characterized by the
following relations:

F ¼ I þ ru; E ¼ 12ðF T F  IÞ;

R ¼ kTrðEÞI þ 2lE;

T ðF Þ ¼ F R:

This constitutive law does not satisfy any type of convexity (see Remark 14) so it is not clear that the
elasticity Eq. (28) are well posed. Nevertheless, we follow the usual practice in structural optimization [13]
and adopt this Saint Venant–Kirchhoff law. The intensity of the force is now varying between 1 and 3, while
the Lagrange multiplier for the weight constraint is adjusted to maintain a given global material proportion
of 40% of the total volume. We first perform 10 iterations in the setting of linear elasticity to get a rough

Fig. 26. Optimal shape of a cantilever in nonlinear elasticity for different intensities of the force compared to the optimal shape in
linear elasticity. Fine mesh (160  80), initialized as in Fig. 1 (middle).
 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 391

Fig. 27. Deformed configuration for the optimal cantilever in nonlinear elasticity: the intensity of the force is 3 and the maximum
displacement norm is equal to 0.4 (2  1 domain).

approximation of the shape, and then 100 in nonlinear elasticity. The resulting optimal shape in Fig. 26
displays thicker bars under compression than under traction as it should be, the disymmetry being am-
plified by the intensity of the force (Fig. 27).

9. Conclusion

We have proposed a method for shape and topology optimization in two and three space dimensions
which has the following advantages:
1. it allows for drastic topology changes during the optimization process;
2. its cost is moderate in terms of CPU time since this is an Eulerian shape capturing method;
3. it can handle very general objective functions and mechanical models, including nonlinear elasticity and
design-dependent loads;
4. with a good initialization it is as efficient as the homogenization method (when the latter one is avail-
able).
However, unlike the homogenization method [1,9,19], it is not a relaxation method, which means that
local minima have not been fully eliminated in favor of global minima. It is clear from the numerical ex-
amples that there still exist local (and non-global) minima to which the method may converge if the ini-
tialization is too far from a global minimum. In practice, this level-set method behaves differently in 2-d and
3-d. In 2-d the method works as if it can reduce the topology (i.e. the number of holes) but it cannot create
new holes since there is no nucleation mechanism in our algorithm. Therefore, one must be careful in the
choice of the initialization that should contain a large number of holes if one seek a non-trivial topology. In
3-d, there are less topological restrictions and the algorithm can easily create holes or non-simply connected
shapes even if the initial guess is convex. The mechanism of 3-d topology changes is the possible crossing of
two separate zero level-set surfaces without breaking the connectivity of the shape or of the void. Fre-
quently the final shape has a much more complex topology than the initial shape. Nevertheless, different
initializations can still drive the algorithm to different optimal shapes in 3-d.
The problem of choosing an adequate initialization can be managed in two ways, at least. First, one can
run the homogenization method on a simplified problem (say, linear elasticity with compliance objective
function) as a pre-processor in order to find a correct initial topology. Second, one can couple our method
with the so-called bubble method, or topological derivative, proposed by [21,22,32], which is precisely a
criteria for hole nucleation.
392 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393

As any shape-capturing algorithm our level-set method cannot attain high precision for the shape
boundary. This is not a problem for topology optimization since these methods are usually post-processed
by a classical shape optimization method based on mesh deformation. Finally, the present level-set method
of shape and topology optimization can easily (in principle) be extended to multi-physics applications.

References

[1] G. Allaire, Shape Optimization by the Homogenization Method, Springer Verlag, New York, 2001.
[2] G. Allaire, E. Bonnetier, G. Francfort, F. Jouve, Shape optimization by the homogenization method, N€ umerische Mathematik 76
(1997) 27–68.
[3] G. Allaire, F. Jouve, Optimal design of micro-mechanisms by the homogenization method, Eur. J. Finite Elements 11 (2002) 405–
416.
[4] G. Allaire, F. Jouve, A.-M. Toader, A level-set method for shape optimization, C.R. Acad. Sci. Paris, S
erie I 334 (2002) 1125–
1130.
[5] G. Allaire, R.V. Kohn, Optimal design for minimum weight and compliance in plane stress using extremal microstructures, Eur. J.
Mech. A/Solids 12 (6) (1993) 839–878.
[6] F. Almgren, J. Taylor, Optimal geometry in equilibrium and growth. Symposium in Honor of Benoit Mandelbrot (Curacßao,
1995), Fractals 3 (1995) 713–723.
[7] L. Ambrosio, G. Buttazzo, An optimal design problem with perimeter penalization, Calc. Var. 1 (1993) 55–69.
[8] J.M. Ball, Convexity conditions and existence theorems in nonlinear elasticity, Arch. Rational Mech. Anal. 63 (1977)
337–403.
[9] M. Bendsoe, Methods for Optimization of Structural Topology Shape, and Material, Springer Verlag, New York, 1995.
[10] M. Bendsoe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comp. Meth. Appl.
Mech. Eng. 71 (1988) 197–224.
[11] M. Bendsoe, O. Sigmund, Topology Optimization. Theory, Methods, and Applications, Springer Verlag, New York, 2003.
[12] B. Bourdin, A. Chambolle, Design-dependent loads in topology optimization, ESAIM Control Optim. Calc. Var. 9 (2003) 19–48.
[13] T. Buhl, C. Pedersen, O. Sigmund, Stiffness design of geometrically nonlinear structures using topology optimization, Struct.
Multidisciplinary Optim. 19 (2) (2000) 93–104.
[14] M. Burger, A framework for the construction of level set methods for shape optimization and reconstruction, Interf. Free
Boundaries 5 (2003) 301–329.
[15] A. Chambolle, A density result in two-dimensional linearized elasticity and applications, Arch. Rational Mech. Anal. 167 (2003)
211–233.
[16] A. Chambolle, C.J. Larsen, C 1 -regularity of the free boundary for a two-dimensional optimal compliance problem, Calc. Var.
Partial Diff. Eq., in press.
[17] J. C
ea, Conception optimale ou identification de formes, calcul rapide de la d
eriv
ee directionnelle de la fonction co^
ut, Math.
Model. Numer. Anal. 20 (3) (1986) 371–402.
[18] D. Chenais, On the existence of a solution in a domain identification problem, J. Math. Anal. Appl. 52 (1975) 189–289.
[19] A. Cherkaev, Variational Methods for Structural Optimization, Springer Verlag, New York, 2000.
[20] P.G. Ciarlet, Mathematical Elasticity, Three-Dimensional Elasticity, vol. I, North-Holland, Amsterdam, 1988.
[21] H. Eschenauer, A. Schumacher, Bubble method for topology and shape optimization of structures, Struct. Optim. 8 (1994) 42–51.
[22] S. Garreau, P. Guillaume, M. Masmoudi, The topological asymptotic for PDE systems: the elasticity case, SIAM J. Control
Optim. 39 (6) (2001) 1756–1778.
[23] B. Mohammadi, O. Pironneau, Applied Shape Optimization for Fluids, Clarendon Press, Oxford, 2001.
[24] F. Murat, S. Simon, Etudes de problemes dÕoptimal design. Lecture Notes in Computer Science 41, Springer Verlag, Berlin, 1976,
pp. 54–62.
[25] S. Osher, F. Santosa, Level-set methods for optimization problems involving geometry and constraints: frequencies of a two-
density inhomogeneous drum, J. Comp. Phys. 171 (2001) 272–288.
[26] S. Osher, J.A. Sethian, Front propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations, J.
Comp. Phys. 78 (1988) 12–49.
[27] O. Pironneau, Optimal Shape Design for Elliptic Systems, Springer-Verlag, New York, 1984.
[28] J.A. Sethian, Level-Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics,
Computer Vision and Materials Science, Cambridge University Press, Cambridge, MA, 1999.
[29] J. Sethian, A. Wiegmann, Structural boundary design via level-set and immersed interface methods, J. Comp. Phys. 163 (2000)
489–528.
[30] O. Sigmund, On the design of compliant mechanisms using topology optimization, Mech. Struct. Mach. 25 (1997) 493–524.
 G. Allaire et al. / Journal of Computational Physics 194 (2004) 363–393 393

[31] J. Simon, Differentiation with respect to the domain in boundary value problems, Numer. Funct. Anal. Optim. 2 (1980) 649–687.
_
[32] J. Sokołowski, A. Zochowski, Topological derivatives of shape functionals for elasticity systems, Mech. Struct. Mach. 29 (3)
(2001) 331–349.
[33] J. Sokołowski, J.P. Zolesio, Introduction to shape optimization: shape sensitivity analysis, Springer Series in Computational
Mathematics, vol. 10, Springer, Berlin, 1992.
[34] M.Y. Wang, X. Wang, D. Guo, A level-set method for structural topology optimization, Comput. Meth. Appl. Mech. Engrg. 192
(2003) 227–246.