\begin{document}

\title{Distributed shape derivative via averaged adjoint method and applications}

\author[1]{Antoine Laurain\thanks{Universidade de S\~{a}o Paulo, Instituto de Matem\'atica e Estat\'istica, Departamento de Matem\'atica Aplicada, Rua do Mat\~{a}o, 1010
		Cidade Universitaria, CEP 05508-090, S\~{a}o Paulo, SP, Brazil}}
\author[2]{Kevin Sturm\thanks{Universit\"at Duisburg-Essen, Fakult\"at f\"ur Mathematik,
		Thea-Leymann-Stra\ss e 9,  D-45127 Essen, Germany;  {\bf email}: kevin.sturm@uni-due.de}}
\affil[1]{Department of Mathematics, University of S\~{a}o Paulo}
\affil[2]{Faculty of Mathematics, University of Essen-Duisburg}

\maketitle

\begin{abstract} 
	The structure theorem of Hadamard-Zol\'esio states that the derivative of a shape functional is a distribution on the boundary of the domain depending only on the normal perturbations of a smooth enough boundary. 
	Actually the domain representation, also known as distributed shape derivative, is more general than the boundary expression as it is well-defined for shapes having a lower regularity.
	It is customary in the shape optimization literature  to assume regularity of the domains and use the boundary expression of the shape derivative for numerical algorithms. 
	In this paper we describe several advantages of the distributed shape derivative in terms of generality, easiness of computation and numerical implementation. 
	We identify a tensor  representation of the distributed shape derivative, study its properties and show how it allows to recover the boundary expression directly. 
	We  use a novel Lagrangian approach, which is applicable to a large class of shape optimization problems, to compute the distributed shape derivative. 
	We also apply the technique to retrieve the distributed shape derivative for  electrical impedance tomography.  
	Finally we explain how to adapt the level set method to the distributed shape derivative framework and present numerical results. 
\varepsilonnd{abstract}

\subjclass{49Q10, 35Q93, 35R30, 35R05}\newline

\text{ker}ywords{Shape optimization, distributed shape gradient, electrical impedance tomography, Lagrangian method, level set method}

\maketitle

\section*{Introduction}

In his research on elastic plates \cite{Ha} in 1907, Hadamard 
showed how to obtain the derivative of a shape functional $J(\mathcal{B}ega)$ 
by considering normal perturbations of the  boundary $\partial\mathcal{B}ega$ of a smooth set  $\mathcal{B}ega$. This fundamental result of shape optimization  was  made rigorous later by Zol\'esio \cite{MR2731611} in the so-called ``structure theorem''. When $J(\mathcal{B}ega)$  and the domain are smooth enough, one may also write the shape derivative as an integral over $\partial\mathcal{B}ega$, which is the canonical form in the shape optimization literature. 

However, when $\mathcal{B}ega$ is less regular, the shape derivative can  often be written as a domain integral even when the boundary expression is not available. The {\it domain expression} also known as {\it distributed shape derivative} has been generally ignored in the shape optimization literature for several reasons:  firstly the boundary representation provides a straightforward way of determining an explicit descent direction since it depends linearly on the boundary perturbation $\theta$ and not on its gradient, secondly this descent direction only needs to be defined on the boundary. When considering the domain expression, these two advantages disappear as the shape derivative is defined on $\mathcal{B}ega$ and depends on the gradient of $\theta$, so that a partial differential equation needs to be solved to obtain a descent direction $\theta$ on  $\mathcal{B}ega$. 

It seems that these drawbacks  would definitely rule out the distributed shape derivative, however they turn out to be less dramatic than expected in many situations and the domain formulation has other less foreseeable advantages  over the boundary representation. In this paper we advocate for the use of the distributed shape derivative and discuss the advantages of this formulation. 

The boundary representation has the following drawbacks. First of all  if the data is not smooth enough the integral representation does not exist so that the more general domain representation is the only rigorous alternative. Even when the boundary representation exists and has the form
$\int_{\partial\mathcal{B}ega} g\, \theta\cdot n $, it is usually not legitimate to choose $\theta\cdot n = -g$ on $\partial\mathcal{B}ega$ for a descent direction if $g$ is not smooth enough, for instance if $g\in L^1(\partial\mathcal{B}ega)$. Therefore, a smoother $\theta$ must be chosen, which requires to solve a partial differential equation on the boundary $\partial\mathcal{B}ega$. When taking $\theta\cdot n = -g$ is legitimate, it might still not be desirable as this may yield a $\theta$ with  low regularity, in which case one needs to regularize $\theta$ on the boundary as well. 

In these cases the first advantage of the boundary representation disappears. The second advantage of the boundary 
representation is that the perturbation field only needs to be defined on the boundary instead of on the whole domain, reducing the cost of the computation. Actually, the distributed shape derivative also has its support on the boundary, and may be computed in a small neighborhood of the boundary so that the additional cost is minimal. 
In addition, in most shape optimization applications, $g$ is the restriction of a function defined in a neighborhood of the boundary and not a quantity depending only on the boundary such as the curvature. Therefore from a practical point of view, $g$ must be evaluated in a neighborhood of $\partial\mathcal{B}ega$ anyway. 
Also, in many numerical applications, $\theta$ must be extended to a neighborhood of $\Gamma$ or even to the entire domain $\mathcal{B}ega$. This is the case for level set methods for instance, where the level set function must be updated on $\mathcal{B}ega$, or when one wishes to update the mesh along with the domain update, to avoid re-meshing the new domain. The distributed shape derivative then directly gives an extension of $\theta$ well-suited to the optimization problem. 

Recent results have shown that the distributed shape derivative is also more accurate than the boundary representation from a numerical point of view; see \cite{MR3348199} for a comparison. Indeed functions such as gradients of the state and adjoint state appearing in the distributed shape derivative only need to be defined at grid points and not on the interface. Therefore one avoids interpolation of these irregular terms.   This is particularly useful for transmission problems where the boundary representation requires to compute the jump of a function over the interface, a delicate and error-prone operation from the numerical point of view. 

Having considered these equivalent expressions of the shape derivative (i.e. boundary and domain expression) leads to a general form of the shape derivative using tensors. We introduce such a tensor representation in Section  \ref{sec:tensor_representation} which covers a large class of problems and in particular contains the boundary and domain expression. We show how this abstract form allows to identify simple relations between the domain and boundary expressions of the shape derivative.  

In this paper we also extend and simplify the averaged adjoint method from \cite{sturm13}, a  Lagrangian-type method which is well-suited to compute the shape derivative of a cost function in an efficient way.

Lagrangian methods are commonly used in shape optimization and have the advantage of providing the shape derivative without the need to compute the material derivative of the state; see \cite{MR862783,MR948649,MR2166150,sturm13,MeftBelh13,SturmHeomHint13}.    
Compared to these known shape-Lagrangian methods, the averaged adjoint method is fairly general due to minimal required conditions. The assumptions are for instance less restrictive than those required for the theorem of Correa-Seeger \cite{MR948649}, therefore it can be applied to more general situations such as non-convex functionals.  As the direct approach our method can also be applied for problems depending on nonlinear partial differential equations. In this paper we give an example of application to a transmission problem (in electrical impedance tomography - see Section \ref{section4}).  Our method  provides the domain expression of the shape derivative and the boundary expression can be computed easily from the tensor representation of the domain expression.

To complete the numerical implementation aspect, we also show  how the domain expression of the shape derivative can be used in the level set method framework \cite{MR2033390,FLS,MR965860,MR2806573,HLsegmentation,MR2472382,MR2459656,MR2356899}. 
The level set method can be modified to use the domain expression which leads to a method which is actually easier to implement.
Combining all these techniques, we obtain a straightforward and general  way of solving the shape optimization problem, from the rigorous theoretical computation of the shape derivative to the numerical implementation.

In Section \ref{section1} we recall the concept of shape derivative and the structure theorem on an abstract level. 
In Section \ref{section2} a shape-Lagrangian method, the averaged adjoint method, is described. 
In Section \ref{sec:tensor_representation} we identify a general tensor representation of the shape derivative, establish some of its properties, and give a few examples.
In Section \ref{section3} we explain how to compute descent directions for the distributed shape derivative for use in gradient methods. 
In Section \ref{section4} we apply the results of Sections \ref{section2} and  \ref{sec:tensor_representation} to the inverse problem of electrical impedance tomography. 
In Section \ref{section5} we extend the level set method to the case of the distributed shape derivative and finally  in Section \ref{section6} we show numerical results for various problems including the  problem of electrical impedance tomography.

\section{The structure theorem revisited}\label{section1}

Our aim in this section is to describe properties of the shape derivative on an abstract level and to emphasize that all representations of the shape derivative satisfy the same {\it structure theorem}.

Let $\mathcal{P}(D)$ be the set of subsets of $D\subset\mathbf{R}^d$ compactly contained in $D$, where the so-called ``universe''
$D\subset \mathbf{R}^d$ is assumed to be open and bounded.
Define for $k\geq 0$ and $0\leq \alpha \leq 1$,
\begin{align}
C^{k,\alpha}_c(D,\mathbf{R}^d) &:=\{\theta\in C^{k,\alpha}(D,\mathbf{R}^d)|\quad \theta\text{ has compact support in } D\}.
\varepsilonnd{align}
Also for given domain $\mathcal{B}ega\subset D$ with at least a $C^1$ boundary we introduce the space of vector field
\begin{align}
C^{k,\alpha}_{\partial \mathcal{B}ega}(D ,\mathbf{R}^d) &:=\{\theta\in C^{k,\alpha}_c( D,\mathbf{R}^d)|\quad \theta\cdot n =0 \text{ on } \partial \mathcal{B}ega\}
\varepsilonnd{align}
where $n$ is the outward unit normal vector to $\mathcal{B}ega$.

Consider a vector field $\theta\in C^{0,1}_c(D,\mathbf{R}^d)$ and the associated flow
$\Phi_t^{\theta}:\overline{D}\rightarrow \mathbf{R}^d$, $t\in [0,\tau]$ defined for each $x_0\in \overline{D}$ as $\Phi_t^{\theta}(x_0):=x(t)$, where $x:[0,\tau]\rightarrow \mathbf{R}$ solves 
\begin{align}\label{Vxt}
\begin{split}
\dot{x}(t)&= \theta(x(t))    \quad \text{ for } t\in (0,\tau),\quad  x(0) =x_0.
\varepsilonnd{split}
\varepsilonnd{align}
We will sometimes use the simpler notation $\Phi_t=\Phi_t^{\theta}$ when no confusion is possible.
Since $\theta\in C^{0,1}_c(D,\mathbf{R}^d)$ we have by Nagumo's theorem \cite{MR0015180} that for fixed $t\in [0,\tau]$ the flow $\Phi_t$ is a homeomorphism from $D$ into itself and maps boundary onto boundary and
interior onto interior. Further, we consider the family 
\begin{equation}\label{domain}
\mathcal{B}ega_t := \Phi_t^{\theta}(\mathcal{B}ega) 
\varepsilonnd{equation}
of perturbed domains. 

In the following let $J : \mathfrak P \rightarrow \mathbf{R}$ be a shape function defined on
some admissible set $\mathfrak P\subset \mathcal P(D)$.

\begin{definition}\label{def1} 
	The Eulerian semiderivative of $J$ at $\mathcal{B}ega$ in direction $\theta \in C^{0,1}_c(D,\mathbf{R}^d)$, when the limit exists,
	is defined by
	\begin{equation}
	dJ(\mathcal{B}ega)(\theta):= \lim_{t \searrow 0}\textnormal{fr}ac{J(\mathcal{B}ega_t)-J(\mathcal{B}ega)}{t}.
	\varepsilonnd{equation}
	\begin{itemize}
		\item[(i)] $J$ is said to be \textit{shape differentiable} at $\mathcal{B}ega$ if it has a Eulerian semiderivative at $\mathcal{B}ega$ for all $\theta \in C^\infty_c(D,\mathbf{R}^d)$ and the mapping
		\begin{align*}
		dJ(\mathcal{B}ega):  C^\infty_c(D,\mathbf{R}^d)  &  \to \mathbf{R},\; \theta     \mapsto dJ(\mathcal{B}ega)(\theta)
		\varepsilonnd{align*}
		is linear and continuous, in which case $dJ(\mathcal{B}ega)(\theta)$ is called the \textit{shape derivative} at $\mathcal{B}ega$.
		\item[(ii)] 
		The shape derivative $dJ(\mathcal{B}ega)$ is of finite order if there is an integer $l\ge 0$ and a constant $c>0$ such that for each compact $K\subset D$ 
		$$|dJ(\mathcal{B}ega)(\theta)|  \le c \|\theta\|_{l} \quad \forall \theta\in C^\infty_c(K,\mathbf{R}^d),$$
		where $\|\theta\|_{l} := \sum_{|\alpha| \le l} |D^\alpha \theta|_\infty$.
		The smallest such integer $l\ge 0$ is called order of $dJ(\mathcal{B}ega)$.
	\varepsilonnd{itemize}
\varepsilonnd{definition}

The shape derivative from Definition \ref{def1} has a particular structure. 
Intuitively, it is clear that the form functional stays constant 
for a transformation $\Phi$ that leaves $\mathcal{B}ega$ unchanged, that is $\Phi(\mathcal{B}ega)=\mathcal{B}ega$,  even if some points inside $\mathcal{B}ega$ move  and consequently the shape derivative is zero in this case. This property is valid  when $\mathcal{B}ega$ is open or closed; cf. \cite{MR2731611}. Mathematically, this is expressed in the following basic theorem 
proved in \cite{zolesio1979identification}.

\begin{theorem}\label{lem:1}
	Let  $\mathcal{B}ega \in \mathfrak P$ be open or closed. Let $\theta \in C^{0,1}_c(D,\mathbf{R}^d)$ be a vector field with 
	compact support in $\mathcal{B}ega$ and denote by
	$\Phi_t$ its flow defined in \varepsilonqref{Vxt}.  

	Then we have 
	$$  dJ(\mathcal{B}ega)(\theta)=0. $$
\varepsilonnd{theorem}

Note that the shape derivative of $J(\mathcal{B}ega)$ always exists for vector fields with compact support in $\mathcal{B}ega$, even if it does not exist for other vector fields. An important consequence of Theorem \ref{lem:1}, also for numerical methods, is that {\it independently of the representation} of the shape derivative and the regularity of the domain $\mathcal{B}ega$, the values of $\theta$ outside the boundary of $\mathcal{B}ega$ have no influence on the shape derivative. 

\begin{corollary}
	Let $\mathcal{B}ega \in \mathfrak P$ be a set with $C^{1}$-boundary. 
	Assume that $J$ is shape differentiable on $\mathfrak P$.
	Let $\theta \in C^{0,1}_{\partial\mathcal{B}ega}(D,\mathbf{R}^d)$.
	Then we have
	$$ dJ(\mathcal{B}ega)(\theta)=0. $$
\varepsilonnd{corollary}

The previous discussion immediately yields the following fundamental result of shape optimization.
\begin{theorem}[Structure Theorem]\label{thm:structure_theorem}
	Assume $\Gamma :=\partial \mathcal{B}ega $ is compact and $J$ is shape differentiable. Denote the \textit{shape derivative} by
	\begin{equation}
	dJ(\mathcal{B}ega ):C^\infty_c(D,\mathbf{R}^d)\rightarrow \mathbf{R},\quad \theta \mapsto dJ(\mathcal{B}ega)(\theta).
	\varepsilonnd{equation} 
	Assuming $dJ(\mathcal{B}ega )$ is of order $k\ge 0$ and $\Gamma$ of class $C^{k+1}$, then there exists a linear and continuous functional $g: C^k(\Gamma)\rightarrow \mathbf{R}$ such that
	\begin{equation}\label{volume}
	dJ(\mathcal{B}ega)(\theta)=g(\theta_{|\Gamma}\cdot n),
	\varepsilonnd{equation}

\varepsilonnd{theorem}

\begin{proof}
See \cite[pp. 480-481]{MR2731611}.
\varepsilonnd{proof}

\section{Shape derivatives via averaged adjoint method}\label{section2}

Lagrangian methods in shape optimization allow to compute the shape derivative of functions depending on the solution of partial differential equations without the need to compute the material derivative of the partial differential equations; see \cite{MR2731611} for a description of such a method in the linear case. 
Here we extend  and simplify the averaged adjoint method, a Lagrangian-type method introduced in \cite{sturm13}. 
With this approach the computation of the domain representation of the shape derivative is fast, the retrieval of the boundary form is convenient, 
and no saddle point assumptions is required unlike in \cite{MR2731611}.

Let two vector spaces $E = E(\mathcal{B}ega), F=F(\mathcal{B}ega)$ and  $\tau>0$ be given, and consider a parameterization 
	$\mathcal{B}ega_t = \Phi_t(\mathcal{B}ega)$ for $t\in [0,\tau]$. 
	Ultimately, our goal is to differentiate shape functions of the type 
	$J(\mathcal{B}ega_t)$ which can be written using a Lagrangian as $J(\mathcal{B}ega_t) = \mathcal{L}(\mathcal{B}ega_t, u^t,\hat\psi)$,
	where $u^t\in E(\mathcal{B}ega_t)$ and $\hat\psi\in F(\mathcal{B}ega_t) $.
	The main appeal of the Lagrangian is that 
	we actually only need to compute the derivative with respect to $t$ of $\mathcal{L}(\mathcal{B}ega_t,\hat\varphi,\hat\psi)$ to compute the derivative of $J(\mathcal{B}ega_t)$,  indeed this is the main result of Theorem \ref{thm:sturm}, but this requires a few explanations.

	Since $\mathcal{L}(\mathcal{B}ega_t,\hat\varphi,\hat\psi)$,  is often constituted of integrals on $\Phi_t(\mathcal{B}ega)$, using a change of variable we can rewrite these integrals to integrals on the fixed domain $\mathcal{B}ega$, and consequently transfer the dependence on $t$ to the integrand. 
	However, in the process appear the composed functions $\hat\varphi\circ\Phi_t\in E(\mathcal{B}ega)$ and $\hat\psi\circ\Phi_t\in F(\mathcal{B}ega)$, whose derivatives are not straightforward to compute since $\hat\varphi$ and $\hat\psi$ are defined on the moving spaces $E(\mathcal{B}ega_t)$ and $F(\mathcal{B}ega_t)$.

	Fortunately, and this is the crucial point of the shape-Lagrangian approach,  to compute the shape derivative we can reparameterize the problem by considering
	$\mathcal{L}(\mathcal{B}ega_t, \Psi_t\circ \varphi, \Psi_t\circ \psi)$ instead of $\mathcal{L}(\mathcal{B}ega_t,\hat\varphi,\hat\psi)$, where 
	$\Psi_t$ is an appropriate bijection between
	$E(\mathcal{B}ega)$ and $E(\mathcal{B}ega_t)$, and $\varphi\in E(\mathcal{B}ega)$, $\psi\in F(\mathcal{B}ega)$.
	Now the change of variable in the integrals yields functions $\varphi$ and $\psi$ in the integrand, which are defined on fixed spaces.
	In this paper $E$ and $F$ are $H^1$-spaces,  and in this case we may consider the particular reparameterization
	$\mathcal{L}(\mathcal{B}ega_t,\varphi\circ\Phi_t^{-1},\psi\circ\Phi_t^{-1})$. 
	For spaces such as 
	$H(\text{curl};\mathcal{B}ega)$, other transformations $\Psi_t$ can be used;  see \cite{HLY,MR2971616,MR2002150}. 

Thus we are led to consider general functions of the type $G:[0,\tau]\times E\times F \rightarrow \mathbf{R}$ with
	$$ G(t,\varphi,\psi):
	=\mathcal{L}(\Phi_t(\mathcal{B}ega),\varphi\circ\Phi_t^{-1},\psi\circ\Phi_t^{-1}).$$
	This is precisely what we do in \varepsilonqref{eq:G_on_the_moved_domain} when showcasing an application of the method.
	The main result of this section, Theorem \ref{thm:sturm}, shows that to obtain the shape derivative of $\mathcal{L}$, it is enough to compute the derivative with respect to $t$ of $G$ while assigning the proper values to $\varphi$ and $\psi$. The main ingredient is the introduction of the averaged adjoint equation. 

	In addition, in this paper we consider the following specific form
	\begin{equation}
	\label{G_lag}
	G(t,\varphi,\psi):= a(t,\varphi,\psi) + b(t,\varphi),
	\varepsilonnd{equation}
	where
	$$ a:[0,\tau]\times E \times  F \rightarrow \mathbf{R}, \qquad b:[0,\tau]\times E \rightarrow \mathbf{R} ,$$ 
	are functions such that $\psi \mapsto a(t,\varphi,\psi)$ is linear for all 
	$t\in [0,\tau]$ and $\varphi\in E$. The function $G$ is commonly called {\it Lagrangian}, hence the name of the method.  
	In the applications we have in mind, the function 
	$b$ arises from the objective function while $a$ corresponds to the constraint, after transporting back  to the fixed domain $\mathcal{B}ega$.

Throughout the paper, the Greek letters $\varphi$ and $\psi$ are used for variables, while the roman letters $u,p$ are used for the solutions of the state and adjoint states, respectively.

Let us  assume that for each $t\in [0,\tau]$ the equation
	\begin{equation}\label{eq:state_G}
	d_\psi G(t,u^t,0;\hat\psi)= a(t,u^t, \hat\psi) = 0\;\text{ for all } \hat\psi \in  F.
	\varepsilonnd{equation}
	admits a unique solution $u^t\in E$. 
Further we make the following assumptions for  $G$.

\begin{thmx}
	\label{amp:gateaux_diffbar_G}
	\label{amp:affine-linear}
	For every $(t,\psi)\in [0,\tau]\times F$
	\begin{enumerate}
		\item[(i)]  $[0,1]\ni s\mapsto G(t,su^t+s(u^t-u^0),\psi)$ is absolutely continuous. 
		\item[(ii)] $[0,1]\ni s\mapsto d_\varphi G(t,su^t+(1-s)u^0,\psi;\hat{\varphi})$ belongs to $L^1(0,1)$ for all $\hat{\varphi}\in E$.
	\varepsilonnd{enumerate}
\varepsilonnd{thmx}

When Assumption \ref{amp:affine-linear} is satisfied,  for $t\in [0,\tau]$ we introduce the \textit{averaged adjoint equation} associated with $u^t$ and $u^0$: 
Find $p^t\in F$ such that 

\begin{equation}\label{averated_}
\int_0^1 d_\varphi G(t,su^t+(1-s)u^0,p^t;\hat{\varphi})\, ds =0 \quad \text{ for all } \hat{\varphi}\in E.
\varepsilonnd{equation}
Notice that, in view Assumption \ref{amp:affine-linear}, for all $t\in[0,\tau]$,
\begin{equation}\label{eq:main_averaged}
G(t,u^t,p^t)-G(t,u^0,p^t) = \int_0^1 d_\varphi G(t,su^t+(1-s)u^0,p^t;u^t-u^0)\, ds =0.
\varepsilonnd{equation}

We can now state the main result of this section.

	\begin{thmy}
		We assume that
		$$ \lim_{t\searrow 0} \textnormal{fr}ac{G(t,u^0,p^t)-G(0,u^0,p^t)}{t}=\partial_tG(0,u^0,p^0).$$
	\varepsilonnd{thmy}

	\begin{theorem}
		\label{thm:sturm}
		Let  (H0) and (H1) be satisfied and assume there exists a unique solution $p^t$  of the averaged adjoint equation \varepsilonqref{averated_}.

		Then for $\psi \in F$  we obtain
		\begin{equation}\label{eq:dt_G_single}
		\dt b(t,u^t) |_{t=0} = \dt(G(t,u^t,\psi))|_{t=0}=\partial_t G(0,u^0,p^0).
		\varepsilonnd{equation}
	\varepsilonnd{theorem}

\begin{proof}
	Put $g(t) := G(t,u^t,0)-G(0,u^0,0)$, and note that $g(t) = G(t,u^t,\psi)-G(0,u^0,\psi)$ for all $\psi \in F$ and $g(0)=0$. We have to show that 
	$$g'(0):= \lim_{t\searrow 0}\textnormal{fr}ac{G(t,u^t,0)-G(0,u^0,0)}{t} \quad \text{ exists. }$$  
	Thanks to Assumption \ref{amp:affine-linear} we can define the averaged adjoint $p^t$ and using that $G$ is affine with respect 
	to the third argument, we obtain 
	$$g(t) =  \underbrace{G(t,u^t,p^t)-G(t,u^0,p^t)}_{=0 \mbox{ in view of } \varepsilonqref{eq:main_averaged}} + G(t,u^0,p^t)-G(0,u^0,p^t) .$$
	Dividing by $t>0$ and using Assumption (H1) yields
	$$ g'(0) = \lim_{t\searrow 0} \textnormal{fr}ac{g(t) - g(0)}{t} =  \lim_{t\searrow 0}\textnormal{fr}ac{G(t,u^0,p^t)-G(0,u^0,p^t)}{t}  = \partial_t G(0,u^0,p^0) $$
	which concludes the proof.
\varepsilonnd{proof}

\begin{remark}
	In terms of $a$ and $b$, equation \varepsilonqref{averated_} reads:
	$$ \int_0^1 d_\varphi a(t,su^t+(1-s)u^0,p^t;\hat{\varphi})\, ds = - \int_0^1 d_\varphi b(t,su^t+(1-s)u^0;\hat{\varphi})\, ds $$
for all $\hat{\varphi}\in E$.
	If $\varphi \mapsto a(t,\varphi,\psi)$ is in addition linear, then  \varepsilonqref{averated_}  becomes 
	$$        a(t,\hat{\varphi},p^t) = - \int_0^1 d_\varphi b(t,su^t+(1-s)u^0;\hat{\varphi})\, ds $$
 for all $\hat{\varphi}\in E$.
\varepsilonnd{remark}

\section{Tensor representation of the shape derivative}\label{sec:tensor_representation}

In this section we identify tensor representations of the shape derivative that correspond to a large class of problems studied in the literature for PDE-constrained shape optimization. This tensor representation has several interesting properties that we investigate. In particular we exhibit the link between this tensor representation and the usual boundary expression of the shape derivative.

\subsection{Definition and properties}\label{def_prop}

\begin{definition}\label{def:tensor}
Let $\mathcal{B}ega \in \mathfrak P$ be a set with $C^k$-boundary, $k\ge 1$. 
	A shape differentiable function $J$ of order $k$ is said to admit a tensor representation if
	there exist tensors $\mathbf{S}_l\in  L^1(D, \mathcal L^l(\mathbf{R}^d,\mathbf{R}^d))$  and $\mathfrak{S}_l\in L^1(\partial \mathcal{B}ega; \mathcal L^l(\mathbf{R}^d,\mathbf{R}^d))$, $l=0,..,k$,
	such that
	\begin{equation}\label{ea:volume_from}
	dJ(\mathcal{B}ega)(\theta) = \sum_{l=0}^k \int_D \mathbf{S}_l\cdot D^l\theta \,dx + \int_{\partial \mathcal{B}ega} \mathfrak{S}_l\cdot D^l_\Gamma \theta\, ds \quad \text{ for all } \theta\in C^k_c(D,\mathbf{R}^d),
	\varepsilonnd{equation}
	where $D_\Gamma \theta: = D\theta -(D\theta n)\otimes n$ is the tangential derivative of $\theta$ along $\partial\mathcal{B}ega$. Here 
		$\mathcal L^l(\mathbf{R}^d, \mathbf{R}^d) $ denotes the space of multilinear maps from $\mathbf{R}^d \times \cdots  \times \mathbf{R}^d$ to $\mathbf{R}^d$.
\varepsilonnd{definition}

Most if not all examples involving PDEs from classical textbooks \cite{MR2731611,MR2512810,SokZol92} can be written in the form \varepsilonqref{ea:volume_from}.

\begin{remark}
	\begin{itemize}
		\item[(a)] A particular case of the tensor representation \varepsilonqref{ea:volume_from}  is the Eshelby energy momentum tensor in continuum mechanics introduced in \cite{MR0489190}; see also \cite{MR3013681}. In this case only $\mathbf{S}_1$ is not zero.
		\item[(b)] When $J$ is is shape differentiable in $\mathcal{B}ega$ then by definition  $\theta \mapsto dJ(\mathcal{B}ega)(\theta)$ is a distribution, and if  $\partial \mathcal{B}ega$ is compact, the distribution $\theta \mapsto dJ(\mathcal{B}ega)(\theta)$  is of finite order.
		\item[(c)] If $dJ(\mathcal{B}ega)$ is of order $k=1$ and $ |dJ(\mathcal{B}ega)(\theta)|  \le C\|\theta\|_{H^1(D,\mathbf{R}^d)} $ for all $\theta\in C^\infty_c(D,\mathbf{R}^d)$
		then by density of $C^\infty_c(D,\mathbf{R}^d)$ in $H^1_0(D,\mathbf{R}^d)$ the derivative $dJ(\mathcal{B}ega)$ 
		extends to a continuous functional on $H^1_0(D,\mathbf{R}^d)$, that is,
		$$ |\widehat{dJ(\mathcal{B}ega)}(\theta)|  \le c\|\theta\|_{H^1(D,\mathbf{R}^d)}  \quad \text{ for all } \theta\in H^1_0(D,\mathbf{R}^d).  $$ 
		Therefore by the  theorem of Riesz, we obtain a vector field $W$ in $H^1_0(D,\mathbf{R}^d)$ 
		such that
		$$ \forall \theta\in H^1_0(D,\mathbf{R}^d),\quad  \widehat{dJ(\mathcal{B}ega)}(\theta) = \int_D DW \cdot D \theta  + W \cdot \theta \, dx  $$
		and this defines a tensor representation with $\mathbf{S}_1= DW$, $\mathbf{S}_0=W$, $\mathfrak{S}_1=0$ and $\mathfrak{S}_0=0$.
		\item[(d)] The assumption that $\mathcal{B}ega$ be a set of class $C^k$ can be reduced if $\mathfrak{S}_l\varepsilonquiv 0$ for all $0\leq k_0\leq l\leq k$.
	\varepsilonnd{itemize}
\varepsilonnd{remark}

The tensor representation \varepsilonqref{ea:volume_from} is not unique in the sense that there might be
several ways to choose the tensors $\mathbf{S}_l$ and $\mathfrak{S}_l$. This is expressed by the fact that 
these tensors are correlated. We describe these relations below in the case $k=1$ in Proposition 
\ref{tensor_relations}, which also describes the link between the tensor representation and the usual boundary representation \varepsilonqref{volume} of the shape derivative. 

\begin{proposition}\label{tensor_relations}
	Let $\mathcal{B}ega $ be a subset of $D$ with $C^1$-boundary. Suppose
	that the derivative  $dJ(\mathcal{B}ega)$ has the representation
	\begin{equation}\label{eq:first_order_tensor}
	dJ(\mathcal{B}ega)(\theta) = \int_D \mathbf{S}_1\cdot D\theta +  \mathbf{S}_0\cdot \theta\,dx + \int_{\partial \mathcal{B}ega} \mathfrak{S}_1\cdot D_\Gamma \theta + \mathfrak{S}_0\cdot \theta\, ds .
	\varepsilonnd{equation}
	If $\mathbf{S}_l$ is of class $W^{1,1}$ in $\mathcal{B}ega$ and $D\setminus \overline \mathcal{B}ega$ then indicating by $+$ and $-$ the restrictions of the tensors to $\mathcal{B}ega$ and $D\setminus\overline \mathcal{B}ega$, respectively, we get
	\begin{equation} \label{eq:equvilibrium_strong}
	\begin{split}
		-\operatorname{div}(\mathbf{S}_1^+) + \mathbf{S}_0^+ &= 0 \quad \text{ in } \mathcal{B}ega \\
		-\operatorname{div}(\mathbf{S}_1^-) + \mathbf{S}_0^- &= 0 \quad \text{ in } D\setminus\overline \mathcal{B}ega.
	\varepsilonnd{split}
	\varepsilonnd{equation}
	Moreover, 

	we can rewrite the tensor representation  as a distribution on the boundary:
	$$dJ(\mathcal{B}ega)(\theta) = \int_{\partial \mathcal{B}ega} [(\mathbf{S}_1^+-\mathbf{S}_1^-)n] \cdot\theta + \mathfrak{S}_1\cdot D_\Gamma \theta + \mathfrak{S}_0\cdot \theta \, ds $$
	where $n$ denotes the outward unit normal vector to $\mathcal{B}ega$.
	
	If the boundary $\partial \mathcal{B}ega$ is $C^2$ and $\mathfrak{S}_1\in W^{1,1}(\partial \mathcal{B}ega; \mathcal L^1(\mathbf{R}^d,\mathbf{R}^d))$, then we obtain a more regular distribution, the so-called boundary expression of the shape derivative: 
	\begin{equation}\label{eq:general_boundary_exp}
	dJ(\mathcal{B}ega)(\theta) = \int_{\partial \mathcal{B}ega} g_1\, \theta\cdot n\, ds,
	\varepsilonnd{equation}
	where
	\begin{equation}
	\label{g1} 
	g_1 :=  [(\mathbf{S}_1^+-\mathbf{S}_1^-)n]\cdot n + \mathfrak{S}_0\cdot n + \mathfrak{S}_1\cdot D_\Gamma n - \operatorname{div}_\Gamma (\mathfrak{S}_1^T n) + \mathcal{H}(\mathfrak{S}_1^T n\cdot n). \varepsilonnd{equation}
	and $\mathcal{H}=\operatorname{div}_\Gamma n$ denotes the mean curvature\footnote{We define the mean curvature as the sum of the principal curvatures $\kappa_i$, that is, $\mathcal H := \sum_{i=1}^d \kappa_i$.} of $\partial\mathcal{B}ega$ while $\operatorname{div}_\Gamma := \operatorname{tr}(D_\Gamma)$ is the tangential divergence.  
\varepsilonnd{proposition}

\begin{proof}
Applying Theorem \ref{lem:1} we have
$$dJ(\mathcal{B}ega)(\theta) = \int_D \mathbf{S}_1\cdot D\theta +  \mathbf{S}_0\cdot \theta\,dx + \int_{\partial \mathcal{B}ega} \mathfrak{S}_1\cdot D_\Gamma \theta + \mathfrak{S}_0\cdot \theta\, ds =0 \quad \text{ for all }\theta\in C^1_c(\mathcal{B}ega\cup (D\setminus \overline \mathcal{B}ega), \mathbf{R}^d).  $$
An integration by parts shows \varepsilonqref{eq:equvilibrium_strong}.

Then, when $\partial\mathcal{B}ega$ is $C^1$, replacing \varepsilonqref{eq:equvilibrium_strong} in the expression of the shape derivative and using Green's formula we obtain  
\begin{equation}\label{eq:domain}
\begin{split}
	dJ(\mathcal{B}ega)(\theta) =& \int_{\partial \mathcal{B}ega} \mathfrak{S}_1\cdot D_\Gamma \theta + \mathfrak{S}_0\cdot \theta\,\,ds  + \int_{\partial \mathcal{B}ega} [(\mathbf{S}_1^+-\mathbf{S}_1^-)n]\cdot \theta \, ds\\
	&+ \int_\mathcal{B}ega( \underbrace{-\operatorname{div}(\mathbf{S}_1^+) + \mathbf{S}_0^+}_{=0})\cdot \theta \, dx + \int_{D\setminus \overline \mathcal{B}ega} (\underbrace{-\operatorname{div}(\mathbf{S}_1^-) + \mathbf{S}_0^-}_{=0})\cdot \theta \, dx\\
	\stackrel{\varepsilonqref{eq:equvilibrium_strong}}{=} & \int_{\partial \mathcal{B}ega} [(\mathbf{S}_1^+-\mathbf{S}_1^-)n]\cdot \theta +\mathfrak{S}_1\cdot D_\Gamma \theta + \mathfrak{S}_0\cdot \theta\, \, ds\quad \text{ for all } \theta\in C^1_c(D, \mathbf{R}^d).
\varepsilonnd{split}
\varepsilonnd{equation}
With a slight abuse of notation we keep the same notation $n$ for the extension of the normal to a neighborhood of $\partial \mathcal{B}ega$. Let $\theta\in C^1(\overline D,\mathbf{R}^d)$ and define
$\theta_\tau  :=\theta - (\theta\cdot n) n$ the tangential part of $\theta$. Then $\theta_\tau \cdot n =0$  on 
$\partial \mathcal{B}ega$ and hence  from the structure theorem we get $dJ(\mathcal{B}ega)(\theta_\tau)=0$ which yields in view of \varepsilonqref{eq:domain}:
\begin{equation}\label{eq:boundary_exp}
\begin{split}
	dJ(\mathcal{B}ega)(\theta) &=  dJ(\mathcal{B}ega)((\theta\cdot n)n) \\
	&= \int_{\partial \mathcal{B}ega} ((\mathbf{S}_1^+-\mathbf{S}_1^-)n\cdot n) (\theta\cdot n) 
	+\mathfrak{S}_1\cdot D_\Gamma (n(\theta\cdot n)) + (\mathfrak{S}_0\cdot n) (\theta\cdot n) \, ds\\
	& =  \int_{\partial \mathcal{B}ega} ((\mathbf{S}_1^+-\mathbf{S}_1^-)n\cdot n) (\theta\cdot n) + (\mathfrak{S}_0\cdot n) (\theta\cdot n)\, ds\\
	& + \int_{\partial \mathcal{B}ega} \mathfrak{S}_1\cdot D_\Gamma n (\theta\cdot n) + n\cdot \mathfrak{S}_1\nabla_\Gamma (\theta\cdot n) \, ds,
\varepsilonnd{split}
\varepsilonnd{equation}
where we used that for all functions $f\in C^1(\mathbf{R}^d,\mathbf{R}^d)$ and $g\in C^1(\mathbf{R}^d)$ we have
\begin{equation}
\begin{split}
	D(gf) & = g Df + f\otimes \nabla g. 

\varepsilonnd{split}
\varepsilonnd{equation}

Finally using $\mathfrak{S}_1\in W^{1,1}(\partial \mathcal{B}ega; \mathcal L^1(\mathbf{R}^d,\mathbf{R}^d))$ we integrate by parts on the boundary $\partial \mathcal{B}ega$ to transform the last term in \varepsilonqref{eq:boundary_exp}
$$
\int_{\partial \mathcal{B}ega} n\cdot \mathfrak{S}_1\nabla_\Gamma (\theta\cdot n) \, ds = \int_{\partial \mathcal{B}ega} (- \operatorname{div}_\Gamma (\mathfrak{S}_1^T n) + \mathcal{H}(\mathfrak{S}_1^T n\cdot n) )(\theta\cdot n) ds.
$$
Therefore \varepsilonqref{eq:boundary_exp} reads
\begin{equation}\label{eq:boundary_exp_2}
\begin{split}
	dJ(\mathcal{B}ega)(\theta) 
	= & \int_{\partial \mathcal{B}ega} ((\mathbf{S}_1^+-\mathbf{S}_1^-)n\cdot n) (\theta\cdot n) + (\mathfrak{S}_0\cdot n) (\theta\cdot n)\, ds\\
	& + \int_{\partial \mathcal{B}ega} \mathfrak{S}_1\cdot D_\Gamma n (\theta\cdot n) + (- \operatorname{div}_\Gamma (\mathfrak{S}_1^T n) + \mathcal{H}(\mathfrak{S}_1^T n\cdot n) )(\theta\cdot n)\, ds,
\varepsilonnd{split}
\varepsilonnd{equation}
which we can rewrite as \varepsilonqref{eq:general_boundary_exp}.
\varepsilonnd{proof}

\begin{remark}
In Proposition \ref{tensor_relations}, if $\mathfrak{S}_1\varepsilonquiv 0$, one can still obtain \varepsilonqref{eq:equvilibrium_strong} when $\mathcal{B}ega$ is only Lipschitz instead of $C^1$. 
\varepsilonnd{remark}

\begin{corollary}\label{corollary2}
	Let the assumptions of Proposition~\ref{tensor_relations} be satisfied. Suppose that the tensor 
	$\mathfrak{S}_1:\partial\mathcal{B}ega\rightarrow \mathcal L(\mathbf{R}^d,\mathbf{R}^d)$ has the form $\mathfrak{S}_1 = \alpha(I -n\otimes n)$, where $\alpha\in C^0(\partial\mathcal{B}ega)$. Then 
	\varepsilonqref{eq:general_boundary_exp} simplifies to 
	\begin{equation}
	\label{shape_der_simplified}
	dJ(\mathcal{B}ega)(\theta) = \int_{\partial \mathcal{B}ega} g_1 \, \theta\cdot n\, \, ds, 
	\varepsilonnd{equation}
	where $g_1$ is given by
	$$ g_1 :=  [(\mathbf{S}_1^+-\mathbf{S}_1^-)n]\cdot n + \mathfrak{S}_0\cdot n + \alpha\mathcal{H}. $$
\varepsilonnd{corollary}

\begin{proof}
First  $\mathfrak{S}_1^T = \alpha(I -n\otimes n)^T = \mathfrak{S}_1$ and $ \mathfrak{S}_1^T n = \alpha(I - n\otimes n)n = \alpha(n- (n\cdot n)n) = 0$, thus the two last terms in \varepsilonqref{g1} vanish. Concerning the third term in \varepsilonqref{g1} we write
\begin{align*} 
\mathfrak{S}_1\cdot D_\Gamma n = \alpha (I -n\otimes n)\cdot D_\Gamma n 
& = \alpha(\operatorname{tr}(D_\Gamma n) - (n\otimes n)\cdot D_\Gamma n )=\alpha(\operatorname{div}_\Gamma n - (D_\Gamma n n)\cdot n)= \alpha\mathcal{H}.
\varepsilonnd{align*}
where we have used $(D_\Gamma n n)\cdot n=0$. 
\varepsilonnd{proof}

\begin{remark}
	The particular tensor $\mathfrak{S}_1 = \alpha(I -n\otimes n)$ in Corollary \ref{corollary2} is commonly encountered in shape optimization problems. In fact, \varepsilonqref{shape_der_simplified} corresponds to a standard formula that can be found in most textbooks on shape optimization.
\varepsilonnd{remark}

\begin{remark}
	Recall that for given vector fields $\theta, \zeta$, the second order shape derivative is defined 
	by $$d^2J(\mathcal{B}ega)(\theta)(\zeta) := \dt dJ(\Phi_t^\zeta(\mathcal{B}ega))(\theta)|_{t=0}.$$
	Once we have identified a tensor representation \varepsilonqref{ea:volume_from} for the shape derivative $dJ(\mathcal{B}ega)(\theta)$ for fixed $\theta$, it is convenient to differentiate it to also obtain a tensor representation for the second-order shape derivative. Further, Proposition \ref{tensor_relations} or Corollary \ref{corollary2} can also be applied to obtain a boundary expression for the second order shape derivative.
\varepsilonnd{remark}

Similar relations as in Proposition \ref{tensor_relations} could be obtained for any tensor representation of order $k$. For instance in the case $k = 2$ we obtain the relations

\begin{equation} \label{eq:equvilibrium_strong_second}
\begin{split}
	\mathcal A\mathbf \mathbf{S}_2^+-\operatorname{div}(\mathbf{S}_1^+) + \mathbf{S}_0^+ &= 0 \quad \text{ in } \mathcal{B}ega, \\
	\mathcal A\mathbf \mathbf{S}_2^- -\operatorname{div}(\mathbf{S}_1^-) + \mathbf{S}_0^- &= 0 \quad \text{ in } D\setminus\overline \mathcal{B}ega,
\varepsilonnd{split}
\varepsilonnd{equation}
where $(\mathcal A\mathbf{S}_2)_l = \sum_{i,j=1}^d \partial_{x_ix_j}(\mathbf{S}_2)_{ijl}$. 

Using the averaged adjoint approach from Theorem  \ref{thm:sturm} yields  the tensor representation \varepsilonqref{ea:volume_from}  of the shape derivative.
	Then Proposition \ref{tensor_relations} can be used to immediately derive the standard boundary expression of the shape gradient from this  tensor representation.

\subsection{Examples of tensor representations}

In this section we present several examples of representations corresponding to Definition \ref{def:tensor} and apply the observations from  Section \ref{def_prop}.

\subsubsection*{First order tensor representation}
A basic example of a first order tensor representation of the shape derivative is for 
$$J(\mathcal{B}ega) = \int_\mathcal{B}ega f \, dx + \int_{\partial \mathcal{B}ega} g \, ds$$ 
with $f,g\in C^2(\mathbf{R}^d)$. Then one easily computes
$$ dJ(\mathcal{B}ega)(\theta) = \int_\mathcal{B}ega \nabla f \cdot \theta + f\operatorname{div}(\theta)\,dx + \int_{\partial \mathcal{B}ega} \nabla g\cdot \theta + g\operatorname{div}_\Gamma \theta\, ds . $$
The corresponding tensor representation \varepsilonqref{ea:volume_from} is
\begin{align*}
\mathbf{S}_1^+:= f  &I,\quad \mathbf{S}_1^-:= 0 , \quad \mathbf{S}_0^+:= \nabla f, \quad \mathbf{S}_0^-:= 0, \quad
\mathfrak{S}_1:= g  (I-n\otimes n),\quad  \mathfrak{S}_0 := \nabla g.
\varepsilonnd{align*}
Note that $\mathfrak{S}_1$ has the form assumed in Corollary \ref{corollary2}. Applying this Corollary, assuming the domain has enough regularity, we obtain in view of \varepsilonqref{shape_der_simplified} the classical formula:
$$ dJ(\mathcal{B}ega)(\theta) =\int_{\partial\mathcal{B}ega} g_1 \, \theta\cdot n\, \, ds, 
$$
where $g_1$ is given by
$$ g_1 :=  f + \partial_n g + g\mathcal{H}. $$
Note that in the particular case $f=0$ we have obtained as a byproduct the formula
\begin{equation}\label{tangential_green}
\int_{\partial \mathcal{B}ega} \nabla g\cdot \theta + g\operatorname{div}_\Gamma \theta\, ds  =\int_{\partial\mathcal{B}ega} (\partial_n g + g\mathcal{H}) \, \theta\cdot n\, \, ds, 
\varepsilonnd{equation}
and when in addition $\partial_n g=0$ or $g$ is defined only on $\partial\mathcal{B}ega$, \varepsilonqref{tangential_green} becomes the classical {\it tangential Green's formula}; see for instance \cite[proposition 5.4.9]{MR2512810}. 

\subsubsection*{Non-homogeneous Dirichlet problem}

The following problem was already considered for instance in \cite{MR2731611}. Here we present a fairly 
easy way 

to compute the shape derivative. Let $\mathcal{B}ega$ be an open and bounded subset of $\mathbf{R}^d$ 
that is contained in an open and bounded set $D$. Consider
\begin{align}\label{eq:state_non_homo}
-\Delta v &= f\mbox{ in }\mathcal{B}ega,\\
v & = g \mbox{ on }\partial\mathcal{B}ega,\label{eq:state_non_homo_}
\varepsilonnd{align}
where $f\in L^2(D)$ and $g\in H^2(D)$. Introducing the 
variable $u:=v-g$, we observe that \varepsilonqref{eq:state_non_homo}-\varepsilonqref{eq:state_non_homo_}
is equivalent to the homogeneous Dirichlet problem
\begin{align}\label{eq:state_non_homo-2}
-\Delta u &= \Delta g + f\mbox{ in }\mathcal{B}ega,\\\label{eq:state_non_homo-2_}
u & = 0 \mbox{ on }\partial\mathcal{B}ega.
\varepsilonnd{align}
Consider the cost function
\begin{equation}\label{eq:cost_dirichlet_non}
J(\mathcal{B}ega)=\int_\mathcal{B}ega |v-u_d|^2\,dx=\int_\mathcal{B}ega |u+g-u_d|^2\,dx.
\varepsilonnd{equation}
The weak formulation of \varepsilonqref{eq:state_non_homo-2},\varepsilonqref{eq:state_non_homo-2_} reads:  
\begin{equation}\label{eq:non_dir}
\mbox{Find } u\in H^1_0(\mathcal{B}ega):  \int_\mathcal{B}ega \nabla u\cdot \nabla \psi \, dx = \int_\mathcal{B}ega  -\nabla g\cdot \nabla \psi +  f\psi\, dx \quad \text{ for all } \psi \in H^1_0(\mathcal{B}ega).
\varepsilonnd{equation}
Note that the previous weak formulation is already well-defined for arbitrary open and bounded 
set $\mathcal{B}ega$. We do not need to impose any regularity on $\mathcal{B}ega$.
The perturbed problem of the previous equation, which is obtained by considering \varepsilonqref{eq:non_dir} on $\Phi_t(\mathcal{B}ega)$ and
performing a change of variables, reads: find $u^t\in H^1_0(\mathcal{B}ega)$ such that 
\begin{equation}\label{eq:non_dir_per}
\int_\mathcal{B}ega A(t)\nabla u^t\cdot \nabla \psi \, dx = \int_\mathcal{B}ega - A(t)\nabla g^t\cdot \nabla \psi +  \xi(t)f^t \psi\, dx \quad \text{ for all } \psi \in H^1_0(\mathcal{B}ega),
\varepsilonnd{equation}
where $\xi(t) := \det(D\Phi_t)$ and $A(t) :=\xi(t)D\Phi_t^{-1}D\Phi_t^{-T}$. The following continuity result is standard:

\begin{lemma}
	There exists a constant $c>0$ such that $ \|u^t-u^0\|_{H^1_0(\mathcal{B}ega)} \le ct$ for all $ t\in [0,\tau]. $
\varepsilonnd{lemma}

Introduce 
\begin{align*}
a(t,\varphi,\psi) & :=  \int_\mathcal{B}ega A(t)\nabla \varphi\cdot \nabla \psi \, dx + \int_\mathcal{B}ega A(t)\nabla g^t\cdot \nabla \psi  - \xi(t)f^t \psi\, dx  \\
b(t,\varphi) &  :=\int_\mathcal{B}ega \xi(t) |\varphi +g^t-u_d^t|^2\,dx.
\varepsilonnd{align*}
where $g^t:=g\circ \Phi_t$ and $u_d^t:=u_d\circ \Phi_t$.

Recall that the associated Lagrangian \varepsilonqref{G_lag} is
$G(t,\varphi,\psi) = a(t,\varphi,\psi) + b(t,\varphi).
$
The averaged adjoint equation \varepsilonqref{averated_} reads 
$$\int_\mathcal{B}ega A(t)\nabla \varphi\cdot \nabla p^t \, dx = \int_\mathcal{B}ega (u^t+u^0 + 2g^t - 2u_d^t)\varphi\,dx \quad \text{ for all } \varphi\in H^1_0(\mathcal{B}ega). $$

The following continuity result for the adjoint is standard:
\begin{lemma}
	There exists a constant $c>0$ such that 
	$$ \|p^t-p^0\|_{H^1(\mathcal{B}ega)} \le ct \quad \text{ for all } t\in [0,\tau].$$
\varepsilonnd{lemma}
One readily verifies that all assumptions of Theorem~\ref{thm:sturm} are satisfied, except for (H1) which we now prove.
	Indeed using $p^t\rightarrow p^0$ in $H^1_0(\mathcal{B}ega)$ as $t$ goes  to zero and the strong differentiability of $t\mapsto A(t)$ and $t\mapsto \xi(t)$, we get 
	\begin{align*}
	&\lim_{t\searrow 0 }  \textnormal{fr}ac{ G(t,u^0,p^t) - G(0,u^0,p^t)}{t}   \\
	&\hspace{1cm} 
	=\lim_{t\searrow 0} \bigg(\int_{\mathcal{B}ega} \left(\textnormal{fr}ac{ A(t) - I }{t} \right) \nabla u^0 \cdot \nabla p^t \, dx + \int_{\mathcal{B}ega} \left(\textnormal{fr}ac{ A(t)\nabla g^t - \nabla g }{t}\right) \cdot \nabla p^t {\color{red}\, dx} - 
	\left(\textnormal{fr}ac{\xi(t)f^t - f}{t} \right) p^t \, dx \\
	& \hspace{1cm} + \int_{\mathcal{B}ega} \textnormal{fr}ac{ \xi(t)| u^0 + g^t - u_d^t |^2  - | u^0 + g - u_d  |^2}{t}\bigg)\\
	&\hspace{1cm} = \partial_t G(0,u^0,p^0),
	\varepsilonnd{align*}
	which shows that (H1) is satisfied.

Hence, applying Theorem~\ref{thm:sturm} yields
$$ dJ(\mathcal{B}ega)(\theta) = \partial_t a(0,u,p) + \partial_tb(0,u) ,$$
which is by definition equivalent to 
\begin{equation}
\begin{split}
	dJ(\mathcal{B}ega)(\theta)    = &\int_\mathcal{B}ega A'(0)(\nabla u \cdot \nabla p + \nabla g\cdot \nabla p) \, dx + \int_\mathcal{B}ega \nabla(\nabla g\cdot \theta)\cdot \nabla p -  \operatorname{div}(f\theta) p \, dx\\
	& + \int_\mathcal{B}ega \operatorname{div}(\theta) |u+g-u_d|^2 + (\nabla (g-u_d)\cdot \theta))(u+g-u_d)\,dx.
\varepsilonnd{split}
\varepsilonnd{equation}

Since $A'(0) = (\operatorname{div} \theta) I - D\theta^T - D\theta$ we obtain the tensor representation \varepsilonqref{eq:first_order_tensor} with:
\begin{align*} 
\mathbf{S}_1 &= I(\nabla u\cdot \nabla p  + \nabla g \cdot \nabla p - fp + |u+g-u_d|^2)  -\nabla u\otimes \nabla p - \nabla p \otimes \nabla u - \nabla p \otimes \nabla g, \\
\mathbf{S}_0 &= D^2 g\nabla p - p\nabla f  
+(u+g-u_d)  \nabla (g-u_d),   \\
\mathfrak{S}_1 &= 0, \quad \mathfrak{S}_0 = 0. 
\varepsilonnd{align*}
Now applying  \varepsilonqref{shape_der_simplified} we get immediately
\begin{align*}
g_1 = \nabla u\cdot \nabla p +\nabla  g \cdot \nabla p - fp + |u+g-u_d|^2  -2\partial_n u\partial_n p - \partial_n p \partial_n g.
\varepsilonnd{align*}
Using the definition of the tangential gradient and $p=0,u=0$ on $\Gamma$ implies $\nabla_\Gamma u=\nabla_\Gamma p=0$, so we obtain the simpler expression 
\begin{align*}
g_1 = -\partial_n u\partial_n p  + |u+g-u_d|^2 = -\partial_n (v-g)\partial_n p  + |v-u_d|^2 .
\varepsilonnd{align*}
Finally, substituting back $u=v-g$ we obtain the formula
$$dJ(\mathcal{B}ega)(\theta) = \int_{\partial \mathcal{B}ega} (-\partial_n (v-g)\partial_n p  + |v-u_d|^2)\, \theta\cdot n\, ds  .$$
This formula can be found for instance in \cite[p. 566, Formula 6.38]{MR2731611}, where
the adjoint has the sign opposite to our case.

\subsubsection*{Elliptic problem: first order tensor representation}

Suppose that $\mathcal{B}ega \subset D\subset \mathbf{R}^d$ is a smooth bounded domain, where $D\subset \mathbf{R}^d$ is the smooth ``universe''. 
Let us  consider the Dirichlet problem:
\begin{align}\label{eq:state}
\begin{split}
-\operatorname{div}(M\nabla u) + u&= f\mbox{ in }\mathcal{B}ega,\\
u & = 0 \mbox{ on }\partial\mathcal{B}ega,
\varepsilonnd{split}
\varepsilonnd{align}
where $M\in \mathbf{R}^{d,d}$ is a positive definite matrix.
Consider the cost function
\begin{equation}\label{eq:cost_dirichlet}
J(\mathcal{B}ega)=\int_\mathcal{B}ega |u-u_d|^2\,dx,
\varepsilonnd{equation}
where $u_d\in C^1(\mathbf{R}^d)$. Let us introduce 
\begin{align*}
a(t,\varphi,\psi) & :=  \int_\mathcal{B}ega (M Q^t\nabla \varphi \cdot Q^t\nabla \psi   + \varphi \psi )\xi(t) \, dx - \int_\mathcal{B}ega \xi(t) f^t\psi \, dx \\
b(t,\varphi) &  :=\int_\mathcal{B}ega \xi(t) |\varphi-u_d^t|^2\,dx,
\varepsilonnd{align*}
where $Q^t:=D\Phi_t^{-T}$ and $\xi(t) := \det(D\Phi_t)$.
Then the weak formulation of \varepsilonqref{eq:state} on the perturbed domain $\mathcal{B}ega_t$, once transported back to $\mathcal{B}ega$ is
\begin{equation}n
a(t,u^t,\psi) = 0 \quad \text{ for all } \psi\in H^1_0(\mathcal{B}ega).
\varepsilonnd{equation}n
The Lagrangian corresponding to the minimization of $J(\mathcal{B}ega)$ and the PDE constraint \varepsilonqref{eq:state} is 
\begin{equation}\label{eq:lagrange}
G(t,\varphi,\psi) = b(t,\varphi) + a(t,\varphi,\psi). 
\varepsilonnd{equation}

 It can be shown using Theorem \ref{thm:sturm}  that $dJ(\mathcal{B}ega)(\theta) =  \partial_t G(0,u,p)$, where $p\in H^1_0(\mathcal{B}ega)$ denotes the adjoint state:
\begin{equation}\label{eq:adjoint_state_simple_dirichlet}
\begin{split}
	\int_\mathcal{B}ega M\nabla \psi\cdot \nabla p +p\psi\, dx=-\int_\mathcal{B}ega 2 (u-u_d)\psi\,  dx\quad \text{ for all } \psi \in H^1_0(\mathcal{B}ega) .
\varepsilonnd{split}
\varepsilonnd{equation}
\\
The tensor representation \varepsilonqref{ea:volume_from} of the shape derivative of $J(\mathcal{B}ega)$ in direction $\theta\in C^2_c(D,\mathbf{R}^d)$ is  given by
\begin{equation}\label{eq:shape_derivative}
dJ(\mathcal{B}ega)(\theta) = \int_\mathcal{B}ega \mathbf S_1\cdot  D\theta +  \mathbf{S}_0\cdot \theta \, dx.
\varepsilonnd{equation}
where we use the relation $(\nabla p \otimes \nabla u)\cdot D\theta = D\theta \nabla u\cdot \nabla p$ to get the tensors 
\begin{align}
\label{f_mathbf}
\mathbf{S}_0 & = -2(u-u_d)\nabla u_d - p\nabla f ,\\
\mathbf{S}_1 & = - \nabla p \otimes M \nabla u -  \nabla u \otimes M^T\nabla p + (M\nabla u\cdot \nabla p+ up-fp+(u-u_d)^2)I.  
\varepsilonnd{align}

In the simple case where $M=I$, assuming $u,p\in C^2(\overline \mathcal{B}ega)$, we know from the previous discussion that \varepsilonqref{eq:equvilibrium_strong} 
 is satisfied. 
Noting that 

\begin{align*}
\operatorname{div}(\nabla p\otimes \nabla u) &= \Delta u \nabla p + (D^2p)^T\nabla u, \\
\operatorname{div}(\nabla u\otimes \nabla p) &= \Delta p \nabla u + (D^2u)^T\nabla p, \\ 
\nabla (\nabla u\cdot \nabla p) &= D^2u\nabla p +  D^2 p \nabla u,
\varepsilonnd{align*}
the relation
\begin{equation}\label{nes}
- \operatorname{div}(\mathbf{S}_1) + \mathbf{S}_0 = 0 \quad \text{ in } \mathcal{B}ega
\varepsilonnd{equation}
is equivalent to
$$ (-\Delta u + u -f)\nabla p + (-\Delta p + p + 2(u-u_d))\nabla u = 0  \quad \text{ in } \mathcal{B}ega.$$ 
Therefore, we observe that the fundamental relation \varepsilonqref{eq:equvilibrium_strong} between the tensors $\mathbf{S}_1$ and $\mathbf{S}_0$  corresponds to the  strong solvability of the state and adjoint state equation.

\section{Descent directions}\label{section3}

In this paper we are interested in numerical methods for shape optimization problems of the type
\begin{equation}
\min_{\mathcal{B}ega \in \mathfrak{P}}J(\mathcal{B}ega),
\varepsilonnd{equation}
where $\mathfrak{P}\subset\mathcal{P}(D)$ is the admissible set. Assume  $J:\mathfrak{P}\rightarrow \mathbf{R}$ is shape differentiable at $\mathcal{B}ega\subset D\subset \mathbf{R}^d$.

\begin{definition}[descent direction]\label{def:descent}
	The vector field $\theta\in C^{0,1}_c(D,\mathbf{R}^d)$ is called a {\it descent direction} for $J$ at $\mathcal{B}ega$ if there exists an ${\varepsilon}$ such that
	$$ J(\Phi_t^\theta(\mathcal{B}ega))< J(\mathcal{B}ega)\mbox{ for all } t\in (0,{\varepsilon}).$$
	If the Eulerian semiderivative of $J$  at $\mathcal{B}ega$ in direction $\theta$ exists and if it is a descent direction then by definition
	\begin{equation}\label{eq:descent}
	dJ(\mathcal{B}ega)(\theta) < 0.
	\varepsilonnd{equation}
\varepsilonnd{definition}

Descent directions are used in iterative methods for finding approximate (possibly local) minimizers of  $J(\mathcal{B}ega)$. Typically, at a given starting point $\mathcal{B}ega$, one determines a descent direction $\theta$ and proceeds along this direction as long as the cost functional $J$ reduces sufficiently using a step size strategy.

In this section we give a general setting for computing descent directions in the framework of gradient methods using the domain and boundary representations of the shape derivative according to Theorem \ref{thm:structure_theorem}. We show how a descent direction $\theta$ with any regularity $H^s$, $s\geq 1$ can be obtained by solving an appropriate partial differential equation. 
We also show how to deal with bound constraints on $\theta$. 

In order to develop a setting allowing to define general descent directions, we recall sufficient conditions for the solvability of the following operator equation
\begin{equation}n
A\theta=f,
\varepsilonnd{equation}n
where $A:\textbf{E}\rightarrow \textbf{E}'$ is an operator between a Banach space $\textbf{E}$ and its dual $\textbf{E}'$.  Sufficient conditions
for the bijectivity of $A$  are given by the theorem of Minty-Browder \cite[p.364, Theorem 10.49]{MR2028503}.

\begin{theorem}[Minty-Browder]\label{MB}
	Let $(\textbf{E};\|\cdot\|_{\textbf{E}})$ be a reflexive separable Banach space and $A:\textbf{E}\rightarrow \textbf{E}'$ a bounded, hemi-continuous, monotone and coercive operator. Then $A$ is surjective, i.e. for each $f\in \textbf{E}'$ there exists $\theta\in \textbf{E}$ such that $A\theta=f$.
	Moreover if $A$ is strictly monotone then it is bijective.
\varepsilonnd{theorem}

Let $A:\textbf{E}\rightarrow \textbf{E}'$ be an operator on a reflexive, separable Banach space $\textbf{E}$ satisfying the assumptions of  Theorem \ref{MB} with $A(0)\theta\ge 0$ for $v\in \textbf{E}$. 
Assume $dJ(\mathcal{B}ega)$ can be extended to $\textbf{E}'$ if necessary; for simplicity we keep the same notation for the extension. Introduce the bilinear form

	\begin{equation}\label{eq:bilinear}
	\mathcal B:\textbf{E}\times \textbf{E}\rightarrow \mathbf{R},\qquad \mathcal B(\theta,\zeta):=\langle A\theta,\zeta \rangle_{\textbf{E}',\textbf{E}}.
	\varepsilonnd{equation}

Consider the variational problem:
\begin{equation}
\label{VP_1}
(\text{VP})\qquad
\mbox{Find  } \theta_1\in \textbf{E}\mbox{  such that  }  \mathcal B(\theta_1,\zeta)=- dJ(\mathcal{B}ega)(\zeta) \mbox{  for all  } \zeta\in \textbf{E} ,
\varepsilonnd{equation}
Then the solution $\theta_1$ of (VP) is a descent direction since $dJ(\mathcal{B}ega)(\theta_1) = -\mathcal B(\theta_1,\theta_1)\leq 0$.

In certain situations it is desirable to have bound constraints on the shape perturbation. This may be handled by considering the more general case of a variational inequality. Given a subset $K\subset \textbf{E}$ with $0\in K$, consider the variational inequality:
\begin{equation}n
(\text{VI})\qquad
\mbox{Find  }\theta_2\in K\mbox{  such that  } \mathcal B(\theta_2,\theta_2-\zeta) \le  dJ(\mathcal{B}ega)(\zeta-\theta_2) \mbox{  for all  } \zeta\in K .
\varepsilonnd{equation}n
The solution $\theta_2$ of  (VI) yields a descent direction for $J$ at $\mathcal{B}ega$ since taking $\zeta = 0\in K$ we get 
$$dJ(\mathcal{B}ega)(\theta_2) \le -\mathcal B(\theta_2,\theta_2)\leq 0.$$ In view of Theorem \ref{thm:structure_theorem}, we choose $\textbf{E}\subset H^{s}(D)$ where
$s$ is such that $dJ(\mathcal{B}ega):H^s(D,\mathbf{R}^d)\rightarrow \mathbf{R}^d$ is continuous. When $\textbf{E}$ is a Hilbert space, one may identify $\textbf{E}'$ with $\textbf{E}$. Therefore if $\mathcal B$ is bilinear, coercive, and continuous, then Lax Milgram's lemma ensures that (VP) has a unique solution. For all other cases we may have to use Theorem \ref{MB} or similar results.

\begin{remark}
		\begin{itemize}
			\item[(a)]
			Let $\textbf{E}:= H^1_0(D,\mathbf{R}^d)$, 
			$\mathcal B(\theta, \zeta) := \int_D D\theta:D\zeta \,dx $ and $\mathcal{B}ega^+\Subset D$. Then 
			\varepsilonqref{VP_1} reads: find $\theta \in H^1_0(D,\mathbf{R}^d)$, such that
			$ \mathcal B(\theta, \zeta) = - dJ(\mathcal{B}ega^+)(\zeta)$ for all $\zeta\in H^1_0(D,\mathbf{R}^d)$. Under the assumption that $\partial \mathcal{B}ega^+ \in C^2$, $\theta|_{\mathcal{B}ega^+} \in H^2(\mathcal{B}ega^+)$,
			and $\theta|_{D\setminus \overline{\mathcal{B}ega^+}} \in H^2(D\setminus \overline{\mathcal{B}ega^+})$, Proposition \ref{tensor_relations} yields
			$$ \int_{\partial \mathcal{B}ega^+} g \; \zeta\cdot n \, ds = dJ(\mathcal{B}ega^+)(\zeta) \quad \text{ for all }\zeta\in H^1_0(D,\mathbf{R}^d), \quad \text{ where } \quad g= -(D\theta^+ n -D\theta^- n)\cdot n. $$ 
			This shows that the  restriction to $\partial \mathcal{B}ega^+$  of the obtained descent direction $\theta$ is 
			more regular than the function $g$.	

			\item[(b)]  Let $n$ be an extension of the unit normal to $\mathcal{B}ega^+$ in  $D$. If $\theta$ defined on $D$ is a descent direction then $(\theta\cdot n)n|_{\partial\mathcal{B}ega^+}$ is also a descent direction, for the tangential part of $\theta$ does not influence the derivative. Indeed define  $\theta_{\tau}:= \theta - (\theta\cdot n)n$, then by Nagumo's theorem $dJ(\mathcal{B}ega^+)(\theta_{\tau}) = 0$ and thus $dJ(\mathcal{B}ega^+)(\theta) = dJ(\mathcal{B}ega^+)((\theta\cdot n)n)$.
			However, $\theta$ and $(\theta\cdot n)n$ lead  to different transformations of the domains in general, indeed the tangential term actually has an influence for large deformations,  which means $\Phi_t^\theta(\mathcal{B}ega^+)  \ne \Phi_t^{(\theta \cdot n)n}(\mathcal{B}ega^+)$. This influence appears for instance in the shape Hessian.
		\varepsilonnd{itemize}
	\varepsilonnd{remark}
	
	
	
	
	
	

 	\section{Electrical impedance tomography}\label{section4}

	We consider an application of the results above to a typical and important interface problem: the inverse problem of electrical impedance tomography (EIT) also known as the inverse conductivity or Calder\'on's problem \cite{MR590275} in the mathematical literature. It is an active field of research with an extensive literature; for further details we point the reader toward the survey papers \cite{MR1955896,MR1669729} as well as \cite{MR2986262} and the references therein. We consider the particular case where the objective is to reconstruct a piecewise constant  conductivity $\sigma$ which amounts to determine an interface ${\Gamma^+}$ between some inclusions and the background. We refer the reader to \cite{MR2329288, MR2132313, MR2536481,HintLauNov,MR2886743, 0266-5611-31-7-075009,Canelas2013} for more details on this approach. 
	
	The main interest of studying EIT is to apply the approach developed in this paper to a problem which epitomizes general interface problems and simultaneously covers the entire spectrum of difficulties encountered with severely ill-posed inverse problem.

	\subsection{Problem statement}\label{sec:eit}

	Let $D\subset \mathbf{R}^d$ be a Lipschitz domain, and $\mathcal{B}ega^+,\mathcal{B}ega^-\subset D$ open sets such that $D=\mathcal{B}ega^+\cup \mathcal{B}ega^-\cup \Gamma^+$, where ${\Gamma^+}=\partial \mathcal{B}ega^+=\overline{\mathcal{B}ega^+}\cap \overline{\mathcal{B}ega^-}$ and $\Gamma = \partial D=\partial\mathcal{B}ega^-\setminus\Gamma^+$; see Figure \ref{partition}. In this section $n$ denotes either the outward unit normal vector to $D$ or the outward unit normal vector to $\mathcal{B}ega^+$. Decompose $\Gamma$ as $\Gamma = \Gamma_d\cup\Gamma_n$.  Let $\sigma=\sigma^+\chi_{\mathcal{B}ega^+}+\sigma^-\chi_{\mathcal{B}ega^-}$ where $\sigma^\pm$ are scalars and $f=f^+\chi_{\mathcal{B}ega^+}+f^-\chi_{\mathcal{B}ega^-}$ where $f^\pm\in H^1(D)$.

	\begin{figure}
		\begin{center}
			\includegraphics[width=0.3\textwidth]{./EIT_sketch_small2}
			\caption{Partition $D=\mathcal{B}ega^+\cup \mathcal{B}ega^-\cup \Gamma$. }\label{partition}
		\varepsilonnd{center}
	\varepsilonnd{figure}

	Consider the following problems: find $u_n\in H^1_d (D)$ such that
	\begin{equation}\label{varun}
	\int_D \sigma \nabla u_n \cdot \nabla z = \int_D fz + \int_{\Gamma_n} gz\ \mbox{ for all }\ z\in H^1_d (D) 
	\varepsilonnd{equation}
	and find $\dot{\mathbf{u}}i\in H^1_{dn}(D)$ such that 
	\begin{equation}\label{varud}
	\int_D \sigma \nabla \dot{\mathbf{u}}i \cdot \nabla z = \int_D fz\  \mbox{ for all }\ z\in H^1_0(D) 
	\varepsilonnd{equation}
	where
	\begin{align*}
	H^1_d (D)&:=\{v\in H^1(D)\ |\ \; v=0\mbox{ on }\Gamma_d\}, \\
	H^1_{dn}(D)&:=\{v\in H^1(D)\ |\ \; v=0\mbox{ on }\Gamma_d, v=h\mbox{ on }\Gamma_n\},\\ 
	H^1_{0}(D)&:=\{v\in H^1(D)\ |\ \; v=0\mbox{ on }\Gamma\}
	\varepsilonnd{align*}
	and $g\in H^{-1/2}(\Gamma_n)$ represents the input, in this case the electric current applied on the boundary and $h\in H^{1/2}(\Gamma_n)$ is the measurement of the potential on $\Gamma_n$, or the other way around, i.e. $h$ can be the input and $g$ the measurement.  Define also the space
	$$PH^k(D) := \{ u = u^+\chi_{\mathcal{B}ega^+} + u^-\chi_{\mathcal{B}ega^-}|\ u^+\in H^k(\mathcal{B}ega^+),\ u^-\in H^k(\mathcal{B}ega^-) \}. $$
	Consider the following assumption which will be used only to derive the boundary expression of the shape derivative but is not required for the domain expression:
	\begin{assumption}\label{assump1}
		The domains $D, \mathcal{B}ega^+,\mathcal{B}ega^-$ are of class $C^k$, $f\in PH^{\max(k-2,1)}(D)$, $g\in H^{k-\textnormal{fr}ac{3}{2}}(D)$ and  $h\in H^{k-\textnormal{fr}ac{1}{2}}(D)$ for $k\geq 2$. 
	\varepsilonnd{assumption}
	Applying Green's formula under Assumption \ref{assump1},  equations \varepsilonqref{varun} and \varepsilonqref{varud}  are equivalent to the following transmission problems where $u_n = u_n^+\chi_{\mathcal{B}ega^+} + u_n^-\chi_{\mathcal{B}ega^-}$ and $\dot{\mathbf{u}}i = \dot{\mathbf{u}}i^+\chi_{\mathcal{B}ega^+} + \dot{\mathbf{u}}i^-\chi_{\mathcal{B}ega^-}$:
	\begin{align}
	\label{eit1.1_}-\sigma^+\Delta u_n^+ & = f \mbox{ in }\mathcal{B}ega^+, \quad  -\sigma^-\Delta u_n^- = f \mbox{ in }\mathcal{B}ega^-, \\
	\label{eit1.2_}  u_n^- & = 0 \mbox{ on }\Gamma_d ,\\
	\label{eit1.3_}\sigma^- \partial_n u_n^- & = g \mbox{ on }\Gamma_n ,
	\varepsilonnd{align}
	\begin{align}
	\label{eit2.1_}-\sigma^+\Delta \dot{\mathbf{u}}i^+ & = f \mbox{ in }\mathcal{B}ega^+,\quad  -\sigma^-\Delta \dot{\mathbf{u}}i^- = f \mbox{ in }\mathcal{B}ega^-,\\
	\label{eit2.2_}\dot{\mathbf{u}}i^- & = 0 \mbox{ on }\Gamma_d, \\
	\label{eit2.3_}\dot{\mathbf{u}}i^- & = h \mbox{ on }\Gamma_n, 
	\varepsilonnd{align}
	with the transmission conditions
	\begin{equation}
	\begin{split}
		\sigma^+\partial_nu_n^+=\sigma^-\partial_nu_n^-, &  \qquad \sigma^+\partial_n\dot{\mathbf{u}}i^+=\sigma^-\partial_n\dot{\mathbf{u}}i^- \quad \text{ on } {\Gamma^+},\\ 
		u_n^+=u_n^-,     &\qquad  \dot{\mathbf{u}}i^+=\dot{\mathbf{u}}i^-\quad \text{ on } {\Gamma^+}. \varepsilonnd{split}
	\varepsilonnd{equation}

	On $\Gamma_d$ we impose homogeneous Dirichlet conditions, meaning that the voltage is fixed and no measurement is performed. One may take $\Gamma_d =\varepsilonmptyset$, in which case   \varepsilonqref{varun} becomes a pure Neumann problem and  additional care must be taken for the uniqueness and existence of a solution. The situation $\Gamma_d \neq \varepsilonmptyset$ corresponds to partial measurements. Alternatively, it is also possible to consider a slightly different problem where each function $u_n$ and $\dot{\mathbf{u}}i$ has both the boundary conditions \varepsilonqref{eit1.3_} and \varepsilonqref{eit2.3_} on different parts of the boundary. 
	
	Several measurements can be made by choosing sets of functions $\{g_i\}_{i=1}^I$ and $\{h_i\}_{i=1}^I$. 
	Writing $u_{n,i}$ and $u_{d,i}$ for the corresponding states, the problem of electrical impedance tomography is: 
		\begin{equation}\label{EIT_pb}
		\mbox{(EIT): Given $\{g_i\}_{i=1}^I$ and $\{h_i\}_{i=1}^I$,  find $\sigma$ such that $u_{n,i} = u_{d,i}$ in $D$ for $i=1,..,I$.}
		\varepsilonnd{equation}
		Note that $u_{n,i} = u_{n,i}(\mathcal{B}ega^+)$ and $u_{d,i} = u_{d,i}(\mathcal{B}ega^+)$ actually depend on $\mathcal{B}ega^+$ through $\sigma=\sigma(\mathcal{B}ega^+)$, however we  often write $u_{n,i}$ and $u_{d,i}$ for simplicity. In this section, we assume that the conductivities $(\sigma^+,\sigma^-)$ are known, therefore the EIT problem \varepsilonqref{EIT_pb} reduces to the following shape optimization problem where $\mathcal{B}ega^+$ is the unknown. 

		\begin{align}\label{EIT-SO}
		\begin{split}
		\mbox{Given }& \mbox{$\{g_i\}_{i=1}^I$, $\{h_i\}_{i=1}^I$ and $(\sigma^+,\sigma^-)$ with $\sigma=\sigma^+\chi_{\mathcal{B}ega^+}+\sigma^-\chi_{\mathcal{B}ega^-}$,} \\
		& \mbox{find $\mathcal{B}ega^+$ such that $u_{n,i} = u_{d,i}$ in $D$ for $i=1,..,I$.}
		\varepsilonnd{split}
		\varepsilonnd{align}

		Actually, the result for several measurements can be straightforwardly deduced from the case of one measurement by summing  the cost functionals corresponding to each measurement, therefore in this section we stick to the  case $I=1$ of one measurement $g$ for simplicity of presentation. In section \ref{section6} we consider several measurements for the numerics.
	
	
	The notion of well-posedness due to Hadamard requires the existence and uniqueness of a solution and the continuity of the inverse mapping. The severe ill-posedness of EIT is well-known: uniqueness and continuity of the inverse mapping depend on the regularity of $\sigma$, the latter being responsible for the instability of the reconstruction process. Additionally, partial measurements often encountered in practice render the inverse problem even more ill-posed. We refer to the reviews \cite{MR1955896,MR1669729, MR2986262} and the references therein for more details. A standard cure against the ill-posedness is to regularize the inverse mapping. In this paper the regularization is achieved by considering smooth perturbations of the domains $\mathcal{B}ega^+$.
	
	To solve problem \varepsilonqref{EIT-SO}, we use an optimization approach by minimizing the shape functionals
	\begin{align}
	\label{eit3.1}J_1(\mathcal{B}ega^+) &= \textnormal{fr}ac{1}{2}\int_D (\dot{\mathbf{u}}i(\mathcal{B}ega^+) - u_n(\mathcal{B}ega^+))^2 ,\\
	\label{eit3.2}J_2(\mathcal{B}ega^+) &= \textnormal{fr}ac{1}{2}\int_{\Gamma_n} (u_n(\mathcal{B}ega^+) - h)^2. 

	\varepsilonnd{align}
	Since $\dot{\mathbf{u}}i,u_n\in H^1(D)$ and $h\in H^{1/2}(\Gamma_n)$, $J_1$ and $J_2$ are well-defined. Note that $J_1$ and $J_2$ are redundant for the purpose of the reconstruction but our aim is to provide an efficient way of computing the shape derivative of two functions which are often encountered in the literature.  To compute these derivatives we use the approach described in Section \ref{section2}. First of all introduce

	\begin{align}
	\label{eit4.1}F_1(\varphi_d,\varphi_n) & := \textnormal{fr}ac{1}{2}\int_D (\varphi_d - \varphi_n)^2, \\
	\label{eit4.2}F_2(\varphi_n) & := \textnormal{fr}ac{1}{2}\int_{\Gamma_n} (\varphi_n - h)^2. 
	\varepsilonnd{align}
	Note that $J_1(\mathcal{B}ega^+) = F_1(\dot{\mathbf{u}}i(\mathcal{B}ega^+),u_n(\mathcal{B}ega^+))$ and $J_2(\mathcal{B}ega^+) =F_2(u_n(\mathcal{B}ega^+))$. 
	Next consider $\mathfrak{P}$ a subset of $\mathcal{P}(D)$ and the Lagrangian $\mathcal{L}:\mathfrak{P}\times H^1_d(D)\times H^1_d (D)\times H^1_0 (D)\times H^1_d (D)\rightarrow \mathbf{R}$: 

\begin{align}
\label{G}
\begin{split}		
\mathcal{L}(\mathcal{B}ega^+,\mathbf{{\varphi}},\mathbf{{\psi}}) & := \alpha_1 F_1(\varphi_d,\varphi_n) + \alpha_2 F_2(\varphi_n)\\ 
& + \int_D \sigma \nabla \varphi_d\cdot\nabla \psi_d - f\psi_d  
+ \int_{\Gamma_n} -\sigma^-\partial_n\psi_d (\varphi_d - h)\\ 
&+ \int_D \sigma \nabla \varphi_n\cdot\nabla \psi_n - f\psi_n - \int_{\Gamma_n} g\psi_n,
\varepsilonnd{split}		
\varepsilonnd{align}
	where $\mathbf{{\varphi}} :=(\varphi_d,\varphi_n)$ and $\mathbf{{\psi}} :=(\psi_d,\psi_n)$. 
	The term $-\sigma^-\partial_n\psi_d$ in the second integral of \varepsilonqref{G} is used to enforce the boundary condition \varepsilonqref{eit2.3_}.

	Introduce the objective functional
	\begin{align*}
	J(\mathcal{B}ega^+) := \alpha_1 J_1(\mathcal{B}ega^+) + \alpha_2 J_2(\mathcal{B}ega^+).
	\varepsilonnd{align*}

To compute the derivative of the Lagrangian depending on \varepsilonqref{eit1.1_}-\varepsilonqref{eit2.3_} we apply the averaged adjoint method from Section \ref{section2}. 
	

	
	

	\subsection{State and adjoint equations}

	The state $\mathbf{u}:=(\dot{\mathbf{u}}i,u_n)$ and adjoint state $\mathbf{p}:=(p_d,p_n)$ are solutions of the equations: 
	\begin{align}
		\label{lag1} \partial_{\mathbf{{\psi}}}  \mathcal{L}(\mathcal{B}ega^+,\mathbf{u}, \mathbf{p})(\hat \mathbf{{\psi}}) &= 0\mbox{ for all } \hat \mathbf{{\psi}} \in H^1_0(D)\times H^1_d(D),\\
		\label{lag2} \partial_{\mathbf{{\varphi}}}   \mathcal{L}(\mathcal{B}ega^+, \mathbf{u} ,\mathbf{p})(\hat \mathbf{{\varphi}}) &= 0\mbox{ for all } \hat \mathbf{{\varphi}} \in H^1_d(D)\times H^1_d(D) .
		\varepsilonnd{align}
	
	Writing \varepsilonqref{lag1} explicitely, one can obtain easily the state equations \varepsilonqref{varun} and \varepsilonqref{varud}. 
	Then \varepsilonqref{lag2} yields the equation for the adjoint $p_d$:

	$$\partial_{\varphi_d} \mathcal{L}(\mathcal{B}ega^+,\mathbf{u},\mathbf{p})(\hat{\varphi}_d) = 0,\mbox{ for all } \hat{\varphi}_d \in H^1_d (D),$$
	which leads to 
	\begin{align}\label{eq:lambda}
	&\int_D \sigma \nabla  p_d\cdot \nabla \hat{\varphi}_d \, dx  = - \alpha_1 \int_D (\dot{\mathbf{u}}i - u_n) \hat{\varphi}_d \, dx -\int_{\Gamma_n} -\sigma^-\partial_np_d\hat{\varphi}_d \, ds \quad \text{ for all } \hat{\varphi}_d \in H^1_d (D)
	\varepsilonnd{align}
	which is the variational formulation for the adjoint state $p_d$. This yields the following variational formulation when test functions are restricted to $H^1_0(D)$:
	\begin{align}\label{eit4b.0}
	& \int_D \sigma \nabla  p_d\cdot \nabla \widetilde\varphi \, dx = -\alpha_1 \int_D (\dot{\mathbf{u}}i - u_n) \widetilde\varphi\, dx\quad  \text{ for all } \widetilde\varphi\in  H^1_0(D).
	\varepsilonnd{align}
	If we use Assumption~\ref{assump1}, we get $p_d\in PH^k(D)$ and using Green's formula in $\mathcal{B}ega^+$ and $\mathcal{B}ega^-$ with  $\widetilde\varphi\in C_c^\infty(\mathcal{B}ega^+)$ 
	and $\widetilde\varphi \in C_c^\infty(\mathcal{B}ega^-)$, we obtain the strong form
	\begin{align}
	\label{eit4b.1}-\operatorname{div}(\sigma \nabla p_d) & = -\alpha_1(\dot{\mathbf{u}}i - u_n)\mbox{ in }\mathcal{B}ega^+ \mbox{ and } \mathcal{B}ega^-.
	\varepsilonnd{align}
	Hence using now Green's formula in \varepsilonqref{eq:lambda} and using 
	\varepsilonqref{eit4b.1} gives 
	\begin{align*}
	& \int_{\Gamma^+} [\sigma \partial_n p_d]_{\Gamma^+} \hat{\varphi}_d \, ds  + \int_{\Gamma_n} (\sigma\partial_n p_d  - \sigma^-\partial_np_d) \hat{\varphi}_d\, ds =0 \quad \text{ for all } \hat{\varphi}_d\in H^1_d (D),
	\varepsilonnd{align*}
	where $[\sigma \partial_n p_d]_{\Gamma^+} = \sigma^+ \partial_n p_d^+ - \sigma^- \partial_n p_d^-$ is the jump of $\sigma \partial_n p_d$ across $\Gamma^+$. Since 
	the integral on $\Gamma_n$ above vanishes and $p_d\in H^1_0(D)$, we obtain 
	\begin{align}

	\label{eit4b.3} p_d & = 0 \mbox{ on }\Gamma,\\
	\label{eit4b.4} \sigma^+ \partial_n p_d^+ & = \sigma^- \partial_n p_d^-  \mbox{ on }\Gamma^+.
	\varepsilonnd{align}

	

	Finally solving 
	$$\partial_{\varphi_n} \mathcal{L}(\mathcal{B}ega^+,\mathbf{u},\mathbf{p})(\hat{\varphi}_n) = 0,  \mbox{ for all } \hat{\varphi}_n\in H^1_d (D),$$
	leads to the variational formulation
	\begin{align}\label{eit5}
	& \int_D -\alpha_1(\dot{\mathbf{u}}i - u_n)\hat{\varphi}_n +\sigma \nabla p_n\cdot \nabla \hat{\varphi}_n  + \int_{\Gamma_n}  \alpha_2(u_n - h )\hat{\varphi}_n =0
	\varepsilonnd{align} 
	for all $\hat{\varphi}_n \in H^1_d (D)$. 
	
	
	Similarly as for $p_d$ we get, under Assumption~\ref{assump1}, $p_n\in PH^k(D)$ and the strong form 
	\begin{align}
	\label{eit4.3} -\operatorname{div}(\sigma \nabla p_n) & = \alpha_1(\dot{\mathbf{u}}i - u_n)\mbox{ in }\mathcal{B}ega^+ \mbox{ and } \mathcal{B}ega^-,\\
	\label{eit4.4} \sigma\partial_n p_n & = -\alpha_2(u_n - h )\mbox{ on }\Gamma_n,\\
	\label{eit4.5} p_n & = 0\mbox{ on }\Gamma_d,\\
	\sigma^+\partial_np_n^+ & =\sigma^-\partial_n p_n^- \text{ on } {\Gamma^+}, \quad  p_n^+=p_n^- \text{ on } {\Gamma^+}.
	\varepsilonnd{align}

 	\subsection{Shape derivatives}

	
	Let us consider a transformation $\Phi_t^{\theta}$ defined by \varepsilonqref{Vxt} with $\theta\in C^1_c(D,\mathbf{R}^d)$. Note that $\Phi_t^{\theta}(D) = D$ but in general $\Phi_t^{\theta}(\mathcal{B}ega^+) \neq \mathcal{B}ega^+$. We use the notation $\mathcal{B}ega^+(t):= \Phi_t^{\theta}(\mathcal{B}ega^+)$. Our aim is to show the shape differentiability of $J(\mathcal{B}ega^+)$ with the
	help of Theorem~\ref{thm:sturm}. 
	Following the methodology described in Section \ref{section2} we introduce
	\begin{equation}\label{eq:G_on_the_moved_domain}
	G(t,\mathbf{{\varphi}},\mathbf{{\psi}}):=\mathcal{L}(\mathcal{B}ega^+(t),\mathbf{{\varphi}}\circ\Phi_t^{-1},\mathbf{{\psi}}\circ\Phi_t^{-1}).
	\varepsilonnd{equation}
We proceed to the  change of variables $\Phi_t(x)=y$ in \varepsilonqref{eq:G_on_the_moved_domain}  to get the canonical form \varepsilonqref{G_lag}. First of all let us denote $f_{\mathcal{B}ega^+(t)}=f^+\chi_{\mathcal{B}ega^+(t)}+f^-\chi_{\mathcal{B}ega^-(t)}$ and $\sigma_{\mathcal{B}ega^+(t)}=\sigma^+\chi_{\mathcal{B}ega^+(t)}+\sigma^-\chi_{\mathcal{B}ega^-(t)}$; recall that $\sigma^\pm$ are scalars but $f^{\pm}$ are functions. Then note that   the  change of variables $\Phi_t(x)=y$ leads to considering the following functions inside the integrals:
\begin{align*}
\sigma_{\mathcal{B}ega^+(t)}\circ \Phi_t
& =\sigma^+\chi_{\mathcal{B}ega^+(t)}\circ \Phi_t+\sigma^-\chi_{\mathcal{B}ega^-(t)}\circ \Phi_t = \sigma^+\chi_{\mathcal{B}ega^+}+\sigma^-\chi_{\mathcal{B}ega^-}= \sigma,\\
f_{\mathcal{B}ega^+(t)}\circ \Phi_t
& =f^+\circ \Phi_t\, \chi_{\mathcal{B}ega^+(t)}\circ \Phi_t
+f^-\circ \Phi_t\, \chi_{\mathcal{B}ega^-(t)}\circ \Phi_t
=f^+\circ \Phi_t\, \chi_{\mathcal{B}ega^+}
+f^-\circ \Phi_t\, \chi_{\mathcal{B}ega^-}.
\varepsilonnd{align*}
Thus we introduce the function $\tilde f_t := f^+\circ \Phi_t\, \chi_{\mathcal{B}ega^+} +f^-\circ \Phi_t\, \chi_{\mathcal{B}ega^-}$.
Now we  obtain the canonical form  \varepsilonqref{G_lag} for the Lagrangian:

	\begin{align}
	\label{G_Omega}
	G(t,\mathbf{{\varphi}},\mathbf{{\psi}}) & = a(t,\mathbf{{\varphi}},\mathbf{{\psi}}) + b(t,\mathbf{{\varphi}}),
	\varepsilonnd{align}
	with

	\begin{align*}
	a(t,\mathbf{{\varphi}},\mathbf{{\psi}}) :=& \int_D \sigma A(t)\nabla \varphi_d\cdot\nabla \psi_d
	-\tilde f_t \psi_d \xi(t)  - \int_{\Gamma_n} \sigma^{-1} \partial_n\psi_d(\varphi_d - h) \\
	\notag & + \int_D \sigma A(t)\nabla \varphi_n\cdot\nabla \psi_n 
	-\tilde f_t  \psi_n \xi(t)  - \int_{\Gamma_n} g\psi_n,\\
	b(t,\mathbf{{\varphi}}) :=&  \textnormal{fr}ac{\alpha_1}{2}\int_D (\varphi_d - \varphi_n)^2 \xi(t) + \textnormal{fr}ac{\alpha_2}{2}\int_{\Gamma_n} (\varphi_n - h)^2
	\varepsilonnd{align*}

	where the Jacobian $\xi(t)$ and $A(t)$ are defined as
	$
	\xi(t)  := \det (D\Phi_t)$ and  
	$A(t)  := \xi(t)D\Phi_t^{-1}D\Phi_t^{-T}$.

	In the previous expression \varepsilonqref{G_Omega}, one should note that the integrals on subsets of $\Gamma$ are unchanged since $\Phi_t^{-1} = I$ on $\Gamma$.  Thus we  have $\Phi_t^{\theta}(D) = D$, however the terms inside the integrals on $D$ are modified by the change of variable since $\Phi_t^{-1}\neq I$ inside $D$.
	Note that
	\begin{align*}
	J(\mathcal{B}ega^+(t))=G(t,\mathbf{u}^t,\mathbf{{\psi}}),\text{ for all } \mathbf{{\psi}}\in H^1_0 (D)\times H^1_d (D), 
	\varepsilonnd{align*}
	where $\mathbf{u}^t=(u_n^t,\dot{\mathbf{u}}i^t):=(u_{n,t}\circ\Phi_t,u_{d,t}\circ\Phi_t)$ and $u_{n,t},u_{d,t}$ solve \varepsilonqref{varun},\varepsilonqref{varud}, respectively, with the domains $\mathcal{B}ega^+$ and $\mathcal{B}ega^-$ replaced by $\mathcal{B}ega^+(t)$ and $\mathcal{B}ega^-(t)$. 
	As one can verify by applying a change of variables to \varepsilonqref{varun} and \varepsilonqref{varud} on the domain $\mathcal{B}ega^+(t)$ the functions $u_n^t,\dot{\mathbf{u}}i^t$ satisfy
	\begin{align}\label{varun_t}
	\int_D \sigma A(t)\nabla u_n^t \cdot \nabla \hat{\psi}_n &= \int_D \tilde f_t\hat{\psi}_n + \int_{\Gamma_n} g\hat{\psi}_n \mbox{ for all } \hat{\psi}_n\in H^1_d (D), \\
	\label{varud_t}
	\int_D \sigma A(t) \nabla \dot{\mathbf{u}}i^t \cdot \nabla \hat{\psi}_d &= \int_D \tilde f_t \hat{\psi}_d  \mbox{ for all } \hat{\psi}_d\in H^1_0(D) .
	\varepsilonnd{align}

	Applying standards estimates for elliptic partial differential equations and the fact that $A(t)$ is uniformly bounded from below and above for $t$ small enough, we infer from equations \varepsilonqref{varun_t},\varepsilonqref{varud_t} the existence of constants  $C_1,C_2>0$ independent of $t$ and $\tau>0$ such that for all $t\in [0,\tau]$:
	\begin{equation}\label{apriori_u_D_u_N}
	\|\dot{\mathbf{u}}i^t\|_{H^1(D)}\le C_1,\quad \text{ and } \|u_n^t\|_{H^1(D)}\le C_2.
	\varepsilonnd{equation} 
	From these estimates, we get
	$\dot{\mathbf{u}}i^t\rightharpoonup w_d \text{ and } u_n^t\rightharpoonup w_n \text{ in } H^1(D)\mbox{ as }t\to 0.$
	Passing to the limit in \varepsilonqref{varun_t} and \varepsilonqref{varud_t} yields  $w_d=\dot{\mathbf{u}}i$ and $w_n=u_n$ by uniqueness.
	
	Let us now check Assumption \ref{amp:gateaux_diffbar_G} and the conditions of Theorem~\ref{thm:sturm} for the function $G$ given by \varepsilonqref{G_Omega} and the Banach spaces $E = H^1_d (D)\times H^1_d (D)$ and $F = H^1_0 (D)\times H^1_d (D)$. First of all equation \varepsilonqref{eq:state_G} admits a unique solution $\mathbf{u}^t := (u_n^t,\dot{\mathbf{u}}i^t)$ for each $t\in[0,\tau]$.
	The  conditions of Assumption \ref{amp:gateaux_diffbar_G} are readily satisfied and also the
	function $G$ is affine with respect to $\psi = (\psi_d,\psi_n)$.
	
	Regarding the conditions of Theorem~\ref{thm:sturm}, first note that applying   Lax-Milgram's lemma, we check that both equations \varepsilonqref{eq_bar_p_t_1} and \varepsilonqref{eq_bar_p_t_2} have indeed a unique solution in $F = H^1_0 (D)\times H^1_d (D)$: 
	\begin{align}
	\label{eq_bar_p_t_1}  &\int_D \sigma A(t)\nabla p_d^t\cdot\nabla \hat{\varphi}_d +\textnormal{fr}ac{\alpha_1}{2}\int_D \xi(t) (\dot{\mathbf{u}}i^t+\dot{\mathbf{u}}i - (u_n^t+u_n)) \hat{\varphi}_d
	- \int_{\Gamma_n} \sigma^{-1}\partial_n p_d^t\hat{\varphi}_d =0,\\ 

	\label{eq_bar_p_t_2} &\int_D \sigma A(t)\nabla p_n^t\cdot\nabla \hat{\varphi}_n-\textnormal{fr}ac{\alpha_1}{2}\int_D \xi(t)(\dot{\mathbf{u}}i^t+\dot{\mathbf{u}}i - (u_n^t+u_n)) \hat{\varphi}_n   + \textnormal{fr}ac{\alpha_2}{2}\int_{\Gamma_n} (u_n^t+u_n - 2h)\hat{\varphi}_n =0,
	\varepsilonnd{align}
	for all $\hat{\varphi}_d$, $\hat{\varphi}_n$ in $H^1_d (D)$. 
	Therefore there exists a unique solution $\mathbf{p}^t= (p_n^t,p_d^t)$ of the averaged adjoint equation \varepsilonqref{averated_}.
	
	Now we check Assumption (H1). Testing  \varepsilonqref{eq_bar_p_t_1} with 
	$\hat{\varphi}_d=p_d^t$ and \varepsilonqref{eq_bar_p_t_2} with $\hat{\varphi}_n=p_n^t$, we conclude by an application of H\"older's inequality together with \varepsilonqref{apriori_u_D_u_N} the existence of constants $C_1,C_2$ and $\tau>0$ such that for all $t\in [0,\tau]$
	\begin{equation}n
	\|p_d^t\|_{H^1(D)}\le C_1,\quad \text{ and } \|p_n^t\|_{H^1(D)}\le C_2.
	\varepsilonnd{equation}n 
		We get that for each sequence $t_k$ converging to zero, there exists a subsequence also denoted $t_k$ such that 
		$p_d^{t_k}\rightharpoonup q_d$ and $p_n^{t_k}\rightharpoonup q_n$ for two elements $q_d,q_n\in H^1(D)$.
		Passing to the limit in 
		\varepsilonqref{eq_bar_p_t_1} and \varepsilonqref{eq_bar_p_t_2} yields  $q_d=p_d$ and $q_n=p_n$ by uniqueness, where $p_d$ and $p_n$ are solutions of the adjoint equations. Since the limit is unique, we have in fact $p_d^{t}\rightharpoonup p_d$ and $p_n^{t}\rightharpoonup p_n$ as $t\to 0$.
	Finally, differentiating $G$ with respect to $t$ yields
	\begin{align*}
	\partial_tG(t,\mathbf{{\varphi}},\mathbf{{\psi}}) & = \textnormal{fr}ac{\alpha_1}{2}\int_D (\varphi_d - \varphi_n)^2  \xi(t)\operatorname{tr}(D\theta_t D \Phi_t^{-1})\\
	&\hspace{-1.5cm} + \int_D \sigma A'(t)\nabla \varphi_d\cdot\nabla \psi_d -\tilde f_t \psi_d \xi(t)\operatorname{tr}(D\theta_t D \Phi_t^{-1}) -\psi_d  \widetilde{\nabla} f_t\cdot \theta_t  \xi(t)\\
	&\hspace{-1.5cm} + \int_D \sigma A'(t)\nabla \varphi_n\cdot\nabla \psi_n - \tilde f_t\psi_n \xi(t)\operatorname{tr}(D\theta_t D \Phi_t^{-1}) -\psi_n \widetilde{\nabla} f_t\cdot \theta_t \xi(t).
	\varepsilonnd{align*} 
	where 
	$$\widetilde{\nabla} f_t := \nabla f^+\circ \Phi_t \, \chi_{\mathcal{B}ega^+} + \nabla f^-\circ \Phi_t\, \chi_{\mathcal{B}ega^-},\qquad \theta_t = \theta\circ\Phi_t$$
	$$A'(t) = \operatorname{tr}(D\theta^t D\Phi_t^{-1}) A(t) - D\Phi_t^{-T} D\theta_t A(t) -(D\Phi_t^{-T} D\theta_t A(t) )^T$$ 
	and $D\theta_t$ is the Jacobian matrix of $\theta_t$.
	In view of $\theta\in C^1_c(D,\mathbf{R}^d)$, the functions $t\mapsto D\theta_t$ and $t\mapsto \operatorname{tr}( D\theta_t \Phi_t^{-1})= \operatorname{div}(\theta)\circ \Phi_t$ are continuous on $[0,T]$. Moreover $\varphi_d,\varphi_n,\psi_d,\psi_n$ are in $H^1(D)$, $f\in PH^1(D)$ so that $\partial_tG(t,\mathbf{{\varphi}},\mathbf{{\psi}}) $ is well-defined for all $t\in [0,T]$. 
		Using the weak convergence of $\mathbf{p}^t$ and the strong differentiability of 
		$t\mapsto A(t)$ and $t\mapsto \xi(t)$  it follows 
		\begin{equation}
		\lim_{t \searrow 0 }  \textnormal{fr}ac{G( t,\mathbf{u}^0,\mathbf{p}^t ) - G(0,\mathbf{u}^0, \mathbf{p}^t) }{t}=\partial_t G(0,\mathbf{u}^0,\mathbf{p}^0).
		\varepsilonnd{equation}
	
	Thus we have verified all assumptions from Theorem \ref{thm:sturm}. This yields
	\begin{equation}n
	dJ(\mathcal{B}ega^+)(\theta)=\dt \left(G(t,\mathbf{u}^t,\mathbf{{\psi}})\right)|_{t=0}=\partial_t G(0,\mathbf{u}^0,\mathbf{p}^0)\text{ for all } \mathbf{{\psi}}\in F= H^1_0 (D)\times H^1_d (D) , 
	\varepsilonnd{equation}n 
	and therefore we have proved the following result.

	\begin{proposition}[distributed shape derivative]\label{domexpr} 
		Let $D\subset\mathbf{R}^d$ be a Lipschitz domain,  $\theta\in C^1_c(D,\mathbf{R}^d)$, $f~\in~PH^1(D)$, $g\in H^{-1/2}(\Gamma_n)$, $h\in H^{1/2}(\Gamma_n)$, $\mathcal{B}ega^+\subset D$ is an open set, then the  shape derivative of $J(\mathcal{B}ega^+)$ is given by 

		\begin{align}\label{volexpr}
		\begin{split}
		dJ(\mathcal{B}ega^+)(\theta) & = \int_D \left(\textnormal{fr}ac{\alpha_1}{2}(\dot{\mathbf{u}}i - u_n)^2 -f(p_n+p_d)\right)\operatorname{div} \theta \\
		&\hspace{-1cm}+ \int_D  - (p_d+p_n)\widetilde{\nabla} f\cdot \theta + \sigma A'(0)(\nabla \dot{\mathbf{u}}i\cdot\nabla p_d  + \nabla u_n\cdot\nabla p_n)  ,
		\varepsilonnd{split}
		\varepsilonnd{align} 

		where $\widetilde{\nabla} f := \nabla f^+ \, \chi_{\mathcal{B}ega^+} + \nabla f^-\, \chi_{\mathcal{B}ega^-}$, $A'(0) = (\operatorname{div} \theta) I - D\theta^T - D\theta$,  $u_n,\dot{\mathbf{u}}i$ are solutions of  \varepsilonqref{varun},\varepsilonqref{varud} and  $p_n,p_d$ of \varepsilonqref{eit5}, \varepsilonqref{eq:lambda}. 
		
		The shape derivative \varepsilonqref{volexpr} also has the tensor representation corresponding to \varepsilonqref{ea:volume_from}:
		\begin{align}\label{volexpr2}
		\begin{split}
		dJ(\mathcal{B}ega^+)(\theta) & = \int_D \mathbf{S}_1 \cdot D\theta + \mathbf{S}_0\cdot \theta  ,
		\varepsilonnd{split}
		\varepsilonnd{align} 
		where 
		\begin{align*}
		\mathbf{S}_1 & = -\sigma (\nabla \dot{\mathbf{u}}i \otimes \nabla p_d + \nabla p_d \otimes \nabla \dot{\mathbf{u}}i + \nabla u_n \otimes \nabla p_n + \nabla p_n \otimes \nabla u_n)   \\
		&\quad  +\sigma(\nabla \dot{\mathbf{u}}i\cdot\nabla p_d  + \nabla u_n\cdot\nabla p_n)I + \left(\textnormal{fr}ac{\alpha_1}{2}(\dot{\mathbf{u}}i - u_n)^2 -f(p_n+p_d)\right)I,\\
		\mathbf{S}_0 & = - (p_d+p_n)\widetilde{\nabla} f.

		\varepsilonnd{align*}
	\varepsilonnd{proposition}

	Note that the volume expression of the shape gradient in  Proposition \ref{domexpr}  has been obtained without any regularity assumption on $\mathcal{B}ega^+$. In order to obtain a boundary expression on the interface $\Gamma^+$  we need more regularity of  $\mathcal{B}ega^+$ provided by Assumption \ref{assump1}. If it is satisfied,  we can apply Corollary \ref{corollary2} to obtain directly the boundary expression of the shape derivative, using mainly the standard tensor relation $(\nabla \dot{\mathbf{u}}i \otimes \nabla p_d) n = (\nabla p_d\cdot n)\nabla \dot{\mathbf{u}}i$, which yields Proposition \ref{prop2}.

	\begin{proposition}[boundary expression]\label{prop2}
		Under Assumption \ref{assump1} and $\theta\in C^1_c(D,\mathbf{R}^d)$ the shape derivative of $J(\mathcal{B}ega^+)$ is given by
   \begin{align*}
   dJ(\mathcal{B}ega^+)(\theta)  =&  \int_{\Gamma^+} \left[ \sigma(-\partial_n \dot{\mathbf{u}}i\partial_n p_d - \partial_n u_n\partial_n p_n )\right]_{\Gamma^+} \theta\cdot n
   \\
   & + \int_{\Gamma^+}  ( [\sigma]_{\Gamma^+}(\nabla_{\Gamma^+} \dot{\mathbf{u}}i\cdot\nabla_{\Gamma^+} p_d + \nabla_{\Gamma^+} u_n\cdot \nabla_{\Gamma^+} p_n ) -[f]_{\Gamma^+}(p_n+p_d)) \theta\cdot n.
   \varepsilonnd{align*}
	\varepsilonnd{proposition}
	Note that our results cover and generalize several results that can be found in the literature of shape optimization approaches for EIT, including \cite{MR2329288,MR2536481}. For instance taking $\alpha_2 = 1$, $\alpha_1 = 0$ in Proposition \ref{prop2} we get $p_D \varepsilonquiv 0$ which yields the same formula
 as the one obtained in \cite[pp. 533]{MR2329288}. 
	
	Note also that from a numerical point of view, the boundary expression in Proposition \ref{prop2} is delicate to compute compared to the domain expression in Proposition \ref{domexpr} for which the gradients of the state and adjoint states can be straightforwardly computed at grid points when using the finite element method for instance. The boundary expression, on the other hand, needs here the computation of the normal vector and the interpolation of the gradients on the interface $\Gamma^+$ which requires a precise description of the boundary and introduces an additional error.

	
	

	\section{Level set method}\label{section5}

	The level set method, originally introduced in \cite{MR965860}, gives a general framework for the computation of evolving interfaces using an implicit representation of these interfaces.  The core idea of this method is to represent the boundary of the moving domain $\mathcal{B}ega^+(t)\subset D \in \mathbf{R}^N$ as the level set of a continuous function $\phi(\cdot,t): D\to\mathbf{R}$.
	
	Let us consider the family of domains $\mathcal{B}ega^+(t)\subset D$ as defined in \varepsilonqref{domain}. Each domain $\mathcal{B}ega^+(t)$ can be defined as
	\begin{equation}
	\mathcal{B}ega^+(t) :=\{x\in D,\ \phi(x,t) < 0\}
	\varepsilonnd{equation}
	where $\phi: D\times \mathbf{R}^+ \to \mathbf{R}$ is continuous and called {\it level set function}. Indeed, if we assume $|\nabla\phi(\cdot,t)|\neq 0$ on the set  $\{x\in D,\ \phi(x,t) = 0\}$ then we have
	\begin{equation}
	\partial \mathcal{B}ega^+(t) = \{x\in D,\ \phi(x,t) = 0\},
	\varepsilonnd{equation}
	i.e. the boundary $\partial \mathcal{B}ega^+(t)$ is the zero level set of  $\phi(\cdot,t)$.
	
	
	Let $x(t)$ be the position of a particle on the boundary
	$\partial \mathcal{B}ega^+(t)$ moving with velocity $\dot{x}(t)=\theta(x(t))$ according to \varepsilonqref{Vxt}.
	Differentiating the relation $\phi(x(t),t)=0$ with respect to $t$
	yields the Hamilton-Jacobi equation:
	\begin{equation*}
	\partial_t\phi (x(t),t)+ \theta(x(t))\cdot \nabla \phi(x(t),t) = 0  \quad \mbox{ in } \partial \mathcal{B}ega^+(t)\times\mathbf{R}^+,
	\varepsilonnd{equation*}
	which is then extended to all of $D$ via the equation
	\begin{equation}
	\partial_t\phi(x,t) + \theta(x)\cdot \nabla \phi(x,t) = 0  \quad \mbox{ in } D\times\mathbf{R}^+,
	\label{eq:transport}
	\varepsilonnd{equation}
	or alternatively to $U\times\mathbf{R}^+$ where $U$ is a neighbourhood of $\partial \mathcal{B}ega^+(t)$.
	
	Traditionally, the level set method has been designed to track smooth interfaces moving along the normal direction to the boundary. Theoretically, this is supported by Theorem \ref{thm:structure_theorem}, i.e. if the domain $\mathcal{B}ega^+(t)$ and the shape gradient are smooth enough then the shape derivative only depends on $\theta\cdot n$ on $\partial\mathcal{B}ega^+(t)$. In this case, we may choose for the optimization a vector field $\theta = \vartheta_n n$ on $\partial\mathcal{B}ega^+(t)$. Then, noting that an extension to $D$ of the unit outward normal vector $n$ to $\mathcal{B}ega^+(t)$ is given 
	by $n =\nabla \phi /|\nabla \phi|$, and extending $\vartheta_n$ to all of $D$, one obtains from \varepsilonqref{eq:transport} the level set equation 
	\begin{equation}
	\partial_t \phi + \vartheta_n |\nabla \phi| = 0 \quad \mbox{ in } D\times\mathbf{R}^+.
	\label{eq:HJ1}
	\varepsilonnd{equation}
	The initial data $\phi(x,0)=\phi_0(x)$ accompanying the Hamilton-Jacobi equation \varepsilonqref{eq:transport} or \varepsilonqref{eq:HJ1} is chosen as the signed distance function to the
	initial boundary $\partial \mathcal{B}ega^+(0)$ in order to satisfy the condition $|\nabla u|\neq 0$ on $\partial\mathcal{B}ega^+$, i.e.
	\begin{equation}
	\phi_0(x) = 
	\left\{\begin{array}{rl}
	d(x,\partial \mathcal{B}ega^+(0)), & \mbox{ if } x\in (\mathcal{B}ega^+(0))^c, \\
	-d(x,\partial \mathcal{B}ega^+(0)), & \mbox{ if } x\in \mathcal{B}ega^+(0) .
	\varepsilonnd{array}
	\right.
	\varepsilonnd{equation}
	
	

	\subsection{Level set method and domain expression}

	
	In the case of the distributed shape derivative, for instance  \varepsilonqref{volexpr} or \varepsilonqref{volexpr2}, $\phi$ is not governed by \varepsilonqref{eq:HJ1} but rather by the Hamilton-Jacobi equation \varepsilonqref{eq:transport}. Indeed we obtain a descent direction $\theta$ defined in $D$ by solving \varepsilonqref{VP_1}, where $dJ(\mathcal{B}ega^+)$ is given by Proposition \ref{domexpr} which can subsequently be used in \varepsilonqref{eq:transport} to compute the evolution of $\phi$. On the other hand, in the usual level set method, one solves a PDE {\it on the  boundary} $\partial\mathcal{B}ega^+$ in an analogous way as for   \varepsilonqref{VP_1} (for instance using a Laplace-Beltrami operator), and uses the boundary expression from Proposition \ref{prop2} to obtain $\vartheta_n = \theta\cdot n$ on $\partial\mathcal{B}ega^+$. 
	
	Numerically it is actually more straightforward in many cases to use \varepsilonqref{eq:transport} instead of   \varepsilonqref{eq:HJ1}. Indeed, when using \varepsilonqref{eq:HJ1}, $\vartheta_n$ is initially only given on $\partial \mathcal{B}ega^+(t)$ and must be extended to the entire domain $D$ or at least to a narrow band around $\partial \mathcal{B}ega^+(t)$. 
	Therefore it is convenient to use \varepsilonqref{eq:transport} with  $\theta$  already defined in $D$ as is the case of the distributed shape derivative, which provides an extension to $D$ or to a narrow band around $\partial \mathcal{B}ega^+(t)$.
	
	In shape optimization,  $\vartheta_n$ usually depends on the solution of one or several PDEs and their gradient. Since the boundary $\partial\mathcal{B}ega^+(t)$  in general does not match the grid nodes where $\phi$ and the solutions of the  partial differential equations are defined in the numerical application, the computation of $\vartheta_n$ 
	requires the  interpolation on $\partial\mathcal{B}ega^+(t)$ of functions defined at the grid points only, complicating the numerical implementation  and introducing an 
	additional interpolation error. This is an issue in particular for interface problems where $\vartheta_n$ is the jump of a function across the interface, as in Proposition \ref{prop2}, which requires multiple interpolations and is error-prone. In the distributed shape derivative framework  $\theta$ only needs to be defined at grid nodes.  
	
	
	

	\subsection{Discretization of the level set equation}

	Let $D$ be the unit square $D=(0,1)\times(0,1)$ to fix ideas.
	For the discretization of the Hamilton-Jacobi equation \varepsilonqref{eq:transport},
	we first define the mesh grid corresponding to $D$. We introduce the nodes $P_{ij}$
	whose coordinates are given by $(i\Delta x, j \Delta y)$, $1\leq i,j\leq N$ where $\Delta x$ 
	and $\Delta y$ are the steps discretization in the $x$ and $y$
	directions respectively. Let us also write $t^k = k\Delta t$
	the discrete time for $k\in \mathbb{N}$, where $\Delta t$ is the time
	step. We are seeking for an approximation
	$\phi_{ij}^k \simeq \phi(P_{ij},t^k)$.
	
	In the usual level set method, the level set equation \varepsilonqref{eq:HJ1} is discretized using an explicit upwind scheme proposed by Osher and Sethian \cite{MR1939127,MR965860,MR1700751}. This scheme applies to the specific form \varepsilonqref{eq:HJ1} but is not suited to discretize \varepsilonqref{eq:transport} required for our application. Equation \varepsilonqref{eq:transport} is of the form
	\begin{equation}
	\partial_t\phi + H(\nabla \phi)= 0  \quad \mbox{ in } D\times\mathbf{R}^+.
	\varepsilonnd{equation}
	where $ H(\nabla \phi) := \theta\cdot \nabla \phi $ is the so-called Hamiltonian. We use the Local Lax-Friedrichs flux originally conceived in \cite{MR1111446} and which reduces in our case to:
	$$\hat H^{LLF}(p^-,p^+,q^-,q^+)  = H\left(\textnormal{fr}ac{p^- + p^+}{2},\textnormal{fr}ac{q^- +q^+}{2}\right) -\textnormal{fr}ac{1}{2} (p^+-p^-)\alpha^x -\textnormal{fr}ac{1}{2}(p^+-p^-)\alpha^y  $$
	where $\alpha^x = |\theta_x|$, $\alpha^y = |\theta_y|$, $\theta = (\theta_x,\theta_y)$ and
	\begin{align*}
	p^- &= D_x^{-}\phi_{ij}=\dfrac{\phi_{ij}-\phi_{i-1,j}}{\Delta x}, &
	p^+  = D_x^{+}\phi_{ij}=\dfrac{\phi_{i+1,j}-\phi_{ij}}{\Delta x},\\
	q^- &= D_y^{-}\phi_{ij}=\dfrac{\phi_{ij}-\phi_{i,j-1}}{\Delta y}, &
	q^+  = D_y^{+}\phi_{ij}=\dfrac{\phi_{i,j+1}-\phi_{ij}}{\Delta y}
	\varepsilonnd{align*}
	are the backward and forward approximations of the $x$-derivative and $y$-derivative of 
	$\phi$ at $P_{ij}$, respectively. Using a forward Euler time discretization, the numerical scheme corresponding to \varepsilonqref{eq:transport} is 
	\begin{equation}
	\phi_{ij}^{k+1}=\phi_{ij}^k - \Delta t \ \hat H^{LLF}(p^-,p^+,q^-,q^+)
	\varepsilonnd{equation}

	

	
	
	For numerical accuracy, the solution of the level set equation 
	\varepsilonqref{eq:transport} should not be too flat or too steep. 
	This is fulfilled for instance if $\phi$ is the distance function
	i.e. $|\nabla \phi| = 1$. Even if one initializes $\phi$ using a
	signed distance function, the solution
	$\phi$ of the level set equation \varepsilonqref{eq:transport} does not generally
	remain close to a distance function. We may occasionally perform a reinitialization
	of $\phi$ by solving a parabolic equation up to the
	stationary state; see \cite{MR2356899,MR2459656,MR1723321}.
	Although in the level set method this reinitialization is standard, in the case of the distributed shape gradient, we observe experimentally that the level set function $\phi$ stays close to a distance function during the iterations and we do not need to reinitialize. 
		The regularization of the shape gradient could explain this observed stability of the level set function.
	
	The computational efficiency of the level set method can be improved by using the so-called ``narrow band'' approach introduced in \cite{MR1329634}, which consists in computing and updating the level set function only on a thin region around the interface. This allows to reduce the complexity of the problem to $N\log(N)$ instead of $N^2$ in two dimensions.
	In this paper we do not implement this approach but we mention that it could also be applied to the distributed shape derivative approach and equation \varepsilonqref{eq:transport} by taking $\theta$ with a support in a narrow band around the moving interface, which can be achieved by choosing the appropriate space $\textbf{E}$ in \varepsilonqref{eq:bilinear}. 
	
	
	
	

	\section{Application and numerical results}\label{section6}

	
	

	\subsection{Electrical impedance tomography}
	
	In this section we give numerical results for the problem of electrical impedance tomography presented in Section \ref{sec:eit}, precisely we look for an approximate solution of the shape optimization problem \varepsilonqref{EIT-SO}.
	Using the notations of Section  \ref{sec:eit} we take  $D = (0,1)\times (0,1)$ and $\Gamma_d = \varepsilonmptyset$, i.e. we have measurements on the entire boundary $\Gamma$.  For easiness of implementation, we consider a slightly different problem than the one in Section \ref{sec:eit}. Denote $\Gamma_{t}$, $\Gamma_{b}$, $\Gamma_{l}$ and $\Gamma_{r}$ the four sides of the square, where the indices $t,b,l,r$ stands for top, bottom, left and right, respectively. We consider the following problems: find $u_n\in H^1_{tb}(D)$
	\begin{equation}\label{varun-b}
	\int_D \sigma \nabla u_n \cdot \nabla \varphi = \int_D f\varphi + \int_{\Gamma_{l}\cup\Gamma_{r}} g\varphi\ \mbox{ for all }\ \varphi\in H^1_{0,tb}(D) 
	\varepsilonnd{equation}
	and find $\dot{\mathbf{u}}i\in H^1_{lr}(D)$ such that 
	\begin{equation}\label{varud-b}
	\int_D \sigma \nabla \dot{\mathbf{u}}i \cdot \nabla \varphi = \int_D f\varphi + \int_{\Gamma_{t}\cup\Gamma_{b}} g\varphi\ \mbox{ for all }\ \varphi\in H^1_{0,lr}(D)  
	\varepsilonnd{equation}
	where
	\begin{align*}
	H^1_{tb}(D)&:=\{v\in H^1(D)\ |\ \; v=h\mbox{ on }\Gamma_t \cup \Gamma_b\}, \\
	H^1_{lr}(D)&:=\{v\in H^1(D)\ |\ \; v=h\mbox{ on }\Gamma_l \cup \Gamma_r\}, \\
	H^1_{0,tb}(D)&:=\{v\in H^1(D)\ |\ \; v=0\mbox{ on }\Gamma_t \cup \Gamma_b\}, \\
	H^1_{0,lr}(D)&:=\{v\in H^1(D)\ |\ \; v=0\mbox{ on }\Gamma_l \cup \Gamma_r\}.
	\varepsilonnd{align*}
	In our experiments we choose $f\varepsilonquiv 0$. The results of Section \ref{sec:eit} can be straightforwardly adapted to equations \varepsilonqref{varun-b}, \varepsilonqref{varud-b}.

	
	We use the software package FEniCS for the implementation; see \cite{fenics:book}. The domain $D$ is meshed using a regular grid of $128\times 128$ elements and we describe the evolution of the interface $\Gamma^+$ using the level set method from Section \ref{section5}. The conductivity values are set to $\sigma_0 = 1$ and $\sigma_1 = 10$. 
	
	We obtain measurements $h_i$ corresponding to fluxes $g_i$, $i=1,..,I$, by taking the trace on $\Gamma$ of the solution of a Neumann problem where the fluxes are equal to $g_i$. To simulate real noisy EIT data, the measurements $h_i$ are corrupted by adding a normal Gaussian noise with mean zero and standard deviation $\delta*\|h_i\|_\infty$, where $\delta$ is a parameter. The noise level is computed as
	\begin{equation}\label{noise_compute}
	noise = \textnormal{fr}ac{\sum_{i=1}^I\|h_i - \tilde h_i\|_{L^2(\Gamma)}}{\sum_{i=1}^I\|h_i\|_{L^2(\Gamma)}} 
	\varepsilonnd{equation}
	where $\tilde h_i$ is the noisy measurement and $h_i$ the synthetic measurement without noise on $\Gamma$.
	
	
	We use a variation of the functional \varepsilonqref{eit3.1}, i.e. in our context:
	\begin{align}\label{funcJ_num}
	J(\mathcal{B}ega^+) &=  \sum_{i=1}^I \mu_i  \int_\mathcal{B}ega   \textnormal{fr}ac{1}{2} (u_{d,i}(\mathcal{B}ega^+) - u_{n,i}(\mathcal{B}ega^+))^2 ,
	\varepsilonnd{align}
	where $u_{d,i}$ and $u_{n,i}$ correspond to the different fluxes $g_i$.

	Here the coefficients $\mu_i$ are weights associated to the fluxes $g_i$. In our experiments we choose the weights $\mu_i$ such that each term of the sums in \varepsilonqref{funcJ_num} are equal to $1$ on initialization in order to have a well-distributed influence of each term. Practically, the $\mu_i$ are thus calculated during the first iteration.

	We use the distributed shape derivative $dJ(\mathcal{B}ega^+)$ from Proposition \varepsilonqref{domexpr}. 
	We obtain a descent direction by solving \varepsilonqref{VP_1} with $\textbf{E}$ a finite dimensional subspace of $H^1_0(D)$ and
		$\mathcal B(v,w) = \int_D D v \cdot D w.$

		We choose $\textbf{E}$ to be  the space of linear Lagrange elements. 
	
	
	
	Since we use a gradient-based method we implement an Armijo line search to adjust the time-stepping. The algorithm is stopped when the decrease of the functional becomes insignificant, practically when the following stopping criterion is repeatedly satisfied: 
	$$J(\mathcal{B}ega^+_{k}) - J(\mathcal{B}ega^+_{k+1}) < \gamma (J(\mathcal{B}ega^+_{0}) - J(\mathcal{B}ega^+_{1})) $$  
	where $\mathcal{B}ega^+_{k}$ denotes the $k$-th iterate of $\mathcal{B}ega^+$. We take $\gamma = 5.10^{-5}$ in our tests. 
	
	

	
	\begin{figure}[ht]
		\begin{center}
			\includegraphics[width=0.32\textwidth]{./initialization01}
			\includegraphics[width=0.32\textwidth]{./reconstruction00}
			\includegraphics[width=0.32\textwidth]{./reconstruction003}\\
			\includegraphics[width=0.32\textwidth]{./reconstruction01}
			\includegraphics[width=0.32\textwidth]{./reconstruction02}
			\includegraphics[width=0.32\textwidth]{./reconstruction05}
		\varepsilonnd{center}
		\caption{Reconstruction (continuous contours) of two ellipses (dashed contours) with different noise levels and using three measurements. From left to right and top to bottom: initialization (continuous contours - top left), $0\%$ noise (367 iterations), $0.43\%$ noise (338 iterations), $1.44\%$ noise (334 iterations), $2.83\%$ noise (310 iterations), $7\%$ noise (356 iterations).}\label{rec_twoellipses_noise}
	\varepsilonnd{figure}
	

	In Figure \ref{rec_twoellipses_noise} we compare the reconstruction for different noise levels computed using \varepsilonqref{noise_compute}. We take in this example $I=3$, i.e. we use three fluxes $g_i$, $i=1,2,3$, defined as follows:
	\begin{align*}
	g_1 &= 1\mbox{ on }\Gamma_l\cup\Gamma_r \mbox{ and } g_1 = -1\mbox{ on }\Gamma_t\cup\Gamma_b, \\ 
	g_2 &= 1\mbox{ on }\Gamma_l\cup\Gamma_t \mbox{ and } g_2 = -1\mbox{ on }\Gamma_r\cup\Gamma_b,\\
	g_3 &= 1\mbox{ on }\Gamma_l\cup\Gamma_b \mbox{ and } g_3 = -1\mbox{ on }\Gamma_r\cup\Gamma_t. 
	\varepsilonnd{align*}
	Without noise, the reconstruction is very close to the true object and  degrades as the measurements become increasingly noisy,   as is usually the case in EIT. However, the reconstruction is quite robust with respect to noise considering that the problem is severely ill-posed. We reconstruct two ellipses and initialize with two balls placed at the wrong location. The average number of iterations until convergence is around $340$ iterations.  
	
	In Figure \ref{rec_three_inclusions} we reconstruct three inclusions this time using $I=7$ different measurements, with $1.55\%$ noise. The reconstruction is close to the true inclusion and is a bit degraded due to the noise. Figure \ref{rec_three_inclusions2} shows the convergence history of the cost functional in log scale for this example. 
	
	Our algorithm gives good results in comparison to existing results in the literature using level set methods to solve the EIT problem.
	In \cite{MR2132313} the EIT problem has been treated numerically using a level set method, which is not based on the use of shape derivatives but on the differentiation of a smooth approximation of the Heaviside function to represent domains. 
	In \cite{MR2536481} the level set method using the boundary expression of the shape derivative is used based on equation \varepsilonqref{eq:HJ1}.
	
	Our algorithm converges fast in comparison to \cite{MR2132313,MR2536481}:
	convergence occurs after around 300 iterations. In \cite{MR2132313} convergence occurs between 200 iterations for one inclusion and up to 50000 iterations for two inclusions. 
	In \cite{MR2536481} convergence occurs after 2000 or 10000 iterations on two examples with three inclusions. 
	Concerning measurements we obtain good reconstruction of two inclusions with $I=3$ and three inclusions with $I=7$, while in \cite{MR2132313} sets of $4,12,28$ and $60$ measurements are used but usually $60$ measurements are required for complicated shapes such as two inclusions. In \cite{MR2536481}, $60$ measurements are used.
	Nevertheless, our results are not directly comparable since the conductivities are unknown in \cite{MR2132313, MR2536481}, which makes the inverse problem harder and might explain the slower convergence.
	Also, the reconstructed shapes are not the same  although the complexity of the unknown shapes is comparable since we also consider two and three inclusions as in \cite{MR2132313, MR2536481}.
	Only an exact comparison using the same problem, test case, initialization, noise level, number and type of measurements could allow to conclude.
	
	
	\begin{figure}[ht]
		\begin{center}
			\includegraphics[width=0.4\textwidth]{./initialization_three}
			\includegraphics[width=0.4\textwidth]{./reconstruction_three}
		\varepsilonnd{center}
		\caption{Initialization (continuous contours - left) and reconstruction (continuous contours - right) of two ellipses and a ball (dashed contours) with $1.55\%$ noise (371 iterations) and using seven measurements.}\label{rec_three_inclusions}
	\varepsilonnd{figure}
	
	\begin{figure}[ht]
		\begin{center}
			\includegraphics[width=0.5\textwidth]{./cost_history}
		\varepsilonnd{center}
		\caption{History of cost functional corresponding to Figure \ref{rec_three_inclusions} in logarithmic scale.}\label{rec_three_inclusions2}
	\varepsilonnd{figure}
	

	
	{\bf Acknowledgments.} 
	The authors would like to thank the reviewers for their
	helpful comments. Antoine Laurain acknowledges support from the DFG research center MATHEON (MATHEON - Project C37, Shape/Topology optimization methods for inverse problems). Kevin Sturm acknowledges support from the  DFG research center MATHEON, Project C11.
	

	
	

	

\varepsilonnd{document}