On second order shape optimization methods for electrical impedance
  tomography
On second order shape optimization methods for electrical impedance
tomography.
L. Afraites, M. Dambrine and D. Kateb∗
October 29, 2018
Abstract
This paper is devoted to the analysis of a second order method for recovering the a priori un-
known shape of an inclusion ω inside a body Ω from boundary measurement. This inverse problem
- known as electrical impedance tomography - has many important practical applications and hence
has focussed much attention during the last years. However, to our best knowledge, no work has yet
considered a second order approach for this problem. This paper aims to fill that void: we investigate
the existence of second order derivative of the state u with respect to perturbations of the shape of
the interface ∂ω, then we choose a cost function in order to recover the geometry of ∂ω and derive
the expression of the derivatives needed to implement the corresponding Newton method. We then
investigate the stability of the process and explain why this inverse problem is severely ill-posed by
proving the compactness of the Hessian at the global minimizer.
Keywords: inverse problems, identification of inhomogenities, shape calculus, order two methods.
1 Introduction and statement of the results.
Let Ω be a bounded open set with smooth boundary in R2 or R3. Consider a L∞ function σ such
that there exists a real c with σ(x) ≥ c > 0. Consider the elliptic equation
− div
σ(x)∇u
= 0 in Ω, (1)
with the Dirichlet boundary condition
u = f on ∂Ω. (2)
Define the Dirichlet-to-Neumann map as
Λσ : f 7→ σ (∂nu)|∂Ω ,
where u solves (1),(2) and n is the outer unit normal vector to ∂Ω. The inverse conductivity problem
of Calderón is to determine σ from Λσ. Electrical impedance tomography aims to form an image
of the conductivity distribution σ from the knowledge of Λσ. When σ is smooth enough, one can
reconstruct σ from Λσ (see the works of Sylvester and Uhlmann [21], Nachmann [15, 16] and Novikov
∗Laboratoire de Mathématiques Appliquées de Compiègne. Université de Technologie de Compiègne.
http://arxiv.org/abs/0704.0708v1
[17]). When the conductivity distribution is only L∞, Astala and Päivärinta have recently shown in
[3] that, in dimension two, the map Λσ determines σ ∈ L
∞(Ω).
We are interested in a particular case of that problem: when a body is inserted inside a given
object with a distinct conductivity, the question of determining its shape from boundary measurement
arises in many fields of modern technology. In the context of the inverse problem of conductivity
of Calderón, we restrict the range of admissible conductivity distributions to the family of piecewise
constant functions which take only two distinct values σ1, σ2 > 0 which are assumed to be known.
The conductivity distribution is then defined by an open subset ω as
σ = σ1χΩ\ω + σ2χω. (3)
Here, the only unknown of the problem is ω a subdomain of Ω with a smooth boundary ∂ω; its outer
unit normal vector is denoted by n. The notation χω (respectively. χΩ\ω) denotes the characteristic
function of ω (respectively. Ω \ ω). The second main difference arises from practical considerations:
it is unrealistic from the point of view of applications to know the full graph of Dirichlet-to-Neumann.
Therefore, we will assume that one has access to a single point in that graph. This non destructive
testing problem is usually written from a numerical point of view as the minimization of a cost
function: typically a least-square matching criterion. Many authors have investigated the steepest
descent method for this problem [13, 7, 10, 18, 1] with the methods of shape optimization since the
unknown parameter is a geometrical domain.
This work is devoted to the study of second order methods for this problem that has only be
considered before for simplified models in [5, 2]. By introducing second order methods, one aims to
reach two distinct objectives.
• On one hand, we provide all the needed material to design a Newton algorithm. We will give
differentiability results for the state function and for the objective that we have chosen to study
in this work. Nevertheless, we point out that the discretization of a Newton method for this
problem turns out to be very delicate; this is why, in the present paper, we will neither discuss
about this problem nor present numerical examples. This topic is actually the main objective
of a work in progress.
• On the other hand, we analyze rigorously the well-posedness of the optimization method. This is
justified by the huge numerical literature devoted to the numerical study of this question in the
field of inverse problems; the numerical experiments insist on the ill-posedness of this problem.
We will explain the instability in the continuous settings in terms of shape optimization. We
show that the shape Hessian is not coercive -in fact its Riesz operator is compact – and this
explains the unstability of the minimization process.
Let us describe the precise problem under consideration and the notations. We consider a bounded
domain Ω ⊂ Rd (d = 2 or 3) with a C2 boundary. It is filled with a material whose conductivity is σ1
and with an unknown inclusion ω in Ω of conductivity σ2 6= σ1. We search to reconstruct the shape
of ω by measuring on ∂Ω, the input voltage and the corresponding output current. In the sequel, we
fix d0 > 0 and consider inclusions ω such that ω ⊂⊂ Ωd0 = {x ∈ ω, d(x, ∂Ω) > d0}. We also assume
that the boundary ∂ω is of class C4,α. The inverse problem arises when one has access to the normal
vector derivative of the potential u that solves (1)-(2) when the conductivity distribution is defined
by (3) . Assume that ones knows
σ1∂nu = g on ∂Ω, (4)
then the problem (1)-(2)-(4) is overdetermined. The electrical impedance tomography problem we
consider is to recover the shape of ω from the knowledge of the single Cauchy pair (f, g).
In order to recover the shape of the inclusion ω, an usual strategy is to minimize a cost function.
Many choices are possible; however it turns out that a Kohn and Vogelius type objective leads to a
minimization problem with nicer properties than the least squares fitting approaches (we refer to [1]
for a comparison of different objectives with order one methods and to [2] for the case of a perfectly
insulated inclusion). Therefore, we study such a cost function in this work.
Let us define this criterion. Its distinctive feature is to involve two state functions ud and un: the
state ud solves (1)-(2) while un solves (1)-(4). The Kohn -Vogelius objective JKV is then defined as:
JKV (ω) =
σ|∇(ud − un)|
2 (5)
Let us sum up the results of this paper concerning the minimization of this objective. We first prove
differentiability results for the state ud. In the sequel, we use the convention that a bold character
denotes a vector. If h denotes a deformation field, it can be written as h = hτ + hnn on ∂ω. Note
also that in the following lines, n denotes the outer normal field to ∂ω pointing into Ω\ω. Hence, for
x ∈ ∂ω, we define, when the limit exists, u±(x) (resp. (∂nu)
±(x)) as the limit of u(x± tn(x)) (resp.
〈∇u(x± tn(x),n(x))) when t > 0 tends to 0. Note that hτ is a vector while hn is a scalar quantity.
The admissible deformation fields have to preserve ∂Ω and the regularity of the boundaries:
therefore the space of admissible fields is
H = {h ∈ C4,α(Rd,Rd), Supp(h) ⊂ Ωd0}.
The following result concerns the first order derivative of the state functions ud and un. It was
derived in [7, 18, 1].
Theorem 1 Let Ω be an open smooth subset of Rd (d = 2 or 3) and let ω be an element of Ωd0 with
a boundary of class C4,α. Then the state functions ud and un are shape differentiable; furthermore
their shape derivative u′d and u
n belongs to H
1(Ω \ ω) ∪H1(ω) and satisfy


∆u′d = 0 in Ω \ ω and in ω,
d on ∂ω,[
= [σ]divτ (hn∇τud) on ∂ω,
u′d = 0 on ∂Ω.


∆u′n = 0 in Ω \ ω and in ω,
n on ∂ω,[
= [σ]divτ (hn∇τun) on ∂ω,
∂u′n = 0 on ∂Ω.
The main result of this work concerns the second order derivative. It is given is the following theorem.
Theorem 2 Let Ω be an open smooth subset of Rd (d = 2 or 3) and let ω be an element of Ωd0
with a C4,α boundary. Let h1 and h2 be two deformation fields in H. Then the state ud has a second
order shape derivative u′′d ∈ H
1(Ω \ ω) ∪H1(ω) that solves


∆u′′d = 0 in Ω \ ω and in ω,[
h1,nh2,nH − h1τ .(Dnh2τ )
[∂nud]−
h1,n[∂n(ud)
2] + h2,n[∂n(ud)
h1τ .∇h2,n + h2τ .∇h1,n
[∂nud] on ∂ω,[
= divτ
σ∇τ (ud)
+ h1,n
σ∇τ (ud)
+ h1τ .(Dnh2τ )[σ∇τud]
− divτ
(h1τ .∇τh2,n +∇τh1,n.h2τ ) [σ∇τud]
+ divτ
h2,nh1,n(2Dn−HI) [σ∇τud]
on ∂ω,
u′′d = 0 on ∂Ω.
Here, (ud)
i denotes the first order derivative of u in the direction of hi as given in (6), Dn stands for
the second fundamental form of the manifold ∂ω and H stands for the mean curvature of ∂ω. The
twin result concerning un is an easy adaption of Theorem 2. Once the differentiability of the state
function has been established, one can consider the objectives. In [1], we have shown the first order
result.
Theorem 3 Let Ω be an open smooth subset of Rd (d = 2 or 3) and let ω be an element of Ωd0
with a C4,α boundary. Let h1 and h2 be two deformation fields in H. The Kohn-Vogelius objective is
differentiable with respect to the shape and its derivative in the direction of a deformation field h is
given by:
DJKV (ω)h = [σ]
2 − |∂nu
+ |∇τud|
2 − |∇τun|
hn. (9)
We now give the second-order derivative of the Kohn and Vogelius criterion.
Theorem 4 Let Ω be an open smooth subset of Rd (d = 2 or 3) and ω be an element of Ωd0 with a
C4,α boundary. Let h1 and h2 be two deformation fields in H. The Kohn-Vogelius objective is twice
differentiable with respect to the shape and its second derivative in the directions h1 and h2 is given
D2JKV (ω)(h1,h2) =
σ|∇v|2
h1τ .∇τ (h2,n) + h2τ .∇τ (h1,n)− h2τ .(Dnh1τ )
σ|∇v|2
h1,nh2,n + 2
σ∇v.(h1,n∇v
2 + h2,n∇v
∂n(un)
2 + ∂n(un)
1 − ∂nv
1(ud)
2 − ∂nv
2(ud)
σ∂n(un)
− σ1∂nv
where we have set v = ud − un.
To investigate the properties of stability of this cost function, we are led to consider an admissible
inclusion ω∗ to solve both (1)-(2) and (1)-(4) in order to obtain the corresponding measurements f∗
and g∗. It is obvious that the domain ω∗ realizes the absolute minimum of the criterion JKV since,
by construction, we can write ud = un in Ω and hence JKV (ω
∗) = 0. We will check that the Euler
equation
DJKV (ω
∗)(h) = 0,
holds. We will also prove that
D2JKV (ω
∗)(h,h) =
σ|∇v′|2. (11)
Moreover, if hn 6= 0, then D
2JKV (ω
∗)(h,h) > 0 holds. Nevertheless, (11) does not mean that the
minimization problem is well-posed. In fact, it is the following theorem that explains the instability
of standard minimization algorithms.
Theorem 5 Assume that ω∗ is a critical shape of JKV for which the additional condition un = ud
holds. Then the Riesz operator corresponding to D2JKV (ω
∗) defined from H1/2(∂ω∗) with values in
H−1/2(∂ω∗) is compact. Moreover, the minimization problem is severely ill-posed in the following
sense: if the target domain is C∞ and if λn denotes the n
th eigenvalue of D2JKV (ω
∗), then λn =
o(n−s) for all s > 0.
Theorem 5 has two main consequences. First, the shape Hessian at the global minimizer is
not coercive. This means that this minimizer may not be a local strict minimum of the criterion.
Moreover, the criterion provides no control of the distance between the parameter ω and the target
ω∗. The second consequence concerns any numerical scheme used to obtain this optimal domain ω∗.
One has to face this difficulty and this explains why frozen Newton or Levenberg-Marquard schemes
have been used to solve numerically this problem [7, 1].
The paper is organized as follows. In a first section, we state some preliminary results. Some are
well known facts in shape optimization and will be recalled without proof for the sake of readability.
Some of them (e.g the derivatives of a Laplace-Beltrami operator and the tangential regularity of the
solution to (1)-(2) along the discontinuity of the conductivity distribution) are less known and will
be proved thanks to potential layer methods. Hence we will tackle the computations in Section 3
that we consider as the core of this work : it is essentially devoted to prove Theorem 2. After a first
part where we prove the existence of a second order derivative for the state, we propose two distinct
methods to find the boundary value problem solved by this second order derivative. The first method
(subsection 3.3) follows the lines of classical proofs of shape differentiability by differentiating the
weak formulation of problem (1)-(2) and interpreting the result in terms of differential operator and
boundary conditions. The alternative method (subsection 3.4) consists in a direct differentiation of
the boundary conditions. Finally, Section 4 is devoted to the analysis of the criterion, we establish
Theorem 4 and Theorem 5. We will present their consequences on the stability of critical shapes.
2 Preliminary results.
2.1 Elements of shape calculus
Before entering the proof of Theorem 2, we recall without proof some basic facts from shape opti-
mization (see [6] for references). Let h be a deformation field in C2(Ω,Rd) with ‖h‖C2 < 1. We set
Tt(h, .) = Id + th and denote by Ωt the transported domain Ωt = Tt(Ω). To avoid heavy notations,
we will misuse the notation Tt instead of Tt(h, .).
Material and shape derivatives. Classically, in mechanics of continuous media, the material
derivative is defined as being a positive limit. In our context, for any vector field h ∈ H, we define
the material derivative of the domain functional y = y(Ω) at Ω in an admissible direction h as the
limit
ẏ(Ω;h) = lim
y(Ωt) ◦ Tt − y(Ω)
, (12)
Similarly, one can define the material derivative ẏ(∂Ω,h) for any domain functional y = y(∂Ω) which
depends on ∂Ω. Another kind of derivative occurs : it is called the shape derivative of y(Ω,h). It is
viewed as a first local variation. Its definition is given by the following
Definition 1 The shape derivative y′ = y′(Ω;h) of a functional y(Ω) at Ω in the direction of a
vector field h is given by
y′ = ẏ − h.∇y. (13)
For more details on these derivations, the reader can consult [20, 6].
Elements of tangential derivatives. We will need in the sequel to manipulate the tangential
differential operators on a manifold. For the reader’s convenience, we recall from [4, 6] some definitions
and also some useful rules of calculus.
Definition 2 The tangential divergence of a vector field V ∈ C1(Rd,Rd) is given by
divτ (V) = div (V)−DV.n.n, (14)
where the notation DV denotes the Jacobian matrix of V. When the vector V ∈ C1(∂Ω,Rd) is
defined on ∂Ω, then the following notation is used to define the tangential divergence
divτ (V) = div
− (DṼ.n).n, (15)
where Ṽ stands for an arbitrary C1 extension of V on an open neighborhood of ∂Ω.
We introduce now, the notion of tangential gradient ∇τ of any smooth scalar function f in
C1(∂Ω,Rd).
Definition 3 Let an element f ∈ C1(∂Ω,Rd) be given and let f̃ be an extension of f in the sense
that f̃ ∈ C1(U) and f̃ |∂Ω = f and where U is an open neighborhood of ∂Ω. Then the following
notation is used to defined the tangential gradient
∇τf = ∇f̃ |∂Ω −∇f̃ .n n on ∂Ω. (16)
The details for the existence of such an extension can be found in [4]. Let us remark that these
definitions do not depend on the choice of the extension. Furthermore, one can show the important
relation ∫
∇τf.F = −
f divτ (F) , (17)
for all elements f ∈ C1(∂Ω) and all vector fields F ∈ C1(∂Ω,Rd) satisfying Fn = 〈F, n〉 = 0.
Integration by parts on ∂Ω. In general, the condition above Fn = 0 is not always satisfied.
We are then led to find another formula to extend the formula in the general case. The extension of
this integration by parts formula to fields with a normal vector component involves curvature.
First, we point out that the curvature is connected to the normal vector via the tangential
divergence operator. Recall that the mean curvature of ∂Ω is defined as H = divτ (n). Making use
of the form of divτ (n) on the boundary, one shows straightforwardly the following statement.
Proposition 1 Let Ω be an open subset of R3 with a C2 boundary. For any unitary extension N of
n on a neighborhood of ∂Ω, one has
div (N ) = H on ∂Ω.
Assume that the manifold ∂Ω has no borders. If F ∈ H2(∂Ω)3 and f ∈ H2(∂Ω), then we have
∇f.F+ fdivτ (F) =
(∇f.n+Hf)F.n. (18)
We assume now that the domain Ω has a C3 boundary. The simplest second-order derivative is the
Laplace Beltrami operator; it is defined as follows (see [20, 4, 6]) thanks to the following usual chain
rule.
Definition 4 Let f ∈ H2(∂Ω). The Laplace-Beltrami ∆τ of f is defined as follows
∆τf = divτ (∇τf) . (19)
There is a relation connecting the Laplace operator and the Laplace-Beltrami operator. Let us denote
by ∂2nnf = (D
2f.n).n where D2f stands for the Hessian of f .
Proposition 2 Let Ω be a domain with a boundary ∂Ω of class C3. For all functions f ∈ H3(Ω), it
holds
∆f = ∆τf +H∂nf + ∂
nnf, on ∂Ω. (20)
We need to compute shape and material derivative of special vector fields: the outer unit normal
vector n, the tangential gradient and the Laplace-Beltrami operator applied to a function. While the
derivative of the normal vector is obtained by a straightforward calculus, we have to transport from
∂Ωt to ∂Ω the Laplace-Beltrami operator and the tangential gradient in order to compute the other
derivatives.
Derivatives of the normal vector. We describe the material and shape derivatives of the normal
vector. We will denote by n the gradient of the signed distance to ∂Ω. This is an unitary extension
of the unitary normal vector n at ∂Ω which is smooth in the vicinity of ∂Ω. This extension furnishes
a symmetric Jacobian Dn that satisfies Dnn = 0 on ∂Ω. The direction h will be supposed to be in
C2(Rd,Rd) or in C2(∂Ω,Rd).
Proposition 3 The material derivative ṅ of the normal vector n at Ω in the direction of a vector
field h ∈ C1(Rd,Rd) is given by
ṅ = −∇τ (h.n) +Dnhτ ,
where hτ = h− h.n n.
Concerning its shape derivative defined as n′ = (∂tnt)|t=0 where nt is any smooth unitary extension
of n to ∂Ωt, we obtain.
Proposition 4 The shape boundary n′ in the direction of h is given by
n′ = −∇τ (h.n).
Derivative of the tangential gradient. For f ∈ H3(∂Ω), we compute the material derivative of
∇τf . We first compute the difference
∇τf −∇ḟ .
Proposition 5 For all functions f ∈ C2(R3) and directions h ∈ C2(∂Ω,R3), one has
˙∇τf = ∇ḟ + (D
2fh)τ −∇f.n ṅ−∇f.ṅ n
Proof of Proposition 5. We differentiate ∇f and ∇f.n n and obtain
∇̇f = ∇f ′ +D2fh
while
˙∇f.n n = ∇f.ṅ n+∇f.n ṅ+∇f ′.n n+ (D2fh).n n.
The two former equations give the desired result. �
Derivative of the Laplace-Beltrami operator. Now, we want to compute the material derivative
∆τf . We begin to study how to transport the Laplace-Beltrami operator when one works on ∂Ωt.
Let ∆τ,t denote the Laplace-Beltrami operator on the manifold ∂Ωt. To compute the derivative of a
Laplace-Beltrami operator, we need the following proposition that we quote from [20].
Proposition 6 Let f ∈ H5/2(Rd), then
∆τ,tf
◦ Tt γτ (t)
φ = −
∇(f ◦ Tt)− (B(t) n).∇(f ◦ Tt)
.∇φ, ∀φ ∈ D(Rd).
In the former proposition, we set
γ(t) = detDTt,
γτ (t) = γ(t)‖(DT
T .n‖Rd ,
B(t) =
D(T−1t )(D(T
‖(DT−1t )
T .n‖2
C(t) = γτ (t)D(T
t )(D(Tt)
−1)T .
. (22)
A straightforward computation gives
γ′(0) = divτ (h) ,
γ′τ (0) = divτ (h) = divτ (hτ ) +Hhn,
B′(0) = 2(Dhn).nI − (Dh+ (Dh)T ),
C′(0) = divτ (h) I − (Dh+ (Dh)
. (23)
Theorem 6 Let f ∈ D(Rd). The material derivative of ∆τf in the direction h is given by
∆τf = ∆τ ḟ+∇τf.∇τ
divτ (hτ )
+∇τ (Hhn).∇τf − divτ
Dh+ (Dh)T
Proof of Theorem 6 : Formula (24) is shown in a weak sense. For each test function φ ∈ C∞(∂Ω),
there exists an extension φ̃ ∈ D(Rd) such that ∂nφ̃ = 0; this can be done by extending φ as a constant
along the orbits of the gradient of the signed distance function to ∂Ω and the use of a cut-off function.
For f ∈ D(Rd), we set
A(t) =
(∆τ,tf) ◦ Tt −∆τf
γτ (t) φ.
After an integration by parts on ∂Ω, we obtain:
A(t) =
1− γτ (t)
(∆τ,tf) ◦ Tt φ+
(∆τ,tf) ◦ Ttφ+
∇τf.∇τφ
1− γτ (t)
(∆τ,tf) ◦ Tt φ
∇τf − C(t)∇ (f ◦ Tt)
.∇φ̃+
(B(t)n.∇(f ◦ Tt)
C(t) n.∇φ̃
Since ∂nφ̃ = 0 and C(0) = I, we get
A(t) =
1− γτ (t)
(∆τ,tf) ◦ Tt φ+
∇τ (f − f ◦ Tt)
.∇τ φ̃+
C(0)− C(t)
∇(f ◦ Tt).∇τ φ̃.
When t→ 0, it then comes
∆τfφ = −
γ′τ (t)∆τ fφ+∇τ ḟ .∇τφ+
C′(0).∇f
.∇τφ,
∆τ ḟ − divτ (h)∆τf
Dh+ (Dh)T − divτ (h) I
∇f.∇τφ,
∆τ ḟ − divτ (h)∆τf + divτ
divτ (h)∇τf
− divτ
Dh+ (Dh)T
Expanding the double divergence term, we obtain:
∆τf = ∆τ ḟ +∇τf.∇τdivτ (h)− divτ
Dh+ (Dh)T
In order to explicit these derivatives, we let appear the curvatures of ∂Ω by means of
∇τf.∇τdivτ (h) = ∇τf.∇τ
divτ (hτ ) +Hhn
and this ends the proof of the theorem (24). �
3 Existence of the second order derivative of the state. Proof
of Theorem 2.
The section is devoted to prove Theorem 2. We follow the usual strategy to derive existence in shape
optimization. In section 3.2, we will write the weak formulation of the problem, then transport it on
the reference domain, pass to the limit and obtain existence of the material derivative. In a second
time, we will seek a boundary value problem solved by the material derivative. This will provide
a characterization of the second order shape derivative. Two strategies, that we will detail, are
possible: the first one explored in section 3.3 consists in working on the variational formulation while
the second one uses the tangential differential calculus by differentiating the boundary conditions.
This last approach will be presented in section 3.4. The computations that will be made in subsections
3.3 and 3.4 require some regularity of the traces of the state ud on the interface of discontinuity ∂ω.
For the sake of readability, we postponed in subsection 3.5 all the needed justifications.
3.1 Preliminary results.
In the sequel, we will use some technical formulae. To preserve the readability of the proof of the
main result, we state them in this paragraph. The tools needed for proving these results can be found
in [20]. Given a smooth vector field h, we denote
Ah = Dh+Dh
T − div (h) I
We begin with the following formula.
Lemma 1 It holds:
∇u.Ah∇v = ∇(h.∇u).∇v +∇(h.∇v)∇u − div
(∇u.∇v)h
. (25)
Given two smooth vector fields h1 and h2, we set
A = Dh2Ah1 +Ah1Dh2
T −Ah1div (h2)− (Ah1)
′(h2), (26)
b = (h2.∇u)Ah1∇v + (h2.∇v)Ah1∇u − ((Ah1∇u).∇v)h2.
Here, the notation (Ah1)
′(h2) stands for the matrix defined by its elements
((Ah1)
′(h2))k,l = ∇(((Ah1 )
′)k,l).h2
Lemma 2 One has:
∇u.A∇v = div (b)− (h2.∇u)div
(Ah1∇v)
− (h2.∇v)div
(Ah1∇u)
. (27)
We need the following crucial result
Lemma 3 If u is harmonic then
Ah1∇u
= ∆(h1.∇u). (28)
Proof of Lemma 3 For any harmonic function u in Ω and for every test function φ ∈ D(Ω), we
can write ∫
∇u̇∇φ =
Ah∇u∇φ
then ∫
∆u̇ φ =
div (Ah∇u) φ
Since u̇ = u′ + h.∇u and since u′ is harmonic in Ω, we obtain the result. �
3.2 Proof of existence of the second order derivative.
We follow Hettlich and Rundell [8] and Simon [19] to define the second order derivative of an op-
erator with respect to a domain. We compute the second derivative by considering two admissible
deformations h1,h2 ∈ H that will describe the small variations of ∂ω. Simon shows that the second
derivative F ′′(∂ω;h1,h2) of F (∂ω) is defined as a bounded bilinear operator satisfying
F ′′(∂ω;h1,h2) =
F ′(∂ω;h1)
h2 − F
′(∂ω;Dh1 h2) (29)
For more details, the reader can consult the appendix in page 613 of [8].
Let us begin the proof. Let h1,h2 ∈ H be two vector fields. The direction h1 being fixed, we
consider u̇1,h2 the variation of u̇1 with respect to the direction h2. We recall from [1] that the material
derivative u̇1 of u in the direction h1 satisfies
∀v ∈ H10 (Ω),
σ∇u̇1.∇v =
σ∇u.Ah1∇v.
Let φ2 : Ω 7→ Ω be the diffeomorphism defined by φ2(x) = x+h2(x) and we set ψ2 = φ
2 . Setting
ωh2 =
x+ h2(x), x ∈ ω
, Ωh2 =
x+ h2(x), x ∈ Ω
= Ω and σh2 = σ ◦ φ2, we get
σh2∇u̇1,h2 .∇v =
σh2∇uh2 .Ah1∇v (30)
where uh2 is the solution of the original problem with ωh2 instead of ω. Making the change of
variables x = φ2(X), we get the integral identity on the fixed domain Ω :
σ∇˜̇u1,h2 .
Dψ2(Dψ2)
T det(Dφ2)
σ∇ũh2 .
Dψ2Ãh1(Dψ2)
T det(Dφ2)
∇v (31)
with the notations ũ = u ◦φ2 and Ãh1 = Ah1 ◦φ2. Since the material derivative u̇1 of u with respect
to the direction h1 satisfies ∫
σ∇u̇1.∇v =
σ∇u.Ah1∇v,
the difference of (30) and (31) gives
˜̇u1,h2 − u̇1
.∇v =
σ∇˜̇u1,h2 .
I −Dψ2(Dψ2)
T det(Dφ2)
σ∇ũh2 .
Dψ2Ãh1(Dψ2)
T det(Dφ2)−Ah1
(∇ũh2 −∇u).Ah1∇v.
We quote from [13] and [8] the following asymptotic formulae
‖ − div (hi) ‖∞ = O(‖hi‖
‖Dψi(Dψi)
T det(Dφi)− I +Ahi‖∞ = O(‖hi‖
‖Dψ2Ãh1(Dψ2)
T det(Dφ2)−Ah1 +Dh2Ah1 +Ah1(Dh2)
T − div (h2)Ah1 − (Ah1)
′(h2)‖∞ = O(‖h2‖
Making the adequate substitutions, we easily check that the material derivative of u̇1 with respect
to h2 exists. This derivative, denoted by ü1, satisfies
σ∇ü1.∇v dx =
∇u̇1.Ah2∇v +∇u̇2.Ah1∇v −∇u.A∇v
. (32)
where A is defined in (26).
3.3 Derivation of (8) from the weak formulation.
We want to make explicit the problem solved by (u′)′. To achieve this, we should write the right
hand side
∇u̇1.Ah2∇v +∇u̇2.Ah1∇v −∇u.A∇v
as the sum of an integral with ∇v in factor and an integral of a divergence to identify the jump
conditions on ∂ω. To that end, we will use algebraic identities that involve second order derivatives
of u, u̇i and of the test function v ∈ D(Ω). Using Lemma 1, we obtain:
σ∇u̇1.Ah2∇v =
∇(h2.∇u̇1).∇v +∇(h2.∇v)∇u̇1 − div
(∇u̇1.∇v)h2
σ∇u̇2.Ah1∇v =
∇(h1.∇u̇2).∇v +∇(h1.∇v)∇u̇2 − div
(∇u̇2.∇v)h1
Concerning the remaining terms, we use Lemma 2 to get
σ∇u.A∇v =
σ div
(h2.∇u)Ah1∇v + (h2.∇v)Ah1∇u− (Ah1∇u.∇v)h2
(h2.∇u)div
Ah1∇v
+ (h2.∇v)div
Ah1∇u
We apply Lemma 3 and gather the expressions obtained for F .
∇ (h1.∇u̇2 + h2.∇u̇1) .∇v +∇(h2.∇v).∇u̇1 +∇(h1.∇v).∇u̇2
σ div
(Ah1∇u.∇v −∇u̇1.∇v)h2 − (∇u̇2.∇v)h1
(h2.∇v)∆(h1.∇u)− div
(h2.∇v)Ah1∇u
−∇(h2.∇u).Ah1∇v
Using (25), we remove the dependency on Ah1∇v:
∇(h2.∇u).Ah1∇v = ∇(h1.∇(h2.∇u)).∇v +∇(h1.∇v)∇(h2.∇u)− div
(∇(h2.∇u).∇v)h1
Therefore, we write F = F1 + F2 where
∇ (h1.∇u̇2 + h2.∇u̇1)−∇(h1.∇(h2.∇u))
.∇v, (34)
∇(h1.∇v).∇(u̇2 − h2.∇u) +∇(h2.∇v).∇u̇1 + (h2.∇v)∆(h1.∇u)
σ div
(Ah1∇u.∇v −∇u̇1.∇v)h2 +
∇(h2.∇u).∇v −∇u̇2.∇v
h1 − (h2.∇v)Ah1∇u
The connection between second order material and shape derivatives is given by:
ü1 = (u
2 + h1.∇u̇2 + h2.∇u̇1 − h1.∇(h2.∇u),
incorporating this expression in (34), we rewrite (32) as:
∀v ∈ H10 (Ω),
σ∇(u′1)
2.∇v = F2. (35)
Testing it against v ∈ D(Ω \ ∂ω), we get ∆(u′1)
2 = 0 in Ω \ ω and in ω. We now deduce the jump
conditions for (u′1)
2. To obtain the jump of the potential, we simply write that ü1 ∈ H
0(Ω), hence
[ü1] = 0 on ∂ω and then
[(u′1)
2] = −h1.∇u
2 − h2.∇u̇1.
To express the jump of the flux, we then apply the Gauss formula in (35) to get
[σ∂n(u
2]v = F2. (36)
The second term F2 contains all the jumps of the flux on the interface ∂ω.
A simplified expression of F2. To get a simplified formula for F2 under a boundary integral,
some lengthy but straightforward calculations are needed. We summarize the result by means of the
following lemma
Lemma 4 One has:
2h2,nh1,nDn [σ∇τu]− h2,nn.∇h1,n [σ∇τu] + h2,nh1τ .Dnn [σ∇τu]
h1τ .∇τ (h2,n) [σ∇τu]− h1,nh2,nH [σ∇τu]
+ divτ
Proof of lemma First, write :
σ∇(h1.∇v).∇(u̇2 − h2.∇u) = σ1
∇(h1.∇v).∇u
2 + σ2
∇(h1.∇v).∇u
[σ∂nu
2](h1.∇v)
Note that the normal vector is oriented from ω to Ω \ ω. In the same spirit, we write
∇(h2.∇v).∇u̇1 + (h2.∇v)∆(h1.∇u) = ∇(h2.∇v).∇(u̇1 − h1.∇u) + div
(h2.∇v).∇(h1.∇u)
By a argument of symmetry, we then can write:
σ∇(h2.∇v).∇(u̇1 − h1.∇u) = −
[σ∂nu
1](h2.∇v).
To drop the dependency in Ah1 , we use (25) and get after expansion:
Ah1∇u.∇v
= div
∇(h1.∇v).∇u +∇(h1.∇u)∇v
− div
(∇u.∇v)h1
(h2.∇v)Ah1∇u
= ∇(h2.∇v).Ah1∇u + (h2.∇v)div
Ah1∇u
= ∇(h1.∇(h2.∇v)).∇u +∇(h1.∇u)∇(h2.∇v) + (h2.∇v)∆(h1.∇u)
− div
∇(h2.∇v).∇u
h1.∇(h2.∇v)
.∇u+ div
(h2.∇v)∇(h1.∇u)−
∇(h2.∇v).∇u
After integrating by parts, we conclude thanks to the state equation and obtain
h1.∇(h2.∇v)
.∇u = −
h1.∇(h2.∇v)
div (σ∇u) = 0
We substitute the shape derivative u′ to the material one u̇:
F2 = −
[σ∂nu
1](h2.∇v) + [σ∂nu
2](h1.∇v)−
σ div
(∇u.∇v)h1
σ div
∇(h1.∇v).∇u)h2 + (∇(h2.∇v).∇u
(∇u′2.∇v)h1 + (∇u
1.∇v)h2
First, we use the continuity of the flux on ∂ω, then we integrate by parts on ∂ω and finally we
incorporate the expressions of the jumps of the shape derivatives u′ to obtain
∇(h2.∇v).∇u
σ∇u.∇(h2.∇v)
h1,n = −
[σ∇τu]h1,n∇τ (h2.∇v)
[σ∇τu]h1,n
h2.∇v =
h2.∇v.
This leads to a simplified expression for F2:
F2 = −
σ div
(∇u.∇v)h1
(∇u′1.∇v)h2 + (∇u
2.∇v)h1
Let us study each term of this sum. Using Gauss formula and integrating by parts on the manifold
∂ω, we obtain
σ div
∇u′1.∇v)h2
σ∇u′1.∇v
∂nv −
∂nv +
By symmetry, we also get:
σ div
∇u′2.∇v)h1
∂nv +
We now turn to the term with a double divergence. We first write it as a boundary integral thanks
to Gauss formula as
σ div
(∇u.∇v)h1
h2,ndiv
σ(∇u.∇v)
then, we use (14) to introduce the tangential operators
σ div
(∇u.∇v)h1
h2,ndivτ
σ(∇u.∇v)
h2,nD(h1
σ(∇u.∇v)
)n.n.
We study each of these terms. We start with the one involving tangential derivatives: we expand the
tangential divergence to incorporate the jump relation for the state u.
σ(∇u.∇v)
= divτ (h1)
σ(∇u.∇v)
+ h1.∇τ [σ∇u.∇v]
= divτ (h1) [σ∇τu] .∇τv + h1.∇τ [σ∇τu.∇τv] .
Then, the first term becomes:
h2,ndivτ
σ(∇u.∇v)
h2,ndivτ (h1) [σ∇τu] .∇τv +
h2,nh1.∇τ [σ∇τu∇τv] .
We use the integration by parts formula (18) to get:
h2,ndivτ
σ(∇u.∇v)
h1,nh2,nH [σ∇τu] .∇τv − divτ
divτ (h1)h2,n [σ∇τu]
v − divτ
h1h2,n
[σ∇τu] .∇τv
h1h2,n
− divτ (h1)h2,n − h1,nh2,nH
[σ∇τu]
Expanding
h1h2,n
[σ∇τu]
v = divτ
divτ (h1) h2,n [σ∇τu] + h1.∇τ (h2,n) [σ∇τu]
= divτ
divτ (h1) h2,n [σ∇τu]
v + divτ
h1τ∇τh2,n [σ∇τu]
we obtain the new expression:
h2,ndivτ
σ(∇u.∇v)
h1τ∇τh2,n − h1,nh2,nH
[σ∇τu]
v. (38)
Now, we consider the term involving normal components. We have
n.D(h1
σ∇(∇u.∇v)
)n = n.∇(h1,n [σ∇u.∇v])− [σ∇u.∇v]h1τ .Dnn
= n.∇(h1,n) [σ∇τu]∇τv + h1,nn.∇([σ∇u.∇v]).
Then, we get
h2,nD(h1
σ(∇u.∇v)
)n.n =
h2,nn.∇(h1,n) [σ∇τu] .∇τv + h2,nh1,nn.∇([σ∇u.∇v])
−divτ
h2,nn.∇(h1,n) [σ∇τu]
v + h2,nh1,nn.∇([σ∇u.∇v]).
A straightforward calculus leads to
n.∇([σ∇u.∇v]) = n.
σD2u∇v
+D2v [σ∇u]
∇τv +D
2v [σ∇τu]
= ∂nv
∇τv + n.D
2v [σ∇τu] .
where D2u is the Hessian matrix of u. From (20) and from the jump conditions for the state u, we
deduce that [
= − [σ∆τu] .
When one differentiates the relation expressing the continuity of the flux for the state along the
tangential direction ∇τv, one gets ([6], p 235):
0 = ∇[σ∂nu].∇τv = [σD
2u]∇τv.n+ [σ∇u].(Dn∇τv).
In the same spirit, it comes that
∇∂nv.[σ∇τu] = D
2v[σ∇τu].n+∇v.(Dn[σ∇τu]). (40)
Since Dn is a symmetric matrix and Dnn = 0, one checks ∇v.(Dn[σ∇τu]) = [σ∇u].(Dn∇τv).
n.∇([σ∇u.∇v])) = − [σ∆τu] ∂nv − 2Dn [σ∇τu] .∇τv + [σ∇τu]∇τ∂nv
We integrate this expression on ∂ω and obtain after some integration by parts:
h2,nh1,nn.∇([σ∇u.∇v])
h2,nh1,n [σ∆τu]∂nv +
h2,nh1,n [σ∇τu]∇τ∂nv − 2
h2,nh1,nDn [σ∇τu] .∇τv,
h2,nh1,n [σ∆τu] + divτ
h2,nh1,n [σ∇τu]
∂nv + 2
h2,nh1,nDn [σ∇τu]
Hence
h2,nD(h1
σ(∇u.∇v)
)n.n = −
h2,nh1,n [σ∆τu] + divτ
h2,nh1,n [σ∇τu]
2h2,nh1,nDn [σ∇τu]− h2,nn.∇τh1,n [σ∇τu]
(∇u.∇v)h1
h2,nh1,n [σ∆τu] + divτ
h2,nh1,n [σ∇τu]
2h2,nh1,nDn [σ∇τu]− h2,nn.∇h1,n [σ∇τu]
h1τ .∇τ (h2,n) [σ∇τu]− h1,nh2,nH [σ∇τu]
Gathering all the terms, we write F2 as:
2h2,nh1,nDn [σ∇τu] +
h1τ .∇τ (h2,n)− h2,nn.∇h1,n − h1,nh2,nH
[σ∇τu]
+ divτ
h2,nh1,n [σ∆τu] + divτ
h2,nh1,n [σ∇τu]
h1,ndivτ
h2,n [σ∇τu]
+ h2,ndivτ
h1,n [σ∇τu]
We end the proof after expanding the tangential divergence of the last term of F2. �
Let us return to the weak formulation (36) of the derivative. By identification, we get
[σ∂n(u
2] = divτ
+ divτ
− divτ
h2,nh1,n(2Dn−HI) [σ∇τu]
− divτ
h1τ .∇τ (h2,n) [σ∇τu]− h2,nn.∇h1,n [σ∇τu] + h2,nh1τ .Dnn [σ∇τu]
It remains to compute the jump of the flux for the second order derivative. Since
u′′1,2 = (u
2 − u
Dh1 h2
where u′Dh1 h2 is the first shape derivative of u in the direction of the vector field Dh1 h2. Thanks
to (6), we can write the jump under the form
[σ∂nu
1,2] = [σ∂n(u
2]− [σ∂nu
Dh1 h2
] = [σ∂n(u
2]− divτ
Dh1h2.n[σ∇τu]
. (42)
Let us split the field h2 in two parts: Dh1 h2.n = h2,nn.Dh1 n +Dh1 h2τ .n. In the spirit of (40),
we obtain
Dh1h2τ .n = ∇τh1,n.h2τ − h1τ .Dnh2τ . (43)
Thanks to (39), the jump [σ∂nu
Dh1 h2
] then can be written under the form
[σ∂nu
Dh1 h2
] = divτ
(h2,nn.∇h1,n +∇τh1,n.h2τ − h1τ .Dnh2τ )[σ∇τu]
Gathering all the terms, simplifications occur and we get:
[σ∂nu
1,2] =divτ
+ h1,n
− divτ
(h1τ .∇τh2,n +∇τh1,n.h2τ ) [σ∇τu]
− divτ
h2,nh1,n(2Dn−HI) [σ∇τu]
+ divτ
h1τ .Dnh2τ )[σ∇τu]
To get the jumps of the potential, we use (41) and obtain
u′′1,2
(u′1)
u′Dh1 h2
= −h1.
− h2. [∇u̇1]−
u′Dh1 h2
= −h2,n
− h1,n
− h2τ .
− h1τ .
− h2,nn.
∇(h1.∇u)
− h1,nn.
∇(h2.∇u)
u′Dh1 h2
Thanks to the jump of the potential for the first order shape derivative given in (6), it comes that
h2τ .
= −h2τ .
∇(h1.∇u)
and h1τ .
= −h1τ .
∇(h2.∇u)
and then:
u′′1,2
= −h2,n
− h1,n
− h2,nn.
∇(h1.∇u)
+ h1τ .
∇(h2.∇u)
u′Dh1 h2
Computing the other jumps that appeared in the former expression, we get
∇(h2.∇u)
= (Dh2)
T [∇u] +
h1τ .
∇(h2.∇u)
= n.Dh2h1τ [∂nu] + h2,nh1τ .
n+ h1τ .
h2τ .
h2,nn.
∇(h1.∇u)
= h2,n [∂nu]n.Dh1n+ h2,nh1,nn.
n+ h2,nn.
h1τ .
u′Dh1 h2
= −Dh1 h2.n [∂nu] = −
h2,nn.Dh1n+ n.Dh1h2τ
[∂nu]
With the help of formula (43), we obtain:
−h2,nn.
∇(h1.∇u)
+ h1τ .
∇(h2.∇u)
u′Dh1 h2
∇τh1,n.h2τ +∇τh2,n.h1τ
[∂nu]
− 2h1τ .Dnh2τ [∂nu] + h1τ .
h2τ − h2,nh1,nn.
n = − [∆τu]−H [∂nu] = −H [∂nu] ,
h1τ .
h2τ = h1τ .D([∇u])h2τ = h1τ .D([∂nu]n)h2τ = [∂nu]h1τ .Dnh2τ .
Finally, we gather the results of these computations to write
u′′1,2
+ h1,n
∇τh1,n.h2τ +∇τh2,n.h1τ
[∂nu]
h2,nh1,nH − h1τ .Dnh2τ
[∂nu])
3.4 How to recover (8) by formal differentiation of the boundary condi-
tions.
The aim of this section is to retrieve the expression of the flux jump [σ∂nu
′′] by computing the normal
derivatives of each of the expressions
[σ∇u′].n and
h1,n[σ∇τu]
. Since
[σ∇u′].n = divτ
h1,n[σ∇τu]
= h1,n[σ∆τu] +∇τh1,n.[σ∇τu],
then, we get
[σ∇u′].n = ˙h1,n[σ∆τu] + h1,n
[σ∆τu] +
˙∇τh1,n.[σ∇τu] +∇τh1,n.
[σ∇τu]. (46)
In order to avoid lengthy computations, we shall concentrate on each normal derivative appearing in
the above formula. Some of the results are straightforward and their proof will be left to the reader.
Combining propositions (3) and (5), we conclude that
˙∇τh1,n = −∇τ (h1.∇τh2,n) + (D
2h1,n.h2)τ −∇h1,n.ṅ n−∇h1,n.n ṅ.
In the same manner, we also get
[σ∇τu] = [σ∇τu
2] + ([σD
2u].h2)τ − [σ∇τu].ṅ n− [σ∇τu]n ṅ.
Hence, we can write
h1.n = h2.∇h,n −∇τh2,n.h1τ .
It remains to simplify the terms A = (D2u.h2)τ .∇τh1,n and B = [σ∇τu].(D
2h1,n.h2)τ . We obtain:
A = −[σ∇τu].(Dn∇τh1,n)h2,n + [σ∆τu]∇τh1,n.h2τ ,
B = (D2h2,n.h2τ ).[σ∇τu] +∇τ (∂nh1,n).[σ∇τu]h2,n − [σ∇τu].(Dn∇τh1,n)h2,n.
We tackle the computation of (∂nu
′)′. For the sake of clearness, we subdivide the work in several
steps.
First step. We compute
h1,n[σ∇τu]
. We expand:
h1,n[σ∇τu]
h1,n[σ∆τu] +
∇τh1,n.[σ∇τu],
h1,n[σ∆τu] + h1,n
[σ∆τu] +
∇τh1,n.[σ∇τu] +∇τh1,n.
[σ∇τu].
Hence, after substitution, one gets
h1,n[σ∇τu]
= divτ
h1,n[σ∇τu
2] + (h2,n∂nh1,n −∇τh2,n.h1τ )[σ∇τu]
+ 2[σ∆τu]∇τh1,n.h2τ − ∂nh1,[σ∇τu].(Dnh2τ ) + [σ∇τu].(D
2h1,n.h2τ )
− 2h2,n[σ∇τu].(Dn∇τh1,n) + h1,n
σ∆τu− [σ∆τu
. (47)
Second step. We compute
[σ∂nu
1]. From the expression of ṅ, we get after some straightforward
computations:
[σ∂nu
1] = [σ∂n(u
2] + ([σD
2u′1]h2).n+ [σ∇τu
1].(Dnh2τ −∇τh2,n). (48)
Third step. We compute σ∂n(u
2. From the jump condition on the flux of the derivative (6) and
(47) and (48), we obtain:
[σ∂(u′1)
2] = divτ
h1,n[σ∇τu
h2,n∂nh1,n −∇τh2,n.h1τ
[σ∇τu]
+ 2∇τh1,n.h2τ [σ∆τu]
−([σD2u′1]h2).n+ [σ∇τu
∇τh2,n −Dnh2τ
− ∂nh1,n[σ∇τu].(Dn h2τ )
+(D2h1,n h2τ ).[σ∇τu]− h2,n(Dn [σ∇τu])£.∇τh1,n.
Taking account of the following calculation,
−([σD2u′1]h2).n+ [σ∇τu
1].∇τh2, n = −
h2,n[σD
2u′1]n+ [σD
2u′1]h2τ
.n+ [σ∇τu
1].∇τh2, n,
= h2,n
[σ∆τu
1] +H [σ∂nu
+ [σu′1].∇τh2,n − ([σD
2u′1]h2τ ).n,
= divτ
h2,n[σ∇τu
+Hh2,n[σ∂nu
1]− ([σD
2u′1]h2τ ).n;
it comes
[σ∂n(u
2] = divτ
h1,n[σ∇τu
2] + h2,n[σ∇τu
h2,n∂nh1,n −∇τh2,n.h1τ
[σ∇τu]
+2[σ∆τu]∇τh1,n.h2τ +Hh2,n[σ∂nu
1]− ([σD
2u′1]h2τ ).n
[σ∇τu
1] + ∂nh1,n[σ∇τu]
.(Dnh2τ ) + (D
2h1,n h2τ ).[σ∇τu]
−2h2,n∇τh1,n.(Dn [σ∇τu]) + h1,n
[σ∆τu]− [σ∆τu
. (49)
This formula remains hard to handle. To get a more convenient one, we decide to derive tangentially
to the direction h2 the boundary identity
[σ∂nu
1] = h1,n[σ∆τu] +∇τh1,n.[σ∇τu].
This leads to:
([σD2u′1]h2τ ).n+(Dnh2τ ).[σ∇τu
1] = ∇τh1,n.h2τ [σ∆τu] + h1,n∇τ [σ∆τu].h2τ
+ (D2h1,n h2τ ).[σ∇τu]− ∂nh1,n[σ∇τu].(Dnh2τ ) + [σ∆τu]h2τ .∇τh1,n.
From (24) and subtracting (50) from (49), we can write
[σ∂n(u
2] = divτ
h1,n[σ∇τu
2] + h2,n[σ∇τu
h2,n∂nh1,n −∇τh2,n.h1τ
[σ∇τu]
+divτ
h1,nh2,n(HI − 2Dn).[σ∇τu]
− h1,n
∇τ [σ∆τu].h2τ +∆τ [σ∇τu].h2τ
+h1,n
∇τdivτ (h2τ ) .[σ∇τu]− divτ
Dh2 + (Dh2)
[σ∇τu]
From (24), we obtain
[σ∆τu] = [σ∆τ u̇] +∇τdivτ (h2τ ) .[σ∇τu] +∇τ (Hh2,n).[σ∇τu]
− divτ
Dh2 + (Dh2)
[σ∇τu]
and using the relation between the material and shape derivative, we get
[σ∆τu] = [σ∆τu
′] +∇
[σ∆τu]
.h2 and [σ∆τ u̇] = [σ∆τu
′] + ∆τ
[σ∇u].h2
Injecting these relations in (51) and applying them for h2τ , we get
∆τ ([σ∇τu].h2τ ) +∇τdivτ (h2τ ) .[σ∇τu] = ∇τ [σ∆τu].h2τ + divτ
Dh2 + (Dh2)
[σ∇τu]
This last fact allows us to conclude.
3.5 Justification of the formal computations.
We have to justify rigorously that the right-hand sides of (6),(7),(8) make sense. They involve
tangential derivatives of un and ud along the interface ∂ω up to the order three. The existence of
these derivatives is not clear a priori since the gradient of the solution has a discontinuity along this
interface. Our first aim is to precise the tangential regularity along the interface ∂ω of the solution
u of (1) with either Dirichlet or Neumann boundary conditions.
We should access to the trace of u on the interface ∂ω. Any numerical discretization needs also
to compute the state, its derivatives with respect to the shape and the normal derivatives along the
interface ∂ω. To that end, we introduce for any α ∈ H1/2(∂ω) and β ∈ H−1/2(∂ω) the following
boundary value problems
∆v = 0 in Ω \ ω and in ω,
[v] = α on ∂ω,
[σ∂nv] = β on ∂ω,
v = f1 on ∂Ω.
and (N)
∆v = 0 in Ω \ ω and in ω,
[v] = α on ∂ω,
[σ∂nv] = β on ∂ω,
∂nv = g1 on ∂Ω,
where (f1, g1) ∈ H
1/2(∂Ω)×H−1/2(∂Ω). Note that for α = 0, β = 0 and (f1, g1) = (f, g) then (ud)
and un solve respectively (D) and (N); furthermore the choice of
hn∂nu
+ and β = [σ] divτ (hn∇τu) (53)
leads to (6) and (7) when we take (f1, g) = (0, 0).
Existence of solutions to (D) and (N). To study these problems, we use the integral rep-
resentation in terms of layer potentials. In a first step, we recall some definitions. The Newtonian
potential Γ is defined as:
Γ(x, y) =
ln(|x− y|) if n = 2,
|x− y|
if n = 3.
The integral equations applying to direct problem will be obtained from a study of the classical
single- and double-layer potentials. We begin to introduce the following operators
S∂Ω∂ω : u 7→ S∂Ω∂ωu(x) :=
Γ(x, y)u(y) dσ(y);
S∂ω∂Ω : u 7→ S∂ω∂Ωu(x) :=
Γ(x, y)u(y) dσ(y);
K∂Ω∂ω : u 7→ K∂Ω∂ωu(x) :=
∂nΓ(x, y)u(y) dσ(y) ;
K∂ω∂Ω : u 7→ K∂ω∂Ωu(x) :=
∂nΓ(x, y)u(y) dσ(y)
Note that all these operators have a smooth kernel since the boundaries ∂ω and ∂Ω are assumed to
have no common point. We also denote
SΩ : u 7→ SΩu(x) :=
Γ(x, y)u(y) dσ(y);
KΩ : u 7→ KΩu(x) :=
∂nΓ(x, y)u(y) dσ(y);
Sω : u 7→ Sωu(x) :=
Γ(x, y)u(y) dσ(y);
Kω : u 7→ Kωu(x) :=
∂nΓ(x, y)u(y) dσ(y).
We now obtain some systems of integral equations to compute the state function and their shape
derivatives. Since v is harmonic in Ω \ ω and for all x ∈ ∂Ω ∪ ∂ω, it has the classical boundary
representation:
v(x) =
∂nΓ(x, y)v(y) −
∂nΓ(x, y)v(y)−
Γ(x, y)∂nv(y) +
Γ(x, y)∂nv(y). (54)
Similarly since v harmonic in ω, for all x ∈ ∂ω we can write
v(x) =
∂nΓ(x, y)v(y) −
Γ(x, y)∂nv(y). (55)
Let us denote by vd the solution of the boundary values problem (D) in (52). Let us show how
to compute their restrictions and also their normal vector derivatives on the boundaries. Incorpo-
rating the jump conditions, a straightforward computation leads to the following boundary integral
equations
I + µKω
σ2 + σ1
S∂Ω∂ω
µK∂ω∂Ω
σ2 + σ1
(v+d )|∂ω
(∂nvd)|∂Ω
σ1 + σ2

I −Kω
−σ2K∂ω∂Ω S∂ω∂Ω

σ1 + σ2

K∂Ω∂ωf1

where µ = [σ]/(σ1 + σ2). Thanks to (55), the quantity (∂nvd)
+ is then given by
Sω(∂nvd)
I +Kω
v+d (x)|∂ω − α
Concerning vn, the solution of the Neumann problem (N) in (52), the same kind of computations
gives

I + µKω −
σ2 + σ1
K∂Ω∂ω
µK∂ω∂Ω −
σ2 + σ1
I +KΩ

(v+n )|∂ω
(vn)|∂Ω
σ1 + σ2

I −Kω
−σ2K∂ω∂Ω S∂ω∂Ω

σ1 + σ2
S∂Ω∂ωg1
Finally, the computation of (∂nvn)
is given by
Sω(∂nvn)
I +Kω
v+n (x)|∂ω − α
Concerning the well-posedness of (56), we can state the following result.
Theorem 7 The linear system of integral equation (56) has an unique solution in H1/2(∂ω) ×
H−1/2(∂Ω).
Proof of Theorem 7 Let A be the matricial operator defined on H1/2(∂ω)×H−1/2(∂Ω) as
I + µKω
σ2 + σ1
S∂Ω∂ω
µK∂ω∂Ω
σ1 + σ2
The main argument of the proof is based on the Fredholm alternative. In a first step, we have to show
that the adjoint operator A∗ is injective. Since the boundaries are bounded, the adjoint operator A∗
defined on H−1/2(∂ω)×H1/2(∂Ω) can be written under the form
I + µK∗ω µK
σ2 + σ1
S∂ω∂Ω
σ2 + σ1
. (59)
Let (u, v) ∈ H−1/2(∂ω)×H1/2(∂Ω) be in the kernel of A∗. Consider the potential W defined for each
x ∈ Rd by
W (x) =
σ2 + σ1
Γ(x, y)u(y) +
Γ(x, y)v(y)
. (60)
In a first step, we show that W = 0. The function W satisfies ∆W = 0 in Rd \ (∂ω ∪ ∂Ω) by
construction. We check that W |∂Ω = 0 from the equation corresponding to the second line of A
By the properties of the single layer potential, [W ] = 0 on ∂ω. Furthermore, it holds [σ∂nW ] = 0 on
∂ω. Indeed, we can have ([11])
σ1 + σ2
u+K∗∂Ω∂ωv
σ1 + σ2
v +K∗∂Ω∂ωv
hence,
σ1∂nW
+ − σ2∂nW
− = σ1
I + µK∗ω)u+ µK
∂Ω∂ωv
This corresponds to the first line of A∗(u, v). Then, W solves the Laplace equation (1) with ho-
mogeneous Dirichlet boundary conditions. By the uniqueness of the solution, we get W = 0 in
In a second step, we deduce that u = v = 0. Since W = 0 in Ω, we see that [∂nW ] = 0 on ∂ω.
Since [∂nW ] = σ1u/(σ1 + σ2) on ∂ω , we deduce u = 0. From the second line of A
∗(u, v) = 0, we
see that SΩv = 0 on ∂Ω. Since the single layer potential operator SΩ : H
−1/2(∂Ω) 7→ H1/2(∂Ω) is an
isomorphism, v = 0 holds. The injectivity of A∗ is proved. Since 2A = I + C where C is a compact
operator, we conclude that A has a continuous inverse thanks to the Fredholm alternative. �
In a similar way, the problem (57) is well-posed under some additional assumptions. We define
the adequate space
♦ (∂Ω) =
φ ∈ H1/2(∂Ω) :
φ = 0
We can state the following result.
Theorem 8 If we impose the normalizing condition
then there exists one unique couple ((vn)|∂ω, (vn)|∂Ω) ∈ H
1/2(∂ω)×H
♦ (∂Ω) solution of (57) .
Proof of Theorem 8 Set

I + µKω −
σ2 + σ1
K∂Ω∂ω
µK∂ω∂Ω −
σ1 + σ2
I +KΩ

the operator defined on H1/2(∂ω)×H
♦ (∂Ω). The adjoint B
∗ can be written under the form
I + µK∗ω µK
σ1 + σ2
K∗∂ω∂Ω −
σ1 + σ2
I +K∗Ω
In a first step, we begin to show that B∗ is injective. Let (u, v) ∈ H1/2(∂ω) × H1/2(∂Ω) be in the
kernel of B∗. We introduce the potential
Z(x) = −
σ1 + σ2
Γ(x, y)u(y) +
Γ(x, y)v(y)
, x ∈ Rd.
We can see that Z is a harmonic function in Rd\(∂ω ∪ ∂Ω), satisfying ∂nZ = 0 on ∂Ω. By the
properties of the single layer potential, [Z] = 0 Furthermore, a straightforward calculation shows
that [σ∂nZ] = 0 on ∂ω. Hence, Z solves the boundary value problem
−div (σ∇Z) = 0 in Ω,
∂nZ = 0 on ∂Ω.
The function is therefore constant in Ω. Writing [∂nZ] = 0 on ∂ω, we get easily u = 0 and then
+K∗Ω)v = 0. Since the operator λI −K
Ω is one to one on H
♦ (∂Ω), we deduce that v = 0. We
conclude the proof thanks to the Fredholm alternative. �
Tangential regularity results. Let us consider now the particular case where both α and β
are the zero function and (f1, g1) = (f, g) where f and g are respectively the Dirichlet and Neumann
boundary data. To recover the tangential regularity of the solution u along ∂ω, we look at the first
line of (56) to deduce that
I + µKω
(ud)|∂ω = −
σ2 + σ1
S∂Ω∂ω∂nud|∂Ω +
σ1 + σ2
K∂Ω∂ωf ; (63)
Sω(∂nud)
I +Kω
u+d (x)|∂ω (64)
It is easy to deduce that (ud)|∂ω ∈ C
3,α(∂ω). Indeed, from (63) that we consider as an equation in
(ud)|∂ω with data f and (∂nud)|∂Ω = g, we see that (f, (∂nvd)|∂Ω) belongs to H
1/2(∂Ω)×H−1/2(∂Ω),
thanks to Theorem 7.
In order to give a sense to the jump conditions arising in (6),(7),(8), we need to work in space of
functions of higher regularity. We choose the framework of Hölder spaces. We quote [12] to precise
the behavior of the layer potentials on these spaces.
Theorem 9 (Kirsch [12])
1. If ∂ω is of class C2,α, 0 < α < 1 then the operators Sω and Kω map C
β(∂ω) continuously into
C1,β for all 0 < β ≤ α.
2. Let k ∈ N with k 6= 0. If ∂ω is of class Ck+1,α with 0 < α < 1, then the operators Sω and Kω
map Ck,β(∂ω) continuously into Ck+1,β(∂ω) for all 0 < β ≤ α.
3. Let k be an integer. If ∂ω is of class Ck+2,α, then K∗ω maps C
k,β continuously into Ck+1,β(∂ω)
for all 0 < β ≤ α.
We go back to the proof. Since the two boundaries have no intersection point and since ∂ω is of
class C4,α, it follows that the right hand side of the former equation is of class C3,α(∂ω). We then
conclude the solution of (63) will be of class C3,α since the operator 1/2I + µKω is an isomorphism
from C3,α(∂ω) into itself. With the same arguments, we show straightforwardly that (∂nun)
∈ C2,α.
About the regularity of the jumps of the second derivative. The equations giving the jump
conditions [u′d] and [∂u
d] show obviously that [u
d] and [∂nu
d] belong respectively to C
2,α(∂ω) and
C1,α(∂ω). Hence, it comes straightforwardly that [u′′d ] ∈ C
1,α. With the same arguments, we show
that [∂nu
d ] ∈ C
0,α(see [20] for more details) and then that all the formal computations to get the
equations describing the second derivative have a sense.
Remark 1 In a view of a numerical discretization of the state equation, one has to emphasize that
the choice of a finite elements method seems inappropriate: one should extract tangential derivative
of high order on the interface ∂ω. The obtained numerical accuracy is not sufficient to incorporate
the results in an optimization scheme. On the converse, the systems of boundary integral equations
(56) and (57) are well-suited for this kind of computation. Nevertheless, a discussion of adapted
schemes should be precise and is out of the scope of this manuscript.
3.6 Case of Neumann boundary conditions.
Since the admissible deformation fields have a support with no intersection points with the outer
boundary, it is a straightforward application of the preceding computations to show that un solution
to (1)-(4) is twice differentiable with respect to the shape. Furthermore, its second order derivative
u′′n belongs to H
1(Ω \ ω) ∪H1(ω) and solves


∆u′′n = 0 in ω \ ω and in ω,[
h1,nh2,nH − h1τ .Dnh2τ
[∂nun]−
h1,n[∂n(un)
2] + h2,n[∂n(un)
h1τ .∇h2,n + h2τ .∇h1,n
[∂nun] on ∂ω,[
= divτ
σ∇τ (un)
+ h1,n
σ∇τ (un)
+ h1τ .Dn.h2τ )[σ∇τun]
−divτ
(h1τ .∇τh2,n +∇τh1,n.h2τ + h2,nh1,n(2Dn−HI)) [σ∇τun]
on ∂ω,
n = 0 on ∂Ω;
where we use the notations of Theorem 2.
4 Second order derivatives for the criterion.
4.1 Proof of Theorem 4.
The differentiability of the objective is a direct application of Theorem 2. The computation we make
here is based on the relation
D2JKV (ω)(h1,h2) = D
DJKV (w)h1
h2 −DJKV (w)Dh1h2. (66)
To obtain (10), we compute in a first step the shape gradient in the direction h1. Then, in a second
step, we differentiate the obtained expression in the direction of h2. In the sequel, we adopt the
notation v = ud − un to obtain concise expressions.
DJKV (ω)h1 = σ1
|∇v|2h1
+ 2∇v.∇v′1 + σ2
|∇v|2h1
+ 2∇v.∇v′1
= σ1(A1 + 2B1) + σ2(A2 + 2B2),
where
|∇v|2h1
∇v.∇v′1
|∇v|2h1
∇v.∇v′1.
Now, we use the classical formulae to differentiate a domain integral to get
DA1(ω)h2 =
|∇v|2h1
+ 2div
∇v.∇v′2 h1
|∇v+|2h1
h2,n + 2∇v
+.∇(v′2)
+ h1,n;
DA2(ω)h2 =
|∇v−|2h1
h2,n + 2∇v
−.∇(v′2)
− h1,n.
The terms DBi, i = 1, 2 require more precisions. First, we write
DB1(ω)h2 =
∇v.∇v′1)h2
+∇v′1.∇v
2 +∇v.∇(v
∇v+.∇(v′1)
+h2,n + ∂nv
+((v′1)
+)′2 +
+(v′2)
+ + ∂n(v
+(v′1)
∂nv((un)
∂n(ud)
1(un)
2 + ∂n(ud)
2(un)
Note that we used the Green formula twice to keep the symmetry in h1 and h2. We also use the fact
that the derivatives (ud)
i are harmonic in Ω \ ω to transform the boundary integral on the exterior
boundary into an integral on the moving boundary. We obtain
DB1(ω)h2 = −
∇v+.∇(v′1)
+h2,n + ∂nv
+(((ud)
+ − v∂n(((un)
+((ud)
+ + ∂n(v
+((ud)
+ − ∂n(un)
+(v′2)
+ − ∂n((un)
+(v′1)
By the same methods, we get
DB2(ω)h2 =
∇v−.∇(v′1)
−h2,n + ∂nv
−((v′1)
−(v′2)
− + ∂n(v
−(v′1)
We regroup the different terms and after some straightforward computations, we obtain:
DJKV (ω)h1
(ω)h2 = −
σ|∇v|2h1
h1,n∇v
2 + h2,n∇v
1 + (ud)
2 − ∂n(un)
1 − ∂n(un)
σ∂n((un)
− σ1∂nv
((ud)
In order to compute D2JKV (ω)(h1,h2), the first order derivative of the Kohn-Vogelius objective is
needed. It can be written as follows:
DJKV (w)h = −
σ|∇v|2
hn + 2
σ∂n(un)
− σ1∂nv
Gathering (66),(41) and (42), we write the second derivative of the Kohn-Vogelius criterion as:
D2JKV (ω)(h1,h2) = −
σ|∇v|2h1
σ|∇v|2
(Dh1h2).n
1 + (ud)
2 − ∂n(un)
1 − ∂n(un)
h1,n∇v
2 + h2,n∇v
σ∂n(un)
− σ1∂nv
Let us give a more simplified version for the first term. We decompose the field h2 into normal vector
and tangential parts and we use (43). After some elementary computations, we obtain
σ|∇v|2h1
σ|∇v|2
(Dh1h2).n
σ|∇v|2
h1τ .∇τh2,n + h2τ .∇τh1,n − h2τ .Dnh1τ
σ|∇v|2
h1,nh2,n.
Finally, the second order derivative of the Kohn-Vogelius objective becomes:
D2JKV (ω)(h1,h2) =
σ|∇v|2
h1τ .∇τh2,n + h2τ .∇τh1,n − h2τ .Dnh1τ
σ|∇v|2
h1,nh2,n + 2
h1,n∇v
2 + h2,n∇v
1 + (ud)
2 − ∂n(un)
1 − ∂n(un)
σ∂n(un)
− σ1∂nv
4.2 Analysis of stability. Proof of Theorem 5
Now, we specify the domain ω that is assumed to be a critical shape for JKV . Moreover, we assume
that the additional condition ud = un holds. To emphasize that we deal with such a special domain,
we will denote it ω∗. The assumptions mean that the measurements are compatible and that ω∗ is
a global minimum of the criterion. From the necessary condition of order two at a minimum, the
shape Hessian is positive at such a point.
Let us notice that only the normal component of h appears. Let us also emphasize that there
is no hope to get h = 0 from the structure theorem for second order shape derivative ([6]). The
deformation field h appears in D2JKV (ω
∗)(h,h) only thought its normal component hn since ω
a critical point for JKV . This remark explains why we consider in the statement of Theorem 5 the
scalar Sobolev space corresponding to the normal components of the deformation field.
We now prove Theorem 5. From (67), we deduce
DJ2KV (ω
∗)[h, h] = −2
u′d∂nv
′ − ∂nu
= 2 [σ]
(u′+d − u
n )divτ (hn∇τud)−
d hn∂n(u
d − u
= 2 [σ]
u′+d − u
n , divτ (hn∇τud)
∂nudhn, ∂n
u′d − u
where 〈, 〉 denotes the duality between H1/2(∂ω∗)×H−1/2(∂ω∗) . Let us introduce the operators
T1 : H
1/2(∂ω∗) → H−1/2(∂ω∗) M1 : H
1/2(∂ω∗) → H1/2(∂ω∗)
h 7→ divτ (hn∇τud) h 7→ u
d − u
T2 : H
1/2(∂ω∗) → H1/2(∂ω∗) M2 : H
1/2(∂ω∗) → H−1/2(∂ω∗)
h 7→ hn∂nu
d h 7→ ∂n
u′+d − u
The Hessian can then be written under the form :
D2JKV (ω
∗)(h,h) = 2 [σ]
M1(h), T1(h)
T2(h),M2(h)
From the classical results of Maz’ya and Shaposhnikova on multipliers ([14], [22]), we get easily that
T1 and T2 are continuous operators. In fact, the compactness of the Hessian is a consequence of the
fact that both operators M1 and M2 are compact. We use a regularity argument : we remark that
M1 is the composition of the operators:
R1 : H
1/2(∂ω∗) → H
⋄ (∂Ω) and R2 : H
⋄ (∂Ω) → H
1/2(∂ω∗)
h 7→ −u′n φ 7→ ψ
where ψ is the trace on ∂ω∗ of Ψ solution of
−∆Ψ = 0 in Ω \ ω∗ and in ω∗, ,
[Ψ] = 0 on ∂ω∗,
[σ∂nΨ] = 0 on ∂ω
Ψ = φ on ∂Ω.
While R1 is a continuous operator, we prove that R2 is compact. Let us express u|∂ω∗ = ψ. We use
the integral formula of u to obtain:
I + µKω∗
σ2 + σ1
S∂Ω∂ω∗
µK∂ω∗∂Ω
σ2 + σ1
(u)|∂ω∗
(∂nu)|∂Ω
σ1 + σ2

K∂Ω∂ω∗φ

The matricial operator arising in this equation appeared also in (56). It has a continuous inverse
thanks to Theorem 7. Let us express u|∂ω∗ = ψ:
I + µKω∗)− µS∂Ω∂ω∗S
Ω K∂ω∗∂Ω
σ1 + σ2
K∂Ω∂ω∗ − S∂Ω∂ω∗S
I −KΩ)
φ. (69)
Since the operators K∂Ω∂ω∗ and S∂Ω∂ω∗ are compact, the operator R2 is compact, hence M1 is
compact. The proof of compactness of M2 is similar. Let us mention that a similar strategy of proof
can be found in [5].
The natural question is then to quantify how is this optimization problem ill-posed. This question
is directly in related to the rate at which the singular values of the Hessian operator are decreasing.
Equation (69) shows that this rate is the one of the operators K∂Ω∂ω∗ and S∂Ω∂ω∗ . Now, since for
every u ∈ H1/2(∂Ω), the functions K∂Ω∂ω∗u and S∂Ω∂ω∗u are harmonic outside of ∂Ω and therefore in
Ω, their restrictions on ∂ω∗ are as smooth as ∂ω∗. We conclude that if ∂ω∗ is C∞ then the restriction
belongs to each Hs(∂ω∗) for s > 1/2 then that if λn denotes the n
th eigenvalue of D2JKV (ω
∗), then
λn = o(n
−s) for all s > 0.
References
[1] L Afraites, and M. Dambrine and D. Kateb. Conformal mappings and shape derivatives for the
transmission problem with a single measurement. Preprint HAL 2006 to appear in Numerical
Functional Analysis and Optimization.
[2] L Afraites, and M. Dambrine, and K. Eppler and D. Kateb. Detecting perfectly insulated obstacles
by shape optimization techniques of order two. Preprint HAL 2006.
[3] K. Astala, and L. Pävärinta. Calderón’s inverse conductivity problem in the plane. Ann. of Math.,
(163), 265-299.
[4] M. Delfour, and J.-P. Zolesio. Shapes and Geometries: Analysis, Differential Calculus, and
Optimization SIAM, (2001).
[5] K. Eppler, and H. Harbrecht. A regularized Newton method in electrical impedance tomography
using Hessian information, Control and Cybernetics (34), 203-225.
[6] A. Henrot, and M. Pierre. Variation et optimization de formes. Springer Mathématiques et
Applications, volume 48 (2005).
[7] F. Hettlich, and W. Rundell. The determination of a discontinuity in a conductivity from a single
boundary measurement, Inverse Problems 14 (1998), 67-82.
[8] F. Hettlich, and W. Rundell. A Second Degree Method for Nonlinear Inverse Problems, SIAM J.
Numer. Anal., 37, No.2, (1999), 587–620.
[9] B. Hofmann. Approximation of the inverse electrical impedance tomography problem by an inverse
transmission problem, Inverse problems 14 (1998), 1171-1187.
[10] K. Ito, K. Kunisch, and Z. Li. Level-set function approach to an inverse interface problem,
Inverse Problems 17 (2001), 1225-1242.
[11] R. Kress. Linear Integral Equations. Springer-Verlag, Applied Mathematical Sciences (82).
[12] A. Kirsch. Surface Gradients and Continuity Properties for some integral operators in Classical
Scattering Theory. Mathematical Methods in the Applied Sciences, Vol 11 (1989), 789-804.
[13] A. Kirsch. The Domain Derivative and Two Applications in Inverse Scattering Theory, Inverse
Problems 9 (1993), 81-96.
[14] V.G. Maz’ya and T.O. Shaposhnikova. Theory of multipliers in spaces of differentiable functions,
Pitman, Boston, Monographs and Studies in Mathematics, 23, (1985).
[15] A.I. Nachmann. Reconstruction from boundary measurements, Ann. of Math., 128 (1988),
531-576
[16] A.I. Nachmann. Global uniqueness for a two dimensional inverse boundary value problem, Ann.
of Math., 143 (1996), 71-96
[17] R.G. Novikov. A multidimensional inverse spectral problem for the equation ∆ψ + (v(x) −
Eu(x))ψ = 0, Funktsional. Anal. i Prilozhen. 22 (1988), no 4, 11-22, 96; translation in Funct.
Anal. Appl. 22, 263-272.
[18] O. Pantz. Sensibilité de l’équation de la chaleur aux sauts de conductivité, C.R. Acad. Sci. Paris,
Ser. I 341-5 (2005), 333-337.
[19] J. Simon. Second variation for domain optimization problems, In Control and estimation of
distributed parameter systems, F. Kappel, K. Kunisch and W. Schappacher ed., International
Series of Numerical Mathematics, no 91, Birkhäuser, 361-378.
[20] J. Sokolowski and Jean-Paul Zolesio. Introduction to Shape Optimization: Shape Sensitivity
Analysis, Springer-Verlag (1992).
[21] J. Sylvester, and G. Uhlmann. A global uniqueness for an inverse boundary value problem Ann.
of Math. 125 (1987), 153-169.
[22] H. Triebel. Theory of Function Spaces, Birkhaüser (1983).
	Introduction and statement of the results.
	Preliminary results. 
	Elements of shape calculus
	Existence of the second order derivative of the state. Proof of Theorem ??.
	Preliminary results.
	Proof of existence of the second order derivative.
	Derivation of (??) from the weak formulation.
	How to recover (??) by formal differentiation of the boundary conditions.
	Justification of the formal computations.
	Case of Neumann boundary conditions.
	Second order derivatives for the criterion.
	Proof of Theorem ??.
	Analysis of stability. Proof of Theorem ??
