BOOK2019 Wyklad

Guenter Leugering, Jan Sokolowski and Antoni Zochowski
The Shape Derivative Method. Theory and

Applications
January 2019
Springer
Contents
Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
1 Preliminaries and notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
2 Shape optimization for problems on metric graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.1 Problems on graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 Planar graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2.1 Steklov-Poincaré map for metric graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.2 Problems on general metric graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.3 Shape variations for graph problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.4 Shape sensitivity for the first eigenvalue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2.5 Transmission problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2.6 Topological derivative with respect to material properties . . . . . . . . . . . . . . . . . . . . . . . 23
2.2.7 Shape-topological derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3 Stars with a hole . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.3.1 Homogeneous networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.4 The topological derivative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.4.1 Homogeneous graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.4.2 Sensitivity with respect to edge inclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3 Mixed formulation of elliptic boundary value problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2 Mixed formulation in one space dimension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.3 Generalized Lax-Milgram Lemma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.4 Mixed formulation of elliptic boundary value problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.4.1 Stokes problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.4.2 Linearized elasticity (Navier or Lamé problem) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.5 Piola transformation for finite elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.6 Piola transformation for Stokes problem and the shape derivative calculations . . . . . . . . . . . . 45
3.7 Mixed formulation and finite elements for elasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4 Lagrangian formalism for weak solutions and shape optimization . . . . . . . . . . . . . . . . . . . . . . . 47

4.1 Lagrangian formalism in shape optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.1.1 Distributed shape gradients for Dirichlet problem by velocity method . . . . . . . . . . . . 48
4.1.2 Shape optimization for homogeneous Dirichlet problem in Ω ⊂ R2 . . . . . . . . . . . . . . 49
5 The second-order shape differentiability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.2 Shape Hessian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.3 Numerical evaluation of the shape Hessian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
v
vi Contents
6 Shape optimization with MATLAB in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

6.1 Shape optimization with Poisson equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.1.1 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.1.2 Matlab Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
6.2 Interface problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
6.2.1 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
6.2.2 Matlab Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
6.3 Shape optimization with linear elasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.3.1 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.3.2 Matlab Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
6.4 Example Shape Hessian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.4.1 Setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.4.2 Adjoint equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.4.3 Shape Gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.4.4 Shape Hessian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
6.4.5 Matlab program . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.5 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
A Program Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
A.0.1 Shared functions for Optimization Examples Chapter (??) . . . . . . . . . . . . . . . . . . . . . . 93
Preface
We present the shape optimization problems in the framework of the shape derivative method. First, the
elementary problems are solved in Matlab using this technique. In the second part we propose some research
project for students interested in shape optimization.
This is a book for postgraduated and PhD students. The contents and the purposes the book can be briefly
described as follows.
• Modeling using Partial Differential Equations (PDEs) in variables domains;
• Local variations of solutions to PDEs on variable domains by means of its Shape Derivatives;
• Simple gradient or Newton techniques of numerical solution for shape optimization problems;
• Applied sensitivity analysis for 1D and 2D boundary value problems including:
– Mathematical sensitivity analysis of Elliptic Boundary Value Problems;
– Matlab programs for the shape sensitivity analysis in 1D and 2D;
– Topological sensitivity analysis for graphs;
• Shape sensitivity analysis for 2D Elliptic Boundary Value Problems (EBP): Laplace, Stokes, Elasticity;
• Research seminar problems for compressible Navier-Stokes stationary problems.
In another words a part of the book can be used for the introductory course on Shape Optimization. There
are Matlab programs for selected shape optimization problems in 1D and 2D. This part can be used for the
laboratory on numerical methods for variatonal problems using e.g., the finite element method. The most
advanced part can be used for the Master Projects on the theory and applications of shape optimization.
There is a link of our study to the book on the topological derivative method [?] since the topological
derivative is obtained as a singular limit of shape derivatives. We present an original contribution to the
theory of topological derivatives on graphs. The topological derivatives of solutions to the discrete structures
on graphs are obtained by nonsingular perturbation technique.
Chapter 1
Preliminaries and notation
Chapter 2
Shape optimization for problems on metric graphs
2.1 Problems on graphs
Differential equations on graphs or, more precisely, on metric graphs are ubiquitous. Transportation net-
works conveying fluids, gas, electrical power or traffic can be modelled by hyperbolic systems on networks.
Mechanical networks involving strings, cables and beams are important in civil engineering. Heat flow in
technical devices as well as in buildings can be regarded as parabolic equations on networks. The underly-
ing quasi-static model or one that is obtained after suitable time-discretization is that of an elliptic system
on networks.
2.2 Planar graphs
We consider a simple planar graph (V, E) = G in R2 , with vertices V = {vJ |J ∈ J } and edges E = {ei |i ∈
I }. Let m = |J |, n = kI k be the numbers of vertices and edges, respectively. In general the edge-set
may be a collection of smooth curves in R2 , parametrized by their arc lengths. The restriction to planar
graphs and straight edges is for the sake of simplicity only. The more general case, which is of course also
interesting in the combination of shape and topology optimization, can also be handled. However this is
beyond these notes.
We associate to the edge ei the unit vector ei aligned along the edge. e⊥ i denotes the orthogonal unit
vector. Given a node vJ we define
IJ := {i ∈ I |ei is incident at vJ }
the incidence set, and dJ = |IJ | the edge degree of vJ . The set of nodes splits into simple nodes JS and
multiple nodes JM according to dJ = 1 and dJ > 1, respectively. On G we consider a vector-valued function
r representative of the displacement of the network (see Figure 2.1)
pi
y : G → Rnp := Πi∈I R, pi ≥ 1 ∀i. (2.2.1)
The numbers pi represent the degrees of freedom of the physical model used to describe the behavior of
the edge with number i. For instance, p = 1 is representative of a heat problem, whereas p = 2, 3 is used
in an elasticity context on graphs in 2 or 3 dimensions. The p0i s may change in the network in principle.
However, in this paper we insist on pi = p = 2, ∀i. See Lagnese, Leugering and Schmidt[?] and Lagnese
and Leugering [210] for details on the modeling.
Once the function r is understood as being representative of, say, a deformation of the graph, we may
localize it to the edges
yi := y|ei : [αi , βi ] → R p , i ∈ I , (2.2.2)
where ei is parametrized by x ∈ [αi , βi ] =: Ii ,0 ≤ αi < βi , ì := βi − αi . See Figure 2.1
We introduce the incidence relation
(
1 if ei ends at vJ
diJ :=
−1 if ei starts at vJ
3
4 2 Shape optimization for problems on metric graphs
X
Fig. 2.1 Representation of planar displacement
Accordingly, we define (
0 if diJ = −1
xiJ :=
ì if diJ = 1
We will use the notation yi (vJ ) instead of yi (xiJ ). In order to represent the material considered on the graph,
we introduce stiffness matrices
1 1
Ki := hi [(1 − )I + ei eTi ] (2.2.3)
si si
Obviously, the longitudinal stiffness is given by hi , whereas the transverse stiffness is given by hi (1 −
1
si ). This can be related to 1-d analoga of the Lamé parameters. We introduce Dirichlet and Neumann
simple nodes as follows. As the displacements and, consequently, the forces are vectorial quantities, we
may consider nodes, where the longitudinal (or tangential) displacement or forces are kept zero, while the
transverse displacements of forces are not constrained, and the other way round. We thus define
JDt := {J ∈ JS |yi (vJ ) · ei = 0}
JDn := {J ∈ JS |yi (vJ ) · e⊥

i = 0}
d i
JNt := {J ∈ JS |diJ Ki y (vJ ) · ei = 0}
dx
d
JNn := {J ∈ JS |diJ Ki yi (vJ ) · e⊥ i = 0}
dx
Notice that these sets are not necessarily disjoint. Obviously, the set of completely clamped vertices can be
expressed as
JD0 := JDt ∩ JDn (2.2.4)
Similarly, a vertex with completely homogenous Neumann conditions is expressed as JNn ∩ JNt . At tan-
gential Dirchlet nodes in JDt we may, however, consider normal Neumann-conditions as in JNn and so
on. In particular, in this paper we will consider bilateral contact conditions for the displacements at simple
Dirichlet nodes, see Figure 2.2. For the sake of simplicity we concentrate on such obstacles with respect to
the transverse displacement only.
J̃Dc := {J ∈ JS |yi (vJ ) · e⊥

i ∈ [ai , bi ]}, (2.2.5)
where ai ≤ bi for all i ∈ ID , D ∈ JDc .

We may then consider bilaterally constrained vertices where the the tangential force is zero, i.e. JDc ∩
JNt or those where with zero longitudinal displacement, i.e. JDc ∩ JDt . The most general treatment would
obscure the presentation, and we thus restrict ourselves to the latter case. Thus, we always assume that a
simple vertex under bilateral constraints admits only zero tangential forces. We may therefore define

⊥ d i
JD := J ∈ JS |y (vJ ) · ei ∈ [ai , bi ], diJ Ki y (vJ ) · ei = 0,
c i
(2.2.6)
dx
2.2 Planar graphs 5
Fig. 2.2 A general graph
We do not consider constraints around multiple joints which would restrict the motion of such a joint,
say, to a box. Again, the more general situation can be handled with the analysis presented here. The basic
assumption at a multiple node is that the deformation y is continuous across the joint. In truss design this
is not the case, and consequently pin-joints are considered, however on a discrete level. One may consider
pin-joints for networks of beams also on the continuous level, as in Lagnese, Leugering and Schmidt[?] and
[210]. In this paper we restrict ourselves to ’rigid’ joints in the sense that the angles between edges in their
reference configuration remain fixed. The continuity is expressed simply as
yi (vJ ) = y j (vJ ), i, j ∈ IJ , J ∈ JM
We consider the energy of the system
Zì
1 d i d i
E0 := ∑ Ki y · y + ci yi · yi dx (2.2.7)
2 i∈I dx dx
0
where the primes denote the derivative with respect to the running variable xi , ci represents an additional
spring stiffness term or an elastic support.
In order to analyze the problem, we need to introduce a proper energy space
V := {y : G → Rnp |yi ∈ H 1 (Ii ) (2.2.8)

yi (vD ) = 0, i ∈ ID , D ∈ JD0 (2.2.9)
y (vJ ) = y (vJ ), ∀i, j ∈ IJ , J ∈ JM }
i j
(2.2.10)
V is clearly a Hilbert space in

H := L2 (0, ì )np (2.2.11)
We introduce the bilinear form on V × V
Zì
d i d i
a(r, ϕ) := ∑ [Ki y · ϕ + ci yi · ϕ i ]dx. (2.2.12)
i∈I 0 dx dx
Let now distributed and boundary forces, f i , gJ be given along the edge ei and at the node vJ , respectively,
which define a continuous linear functional in V
Zì
L(ϕ) := ∑ f i · ϕ i dx + ∑ g gJ · ϕ î (vJ ), (2.2.13)
i∈I 0 f J∈JN
where î indicates that the simple nodes have just one incident edge, and where f i ∈ H 1 (0, ì )∗ . We now
consider minimizing the energy over the set of constrained displacements. To this end we introduce the
convex and closed (and hence weakly closed) set
K := V ∩ {(yi )ni=1 |yi (vD ) · e⊥

i ∈ [ai , ϕ ], i ∈ ID , D ∈ JD }
î c
(2.2.14)
The Ritz-approach to deriving the problem now can be stated as follows

1
min a(y, y) − L(y) (2.2.15)
y∈K 2
That this convex optimization problem admits a unique solution is then proved by standard arguments. The
classical first order necessary optimality conditions then read as follows.
n Zì n Z ì
d i d i d i
∑ [Ki y · ( ŷ − y ) + ci yi · (ŷi − yi )]dx − ∑ f i · (ŷi − yi )dx ≥ 0, ∀ŷ ∈ K (2.2.16)
i=1 dx dx dx i=1 0
0
In order to explore this variational inequality, we introduce active and inactive sets with respect to the
bilateral obstacles.
A u := {i|i ∈ ID , D ∈ JDc , yi (vD ) · e⊥
i = bi }
A ` := {i|i ∈ ID , D ∈ JDc , yi (vD ) · e⊥
i = ai } (2.2.17)
A 0 := {i|i ∈ ID , D ∈ JDc , ai < yi (vD ) · e⊥
i < ai }
In order to define proper variations in (2.2.16), we introduce the Hilbert space
V 0 = {ϕ ∈ V |ϕ i (vD ) · e⊥
i = 0, i ∈ ID , D ∈ JD }
c
(2.2.18)
Obviously, if r ∈ K then ŷ = y + ϕ ∈ K , ∀ϕ ∈ V 0 . Taking these variations we obtain from (2.2.16) the

following variational equation
n Zì
d i d i
∑ [Ki y · ϕ + ci yi · ϕ i − f i · ϕ i ]dx = 0, ∀ϕ ∈ V 0 (2.2.19)
i=1 dx dx
0
This variational problem, in turn, can be further analyzed by integration by parts (if additional H 2 -regularity
holds) in order to obtain
d
∑ c dîJ [Ki dx yi (vJ ) · ei ][ϕ i (vJ ) · ei ]
J∈JD
d i
+ ∑ dîJ Ki y (vJ )ϕî (vJ ) − ∑ gJ · ϕî (vJ )
J∈JN dx J∈J N
d (2.2.20)
+ ∑ ∑ diJ Ki yi (vJ ) · ϕ i (vJ )
J∈JM i∈IJ dx
n Zì
d2 i
+∑ {−Ki y + ci yi − f i } · ϕ i dx = 0, ∀ϕ ∈ V 0
i=1 dx2
0
Now (2.2.20) clearly implies the strong formulation of the problem:
d2 i
−Ki y + ci yi = f i in (0, ì )
dx2
d
∑ ∑ diJ Ki dx yi (vJ ) = 0, J ∈ JM
J∈JM i∈I J (2.2.21)
d
diN Ki yi (vN ) = gN , i ∈ IN , N ∈ JN
dx
d
Ki yi (vD ) · ei = 0, i ∈ ID , D ∈ JDc
dx
2.2 Planar graphs 7
We now concentrate on the active and inactive sets. We may take variations in (2.2.16) as follows:
ŷi = yi + ψi , ψ ∈ V , ψi (vD ) = 0, i ∈ ID , D ∈ JDc

i ∈ A u : ψi (vD ) · e⊥
i ≤0 (2.2.22)
i ∈ A l : ψi (vD ) · e⊥
i ≥0
i ∈ A o : ψi (vD ) · e⊥
i = ±ε, ε small
d i
Obviously, taking variations in the inactive case, we obtain Ki dx y (vD )·e⊥
i = 0 which together with (2.2.21)4
gives
d
Ki yi (vD ) = 0 i ∈ A o (2.2.23)
dx
In the active cases we get
d i
diD Ki y (vD ) · e⊥ ⊥
i [ψi (vD ) · ei ] ≥ 0, i ∈ A ∪ A
u l
(2.2.24)
dx
and hence
d i
diD Ki y (vD ) · e⊥
i ≤ 0, i ∈ A
u
dx (2.2.25)
d
diD Ki yi (vD ) · e⊥ i ≥ 0, i ∈ A
l
dx
Putting all together ((2.2.21),(2.2.25), (2.2.23) and the conditions involved in V ) we obtain the strong
formulation of (2.2.16)
d2 i


 −K i y + ci yi = f i , ∀i ∈ I
dx2




yi (vD ) = 0, i ∈ ID , D ∈ JD





d


diD Ki yi (vN ) · ei = 0, i ∈ ID , D ∈ JDc






 dx
 d
 diN Ki yi (vN ) = gJ , i ∈ IN , N ∈ JN




 dx
) = y j (vJ ), ∀i, j ∈ IJ , J ∈ JM
 i
y (v

J



d (2.2.26)
 ∑ diJ Ki yi (vJ ) = 0, J ∈ JM



 i∈IJ dx

ai ≤ yi (vD ) · e⊥ i ≤ bi , i ∈ ID , D ∈ JD
c






d

Ki yi (vD ) · e⊥ i =0 i∈A

 o




 dx
d


diD Ki yi (vD ) · e⊥ i ≤ 0, i ∈ A
u





 dx
d


diD Ki yi (vD ) · e⊥ i ≥ 0, i ∈ A
l


dx
where f i = 0, i ∈ I \ I f , gN = 0, J ∈ JN \ JNg . Notice that (2.2.26) line 6 is an example of the classical
Kirchhoff condition known from electrostatics. See Lagnese, Leugering and Schmidt [?],[210] for the case
without obstacles.
In order to follow this cycle of ideas bottom up, we first consider very simple such equations on networks,
in fact on star-graphs. We provide some first hand information before we embark on shape optimization
problems related to problems on metric graphs. We consider the following star-graph with m edges, m
external nodes at x = ì and one multiple node at x = 0. On this graph we consider the Laplacean with
Dirichlet conditions.
d2 i
− y + ci yi = f i in (0, ì )
dx2
yi (ì ) = ui
(2.2.27)
yi (0) = y j (0), i 6= j = 1, . . . , m
m
d
∑ dx yi (0) = 0
i=1
x = `1
x=0
x = `3 x = `2
Fig. 2.3 The three-star
Example 2.2.1 Let f i = 0, ci = 0, i = 1, . . . , m then

d i
yi (x) = ai x + bi , y (x) = ai , i = 1, . . . , m (2.2.28a)
dx
yi (ì ) = ai ì + bi = ui , i = 1, . . . , m (2.2.28b)
i j
y (0) = y (0), i 6= j = 1, . . . , m ⇔ bi = b j = b, i = 1, . . . , m (2.2.28c)
m m
d
∑ dx yi (0) = 0 ⇔ ∑ ai = 0 (2.2.28d)
i=1 i=1
Samstag, 15. Dezember 12
But
1 i
ai = (u − b) (2.2.28e)
ì
and hence
m m m
1 1
∑ ai = ∑ ì ui − b ∑ ì = 0 (2.2.28f)
i=1 i=1 i=1
m
1 1
b= m
1 i=1
∑ ì ui (2.2.28g)
∑ ì
i=1
This implies
 
1 i 1 1
ai = u − ∑ m ` j u j (2.2.28h)

ì ∑ m `1j j=1
j=1
Hence, the explicit solution is given by

 
m
1 i 1 1 1 1
yi (x) = u − ∑ m ` j u j x + ∑ ì ui (2.2.29)

m
ì ∑ m `1j j=1 ∑ 1 i=1
j=1 ì
i=1
Example 2.2.2 Let f i = 0, ci > 0, i = 1, . . . , m. Then the solutions are given by

√ √
yi (x) = ai sinh( ci x) + bi cosh( ci x), i = 1, . . . , m. (2.2.30)
The boundary and node conditions give

2.2 Planar graphs 9
√ √
yi (ì ) = ai , sinh( ci ì ) + bi cosh( ci ì ) = ui (2.2.31a)
i j
y (0) = y (0), i 6= j = 1, . . . , m ⇔ bi = b j = b, i 6= j = 1, . . . , m (2.2.31b)
m m
d √
∑ dx yi (0) = ∑ ai ci = 0 ⇒ (2.2.31c)
i=1 i=1
1 √
ui − b cosh( ci ì )

ai =
√ (2.2.31d)
sinh( ci ì )
√ !
m
√ m ui ci m
√ √
0 = ∑ ai ci = ∑ √ − ∑ ci coth( ci ì ) b ⇒ (2.2.31e)
i=1 i=1 sinh( ci ì ) i=1
√
1 m ui ci
b= m ∑ √ (2.2.31f)
√ √ sinh( ci ì )
∑ ci coth( ci ì ) i=1
i=1
 
√
ui  1 ui ci   coth(√ci ì )
ai = √ − m √ √ (2.2.31g)
sinh( ci ì )  √ sinh( ci ì ) 
∑ ci coth( ci ì )
i=1
The explicit solution is then assembled using ai , bi form (2.2.31) in (2.2.30). We have
ui
yi (x) = √ (2.2.32)
sinh( ci ì )
 
√
 1 ui ci   coth(√ci ì ) sinh(√ci x)
− √
 m √ √ sinh( ci ì ) 
∑ ci coth( ci ì )
i=1
√
1 m ui ci √
+ ∑ sinh(√ci ì ) cosh( ci x), i = 1, . . . , m.
m√ √
∑ ci coth( ci ì ) i=1
i=1
Example 2.2.3 A Helmholtz-type example is given by f i = 0, ci = −ω 2 ki < 0, where ω denotes a frequency.

In this case the solutions have the following form
p p
yi (x) = ai sin(ω ki x) + bi cos(ω ki x), i = 1, . . . , m. (2.2.33)
The same procedure asin the previous examples provides

m
√
1 ui ki
b= m √ √ ∑ sin(ω √ki ì ) (2.2.34)
∑ ω ki cot(ω k `
i i ) i=1
i=1
1
ai = √ ui
sin(ω ki ì )

uj
p
1 kj p
− √ ∑m √ cos(ω ki ì ) (2.2.35)

p
∑ mω k j cot(ω ki ` j ) j=1 sin(ω ki ` j )
j=1
Exercise 2.1. Let f i = 1, i = 1, . . . , m in (2.2.27). Compute the solutions yi (x) in the examples (2.2.28),(2.2.30)
and (2.2.33).
Exercise 2.2. Write a Matlab routine to solve the problem (2.2.27).
We now look at the same problem, however, with different nodal conditions at the multiple center node.
d2 i
− y + ci yi = f i in (0, ì )
dx2
yi (ì ) = ui
d i d (2.2.36)
y (0) = y j (0), i 6= j = 1, . . . , m
dx dx
m
∑ yi (0) = 0
i=1
Exercise 2.3. Repeat Excercise (2.1), (2.2) for (2.2.36).
Remark 2.2.4 If we introduce the vectors

 1  d 1
y dx y
 y2  d 2
d  dx y 

Y (0) =  .  (0), Y (0) =  .  (0)
 
 . . dx  .. 
y m d m
dx y
and the matrices

 
1 0 0 . . . −1  
 0 1 0 . . . −1 0 0 0 ... 0
  0 0 0 . . . 0
 0 0 1 . . . −1  
A0 =  0 0 0 . . . 0 ,
 , A1 =   (2.2.37)
. . . . . . . . . . . . −1

0 . . . . . . . . . 0
 0 . . . . . . 1 −1
1 1 1 1 1
0 0 0 0 0
then we can write down the nodal conditions for (2.2.27) as follows:
d
A0Y (0) + A1 Y (0) = 0.
dx
Accordingly, for (2.2.36), we have
d
A0 Y (0) + A1Y (0) = 0.
dx
We can, thus identify the nodal conditions as vector-values Robin-type conditions for the vectorial problem.
In particular, if one normalizes the lengths to 1, which can always be achieved by introducing a diagonal
matrix in the equations (see below), one can rewrite the star-graph problems as two point boundary value
problems for a vectorial elliptic equation or an elliptic system.
d2
AY − Y =F (2.2.38)
dx2
d
A0Y (0) + A1 Y (0) = 0, Y (1) = U, (2.2.39)
dx
where A can be regarded as yet another coupling matrix, this time the coupling occurs in the domains,
rather than on the boundary. In a sense, this system can be taken as one-dimensional prototype of a two-
point boundary value problem for elliptic systems. This remark can be extended to general Strum-Liouville
type problems on star-graphs and, in fact, on general metric graphs.
Remark 2.2.5 If the edges have individual lengths ì , we can use wi (x) = yi (xì ), x ∈ (0, 1). We define
diagonal matrices L := (ì aik )ik where aik is the adjacency matrix: aik = 1 if edge i is connected to edge
k (∃ j ∈ J S : di j 6= 0, dk j 6= 0) and 0 else. In the same way L 2 := (`2i aik )ik . Then the system (2.2.38) is
replaced with
d2
AY − L −2 Y =F
dx2
d
A0Y (0) + A1 L −1 Y (0) = 0, Y (1) = U.
dx
2.2 Planar graphs 11
Thus, a problem with varying lengths can be reinterpreted as a problem with normalized lengths but with
varying coefficients. In a more general context, this is precisely the basis for the method of the transforma-
tion on a fixed domain which is the basic tool in this book.
2.2.1 Steklov-Poincaré map for metric graphs
We consider again the general star-graph. We also assume that the data f i , i = 1, . . . , m are constant. In
a small neighbourhood of the central node this may be feasible. In this case the solution on each edge is
clearly quadratic. In fact, simple calculations show that the solutions yi look as follows:
1
yi (x) = − x2 f i + ai x + bi , i = 1, . . . , m (2.2.40)
2 
!
m
ui 1 i

1 i 1 1
+ ì f + `2i f i 

ai =  u − m ∑
ì ` 2 2
∑ `1i i=1 i
 
i=1
m i !
1 u 1 i
bi = b = m ∑ ` + 2 ì f , i = 1, . . . , m. (2.2.41)
∑ `1i i=1 i
i=1
We have then
1
yi (x) = − x2 f i
2
 
!
m
ui 1 i

1 1 1
ui − + `2i f i 

+ m ∑ + ì f x
ì  1 i=1 ì 2 2
`∑
i=1 i
m !
ui 1 i

1
+ m ∑ + ì f , i = 1, . . . , m.
1 i=1 ì 2
∑ ì
i=1
The corresponding Steklov-Poincaré map is then given by
d i
Si (u) = y (ì ) (2.2.42)
dx !
m
ui 1 i

1 1
= ui − m ∑ + ì f + `2i f i − ì f i (2.2.43)
ì 2 2
∑ `1i i=1
i=1
m m
1 ui 1 1 1
= ui − m
1 i=1
∑ ì − m ∑ 2 ì f i − 2 `2i f i − ì f i .
1 i=1
∑ ì ∑ ì
i=1 i=1
If we now assume that the connecting edges ei+m , i = 1, . . . , m start at the nodes ni , i = 1, . . . , m, then the
Steklov-Poincaré-based boundary conditions there read as follows:
m
d i+m 1 yi (ì )
y (0) − yi (ì ) + m ∑ (2.2.44)
dx
∑ 1 i=1 ì
ì
i=1
m
1 1 1
=− m
1 i=1
∑ 2 ì f i − 2 `2i f i − ì f i , i = 1, . . . , m. (2.2.45)
∑ ì
i=1
For the second transmission condition according to (2.2.36), we obtain:

1
yi (x) = − x2 f i + ax + bi , i = 1, . . . , m (2.2.46)
2
m !
1 i 1 2 i
a= m ∑ u + 2 ì f
∑ ì i=1
i=1
m !
1 1 2 i 1
i
bi = u − m ∑ u + 2 ì f ì + 3 `2i f i , i = 1, . . . , m.
i
(2.2.47)
∑ ì i=1
i=1
The corresponding Steklov-Poincaré map is then given by
d i
Si (u) = y (ì ) (2.2.48)
dx !
m m
1 1
= m ∑ u j + 2 ∑ `2j f j − ì f i .
∑ `j j=1 j=1
j=1
If we now assume that the connecting edges ei+m , i = 1, . . . , m start at the nodes ni , i = 1, . . . , m, then the
Steklov-Poincaré-based boundary conditions there read as follows:
!
m
d i+m 1 j
y (0) − m ∑ y (ì ) (2.2.49)
dx
∑ ` j j=1
j=1
!
m
1 1
= m ∑ `2j f j − ì f i , i = 1, . . . , m. (2.2.50)
2
∑ `j j=1
j=1
Similar calculations can be done for reverse Steklov-Poincaré problem, namely, where one starts with Neu-
mann data and solves for the Dirichlet data at the same nodes.
Exercise 2.4. Compute the Neumann-Dirchlet maps for the two star-graph problems (2.2.27), (2.2.36).
2.2.2 Problems on general metric graphs
The Steklov-Poincaré or Dirichlet-Neumann maps, together with the corresponding Neumann-Dirichlet

maps can be used to decompose any given problem on a metric graph into star-graphs. Let us briefly consider
this method. Let G = (N, E) denote the graph of the gas network with nodes N = {n1 , n2 , . . . , n|N| } an edges
E = {e1 , e2 , . . . , e|E| }. For the sake of uniqueness, we associate to each edge a direction.

 − 1, if node n j is the left node of the edge ei ,

di j = + 1, if node n j is the right node of the edge ei ,

0, else.

The edges incident at n j are labelled as follows: i ∈ I j := {i ∈ 1, . . . |E||di j 6= 0}. We introduce the edge
degree at a given node n j as d j := ∑ 1 = |{di j 6= 0|. The set N of nodes is labelled by J = { j = 1, . . . , |N|}
i∈I j
which is decomposed into those nodes with d j = 1, J S , the simple nodes, and d j ≥ 2, J M , the multi-
ple nodes. The simple nodes are split into Dirichlet and Neumann nodes with index sets JDS and JNS ,
respectively. The complete system of type (2.2.27) reads as follows:
d2 i
yi − γi y = f i , i ∈ I , x ∈ (0, 1)
dx2
yi (n j ) = yk (n j ), ∀i, k ∈ I j , j ∈ J M
d
∑ di j γi dx yi (n j ) = 0, j ∈ J M (2.2.51)
i∈I j
yi (n j ) = u j , i ∈ I j , j ∈ JDS
d i
y (n j ) = u j , i ∈ I j , j ∈ JNS ,
dx
whereas the complete elliptic network problem of type (2.2.36) is written as:
d2 i
yi − γi y = f i , i ∈ I , x ∈ (0, 1)
dx2
d d
γi yi (n j ) = γk yk (n j ), ∀i, k ∈ I j , j ∈ J M
dx dx
∑ di j yi (n j ) = 0, j ∈ J M (2.2.52)
i∈I j
yi (n j ) = u j , i ∈ I j , j ∈ JDS
d i
y (n j ) = u j , i ∈ I j , j ∈ JNS .
dx
The problem of well-posedness (2.2.52) is considered next. The analysis is similar for the even more clas-
sical system (2.2.51) which has been analyzed in the literature. See e.g. [?]. It turns out that, in this case,
for any graph G, such that each node is connected to a Dirichlet node by a simple path, the underlying
elliptic operator is self-adjoint and positive definite, thus, the problem admits a unique solution. Indeed, the
Lagrange identity
Z1
d2 i

hAy, ϕi := ∑ yi − γi y − f i
ϕ i dx (2.2.53)
i∈I 0 dx2
Z1
d i d id i
= ∑ ∑ di j γi y (n j )ϕ i (n j ) + ∑ yi ϕ i + γi y ϕ − f i ϕ i dx
j∈J i∈I j dx i∈I dx dx
0
d d
= ∑ ∑ di j γi yi (n j )ϕ i (n j ) − ∑ ∑ di j γi yi (n j ) ϕ i (n j ) + hq, Aϕi
j∈J i∈I j dx j∈J i∈I dx
j
= hq, Aϕi, ∀q, ϕ ∈ D(A),
where the domain D(A) is given by
D(A) := y = (yi )|i∈I ∈ Πi∈I H 2 (0, 1)|yi (n j ) = 0, i ∈ I j , j ∈ JDS

d i
y (n j ) = 0, i ∈ I j , j ∈ JNS (2.2.54)
dx
d d
γi q(n j ) = γk yk (n j ), ∀i, k ∈ I j , j ∈ J M
dx dx )
∑ di j yi (n j ) = 0, j ∈ J M
i∈I j
Thus, we can define the unbounded, self-adjoint operator
d2

Aq := yi − γi 2 yi (2.2.55)
dx i∈I
with domain D(A) in the Hilbert space H := Πi∈I L2 (0, 1) which is positive definite. We may also intro-
duce the energy space
V := y = (yi )|i∈I ∈ Πi∈I H 1 (0, 1)|yi (n j ) = 0, i ∈ I j , j ∈ JDS

d i
y (n j ) = 0, i ∈ I j , j ∈ JNS (2.2.56)
dx )
∑ di j yi (n j ) = 0, j ∈ J M
i∈I j
Theorem 2.2.6 Let a gas network G = (N, E) be given such that each node is connected to a Dirichlet
node by a simple path. Moreover, let controls u j , j ∈ J S and right hand sides f ∈ H be given. Then there
exists a unique solution y ∈ D(A) to Ay = f (2.2.52). If f ∈ V ∗ , then there is a unique solution y ∈ V .
Exercise 2.5. Provide a more detailed proof of theorem (2.2.6). Argue in a similar way for the system
(2.2.51).
We continue the discussion for the use of the Dirichlet-Neumann (Steklov-Poincaré) and Neumann-
Dirichlet maps for the purpose of decomposing problems on a general graph G into star-graph problems.
Assume we are given a multiple node n j , j ∈ J M with edge degree d j . Then d j edges emanate from n j .
Denote these edges i01 , . . . i0d j . If we cut these edges in the middle, then we obtain d j serial nodes n j1 , . . . , n jd j ,
where we apply the canonical transmission conditions and another set of d j edges with labels i11 , . . . i1d j . Then
we impose
0 1 i0 i1
yi1 (n j1 ) = yi1 (n j1 ), . . . , y d j (n jd j ) = y d j (n jd j ),
and
d i0 d 1 d i0 d i1
y 1 (n j1 ) = yi1 (n j1 ), . . . , y d j (n jd j ) = y d j (n jd j ).
dx dx dx dx
Then we use the Dirichlet-Neumann (D2N) or the Neumann-Dirichlet (N2D) mappings in order to cut out
the star-graph with center node n j and the edges with indices i01 , . . . i0d j as indicated above. Therefore, for a
first approach, it is sufficient to consider star-graphs as the contain the problem of multiple nodes. We do
not dwell further on problems on general graphs here and instead refer to [?, 210].
Exercise 2.6. Use the D2N and N2D mappings first for the decomposition of the serial problem
d2 1 d2
α1 y1 − 2
y = f 1 , x ∈ (0, ε), α2 y2 − δ 2 y2 = f 2 , x ∈ (ε, 1)
dx dx
d 1 d 2
y (0) = 0, y (ε) = y (ε), y (ε) = δ y (ε), y2 (1) = u,
1 1 2
dx dx
where δ > 0, 1 > ε > 0. That is to say, replace the coupling transmission conditions by inhomogeneous
Neumann conditions for, say, the first problem, using the D2N map.
Exercise 2.7. Use the last exercise in order to decompose two star-graphs connected by one, two or m − 1
edges.
2.2.3 Shape variations for graph problems
We are now interested in moving the center node of the graph 2.2 and to calculate the sensitivities with
respect that movement. We introduce the following notation. Let ai , i = 1, . . . , m and a be vectors in the
plane, where the ai , i = 1, . . . , m denote the endpoints of the edges labelled 1, . . . , m and a denotes the center
node in the reference configuration. We are going to move a in the direction v. To this end we introduce the
normalized edge vectors
a − ai a + tv − ai
ei := ei (0, 0) := , ei (t, v) := , i = 1, . . . , m (2.2.57)
ka − ai k ka + tv − ai k
See figure 2.2.3. The edges have lengths ì := ka − ai k, i = 1, . . . , m

There are several ways to introduce a vector field that transports the point a to the point a + tv, t > 0.
One natural choice is
Tt (v) : R2 → R2
(2.2.58)
Tt (v)(ai ) = ai , Tt (v)(a) = a + tv
a1
e1
a + tv
a
e3 e2
a3 a2
Fig. 2.4 The deformed three-star
We can realize this field by introducing a velocity V (ai ) = 0, i = 1, . . . , m and V (a) = v, hence,
kx − a1 kkx − a2 k . . . kx − am k
V (x) = v (2.2.59)
ka − a1 kka − a2 k . . . ka − am k
Note, however, that V deforms the graph in that the images of the edges Tt (v)(ai + se1 ) = ai + sei +tV (a1 +
sei ) are no longer straight lines. Alternatively, we can define m linear fields as follows.
sv
Tti (v)(ai + sei ) = ai + sei + t (2.2.60)
ì

which map the straight edges ai + sei , s ∈ [0, ì ] onto the straight edges ai + sei + t `vi , s ∈ [0, ì ]. Let
us introduce the lengths of the edges in the moving domain ì (t) := ka + tv − ai k, i = 1, . . . , m. We also
denote by ì (â) the individual length with a given center multiple node â. In case of â = a +tv − ai we have
ì (â) = ì (t).We can now pose the following optimization problem, where we look for the optimal position
in an ε neighbourhood of the given positon of the multiple node.
`Zi (a)
1 m
min ∑ |yi (x) − yd (x)|2 dx s.t.
â∈Bε(a) 2 i=1
0
d2 i
y + ci yi = f i in (0, ì (â))
dx2 (2.2.61)
yi (ì (â)) = ui
yi (0) = y j (0), i 6= j = 1, . . . , m
m
d
∑ dx yi (0) = 0
i=1
Example 2.2.7 We let ci = 0, i = 1, . . . , m and yd (x) = 0, x ∈ [0, ì ], i = 1, . . . , m. According to Fig. 2.2.3,

we consider the following transformations. See (2.2.58).
v
Tti (v)(x) := ai + xei + t x, i = 1, . . . , m (2.2.62)
ì
Then, clearly, Tti (v)(0) = ai , i = 1, . . . , m and Tti (v)(ì ) = a + tv. We can reformulate problem (2.2.61) as
follows.
d2 i
y = 0 in (0, ì (â))
dx2
yi (0) = ui
(2.2.63)
yi (ì (t)) = y j (` j (t)), i 6= j = 1, . . . , m
m
d
∑ dx yi (ì (t)) = 0
i=1
Let us observe that the unique solution is given by

 
1  i 1 uj  1 uj
yti (x) = u − ∑ m ` j (t)  (ì (t) − x) + ∑ m ` j (t) . (2.2.64)
ì (t) ∑ m ` j1(t) j=1 ∑ m ` j1(t) j=1
j=1 j=1
On the fixed graph we have
(yi )t (x) = yti (Tt (v)(x) (2.2.65)

 
1  i 1 uj  ì (t)
= u − ∑ m ` j (t)  (ì (t) − x)
ì (t) ∑ m ` j1(t) j=1 ì
j=1
1 uj
+ ∑ m ` j (t) .
∑ m ` j1(t) j=1
j=1
From this it is possible to derive the material and shape derivative, ẏi , (yi )0 , directly. We obtain for the
material derivatives
ẏi (0) = 0, i = 1, . . . , m, (2.2.66a)

m `0 (0) m u j `0j (0)
1 j uj 1
ẏi (ì ) = !2 ∑ 2 ∑ `j − ∑m (2.2.66b)
m j=1 ` j i=1 ∑ m `1j j=1 `2j
∑ `1 j=1
j=1 j
which is independent of i and

m
d `0i (0) d i
m
∑ dx ẏi (ì ) = − ∑ y (ì ). (2.2.66c)
i=1 i=1 ì dx
For the shape derivatives we get
d2 i 0
(y ) = 0, in (0, ì (â))
dx2
(yi )0 (0) = 0, i = 1, . . . , m, (2.2.67a)
d d
(yi )0 (ì ) − (y j )0 (` j ) = `0j (0) y j (` j ) − `0i (0) yi (ì ) (2.2.67b)
dx dx
m
d i 0
∑ (y ) (ì ) = 0.
i=1 dx
These relations can be verified using the system (2.2.63) directly, without using the explicit solution (2.2.64).
Exercise 2.8. Verify (2.2.66) and (2.2.67)!
Example 2.2.8 Let f i = ci = 0, i = 1, . . . , m and yd (x) = 0, x ∈ [0, `(â)] Then the silution on the moving
graph is given by  
m m
1  1 uj  1 uj
yti (x) = ui −
m x+
∑ ` j (t)  m ∑ ` j (t)
ì (t) 
∑ ` 1(t) j=1 ∑ ` 1(t) j=1
j=1 j j=1 j
We can equally well keep the center and move the the three edges. It turnes out, in fact, that for such scalar
a1
e1 a1 + tv1
a
e3 e2 a2 + tv2
a3
a3 + tv3 a2
Fig. 2.5 The three-star with transformed edges
problems, the solution does not depend on the angles between the edges. In scalar problems the variable
yi (x) is seen as either a temperature, in case we are dealing with a heat transport problem on the graph,
or it is viewed as an out-of-the-plane displacement, in case the problem is considered as representing a
network of strings with displacements pointing out of the plane. For problems on graphs in the plane, where
the variables yi (x) ∈ R2 are vectros in the plane, the situation is different. This case will be treated later. It
is, thus, sufficient to treat the case, where
ì (t)
Tti (v)(s) := s, s ∈ [0, ì ], (2.2.68)
ì
where,
ì (t) = kai − a + tvi k. (2.2.69)
Clearly,
(ai − a + tvi , vi ) 0 (ai − a, vi )
`0i (t) = , ì (0) = (2.2.70)
kai − a + tvi k kai − ak
We first transform the cost functional onto the fixed graph
`Zi (t) Zì
m m
ì (t) ì (t) 2
∑ |yti (x)|2 dx =∑ |yti ( s)| ds.
0 0 ì ì
0 0
We differentiate with respect to t and obtain

`Zi (t) Zì
d m d m ì (t) ì (t) 2
|yti (x)|2 dx = ∑ |yti ( s)| ds
dt ∑
0 dt 0 ì ì
0 0
Zì
m
`0i (t) ì (t) 2
=∑ |yti ( s)| ds+ (2.2.71)
i=1 ì ì
0
Zì
m
ì (t) ì (t)s ∂ i ì (t)s `0i (t)
∑ 2yti ( ) y( ) sds
0 ì ì ∂t t ì ì
0
We evaluate (2.2.71) at t = 0 and obtain

`Zi (t) Zì Zì
d m `0 (0)
m m
`0i (0)
|yti (x)|2 dx|t=0 =∑ i i 2
|y (s)| ds + 2 ∑ ẏi (s)yi (s) sds, (2.2.72)
dt ∑
0 i=1 ì 0 ì
0 0 0
˙ is the material derivative of yi (s). We use the relation between the shape derivative (yi )0
where now yi (s)
and the material derivative ẏi , namely,
d i `0i (0)
ẏi (x) = (yi )0 (x) + y (x) x.
dx ì
Then (2.2.72) takes the form
`Zi (t) Zì
d m m
`0i (0)
|yti (x)|2 dx|t=0 = ∑ |yi (s)|2 ds
dt ∑
0 i=1 ì
0 0
`
m Zi
d i `0i (0) i `0i (0)
+2 ∑ (yi )0 (x) + y (x) xy (s) sds
0 dx ì ì
0
(2.2.73)
m m Zì
= ∑ `0i (0)yi (ì )2 + 2 ∑ (yi )0 (x)yi (x)dx
i=1 i=1
0
m m Zì
= ∑ (ei , vi )yi (ì )2 + 2 ∑ (yi )0 (x)yi (x)dx
i=1 i=1
0
We consider the problem that the shape derivatives (yi )0 (·), i = 1, . . . , m have to solve. This problem can be
derived by direct differentiation.
d2 i 0
− (y ) = 0 in (0, ì )
dx2
d
(yi )0 (ì ) = − yi (ì )`0i (0)
dx (2.2.74)
(yi )0 (0) = (y j )0 (0), i 6= j = 1, . . . , m
m
d
∑ dx (yi )0 (0) = 0
i=1
Furthermore, we introduce the adjoint problem for the adjoint variable pi (·) as follows.
d2 i
− p = yi in (0, ì )
dx2
pi (ì ) = 0
(2.2.75)
pi (0) = p j (0), i 6= j = 1, . . . , m
m
d
∑ dx pi (0) = 0
i=1
We use the adjoint in (2.2.73) and integrate by parts and use the boundary conditions for pi in (2.2.75). We
have
m Zì m m
d2 i i 0 d i ì d i 0 i ì
−∑ (yi )0 p dx = − ∑ (y ) p |0 + ∑ (y ) p |0
i=1 dx2 i=1 dx i=1 dx
0
m
d i d
=∑ y (ì ) pi (ì )`0i (0)
i=1 dx dx
This gives the shape derivative of the cost functional as follows.
m
1 i 2 d i d
dJ(y, v) = ∑ y (ì ) + y (ì ) pi (ì ) (ei , vi ) (2.2.76)
i=1 2 dx dx
m
1 i 2 d i d i
∑ (u ) + dx y (ì ) dx p (ì ) (ei , vi ) (2.2.77)
i=1 2
We can for instance consider the following scenario: Let ui = u, ì = ` i = 1, . . . , m. It is clear from (2.2.68)
d i
evaluated at t = 0 that yi0 (ì ) = u 6= 0, i = 1, . . . , m and dx y (ì ) = 0 i = 1, . . . , m. Thus,
m
1
dJ(y, v) = (ui )2 ∑ (ei , vi ).
2 i=1
m
If we take the velocity vectors equal vi = v, i = 1, . . . , m, then ∑ ei = 0 is a condition which leads to
i=1
2π
stationarity. This is the configuration where all edges have the same angle of m.
2.2.4 Shape sensitivity for the first eigenvalue
We consider the following system on the star-graph
d2 i
− z = λ zi , i ∈ I , x ∈ (0, 1)
dx2
zi (0) = zk (0), ∀i, k = 1, . . . , m
m
d
∑ dx zi (0) = 0, (2.2.78)
i=1
zi (1) = 0, i = 1, . . . , m. (2.2.79)
We consider the smallest eigenvalue λ . To this end, we define
m Zì m Zì
d id i
a(z, ϕ) := ∑ y ϕ dx, (z, ϕ) := ∑ yi ϕ i dx (2.2.80)
i=1 dx dx i=1
0 0
m
V := {y ∈ Πi=1 H 1 (0, ì )|yi (0) = y j (0), i, j = 1, . . . , m, yi (ì ) = 0}
We pose the eigenvalue problem in variational form as follows.
a(z, ϕ) = λ (z, ϕ), ∀ϕ ∈ V. (2.2.81)
The smallest eigenvalue satisfies
λ (G) = min{a(ϕ, ϕ)|kϕk = 1}. (2.2.82)

ϕ∈V
Clearly, λ (G) > 0. We now consider Gt , where the edges are dependent on t: [0, ì (t)]. Consequently, we
look for the smallest eigenvalue λ (Gt ):
at (ϕ, ϕ) 2
λ (Gt ) = min {at (ϕ, ϕ)|kϕk = 1} = min { k }, (2.2.83)
ϕ∈Vt ϕ∈Vt \{0} kϕ
where
m
Vt := {y ∈ Πi=1 H 1 (0, ì (t))|yi (0) = y j (0), i, j = 1, . . . , m, yi (ì (t) = 0}. (2.2.84)
We introduce
at (ψ, ψ)
F(t, ψ) := p ,
k γ(t, x)ψk2
with γ i (t, x) := ∂x Tti (x). We obtain
∂ 1
F(0, ψ) = lim (F(t, ψ) − F(0, ψ))
∂t t→0+ t
 
 m Zì d d m Zi
` 
= lim ∑ yi ϕ i ∂x Tti (x)dx/ ∑ |ψ i |2 ∂x Tti (x)dx
t→0+ i=1 dx dx i=1

0 0
a0 (ψ, ψ)kψk2 − a(ψ, ψ)k γ 0 (0)ψk2
p
= p (2.2.85)
k γ(0)ψk4
m Zì
0
= a (ψ, ψ) − λ ∑ |ψ i |2 ∂xV i (0, x)dx, (2.2.86)
i=1
0
where we used that kψk = 1, γ(0) = 1. On the other side, by direct computation
d2 i 0
− (z ) = λ (zi )0 + λ 0 zi , i ∈ I , x ∈ (0, 1)
dx2
(zi )0 (0) = (zk )0 (0), ∀i, k = 1, . . . , m
m
d
∑ dx (zi )0 (0) = 0, (2.2.87)
i=1
d i i
(zi )0 (1) = − z V (0, ì ), i = 1, . . . , m. (2.2.88)
dx
We multiply by the eigenelement zi and obtain
m Zì m Zì
d2 i 0 i
−∑ (z ) z dx = λ ∑ (zi )0 zi dx + λ 0 . (2.2.89)
i=1 dx2 i=1
0 0
We come back to to the constrained optimization problem and introduce the Lagrangean:
m Zì m Zì
d2
L (G; z, β ) := ∑ | 2 (zi )|2 dx + β ( ∑ |zi |2 dx − 1) (2.2.90)
i=1 dx i=1
0 0
Shape variation in the direction of V yields at the optimum
m Zì
d i 0d i d d i 2 i
0 = dL (G; z, β ;V ) = ∑ 2 (z ) z + (( z ) V (0, x) dx
i=1 dx dx dx dx
0
m Zì
i 0 id i 2 i
+β ∑ 2(z ) z + ((z ) V (0, x) dx
i=1 dx
0
m Zì m m Zì
d d d
=∑ 2 (zi )0 zi dx + ∑ ( zi )2 (ì )V i (0, ì ) + 2β ∑ (zi )0 zi dx.
i=1 dx dx i=1 dx i=1
0 0
This implies
m Zì Zì
d i 0d i 1 m d m
∑ (z ) z dx + ∑ ( zi )2 (ì )V i (0, ì ) = −β ∑ (zi )0 zi dx. (2.2.91)
i=1 dx dx 2 i=1 dx i=1
0 0
Together with (2.2.89), (2.2.91) implies
m Zì Zì
1 m d i 2 m
λ∑ (zi )0 zi dx + λ 0 + ∑ ( z ) (ì )V i (0, ì ) = −β ∑ (zi )0 zi dx. (2.2.92)
i=1 2 i=1 dx i=1
0 0
We conclude β = −λ and
1 m d i 2 1 m d
dλ (G;V ) = − ∑ ( z ) (ì )V i (0, ì ) = − ∑ ( zi )2 (ì )`0i (0). (2.2.93)
2 i=1 dx 2 i=1 dx
We may now add a volume constraint

m Zì
∑ 1dx = M
i=1
0
together with a new Lagrange multiplier µ and rewrite the Lagrangean form as
Zì
 
m
L (G; z, µ) := λ (G) + µ  ∑ 1dx − M  . (2.2.94)
i=1
0
The shape derivative finally amounts to

m
dλ (G;V ) + µ ∑ V i (0, ì ), (2.2.95)
i=1
which gives
m
1 d
∑ (µ − 2 ( dx zi )2 (ì ))V i (0, ì ) = 0, ∀V i . (2.2.96)
i=1
d i
√
This clearly shows, that the boundary conditions for
√ dx z (ì ) are independent of i. But, as λ ì 6= kπ, we
d i √λ M
have dx z (ì ) = −a sin( λ elli ) . But this implies that the lengths are all equal: ì = m.
2.2.5 Transmission problem
We again consider the star-graph problem, however, for the sake of simplicity, with uniform lengths of the
edges, i.e. ì = `, i ∈ I . We introduce 0 < r < ` and a stiffness parameter δ > 0. We consider uniform
distributed forces f i = 1, i ∈ I , I = {i = 1, . . . , m}. Let χ(0,r)
δ be the function
(
δ δ x ∈ [0, r)
χ(0,r) (x) :=
1 x ∈ [r, `]
The transmission problem on the star-graph then reads as follows.
d2 i
δ
−χ(0,r) y = 1, x ∈ (0, ì ), i ∈ I
dx2
yi (r− ) = yi (r+ ), i ∈ I (2.2.97)
d d
δ yi (r− ) = yi (r+ ), i ∈ I
dx dx
d i
yi (0) = y j (0), ∑ y (0) = 0
i∈I dx
yi (`) = 0, i ∈ I .
We consider the optimization problem
Zr Z`
1 m 1 m
min J (Gδ , y) := ∑ |yi − yid1 |2 dx + ∑ |yi − yid2 |2 dx
2 i=1 2 i=1
0 r
such that (y, δ ) satisfies (2.2.97). (2.2.98)
We will specify the target functions yidk , k = 1, 2 in the example following the analysis. As of now, we
assume yid1 (r− ) = yid2 (r+ ), i = 1, . . . , m. We first consider the shape derivative (yi )0 (x) of yi at x.
d2 i 0
− (y ) = 0, x ∈ (0, ì ), i ∈ I
dx2
(yi )0 (r− ) = (yi )0 (r+ ), i ∈ I (2.2.99)
d d d
δ (yi )0 (r− ) − (yi )0 (r+ ) = (δ − 1) yi (r− )r0 (0), i ∈ I (2.2.100)
dx dx dx
d i
yi (0) = y j (0), ∑ y (0) = 0
i∈I dx
yi (`) = 0, i ∈ I .
In order to establish the gradient of the cost function, we note
m Zr Z`
dJ (r, y) = ∑ (y i
− yid1 )(yi )0 dx + (yi − yid2 )(yi )0 dx. (2.2.101)
i=1 r
0
In order to establish the shape gradient, we can, of course, either solve (2.2.99), or, which is more common
and most often more convenient, introduce the adjoint problem
d2 i
δ
−χ(0,r) p = yi − yid1 , x ∈ (0, r), i ∈ I
dx2
d2
− 2 pi = yi − yid2 , x ∈ (r, `), i ∈ I
dx
yi (r− ) = yi (r+ ), i ∈ I (2.2.102)
d d i +
δ pi (r− ) = p (r ), i ∈ I
dx dx
d i
pi (0) = p j (0), ∑ p (0) = 0
i∈I dx
pi (`) = 0, i ∈ I .
After some calculus, we arrive at

m
d i − d i −
dJ (r, y) = (1 − δ )r0 (0) ∑ p (r ) y (r ). (2.2.103)
i=1 dx dx
Then, depending of the sign of the sum, r0 (0) needs to be positive or negative in order to decrease the cost
function.
Example 2.2.9 We compute the solutions of (2.2.97). We obtain
(
1 1 1 2
i (( − 1)r2 − 2δ x x ∈ [0, r)
y (x) = 21 2δ 2
2 (` − x ) x ∈ [r, `].
We now specify the target functions. In order to provide simple evidence, we introduce 0 < r0 < r and take
(
1 1 1 2
i (( − 1)r2 − 2δ x x ∈ [0, r0 )
yd1 (x) = 12 2δ 2
2 (` − x ) x ∈ [r0 , r), yid2 = 12 (`2 − x2 ), x ∈ [r, `].
Clearly, in order to reduce the gradient to zero, r has to tend to to r0 and, therefore, r0 (0) < 0. In order to
derive this from (2.2.103), we calculate the adjoint explicitly. We obtain
pi (x) = (2.2.104)

 1 r2 r2 − 1 r4 − 1 r4 + (r2 r0 − 13 r03 − 23 r3 )` x ∈ [0, r0 )
δ − 1  2 2 0 14 30 24 3
− (r r0 − 3 r0 − 3 r )(x − `) + 23 r3 x − 12 r2 x2 + 12
1 4
x − 14 r4 x ∈ [r0 , r) (2.2.105)
2δ 2 
 2 1 3 2 3
(r r0 − 3 r0 − 3 r )(x − `) x ∈ [r, `]
This implies
(1 − δ )2 2 1
dJ (r, y) = mr( r3 + r03 − r2 r0 )r0 (0). (2.2.106)
2δ 2 3 3
As the factor in front of r0 (0) is positive, in order to have decent, we have to take r0 (0) < 0. Clearly, the
gradient is zero for y = r0 .
2.2.6 Topological derivative with respect to material properties
We now consider the sensitivity of states with respect to the inclusions. That is to say, we introduce a
different stiffness parameter close to the origin at the center node and establish a sensitivity result. This
time we consider inhomogeneous Dirichlet conditions at the simple nodes.
d2 i
− y = 0, x ∈ (0, 1), i ∈ I
dx2
yi (r− ) = yi (r+ ), i ∈ I (2.2.107)
d d
δ yi (r− ) = yi (r+ ), i ∈ I
dx dx
i j d i
y (0) = y (0), ∑ y (0) = 0
i∈I dx
yi (1) = gi , i ∈ I .
We compute the following solution:
yi (x) = (2.2.108)
m m
1 1 1
− (gi − ∑ g j )x + m ∑ g j , x ∈ [0, r)
δ (r − 1) − r m j=1 j=1
m m
1 1 1
− δ (gi − ∑ g j )x + m ∑ g j + gi
δ (r − 1) − r m j=1 j=1
m
δ 1
+ (gi − ∑ g j ), x ∈ [r, 1].
δ (r − 1) − r m j=1
With this, it is easy to calculate

!
1 i 1 1 − δ1 1 m
yr (x) − yi (x)0 = gi − ∑ gj x (2.2.109)
r δ 1 − r − δr m j=1
and, therefore, !
m
d i 1−δ 1
y (x)|r=0 − gi − ∑ gj x. (2.2.110)
dr r δ m j=1
2.2.7 Shape-topological derivatives
We now embark on situations where the solutions are not differentiable with respect to shape and topology
parameters in general. The generic example of such problems is provided by boundary obstacle problems.
2.2.7.1 Single link case
To this end, we first consider the single link case. For two variable parameters ε,t ∈ R+ , we start with the
description of geometry. We define a disconnected set joining two segments x ∈ (0, ε) ∪ (ε, r + t) such that
0 < r0 < r < r1 , r0 − r < t < r1 − r, and 0 < ε < ε0 < r0 . One geometric parameter ε associates the size
of inhomogeneity ωε = (0, ε) in the domain, and the other geometric parameter t defines position of the
moving boundary Γt = {x : x = r +t}. We thus look at a transmission problem with transmission boundary at
x = ε and consider unilateral inequality constraints at the end r +t, where t describes the moving boundary.
δ
The material stiffness represented with the help of the characteristic function such that χ(0,ε) (x) = δ for
x < ε, otherwise χ(0,ε)
δ (x) = 1 for x > ε, where δ ∈ R+ stands for a given stiffness parameter. Its two limit
cases correspond to the hole as δ & 0+ and to the rigid inclusion as δ % +∞.
For fixed g ∈ R, the space of functions is given by
Ht := {v ∈ H 1 (0, r + t) : v(0) = 0},
the admissible set is represented by the inequality constraint
Kt := {v ∈ Ht : v(r + t) + g ≥ 0},
and the variational inequality describing the unilateral problem takes the specific form: Find u(ε,t) − g ∈ Kt
such that
Z r+t
δ
χ(0,ε) (u(ε,t) )0 (v − u(ε,t) )0 dx ≥ 0 for all v − g ∈ Kt , (2.2.111)
0
where denoted the dependence of the solution on ε and t.

The variational inequality (2.2.111) implies the boundary value problem:
−(u(ε,t) )00 (x) = 0 for x ∈ (0, ε) ∪ (ε, r + t), (2.2.112)

u(ε,t) (0) = g, (2.2.113)
u(ε,t) (ε + ) − u(ε,t) (ε − ) = 0, (u(ε,t) )0 (ε + ) − δ · (u(ε,t) )0 (ε − ) = 0, (2.2.114)
u(ε,t) (r + t) ≥ 0, (u(ε,t) )0 (r + t) ≥ 0, (2.2.115)
(u(ε,t) )0 (r + t) · u(ε,t) (r + t) = 0, (2.2.116)
where ε − and ε + approach ε from below and above, respectively. It is derived from (2.2.111) in the standard
way by applying integration by parts for all v − g ∈ Kt :
Z r+t
(u(ε,t) )00 (v − u(ε,t) ) dx + (u(ε,t) )0 (r + t) v(r + t) − u(ε,t) (r + t)
δ

− χ(0,ε)
0
− (u(ε,t) )0 (ε + ) − δ (u(ε,t) )0 (ε − ) · v(ε) − u(ε,t) (ε) ≥ 0.

We establish the solution to (2.2.112)–(2.2.116) explicitly. Indeed, for arbitrary c(ε,t) ∈ R the relations
(2.2.112)–(2.2.114) can be solved by
c c
(
u(ε,t) (x) = g + (ε,t)
δ x, (u(ε,t) )0 = (ε,t)
δ , x ∈ (0, ε)
u(ε,t) (x) = g + c(ε,t) x + ε 1−δ (ε,t) )0 = c

δ , (u (ε,t) , x ∈ (ε, r + t)
implying the piecewise-linear continuous function
u(ε,t) (x) = g + c(ε,t) x + 1−δ

δ min{ε, x} .
(2.2.117)
With (2.2.117) the complementarity condition (2.2.116) takes the form
c(ε,t) · g + c(ε,t) (r + t + 1−δ

δ ε) = 0.
Hence, due to (2.2.115), the nonnegative constant c(ε,t) can be found uniquely
−1
c(ε,t) = max 0, −g r + t + 1−δ

δ ε . (2.2.118)
As ε & 0+ , from (2.2.117) and (2.2.118) we have the reference state
u(0,t) (x) = g + c(0,t) x, (2.2.119)

c(0,t) = max 0, −g(r + t)−1 ,

(2.2.120)
which solves the reference variational inequality corresponding to (??): Find u(0,t) − g ∈ Kt such that
Z r+t
(u(0,t) )0 (v − u(0,t) )0 dx ≥ 0 for all v − g ∈ Kt . (2.2.121)
0
In order to study the asymptotic properties of the solution, we need to consider a layer near the interface
point x = ε. To this end, we consider mapping (0, ε) 7→ (0, 1), x 7→ εy on the fixed domain. We pose the
auxiliary transmission problem: Find w ∈ H 1 (R+ ) such that
Z ∞
δ
χ(0,1) w0 (y)v0 (y) dy = (1 − δ )v(1) for all v ∈ H 1 (R+ ). (2.2.122)
0
Using integration by parts, variational equation (2.2.122) implies the boundary value problem:
−w00 (y) = 0 for y ∈ (0, 1) ∪ (1, ∞),

w(x) → 0 as x % ∞,
w(1+ ) − w(1− ) = 0, w0 (1+ ) − δ · w0 (1− ) = −(1 − δ ),
where 1− and 1+ approach 1 from below and above, respectively. The unique solution of this problem is
given by the piecewise linear continuous function
1−δ
w(y) = δ min{0, y − 1}. (2.2.123)
x
After stretching the coordinates y = ε in (2.2.123), we get the boundary layer
εw( εx ) = 1−δ
δ min{0, x − ε}, kεw( εx )k = O(ε 1/2 ) in H 1 (R+ ), (2.2.124)
where the square root asymptotic order is due to the semi-norm estimate
rZ
√
rZ
∞ 2 ε
1−δ 2
εw( εx )0 dx =

δ dx = O( ε).
0 0
We now provide an asymptotic expansion in a non-standard format:

Theorem 2.2.10 The solutions u(ε,t) and u(0,t) of variational inequalities (2.2.111) and (2.2.121) admit the
following residual estimate as ε & 0+ :
u(ε,t) = u(0,t) + ε q̃(ε,t) + O(ε 3/2 ) in H 1 (0, r + t) (2.2.125)
with the principal asymptotic term defined in Ht by
ε q̃(ε,t) (x) := (u(0,t) )0 (0) · εw( εx ) + ε 1−δ x

δ (1 − r+t ) .
(2.2.126)
proof:
Indeed, for sufficiently small ε we have (r + t)(r + t + 1−δ −1 > 0, hence from (2.2.118) it follows that
δ ε)
−1
c(ε,t) = 1 + δ1−δ · max 0, −g(r + t)−1 ,

(r+t) ε
and together with (2.2.120) this results in the expansion

2
c(ε,t) = c(0,t) 1 − δ1−δ

(r+t) ε + O(ε) . (2.2.127)
Substituting (2.2.127), (2.2.124), and (2.2.119) in (2.2.117) we get
u(ε,t) (x) = u(0,t) (x) − c(0,t) x + c(0,t) 1 − δ1−δ 1−δ x 2

(r+t) ε x + δ ε + εw( ε ) + O(ε )
and derive iteratively the following uniform estimates:

u(ε,t) (x) = u(0,t) (x) + O(ε 1/2 ),

u(ε,t) (x) = u(0,t) (x) + (u(0,t) )0 (0) · εw( εx ) + O(ε),
u(ε,t) (x) = u(0,t) (x) + (u(0,t) )0 (0) · εw( εx ) + ε 1−δ x 3/2 ),

δ (1 − r+t ) + O(ε
where we have used c(0,t) = (u(0,t) )0 (0). The latter equality enforces (2.2.125) with notation (2.2.126), thus
completing the proof.
Remark 2.2.11 We remark that the expansion (2.2.125) is non-standard in the sense that the ’first-order
term’ (2.2.126) still contains the small parameter, but in a uniform way. We will see below that this notation
is useful.
2.2.7.2 Uniform asymptotic expansion in the problem
Provided with the asymptotic formula (2.2.125), we embark on constrained optimization problems, where
the constraint is given by the variational inequality (2.2.111). We discuss examples for various objectives
J(u(ε,t) ) subject to the optimal state u(ε,t) . The state-constrained optimization problem takes the specific
form:
min J(u(ε,t) ) (2.2.128)

(ε,t)∈(0,ε0 )×(r0 −r,r1 −r)
subject to Π (u(ε,t) ) = min Π (v),

v−g∈Kt
and the strain energy (SE) functional Π : Ht 7→ R is

Z r+t
Π (v) := 1
2
δ
χ(0,ε) (v0 (x))2 dx. (2.2.129)
0
The variational inequality (2.2.111) implies the first order optimality condition for the constrained mini-
mization of Π over v − g ∈ Kt .
Remark 2.2.12 Problem (2.2.128) is an example of an MPEQ, i.e. a minimization problem with equilibrium
constraints. These problems in higher dimensions are subject to current research, in particular in the context
of optimization problems for contact and crack problems in mechanics. See e.g. [202, 204, 205, 206, 207,
208, 211, 212, 214, 215, ?]. We also remark that, for fixed ε ∈ (0, ε0 ), variations of the parameter t ∈
(r0 − r, r1 − r) describe regular perturbations of the moving boundary of the domain (0, ε) ∪ (ε, r + t), thus,
shape variations. In contrast, the limiting procedure ε & 0+ implies diminishing of the inhomogeneity
ωε = (0, ε), and, hence, the topology change from the disconnected set to the 1-connected set (0, r + t).
We are going to consider two examples, first we take the strain energy itself as cost function and in the
second example we take its shape derivative.
Example 2.2.13 We start with the control the optimal value function JSE = Π of the strain energy (2.2.129)
with respect to the topology change as ε & 0+ . Relying on small ε, we substitute the optimal state u(ε,t)
with its asymptotic model (2.2.125) and (2.2.126), thus calculating the approximation of the optimal value
function
Π (u(ε,t) ) = Π u(0,t) + c(0,t) εw( εx ) + ε 1−δ x 3/2

δ (1 − r+t ) + O(ε )
Z r+t
) 2
= 12 (u(0,t) )0 + c(0,t) εw0 ( εx ) − ε(1−δ dx + O(ε 3/2 )
δ

χ(0,ε) δ (r+t)
0
Z r+t
c2(0,t) h ε
Z i
) 1−δ 2 )
δ 1 − 2ε(1−δ 1 − 2ε(1−δ

= 2 δ (r+t) + δ 2 dx + δ (r+t) dx + o(ε) (2.2.130)
0 ε
c2(0,t) ε(1−δ )
= 2 (r + t − δ ) + o(ε),
c2
Π (u(0,t) ) = (0,t)
2 (r + t),
due to (2.2.119), (2.2.124), and (2.2.129). From (2.2.130) it follows that the function (0, ε0 ) 7→ R, ε 7→
Π (u(ε,t) ) is differentiable at ε = 0 with the topological derivative
d
dε Π (u
(ε,t)
)|ε=0 = −c2(0,t) (1−δ )
2δ = −Π (u
(0,t) 1−δ
) δ (r+t) . (2.2.131)
Example 2.2.14 In the second example, we control the objective function JSERR = − dtd Π of the strain
energy release rate, which implies shape variation and associates a Griffith functional used in fracture
mechanics.
In order to calculate − dtd Π , we introduce let a cut-off function η be such that η(x) = 0 as x < ε and
η(x) = 1 as x > ε + β , with some β such that ε + β < r0 . For small s ∈ (r0 − r −t, r1 − r −t), the translation
Φs : (0, r + t) 7→ (0, r + t + s), z = x + sη(x) yields the representation of Π (u(ε,t+s) ) as
Z r+t+s Z r+t
(u(ε,t+s) ◦Φs )0x 2
1
2
δ
χ(0,ε) (u(ε,t+s) )0z )2 dz = 1
2
δ
χ(0,ε) 1+sη 0 (1 + sη 0 ) dx
0 0
Z r+t 2
(ε,t+s)
= Π (u ◦ Φs ) − s
2
δ
χ(0,ε) (u(ε,t+s) ◦ Φs )0 η 0 dx + o(s).
0
Since u(ε,t+s) ◦ Φs − g ∈ Kt , we infer u(ε,t+s) ◦ Φs → u(ε,t) strongly in Ht as s → 0, and conclude, see [202]
for detail, with the asymptotic expansion
Z r+t 2
Π (u(ε,t+s) ) = Π (u(ε,t) ) − 2s δ
χ(0,ε) (u(ε,t) )0 η 0 dx + o(s). (2.2.132)
0
From (2.2.132) it follows directly the explicit formula of the shape derivative
Z r+t 2
JSERR (u(ε,t) ) := − dtd Π (u(ε,t) ) = 1
2
δ
χ(0,ε) (u(ε,t) )0 η 0 dx. (2.2.133)
0
We observe that JSERR depends on u(ε,t) , but not on ε q̃(ε,t) in expansion (2.2.125). The latter fact is in
accordance with the assertion in [211, 212].
For the shape-topological control, now we insert (2.2.125) in (2.2.133), which implies the asymptotic
model
JSERR u(0,t) + c(0,t) εw( εx ) + ε 1−δ x 3/2

δ (1 − r+t ) + O(ε )
Z r+t
) 2 0
= 21 (u(0,t) )0 + c(0,t) εw0 ( εx ) − ε(1−δ η dx + O(ε 3/2 )
δ

χ(0,ε) δ (r+t)
0
c2(0,t) c2(0,t)
Z ε+β
) 0 )
= 2 1 − 2ε(1−δ
δ (r+t) η dx + o(ε) = 2 1 − 2ε(1−δ
δ (r+t) + o(ε)
(2.2.134)
ε
= JSERR (u(0,t) ) − ε c2(0,t) δ1−δ
(r+t) + o(ε),
c2(0,t)
JSERR (u(0,t) ) = 2 .
In particular, (2.2.134) follows formula for the shape-topological derivative

(ε,t)
d
dε JSERR (u )|ε=0 = −c2(0,t) δ1−δ
(r+t) . (2.2.135)
2
Moreover, in view of definition (2.2.133), it implies the mixed second derivative − ∂ ∂ε ∂t Π (u(ε,t) )|ε=0 which
∂2 (ε,t) )| ∂2 (ε,t) )|
is symmetric: ∂ ε ∂t Π (u ε=0 = ∂t ∂ ε Π (u ε=0 . Thus, we have proved the following.
Theorem 2.2.15 For the solutions u(ε,t) and u(0,t) of variational inequalities (2.2.111) and (2.2.121), there
exists the shape-topological derivative
(ε,t) 2 2
d
dε JSERR (u )|ε=0 = − ∂ ∂ε ∂t Π (u(ε,t) )|ε=0 = − ∂t∂∂ ε Π (u(ε,t) )|ε=0
2 1−δ
(2.2.136)
= −c(0,t) δ (r+t) .
2.2.7.3 Serial case
We now look at a serial case as a preparation for the star-graph. To this end, we take two strings. The first
string extends from 0 to ` > 0, the second from 0 to the moving end at x = r + t, r > 0,t > 0. We consider a
subinterval [0, ε) on which the stiffness parameter δ ≥ 0 applies. In particular,
d2 1
− u = 0, x ∈ (0, `),
dx2 ε,t
d2
− 2 u2ε,t = 0, x ∈ (0, r + t)
dx
u1ε,t (`) = g, u2ε,t (r + t) ≥ 0 (2.2.137)
u1ε,t (ε + ) − u1ε,t (ε − ) = 0,
u2ε,t (ε + ) − u2ε,t (ε − ) = 0
d 1 + d
uε,t (e ) − δ u1ε,t (ε − ) = 0,
dx dx
d 2 + d
uε,t (e ) − δ u2ε,t (ε − ) = 0.
dx dx
We find the following explicit solutions for ε > 0,t > 0
(
1 g g≥0
uε,t (x) = δg (2.2.138)
δ (r+t+`)+2(1−δ )ε (x − `) + g g ≤ 0
(
g g≥0
u2ε,t (x) = δg δ (r+t) ,
− δ (r+t+`)+2(1−δ )ε x + δ (r+t+`)+2(1−δ )ε g≤0
whereas for ε = 0, we find

(
g g≥0
u10,t (x) = g (2.2.139)
r+t+` (x − `) + g g≤0
(
g g≥0
u20,t (x) = g r+t .
− r+t+` x + r+t+` g≤0
In contrast to the previous example, where the jump in the stiffness approaches the fixed boundary with
inhomogeneous Dirichlet conditions, no boundary layer appears. We find the asymptotic expansion in case
g < 0 (the other case is trivial):
2(1 − δ )
u1ε,t (x) = u10,t (x) − g ε(x − `) + O(ε 2 ) (2.2.140)
δ (r + t + l)2
2(1 − δ )
u2ε,t (x) = u20,t (x) + g ε(x − (r + t)) + O(ε 2 ). (2.2.141)
δ (r + t + l)2
Thus, in this case, we find shape-topological derivatives of the solutions directly.
d 2(1 − δ ) g
|ε=0 u1ε,t (x) = (ut1 )0 (x) = − g(x − `) (2.2.142)
dε δ (r + t + `)2
d 2(1 − δ ) g
|ε=0 u2ε,t (x) = (ut2 )0 (x) = g(x − (r + t)).
dε δ (r + t + `)2
Exercise 2.9. For fixed ε0 > 0, define ε = ε0 + s, s > 0. Compute the shape derivative of the solutions of
system (2.2.137) at the transmission boundary at x = ε0 . What happens if then ε0 → 0? Find a relation
between the shape-topological derivative (2.2.142) and shape derivative in the limit as ε0 → 0.
Exercise 2.10. Extend the result for the serial case to the case of a star-graph, where the stiffness parameter
δ > 0 is distributed in a circular neighbourhood of the center node at x = 0.
Exercise 2.11. Write a Matlab routine for the case of a single link and the serial link.
2.3 Stars with a hole
We consider a star-graph GJ 0 with m edges and center at the node vJ 0 . In particular, we may without loss
of generality, assume that the edges ei stretch from the center to the simple boundary nodes, which we will
2.3 Stars with a hole 29
label from 1 to m. By this assumption we consider the multiple node at the center as being reached at x = 0
for all outgoing edges. Thus, the data ui are picked up at the ends x = ì . As a slight conceptual variation
with respect to previous discussions on differential equations on networks, we now consider vectorial states
on the graph and, therefore, collect the stiffness information in a symmetric, positive definite constant
matrix Ki . By this, the reader should be able to redo the problems discussed in the previous sections in
the context of vectorial systems that are more relevant in the applications. Indeed, if one considers elastic
strings stretched from node to node in space, the resulting system becomes more relevant and gives rise
to still open questions, in particular for nonlinear systems. See []...[] for further reading. The systems then
reads as follows.
2

 −Ki d yi + ci yi = f i , i ∈ I

dx2




 i
 y (ì ) = ui , i = 1, . . . , m

(2.3.1)

 yi (0) = y j (0), ∀i, j = 1, . . . , m

 m
 d
∑ dx yi (0) = 0.




i=1
Fig. 2.6 Cutting a hole into star-like subgraph
We are going to cut out the center and connect the corresponding cut-nodes via a circuit as seen in the
Figure 2.6. In general we have numbers ρi ∈ [0, ì ), i = 1, . . . , m which are taken to be the lengths of the
edges that are cut out. Thus the remaining edges have lengths ì −ρi . At x = ρi we create a new multiple node
vi . We connect these nodes by edges em+i , i = 1, . . . , m with lengths σ i (ρi ). After that, these nodes receive
a new edge degree. In this paper we assume that all these nodes have the same edge degree di = 3. More
complicated cutting procedures can be introduced, but obscure the ideas of this first paper on topological
derivatives of graph problems.
The problem we have to solve is the following
d2 i

−Ki y + ci yi = f i , i ∈ I


dx2




yi (ì ) = ui , ı = 1, . . . , m






yi (ρi ) = ym+i (0) = ym+i−1 (σ i (ρi )), ∀i = 2, . . . , m



(2.3.2)

 y1 (ρ1 ) = ym+1 (0) = y2m (σ 2m (ρ2m )),

 d d d
−Ki yi (ρi ) − Km+i ym+i (0) + Km+i−1 ym+i−1 (σ m+i−1 (ρm+i−1 )) = 0, i = 2, . . . , m



dx dx dx





 d d d
−K1 y1 (ρ1 ) − Km+1 ym+1 (0) + K2m y2m (σ 2m (ρ2m )) = 0.


dx dx dx
We proceed to derive the solutions to (2.3.1) and(2.3.2), respectively. To this end we look at
d2 i d2
−Ki 2
y + ci yi = f i ⇔ 2 yi = ci Ki−1 (yi ) − Ki−1 f i
dx dx
and define Ai := ci Ki−1 , f i := −ci Ai f i . The general solution of the homogeneous equation ( f i = 0) is given
by
1 1
(yi )H (x) = sinh(Ai2 x)ai + cosh(Ai2 x)bi , (2.3.3)
where the sin- and cos-operators are defined vie spectral resolution
1 p 1
sinh(Ai2 (x))ξ = ∑ sinh(λi2j x)(ξ , ϕi j )ϕi j (2.3.4)
j=1
1
accordingly for cosh(Ai2 (x))
The inhomogeneous equation is then solved by variation of constants as follows
Zx
i I −1 1
(y ) (x) = Ai 2 sinh(Ai2 (x − s) f i (s)ds. (2.3.5)
0
We will treat the case f i = 0 only. The general case is then a matter of additional but straightforward
calculus.
Lemma 2.3.1 The solution y to problem (2.3.1) with f i = 0, i = 1, . . . , m is given by
1 1
(yi )(x) = sinh((Ai2 )(x))ai + cosh(Ai2 (x))b (2.3.6)
with the coefficient-vectors ai , b
1 1
ai = sinh(Ai2 ì )−1 (ui − cosh(Ai2 ì )
m m
− 12 1 − 12 1
·( ∑ ci Ai coth Ai2 ì )−1 ∑ ci Ai sinh(Ai2 ì )−1 ui (2.3.7)
i=1 i=1
m m
− 21 1 − 12 1
b = ( ∑ ci Ai coth(Ai2 ì ))−1 ∑ ci Ai sinh(Ai2 ì )−1 ui . (2.3.8)
i=1 i=1
The Stekov-Poincaré map is provided by

1 1 1
SJi 0 (u) = Ai2 (cosh(Ai2 ì )ai + sinh(Ai2 ì )b) (2.3.9)
with ai , b according to (2.3.7),(2.3.8).

The situation appears to be much more simple in case all material parameters and geometrical data are
equal.
1
ci = 1, Ki = Id = A 2 , ì = `, f i = 0, i = 1, . . . m (2.3.10)
Example 2.3.2 Let assumption (2.3.10) hold true. Then the solution y to (2.3.1) is given by
m
1 1
(yi )(x) = sinh(x)(ui − ∑ u j) (2.3.11)
sinh(`) m j=1
1 1 m
+ cosh(x) ∑ ui
cosh(`) m i=1
The Steklov-Poincaré map is given by

m m
1 1
Si (u)J 0 = coth(`)(ui − ∑ u j ) + tanh(`) m ∑ u j (2.3.12)
m j=1 j=1
We proceed to problem (2.3.2). Again, we will treat the general case first and will then restrict to as-
sumption (2.3.10) in order to better reveal the underlying structure.
We introduce the ansatz for the solution as follows
1 1
(yi )ρ (x) := sinh(Ai2 x)ai + cosh(Ai2 x)bi
ρ ρ
(2.3.13)
From the Dirichlet conditions in (2.3.2)2 we infer

1 1
(yi )ρ (ì ) = sinh(Ai2 (ì ))ai + cosh(Ai2 (ì ))bi = ui , i = 1, . . . , m.
ρ ρ
(2.3.14)
From the continuity requirement in (2.3.2)3,4 we obtain

1 1
(yi )ρ (ρi ) = sinh(Ai2 ρi )ai + cosh(Ai2 ρi )bi = (ym+i )ρ (0) = bm+i
ρ ρ ρ
(2.3.15)
m+i−1 m+i−1
= (y )ρ (σ (ρm+i−1 )), i = 2, . . . m
1 1
(y1 )ρ (ρ1 ) = sinh(A12 ρ1 )a1 + cosh(A12 ρ1 )b1
ρ ρ
(2.3.16)
m+1 ρ 2m 2m
= (y )ρ (0) = bm+1 = (y )ρ (σ (ρ2m ))
The Kirchhoff conditions in (2.3.2) result in

−1 1
ρ
1
ρ −1 ρ
−ci Ai 2 [cosh(Ai2 ρi )ai + sinh(Ai2 ρi )bi ] − cm+i Am+i
2
am+i (2.3.17)
− 21 1
(σ m+i−1 (ρm+i−1 )))am+i−1
ρ
+cm+i−1 Am+i−1 2
[cosh(Am+i−1
1
(σ m+i−1 (ρm+i−1 )))bm+i−1 ] = 0, i = 2, . . . , m
ρ
2
+ sinh(Am+i−1
−1 1
ρ
1
ρ −1 ρ
−c1 A1 2 [cosh(A12 ρ1 )a1 + sinh(A12 ρ1 )b1 ] − cm+1 Am+1
2
am+1 (2.3.18)
− 21 1
+c2m A2m [cosh(A2m (σ 2m (ρ2m )))a2m 2 ρ
1
(σ 2m (ρ2m )))b2m ] = 0.
ρ
2
+ sinh(A2m
ρ ρ
This set of equations (2.3.14)-(2.3.18) constitutes 4m conditions on the 4m unknowns ai , bi , i =
1, . . . , 2m. The problem is as to whether there is an asymptotic expansion of (yi )ρ in terms of ρ for small
ρ := (ρi )i=1,...,m . This problem is a singular perturbation problem. Notice that the graph with ρ = 0 is the
original star-graph with m edges, while for every ρ > 0 (i.e. ρi > 0), the graph has 2m edges and contains
exactly one circuit. We may of course also formally start with a star-graph consisting of 2m edges with
serial joints at xi = 0, xm+i = ρi , i = 1, . . . , m so that the edges ei , i = 1, . . . , m have length ì − ρi to begin
with, while the other edges em+i , i = 1, . . . , m stretch from the center (at xm+i = 0) to the serial nodes at
xm+i = ρi . But still, the perturbation is then singular with respect to the subgraphs spanned by the edges
em+i , i = 1 . . . m.
Our analysis depends on the expansion of the set of equations (2.3.14) to (2.3.18) up to second order
terms. The asymptotic analysis is based on the expansions of sinh(x), cosh(x) on the matrix level. We use
the asymptotic expansions in (2.3.4) as follows
( 1 1
sinh(Ai2 (σ i (ρi )))ξ = σ i (ρi )Ai2 ξ + O(ρi2 )
1 (2.3.19)
cosh(Ai2 (σ i (ρi )))ξ = ξ + O(ρi2 )
By (2.3.14) we have
1 1
ai = (sin(Ai2 (ì ))−1 (ui − cosh(Ai2 (ì ))bi ), i = 1, . . . , m
ρ ρ
(2.3.20)
We expand (2.3.15) and(2.3.16)
1
ρ ρ ρ
Ai2 ρi ai + bi = bm+i (2.3.21)
1
= σ m+i−1 (ρm+i−1 )Am+i−1 am+i−1 + bm+i−1 + O(ρ 2 ), i = 2, . . . m
2 ρ ρ
1 1
A12 ρ1 a1 + b1 = bm+1 = σ 2m (ρ2m )A2m a2m + b2m + O(ρ 2 )
ρ ρ ρ
2 ρ ρ
(2.3.22)
We now proceed to the Kirchhoff conditions at the multiple nodes (2.3.17),(2.3.18)
−1 ρ
1
ρ −1 ρ
−ci Ai 2 [ai + ρi Ai2 bi ] − cm+i Am+i
2
am+i
− 12 ρ
1
ρ (2.3.23)
+cm+i−1 Am+i−1 [am+i−1 + σ m+i−1 (ρm+i−1 )Am+i−1
2
bm+i−1 ]
2
= 0 + O(ρ ), i = 2, . . . , m
and
−1 ρ
1
ρ −1 ρ
−c1 A1 2 [a1 + ρ1 A12 b1 ] − cm+1 Am+1
2
am+1
− 12 1 (2.3.24)
ρ ρ
+c2m A2m [a2m + σ 2m (ρ2m )A2 b2m ] = 0 + O(ρ 2 ) 2
We reformulate the system (2.3.21),(2.3.22),(2.3.23),(2.3.24) as follows

1 1 1 1
2 ρ ρ
Ai−1 2
ρi−1 − tanh(Ai−1 ì−1 ) ai−1 − Ai2 ρi − tanh(Ai2 ì ) ai
1
ρ
+σ m+i−1 (ρm+i−1 )Am+i−1
2
am+i−1
1 1
= cosh(Ai ì )−1 ui −cosh(Ai−1

2 2
ì−1 )−1 ui−1, i = 2, . . . , m (2.3.25)
1 1 1 1
ρ ρ
− A12 ρ1 − tanh(A12 `1 ) a1 + Am2 ρm − tanh(Am2 `m ) am
1 1 1
a2m = cosh(A12 `1 )−1 u1 − cosh(Am2 `m )−1 um + O(ρ 2 )
ρ
+σ 2m (ρ2m )A2m
2

−1 1
ρ
− ci Ai 2 + (σ m+i−1 (ρm+i−1 )cm+i−1 − ρi ci ) tanh(Ai2 ì ) ai
−1 ρ −1 ρ
2
−cm+i Am+i 2
am+i + cm+i−1 Am+i−1 am+i−1
1
= − σ m+i−1 (ρm+i−1 −1

)cm+i−1 − ρi ci cosh(Ai ì ) ui , i = 2, . . . m
2

1 1
(2.3.26)
− ρ
− c1 A1 2 + (σ 2m (ρ2m )c2m − ρ1 c1 ) tanh(A12 `1 )a1
−1 ρ −1 ρ
2
−cm+1 Am+1 am+1 + c2m A2m2 a2m
1
= − σ 2m (ρ2m )c2m − ρ1 c1 cosh(A12 `1 )−1 u1 + O(ρ 2 )

Now, (2.3.25)-(2.3.26) constitute a system of 2m linear asymptotic equations to order 2 in the 2m vari-
ρ
ables ai , i = 1, . . . , 2m.
ρ
Theorem 2.3.3 The system of equations (2.3.23) to (2.3.26) admits a unique solution ai , i = 1, . . . 2m.
Moreover, we have the asymptotic expansion
ρ
ai = ai + O(ρ), i = 1, . . . , m, (2.3.27)
where ai is given by (2.3.7) There exists a function si (·) such that
(yi )ρ (x) = (yi )(x) + O(ρ)si (x), i = 1, . . . , m, (2.3.28)
where (yi ) is the solution of the star-graph problem (2.3.1) ρ = 0.

Proof: Using equations (2.3.21) and (2.3.22), taking appropriate differences, we realize that bi = b̂ + O(ρ).
ρ
This information is inserted into equations (2.3.23) and (2.3.24). If we write all quantities involving ai with
indices i = 1 . . . m on the left and the other terms on the right side, we obtain after summing up, using a
’telescope-sum’, only O(ρ)-terms on the right hand side, i.e. we have
q
− 12 ρ
∑ ci Ai ai = O(ρ) (2.3.29)
i=1
ρ
Then we use the expression (2.3.20) for ai in (2.3.29) to obtain
!
m m
−1 1
−1 −1 1
−1
∑ ci Ai 2 sinh(Ai ì ) ui =2
∑ ci Ai 2 2
coth(Ai ì ) b̂
i=1 i=1
ρ
From this and (2.3.8) we see that up to terms of order O(ρ), b̂ = b. Then ai , up to the order O(ρ), are given
by ai in (2.3.7).
2.3.1 Homogeneous networks
In this subsection we consider the network under the assumption (2.3.10), i.e. all material and geometrical
quantities are the same, and a symmetric hole. Under this assumption the system of equations (2.3.25) to
(2.3.26) reduces to
ai−1 − ai − σ ρ coth(`)am+i−1 = − 1+ρ coth(`)

ρ ρ ρ 2
sinh(`) (ui − ui−1 ) + O(ρ ),
−a1 + am − σ ρ coth(`)a2m = − 1+ρ coth(`)
ρ ρ 2
sinh(`) (u1 − um ) + O(ρ ), (2.3.30)
ρ ρ ρ 1−σ
−(1 + (σ − 1)ρ tanh(`))ai − am+i + am+i−1 = cosh(`) ui + O(ρ 2 )
ρ ρ ρ 1−σ
−(1 + (σ − 1)ρ tanh(`))a1 − am+1 + a2m = cosh(`) u1 + O(ρ 2 ),
where the first and the third equations hold for i = 2, . . . , m, respectively. This system has a very particular
sparse structure which reflects the adjacency structure of the graph. To obtain the direct explicit solution is,
nevertheless, a matter of substantial calculations. Instead we look at an example.
Example 2.3.4 In this example we reduce the graph to a tripod. See figure 2.6. Here we can solve (2.3.30)
analytically an obtain
3
ρ 1 1
ai = sinh(`) (ui − 3 ∑ u j)
(j=1
3
1
+ρ cosh(`) (1 − 31 σ ) coth(`)2 (ui − 13 ∑ u j ) (2.3.31)
j=1
)
3
+ (σ − 1) 13 ∑ u j + O(ρ 2 ),
j=1
3
ρ 1 1
bi = cosh(`) 3 ∑ uj
j=1 (
3
sinh(`)
(1 − 31 σ ) coth(`)2 (ui − 31 ∑ u j )

−ρ cosh(`) 2 (2.3.32)
j=1

3
+ (σ − 1) 13 ∑ ui + O(ρ 2 ),
i=1
where i = 1, 2, 3.
ρ
We also display the coefficients am+i , i = 1, 2, 3 in order to reveal the behavior of the edges introduced
by cutting the hole.
ρ 1
a4 = (u2 − u1 ) (2.3.33)
3 sinh(`)
ρ σ
(1 − ) coth(`))(u2 − u1 ) + O(ρ 2 )

+
3 sinh(`) 3
ρ 1
a5 = (u3 − u2 ) (2.3.34)
3 sinh(`)
ρ σ
(1 − ) coth(`))(u3 − u2 ) + O(ρ 2 )

+
3 sinh(`) 3
ρ 1
a6 = (u1 − u3 ) (2.3.35)
3 sinh(`)
ρ σ
(1 − ) coth(`)(u1 − u3 ) + O(ρ 2 )

+
3 sinh(`) 3
The remaining bm+i , 1 = 1, 2, 3 are of course given by bi , i = 1, 2, 3 according to (2.3.21),(2.3.22). This

completely determines the solution (yi )ρ (x), i = 1, . . . , 6. We list the first three members for easier reference:
3 3
1
(yi )ρ (x) = ui − 13 ∑ u j sinh(x) + cosh(`)
1 1

sinh(`) 3 ∑ u j cosh(x)
( "j=1 j=1
3
1
+ρ cosh(`) (1 − 13 σ ) coth(`)2 (ui − 13 ∑ u j )
j=1
#
3
+ (σ − 1) 31 ∑ u j sinh(x) (2.3.36)
j=1
"
3
sinh(`) 1 2 1
− cosh(`) 2 (1 − 3 σ ) coth(`) (ui − 3 ∑ u j )
j=1
# )
3
+(σ − 1) 13 ∑ u j cosh(x) + O(ρ 2 ), i = 1, 2, 3
j=1
The Steklov- Poincaré-map is then obtained using

3 3
d
( dx (yi ))ρ (`) = coth(`)(ui − 31 ∑ u j ) + tanh(`) 13 ∑ u j
( j=1 j=1
3
+ρ (1 − tanh2 (`))[(1 − 13 σ ) coth2 (`)(ui − 31 ∑ u j )) (2.3.37)
j=1
)
3
+ (σ − 1) 13 ∑ u j ] , i = 1, . . . m.
j=1
It is apparent that (2.3.36),(2.3.37) provide the second order asymptotic expansion we were looking for.
We consider the following experiment: we apply longitudinal forces ui = uei with the same magnitude at
the simple nodes of the network. The (outer) edges ei , 1 = 1, 2, 3 or, respectively the edges of the original
star, are given by √ √
3 1 3 1
e1 = (0, 1), e2 = (− , − ), e3 = ( ,− )
2 2 2 2
which together with the orthogonal complements
√ √
1 3 1 3
e⊥
1 = (−1, 0), e⊥
2 = ( , − ), e⊥
3 = ( , )
2 2 2 2
3
form the local coordinate systems of the edges. Obviously ∑ ei = 0. Thus the solution to the unperturbed
i=1
problem is given by
1
(yi )(x) = u sinh(x)ei (2.3.38)
sinh(`)
This is in agreement with the fact that that particular reference configuration is completely symmetric. Now,
d
the solution (yi )ρ to the perturbed system and ( dx (yi ))ρ (`) are then given by
1
(yi )ρ (x) = sinh(x)uei
sinh(`)
σ 1
+ρ(1 − ) (coth(`) sinh(x) − cosh(x)) uei + O(ρ 2 ) (2.3.39)
3 sinh(`)
d σ
((yi ))ρ (`) = coth(`)uei + ρ(coth(`)2 − 1)(1 − )uei + O(ρ 2 )
dx 3
The energy of the unperturbed system is given by
Z`
1 3 d i d i 3
E0 = ∑ (y ) · (y ) + (yi ) · (yi )dx = coth(`)u2 (2.3.40)
2 i=1 dx dx 2
0
The energy of the perturbed system is given by

`−ρ Zσ ρ
1 3 d d 1 6 d i d i
Z
E ρ
= ∑ [ (yi ) · (yi ) + (yi ) · (yi )]dx + ∑ [ (y ) · (y ) + (yi ) · (yi )]dx (2.3.41)
2 i=1 dx dx 2 i=4 dx dx
0 0
1 1 1 σ
= hSρ u, ui = hS0 u, ui + ρ (1 − ) ((coth(`))2 − 1) u2

(2.3.42)
2 2√ 2 3
1 0 3 √
= hS u, ui + ρ ( 3 − 1) sinh(`)−2 u2 (2.3.43)
2 2
From these experiments we may draw the conclusion, that nodes of edge √ degree 3 under symmetric load,
where the configuration is at 120° between the edges (this amounts to σ = 3) are not going to be replaced
by hole, which would, in turn result in 3 new multiple nodes of edge degree 3. This seems to support the
optimality of such graphs being observed by Buttazzo [?].
Remark 2.3.5 The completely analogous formulae are obtained in the scalar case ((yi )(x) ∈ R), relevant
for instance in problems of heat transfer or electrical currents in networks.
If the loads are not symmetric, and/or if the geometry of the ’hole’ is not uniform, the energy may in
fact drop. A more detailed analysis can be found in []. Suffice it to say here, that nodes with higher edge
degree, according to our analysis, are ’more likely’ to be released by a hole, as even in the symmetric case
the number σ (ρ) which measures the new edge-lengths will be less than 1.
This is true e.g. for a node with edge degree 6 and beyond. Thus, the total length of the new edges is
smaller than the total length of the removed edges. This, in turn, is intuitive with respect to the fact that in
the higher-dimensional problem (in 2- or 3-d, no graphs), digging a hole reduces the amount of mass.
Example 2.3.6 Here we consider the homogeneous situation for a star with edge degree 6 at the multiple
node. In this case σ = 1 for the symmetric situation. See Figure 2.7
9
3
10
5
8
11
7
6 12
Fig. 2.7 Graph with ’critical’ edge degree 6
We calculate
6
ρ 1 1
a1 = (u1 − ∑ u j) (2.3.44)
sinh(`) 6 j=1
cosh(`)
+ρ {(−u5 − u3 − 4u2 − 4u6 + 10u1 )
cosh2 (`) − 1
)
1 6
−7(u1 − ∑ u j )
6 j=1
Notice that the edges 2 and 6 are the ’neighboring’ edges of edge 1 in the original star-graph. The other
ρ ρ
coefficients ai , 1 = 2, . . . , 6 are then obvious. For the sake of brevity, we only display e.g. a12 :
ρ 1
a12 = [5(u1 − u6 ) + 3(u2 − u5 ) + (u3 − u4 )]
12 sinh(`)
cosh(`)
−ρ [25(u1 − u6 ) − 9(u2 − u5 ) − 7(u3 − u4 )]
144(cosh2 (`) − 1)
+O(ρ 2 ) (2.3.45)
Again, observe that edge 12, in terms of the edges of the original graph, has direct neighbors 1 and 6,
the next level is 2 and 5 and finally we have 3 and 4. One realizes a consequent scaling. Also note that
ρ ρ
ai = 0 if ui are all equal. This shows that the coefficients bi in that case are independent of ρ and thus the
energy will not change for this limiting case.
2.4 The topological derivative
We are now in the position to define the topological derivative of an elliptic problem on a graph.
Let G be a graph, and let vJ ∈ JM be a multiple node with edge degree dJ . Let Gρ be the graph obtained
dJ
from G by replacing vJ with a cycle of length ∑ ci ρ with vertices v1J , . . . vdJ J of edge degree 3 each, such
i=1
that the distance from vJ to viJ is equal to ρ. Thus, the number nρ of edges of Gρ is n + dJ . Let J : G → R
be a functional on the edges of G
n Zì
d d i
J(y, y, G) := ∑ F(x, yi , y) (2.4.1)
dx i=1 dx
0
and let
ρ
n+dJ Zì
d d
J(yρ , yρ , Gρ ) := ∑ F(x, yρi , yρ ) (2.4.2)
dx i=1 dx i
0
be its extension to Gρ. Assume we have an asymptotic expansion as follows
d d
J(yρ , yρ , Gρ ) = J(r, y, G) + ρT (vJ ) + O(ρ 2 ) (2.4.3)
dx dx
then we define the topological gradient of J(Gρ ) with respect to ρ for ρ = 0 at the vertex vJ as follows.
d d
J(yρ , dx yρ , Gρ ) − J(r, dx y, G)
T (vJ ) = lim (2.4.4)
ρ→0 ρ
We consider the energy functional or, equivalently, the compliance which is the most natural criterion to
begin with. There are five such functionals relevant for the analysis of this paper: E 0 (y) on the entire graph
0
G , E ρ (yρ ) on the entire graph with the hole Gρ , ECS (y) on the graph G \ S J , where the star-graph without
0
hole S J has been cut out along edges ei , i ∈ IJ 0 , ES0 (y; v) on the star-graph without hole, and ES (y; v) on
ρ
the star-graph with hole. Obviously

1
ES0 (y; u) = hS0 u, ui, (2.4.5)
2
ρ 1
ES (y; u) = hSρ u, ui, (2.4.6)
2
E 0 (y) = ECS (y) + ES0 (y, y), E ρ (y0ρ ) = ECS (yρ ) + ES (yρ , yρ ),
ρ
(2.4.7)
ρ
where it is understood that in ES (yρ , ·) and ES0 (y, ·) we insert ui = yρ (ì ) and ui = y0 (ì ), respectively. Thus
1 1
E ρ (yρ ) − E 0 (y) = hSρ (ỹ), ỹi − hS0 (ỹ), ỹi, (2.4.8)
2 2
2.4 The topological derivative 37
0
where ỹ solves the problem on G \ S J and ui = ỹi (ì ), i ∈ IJ 0 . Thus the asymptotic analysis of the last
section carries over to the entire graph. As we have done the complete asymptotic analysis up to order 2 in
the homogeneous case only, we consequently dwell on this case now,the more general case will be subject
of a forthcoming publication.
2.4.1 Homogeneous graphs
In order to find an expression of the topological gradient in terms of the solutions y at the node vJ 0 , the one
that is cut out, we need to express the solution in terms of the data ui .
Example 2.4.1 We consider the star-graph as above with 3 edges. Obviously

3 3
1 d 1
ui −
3 ∑ u j = sinh(`) dx yi (0), 3 ∑ u j = cosh(`)yi (0). (2.4.9)
j=1 j=1
3 3 3 3
Thus using the fact that ∑ kui − 31 ∑ u j k2 = ∑ kui k2 − 31 (k ∑ u j k)2 we can express the bilinear expression
i=1 j=1 i=1 i=1
d
hS ρ (u), ui in terms of ky0 (0)k2 and k( dx y0 (0)k2 (where we omit the index 0) as follows
ρ
= hSi0 (u), ui
hSi (u), ui
3 3 (2.4.10)
1 d i
+ρ (1 − 3σ) ∑ k dx y (0)k2 + (σ − 1) ∑ kyi (0)k2
i=1 i=1
says that the energy function in the homogeneous case, when cutting out a symmetric hole e.g. σ i =
This √
σ = 3, i = 1, 2, 3, we have
( )
3 3
1 1 d i
TE (y, vJ 0 ) = 2 i
(1 − σ ) ∑ k y (0)k + (σ − 1) ∑ ky (0)k 2
(2.4.11)
2 3 i=1 dx i=1
The situation will be different for such vertices having a higher edge-degree as 6, and those having non-
symmetric holes. We expect that such networks are more likely to be reduced to edge-degree 3 by tearing a
hole. But this has to be confirmed by more detailed studies.
More general functionals will be considered in a forthcoming publication.
2.4.2 Sensitivity with respect to edge inclusion
We now consider a different situation where a node with edge degree dJ = N is released into a node of edge
degree 3 and one of degree N − 1 by introduction of a new edge eN+1 . See the figure 2.8 below.
Fig. 2.8 N-node turns into 3-node plus (N+1)-node
We consider this procedure in an explicit example with edge degree 4.

Let, therefore, vJ be a node with edge degree 4. As visualized in figure 2.8, we will introduce an addi-
ρ ρ ρ
tional new edge e5 of length ρ > 0 which together with the two new edges e1 , e2 given by
ρ `1 e1 − ρeN+1
e1 :=
k`1 e1 − ρeN+1 k
ρ `2 e1 − ρeN+1
e2 := , (2.4.12)
k`2 e2 − ρeN+1 k
where in our case study below N = 4.

ρ ρ
The new lengths ` − σ of the edges e1 , e2 (we consider a symmetric situation where the new additional
edge eN+1 equally divides the angle between e1 , e2 with an inclination α towards the corresponding unit
vectors) can b e computed by elementary trigonometry. The number σ is then found to be
1 1
σ = ρ cos α − ρ 2 (1 − cos2 α) + O(ρ 3 ) (2.4.13)
2` `
It is interesting to notice that for cos α > 12 the new graph has actually a smaller total length. This is
in contrast to the standard situation where cutting out a hole - which in fact implies introducing new the
edges forming that hole - has the opposite effect. For the sake of simplicity, we calculate the sensitivities
with respect to introducing the new edge of length ρ for the Laplacian on the graph only. Thus, we do not
consider an extra stiffening part due to the presence of a term cyi .
d2 i

− y = 0 in Ii , i = 1, . . . , 5


dx2




yi (`) = ui , i = 1, . . . 4,






y1 (σ ) = y2 (σ ) = y5 (ρ),



d 1 d d (2.4.14)

 y (σ ) + y2 (σ ) − y5 (ρ) = 0,


 dx dx dx
3 4
y5 (0),

y (0) = y (0) =





 d 3

d d
y (0) + y4 (0) + y5 (0) = 0


dx dx dx
We perform a similar analysis as in section 2.3 and therefore omit the details. We obtain
1 1 4 1 4
y1ρ (x) = (u1 − ∑ ui )x + ∑ ui
` 4 i=1 4 i=1
( )
ρ 1 1 4
− 2 [ cos α + 1](u2 − u1 ) + (2 − cos α)(u1 − ∑ ui ) (x − `)
2` 2 4 i=1
= y10 (x)
( )
ρ 1 1 4
− 2 [ cos α + 1](u2 − u1 ) + (2 − cos α)(u1 − ∑ ui ) (x − `)
2` 2 4 i=1
y2ρ (x) = y20 (x)

( )
ρ 1 1 4
− 2 ([ cos α + 1](u1 − u2 ) + (2 − cos α(u2 − ∑ ui ) (x − `)
2` 2 4 i=1
y3ρ (x) = r30 (x)

( )
ρ 1 1 4
− 2 [1 − cos α](u2 − u1 ) + (2 − cos α(u2 − ∑ ui ) (x − `)
2` 2 4 i=1
2.4 The topological derivative 39
y4ρ (x) = r40 (x)

( )
ρ 1 1 4
− 2 [1 − cos α](u1 − u2 ) + (2 − cos α)(u1 − ∑ ui ) (x − `)
2` 2 4 i=1
In order to calculate the energy we use the Steklov-Poincaré mapping and multiply by yi (`).
As before, the calculations can be done for scalar problems as well as for vectorial in-plane models. We
dispense with the display of the lengthy formulae. Instead, we give two different scenarios for topological
derivatives.
Example 2.4.2 In the scalar case we may set u1 = u2 and u3 = u4 = 0, i.e., we apply Dirichlet conditions
at the ends of edges 3 and 4 and pull at the end of the edges 1 and 2 by the same amount. This results in:
ρ
hSρ u, ui = hS0 u, ui − (2 − cos α)u2 (2.4.15)
2`2
Obviously, the introduction of a new edge is enhanced. One obtains a decomposition into two multiple nodes
with edge degree 3
Example 2.4.3 In the second example we take the planar model and set u1 = ue1 , u2 = ue2 and again
u3 = 0 = u4 . Now we obtain

3ρ 2 4
hSρ u, ui = hS0 u, ui − 2 cos(α)(cos2 α + cos α − ) u2 (2.4.16)
4` 3 3
For small enough angles α (e.g. 0 < α < π/6) the expression with ρ, i.e. the topological derivative of the
energy becomes negative. This shows that in the planar situation, the opportunity to create an additional
edge depends on the angles between the edges 1 and 2.
Obviously, the examples above can be generalized to more general networks including distributed loads
and obstacles. It is also possible to extend this analysis to 3-d networks.
Exercise 2.12.
Exercise 2.13.
Exercise 2.14.
Chapter 3
Mixed formulation of elliptic boundary value problems
3.1 Introduction
We are going to use the saddle-point variational formulation for elliptic problems. The method is well
known in numerical solution of PDEs.
3.2 Mixed formulation in one space dimension
Gatica - Mixed Finite Element Method: str. 10-14

Denote Ω = (0, 1), and consider the Neumann problem
−u00 + κu = f in Ω , u0 (0) = a, u0 (1) = b.
R1
The weak formulation in H 1 (Ω involves the bilinear form a(u, v) = {u0 v0 + κuv} and linear form F(v) =
0
R1
f v + {bv(1) − av(0)}.
0
Exercice 3.2.1 str. 11 apply the Lax-Milgram lemma to the variational formulation
The mixed formulation is a first order system, and it is obtained with σ := u0 and
−σ 0 + κu = f in Ω , σ (0) = a, σ (1) = b.
Exercice 3.2.2 str. 12 apply the generalized Lax-Milgram lemma to the mixed formulation by using (str 13
the translation technique or the Lagrange multiplier technique) to the integral formulation
3.3 Generalized Lax-Milgram Lemma
Gatica - Mixed Finite Element Method: str. 9
3.4 Mixed formulation of elliptic boundary value problems
Let us consider the mixed formulation of the principal elliptic boundary problems in two and three spatial
dimensions. At the first stage, we are going to write the mixed formulations for the state equations. Then,
the adapted mixed formulation is given for the associated shape derivatives of solutions. Thus, we do not
need to differentiate, with respect to the shape, of the mixed variational problems, this differentation is
performed at the level of variational formulations for the state equations.
41
42 3 Mixed formulation of elliptic boundary value problems
Notation and examples
Given f ∈ L2 (Ω ), the Poisson problem for u ∈ H01 (Ω ) reads
−∆ u = f ,
can be formulated as follows, seek a tuple (p, u) ∈ L2 (Ω ; R2 ) × H01 (Ω ), such that
∇u = p , (3.4.1)
−div p = f . (3.4.2)
It means that instead of one scalar function u we need to compute two functions u, p, where p stands for the
gradient of u. But this approach has some advantages from the numerical point of view, in particular for the
Stokes problem, as we can see later. This approach can be formalized in the following way.
Denote Q := L2 (Ω ; R2 ) and V := H01 (Ω ), as well as Λ u := ∇u and define the bilinear and linear forms
associated with the mixed formulation for all p, q ∈ Q and u, v ∈ V ,
Z Z
a(p, q) := p · qdx , b(q, u) := q · ∇udx , c(u, v) := 0 , (3.4.3)
Ω Ω
Z
lQ (q) := 0 lV (v) := f vdx . (3.4.4)
Ω
Thus the mixed formulation of the weak variational formulation of the Poisson problem can be considered
within the abstract framework.
Problem 3.1. Find (p, u) ∈ Q × V such that
a(p, q) − b(q, u) = lQ (q) ∀q ∈ Q , (3.4.5)

b(p, v) + c(u, v) = lV (v) ∀v ∈ V . (3.4.6)
Proposition 3.4.1 The function u solves the Poisson problem if and only if (u, ∇u) solves the mixed problem.
Proof. Observe that for q ∈ Q,
Z
a(p, q) − b(q, u) = (p − ∇u) · qdx = 0 = lQ (q) ,
Ω
and for v ∈ V , Z Z
b(p, v) + c(u, v) = p · ∇vdx = f vdx = lV (v) .
Ω Ω
We can also have formally a shorter notation for the mixed formulation with the bilinear form
A (U, P) := a(p, q) − b(q, u) + b(p, v) + c(u, v) ∀U := (p, u), P := (q, v) ∈ Q := Q × V ,
and the linear form

l(P) := lQ (q) + lV (v) .
Exercice 3.4.2 Show the existence of the solution to the mixed formulation of the Poisson problem by a
generalization of the Lax-Milgram lemma.
Remark 3.4.3 Let us note that there are α, β > 0 such that for P := (q, v) with q := p − ∇u and v := 2u we
have
α(kUkQ kPkQ ) ≤ α(kpkQ + kukQ )(kpkQ + β kukQ ) ≤ kpk2Q + kuk2Q = A (U, P) .
Mixed formulation of the nonhomogeneous Dirichlet boundary value problem for the shape derivatives
If we assume that the domain Ω moves, and the evolution of the domain within a tube is given by the family
Ωt , −δ < t < δ , then the solution of the solution of Poisson problem in the moving domain denoted by t →
3.4 Mixed formulation of elliptic boundary value problems 43
ut (y), y ∈ Ωt is differentiable, its derivative at t = 0 is denoted by u0 (x), x ∈ Ω solves the nonhomogeneous

Dirichlet problem
∂u
−∆ u0 = 0 in Ω , u0 = − (V · n) on ∂ Ω ,
∂n
where V · n stands for the normal component of the velocity vector field on ∂ Ω .
The problem for the shape derivative w := u0 with the nonhomogeneous Dirichlet boundary conditions
admits even a very week solution in L2 (Ω ) and it can be written in variational form
∂u ∂ϕ
Z Z
w ∈ L2 (Ω ) : w∆ ϕ = − (V · n) ∀ϕ ∈ H 2 (Ω ) ∩ H01 (Ω ) .
∂n ∂n
Ω Ω
Thus, the mixed formulation becomes, with p = −∇w :

Find w ∈ L2 (Ω ), p ∈ L2 (Ω ; div ) such that
Z
(p + ∇w) · q = 0 ∀q ∈ L2 (Ω ; div ) , (3.4.7)
Ω
Z
(div p − f )v = 0 ∀v ∈ L2 (Ω ) . (3.4.8)
Ω
3.4.1 Stokes problem
We denote by Π the projection mapping

1
Z
Π p := p − pdx .
measΩ
Ω
The starting point is the following elliptic boundary value problem written for the velocity u and the pressure
p of the fluid
−∆ u + ∇p = f in Ω , (3.4.9)
div u = 0 in Ω , (3.4.10)
u=0 on ∂ Ω , (3.4.11)
Π p = p. (3.4.12)
The proper choice of the spaces can be
p ∈ Q := {q ∈ L2 (Ω ) : Π q = q} , (3.4.13)
u ∈ V := {v = (v1 , v2 ), vi ∈ H01 (Ω ), i = 1, 2} . (3.4.14)
With this choice we can define the bilinear and linear forms as follows
Z
a(p, q) := − pq dx , (3.4.15)
Ω
Z
b(q, v) := − qdiv v dx , (3.4.16)
Ω
Z
c(u, v) := 2 ε(u) : ε(v) dx , ε(u) := ∇u + ∇u> ∈ L2 (Ω ; R4sym ) , (3.4.17)
Ω
Z
lQ (q) := 0 , lV (v) := f · vdx . (3.4.18)
Ω
44 3 Mixed formulation of elliptic boundary value problems
3.4.2 Linearized elasticity (Navier or Lamé problem)
The boundary value problem for linearized elasticity model is introduced in two spatial dimensions. The
shape optimization theory is well developed within the elasticity theory in three spatial dimensions, however
we restrict ourselves to the simple framework of plane problems.
2×2 reads
The constitutive relation in two spatial dimensions for the stress tensor σ (u) ∈ Rsym
σ (u) := C : ε(u) = λ tr ε(u)I + 2µε(u)
for all strains ε(u) ∈ R2×2

sym written as the symmetric 2 × 2-matrices, with Lamé parameters λ , µ > 0, and the
Hook’s tensor of material properties C, in the tensor or matrix notations. The inverse formula for the strain
tensor ε(u) = C−1 : σ (u) in the vector notation reads
1 λ
ε(u) = σ (u) + tr σ (u)I
2µ 4µ(λ + µ)
for all stress matrices σ ∈ R2×2

sym .
Remark 3.4.4 It is possible to use also the so-called Voigt notation for the strains and the stresses in the
form of column vectors in R3 for the case of two spatial dimensions. This simply means that the symmetric
2 × 2 matrices of strains and stresses which contains three independent components each, in the system of
two linear elliptic equations, are replaced by the vectors of the length 3. In this case the Hook’s tensor is
replace by the corresponding matrix or matrix function.
Let the body forces f = ( f1 , f2 ), fi ∈ L2 (Ω ) be given in the elastic body of the reference configuration Ω .
Assuming that the displacement field u solves the homogeneous elliptic boundary value problem in Ω
−div σ (u) = f in Ω , (3.4.19)

u=0 on ∂ Ω , (3.4.20)
we obtain the mixed formulation with the notation

Z
1
Z
λ
a(σ , τ) := σ : τ+ tr σ tr τ , b(τ, u) := τ : ε(u) , (3.4.21)
2µ 4µ(λ + µ)
Ω Ω
Z
c(u, v) := 0 , lV (v) := f ·v, lQ (τ) := 0 . (3.4.22)
Ω
The mixed formulation reads, find the functions σ ∈ Q and u ∈ V such that
 
 Z 
a(σ , τ) − b(τ, u) = 0 = (C−1 σ − ε(u)) : τ) ∀τ ∈ Q , (3.4.23)
 
Ω
 
 Z Z Z 
b(σ , v) = lV (v) = σ : ε(u) = − σ · v = f · v ∀v ∈ V . (3.4.24)
 
Ω Ω Ω
Remark 3.4.5 In this case the operator Q × V → (Q × V )? associated with the bilinear form A is linear,
bounded and bijective, therefore there is a unique solution of the mixed problem (3.4.23)-(3.4.24).
3.5 Piola transformation for finite elements
Gatica - Mixed Finite Element Method: str. 74-77

3.7 Mixed formulation and finite elements for elasticity 45
3.6 Piola transformation for Stokes problem and the shape derivative calculations
znalezc w literaturze Stokes, sprawdzic to co mam juz napisane
3.7 Mixed formulation and finite elements for elasticity
str. 121
Chapter 4
Lagrangian formalism for weak solutions and shape
optimization
4.1 Lagrangian formalism in shape optimization
The distributed gradient of the shape functional is an important tool in numerical methods of shape opti-
mization, see for example [2, 3]. It turns out that the boundary form of the shape gradient is less efficient in
the gradient type algorithms.
We are going to show that the distributed form of the shape gradient can be obtained by an application
of Lagrangian formalism. The general idea is the same as in the calculus of variations :
Z Z
u ∈ H01 (Ω ) : ∇u · ∇v dx = f v dx ∀v ∈ H01 (Ω )
Ω Ω
is the weak form of the homogeneous Dirichlet boundary value problem in Ω ⊂ D, where f ∈ L2 (D) is
given. The tracking type cost is
1
Z
J(Ω ) = (u − zd )2 dx,
2
Ω
where zd ∈ L2 (D) is given. We introduce the Lagrange function defined on D ⊂ Rd ,
1
Z Z Z
L (Ω , u, v) = (u − zd )2 + ∇u · ∇v dx − f v dx
2
Ω Ω Ω
for Ω ⊂ D, u, v ∈ H01 (Ω ).
The partial derivatives ∂u L (Ω , u, p)(ϕ) = 0, ∀ϕ ∈ H01 (Ω ) and ∂v L (Ω , u, v)(ϕ) = 0, ∀ϕ ∈ H01 (Ω ) de-
fines the adjoint state and the state equations, respectively. To define the derivative of Ω → L (Ω , u, v) we
can use the parametrization of domains by t → Ωt = Tt (V )(Ω ) which results in the function
1
Z Z Z
t → L (Ωt , u, v) := (u − zd )2 dx + ∇u · ∇v dx − f v dx,
2
Ωt Ωt Ωt
however, with u, v ∈ H01 (Ωt ).

Remark 4.1.1 In the case of nonhomogeneous Dirichlet boundary conditions u = ϕ on ∂ Ω the conditions
should be taken into account in definition of the Lagrangian. The result of the shape sensitivity analysis of
the boundary value problem in such a case is presented in Chapter ??.
In order to differentiate the latter function we make the change of variables Ωt = Tt−1 (Ω ) which leads
to new function defined in fixed domain Ω , for (u, v) ∈ H01 (Ω ) × H01 (Ω ),
1
Z Z Z
G (t, u, v) = γ(t)(u − ztd )2 dx + A(t)∇u · ∇v dx − γ(t) f t v dx,
2
Ω Ω Ω
where A(t) = det(DTt )DTt−1 ∗ DTt−1 .

We have the saddle-point property
47
48 4 Lagrangian formalism for weak solutions and shape optimization
1
Z
j(t) := G (t, ut , vt ) = min max G (t, u, v) = γ(t)(ut − ztd )2 dx
u∈H01 (Ω ) v∈H01 (Ω ) 2
Ω
The distributed shape gradient
d
j0 (0) = G (t, ut , pt )|t=0 = dJ(Ω ;V )
dt
depends upon the state and adjoint state but it is independent of the material derivatives of state and adjoint
state, namely, we evaluate the partial derivative of t → G (t, ·, ·),
dJ(Ω ;V ) = ∂t G (t, ut , pt )|t=0
since the material derivatives do exist u̇, ṗ ∈ H01 (Ω ), and we take into account that ∂t L (Ω , ut , v) =
∂u L (Ω , u, v)(u̇) = 0 as well as ∂t L (Ω , u, pt ) = ∂v L (Ω , u, v)( ṗ) = 0.
Remark 4.1.2 Let us note that
t → L (Ωt , ut , pt ) = min max L (Ωt , u, p),

u∈H01 (Ωt ) v∈H01 (Ωt )
however [127]
∂u
u0 = −(V · n) on ∂ Ω ,
∂n
and
∂p
p0 = −(V · n) on ∂ Ω ,
∂n
therefore, the shape derivatives u0 , p0 6∈ H01 (Ω ) and we cannot differentiate the Lagrangian in the variable
domain setting with ut , pt ∈ H01 (Ωt ).
4.1.1 Distributed shape gradients for Dirichlet problem by velocity method
We are going to describe the Lagrangian formalism for the shape differentiability of functionals. We restrict
ourselves to the homogeneous Dirichlet problem and the tracking type shape functional.
The Lagrangian method of differentiation is simple when applicable. All the details on the velocity
method [127] can be found in Chapter ??.
We set [127],
Ωt := Tt (V )(Ω ),
where the velocity field V is associated with the bilipschitzian transformation Tt : Rd → Rd , d = 1, 2, 3,
∂ Tt
V (t, x) = ◦ Tt−1 (x).
∂t
We denote t → γ(t) = det(DTt ) and t → A(t) = det(DTt )DTt−1 ∗ DTt−1 , so that [127], γ 0 (0)(x) = divV (0, x)
and A0 (0)(x) = divV (0, x) − DV (0, x) − ∗ DV (0, x) for d = 2, 3, for d = 1 the results are simpler. We recall
that DTt is the Jacobian matrix function for transfomation Tt and det(DTt ) stands for its determinant; ∗ DV
denotes the transposed of matrix DV .
Remark 4.1.3 In general the integral shape functional in variable domain t → Ωt

Z
t → j(t) := ϕ(t, x)dx
Ωt
admits the right-derivative at t = 0

Z Z
∂ϕ
j0 (0) = (0, x)dx + V (0, x) · n(x)dΓ (x), (4.1.1)
∂t
Ω Γ
4.1 Lagrangian formalism in shape optimization 49
where Γ ⊂ ∂ Ω is the moving portion of the boundary. Such a formula is also useful for the shape functional
depending on solutions to state equation in a variable domain setting. However, the straightforward evalu-
ation of the shape derivative can be performed in the fixed domain setting by using the material derivatives.
Thus, we rewrite Z
t → j(t) := γ(t)ϕ(t, Tt (x))dx
Ω
and differentiate the latter expression to obtain

Z
0 0 ∂ϕ
j (0) = γ (0)ϕ(0, x) + (0, x) + ∇ϕ(0, x) ·V (0, x) dx. (4.1.2)
∂t
Ω
It is easy to see, that (4.1.1) and (4.1.2) are equivalent by the Stokes Theorem.
4.1.2 Shape optimization for homogeneous Dirichlet problem in Ω ⊂ R2
We are going to explain, in the fixed domain setting in R2 , the Lagrangian formalism of the velocity method
in shape differentiability.
Let us assume that there are given zd , f ∈ L2 (D), where D ⊂ R2 is hold-all domain, with Ω ⊂ D, and
consider the shape optimization problem
Problem 4.1. Minimize Z
D 3 Ω → J(Ω ) = (u − zd )2 dx
Ω
over a family of admissible domains Ω ∈ Uad , subject to the state equation
−∆ u = f in Ω , (4.1.3)
u=0 on ∂ Ω . (4.1.4)
Velocity method for variable domain problem
In the variable domain t → Ωt = Tt (Ω ), the cost reads

Z
J(Ωt ) = (ut − zd )2 dx
Ωt
where ut solves Z Z
ut ∈ H01 (Ωt ) : ∇ut · ∇ψdx = f ψdx ∀ψ ∈ H01 (Ωt ).
Ωt Ωt
−1
Taking into account that Ω = Tt (Ωt ) we have the equivalent variational problem for the transported
function ut = ut ◦ Tt ,
Z Z
ut ∈ H01 (Ω ) : A(t)∇ut · ∇ψdx = γ(t) f t ψdx ∀ψ ∈ H01 (Ω ).
Ω Ω
The cost becomes Z

t → j(t) := J(Ωt ) = γ(t)(ut − ztd )2 dx
Ω
the function of one variable for a given velocity field associated with the transformation Tt . Thus, we can
consider the C1 function
Z Z Z
t → G(t, ϕ, ψ) = γ(t)(ϕ − ztd )2 dx + A(t)∇ϕ · ∇ψdx − γ(t) f t ψdx
Ω Ω Ω
50 4 Lagrangian formalism for weak solutions and shape optimization
for ϕ, ψ ∈ H01 (Ω ).
In order to introduce the saddle point formulation for the evaluation of the derivative
1
j0 (0) = lim ( j(t) − j(0)) ,
t→0 t
we introduce the specific adjoint state pt ∈ H01 (Ω ). The adjoint equation
∂ϕ G(t, ut , pt )(ψ) = 0 ∀ψ ∈ H01 (Ω ),
is obtained by the differentiation of ϕ → G(t, ϕ, ψ) at the point t, ϕ = ut , ψ = pt , thus

Z Z
pt ∈ H01 (Ω ) : A(t)∇pt · ∇ψdx = −2 γ(t)(ut − ztd )ψdx ∀ψ ∈ H01 (Ω ).
Ω Ω
Since the variational problems for the state ut ∈ H01 (Ω ) and the adjoint state pt ∈ H01 (Ω ) are given by
∂ϕ G(t, ut , pt )(ϕ) = 0, ∀ϕ ∈ H01 (Ω ) and ∂ψ G(t, ut , pt )(ψ) = 0, ∀ψ ∈ H01 (Ω ), respectively, hence for solu-
tions (ut , pt ) ∈ H01 (Ω ) × H01 (Ω ) of the state and adjoint equations the min-max theorem
G(t, ut , pt ) = min max G(t, ϕ, ψ)

ϕ∈H01 (Ω ) ψ∈H01 (Ω )
is satisfied and the derivative j0 (0) is a volume integral which is evaluated by the partial derivative
∂t G(t, ϕ, ψ) of t → G(t, ϕ, ψ) at the point (0, u, p) ∈ R × H01 (Ω ) × H01 (Ω ),
d
j0 (0) = G(t, ut , pt )|t=0 = ∂t G(0, u, p).
dt
in view of u̇ ∈ H01 (Ω ) and
∂ϕ G(t, u, p)(u̇) = 0
as well as ṗ ∈ H01 (Ω ) and
∂ψ G(t, u, p)( ṗ) = 0.
Distributed shape derivative of shape functional
In order to complete our analysis, we need an explicit formula for the shape derivative
j0 (0) = dJ(Ω ;V ).
It follows that
Z Z
dJ(Ω ;V ) = γ 0 (0)(u − zd )2 dx + 2 (u̇ + ∇zd ·V )(u − zd )dx. (4.1.5)
Ω Ω
and simple calculations show that the material derivative u̇ given by variational problem
Z Z
u̇ ∈ H01 (Ω ) : ∇u̇ · ∇ψ dx + A0 (0)∇u · ∇ψ dx = divV f ψ dx + ∇ f ·V ψ dx ∀ψ ∈ H01 (Ω )
Ω Ω
can be eliminated. As a result

Z
dJ(Ω ;V ) = divV [(u − zd )2 − f p] − 2(u − zd )∇zd ·V + A0 (0)∇u · ∇p − p∇ f ·V dx.
Ω
Now, we evaluate the boundary form of shape gradient. To this end the integration by parts is performed,
Z Z Z
divV [(u − zd )2 − f p] dx = − V · ∇[(u − zd )2 − f p] dx + [(u − zd )2 − f p](V · n)(x) dx,
Ω Ω ∂Ω
4.1 Lagrangian formalism in shape optimization 51
Z Z
A0 (0)∇u · ∇p dx = Vi,i u,l p,l −Vi, j u,i p, j −V j,i u,i p, j dx
Ω Ω
Z
=− Vi [u,l p,l ],i −Vi [u,i p, j ], j −V j [u,i p, j ],i dx
Ω
Z
+ Vi [u,l p,l ]ni −Vi [u,i p, j ]n j −V j [u,i p, j ]ni dx . (4.1.6)
∂Ω
Exercice 4.1.4 Taking into account the homogeneous Dirichlet boundary conditions and the Structure The-
orem of shape gradient [127] which says that the result depends exclusively on the normal component of
the velocity field, show that from the boundary integral in (4.1.6) we obtain
∂u ∂ p
Z
dJ(Ω ;V ) = − (V · n)(x) dx . (4.1.7)
∂n ∂n
∂Ω
Remark 4.1.5 If we introduce the adjoint state pe as a weak solution of the homogeneous Dirichlet problem
−∆ pe = zd − u
i.e., for the Lagrangian of the form

1
Z Z Z
L (Ω , u, v) = (u − zd )2 − ∇u · ∇v dx + f v dx .
2
Ω Ω Ω
Then we change the sign in (4.1.7)
∂ u ∂ pe
Z
dJ(Ω ;V ) = (V · n)(x) dx .
∂n ∂n
∂Ω
Chapter 5
The second-order shape differentiability
5.1 Introduction
The shape Hessian is introduceed in Section ?? in one spatial dimension.
5.2 Shape Hessian
Let us consider the homogeneous Dirichlet problem for the strictly positive right-hand side f > 0,
−∆ u = f in Ω , (5.2.1)
u=0 on ∂ Ω . (5.2.2)
The solution operator for the above boundary value problem is denoted by u := S ( f ).
The associated tracking type shape functional
1
Z
J(Ω ) = (u − zd )2 . (5.2.3)
2
Ω
The given element zd is defined by a solution of the homogeneous Dirichlet problem,
zd ∈ H01 (C(R, r)) : −∆ zd = f in C(R, r),
where, the reference domain C(R, r) := {r < ρ < R}, in polar coordinates (ρ, θ ). This function is extended
by zero outside of the ring C(R, r).
In this way there is a unique minimizer Ω ∗ := C(R, r) for the shape optimization problem under consid-
erations. The proof is left to the reader.
The shape gradient is given by
∂u ∂ p
Z
dJ(Ω ;V ) = V ·n (5.2.4)
∂n ∂n
∂Ω
with the adjoint state given by the homogeneous Dirichlet problem
−∆ p = u − zd in Ω , (5.2.5)
p=0 on ∂ Ω , (5.2.6)
We have in a ususal way the nonhomogeneous Dirichlet problems for the shape derivatives u0 , p0 of the
state u and of the adjoint state p, namely
53
54 5 The second-order shape differentiability
−∆ u0 = 0 in Ω , (5.2.7)
∂u
u0 = −V · n on ∂ Ω , (5.2.8)
∂n
and
−∆ p0 = u0 in Ω , (5.2.9)
∂p
p0 = −V · n on ∂ Ω , (5.2.10)
∂n
It is convenient to introduce the notation DN for the Dirichlet-to-Neumann or Steklov-Poincaré map, and
S for the inverse of negative Dirichlet-Laplacian. Thus, the mapping DN is defined by the the nonhomo-
geneous Dirichlet problem
∆w = 0 in Ω, w=ϕ on ∂Ω ⊂ ∂Ω
and it is DN : ϕ → ∂w
∂n =: DN (ϕ). In this notation we have
u = S [ f ], p = S [u − zd ]; (5.2.11)

∂u ∂u ∂p
u0 = −DN V ·n , p0 = −S ◦ DN V · n − DN V · n (5.2.12)
∂n ∂n ∂n
Now we want to evaluate the shape Hessian. To this end we differentiate with respect to t at t = 0 the shape
gradient written in Ωt
d ∂ ut ∂ pt
Z
d 2 J(Ω ;V,V ) := Vt · nt
dt ∂ nt ∂ nt
∂ Ωt
Proposition 5.2.1 The symmetric part of the shape Hessian is given by the formula
0
∂ u ∂ p ∂ u ∂ p0
Z
2
∂ J(Ω ;V,V ) := V · n + (5.2.13)
∂n ∂n ∂n ∂n
∂Ω

∂u ∂ p ∂ ∂u ∂ p
Z
(V · n)2 κ + (5.2.14)
∂n ∂n ∂n ∂n ∂n
∂Ω
5.3 Numerical evaluation of the shape Hessian
We use the explicit solution for the Laplacian in the ring C(R, r) = {r < ρ < R}
∆ u(x) = 0 in r2 < x12 + x22 < R2 } (5.3.1)

u(x) = f1 (x) on x12 + x22 = r2 , (5.3.2)
u(x) = f2 (x) on x12 + x22 = R2 , (5.3.3)
which is rewritten in polar coordinates
ρ 2 uρρ + ρuρ + uϕϕ = 0 for r < ρ < R, 0 ≤ ϕ < 2π , (5.3.4)

u(r, ϕ) = f1 (ϕ), u(R, ϕ) = f2 (ϕ) for 0 ≤ ϕ < 2π . (5.3.5)
The solution obtained by the separation of variables is a linear combinations of the elementary solutions
1, ln ρ, ρ n cos (nϕ), ρ n sin (nϕ), ρ −n cos (nϕ), ρ −n sin (nϕ), (5.3.6)
thus
∞
an ρ n + bn ρ −n cos (nϕ) + cn ρ n + dn ρ −n sin (nϕ)

u(ρ, ϕ) = a0 + b0 ln ρ + ∑ (5.3.7)
n=1
5.3 Numerical evaluation of the shape Hessian 55
Taking into account the nonhomogeneous Dirichlet boundary conditions at ρ = r combined with the expan-
sion of the function f1 into the Fourier series we obtain
∞
an rn + bn r−n cos (nϕ) + cn rn + dn r−n sin (nϕ)

u(r, ϕ) = a0 + b0 ln r + ∑ (5.3.8)
n=1
Z2π
1
a0 + b0 ln r = f1 (s)ds , (5.3.9)
2π
0
Zπ
1
an rn + bn r−n = f1 (s) cos (ns)ds , (5.3.10)
π
0
Zπ
n −n 1
cn r + dn r = f1 (s) sin (ns)ds , (5.3.11)
π
0
and replacing r by R in the same way we get

∞
an Rn + bn R−n cos (nϕ) + cn Rn + dn R−n sin (nϕ)

u(R, ϕ) = a0 + b0 ln R + ∑ (5.3.12)
n=1
Z2π
1
a0 + b0 ln R = f1 (s)ds , (5.3.13)
2π
0
Zπ
1
an Rn + bn R−n = f1 (s) cos (ns)ds , (5.3.14)
π
0
Zπ
1
cn Rn + dn R−n = f1 (s) sin (ns)ds , (5.3.15)
π
0
Therefore, we can compute the shape Hessian for a ring using this technique, since the normal derivative
∂
on the boundary is just the partial derivative ± .
∂ρ
Chapter 6
Shape optimization with MATLAB in 2D
We are interested in specific problems of shape optimization with the known, smooth solutions. Such solu-
tions are computed using the Matlab programs of optimization combined with the finite element approxi-
mations.
Dirichlet problem
Let D be a smooth hold all domain, and let us consider a simply connected domain Ω ⊂ D of the condenser
type with the boundary ∂ Ω which consists of two parts, the exterior boundary denoted ΓR and the interior
boudary denoted Γr . For example, the domain can be selected in the form of a ring, e.g., C(R, r) := BR \ Br .
We denote by S := D \ Ω the complement of Ω in D, hence S is a double connected compact domain.
For a given f > 0, f ∈ C(D), let us consider the homogeneous Dirichlet boundary value problem in the
simply connected domain Ω
−∆ u = f in Ω , (6.0.1)
u=0 on ∂ Ω . (6.0.2)
The solution of (6.0.1)-(6.0.2) for the ring C(R, r) is denoted by uC , and its extension by zero to D is denoted
by zd .
Problem 6.1. Consider the shape optimization problem
1
Z
inf J(Ω ) := (u − zd )2
2
Ω
over the family of simply connected domains in D with the volume bounded from below by the volume of
the ring |C(R, r)|.
Proposition 6.0.1 The unique solution of 6.1 is the ring C(R, r).
Proof. Existence Let Ωn be the minimizing sequence. We have J(C(R, r)) = 0, hence lim J(Ωn ) = 0.
Denote Sn := D \ Ωn → S in the Hausdorf distance, for a subsequence, with Ω := D \ S simply connected.
By the result [1] of Sverak H01 (Ωn ) → H01 (Ω ) in the sense of Kuratowski-Mosco. The solutions un ∈ H01 (Ωn )
in Ωn converge weakly in H01R(D) to u ∈ H01 (Ω ) ⊂ H01 (D) thus strongly in L2 (D). We have u > 0 in Ω by
the maximum principle, thus (u − zd )2 = 0 implies Ω ≡ C(R, r) and u ≡ zd in C(R, r).
Ω
u2 > 0.
R
Uniqueness For Ω 6= C(R, r) the cost functional J(Ω ) =
Ω \C(R,r)
6.1 Shape optimization with Poisson equation
6.1.1 Analysis
Let Ω ⊂ R2 with one hole. The objective functional reads as
57
58 6 Shape optimization with MATLAB in 2D
1 1
Z Z
J(Ω ) = k(u − zD )2 dx = (u − zd )2 dx (6.1.1)
2 Ω 2 E
where k is the characteristic function of the observation area E ⊂ Ω ⊂ D and u is solution of the state
equation )
−∆ u = f in Ω
(6.1.2)
u=0 on Γ
with a constant function f := b. The target function is denoted by zd and is given in polar coordinates (ρ, θ )
by the solution of the state equation on the ring C(r, R) := {x ∈ R2 : r ≤ |x| ≤ R}
(R2 −r2 ) 2 ln(r)−r2 ln(R)

(
b −ρ 2 + ln(R)−ln(r) ln(x) + R ln(R)−ln(r) if r ≤ ρ ≤ R
zd (ρ) = (6.1.3)
4 0 otherwise.
The weak form of the state equation reads with HΓ1 (Ω ) := {v ∈ H 1 (Ω ) | v|Γ = 0} as
)
Find u ∈ HΓ1 (Ω ) s.t
(6.1.4)
(∇u, ∇ϕ)Ω = ( f , ϕ)Ω ∀ϕ ∈ HΓ1 (Ω ).
The adjoint equation to the state equation follows by the Langrange functional
L(u, w, λ ) := J(Ω ) − (∇u, ∇w)Ω + ( f , w)Ω − (u, λ )Γ (6.1.5)
with w ∈ HΓ1 (Ω ) and its partial derivative with respect to u in any direction ū ∈ H 1 (Ω )
∂L
(u, w, λ )[ū] = δu L(u, w, λ )[ū] = δu J(Ω )[ū] − (∇ū, ∇w)Ω − (ū, λ )Γ . (6.1.6)
∂u
This derivative has to vanish in every direction and the adjoint equation holds
)
Find w ∈ HΓ1 (Ω ) s.t
(6.1.7)
(∇w, ∇ϕ)Ω + (λ , ϕ)Γ = du J(Ω )[ϕ] ∀ϕ ∈ H 1 (Ω ).
with Lagrage multiplier λ = −∇w · n and du J(Ω )[ϕ] = (k(u − zd ), ϕ)Ω . The shape derivative problem of
the Poisson equation reads as
Find v ∈ H 1 (Ω ) s.t



1 (6.1.8)
(∇v, ∇ϕ)Ω = 0 ∀ϕ ∈ HΓ (Ω )

v = −(∇u · n)(V · n)

on Γ
The shape gradient reads with the previous equations as

Z Z
δ J(Ω ;V ) = k(u − zd )v dx = (∇w · n)(∇u · n)(V · n) dω (6.1.9)
Ω Γ
6.1 Shape optimization with Poisson equation 59
6.1.2 Matlab Example
In this Matlab code shape optimization with the Poisson equation is implemented with shape gradient
method.
1 function SO2D_Poisson_spline
2 close all;clc
3
4 % folder for vtu output

5 runstr = fullfile(mfilename,datestr(now,'yy-mm-dd_HH-MM-SS'));
In this section the parameters tau as the initial stepsize, qual the area of the biggest triangle and f the
value of the right hand side are declared.
7 %% === parameters ======================================================== %
8 tau = 1e0; % initial stepsize
9 qual = 1e-3; % mesh quality
10 f = 1; % scalar value for RHS
As next the non-optimization boundary is described by one circles, which is stored in the struct inner.
Further a point is chosen to define a hole in the middle (holes).
12 %% === initial domain ==================================================== %
13 inner.N = 200; % number of points
14 inner.alpha = linspace(0,2*pi,inner.N+1);
15 inner.alpha(end) = [];
16 inner.radius = 1;
17 inner.p = ... % point coordinates
18 inner.radius * [cos(inner.alpha) ; sin(inner.alpha)];
19 inner.e = ... % edges
20 [1:length(inner.p);circshift((1:length(inner.p))',-1)'];
21
22 % define hole
23 holes = [0;0];
The optimization boundary is described by a spline (spline) with control points spline.b. In this
initial setting the control points are given by a disturbed circle (spline.radius). Furthermore the target
shape of the interface is stored in target to compute the Hausdorff distance, the point coordinates of the
spline are returned by the function ShapeGradToSpline(0,spline).
25 %% === Spline ============================================================ %
26 % initial control points
27 spline.N = 60; % number of control points
28 spline.N_inter = 6; % number of intermediate points
29 spline.N_total = spline.N * spline.N_inter;
30 spline.h_max = 0.05;
31 a = linspace(0,2*pi,spline.N+1);
32 a(end) = [];
33 % spline.radius = 2+0.0*sin(4*a)+0.0*cos(7*a);
34 spline.radius = 2+0.5*(2*rand(1,length(a))-1);
35
36 % control point coordinates
37 bx = spline.radius.*cos(a);
38 by = spline.radius.*sin(a);
39 spline.b = [bx;by];
40
41 % target shape
42 target.N = spline.N_total;
43 target.radius = 2;
44 target.alpha = linspace(0,2*pi,target.N+1);
45 target.alpha(end) = [];
46 target.p = target.radius * [cos(target.alpha) ; sin(target.alpha)];
47
48 % discretization points
49 spline = ShapeGradToSpline(0,spline);
Here the analytical solution is stored in a function handle zD in polar coordinates.

51 %% === analytical solution =============================================== %
52 r = inner.radius;
53 R = target.radius;
54 zD = @(x) f*(x.ˆ2 *(log(R) - log(r)) + Rˆ2*(log(r) - log(x)) ...
55 + rˆ2*(log(x) -log(R)))/(4*(log(r)-log(R))).*(x<=R);
Some initial definitions for the objective functional J, shape derivative dJ and Hausdorff distance dh
are added.
57 %% === shape gradient ==================================================== %
58 % init()
59 iter = 1;
60 J(1) = inf;
61 dJ(1) = inf;
62 dh(1) = inf;
63 rejected = false;
64 time(1) = 0;
Here starts the shape gradient algorithm. The content inside the while-loop is repeated until the shape
derivative or the Hausdorff distance reaches a bound. First the Hausdorff distance between the current shape
and the target shape is computed. Then the triangulation of the domain Ω is created with T RIANGLE inside
the function boundary2mesh. This function returns the points p, the edges e and the triangles t for P2
finite elements.
66 while (abs(dJ(end)) >=eps) && (dh(end) >=eps)
67 %% === triangulation ...
===================================================== %
68 dh(iter) = hausdorff(spline.p,target.p);
69
70 % Omega
71 p = [spline.p inner.p];
72 e = [spline.e inner.e+spline.N_total];
73 [p,e,t] = boundary2mesh(p,e,holes,qual);
In this section the primal problem solved with the finite element method. The stiffness matrix for the do-
main is returned by the function AssembleStiffness(p,t). The mass matrix to get a FE-description
of the right hand side is returned by AssembleMass(p,t). Furthermore the boundary conditions are
applied to the system.
75 %% === primal problem ...
==================================================== %
76 % Assemble system
77 A = AssembleStiffness_fast(p,t);
78 M = AssembleMass_fast(p,t);
79 F = f*sum(M,2);
80
81 % apply boundary condition

82 Ind = unique(e(1:3,:));
83 A(Ind,:) = 0;
84 A(Ind,Ind) = speye(length(Ind));
85 F(Ind) = 0;
86
87 % solve
88 u = A\F;
Here the adjoints system is solved. Before appling boundary conditions to the system, the objective
functional is evaluated.
90 %% === adjoint problem ...
=================================================== %
91 zDh = zD(norm2(p))';
92 F = M * (u-zDh);
93
94 J(iter) = 1/2 * (u-zDh)' *F;
95
96 F(Ind) = 0;
97 w = A\F;
To enable a changing step size, a stepsize control (scc(...)) is called. It decides whether a iteration
is accepted or rejected and when it is accepted the stepsize should be increased, decreased or hold.
6.1 Shape optimization with Poisson equation 61
99 %% === step size control ...

================================================= %
100 if iter > 1
101 [tau,iter,rejected] = ssc(tau,dJ,J,iter);
102 sfigure(3)
103 plot(time,J)
104 sfigure(4)
105 plot(time,dh)
106 end
107
108 if rejected
109 spline = old; % take old domain
110 else
111 old = spline; % store current domain
112 end
The shape gradient is evaluted with the solution of the primal and adjoints system at each point with
index stored in e. The functions nablaU_e(p,e,t,u) returns the derivative of u ∇u on the boundary
of the domain (p,t) at the point with index contained in e.
114 %% === shape gradient ...
==================================================== %
115 if ˜rejected
116 % nabla u and nabla w
117 DU = nablaU_e(p,e,t,u);
118 DW = nablaU_e(p,e,t,w);
119
120 idx_ef = 1:(spline.N_total);

121 tangent = p(:,e(2,idx_ef))-p(:,e(1,idx_ef));
122 normal = [tangent(2,:);-tangent(1,:)];
123 normal = bsxfun(@times,normal,1./norm2(normal));
124
125 g = -sum(DU(5:6,idx_ef).*normal,1).*sum(DW(5:6,idx_ef).*normal,1);
126 dJ(iter) = -sum(norm2(p(:,e(1,idx_ef))-p(:,e(2,idx_ef))).*g.ˆ2)*tau;
127 end
The result is written to a vtu-file for use in PARAVIEW with the function WriteBinaryVTK.
129 %% === Write VTK ...
======================================================== %
130 if ˜rejected
131 WriteBinaryVTK(fullfile('Results',runstr,sprintf('2D_Shape_%04d.vtu',iter)),...
132 p,t(1:3,:) ,...
133 'U', full(u') ,...
134 'W', full(w') ,...
135 'error',full(u'-zDh') ,...
136 'J', J(iter) +0*p(1,:) ,...
137 'dh', dh(iter) +0*p(1,:) ,...
138 'zD', zDh' );
139 end
In this last part the computed shape gradient is used to calculate the displacement of the spline control
points. Further the displacement of the control points is restricted to the radial direction, thus one has to
solve the problem
! 2

2
min kg −V · nk = min g − ∑ ci bi ϕi · n

V c
i
where ϕ are the spline basis functions, b the control points and n the outer normal vector. This is done
with the function ShapeGradToSpline, which returnes the new point coordinates, control points and
normals.
141 %% === Shape gradient - SPLINE VERSION ...
================================== %
142 [spline] = ShapeGradToSpline(g*tau,spline);
143
144 iter = iter+1;

145 time(iter) = time(iter-1)+tau;
146 end
147
148 end
6.2 Interface problem 63
6.2 Interface problem
6.2.1 Analysis
Two domains Ωα and Ωβ with interface layer Σ which is the optimization

boundary. 
−aα ∆ uα = 0 in Ωα 


uα = 0 on Γα 




−aβ ∆ uβ = 0 in Ωβ 

(6.2.1)
aβ δn β uβ = f on Γβ 


uα = uβ

on Σ 




aα δnα uα = −aβ δnβ uβ on Σ

The objective functional reads as

1 1
Z Z
J(Σ ) = (uα − ud )2 dx + (uβ − ud )2 dx. (6.2.2)
2 Ωα 2 Ωβ
The target function for ring C(r, R) is given by

(
Rf aβ (ln(ρ) − ln(r)) if ρ ∈ [r, S)
ud (x, y) = ud (ρ(x, y)) = (6.2.3)
aα aβ aβ (ln(S) − ln(r)) + aα (ln(ρ) − ln(S)) if ρ ∈ [S, R]
The adjoint system for the interface problem reads as


−aα ∆ wα = uα − ud in Ωα 


wα = 0 on Γα 




−aβ ∆ wβ = uβ − ud in Ωβ 

(6.2.4)
aβ δn wβ = 0 on Γβ 


wα = wβ

on Σ 




aα δnα wα = −aβ δnβ wβ on Σ

and its weak form is given by


aα (∇wα , ∇ϕα )Ωα + (λ1 , ϕα )Γα + (λ2 , ϕα )Σ = (uα − ud , ϕα ) ∀ϕα ∈ H 1 (Ωα ) 


aβ (∇wβ , ∇ϕβ )Ωβ − (λ2 , ϕβ )Σ = (uβ − ud , ϕβ ) ∀ϕα ∈ H 1 (Ωβ ) (6.2.5)


wα = wβ on Σ 
with Lagrange multipliers

)
λ1 = −aα δnα wα on Γα
(6.2.6)
λ2 = −aα δnα wα = aβ δnβ wβ on Σ
The weak shape derivatives system without simplifications is given by
aα (∇u0α , ∇ϕα )Ωα = aα (δnα uα , ∇ϕα ·V )δ Ωα − aα (∇uα · ∇ϕα ,V · nα )δ Ωα




aβ (∇u0β , ∇ϕβ )Ωβ = aβ (δnβ uβ , ∇ϕβ ·V )δ Ωβ − aβ (∇uβ · ∇ϕβ ,V · nβ )δ Ωβ






+ (δnβ ( f ϕβ ) + f ϕβ κ,V · nβ )δ Ωβ (6.2.7)

u0α

on Γα = −(δnα uα )(V · nα ) 



u0α = u0β + ∇uβ ·V − ∇uα ·V

on Σ 
Thus the shape gradient reads as

δ J(Ω ;V ) = (uα − ud , u0α )Ωα + (uβ − ud , u0β )Ωβ (6.2.8)

= +aα ((δnα uα )(δnα wα ),V · nα )Σ − aα ((δτα uα )(δτα wα ),V · nα )Σ (6.2.9)
− aβ ((δnβ uβ )(δnβ wβ ),V · nα )Σ + aβ ((δτβ uβ )(δτβ wβ ),V · nα )Σ (6.2.10)
In this Matlab code the optimization of the interface is implemented by a shape gradient method.
1 function SO2D_Interface_spline
2 close all;clc;
3

In this section the parameters tau as the initial stepsize, qual the area of the biggest triangle, f the value
of the inhomogeneous Neumann boundary, a_alpha and a_beta the values for the domain parameters
aα , aβ are declared.
10 %% === parameters ======================================================== %
11 tau = 1e-1; % initial stepsize
13 f = 1; % inhomogenouos Neumann boundary value
14 a_alpha = 1; % domain parameters
15 a_beta = 2;
As next the non-optimization boundary is described by two circles , which are stored in two structs
inner and outer. Further a point is chosen to define a hole in the middle (holes).
16 %% === initial domain ==================================================== %
17 % inner boundary
23 inner.radius * [cos(inner.alpha); sin(inner.alpha)];
25 [1:length(inner.p); circshift((1:length(inner.p))',-1)'];
26
27 % outer boundary
28 outer.N = 200;
29 outer.alpha = linspace(0,2*pi,outer.N+1);
30 outer.alpha(end) = [];
31 outer.radius = 2;
32 outer.p = ...
33 outer.radius * [cos(outer.alpha) ; sin(outer.alpha)];
34 outer.e = ...
35 [1:length(outer.p);circshift((1:length(outer.p))',-1)'];
36
37 % define hole
38 holes = [0;0];
39 %% === Spline ============================================================ %
40 % initial interface defined by spline
46 a(end) = [];
47 % spline.radius = 1.5+0.2*sin(12*a)+0.0*cos(7*a);
48 spline.radius = 1.5+0.3*(2*rand(1,length(a))-1);
49
54
55 % target shape
57 target.radius = 1.5;
60 target.p = target.radius * [cos(target.alpha) ; sin(target.alpha)];
61
63 [spline] = ShapeGradToSpline(0,spline);
Here the analytical solution is stored in a function handle zD in polar coordinates.

67 R = outer.radius;
69 S = target.radius;
70
71 zD = @(x) R*f/(a_beta*a_alpha)*( ...
72 (x<S).*(a_beta*(log(x)-log(r))) + ...
73 (x>=S).*((a_beta*(log(S)-log(r)))+a_alpha*(log(x)-log(S))) ...
74 );
are added.
76 %% === shape gradient ==================================================== %
77 % init()
78 iter = 1;
79 J(1) = inf;
80 dJ(1) = inf;
82 time(1) = 0;
83 dh(1) = inf;
derivative or the Hausdorff distance reaches a bound. First the Hausdorff distance between the current
shape and the target shape is computed. Then the triangulation of each subdomain Ωα and Ωβ is created
with T RIANGLE inside the function boundary2mesh. This function returns the points p, the edges e
and the triangles t for P2 finite elements. At last the subdomains are combined and the double entries are
erased. Variable names with capital letters (e.g. P, T) mark the relation to the combined domain.
===================================================== %
88
89 % Omega_alpha
90 p1 = [spline.p inner.p];
91 e1 = [spline.e inner.e+spline.N_total];
92 [p1,e1,t1] = boundary2mesh(p1,e1,holes,qual);
93
94 % Omega_beta
95 p2 = [spline.p outer.p];
96 e2 = [spline.e outer.e+spline.N_total];
97 [p2,e2,t2] = boundary2mesh(p2,e2,holes,qual);
98
99 % combination
100 dblidx = unique(e2(1:3,1:length(spline.e)));
101 p2_idx = 1:length(p2);
102 not_dblidx = setdiff(p2_idx,dblidx);
103 trans = p2_idx;
104 trans(not_dblidx) = (1:length(not_dblidx)) + length(p1);
105
106 P = [p1,p2(:,not_dblidx)];
107 T = [t1,trans(t2)];
108 ID = [ones(1,length(t1)),zeros(1,length(t2))];
In this section the primal problem solved with the finite element method. The stiffness matrix for the
domain represented by points p and triangles t is returned by the function AssembleStiffness(p,t).
Furthermore the boundary conditions and interface conditions are applied to the system.

==================================================== %
112 A1 = a_alpha*AssembleStiffness_fast(p1,t1);
113 F1 = zeros(length(p1),1);
114 A2 = a_beta*AssembleStiffness_fast(p2,t2);
115 F2 = zeros(length(p2),1);
116

118 dirichlet = unique(e1(1:3,(spline.N_total+1):end));
119 neumann = unique(e2(1:3,(spline.N_total+1):end));
120 i_neumann = (spline.N_total+1):length(e2);
121 interface1 = unique(e1(1:3,1:spline.N_total));
122 interface2 = unique(e2(1:3,1:spline.N_total));
123
124 A1(dirichlet,:) = 0;
125 A1(dirichlet,dirichlet) = speye(length(dirichlet));
126 F1(dirichlet) = 0;
127
128 le = norm2(p2(:,e2(1,i_neumann))-p2(:,e2(2,i_neumann)))';
129 F2(e2(1,i_neumann)) = F2(e2(1,i_neumann))+1/6.*le.*f;
132
133 A_flux = [A1(interface1,:),A2(interface2,:)];

134 F_flux = zeros(length(interface1),1);
135
136 A1(interface1,:) = 0;
137 A1(interface1,interface1) = speye(length(interface1));
138 A2(interface2,:) = 0;
139 A2(interface2,interface2) = speye(length(interface2));
140 A_cont = [A1(interface1,:),-A2(interface2,:)];
141 F_cont = zeros(length(interface1),1);
142
143 F1(interface1) = [];
145 A1(interface1,:) = [];
146 A2(interface2,:) = [];
147
148 A = [blkdiag(A1,A2); A_flux; A_cont];

149 F = [F1; F2; F_flux; F_cont];
150
151 % solve
152 u = A\F;
153 u1 = u(1:length(p1));
154 u2 = u(length(p1)+1:end);
155 U = [u1;u2(not_dblidx)];
Here the adjoints system is solved. The mass matrix to get a FE-description of the right hand side is
returned by AssembleMass(p,t). Before appling boundary conditions to the system, the objective
=================================================== %
158 zDh = zD(norm2([p1,p2]))';
159 error = u-zDh;
160 err1 = error(1:length(p1));
161 err2 = error(length(p1)+1:end);
162
163 ERROR = [err1;err2(not_dblidx)];

164 ZDH = zD(norm2(P))';
165
166 M1 = AssembleMass_fast(p1,t1);
167 M2 = AssembleMass_fast(p2,t2);
168 F1 = M1 * err1;
169 F2 = M2 * err2;
170
171 J(iter) = 1/2 * err1'*F1 + 1/2 * err2'*F2;

172
173 F1(dirichlet) = 0;
174 F2(neumann) = 0;
175

178
179 F = [F1;F2;F_flux;F_cont];
180
181 w = A\F;
182 w1 = w(1:length(p1));
183 w2 = w(length(p1)+1:end);
184 W = [w1;w2(not_dblidx)];
To enable a changing step size a stepsize control (scc(...)) is called. It decides whether a iteration is
accepted or rejected and when it is accepted the stepsize should be increased, decreased or hold.
================================================= %
187 if iter > 1
189 sfigure(3)
190 plot(time,J)
191 sfigure(4)
192 plot(time,dh)
193 end
194
195 if rejected
197 else
199 end
==================================================== %
203 if ˜rejected
204 Du1 = nablaU_e(p1,e1,t1,u1);
205 Du2 = nablaU_e(p2,e2,t2,u2);
206 Dw1 = nablaU_e(p1,e1,t1,w1);
207 Dw2 = nablaU_e(p2,e2,t2,w2);
208 x = [1 3 5];
209 y = [2 4 6];
210
211 idx_ef = 1:(spline.N_total); % indices of optimization points
212 tangent = p1(:,e1(2,idx_ef))-p1(:,e1(1,idx_ef));
215 gs = zeros(3,length(idx_ef));
216 for i = 1:3
217 gs(i,:) = ...
218 + a_alpha * (dot(Du1([x(i);y(i)],idx_ef),normal).* ...
219 dot(Dw1([x(i);y(i)],idx_ef),normal)) ...
220 - a_alpha * (dot(Du1([x(i);y(i)],idx_ef),tangent).* ...
221 dot(Dw1([x(i);y(i)],idx_ef),tangent)) ...
222 - a_beta * (dot(Du2([x(i);y(i)],idx_ef),normal).* ...
223 dot(Dw2([x(i);y(i)],idx_ef),normal)) ...
224 + a_beta * (dot(Du2([x(i);y(i)],idx_ef),tangent).* ...
225 dot(Dw2([x(i);y(i)],idx_ef),tangent));
226 end
227 g = -gs(3,:);
228 le = norm2(p1(:,e1(1,idx_ef))-p1(:,e1(2,idx_ef)));
229 dJ(iter) = -sum(le.*g.ˆ2)*tau;
230 end
231 %% === Write VTK ...

======================================================== %
232 if ˜rejected
233 WriteBinaryVTK(fullfile('Results',runstr,sprintf('2D_Shape_%04d.vtu',iter)),...
234 P,T(1:3,:) ,...
235 'U', full(U') ,...
236 'W', full(W') ,...
237 'error',full(ERROR') ,...
238 'J', J(iter) +0*P(1,:) ,...
239 'dh', dh(iter) +0*P(1,:) ,...
240 'zD', ZDH' ,...
241 'ID', ID );
242 end
solve the problem
! 2

2

V c
i
normals.
================================== %
245
246 iter = iter+1;

248 drawnow
249 end
250
251 end
6.3 Shape optimization with linear elasticity
6.3.1 Analysis
Let Ω ⊂ R2 with one hole. The objective functional read as

1 1
Z Z
J(Ω ) = k(u − zD )2 dx = (u − zD )2 dx (6.3.1)
2 Ω 2 E
with u solution of the state equation


−div (σ (u)) = f in Ω 


u=0 on ΓD 



σ (u) · n = 0 on ΓN (6.3.2)

σ (u) = λ tr (ε(u))I + 2µε(u) 




ε(u) = 21 (∇u + ∇uT )

and k characteristic function of the observeration area E ⊂ Ω ⊂ D. The target function is denoted by zD
(r − ρ) − 2R2 (λ + µ) + Rrλ + (2rµ + Rλ )ρ

x
zD (x, y) = − 2 2 2
(6.3.3)
2ρ (λ + 2µ)(R (λ + µ) + r µ) y
which is solution on the domain Ω := BR (0) \ Br (0) with right hand side

1 x p
f (x, y) = − 3 , ρ = x2 + y2 (6.3.4)
ρ y
and extended over BR (0) with contant value. The weak form of the state equation reads with V0 := {v ∈
H 1 (Ω )2 | v|ΓD = 0} as
)
Find u ∈ V0 s.t
(6.3.5)
(σ (u), ε(ϕ))Ω = ( f , ϕ)Ω ∀ϕ ∈ V0 .
The adjoint equation to the state equation follows by the Langrange functional
L(u, w, λ ) := J(Ω ) − (σ (u), ε(w))Ω + ( f , w)Ω − (u, λ )ΓD (6.3.6)
with w ∈ V0 and its derivative with respect to u in direction ū ∈ H 1 (Ω )
δu L(u, w)[ū] = δu J(Ω )[ū] − (σ (ū), ε(w))Ω − (ū, λ )ΓD . (6.3.7)
This derivative has to vanish in every direction and the adjoint equation holds
)
Find w ∈ V0 s.t
(6.3.8)
(σ (w), ε(ϕ))Ω + (λ , ϕ)ΓD = du J(Ω )[ϕ] ∀ϕ ∈ H 1 (Ω ).
with Lagrage multiplier λ = −σ (w) · n and du J(Ω )[ϕ] = (k(u − zD ), ϕ)Ω .

The shape derivative problem of the linear elasticity equation reads as
Find v ∈ H 1 (Ω ) s.t



(σ (v), ε(ϕ))Ω = −((σ (u) · τ) · (∇ϕ · τ), V · n)ΓN


(6.3.9)
+ ( f · ϕ, V · n)Γ ∀ϕ ∈ V0 



v = −(∇u · n)(V · n) on ΓD
because

d
dt (σ (ut ), ∇ϕt )Ωt t=0 = (σ (v), ∇ϕ)Ω + (σ (u), ∇ϕ 0 )Ω + (σ (u) : ∇ϕ, V · n)Γ (6.3.10)

d
dt ( f , ϕt )Ωt = ( f , ϕ 0 )Ω + ( f · ϕ, V · n)Γ (6.3.11)

t=0
6.3 Shape optimization with linear elasticity 71
The shape gradient reads with the previous equations as

Z
dJ(Ω ; V ) = k(u − zD )v dx (6.3.12)
ZΩ
= (σ (w) : ε(u))(V · n) dω (6.3.13)
ΓD
Z
+ (−(σ (w) : ε(u)) + f · w)(V · n) dω (6.3.14)
ΓN
In this Matlab code shape optimization on the Neumann boundary with the linear elasticity is implemented
with shape gradient method.
1 function SO2D_LE_spline_neumann
2 close all;clc
3

In this section the parameters tau as the initial stepsize, qual the area of the biggest triangle, Youngs
modulus E and Poissons ratio nu are declared.
7 %% === parameters ======================================================== %
8 tau = 1e4; % initial stepsize
10 E = 69; % material constants
11 nu = 0.33;
As next the non-optimization boundary is described by one circles, which is stored in the struct inner.
Further a point is chosen to define a hole in the middle (holes).
12 %% === initial domain ==================================================== %
18 inner.radius*[cos(inner.alpha) ; sin(inner.alpha)];
20 [1:length(inner.p);circshift((1:length(inner.p))',-1)'];
21
22 % define holes
23 holes = [0;0];
24 %% === Spline ============================================================ %
25 % initial control points
31 a(end) = [];
32 % spline.radius = 2+0.2*sin(4*a)+0.2*cos(7*a);
33 spline.radius = 2+0.4*(2*rand(1,length(a))-1);
34

39
40 % target shape
42 target.radius = 2;
45 target.p = target.radius*[cos(target.alpha) ; sin(target.alpha)];
46
48 [spline] = ShapeGradToSpline(0,spline);
Here the analytical solution is stored in a function handle zD and the right hand side in f.
51 R = target.radius;
52 lambda = E*nu/((1+nu)*(1-2*nu));
53 mu = E/(2*(1+nu));
54 f = @(p) -bsxfun(@times,p,1./norm2(p).ˆ3);
55 zD = @(p) -bsxfun(@times,p,...
56 ((r-norm2(p)).*(Rˆ2*lambda*-2.0-Rˆ2*mu*2.0+mu*r*norm2(p)*2.0+R*lambda*r...
57 +R*lambda*norm2(p))*(1.0/2.0))...
58 ./(norm2(p).ˆ2*(lambda+mu*2.0)*(Rˆ2*lambda+Rˆ2*mu+mu*rˆ2)).*(norm2(p)<=R)...
59 + ((r-R).*(Rˆ2*lambda*-2.0-Rˆ2*mu*2.0+mu*r*R*2.0+R*lambda*r...
60 +R*lambda*R)*(1.0/2.0))...
61 ./(R.ˆ2*(lambda+mu*2.0)*(Rˆ2*lambda+Rˆ2*mu+mu*rˆ2)).*(norm2(p)>R));
are added.
63 %% === shape gradient ==================================================== %
64 % init()
65 iter = 1;
66 J(1) = inf;
67 dJ(1) = inf;
69 time(1) = 0;
70 dh(1) = inf;
derivative or the Hausdorff distance reaches a bound. First the Hausdorff distance between the current shape
and the target shape is computed. Then the triangulation of the domain Ω is created with T RIANGLE inside
the function boundary2mesh. This function returns the points p, the edges e and the triangles t for P2
finite elements.
===================================================== %
75
76 %\Omega
77 p = [spline.p inner.p];
78 e = [spline.e inner.e+spline.N_total];
79 [p,e,t] = boundary2mesh(p,e,holes,qual);
In this section the primal problem solved with the finite element method. The linear elasticity matrix
for the domain represented by points p and triangles t is returned by the function AssembleLE(p,t).
The mass matrix to get a FE-description of the right hand side is returned by AssembleMass(p,t).
Furthermore the boundary conditions are applied to the system.
==================================================== %
83 fh = f(p);
84 A = AssembleLE_fast(p,t);
85 M = AssembleMass_fast(p,t);
86 F1 = M * fh(1,:)';
87 F2 = M * fh(2,:)';
88 F = [F1;F2];
89
91 outer_d = unique(e(1:3,1:spline.N_total));
92 inner_d = unique(e(1:3,spline.N_total+(1:inner.N)));
93 NoQDoF = size(p,2);
94
95 outer_d = [outer_d;outer_d + NoQDoF];

96 inner_d = [inner_d;inner_d + NoQDoF];
97
98 A(inner_d,:) = 0;
99 A(inner_d,inner_d) = speye(length(inner_d));
100 F(inner_d) = 0;
101 F(outer_d) = 0;
102
103 % solve
104 u = A\F;
105 U = reshape(u,[],2)';
Here the adjoints system is solved. Before appling boundary conditions to the system, the objective
=================================================== %
108 zDh = zD(p);
109 error = U-zDh;
110
111 F1 = M * error(1,:)';
112 F2 = M * error(2,:)';
113 F = [F1;F2];
114
115 J(iter) = 1/2 * reshape(error',[],1)'*F ;
116
117 F(outer_d) = 0;
118 F(inner_d) = 0;
119
120 w = A\F;
121 W = reshape(w,[],2)';
To enable a changing step size, a stepsize control (scc(...)) is called. It decides whether a iteration
is accepted or rejected and when it is accepted the stepsize should be increased, decreased or kept.
================================================= %
124 if iter > 1
126 sfigure(3)
127 plot(time,J)
128 sfigure(4)
129 plot(time,dh)
130 end
131
132 if rejected
134 else
136 end
==================================================== %
139 if ˜rejected
141 Du1 = nablaU_e(p,e,t,U(1,:)');
142 Du2 = nablaU_e(p,e,t,U(2,:)');
143 Dw1 = nablaU_e(p,e,t,W(1,:)');
144 Dw2 = nablaU_e(p,e,t,W(2,:)');
145 x = [1 3 5];
146 y = [2 4 6];
147
148 idx_ef = 1:(spline.N_total); % indices of optimization points

149 tangent = p(:,e(2,idx_ef))-p(:,e(1,idx_ef));
150 tangent = bsxfun(@times,tangent,1./norm2(tangent));
153 gs = zeros(3,length(idx_ef));
154 for i = 1:3
155 DU = [Du1([x(i),y(i)],:);Du2([x(i),y(i)],:)];
156 DP = [Dw1([x(i),y(i)],:);Dw2([x(i),y(i)],:)];
157 epsilon_P = 1/2 *(DP + DP([1 3 2 4],:));
158 epsilon_U = 1/2 *(DU + DU([1 3 2 4],:));
159 sigma = 2*mu*epsilon_U;
160 sigma([1 4],:) = bsxfun(@plus,sigma([1 4],:),...
161 lambda * sum(epsilon_U([1 4],:),1));
162 sigma_t = [sum(sigma([1 2],idx_ef).*tangent,1);...

163 sum(sigma([3 4],idx_ef).*tangent,1)];
164 Dpt = [sum(DP([1 2],idx_ef).*tangent,1);...
165 sum(DP([3 4],idx_ef).*tangent,1)];
166 gs(i,idx_ef) = - dot(sigma_t,Dpt)...
167 + dot(W(:,e(i,idx_ef)),fh(:,e(i,idx_ef)));
168 end
169 g = -gs(3,idx_ef);
170 dJ(iter) = -sum(norm2(p(:,e(1,idx_ef))-p(:,e(2,idx_ef))).*g.ˆ2)*tau;
171 end
173 %% === Write VTK ...
======================================================== %
174 if ˜rejected
175 WriteBinaryVTK(...
176 fullfile('Results',runstr,sprintf('2D_Shape_%04d.vtu',iter)),...
177 p,t(1:3,:) ,...
178 'primal', full(U) ,...
179 'adjoint', full(W) ,...
180 'analytical', zDh ,...
181 'error', error ,...
182 'J', J(iter) +0*p(1,:) ,...
183 'dh', dh(iter) +0*p(1,:) );
184 end
solve the problem
! 2

2

V c
i
normals.
================================== %
187
188 iter = iter+1;

190 drawnow
191 end
192
193 end
6.4 Example Shape Hessian
6.4.1 Setting
Let Ω = C(R, r) := {ρ ∈ R2 : r < ρ < R} a ring with the two radii R and r. Futhermore let Ω ⊂ B with a
hold-all B = C(0, RB ). Consider the Poisson equation on Ω with solution u
)
−∆ u = 1 in Ω
u=0 on δ Ω
and the objective functional Z

J(Ω ) = (u − zD )2 dx.
Ω
where the target function is contructed by
(
u∗ in Ω ∗
zD =
0 in B \ Ω ∗
with solution u∗ of the state equation on Ω ∗ = C(r∗ , R∗ ). The analytical solution of the state equation reads
as
u(ρ) = Ψ (ρ) − aΨ (r, R) ln(ρ) + bΨ (r, R)
with the particular solution
1
Ψ (ρ) = − ρ 2 .
4
The coefficients
Ψ (R) −Ψ (r) Ψ (R) ln(r) −Ψ (r) ln(R)
aΨ (r, R) = , bΨ (r, R) =
ln(R) − ln(r) ln(R) − ln(r)
are given by solving the harmonic extension with fourier series of a function Ψ on the domain C(r, R). For
simplification the coefficients are denoted in the following way:
aΨ (r, R) =: aΨ =: a, bΨ (r, R) =: bΨ =: b,
∗ ∗ ∗ ∗
aΨ (r , R ) =: aΨ =: a , bΨ (r∗ , R∗ ) =: bΨ
∗
=: b∗ .
6.4.2 Adjoint equation
The adjoint equation to the Poisson equation reads as:

)
−∆ p = u − zD in Ω
p=0 on δ Ω
The adjoint p is decomposed in the following way:
p = p̃ + η
where p̃ solves )
−∆ p̃ = −zD in Ω
p̃ = 0 on δ Ω
The solution p̃ is obtained by the particular solution P̃ whichs solves
−∆ P̃ = −zD in Ω
and reads
6.4 Example Shape Hessian 77

∗ on Ω ∗
p

P̃ = pr on {ρ < r∗ }
on {ρ > R∗ }

pR

with
1 4 1 2 ∗
p∗ (ρ) = − ρ + ρ (a (1 − ln(ρ)) + b∗ )
64 4
pr (ρ) = p∗ (r∗ ) + r∗ δρ p∗ (r∗ )(ln(ρ) − ln(r∗ ))
pR (ρ) = p∗ (R∗ ) + R∗ δρ p∗ (R∗ )(ln(ρ) − ln(R∗ ))
1 1
δρ p∗ = − ρ 3 + ρ(a∗ ( 21 − ln(ρ)) + b∗ ).
16 2
We get p̃ by harmonic extension of the function P̃:
p̃(ρ) = P̃(ρ) − aP̃ ln(ρ) + bP̃
The second part η of the adjoint solves

)
−∆ η = u in Ω
η =0 on δ Ω .
With the particular solution N, which solves
−∆ N = u in Ω
and reads
1 4 1 2
N(ρ) = ρ − ρ (a(1 − ln(ρ)) + b),
64 4
we have
η(ρ) = N(ρ) − aN ln(ρ) + bN .
6.4.3 Shape Gradient
In the following we use the notation (

−α on Γr
V ·n =
β on ΓR .
The shape gradient reads as:

δu δ p δu δp δu δu
Z
δ J(Ω ;V ) = (V · n) dω = 2π −αr (r) (r) + β R (R) (R)
δΩ δn δn δρ δρ δρ δρ
Thus the shape gradient can be written as

α
δ J(Ω ,V ) = 2π∇F
β
with
δu δp δu δp
(∇F)1 = −r (r) (r) (∇F)2 = R (R) (R)
δρ δρ δρ δρ
and
δu δΨ a δΨ 1
= − =− ρ
δρ δρ ρ δρ 2
δp δ p̃ δ η δ p̃ δ P̃ aP̃
= + = −
δρ δρ δρ δρ δρ ρ
δρ p∗ on Ω ∗

δ P̃  ∗ δ p∗ 1 1
= r δρ p∗ (r∗ ) ρ1 on {ρ < r∗ } = − ρ 3 + ρ(a∗ ( 21 − ln(ρ)) + b∗ )
δρ   ∗ δρ 16 2
R δρ p∗ (R∗ ) ρ1 on {ρ > R∗ }
δη δ N aN δN 1 1
= − = ρ 3 − ρ(a( 12 − ln(ρ)) + b)
δρ δρ ρ δρ 16 2
6.4.4 Shape Hessian
The symmetric part of the shape Hessian is given by the formula:

0
δ u δ p0

δu δp
∂ 2 J(Ω ;V,V ) := − 2παr (r) (r) + (r) (r)
δρ δρ δρ δρ
2 δ2p

1 δu δp δ u δp δu
+ 2πα 2 r (r) (r) − 2 (r) (r) − (r) 2 (r)
r δρ δρ δ ρ δρ δρ δ ρ
0 0

δu δp δu δp
+ 2πβ R (R) (R) + (R) (R)
δρ δρ δρ δρ
δ 2u δ2p

2 1 δu δp δp δu
+ 2πβ R (R) (R) + 2 (R) (R) + (R) 2 (R)
R δρ δρ δ ρ δρ δρ δ ρ
The shape derivative u0 of the state solves the system:
−∆ u0 = 0

in Ω 
δu
u0 = −(V · n) on δ Ω 
δn
The shape derviative u0 is decomposed to
u0 = αQ1 + β Q2 .
with
 
−∆ Q1 = 0 in Ω 
 −∆ Q2 = 0 in Ω 

 
δu Q2 = 0 on Γr
 
Q1 = − on Γr and
δρ  δu 
Q2 = − on ΓR . 

 
Q1 = 0 on ΓR
 
δρ
For this system we get the solutions

δu ln(R) − ln(ρ) δu ln(r) − ln(ρ)
Q1 (ρ) = − (r) and Q2 (ρ) = (R) .
δρ ln(R) − ln(r) δρ ln(R) − ln(r)
The shape derivative p0 of the adjoint solves
−∆ p0 = u0

in Ω 
δp
p0 = −(V · n) on δ Ω 
δn
The shape derviative p0 is decomposed to
p0 = α(Q3 + Q5 ) + β (Q4 + Q6 )
with
 
−∆ Q3 = 0 in Ω 
 −∆ Q4 = 0 in Ω 

 
δp Q4 = 0 on Γr
 
Q3 = − on Γr and
δρ  δp 
Q4 = − on ΓR . 

 

Q3 = 0 on ΓR

δρ
Solutions of this systems are

δp ln(R) − ln(ρ) δp ln(r) − ln(ρ)
Q3 (ρ) = − (r) and Q4 (ρ) = (R)
δρ ln(R) − ln(r) δρ ln(R) − ln(r)
The two further parts of the adjoint shape derivative solve

 
−∆ Q5 = Q1 in Ω   −∆ Q6 = Q2 in Ω 

Q5 = 0 on Γr Q6 = 0 on Γr
 
Q5 = 0 on ΓR Q6 = 0 on ΓR
 
with solutions
Q5 (ρ) = Q̄5 (ρ) − aQ̄5 ln(ρ) + bQ̄5 and Q6 (ρ) = Q̄6 (ρ) − aQ̄6 ln(ρ) + bQ̄6
where the particular solutions are
δ u ρ 2 (ln(R) − ln(ρ) + 1) δu ρ 2 (ln(r) − ln(ρ) + 1)

Q̄5 (ρ) = − (r) and Q̄6 (ρ) = (R) .
δρ 4(ln(R) − ln(r)) δρ 4(ln(R) − ln(r))
The shape hessian with the decompositions of the shape derivative reads as:

2 δ Q1 δ Q2 δp
∂ J(Ω ;V,V ) := − 2παr α (r) + β (r) (r)
δρ δρ δρ

δu δ Q3 δ Q5 δ Q4 δ Q6
− 2παr (r) α (r) + (r) + β (r) + (r)
δρ δρ δρ δρ δρ
δ 2u δ2p

1 δu δp δp δu
+ 2πα 2 r (r) (r) − 2 (r) (r) − (r) 2 (r)
r δρ δρ δ ρ δρ δρ δ ρ

δ Q1 δ Q2 δp
+ 2πβ R α (R) + β (R) (R)
δρ δρ δρ

δu δ Q3 δ Q5 δ Q4 δ Q6
+ 2πβ R (R) α (R) + (R) + β (R) + (R)
δ 2u δ2p

1 δu δp δp δu
+ 2πβ 2 R (R) (R) + 2 (R) (R) + (R) 2 (R)
R δρ δρ δ ρ δρ δρ δ ρ
It can be written as
2
2 α
δ J(Ω ,V,V ) = 2π α β ∇ F
β
with

(∇2 F)11 = r − δρ Q1 (r)δρ p(r) − δρ u(r) δρ Q3 (r) + δρ Q5 (r)

1
+ δρ u(r)δρ p(r) − δρ2 u(r)δρ p(r) − δρ u(r)δρ2 p(r)
r
2
(∇ F)12 = r − δρ Q2 (r)δρ p(r) − δρ u(r) δρ Q4 (r) + δρ Q6 (r)

(∇2 F)21 = R δρ Q1 (R)δρ p(R) + δρ u(R) δρ Q3 (R) + δρ Q5 (R)

(∇2 F)22 = R δρ Q2 (R)δρ p(R) + δρ u(R) δρ Q4 (R) + δρ Q6 (R)

1
+ δρ u(R)δρ p(R) + δρ2 u(R)δρ p(R) + δρ u(R)δρ2 p(R)
R
first order derivatives:

δu δΨ a δΨ 1
= − =− ρ
δρ δρ ρ δρ 2
δp δ p̃ δ η δ p̃ δ P̃ aP̃
= + = −
δρ δρ δρ δρ δρ ρ
δρ p∗ on Ω ∗

δ P̃  ∗ δ p∗ 1 1
= r δρ p∗ (r∗ ) ρ1 on {ρ < r∗ } = − ρ 3 + ρ(a∗ ( 21 − ln(ρ)) + b∗ )
δρ   ∗ δρ 16 2
R δρ p∗ (R∗ ) ρ1 on {ρ > R∗ }
δη δ N aN δN 1 1
= − = ρ 3 − ρ(a( 12 − ln(ρ)) + b)
δρ δρ ρ δρ 16 2
second order derivatives:

δ 2u δ 2Ψ a δ 2Ψ 1
2
= 2
+ 2 2
=−
δ ρ δ ρ ρ δ ρ 2
δ2p δ 2 p̃ δ 2 η δ 2 p̃ δ 2 P̃ aP̃
= 2 + = 2 + 2
δ 2ρ δ ρ δρ δ 2ρ δ ρ ρ
 2 ∗
δρ p on Ω ∗
δ 2 P̃  ∗ δ 2 p∗ 3 1
= −r δρ p∗ (r∗ ) ρ 2 on {ρ < r∗ }
1
= − ρ 2 + (a∗ (− 21 − ln(ρ)) + b∗ )
δ 2ρ   ∗ δ 2ρ 16 2
−R δρ p∗ (R∗ ) ρ12 on {ρ > R∗ }
δ 2η δ 2 N aN δ 2N 3 1
= + = ρ 2 − (a(− 12 − ln(ρ)) + b)
δ 2ρ δ 2ρ ρ 2 δ 2ρ 16 2
shape derivatives:
δ u0 δ Q1 δ Q2
=α +β
δρ δρ δρ
δ Q1 δu 1 1 δ Q2 δu 1 1
= (r) = − (R)
δρ δρ ρ ln(R) − ln(r) δρ δρ ρ ln(R) − ln(r)
δ p0

δ Q3 δ Q5 δ Q4 δ Q6
=α + +β +
δ Q3 δp 1 1 δ Q4 δp 1 1
= (r) = − (R)
δρ δρ ρ ln(R) − ln(r) δρ δρ ρ ln(R) − ln(r)
δ Q5 δ Q̄5 aQ̄5 δ Q6 δ Q̄6 aQ̄6
= − = −
δρ δρ ρ δρ δρ ρ
δ Q̄5 δ u ρ 2 ln(R) − 2 ln(ρ) + 1 δ Q̄6 δu ρ 2 ln(r) − 2 ln(ρ) + 1
= (r) = − (R)
δρ δρ 4 ln(R) − ln(r) δρ δρ 4 ln(R) − ln(r)
6.4.5 Matlab program
This Matlab program computes the optimal shape beginning with the ring with inner radius r and outer
radius R either with gradient descent method (typ=’g’) or Newton method (typ=’h’).
1 function ret = SO2D_Hessian(r,R,typ)

2
3 % optimal radii
4 rs = 1;
5 Rs = 2;
6
7 % initial
8 J = inf;
9 or = r;
10 oR = R;
11 oJ = J;
12 tau = 0.1;
13 gtyp = typ;
14 iter = 1;
15
16 while abs(J) >=eps && iter < 1000

17 % particular solution
18 PSI = @(x) -1/4*x.ˆ2;
19 dPSI = @(x) -1/2*x;
20 ddPSI = @(x) -1/2;
21
22 % fourier coefficients
23 a = @(psi,r,R) (psi(R)-psi(r))/(log(R)-log(r));
24 b = @(psi,r,R) (psi(R)*log(r)-psi(r)*log(R))/(log(R)-log(r));
25
26 % primal
27 u = @(x) PSI(x) - a(PSI,r,R)*log(x) + b(PSI,r,R);
28 du = @(x) dPSI(x) - a(PSI,r,R)./x;
29 ddu = @(x) ddPSI(x)+a(PSI,r,R)./x.ˆ2;
30
31 % target
32 us = @(x) PSI(x) - a(PSI,rs,Rs)*log(x) + b(PSI,rs,Rs);
33 zD = @(x) (x>=rs & x<= Rs).*us(x);
34
35 % adjoint
36 ps = @(x) -1/64*x.ˆ4 + 1/4*x.ˆ2.*(a(PSI,rs,Rs)*( 1-log(x)) + ...
b(PSI,rs,Rs));
37 dps = @(x) -1/16*x.ˆ3 + 1/2*x .*(a(PSI,rs,Rs)*( 1/2-log(x)) + ...
b(PSI,rs,Rs));
38 ddps = @(x) -3/16*x.ˆ2 + 1/2 .*(a(PSI,rs,Rs)*(-1/2-log(x)) + ...
b(PSI,rs,Rs));
39
40 pr = @(x) ps(rs) + rs*dps(rs)*(log(x)-log(rs));
41 dpr = @(x) rs*dps(rs)./x;
42 ddpr = @(x) -rs*dps(rs)./x.ˆ2;
43
44 pR = @(x) ps(Rs) + Rs*dps(Rs)*(log(x)-log(Rs));

45 dpR = @(x) Rs*dps(Rs)./x;
46 ddpR = @(x) -Rs*dps(Rs)./x.ˆ2;
47
48 PT = @(x) (x<rs).* pr(x) + (x>Rs).* pR(x) + (x>=rs & x<= Rs).* ps(x);
49 dPT = @(x) (x<rs).* dpr(x) + (x>Rs).* dpR(x) + (x>=rs & x<= Rs).* dps(x);
50 ddPT = @(x) (x<rs).*ddpr(x) + (x>Rs).*ddpR(x) + (x>=rs & x<= Rs).*ddps(x);
51
52 pt = @(x) PT(x) - a(PT,r,R)*log(x) + b(PT,r,R);

53 dpt = @(x) dPT(x) - a(PT,r,R)./x;
54 ddpt = @(x) ddPT(x) + a(PT,r,R)./x.ˆ2;
55
56 N = @(x) 1/64*x.ˆ4 - 1/4*x.ˆ2.*(a(PSI,r,R)*( 1-log(x)) + b(PSI,r,R));

57 dN = @(x) 1/16*x.ˆ3 - 1/2*x .*(a(PSI,r,R)*( 1/2-log(x)) + b(PSI,r,R));
58 ddN = @(x) 3/16*x.ˆ2 - 1/2 .*(a(PSI,r,R)*(-1/2-log(x)) + b(PSI,r,R));
59
60 n = @(x) N(x) - a(N,r,R)*log(x) + b(N,r,R);
61 dn = @(x) dN(x) - a(N,r,R)./x;

62 ddn = @(x) ddN(x) + a(N,r,R)./x.ˆ2;
63
64 p = @(x) pt(x) + n(x);

65 dp = @(x) dpt(x) + dn(x);
66 ddp = @(x) ddpt(x) + ddn(x);
67
68 % shape derivative: state
69 Q1 = @(x) -du(r)*(log(R)-log(x))/(log(R)-log(r));
70 dQ1 = @(x) du(r)*1./x*1/(log(R)-log(r));
71 Q2 = @(x) du(R)*(log(r)-log(x))/(log(R)-log(r));
72 dQ2 = @(x) -du(R)*1./x*1/(log(R)-log(r));
73
74 % shape derivative: primal

75 Q3 = @(x) -dp(r)*(log(R)-log(x))/(log(R)-log(r));
76 dQ3 = @(x) dp(r)*1./x*1/(log(R)-log(r));
77
78 Q4 = @(x) dp(R)*(log(r)-log(x))/(log(R)-log(r));
79 dQ4 = @(x) -dp(R)*1./x*1/(log(R)-log(r));
80
81 Q5b = @(x) du(r)*x.ˆ2*(log(R)-log(x)+1)/(4*(log(R)-log(r)));

82 dQ5b = @(x) du(r)*x*(2*log(R)-2*log(x)+1)/(4*(log(R)-log(r)));
83
84 Q6b = @(x) -du(R)*x.ˆ2*(log(r)-log(x)+1)/(4*(log(R)-log(r)));

85 dQ6b = @(x) -du(R)*x*(2*log(r)-2*log(x)+1)/(4*(log(R)-log(r)));
86
87 Q5 = @(x) Q5b(x) - a(Q5b,r,R)*log(x) + b(Q5b,r,R);
88 dQ5 = @(x) dQ5b(x) - a(Q5b,r,R)./x;
89
90 Q6 = @(x) Q6b(x) - a(Q6b,r,R)*log(x) + b(Q6b,r,R);
91 dQ6 = @(x) dQ6b(x) - a(Q6b,r,R)./x;
92
93 % hessian
94 H(1,1) = r*(-dQ1(r)*dp(r) - du(r)*(dQ3(r)+dQ5(r)) + 1/r*du(r)*dp(r) - ...
ddu(r)*dp(r) - du(r)*ddp(r));
95 H(1,2) = r*(-dQ2(r)*dp(r) - du(r)*(dQ4(r)+dQ6(r)));
96 H(2,1) = R*(dQ1(R)*dp(R) + du(R)*(dQ3(R)+dQ5(R)));
97 H(2,2) = R*(dQ2(R)*dp(R) + du(R)*(dQ4(R)+dQ6(R)) + 1/R*du(R)*dp(R) - ...
ddu(R)*dp(R) - du(R)*ddp(R));
98 H = 2*pi*H;
99
100 % gradient
101 dJ = 2*pi*[-r*du(r).*dp(r);R*du(R).*dp(R)];
102
103 % objective functional

104 rl = linspace(r,R,101);
105 rc = (rl(1:end-1)+rl(2:end))/2;
106 rh = diff(rl);
107 J = 2*pi*(u(rc)-zD(rc)).ˆ2.*rc*rh';
108
109 % find direction

110 switch lower(typ)
111 case {'hessian','h'}
112 d = -1*H\dJ;
113 tau = 1;
114 fprintf('%s','hessian: ')
115 case {'gradient','g'}
116 d = -dJ;
117 fprintf('%s','gradient: ')
118 end
119
120 fprintf('%3d % e % e %e %e\n',iter,r,R,J,tau)

121
122 ret.r(iter) = r;
123 ret.R(iter) = R;
124 ret.J(iter) = J;
125
126 r = max(r + tau*d(1),eps);

127 R = min(R + tau*d(2),4);
128 iter = iter+1;
129
130 % stepsize control
131 if( iter > 2)
132 if (J <= oJ + 1e-14 && r < R && J>=0 ) || tau < 1e-12
133 oJ = J;
134 or = r;
135 oR = R;
136 tau = tau*1.1;
137 else
138 r = or;
139 R = oR;
140 iter = iter-1;
141 tau = tau/2;
142 end
143 end
144 end
145 ret.typ = gtyp;
146 end
References
1. Sverak, V. On optimal shape design. J. Math. Pures Appl. (9) 72 (1993), no. 6, 537–551.
2. Zhu, Shengfeng ; Gao, Zhiming . Convergence analysis of mixed finite element approximations to shape gradients in the
Stokes equation. Comput. Methods Appl. Mech. Engrg. 343 (2019), 127–150.
3. Zhu, Shengfeng . Effective shape optimization of Laplace eigenvalue problems using domain expressions of Eulerian
derivatives. J. Optim. Theory Appl. 176 (2018), no. 1, 17–34.
4. Ganghoffer, J. F. ; Plotnikov, P. I. ; Sokolowski, J. Nonconvex model of material growth: mathematical theory. Arch. Ration.
Mech. Anal. 230 (2018), no. 3, 839–910.
5. Kazmierczak, Anna ; Sokolowski, Jan ; Zochowski, Antoni . Drag minimization for the obstacle in compressible flow using
shape derivatives and finite volumes. Math. Control Relat. Fields 8 (2018), no. 1, 89–115.
6. Ganghoffer, J. F. ; Goda, I. ; Novotny, A. A. ; Rahouadj, R. ; Sokolowski, J. Homogenized couple stress model of optimal
auxetic microstructures computed by topology optimization. ZAMM Z. Angew. Math. Mech. 98 (2018), no. 5, 696–717.
7. Nowak, Michal ; Sokolowski, Jan ; Zochowski, Antoni . Justification of a certain algorithm for shape optimization in 3D
elasticity. Struct. Multidiscip. Optim. 57 (2018), no. 2, 721–734.
8. S. Nazarov, J. Sokolowski, The topological derivative of the Dirichlet integral under the formation of a thin bridge. (Rus-
sian) ; translated from Sibirsk. Mat. Zh. 45 (2004), no. 2, 410–426, Siberian Math. J. Vol. 45 (2004), no. 2, 341–355.
9. G. Fremiot, W. Horn,A. Laurain,M. Rao, J. Sokolowski, On the analysis of boundary value problems in nonsmooth
domains. Dissertationes Math. Vol. 462 (2009), 149 pp.
10. W. Kosinski, Field singularities and wave analysis in continuum mechanics. Translated from the Polish by Kosiński and
M. Doliński. Translation edited by R. W. Ogden. Ellis Horwood Series: Mathematics and its Applications. Ellis Horwood
Ltd., Chichester; PWN—Polish Scientific Publishers, Warsaw, 1986. 251 pp. ISBN: 0-85312-543-0
11. R. A. Adams Sobolev spaces Academic press, New-York, 1975.
12. D. Adams, On the existence of capacitary strong type estimates in Rn , Arkiv for Matematik, 14 (1976), 125–140.
13. D. Adams and L. Hedberg, Function Spaces and Potential Theory, Springer-Verlag, Berlin, 1996, 366 pp.
14. C. Agmon, A. Douglis, L. Nirenberg, Estimates near the booundary for solutions of elliptic partial differential equations
satisfying general boundary conditions. I. Comm. Pure Apll. Math. 12 (1959), 623–727.
15. L. Ambrosio, Transport equations and Cauchy problem for non-smooth vector fields, Calculus of variations and nonlinear
partial differential equations, Lecture Notes in Math., 1927(2008), Springer, Berlin, 1–41.
16. H. Amann, Linear and quasilinear parabolic problems, Birkhauser Basel Inc., Boston, MA, 1995, 335 pp.
17. A. A. Amosov, A. A. Zlotnik Properties in the large of quasi-averaged equations of the one dimensional motion of a
viscous heat conacting gas, Dokl. Akad. Nauk, 346 (1996), no.2, 151–154.
18. A. A. Amosov, A. A. Zlotnik On quasi-averaged equations of the one-dimensional motion of a viscous barotropic medium
with rapidly oscillating data, Comput. Math. Math. Phys., 36 (1996), no.2, 203–220.
19. S. N. Antontsev, A. V. Kazhikhov, V. N. Monakhov, Boundary value problems of the mechanics of inhomogeneous fluids,
Nauka Sibirsk. Otdel., Novosibirsk, 1983. 320 pp; (English translation: North Holland)
20. D. G. Aronson. Non-negative solutions of linear parabolic equations Annali della Scuola Norm. Sup. Piza 22 (1968),
607–694.
21. S. Axler, P. Bourdon, W. Ramey, Harmonic function theory, second edition, Springer-Verlag, New-York, 2001.
22. N. S. Bakhvalov M. E. Eglit Processes in periodic media not describle in terms of mean characteristics Dokl. Akad.
Nauk SSSR, 268(1983), no.4, 836–840.
23. J. M. Ball, A version of the fundamental theorem for Young measures, in PDEs and Continuum Models of Phase Transi-
tions, Lecture Notes in Phys. 344, Springer, Berlin, (1989), 207–215.
24. V. Barbu, I. Lasiecka, R. Triggiani, Tangential boundary stabilization of Navier-Stokes equations, Memoirs of the AMS,
181, Providence, 2006.
25. H. Beirão da Veiga, Stationary motions and the incompressible limit for compressible viscous limit Houston J. Math., 13
(1987), 527–544.
26. H. Beirão da Veiga An L p - theory for n-dimensional stationary compressible Navier-Stokes equations, and the incom-
pressible limit for compressible fluids. The equilibrium solutions Comm. Math. Phys., 109 (1987), 229–248.
27. H. Beirão da Veiga Existence results in Sobolev spaces for a transport equation Ricerche Mat. 36 (1987), suppl., 173–184.
28. J. A. Bello, E. Fernandez-Cara, J. Lemoine, J. Simon, The differentiability of the drag with respect to variations of a
Lipschitz domain in a Navier-Stokes flow, SIAM J. Control. Optim. 35 (1997), 626–640.
29. J. A. Bello, E. Fernandez-Cara, J. Simon, Optimal shape design for Navier-Stokes flow, in System Modelling and Opti-
mization, Lecture Notes in Control and Inform Sci 180, D. Kall, ed., Springer-Verlag, Berlin, 1998.
30. J. Bergh, J. Löfström, Interpolation spaces. An Introduction, Springer-Verlag, Berlin Heidelberg New-York, 1976.
86 Bibliography
31. Y. Brenier, Resolution d’equations quasiliniaires en dimension N d’espace al’aide d”equations lineaires en dimension
N + 1, Differential equations 50 (1983), 375–390.
32. R. Courant, D. Hilbert, Methods of mathematical physics. Vol. II: Partial differential equations. (Vol. II by R. Courant.)
Interscience Publishers (a division of John Wiley & Sons), New York-London, 1962, 830 pp.
33. M. C. Delfour, J.-P. Zolesio, Shapes and Geometries, Advances in Design and Control. Society for Industrial and Applied
Mathematics (SIAM), Philadelphia, PA, 2001.
34. B. Desjardins, E. Grenier, P.-L. Lions, N. Masmoudi Incompressible limitit for solutions of the isentropic Navier-Stokes
Equationa with Dirichlet boundary data J. Math. Pures. Appl., 78 (1999), 461–471.
35. R. J. DiPerna, P.-L. Lions, Ordinary differential equations, transport theory and Sobolev spaces Invent. Math. 48(1989),
511–547.
36. H. Ebert, Physicalisches Tascenbuch, Friedr, Vieweg & Sohn, Braunschweig, 1957.
37. R. E. Edwards, Functional Analysis. Theory and Applications, Holt, Rinehart, and Wilson, New York, Toronto, London,
1965.
38. L. C. Evans, Partial Differential Equations Am. Math. Soc., 1998.
39. R. Farwig, Stationary solutions to compressible Navier -Stokes equations Comm. Part. Diff. Eq., 14 (1989), 1579–1506.
40. R. Farwig, G. Galdi, H. Sohr A new class of weak solutions of the Navier-Stokes Equations, J. Math. Fluid Mech. 8 (2006),
423–443.
41. G. Fichera Sulle equazioni differenziali lineari elliptico-paraboliche del secondo ordine Atti Accad. Naz. Lincei, Mem.
Cl. Sci. Fis. Mat. Nat. Sez. I. (8)5 (1956) 1–30.
42. E. Feireisl Dynamics of Viscous Compressible Fluids Oxford University Press, Oxford, 2004.
43. E. Feireisl, A. H. Novotný, Singular limits in thermodynamics of viscous fluids, Birkhäuser, Basel, 2009.
44. E. Feireisl, Asymptotic analysis of the full Navier-Stokes-Fourier systems: From compressible to incompressible fluid
flows Uspekhi Mat. Nauk, 62 (2007), 169–192.
45. E. Feireisl, On compactness of solutions to the compressible isentropic Navier-Stokes equations when the density is not
square integrable Comment. Math. Univ. Carolinae 42 (2001), 83–98.
46. E. Feireisl, A. H. Novotný, H. Petzeltová, On the existence of globally defined weak solutions to the Navier-Stokes equa-
tions J. of Math. Fluid Mech. 3 (2001), 358–392.
47. E. Feireisl Shape optimization in viscous compressible fluids Appl. Math. Optim. 47 (2003), 59–78.
48. E. Feireisl, A. H. Novotný, H. Petzeltová, On the domain dependence of solutions to the compressible Navier–Stokes
equations of a barotropic fluid, Math. Methods Appl. Sci., 25 (2002), 1045–1073.
49. W. Feller, An introduction to the Probability Theory and its Applications, 3rd edition, Wiley, New-York, 1968.
50. I. Fonseca, G. Leoni, Modern methods in the calculus of variations: L p spaces, Springer, New York, 2007. 599 pp.
51. J. Frehse, S. Goj, M. Steinhauer, L p -estimates for the Navier-Stokes equations for steady compressible flow, Manuscripta
Math. 116 (2005), 265–275.
52. A. Fursikov, Optimal control of distributed systems. Theory and applications Translation of Mathematical Monographs
187, AMS, Providence, Rhode Island, 2000, 305 pp.
53. G. Galdi, An introduction to the mathematical theory of the Navier-Stokes equations. Vol. I. Linearized steady problems,
Springer-Verlag, Berlin, Heidelberg New-York, 1994, 450 pp.
54. D. Gilbarg, N. S. Trudinger Elliptic Partial Differential Equations of Second Order, Springer-Verlag, 1983, 513 pp.
55. V. Girinon Navier-Stokes equations with nonohomogeneous boundary conditions in a bounded three-dimensional domain,
J. Math. Fluid Mech. (2010) DOI 10.1007/s00021-009-0018x
56. P. Grisward Caracterisation de Quelques Espaces, d’ Interpolation, Arch. Rat. Mech. Anal. 25 (1967), 43–63.
57. L. Hedberg, Spectral synthesis in Sobolev spaces and uniqueness of solutions of the Dirichlet problem, Acta Math., 147
(1981), 237–264.
58. J. G. Heywood, M. Padula, On the uniqueness and existence theory for steady compressible viscous flow, Fundamental
directions in mathematical fluids mechanics, Adv. Math. Fluids Mech., Birkhauser, Basel (2000), 171–189.
59. L. Hörmander Non-elliptic boundary value problems Ann. of Math., 83 (1966), 129–209.
60. S. Ishikawa, Fixed point by new iteration method, Proc. Amer. Math. Soc., (1)44 (1974), 147–150.
61. W. Jager, A. Mikelic Couette flow over a rough boundary and drag reduction, Comm. Math. Phys. 232 (2003), 429–455.
62. H. Kang, H. Koo, Estimates of harmonic Bergman kernel on smooth domains, Functional analysis, 185 (2001), 220–239.
63. B. Kawohl, O. Pironneau, L. Tartar, J. Zolesio, Optimal Shape Design, Lecture Notes in Math. 1740, Springer-Verlag,
Berlin, 2000.
64. R. B. Kellogg Discontinuous solutions of the linearized, steady state, compressible viscous Navier-Stokes equations,
SIAM J. Math. Anal.,19 (1988), 567–579.
65. Jae Ryong Kweon, R. B. Kellogg Compressible Navier-Stokes equations in a bounded domain with inflow boundary
condition, SIAM J. Math. Anal.,28 (1997), 94–108.
66. Jae Ryong Kweon, R. B. Kellogg Regularity of solutions to the Navier-Stokes equations for compressible barotropic flows
on a polygon, Arch. Ration. Mech. Anal., 163 (2000), 36–64.
67. Jae Ryong Kweon, R. B. Kellogg Singularity in the density gradient J. Math. Fluid Mech. 8 (2006), 445–454.
68. J. R. Kweon, R. B. Kellogg, Regularity of solutions to the Navier–Stokes system for compressible flows on a polygon,
SIAM J. Math. Anal., 35 (2004), 1451–1485.
69. J. R. Kweon, R. B. Kellogg, Compressible Stokes problem on nonconvex polygonal domains, J. Differential Equations,
176 (2001), 290–314.
70. J. J. Kohn, L. Nirenberg, Degenerate elliptic-parabolic equations of second order Comm. Pure and Appl. Math,
(4)20(1967), 797–872.
71. A. N. Kolmogorov, S. V. Fomin, Elements of the theory of functions and functional analysis, Mir, Moscow, 1989, 624 pp.
72. L. D. Landau, E. M. Lifshitz, Course of theoretical physics. Vol 6. Fluid mechanics, Pergamon Press, Oxford, 1987.
73. M. A. Krasnosel’skii, Topological methods in the theory of nonlinear integral equations, The Macmillan Co., New York,
1964, 395 pp.
74. O. A. Ladyzenskaja, V. A. Solonnikov, N. N. Ural’ceva, Linear and quasiliniear equation of parabolic type, Americam
Math. Soc., Providence, RI, 1967, 648 pp.
Biblography 87
75. N. S. Landkof, Foundations of modern potential theory, Springer-Verlag, Berlin, Heidelberg, New-York, 1972, 424 pp.
76. C. De Lellis,Notes on hyperbolic systems of conservation laws and transport equations, Handbook of differential equa-
tions: evolutionary equations. Vol. III, 277–382, Hand. Differ. Equ., Elsevier/North-Holland, Amsterdam, 2007.
77. J.-L. Lions, E. Magenes, Problèmes aux limites non homegènes et applications, Vols 1,2, Dunod, Paris,1968, 372 pp.,251
pp.
78. P.-L. Lions, B. Perthame, E. Tudmor, A kinetic formulation of multidimensional conservation laws and related equations,
J. Amer. Math. Soc., 7 (1994), 415–431.
79. P. L. Lions Mathematical topics in fluid dynamics, Vol. 2, Compressible models The Clarendon Press, Oxford University
Press, New York, 1998, 348 pp.
80. P. L. Lions Mathematical topics in fluid dynamics, Vol. 1, Compressible models The Clarendon Press, Oxford University
Press, New York, 1996, 237 pp.
81. P. L. Lions, B. Perthame, P. Souganidis, Existence and stability of entropy solutions for the hyperbolic systems of isentropic
has dynamics in Eulerian and Lagrangian coordinates, Comm. Pure Apll. Math., 59 (1996), 599–638.
82. P. L. Lions, Compacité des solutions des équations de Navier-Stokes compressibles isentropiques, C.R. Acad. Sci. Paris,
Sér. I Math., 317 (1993), 115–120.
83. P. L. Lions, Bornes sur la densité pour les équations de Navier-Stokes compressibles isentropiques avec conditions aux
limites de Dirichlet, C.R. Acad. Sci. Paris, Sér. I Math., 328 (1999), 659–662.
84. P. L. Lions, On some challenging problems in nonlinear partial differential equations, in Mathematics: Frontiers and
Perspectives, V. Arnold et al., Eds., AMS Providence, RI, 2000, 121–135.
85. J. Löfström, Interpolation of boundary value problems of Neumann type on smooth domains J. London Math. Soc. 46
(1992), 499–516.
86. J. Löfström Interpolation of subspaces, Techn. Rep. 10, University of Göterborg (1997).
87. E. Magenes, Interpolational spaces and partial differential equations Uspehi Mat. Nauk, 21 (1966), 169–218.
88. V. G. Maz’ya Sobolev spaces, Springer-Verlag, Berlin, 1985, 486 pp.
89. V. G. Maz’ya, T. O. Shaposhnikova, Theory of Sobolev multipliers. With applications to differential and integral operators,
Springer-Verlag, Berlin, 2009, 609 pp.
90. N. Meyer, An L p -estimate for the gradient of solutions of second order elliptic divergence equations Ann. Scuola Nor.
Super. di Pisa 17 (1963), 189–207.
91. Marwan Moubachir, J.-P. Zolesio Moving Shape Analysis And Control: Applications to Fluid Structure Interactions Chap-
man & Hall/CRC, Boca Raton, 2006.
92. B. Mohammadi, O. Pironneau, Shape optimization in fluid mechanics Ann. Rev. Fluid Mech., 36 (2004), 255–279.
93. A. A. Novotny, J. Sokolowski, Topological derivatives in shape optimization. Interaction of Mechanics and Mathematics.
Springer, Heidelberg, 2013. xxii+412 pp. ISBN: 978-3-642-35244-7; 978-3-642-35245-4
94. A. A. Novotny, J. Sokolowski, A. Zochowski, Applications of the Topological Derivative Method. Studies in Systems,
Decision and Control 188, Springer, 2019. xiv, 212 pp. ISBN: 978-3-030-05431-1;
95. S. Novo Compressible Navier-Stokes Model with Inflow-Outflow Boundary Conditions J. Math. Fluid Mech. 7 (2005),
485-514.
96. S. Novo, A. Novotný, On the existence of weak solutions to the steady compressible Navier-Stokes equations when the
density is not square integrable, J. Math. Kyoto Univ., 42 (2002), 531–550.
97. S. Novo, A. Novotný, On the existence of weak solutions to the steady compressible Navier-Stokes equations in domains
with conical outlets, J. Math. Fluid Mech., 8 (2006), 187–210.
98. A. Novotný, I. Straskraba, Introduction to the mathematical theory of compressible flow Oxford Lecture Series in Math-
ematics and its Applications, 27 Oxford University Press, Oxford, 2004.
99. A. Novotný, M. Padula, Existence and Uniqueness of Stationary solutions for viscous compressible heat conductive fluid
with large potential and small non-potential external forces, Siberian Math. Journal, 34 (1993), 120-146
100. A. Novotný, M. Padula, L p -Approach to Steady flows of Viscous Compressible Fluids in Exterior Domains Arch. Rat.
Mech. Anal, 126 (1994), 243–297.
101. A. Novotný, M. Padula, Physically reasonable solutions to steady compressible Navier-Stokes equations in 3D-Exterior
Domains Math. Ann., 308 (1997), 438–439.
102. A. Novotný, About steady transport equation. I. L p -approach in domains with smooth boundaries Comment. Math. Univ.
Carolin., 37 (1996), 43–89.
103. A. Novotný, About steady transport equation. II. Schauder estimates in domains with smooth boundaries Portugal.
Math., 54 (1997), 317–333.
104. O. A. Oleinik, E. V. Radkevich, Second order equation with non-negative characteristic form American Math. Soc.,
Providence, Rhode Island Plenum Press. New York-London, 1973.
105. M. Padula, Existence and uniqueness for viscous steady compressible motions Arch. Rational Mech. Anal., 97 (1986),
1–20.
106. M. Padula, Steady flows of barotropic viscous fluids, Classical Problems of in Mechanics, 253–345, Quad. Math. 1, Dept.
Math., Seconda Univ. Napoli, Caserta, 1997.
107. P. Pedregal, Parametrized measures and variational principles, Progr. Nonlinear Differential Equations Appl. 30,
Birkhäuser Verlag, Basel, 1997.
108. P. I. Plotnikov, J. Sokolowski, Compressible Navier-Stokes equations. Theory and shape optimization, Instytut Matem-
atyczny Polskiej Akademii Nauk. Monografie Matematyczne (New Series) [Mathematics Institute of the Polish Academy
of Sciences. Mathematical Monographs (New Series)], 73. Birkhäuser/Springer Basel AG, Basel, 2012. xvi+457 pp. ISBN:
978-3-0348-0366-3
109. P. I. Plotnikov, J. Sokolowski, Stationary Boundary Value Problems for Navier-Stokes Equations with Adiabatic Index
ν < 3/2, Doklady Mathematics 70 (2004), 535–538.
110. P. I. Plotnikov, J. Sokolowski, On compactness, domain dependence and existence of steady state solutions to compress-
ible isothermal Navier-Stokes equations J. Math. Fluid Mech. 7 (2005), 529–573.
111. P. I. Plotnikov, J. Sokolowski, Concentrations of solutions to time-discretizied compressible Navier -Stokes equations
Communications in Mathematical Phisics 258 (2005), 567–608.
88 Bibliography
112. P. I. Plotnikov, J. Sokolowski, Domain dependence of solutions to compressible Navier-Stokes equations SIAM J. Con-
trol Optim., 45 (2006), 1147–1539.
113. P. I. Plotnikov, J. Sokolowski, Stationary solutions of Navier-Stokes equations for diatomic gases Russian Mathematical
Surveys, 62 (2007), 561–593.
114. P. I. Plotnikov, J. Sokolowski, Stationary boundary value problems for compressible Navier-Stokes equations Hand-
book of differential equations: stationary partial differential equations. Vol. VI, 313–410, Hand. Diff. Equ., Elsevier/North
Holland, Amsterdam, 2008.
115. P. I. Plotnikov, E. V. Ruban, J. Sokolowski, Inhomogeneous boundary value problems for compressible Navier-Stokes
equations, well-posedness and sensitivity analysis, SIAM J. Math. Analysis 40 (2008), 1152–1200.
116. P. I. Plotnikov, E. V. Ruban, J. Sokolowski, Inhomogeneous boundary value problems for compressible Navier-Stokes
and transport equations, Journal des Mathématiques Pure et Appliquées, 92(2009), 113–162.
117. P. I. Plotnikov, J. Sokolowski, Shape derivative of drag functional, SIAM Journal on Control and Optimization,
(7)48(2010), 4680–4706.
118. P. I. Plotnikov, J. Sokolowski, Inhomogeneous boundary value problem for nonstationary compressible Navier-Stokes
equations, Journal of Mathematical Sciences, (1)170 (2010), 34–130.
119. W. Ramey, H. Yi, Harmonic Bergman functions on half-spaces Trans. Amer. Math. Soc., 348 (1996), 633–660.
120. J.-F. Scheid, J. Sokolowski, Shape optimization for a fluid-elasticity system. Pure Appl. Funct. Anal. 3 (2018), no. 1,
193–217.
121. H. Schlichting, Boundary-layer theory, (McGraw-Hill series in mechanical engineering) New York: McGraw-Hill, 1955.
122. J. Serrin, Mathematical Principles of Classical Fluid Mechanics, Hanbook der Physic VIII/I, Springer, Berlin, 1972.
123. J. Simon, Domain variation for drag in Stokes flow, in Control and Estimation of Distributed Parameter Systems,
Internat. N Ser. Numer. Math. 91, F. Kappel, K. Kuninisch, and W. Schappacher, Eds., Birkhäuser, Basel, 1989, 361–378.
124. T. Slawig ,A formula for the derivative with respect to domain variations in Navier-Stokes flow based on an embedding
domain method SIAM J. Control Optim., 42 (2003), 495–512.
125. T. Slawig, An explicit formula for the derivative of a class of cost functionals with respect to domain variations in Stokes
flow SIAM J. Control Optim., 39 (2000), 141–158.
126. Sokolowski, J.; Zochowski, A. On the topological derivative in shape optimization. SIAM J. Control Optim. 37 (1999),
no. 4, 1251-1272 (electronic).
127. J. Sokolowski, J.-P. Zolésio, Introduction to Shape Optimization. Shape Sensitivity Analysis. Springer Series in Compu-
tational Mathematics 16, Springer Verlag, 1992.
128. E. Stein, Singular integrals and differentiability properties of functions, Princeton Unicersity Press, Princeton, 1970,
290 pp.
129. K. Stroethoff, Harmonic Bergman spaces Holomorphic spaces, 33 (1998), 51–63.
130. L. Tartar, An introduction to Sobolev spaces and interpolation spaces Lecture notes of the Unione Matematica Italiana
3, Springer, Berlin, UMI, Bologna, 2007, 218 pp.
131. Xinfu Chen, Weiqing Xie, Discontinuous solutions of steady state viscous compressible Navier-Stokes equations, J. Diff.
Equations, 115 (1995), 567–579.
132. Zhang Qi S. Gaussian Bounds for the Fundamental Solutions Manuscripta math. 93 (1997), 381–390.
133. W. P. Ziemer, Weakly Differentiable Functions, Springer-Verlag, New York, 1989.
134. Cea, J.: Problems of shape optimal design, in: Optimization of Distributed Parameter Structures, Haug E.J., Cea J. Eds.,
Sijthoff & Noordhoff 1981, pp. 1005-1049.
135. Chenais, D.: On the existence of a solution in a domain identification problem, J. Mathematical Analysis and Applica-
tions, vol. 52, no. 2, 1975.
136. Cowin S.C., Hegedus D.H., Bone remodeling I: A theory of adaptive elasticity, J. Elasticity, 6, 313–326, 1976.
137. Huiskes R. et al.,Effects of mechanical forces on maintenance and adaptation of form in trabecular bone, Nature 404,
704-706, 2000.
138. Jang I., Kim I., Computational study of wolff’s law with trabecular architecture in the human proximal femur using
topology optimization, Journal of Biomechanics, 41(11), pp. 2353–2361, 2008.
139. Klarbring A., Torstenfelt B., Lazy zone bone remodelling theory and its relation to topology optimization, Ann. Solid
Struct. Mech. 4(1), 25–32, 2012.
140. Nazarov S.A., Plamenevsky, B.A.: Elliptic Problems in Domains with Piecewise Smooth Boundaries, de Gruyter Expo-
sitions in Mathematics 13, 1994.
141. Nowak M.,Structural optimization system based on trabecular bone surface adaptation, Journal of Structural and Mul-
tidisciplinary Optimization, Springer Berlin Heidelberg, Volume 32, Number 3 , pp. 241–251, 2006.
142. Nowak M., On Some Properties of Bone Functional Adaptation Phenomenon Useful in Mechanical Design, Acta of
Bioengineering and Biomechanics, Vol. 12, No.2, pp. 49–54, 2010.
143. Pironneau, O.: Optimum design with Lagrangian Finite Elements: Design of an electromagnet, Computational Mathe-
matics in Applied Mechanics and Engineering, vol. 15, 1978, pp. 207-308.
144. Wasiutyński Z., On the congruency of the forming according to the minimum potential energy with that according to
equal strength, Bull. de l’Academie Polonaisedes Sciences, Serie des Sciences Techniques 8(6), pp. 259–268, 1960.
145. Wolff J., Das Gesetz der Transformation der Knochen, Hirschwald, 1892
146. Zolesio, Jean-Paul: The material derivative (or speed) method for shape optimization, in: Optimization of Distributed
Parameter Structures, Haug E.J., Cea J. Eds., Sijthoff & Noordhoff 1981, pp. 1089-1151.
147. P. Plotnikov J. Sokolowski and A. Zochowski, Numerical experiments in drag minimization for compressible Navier-
Stokes flows in bounded domains, 14th International IEEE/IFAC Conference on Methods and Models in Automation and
Robotics MMAR’09, 2009, 4 pages.
148. T. Gallouet R. Herbin and J. Latche, A convergent finite element-finite volume scheme for the compressible Stokes prob-
lem; part I - The isothermal case, to appear in print.
149. R. Eymard T. Gallouet R. Herbin and J. Latche, A convergent finite element-finite volume scheme for the compressible
Stokes problem; part II - The isentropic case, to appear in print.
Biblography 89
150. P. Plotnikov E. Ruban and J. Sokolowski, Inhomogenous boundary value problems for compressible Navier-Stokes equa-
tions: well-posedness and sensitivity analysis, SIAM Journal on Mathematical Analysis, 2008.
151. P.I. Plotnikov, J. Sokolowski Stationary Boundary Value Problems for Navier-Stokes Equations with Adiabatic Index
ν < 3/2, Doklady Mathematics Vol. 70, No. 1, 2004, 535-538.
152. P.I. Plotnikov, J. Sokolowski On compactness, domain dependence and existence of steady state solutions to compressible
isothermal Navier-Stokes equations J. Math. Fluid Mech. 7(2005), no. 4, 529-573.
153. P.I. Plotnikov, J. Sokolowski Concentrations of solutions to time-discretizied compressible Navier -Stokes equations
Communications in Mathematical Physics Volume 258, Number 3, 2005, 567-608.
154. P.I. Plotnikov, J. Sokolowski, Domain dependence of solutions to compressible Navier-Stokes equations , SIAM Journal
on Control and Optimization, 45 (2006), 1165 – 1197.
155. P.I. Plotnikov, J. Sokolowski, Stationary solutions of Navier-Stokes equations for diatomic gases., Russian Mathe-
matical Surveys, 62 (2007), 561– 593, Uspekhi Mat. Nauk 117-148.
156. P.I. Plotnikov, J. Sokolowski, Stationary Boundary Value Problems for Compressible Navier-Stokes Equations, Hand-
book of Differential Equations, Edited by M. Chipot 6 (2008), 313– 410.
157. P.I. Plotnikov, J. Sokolowski, Inhomogeneous boundary value problems for compressible Navier-Stokes equations, well-
posedness and sensitivity analysis, SIAM Journal on Mathematical Analysis 40 (2008), 1152– 1200.
158. P.I. Plotnikov, J. Sokolowski, Compressible Navier-Stokes equations, Discrete and Continuous Dynamical Systems,
Supplement 2009, 602-611.
159. P.I. Plotnikov, J. Sokolowski, Compressible Navier-Stokes equations SIAM J. Control Optim., Volume 48, Issue 7, pp.
4680-4706 (2010).
160. Arora, R.K., Optimization - algorithms and applications CRC Press, (2015)
161. Brampton C., Kim H., Cunningham J., Applications of 3D Level Set Topology Optimization, ASME 2012 Interna-
tional Design Engineering Technical Conferences and Computers and Information in Engineering Conference, paper No.
DETC2012-70870, pp. 875-882, (2012)
162. Bruggi M., On an alternative approach to stress constraints relaxation in topology optimization, Structural and Multidis-
ciplinary Optimization, Volume 36, Issue 2, pp 125–141, (2008)
163. Coelho P., Cardosoa J., Fernandes P., Rodrigues H., Parallel computing techniques applied to the simultaneous design of
structure and material, Advances in Engineering Software, Volume 42, Issue 5, pp. 219–227, (2011)
164. Cowin S.C., Hegedus D.H., Bone remodeling I: A theory of adaptive elasticity, J. Elasticity, 6, 313-326, (1976)
165. Delfour, M. C., Zolésio, J.-P., Shapes and geometries, SIAM, (2011)
166. Fernandes P., Rodrigues H, Jacobs C., A Model of Bone Adaptation Using a Global Optimisation Criterion Based on the
Trajectorial Theory of Wolff, Computer Methods in Biomechanics and Biomedical Engineering, vol. 2, pp. 125-138, (1999)
167. Florio C., Development of a widely applicable gradientless shape optimization based bone adaptation model for com-
parative parametric studies, Structural and Multidisciplinary Optimization, Volume 52, Issue 1, pp. 157–177 (2015)
168. Goda I., Ganghoffer J. F., Maurice, G., Combined bone internal and external remodeling based on Eshelby stress, Inter-
national Journal of Solids and Structures, 94, pp. 138-157, (2016)
169. Haftka, R., Gürdal, Z., Elements of structural optimization, 3rd edition, Kluwer, (1992)
170. Huiskes R., If bone is the answer, then what is the question?, Journal of Anatomy, 197, pp. 145–156, (2000)
171. Huiskes R. et al.,Effects of mechanical forces on maintenance and adaptation of form in trabecular bone, Nature 404,
704-706, (2000)
172. Jang I., Kim I., Computational study of Wolff’s law with trabecular architecture in the human proximal femur using
topology optimization, Journal of Biomechanics, 41(11), pp. 2353-2361, (2008)
173. Klarbring A., Torstenfelt B., Lazy zone bone remodelling theory and its relation to topology optimization, Ann. Solid
Struct. Mech. 4(1), 25-32, (2012)
174. Krog L., Tucker A., Kemp M., and Boyd R., Topology optimization of aircraft wing box ribs, AIAA Paper 2004-4481,
(2004)
175. Le C., Norato J., Bruns T., Ha C., Tortorelli D., Stress-based topology optimization for continua, Structural and Multi-
disciplinary Optimization, Volume 41, Issue 4, pp 605–620, (2010)
176. Lee Y. H., Kim Y., Kim J. J., Jang I. G., Homeostasis-based aging model for trabecular changes and its correlation with
age-matched bone mineral densities and radiographs, European journal of radiology, 84(11), pp. 2261-2268, (2015)
177. Makowski P., Kuś, Optimization of bone scaffold structures using experimental and numerical data, Acta Mechanica,
Volume 227, Issue 1, pp. 139–149, (2016)
178. Marzban A., Nayeb-Hashemi H., Vaziri A., Numerical simulation of load-induced bone structural remodelling using
stress-limit criterion, Computer methods in biomechanics and biomedical engineering, 18(3), pp. 259-268, (2015)
179. Nazarov S.A., Plamenevsky, B.A., Elliptic Problems in Domains with Piecewise Smooth Boundaries, de Gruyter Expo-
sitions in Mathematics 13, (1994)
180. Nowak M.,Structural optimization system based on trabecular bone surface adaptation, Journal of Structural and Multi-
disciplinary Optimization, Springer Berlin Heidelberg, Volume 32, Number 3 , pp. 241-251, (2006)
181. Nowak M., On Some Properties of Bone Functional Adaptation Phenomenon Useful in Mechanical Design, Acta of
Bioengineering and Biomechanics, Vol. 12, No.2, pp. 49-54, (2010)
182. Nutu E., Interpretation of parameters in strain energy density bone adaptation equation when applied to topology opti-
mization of inert structures, Mechanika, Vol. 21(6), pp. 443-449, (2015)
183. Pedersen P., Optimal Designs - Structures and Materials - Problems and Tools. ISBN 87- 90416-06-6, (2003)
184. Sigmund O., On the Optimality of Bone Microstructure. Synthesis in Bio Solid Mechanics, Kluwer, p. 221-234, (1999)
185. Shimoda M., Azegami H., Sakurai T., Traction Method Approach to Optimal Shape Design Problems, SAE Transactions,
Journal of Passenger Cars, 106 (6), pp. 2355-2365, (1998)
186. Shimoda M., Motora S., Ohtani H., Shape Optimization Method for Designing Interface Shapes of Composite Clad
Structures, ECCM15, Venice, Italy, 24-28 June 2012, (2012)
187. Plotnikov, P., Sokolowski, J., Compressible Navier-Stokes equations. Theory and Shape Optimization,
Birkhäuser/Springer Basel AG, Basel, (2012)
90 Bibliography
188. Rodrigues H., Jacobs C., Guedes M., Bendsoe M., Global and local material optimization applied to anisotropic bone
adaptation. In: Perdersen P, Bendsoe MP (eds) Synthesis in bio solid mechanics, Kluwer Academic Publishers, Dordrecht,
pp. 209-220, (1999)
189. Wasiutyński Z., On the congruency of the forming according to the minimum potential energy with that according to
equal strength, Bull. de l’Academie Polonaisedes Sciences, Serie des Sciences Techniques 8(6), pp. 259-268, (1960)
190. Wolff J., Das Gesetz der Transformation der Knochen, Hirschwald, (1892)
191. Wu J., Aage N., Westermann R., Sigmund O., Infill Optimization for Additive Manufacturing - Approaching Bone-like
Porous Structures, IEEE Transactions on Visualization and Computer Graphics, doi:10.1109/TVCG.2017.2655523, (2017)
192. V. Barbu, Optimal control of variational inequalities, Research Notes in Mathematics 100, Pitman (Advanced Publishing
Program), Boston, MA, 1984, iv+298 pages.
193. Avner Friedman, Variational principles and free-boundary problems, Second Edition, Robert E. Krieger Publishing Co.,
Inc., Malabar, FL, 1988, x+710 pages.
194. D. Kinderlehrer, G., Stampacchia, An introduction to variational inequalities and their applications, Classics in Applied
Mathematics 31, Reprint of the 1980 original, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA,
2000, xx+313 pages.
195. G. P. Galdi, M. Kyed, Steady flow of a Navier-Stokes liquid past an elastic body, Arch. Ration. Mech. Anal. 194 (2009),
no. 3, 849–875.
196. V.I. Gorbachuk, M.L. Gorbachuk, Boundary Value Problems for Operator Differential Equations, Kluwer Academic
Publ., Dordrecht, 1991.
197. C. Grandmont, Existence for a three-dimensional steady state fluid-structure interaction problem, J. Math. Fluid Mech.,
4 (2002), no. 1, 76–94.
198. A. Halanay, C.M. Murea, D. Tiba, Existence of a Steady Flow of Stokes Fluid Past a Linear Elastic Structure Using
Fictitious Domain, J. Math. Fluid Mech., 18 (2016), no. 2, 397–413.
199. A. Henrot, M. Pierre, Variation et optimisation de formes: Une analyse géométrique. Mathématiques et Applications,
vol. 48. Springer Berlin Heidelberg, 2005.
200. I. Lasiecka, K. Szulc, A. Żochowski, Linearised coupling of elasticity and Navier-Stokes equations, 20th International
Conference on Methods and Models in Automation and Robotics (MMAR) (2015), 208–210.
201. C. Surulescu, On the stationary motion of a Stokes fluid in a thick elastic tube: a 3D/3D interaction problem, Acta Math.
Univ. Comenianae, Vol. LXXV, 1 (2006), 95–106.
202. M. Hintermüller, V.A. Kovtunenko, From shape variation to topology changes in constrained minimization: a velocity
method based concept. Optim. Methods Softw. 26 (2011), 513–532.
203. A. M. Khludnev, J. Sokolowski, Modelling and control in solid mechanics, International Series of Numerical Mathemat-
ics, Vol. 122, Birkhäuser Verlag, Basel, ISBN 3-7643-5238-8,1997.
204. A.M. Khludnev, V.A. Kovtunenko, Analysis of Cracks in Solids. WIT-Press, Southampton, Boston 2000.
205. A. Khludnev, G. Leugering, On elastic bodies with thin rigid inclusions and cracks. Math. Methods Appl. Sci. 33 (2010),
1955–1967.
206. A. Khludnev, G. Leugering, M. Specovius-Neugebauer, Optimal control of inclusion and crack shapes in elastic bodies.
J. Optim. Theory Appl. 155 (2012), 54–78.
207. A. M. Khludnev, J. Sokolowski, The Griffith’s formula and the Rice–Cherepanov integral for crack problems with uni-
lateral conditions in nonsmooth domains. European J. Appl. Math. 10 (1999), 379–394.
208. V.A. Kovtunenko, G. Leugering, A shape-topological control problem for nonlinear crack - defect interaction: the anti-
plane variational model. SIAM J.Control Optim. 54 (2016), 1329–1351.
209. J.E. Lagnese, G. Leugering, E.J.P.G. Schmidt, Modeling, analysis and control of dynamic elastic multi-link structures.
Systems and Control: Foundations and Applications. Birkhäuser Boston, Inc., Boston, MA, 1994
210. J.E. Lagnese, G. Leugering, Domain decomposition methods in optimal control of partial differential equations. Interna-
tional Series of Numerical Mathematics, 148. Birkhäuser Verlag, Basel, 2004.
211. G. Leugering, J. Sokolowski, A. Zochowski, Shape-topological differentiability of energy functionals for unilateral prob-
lems in domains with cracks and applications. In: Optimization with PDE constraints; ESF Networking Program ’OPTPDE’
(R. Hoppe, editor) (2014), 203–221.
212. G. Leugering, J. Sokolowski, A. Zochowski, Control of crack propagation by shape-topological optimization. Discrete
and Continuous Dynamical Systems - Series A (DCDS-A) 35 (2015), 2625–2657.
213. V.G. Maz’ya, S.A. Nazarov, B.A. Plamenevski, Asymptotic Theory of Elliptic Boundary Value Problems in Singularly
Perturbed Domains. Birkhäuser, Basel, 2000.
214. Zh.O. Oralbekova, K.T. Iskakov, A.L. Karchevsky, Existence of the residual functional derivative with respect to a
coordinate of gap point of medium. Appl. Comput. Math. (2013) 12, 222–233.
215. N. Ovcharova, J. Gwinner, From solvability and approximation of variational inequalities to solution of nondifferentiable
optimization problems in contact mechanics. Optimization 64 (2015), 1683–1702.
6.5 Definitions 91
6.5 Definitions
we add at the end of the book all definitions

In Chapter ?? we use the notation.
n, Ω ,Γ Sobolev spaces W s,p (Ω ) ⊂ L p (Ω )
For boundary value problems in Ω , the solution y(Ω ) is a function in the Sobolev spaces W s,p (Ω ) ⊂
L (Ω ). For the exponent p := 2 and for the integer s := k, the spaces are denoted by H k (Ω ), usually k = 1, 2,
p
1
however we need also the fractional Sobolev space H 2 (Γ ) for the traces on the boundary Γ of the functions
in H 1 (Ω ) .
We need also the family of characteristic functions of the domains in RN , the characteristic function of Ω
is denoted by χ Ω or χ ,

χ Ω (x) = 1, x ∈ Ω (6.5.1)
0, otherwise.
The complement of a domain Ω is denoted by Ω c , it is the complement in the hold all set = B \ Ω (or in
the whole space RN \ Ω , usually it is indicated what we mean by the symbol of complement).
The function χ Ω c or χ c stands for the characteristic function of the complement Ω c , thus χ c equals one
outside of Ω ,
χ c (x) = 0, x ∈ Ω (6.5.2)
1, otherwise.
|Ω | is N-dimensional Lebesgue measure of Ω ,

Z
|Ω | = χΩ .
The set of characteristic functions defined in the hold all domain B is denoted by Char(B)= {χ ∈ L2 (B)
: (1 − χ )χ = 0 a.e. in B}. If the hold all domain is bounded then the set of characteristic functions is
sequentially compact, e.g. in L2 (B), with respect to the strong convergence. However, the set Char(B) is
neither nonconvex, nor closed with respect to the weak-star convergence of characteristic functions. Thus,
a limit of a weakly convergent sequence of characteristic functions belongs to the convex hull of the set
Char(B), the convex hull contains measurable functions such that 0 ≤ ψ ≤ 1.
For the differential expressions on the boundary, the following notation is used,
1. ∇Γ is tangential gradient on Γ ,
2. ∂∂n is normal derivative on Γ ,
3. ∂∂n is conormal derivative on Γ associated to the matrix A,
A
4. divΓ stands for the tangential divergence on Γ ,
5. ∆Γ is the Laplace-Beltrami operator on Γ .
Finally, the symbol y ∗ Tt stands for the distribution y ∈ D 0 (Ωt ) transported to the reference domain Ω .
C0∞ (Ω ), C0∞ (Ω ; Rd ) are the spaces of smooth functions (or of smooth vectorial functions), with compact
supports in Ω .
Appendix A
Program Code
A.0.1 Shared functions for Optimization Examples Chapter (??)
A.0.1.1 Assembling routines
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 function A = H1_D0_S(mesh,DoF,Basis,areaMarker,B)
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 tic
5
6 assert(size(mesh.p,1) == 2)
7 if ˜iscell(B)
8 B = {B};
9 end
10 nDoF = max(DoF(:));
11 nRhs = length(B);
12
13 qweights = permute(Basis.Int.d0.weights,[3 1 2]);
14 qpoints = Basis.Int.d0.points;
15
16 ind = ismember(mesh.ID,areaMarker);
17 if all(˜ind); A = sparse(nDoF,nRhs); return; end
18
19 mesh.t = mesh.t(:,ind);
20 DoF = DoF(:,ind);
21
22 NoDoFpT = size(DoF,1);
23 NoT = size(DoF,2);
24
25 [i2, i1] = meshgrid(1,1:NoDoFpT);

26 it1 = DoF(i1(:),:);
27 it2 = repmat(i2,1,NoT);
28 IT1 = [];
29 IT2 = [];
30 ALOC = [];
31
32 meshTransform = getTransformation(mesh);
33 qpFT = transformQP(meshTransform,qpoints);
34 DT = transformH1d0(meshTransform,Basis);
35 A = sparse(0);
36 for nf = 1:length(B)
37 if ˜isempty(B{nf});
38 it2 = nf + 0*it1(:);
39 D = B{nf}(qpFT(1,:,:),qpFT(2,:,:));
40 Aloc = zeros(NoDoFpT,NoT);
41 for i = 1:NoDoFpT
42 f = DT(:,:,:,i);
43 Aloc(i,:) = ...
permute(sum(bsxfun(@times,sum(f.*D(1,:,:),1),qweights),3).*(meshTransform.absT),[1
3 2]);
44 end
93
94 A Program Code
45 IT1 = [IT1;it1(:)];
46 IT2 = [IT2;it2(:)];
47 ALOC = [ALOC;Aloc(:)];
48 end
49 end
50 A = sparse(IT1,IT2,ALOC,nDoF,nRhs);
51 end
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 function A = H1_D0_M_D0(mesh,DoF,Basis,areaMarker,B)
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 tic
6

9
11 if all(˜ind); A = sparse(nDoF,nDoF); return; end
12
15
18
19 [i2, i1] = meshgrid(1:NoDoFpT,1:NoDoFpT);

20 it1 = DoF(i1(:),:);
21 it2 = DoF(i2(:),:);
22
24 absT = (meshTransform.absT);
25
26 DIM = size(mesh.p,1);
27
28 if nargin == 4 || isempty(B);
29 type = 'unity';
30 else
31 if DIM == 3; b = B(0,0,0); end
32 if DIM == 2; b = B(0,0); end
33 [rows,cols] = size(b);
34 if cols ˜=1; error('check B'); end
35 if rows == 1
36 type = 'scalar';
37 elseif rows == DIM
38 type = 'diagonal';
39 elseif rows == DIM*(DIM+1)/2
40 type = 'symm';
41 elseif rows == DIMˆ2;
42 type = 'matrix';
43 end
45 Dxy = evalAtQuad(qpFT,B);
46
47 end
48
50 Aloc = zeros(NoDoFpT,NoDoFpT,NoT);
51 [i1,i2] = ndgrid(1:DIM,1:DIM);
52 switch type
53 case 'unity'
54 for i =1:NoDoFpT
55 f = DT(:,:,:,i);
56 Aloc(i,i,:) = permute((permute(sum(f.*f,1),[2 3 ...
1])*qweights(:)).'.*absT,[1 3 2]);
57 for j = i+1:NoDoFpT
58 g = DT(:,:,:,j);
A Program Code 95
59 Aloc(i,j,:) = permute((permute(sum(f.*g,1),[2 3 ...

1])*qweights(:)).'.*absT,[1 3 2]);
60 Aloc(j,i,:) = Aloc(i,j,:);
61 end
62 end
63 case 'scalar'
64 for i =1:NoDoFpT
65 f = DT(:,:,:,i);
66 Aloc(i,i,:) = ...
permute(sum(bsxfun(@times,sum(f.*f,1).*Dxy,qweights),3).*absT,[1 ...
3 2]);
68 g = DT(:,:,:,j);
69 Aloc(i,j,:) = ...
permute(sum(bsxfun(@times,sum(f.*g,1).*Dxy,qweights),3).*absT,[1 ...
3 2]);
71 end
72 end
73 case 'diagonal'
74 for i =1:NoDoFpT
75 f = DT(:,:,:,i);
76 Aloc(i,i,:) = ...
permute(sum(bsxfun(@times,sum(f.*f.*Dxy,1),qweights),3).*absT,[1 ...
3 2]);
78 g = DT(:,:,:,j);
79 Aloc(i,j,:) = ...
permute(sum(bsxfun(@times,sum(f.*g.*Dxy,1),qweights),3).*absT,[1 ...
3 2]);
81 end
82 end
83 case 'symm'
84 if DIM == 3; Dxy = Dxy([1 2 3 2 4 5 3 5 6],:,:);end
85 if DIM == 2; Dxy = Dxy([1 2 2 3],:,:);end
86 for i =1:NoDoFpT
87 f = DT(i1(:),:,:,i);
88 g = DT(i2(:),:,:,i);
89 Aloc(i,i,:) = ...
3 2]);
91 g = DT(i2(:),:,:,j);
92 Aloc(i,j,:) = ...
3 2]);
94 end
95 end
96
97 case 'matrix'
98 for i =1:NoDoFpT
99 f = DT(i1(:),:,:,i);
100 for j = 1:NoDoFpT
101 g = DT(i2(:),:,:,j);
102 Aloc(i,j,:) = ...
3 2]);
103 end
104 end
105 otherwise
106 error('No type found');
107 end
108 A = sparse(it2(:),it1(:),Aloc(:),nDoF,nDoF);
109 end
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
96 A Program Code
2 function [A,dA] = H1_D1_M_D1(mesh,DoF,Basis,areaMarker,B)

3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 tic
6

9
12
15
18

20 it1 = DoF(i1(:),:);
21 it2 = DoF(i2(:),:);
22
25
27
28 if nargin == 4 || isempty(B);
29 type = 'unity';
30 else
31 if isa(B,'function_handle')
32 if DIM == 3; b = B(0,0,0); end
33 if DIM == 2; b = B(0,0); end
34 [rows,cols] = size(b);
35 if cols ˜=1; error('check B'); end
36 if rows == 1
37 type = 'scalar';
38 elseif rows == DIM
39 type = 'diagonal';
40 elseif rows == DIM*(DIM+1)/2
41 type = 'symm';
42 elseif rows == DIMˆ2;
43 type = 'matrix';
44 end
46 Dxy = evalAtQuad(qpFT,B);
47
48 else
49 type ='element';
50 Dxy = B;
51 end
52
53 end
54
57 [i1,i2] = ndgrid(1:DIM,1:DIM);
58 switch type
59 case 'unity'
60 for i =1:NoDoFpT
61 f = DT(:,:,:,i);
62 Aloc(i,i,:) = ...
permute(sum(bsxfun(@times,sum(f.*f,1),qweights),3).*absT,[1 ...
3 2]);
64 g = DT(:,:,:,j);
65 Aloc(i,j,:) = ...
permute(sum(bsxfun(@times,sum(f.*g,1),qweights),3).*absT,[1 ...
3 2]);
67 end
A Program Code 97
68 end
69 case 'element'
70 for i =1:NoDoFpT
71 f = DT(:,:,:,i);
72 Aloc(i,i,:) = ...
permute(sum(bsxfun(@times,sum(f.*f,1),qweights),3).*absT,[1 ...
3 2]);
74 g = DT(:,:,:,j);
75 Aloc(i,j,:) = ...
permute(sum(bsxfun(@times,sum(f.*g,1),qweights),3).*absT,[1 ...
3 2]);
77 end
78 end
79 dA = Aloc;
80 Aloc = bsxfun(@times,Aloc,reshape(Dxy,1,1,[]));
81
82 case 'scalar'
83 for i =1:NoDoFpT
84 f = DT(:,:,:,i);
85 Aloc(i,i,:) = ...
permute(sum(bsxfun(@times,sum(f.*f,1).*Dxy,qweights),3).*absT,[1 ...
3 2]);
87 g = DT(:,:,:,j);
88 Aloc(i,j,:) = ...
permute(sum(bsxfun(@times,sum(f.*g,1).*Dxy,qweights),3).*absT,[1 ...
3 2]);
90 end
91 end
92 case 'diagonal'
93 for i =1:NoDoFpT
94 f = DT(:,:,:,i);
95 Aloc(i,i,:) = ...
permute(sum(bsxfun(@times,sum(f.*f.*Dxy,1),qweights),3).*absT,[1 ...
3 2]);
97 g = DT(:,:,:,j);
98 Aloc(i,j,:) = ...
3 2]);
100 end
101 end
102
103 case 'symm'

104 if DIM == 3; Dxy = Dxy([1 2 3 2 4 5 3 5 6],:,:);end
105 if DIM == 2; Dxy = Dxy([1 2 2 3],:,:);end
106 for i =1:NoDoFpT
107 f = DT(i1(:),:,:,i);
108 g = DT(i2(:),:,:,i);
109 Aloc(i,i,:) = ...
3 2]);
111 g = DT(i2(:),:,:,j);
112 Aloc(i,j,:) = ...
3 2]);
114 end
115 progress(i,NoDoFpT)
116 end
117
118 case 'matrix'

119 for i =1:NoDoFpT
120 f = DT(i1(:),:,:,i);
98 A Program Code
122 g = DT(i2(:),:,:,j);
123 Aloc(i,j,:) = ...
3 2]);
124 end
125 end
126 otherwise
127 error('No type found');
128 end
130 end
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 function [A,dA] = ...
H1_D1_M_D0_mixed(mesh,DoF,DoF_test,Basis,Basis_test,areaMarker,B)
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 tic
6 nDoF_test = max(DoF_test(:));
7
8
11
14
17 DoF_test = DoF_test(:,ind);
18
20 NoDoFpT_test = size(DoF_test,1);
22
23 [i2, i1] = ndgrid(1:NoDoFpT_test,1:NoDoFpT);
24 it1 = DoF(i1(:),:);
25 it2 = DoF_test(i2(:),:);
26
29
31 DT_test = transformH1d1(meshTransform,Basis_test);
32 AlocX = zeros(NoDoFpT_test,NoDoFpT,NoT);
33 AlocY = zeros(NoDoFpT_test,NoDoFpT,NoT);
34
35 for i =1:NoDoFpT_test
36 f = DT_test(:,:,:,i);
38 g = DT(:,:,:,j);
39 AlocX(i,j,:) = ...
permute(sum(bsxfun(@times,sum(f(1,:,:).*g,1),qweights),3).*absT,[1 ...
3 2]);
40 IT1(i,j,:) = DoF_test(i,:);
41 IT2(i,j,:) = DoF(j,:);
42 AlocY(i,j,:) = ...
permute(sum(bsxfun(@times,sum(f(2,:,:).*g,1),qweights),3).*absT,[1 ...
3 2]);
43 end
44 end
45 AX = sparse(it2(:),it1(:),AlocX(:),nDoF_test,nDoF);
46 AY = sparse(it2(:),it1(:),AlocY(:),nDoF_test,nDoF);
47 A = cat(1,AX,AY);
48 end
A Program Code 99
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 function A = H1_G0_V(mesh,DoF,Basis,areaMarker,edgeMarker,B)
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 assert(DIM == 2,'Function is only for 2D')
7
8 % restricting to areaMarker
12
13 % adjacency
14 locEdge = [2 3; 1 3 ; 1 2]';
15 idx = ismember(mesh.be,edgeMarker);
16 e = sort(mesh.e(:,idx));
17 Edg = ...
reshape(ismember(sort(reshape(mesh.t(locEdge(:),:),2,[]))',e','rows'),3,[]);
18
19 % restricting to adjacent triangles

20 ind = any(Edg);
21 Edg = Edg(:,ind);
22
23 % check if anything has to be done

25
28
31
32 % modify weights
33 qweights = Basis.Int.g0.weights(:);
34 nqp = length(qweights);
35
37
38 qp = Basis.Int.g0.points;
39 qpFT = transformQP(meshTransform,qp);
40
41 D = B(qpFT(1,:,:),qpFT(2,:,:));
42
43 locEdge = [2 3; 1 3 ; 1 2]';
44 meshTransform.edgeLength = ...
reshape(sqrt(sum((mesh.p(:,mesh.t(locEdge(1,:),:)) - ...
mesh.p(:,mesh.t(locEdge(2,:),:))).ˆ2)),3,[]);
45 edgelen = meshTransform.edgeLength.*Edg;
46
47 el = reshape(repmat(edgelen(:)',nqp/3,1),nqp,[]);
48 qweights = bsxfun(@times,qweights(:),el);
49
50 % transform basis functions

51 DT = transformH1g0(meshTransform,Basis);
52
53 % assembling
54 Aloc = zeros(NoDoFpT,1,NoT);
55 for i =1:NoDoFpT
56 f = DT(:,:,:,i);
57 Aloc(i,1,:) = permute(sum(permute(sum(f.*D,1),[2 3 1]).*qweights',2),[1 ...
3 2]);
58 end
59
60 [i2, i1] = meshgrid(1:NoDoFpT,1);

61 it1 = 1+0*DoF(i1(:),:);
62 it2 = DoF(i2(:),:);
63
64 A = sparse(it2(:),it1(:),Aloc(:),nDoF,1);
65 end
100 A Program Code
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 function A = H1_G0_M_G0(mesh,DoF,Basis,areaMarker,edgeMarker)
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 assert(DIM == 2,'Function is only for 2D')
7
8 % restricting to areaMarker
12
13 % adjacency
14 locEdge = [2 3; 1 3 ; 1 2]';
15 idx = ismember(mesh.be,edgeMarker);
17 Edg = ...
reshape(ismember(sort(reshape(mesh.t(locEdge(:),:),2,[]))',e','rows'),3,[]);
18
19 % restricting to adjacent triangles

20 ind = any(Edg);
21 Edg = Edg(:,ind);
22
23 % check if anything has to be done
25
28
31
32 % modify weights
33 qweights = Basis.Int.g0.weights(:);
34 nqp = length(qweights);
35
37 locEdge = [2 3; 1 3 ; 1 2]';
38 meshTransform.edgeLength = ...
reshape(sqrt(sum((mesh.p(:,mesh.t(locEdge(1,:),:)) - ...
mesh.p(:,mesh.t(locEdge(2,:),:))).ˆ2)),3,[]);
39 edgelen = meshTransform.edgeLength.*Edg;
40
41 el = reshape(repmat(edgelen(:)',nqp/3,1),nqp,[]);
42 qweights = bsxfun(@times,qweights(:),el);
43
44 % transform basis functions

45 DT = transformH1g0(meshTransform,Basis);
46
47 % assembling
49 for i =1:NoDoFpT
50 f = DT(:,:,:,i);
51 Aloc(i,i,:) = permute(sum(permute(sum(f.*f,1),[2 3 1]).*qweights',2),[1 ...
3 2]);
53 g = DT(:,:,:,j);
54 Aloc(i,j,:) = permute(sum(permute(sum(f.*g,1),[2 3 ...
1]).*qweights',2),[1 3 2]);
56 end
57 end
58

60 it1 = DoF(i1(:),:);
61 it2 = DoF(i2(:),:);
62
64 end
A Program Code 101
A.0.1.2 basis transformation routines
1 function [varargout] = getTransformation(mesh)

2 nVpT = size(mesh.t,1);
3 p = mesh.p;
4 t = mesh.t;
5
6 NoT = size(t,2);
7 x = p(1,:);
8 y = p(2,:);
9 T = zeros(2,2,NoT);
10 T(1,1,:) = x(t(2,:))-x(t(1,:));
11 T(1,2,:) = x(t(3,:))-x(t(1,:));
12 T(2,1,:) = y(t(2,:))-y(t(1,:));
13 T(2,2,:) = y(t(3,:))-y(t(1,:));
14 detT = T(1,1,:).*T(2,2,:) - T(2,1,:).*T(1,2,:);
15 absT = abs(detT)/2;
16 absT = (absT(:))';
17
18 adjT(1,1,:) = T(2,2,:);
19 adjT(1,2,:) = -T(1,2,:);
20 adjT(2,1,:) = -T(2,1,:);
21 adjT(2,2,:) = T(1,1,:);
22 invT = bsxfun(@times,adjT,1./detT);
23 detT = (detT(:))';
24 d = p(:,t(1,:));
25 fceArea = 0;
26 meshTransform.absT = absT;
27 meshTransform.invT = invT;
28 meshTransform.T = T;
29 meshTransform.d = d;
30 meshTransform.detT = detT;
31 meshTransform.fceArea = fceArea;
32 if nargout == 1
33 varargout{1} = meshTransform;
34 else
35 varargout = struct2cell(meshTransform);
36 end
37 end
1 function qpFT = transformQP(meshTransform,qpoints)

2 T = meshTransform.T;
3 d = meshTransform.d;
4 DIM = size(T,1);
5 NoT = size(T,3);
6 nqp = size(qpoints,1);
7 FT = reshape(permute(T,[1 3 2]),DIM*NoT,DIM);
8 qpFT = FT * qpoints(:,2:DIM+1)';
9 qpFT = reshape(qpFT,DIM,NoT,nqp);
10 qpFT = bsxfun(@plus,qpFT,d);
11 end
1 function phiT = transformH1d0(meshTransform,Basis)

2 phi = [
3 Basis.Vtx.d0;
4 Basis.Edg.d0;
5 Basis.Fce.d0;
6 Basis.Cll.d0
7 ];
9 NoT = size(T,3);
10 phiT = repmat(permute(phi,[2 4 3 1]),[1,NoT,1,1]);
102 A Program Code
1 function nablaT = transformH1d1(meshTransform,Basis)

2 dphi = [
3 Basis.Vtx.d1;
4 Basis.Edg.d1;
5 Basis.Fce.d1;
6 Basis.Cll.d1
7 ];
8 invT = meshTransform.invT;
9 NoT = size(invT,3);
10 NoDoFpT = size(dphi,1);
11 DIM = size(dphi,2);
12 nqp = size(dphi,3);
13 TT = reshape(permute(invT,[2 3 1]),DIM*NoT,DIM);
14 DD = reshape(permute(dphi,[2 3 1]),DIM,NoDoFpT*nqp);
15 DT = TT*DD;
16 nablaT = reshape(DT,DIM,NoT,nqp,NoDoFpT);
1 function phiT = transformH1g1(meshTransform,Basis)

2 phi = [
3 Basis.Vtx.g0;
4 Basis.Edg.g0;
5 Basis.Fce.g0;
6 Basis.Cll.g0
7 ];
9 NoT = size(T,3);
10 phiT = repmat(permute(phi,[2 4 3 1]),[1,NoT,1,1]);
1 function nablaT = transformH1g1(meshTransform,Basis)

2 dphi = [
3 Basis.Vtx.g1;
4 Basis.Edg.g1;
5 Basis.Fce.g1;
6 Basis.Cll.g1
7 ];
8 invT = meshTransform.invT;
9 NoT = size(invT,3);
10 NoDoFpT = size(dphi,1);
11 DIM = size(dphi,2);
12 nqp = size(dphi,3);
13 TT = reshape(permute(invT,[2 3 1]),DIM*NoT,DIM);
14 DD = reshape(permute(dphi,[2 3 1]),DIM,NoDoFpT*nqp);
15 DT = TT*DD;
16 nablaT = reshape(DT,DIM,NoT,nqp,NoDoFpT);
A.0.1.3 basis generation and degrees of freedom
1 % getLagrangeBasis Lagrange basis function

2 % basis = getLagrangeBasis(order)
3 % returns Lagrange basis functions of polynomial order <order>
4 function basis = getLagrangeBasis(order)
5 [˜,˜,˜] = mkdir(mfilename('fullpath'));
6 switch order
7 case 1
8 load(fullfile(mfilename('fullpath'),'P1Tri'));
9 case 2
10 load(fullfile(mfilename('fullpath'),'P2Tri'));
11 otherwise
12 error('Choose basis order 1 or 2');
13 end
14 end
A Program Code 103
1 function DoF = getDoFs(mesh,basis)

2 t = mesh.t;
3
4 nV = basis.Vtx.nDoF;
5 nE = basis.Edg.nDoF;
6 nF = basis.Fce.nDoF;
7 nC = basis.Cll.nDoF;
8
9 locEdge = [2 3; 1 3; 1 2]';
10 nEpT = 3;
11 nFpT = 0;
12
13
14 maxDoF = 3e6;
15
16 DoF = [];
17 nDoF = 0;
18
19 % vertex DoF
20 if nV > 0
21 Vtx0 = t;
22 Vtx = [];
23 nVtx = 0;
24 for i = 1:nV
25 Vtx = [Vtx; Vtx0 + nVtx];
26 nVtx = max(Vtx(:));
27 end
28 DoF = [DoF;Vtx + nDoF];
30 if (nDoF > maxDoF) DoF = -1, return; end;
31 end
32
33 % edge DoF
34 if nE > 0
35 a1 = locEdge(:);
36 a2 = t(a1,:);
37 a3 = reshape(a2,2,[]);
38 a4 = sort(a3);
39 [˜,˜,a5] = unique(a4','rows');
40 Edg0 = reshape(a5,nEpT,[]);
41
42 Edg = [];
43 nEdg = 0;
44 for i = 1:nE
45 Edg = [Edg; Edg0 + nEdg];
46 nEdg = max(Edg(:));
47 end
48 DoF = [DoF;Edg + nDoF];
50 if (nDoF > maxDoF) DoF = -1, return; end;
51
52 end
53 end
1 function onBoundary = getBoundaryDoFs(mesh,DoF,Basis,faceMarker)

2 nVpT = size(mesh.t,1);
3 nVpE = 2;
4 nVpF = 3;
5
6 DIM = 2;
7 locEdge = [2 3; 1 3; 1 2]';
8 nEpT = 3;
9 nFpT = 0;
10
11
12 onBoundary = struct('Vtx',[],'VtxM',[],'Edg',[],'EdgM',[],'Fce',[],'FceM',[]);
13
104 A Program Code
14 if DIM >=1 && Basis.Vtx.nDoF > 0

15 t = mesh.t;
16 switch DIM
17 case 3
18 idx = ismember(mesh.bf,faceMarker);
19 f = mesh.f(:,idx);
20 Vtx = ismember(t,f);
21 case 2
22 idx = ismember(mesh.be,faceMarker);
23 e = mesh.e(:,idx);
24 Vtx = ismember(t,e);
25 end
26 idx1 = false(size(DoF));
27 for i = 1:Basis.Vtx.nDoF
28 idx1((1:nVpT) + (i-1)*nVpT,:) = Vtx;
29 end
30 onBoundary.Vtx = unique(DoF(idx1));
31 onBoundary.VtxM = Vtx;
32 end
33
34
35 if DIM >= 2 && Basis.Edg.nDoF > 0
36 t = mesh.t;
37 % e = mesh.e;
38 switch DIM
39 case 2
40
41 idx = ismember(mesh.be,faceMarker);
43 case 3
44 idx = ismember(mesh.bf,faceMarker);
45 f = mesh.f(:,idx);
46 locEdge2D = [2 3; 1 3; 1 2]';
47 e = f(locEdge2D(:),:);
48 e = unique(sort(reshape(e,2,[]))','rows')';
49 end
50 a1 = locEdge(:);
51 a2 = t(a1,:);
52 a3 = reshape(a2,nVpE,[]);
53 a4 = sort(a3);
54 e1 = sort(e);
55 a5 = ismember(a4',e1','rows');
56 Edg = reshape(a5,nEpT,[]);
57
58 idx2 = false(size(DoF));
59 for i = 1:Basis.Edg.nDoF
60 idx2((1:nEpT) + (i-1)*nEpT + nVpT*Basis.Vtx.nDoF,:) = Edg;
61 end
62 onBoundary.Edg = unique(DoF(idx2));
63 onBoundary.EdgM = Edg;
64
65 end
66
67 onBoundary.all = [onBoundary.Vtx',onBoundary.Edg',onBoundary.Fce'];
68 onBoundary.allM = ismember(1:max(DoF(:)),onBoundary.all);
69
70 end
1 function mesh = extendP2(mesh)

2 x = mesh.p(1,:);
3 y = mesh.p(2,:);
4
5 locEdge = [2 3; 1 3 ; 1 2]';
6 [edge] = unique(reshape(mesh.t(locEdge(:),:),2,[])','rows')';
7 pEdge = [mean(x(edge));mean(y(edge))];
8 pNode = mesh.p(:,unique(mesh.t(:)));
9 mesh.p = [pNode pEdge];
10 end
A Program Code 105
A.0.1.4 gradient projection
1 % /////////////////////////////////////////////////////////////////////////
2 % |[sdJX,sdJY] = getDescentDirection(mesh,SEL,dJX,dJY,params)
3 % | calculates smoothed descent direction
4 % |
5 % | Input:
6 % | mesh: geometry description
7 % | SEL: boundary selections
8 % | dJX / dJY: shape gradient for x and y deformation
9 % | params: options for smoothing
10 % |
11 % | Output:
12 % | dJX / dJY: smooth descent direction for x and y deformation
13 % \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
14
15 function [sdJX,sdJY] = getDescentDirection(mesh,SEL,dJX,dJY,params)

16 % get DoFs
17 basis = getLagrangeBasis(1);
18 DoF = getDoFs(mesh,basis);
19
20 % smoothing matrix
21 K = H1_D1_M_D1(mesh,DoF,basis,0);
22 M = H1_D0_M_D0(mesh,DoF,basis,0);
23 MOL = K + params.M_fac * M;
24
25 % fixed boundary
26 Vx0 = getBoundaryDoFs(mesh,DoF,basis,SEL.e_Vx_0);
27 Vx0 = Vx0.allM;
28 Vy0 = getBoundaryDoFs(mesh,DoF,basis,SEL.e_Vy_0);
29 Vy0 = Vy0.allM;
30
31 VVx = MOL;
32 VVx(Vx0,:) = 0;
33 VVx(Vx0,Vx0) = speye(sum(Vx0));
34 dJX(Vx0) = 0;
35 M2x = M;
36 M2x(Vx0,:)=0;
37 M2x(Vx0,Vx0)=speye(sum(Vx0));
38
39 VVy = MOL;
40 VVy(Vy0,:) = 0;
41 VVy(Vy0,Vy0) = speye(sum(Vy0));
42 dJY(Vy0) = 0;
43 M2y = M;
44 M2y(Vy0,:)=0;
45 M2y(Vy0,Vy0)=speye(sum(Vy0));
46
47 % smoothing
48 sdJX = VVx\dJX;
49 sdJY = VVy\dJY;
50 for i = 2:params.smooth_steps
51 sdJX = VVx\(M2x*sdJX);
52 sdJY = VVy\(M2y*sdJY);
53 end
54 end

BOOK2019 Wyklad

Uploaded by

Copyright:

Available Formats

BOOK2019 Wyklad

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

BOOK2019 Wyklad

Uploaded by

Copyright:

Available Formats

Guenter Leugering, Jan Sokolowski and Antoni Zochowski

The Shape Derivative Method. Theory and

1 Preliminaries and notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

2 Shape optimization for problems on metric graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

3 Mixed formulation of elliptic boundary value problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

4 Lagrangian formalism for weak solutions and shape optimization . . . . . . . . . . . . . . . . . . . . . . . 47

5 The second-order shape differentiability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

6 Shape optimization with MATLAB in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

2.1 Problems on graphs

2.2 Planar graphs

JDt := {J ∈ JS |yi (vJ ) · ei = 0}

JDn := {J ∈ JS |yi (vJ ) · e⊥

J̃Dc := {J ∈ JS |yi (vJ ) · e⊥

where ai ≤ bi for all i ∈ ID , D ∈ JDc .

Fig. 2.2 A general graph

We consider the energy of the system

V := {y : G → Rnp |yi ∈ H 1 (Ii ) (2.2.8)

V is clearly a Hilbert space in

K := V ∩ {(yi )ni=1 |yi (vD ) · e⊥

The Ritz-approach to deriving the problem now can be stated as follows

In order to define proper variations in (2.2.16), we introduce the Hilbert space

Obviously, if r ∈ K then ŷ = y + ϕ ∈ K , ∀ϕ ∈ V 0 . Taking these variations we obtain from (2.2.16) the

Now (2.2.20) clearly implies the strong formulation of the problem:

ŷi = yi + ψi , ψ ∈ V , ψi (vD ) = 0, i ∈ ID , D ∈ JDc

Fig. 2.3 The three-star

Example 2.2.1 Let f i = 0, ci = 0, i = 1, . . . , m then

Hence, the explicit solution is given by

Example 2.2.2 Let f i = 0, ci > 0, i = 1, . . . , m. Then the solutions are given by

The boundary and node conditions give

Example 2.2.3 A Helmholtz-type example is given by f i = 0, ci = −ω 2 ki < 0, where ω denotes a frequency.

The same procedure asin the previous examples provides

Exercise 2.2. Write a Matlab routine to solve the problem (2.2.27).

Exercise 2.3. Repeat Excercise (2.1), (2.2) for (2.2.36).

Remark 2.2.4 If we introduce the vectors

and the matrices

2.2.1 Steklov-Poincaré map for metric graphs

The corresponding Steklov-Poincaré map is then given by

For the second transmission condition according to (2.2.36), we obtain:

The corresponding Steklov-Poincaré map is then given by

2.2.2 Problems on general metric graphs

The Steklov-Poincaré or Dirichlet-Neumann maps, together with the corresponding Neumann-Dirichlet

= hq, Aϕi, ∀q, ϕ ∈ D(A),

where the domain D(A) is given by

D(A) := y = (yi )|i∈I ∈ Πi∈I H 2 (0, 1)|yi (n j ) = 0, i ∈ I j , j ∈ JDS

Thus, we can define the unbounded, self-adjoint operator

V := y = (yi )|i∈I ∈ Πi∈I H 1 (0, 1)|yi (n j ) = 0, i ∈ I j , j ∈ JDS

2.2.3 Shape variations for graph problems

See figure 2.2.3. The edges have lengths `i := ka − ai k, i = 1, . . . , m

Fig. 2.4 The deformed three-star

Example 2.2.7 We let ci = 0, i = 1, . . . , m and yd (x) = 0, x ∈ [0, `i ], i = 1, . . . , m. According to Fig. 2.2.3,

Let us observe that the unique solution is given by

On the fixed graph we have

(yi )t (x) = yti (Tt (v)(x) (2.2.65)

ẏi (0) = 0, i = 1, . . . , m, (2.2.66a)

which is independent of i and

For the shape derivatives we get

Fig. 2.5 The three-star with transformed edges

We differentiate with respect to t and obtain

We evaluate (2.2.71) at t = 0 and obtain

= cosh(Ai `i )−1 ui −cosh(Ai−1