# Automated computation of topological derivatives with application to
nonlinear elasticity and reaction-diffusion problems

Peter Gangl  TU Graz, Steyrergasse 30/III, 8010 Graz, Austria
gangl(at)math.tugraz.at and  Kevin Sturm  TU Wien, Wiedner Hauptstr. 8-10,
1040 Vienna, Austria  kevin.sturm(at)tuwien.ac.at

###### Abstract.

Purpose While topological derivatives have proven useful in applications of
topology optimisation and inverse problems, their mathematically rigorous
derivation remains an ongoing research topic, in particular in the context of
nonlinear partial differential equation (PDE) constraints.
Design/methodology/approach We present a systematic yet formal approach for
the computation of topological derivatives of a large class of PDE-constrained
topology optimization problems with respect to arbitrary inclusion shapes.
Scalar and vector-valued as well as linear and nonlinear elliptic PDE
constraints are considered in two and three space dimensions including a
nonlinear elasticity model and nonlinear reaction-diffusion problems. The
systematic procedure follows a Lagrangian approach for computing topological
derivatives.
Findings For problems where the exact formula is known, the numerically
computed values show good coincidence. Moreover, by inserting the computed
values into the topological asymptotic expansion, we verify that the obtained
values satisfy the expected behaviour also for other, previously unknown
problems, indicating the correctness of the procedure.
Originality/value We present a systematic approach for the computation of
topological derivatives that is applicable to a large class of problems. Most
notably, our approach covers the topological derivative for a nonlinear
elasticity problem, which has not been reported in the literature.

###### Key words and phrases:

topological derivative, automated differentiation, nonlinear elasticity,
nonlinear diffusion-convection-reaction equation

##  1\. Introduction

The topological derivative concept was first used for finding optimal
locations of holes in mechanical structures in [11] and was later introduced
in a mathematically concise way in the publications [22] and [17]. Given a
shape function $\mathcal{J}$ that maps a shape $\Omega$ to a real number
$\mathcal{J}(\Omega)$, the topological derivative at a spatial point
$z\in\Omega$ measures the sensitivity of $\mathcal{J}$ with respect to a small
topological perturbation of the shape $\Omega$. Denoting the perturbed shape
by $\Omega_{\varepsilon}$, e.g.,
$\Omega_{\varepsilon}:=\Omega\setminus\omega_{\varepsilon}$ with
$\omega_{\varepsilon}=B_{\varepsilon}(z)$, the topological derivative is
defined as

$\displaystyle d\mathcal{J}(\Omega)(z):=\underset{{\varepsilon}\searrow
0}{\mbox{lim
}}\frac{\mathcal{J}(\Omega_{\varepsilon})-\mathcal{J}(\Omega)}{|\omega_{\varepsilon}|},$
(1.1)

thus satisfying a topologically asymptotic expansion of the form

$\displaystyle\mathcal{J}(\Omega_{\varepsilon})=\mathcal{J}(\Omega)+|\omega_{\varepsilon}|d\mathcal{J}(\Omega)(z)+o({\varepsilon})\quad\mbox{as
}{\varepsilon}\searrow 0.$ (1.2)

Since its introduction, the topological derivative concept has found
application mostly in the context of topology optimization for engineering
applications by level set approaches [4, 1, 9], but has also been utilized in
medical applications such as electrical impedance tomography (EIT) [19] or
mathematical image processing [18].

Most practically relevant engineering applications involve a partial
differential equation (PDE) constraint and thus are of the type

$\displaystyle\underset{\Omega}{\mbox{inf }}J(\Omega,u)\qquad\mbox{subject to
}\;e(\Omega;u)=0$ (1.3)

with a PDE operator $e(\Omega;\cdot)$. In many cases (in particular when the
principal part of the PDE operator is subject to a topological perturbation),
the derivation of topological derivatives involves an asymptotic analysis of
the state variable of the form

$\displaystyle
u_{\varepsilon}(z+{\varepsilon}x)=u_{0}(z+{\varepsilon}x)+{\varepsilon}K(x)+o({\varepsilon})\quad\mbox{
as }{\varepsilon}\searrow 0,$ (1.4)

where $u_{\varepsilon}$ is defined by the perturbed PDE operator,
$e(\Omega_{\varepsilon},u_{\varepsilon})=0$, and $K$ is a corrector function
defined as the solution to a transmission problem on the unbounded domain
${\mathbf{R}}^{d}$. Often, also a similar asymptotic analysis of an adjoint
variable is needed. There exist different approaches for the derivation of
topological derivatives for PDE-constrained optimization problems of the type
(1.3). We mention the approaches by Novotny and Sokolowski [20], the approach
by Amstutz [3] as well as the averaged adjoint approach [23]. Finally, we
mention the approach that was introduced by Delfour in [10] and applied to
quasilinear problems in [12, 14], which does not involve an asymptotic
analysis of an adjoint variable.

For many problems involving linear or semilinear PDE constraints, topological
derivatives are rather well-understood. In the cases where the shape of the
inclusion $\omega$ is a disk or ellipse in 2D or a ball or ellipsoid in 3D,
the corrector function $K$ can be determined analytically and the topological
derivative can be obtained in a closed form involving so-called polarisation
tensors [2]. When $\omega$ has a more general shape as well as in the case of
quasilinear PDE constraints, explicit formulas have not been reported in the
literature. Topological derivatives in the context of quasilinear PDE
constraints were first treated theoretically in [5], and later also in [12].
The numerical computation for a quasilinear problem was first discussed in [6]
in the context of two-dimensional nonlinear magnetostatics, and later in the
context of three-dimensional nonlinear magnetostatics in an $H(curl)$ setting
in [14].

In this paper, we make the assumption that an expansion of the form (1.4)
holds and apply the systematic approach of M.C. Delfour [10], see also [12],
to a large class of problems involving linear and nonlinear elliptic PDEs as
constraints. The suggested procedure covers scalar and vector-valued problems
in two and three space dimensions. In particular, we will also treat the
problem of nonlinear elasticity in two space dimensions using a St. Venant-
Kirchhoff material which, to the best of the authors’ knowledge, has not been
reported on in the literature so far.

We stress that our generic automated approach is formal since a rigorous
derivation would require a detailed asymptotic analysis of the variation of
the state (1.4). Usually, this entails to show a sufficiently fast decay of
the corrector $K$ as $|x|\rightarrow\infty$. While the computation is formal,
we show numerically that, for the considered problems, the computed values
satisfy a topologically asymptotic expansion of the form (1.2).

##  2\. Lagrangian approach for computing topological derivatives

In this section, we present the general approach for computing topological
derivatives of a class of model problems based on the Lagrangian approach
introduced in [10]. The considered problems include linear and nonlinear,
scalar and vector-valued elliptic PDE constraints in two or three space
dimensions.

###  2.1. Class of considered problems

We consider a PDE-constrained shape optimization problem in $d$ space
dimensions where the solution to the PDE is ${\mathbf{R}}^{m}$-valued. This
covers the case of scalar quantities, $m=1$ or vector-valued problems such as
elasticity, $m=d$. For a function
$\varphi:{\mathbf{R}}^{d}\rightarrow{\mathbf{R}}^{m}$, we denote by
$\mathcal{D}\varphi\in{\mathbf{R}}^{m\times d}$ its Jacobian,
$(\mathcal{D}\varphi)_{i,k}=\frac{\partial\varphi_{i}}{\partial x_{k}}$ for
$i\in\\{1,\dots m\\},k\in\\{1,\dots d\\}$. We recall the Euclidean vector
product $a\cdot b=\sum_{i=1}^{m}a_{i}b_{i}$ for $a,b\in{\mathbf{R}}^{m}$ and
will also use this notation for $m=1$, noting that the vector product is just
a simple product then. Moreover, we denote by
$A:B=\sum_{i=1}^{m}\sum_{k=1}^{d}A_{i,k}B_{i,k}$ the Frobenius inner product
of two matrices $A,B\in{\mathbf{R}}^{m\times d}$ and again note that the
Frobenius inner product reduces to the Euclidean inner product of two vectors
when $m=1$.

We consider a computational domain $\mathsf{D}$ which is subdivided into two
open disjoint subdomains,
$\mathsf{D}=\Omega^{\text{in}}\cup\Omega^{\text{out}}$. We introduce the
abbreviation $\Omega:=\Omega^{\text{in}}$, such that
$\Omega^{\text{out}}=\mathsf{D}\setminus\overline{\Omega}$. Moreover, we
consider four operators

$\displaystyle A_{1}^{\text{in}},A_{1}^{\text{out}}:$
$\displaystyle\,\mathsf{D}\times{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow{\mathbf{R}}^{m}$ $\displaystyle
A_{2}^{\text{in}},A_{2}^{\text{out}}:$
$\displaystyle\,\mathsf{D}\times{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow{\mathbf{R}}^{m\times d}$

to define the two piecewise defined operators
$A_{1}^{\Omega}:\mathsf{D}\times{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow{\mathbf{R}}^{m}$,
$A_{2}^{\Omega}:\mathsf{D}\times{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow{\mathbf{R}}^{m\times d}$,

$\displaystyle
A_{1}^{\Omega}(x,y_{1},y_{2}):=\chi_{\Omega}(x)A_{1}^{\text{in}}(x,y_{1},y_{2})+\chi_{\mathsf{D}\setminus\Omega}(x)A_{1}^{\text{out}}(x,y_{1},y_{2}),$
(2.1) $\displaystyle
A_{2}^{\Omega}(x,y_{1},y_{2}):=\chi_{\Omega}(x)A_{2}^{\text{in}}(x,y_{1},y_{2})+\chi_{\mathsf{D}\setminus\Omega}(x)A_{2}^{\text{out}}(x,y_{1},y_{2}),$
(2.2)

which will represent the left hand side of an abstract PDE constraint. The
right hand side will comprise
$F_{1}^{\Omega}:\mathsf{D}\rightarrow{\mathbf{R}}^{m}$,
$F_{2}^{\Omega}:\mathsf{D}\rightarrow{\mathbf{R}}^{m\times d}$ with

$\displaystyle
F_{i}^{\Omega}(x):=\chi_{\Omega}(x)F_{i}^{\text{in}}(x)+\chi_{\mathsf{D}\setminus\Omega}(x)F_{i}^{\text{out}}(x),$
(2.3)

for $i=1,2$ with functions
$F_{1}^{\text{in}},F_{1}^{\text{out}}:\mathsf{D}\rightarrow{\mathbf{R}}^{m}$
and
$F_{2}^{\text{in}},F_{2}^{\text{out}}:\mathsf{D}\rightarrow{\mathbf{R}}^{m\times
d}$. Finally, we assume that the boundary $\partial\mathsf{D}$ of $\mathsf{D}$
is subdivided into two subsets $\Gamma_{D},\Gamma_{N}$ and consider a function
$g_{N}:\Gamma_{N}\rightarrow{\mathbf{R}}^{m}$ to represent inhomogeneous
Neumann boundary conditions. For sake of more compact presentation, we only
consider homogeneous Dirichlet conditions on the Dirichlet boundary
$\Gamma_{D}$ and remark that an extension to inhomogenous conditions can be
obtained by minor modifications. For a given admissible subdomain
$\Omega\subset\mathsf{D}$, we consider the PDE constraint to find
$u\in{\mathbf{E}}(\mathsf{D})$ such that

$\displaystyle\int_{\mathsf{D}}A_{1}^{\Omega}(x,u,\mathcal{D}u)\cdot\psi+A_{2}^{\Omega}(x,u,\mathcal{D}u):\mathcal{D}\psi\;\mbox{d}x=\int_{\mathsf{D}}F_{1}^{\Omega}(x)\cdot\psi+F_{2}^{\Omega}(x):\mathcal{D}\psi\;\mbox{d}x+\int_{\Gamma_{N}}g_{N}\cdot\psi\;\mbox{d}S$
(2.4)

for all $\psi\in{\mathbf{E}}(\mathsf{D})$, where ${\mathbf{E}}(\mathsf{D})$ is
the function space on which the PDE is posed and which includes the
homogeneous Dirichlet conditions on $\Gamma_{D}$. Similarly, we consider two
functions
$j^{\text{in}},j^{\text{out}}:\mathsf{D}\times{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow{\mathbf{R}}$ and, for a given admissible set
$\Omega\subset\mathsf{D}$, define the piecewise defined function
$j^{\Omega}:\mathsf{D}\times{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow{\mathbf{R}}$,

$\displaystyle
j^{\Omega}(x,y_{1},y_{2}):=\chi_{\Omega}(x)j^{\text{in}}(x,y_{1},y_{2})+\chi_{\mathsf{D}\setminus\Omega}(x)j^{\text{out}}(x,y_{1},y_{2}).$
(2.5)

For $\varphi\in{\mathbf{E}}(\mathsf{D})$, we will consider the cost function

$\displaystyle
J(\Omega,\varphi,\mathcal{D}\varphi):=J^{\text{vol}}(\Omega,\varphi,\mathcal{D}\varphi)+J^{\text{bnd}}(\varphi,\mathcal{D}\varphi)$
(2.6)

with

$\displaystyle
J^{\text{vol}}(\Omega,\varphi,\mathcal{D}\varphi):=\int_{\mathsf{D}}j^{\Omega}(x,\varphi,\mathcal{D}\varphi)\;dx,\qquad
J^{\text{bnd}}(\varphi,\mathcal{D}\varphi):=\int_{\partial\mathsf{D}}j^{bnd}(x,\varphi,\mathcal{D}\varphi)\;dS$
(2.7)

for a function
$j^{bnd}:\partial\mathsf{D}\times{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow{\mathbf{R}}$. Summarizing, we consider the abstract class of
PDE-constrained topology optimization problem which can be written as

$\displaystyle\begin{aligned} &\underset{\Omega\in\mathcal{A}}{\mbox{min
}}J^{\text{vol}}(\Omega,u,\mathcal{D}u)+J^{\text{bnd}}(u,\mathcal{D}u)\\\
&\mbox{subject to }u\in{\mathbf{E}}(\mathsf{D})\mbox{ solves
}\eqref{eq_abstr_pde}\end{aligned}$ (2.8)

where $\mathcal{A}$ denotes the set of admissible subsets of $\mathsf{D}$.

######  Remark 1.

Although we consider only homogeneous boundary conditions, the formulas for
inhomogeneous Dirichlet boundary conditions will not affect the formula of the
first topological derivative; see, e.g., [8] where mixed (Dirichlet and
Neumann) inhomogeneous boundary conditions are considered.

###  2.2. Domain perturbation

Let $\omega\subset{\mathbf{R}}^{d}$ with $0\in\omega$ represent the shape of
the considered topological perturbation. For a point
$z\in\mathsf{D}\setminus\partial\Omega$ and a given small parameter
${\varepsilon}$, we define $\omega_{\varepsilon}(z):=z+{\varepsilon}\omega$ as
well as the perturbed domain

$\displaystyle\Omega_{\varepsilon}(z):=\begin{cases}\Omega\setminus\overline{\omega}_{\varepsilon}(z),&z\in\Omega,\\\
\Omega\cup\omega_{\varepsilon}(z),&z\in\mathsf{D}\setminus\overline{\Omega}.\end{cases}$

From now on, we assume that a fixed point
$z\in\mathsf{D}\setminus\overline{\Omega}$ is given and set
$\omega_{\varepsilon}:=\omega_{\varepsilon}(z)$,
$\Omega_{\varepsilon}:=\Omega_{\varepsilon}(z)$. Moreover, we use the
abbreviation $\mathsf{D}_{\varepsilon}:=T_{\varepsilon}^{-1}(\mathsf{D})$,
where $T_{\varepsilon}(x):=z+{\varepsilon}x$. We also introduce the
abbreviating notations
$A_{1}^{({\varepsilon})}:=A_{1}^{\Omega_{\varepsilon}}$,
$A_{2}^{({\varepsilon})}:=A_{2}^{\Omega_{\varepsilon}}$,
$F_{1}^{({\varepsilon})}:=F_{1}^{\Omega_{\varepsilon}}$,
$F_{2}^{({\varepsilon})}:=F_{2}^{\Omega_{\varepsilon}}$,
$j^{({\varepsilon})}:=j^{\Omega_{\varepsilon}}$ and define the perturbed
Lagrangian of the abstract optimization problem (2.8)

$\displaystyle
G({\varepsilon},\varphi,\mathcal{D}\varphi,\psi,\mathcal{D}\psi)=$
$\displaystyle\,J(\Omega_{\varepsilon},\varphi,\mathcal{D}\varphi)+\langle
A_{1}^{({\varepsilon})}(\varphi,\mathcal{D}\varphi),\psi\rangle+\langle
A_{2}^{({\varepsilon})}(\varphi,\mathcal{D}\varphi),\mathcal{D}\psi\rangle$
$\displaystyle-\langle F_{1}^{({\varepsilon})},\psi\rangle-\langle
F_{2}^{({\varepsilon})},\mathcal{D}\psi\rangle-\langle g_{N},\psi\rangle,$

where $\langle\cdot,\cdot\rangle$ denotes the duality product in the
corresponding spaces.

###  2.3. Formal variation of the state variable

The following derivation is motivated by the rigorous results obtained in
[12], where a simplified approach to the computation of quasi-linear
derivatives is proposed. Here we generalise these results formally to the
class of non-linear problems introduced in (2.8). We also refer to [6, 5] for
another approach to the computation of topological derivatives for quasi-
linear problems.

Let $\psi\in{\mathbf{E}}(\mathsf{D})$. The perturbed state equation reads:
$u_{\varepsilon}$ solves

$\displaystyle\begin{split}\int_{\mathsf{D}}A_{1}^{({\varepsilon})}(x,u_{\varepsilon},\mathcal{D}u_{\varepsilon})\cdot\psi+A_{2}^{({\varepsilon})}(x,u_{\varepsilon},\mathcal{D}u_{\varepsilon}):\mathcal{D}\psi\;dx=\int_{\mathsf{D}}F_{1}^{({\varepsilon})}(x)\cdot\psi&+F_{2}^{({\varepsilon})}(x):\mathcal{D}\psi\;dx\\\
&+\int_{\Gamma_{N}}g_{N}\cdot\psi\;dS.\end{split}$ (2.9)

Subtracting the state equation (2.9) for ${\varepsilon}>0$ from the state
equation for ${\varepsilon}=0$ leads to

$\displaystyle\begin{split}\int_{\mathsf{D}}[A_{1}^{({\varepsilon})}(x,u_{\varepsilon},\mathcal{D}u_{\varepsilon})-&A_{1}^{(0)}(x,u_{0},\mathcal{D}u_{0})]\cdot\psi+[A_{2}^{({\varepsilon})}(x,u_{\varepsilon},\mathcal{D}u_{\varepsilon})-A_{2}^{(0)}(x,u_{0},\mathcal{D}u_{0})]:\mathcal{D}\psi\;dx\\\
&=\int_{\omega_{\varepsilon}}(F_{1}^{\text{in}}-F_{1}^{\text{out}})(x)\cdot\psi+(F_{2}^{\text{in}}-F_{2}^{\text{out}})(x):\mathcal{D}\psi\;dx.\end{split}$

Using the change of variables $T_{\varepsilon}(x)=z+{\varepsilon}x$, we obtain
with the notation $\mathsf{D}_{\varepsilon}=T_{\varepsilon}^{-1}(\mathsf{D})$:

$\displaystyle\int_{\mathsf{D}_{\varepsilon}}$
$\displaystyle[A_{1}^{({\varepsilon})}(T_{\varepsilon}(x),u_{\varepsilon}\circ
T_{\varepsilon},(\mathcal{D}u_{\varepsilon})\circ
T_{\varepsilon})-A_{1}^{(0)}(T_{\varepsilon}(x),u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ T_{\varepsilon})]\cdot(\psi\circ
T_{\varepsilon})$ $\displaystyle+$
$\displaystyle[A_{2}^{({\varepsilon})}(T_{\varepsilon}(x),u_{\varepsilon}\circ
T_{\varepsilon},(\mathcal{D}u_{\varepsilon})\circ
T_{\varepsilon})-A_{2}^{(0)}(T_{\varepsilon}(x),u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})]:(\mathcal{D}\psi)\circ T_{\varepsilon}\;dx$
$\displaystyle=\int_{\omega}((F_{1}^{\text{in}}-F_{1}^{\text{out}})\circ
T_{\varepsilon})\cdot(\psi\circ
T_{\varepsilon})+((F_{2}^{\text{in}}-F_{2}^{\text{out}})\circ
T_{\varepsilon}):((\mathcal{D}\psi)\circ T_{\varepsilon})\;dx.$

Using $(\mathcal{D}\varphi)\circ
T_{\varepsilon}=\frac{1}{{\varepsilon}}\mathcal{D}(\varphi\circ
T_{\varepsilon})$ and the first variation of the state
$K_{\varepsilon}:=\frac{(u_{\varepsilon}-u_{0})\circ
T_{\varepsilon}}{{\varepsilon}}$, after multiplication with ${\varepsilon}$,
we obtain

$\displaystyle\int_{\mathsf{D}_{\varepsilon}}$
$\displaystyle[A_{1}^{({\varepsilon})}(T_{\varepsilon}(x),u_{0}\circ
T_{\varepsilon}+{\varepsilon}K_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon}+\mathcal{D}K_{\varepsilon})$ $\displaystyle\hskip 113.81102pt-
A_{1}^{(0)}(T_{\varepsilon}(x),u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})]\cdot({\varepsilon}\psi\circ T_{\varepsilon})$
$\displaystyle+$
$\displaystyle\int_{\mathsf{D}_{\varepsilon}}[A_{2}^{({\varepsilon})}(T_{\varepsilon}(x),u_{0}\circ
T_{\varepsilon}+{\varepsilon}K_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon}+\mathcal{D}K_{\varepsilon})$ $\displaystyle\hskip 142.26378pt-
A_{2}^{(0)}(T_{\varepsilon}(x),u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})]:\mathcal{D}(\psi\circ T_{\varepsilon})\;dx$
$\displaystyle=\int_{\omega}((F_{1}^{\text{in}}-F_{1}^{\text{out}})\circ
T_{\varepsilon})\cdot({\varepsilon}\psi\circ
T_{\varepsilon})+((F_{2}^{\text{in}}-F_{2}^{\text{out}})\circ
T_{\varepsilon}):\mathcal{D}(\psi\circ T_{\varepsilon})\;dx.$

Now we make the following assumption:

######  Assumption 1.

  * (i)

We assume that $u_{0}$ is continuously differentiable at $z$.

  * (ii)

Let $K_{\varepsilon}:=\frac{(u_{\varepsilon}-u_{0})\circ
T_{\varepsilon}}{{\varepsilon}}$ for ${\varepsilon}>0$ be the first variation
of the state variable $u_{\varepsilon}$. We assume that $\nabla
K_{\varepsilon}\to\nabla K$ and ${\varepsilon}K_{\varepsilon}\to 0$, where $K$
solves (2.11).

  * (iii)

For all ${\varepsilon}>0$ we have
$\psi\in{\mathbf{E}}(\mathsf{D}_{\varepsilon})$ if and only if $\psi\circ
T_{\varepsilon}\in{\mathbf{E}}(\mathsf{D})$.

Rearranging and replacing $\psi\circ T_{\varepsilon}$ by $\psi$ (using
Assumption 1,item (iii)) yields

$\begin{split}\int_{\mathsf{D}_{\varepsilon}}&[A_{1}^{({\varepsilon})}(T_{\varepsilon}(x),u_{0}\circ
T_{\varepsilon}+{\varepsilon}K_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon}+\mathcal{D}K_{\varepsilon})-A_{1}^{({\varepsilon})}(T_{\varepsilon}(x),u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})]\cdot{\varepsilon}\psi\\\
+&[A_{2}^{({\varepsilon})}(T_{\varepsilon}(x),u_{0}\circ
T_{\varepsilon}+{\varepsilon}K_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon}+\mathcal{D}K_{\varepsilon})-A_{2}^{({\varepsilon})}(T_{\varepsilon}(x),u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})]:\mathcal{D}\psi\;dx\\\
=&\int_{\omega}((F_{1}^{\text{in}}-F_{1}^{\text{out}})\circ
T_{\varepsilon})\cdot({\varepsilon}\psi)+((F_{2}^{\text{in}}-F_{2}^{\text{out}})\circ
T_{\varepsilon}):\mathcal{D}\psi\;dx\\\
&-\int_{\omega}[A_{1}^{\text{in}}-A_{1}^{\text{out}}](T_{\varepsilon}(x),u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})\cdot{\varepsilon}\psi\;dx\\\
&-\int_{\omega}[A_{2}^{\text{in}}-A_{2}^{\text{out}}](T_{\varepsilon}(x),u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon}):\mathcal{D}\psi\;dx\end{split}$ (2.10)

for all $\psi\in{\mathbf{E}}(\mathsf{D}_{\varepsilon})$. With this we can pass
to limit ${\varepsilon}\rightarrow 0$ in (2.10) and get

$\displaystyle\begin{split}\int_{{\mathbf{R}}^{d}}&[A_{2}^{\omega}(z,u_{0}(z),\mathcal{D}u_{0}(z)+\mathcal{D}K)-A_{2}^{\omega}(z,u_{0}(z),\mathcal{D}u_{0}(z))]:\mathcal{D}\psi\;dx\\\
&=\int_{\omega}(F_{2}^{\text{in}}-F_{2}^{\text{out}})(z):\mathcal{D}\psi\;dx-\int_{\omega}[A_{2}^{\text{in}}-A_{2}^{\text{out}}](z,u_{0}(z),\mathcal{D}u_{0}(z)):\mathcal{D}\psi\;dx,\end{split}$
(2.11)

for all $\psi\in{\mathbf{E}}({\mathbf{R}}^{d}):=\lim_{{\varepsilon}\searrow
0}{\mathbf{E}}(\mathsf{D}_{\varepsilon})$ , with the definition
$A_{2}^{\omega}(x,y_{1},y_{2})=\chi_{\omega}(x)A_{2}^{\text{in}}(x,y_{1},y_{2})+\chi_{{\mathbf{R}}^{d}\setminus\omega}(x)A_{2}^{\text{out}}(x,y_{1},y_{2})$.
The limit $\lim_{{\varepsilon}\searrow
0}{\mathbf{E}}(\mathsf{D}_{\varepsilon})$ has to be understood formally and in
practice often is a Beppo-Levi space (see the next remark).

######  Remark 2.

The function space ${\mathbf{E}}({\mathbf{R}}^{d})$ in which (2.11) admits a
solution can often be chosen as a Beppo-Levi space as shown in [23], see also
[12]. For the numerical approximation in the next section we use precisely
that the functions $K_{\varepsilon}$ converge to $K$ and approximate $K$ by
solving an equation on a blown up domain similarly to
$\mathsf{D}_{\varepsilon}$.

###  2.4. Lagrangian theorem for first topological derivatives

In this section we discuss a method proposed by M.C. Delfour in [10, Thm.3.3].
The definite advantage is that it uses the unperturbed adjoint equation and
only requires the asymptotic analysis of the state equation, but it seems to
come with the shortcoming that certain cost functions cannot be treated
rigorously; see [12] and also [8].

Let ${\mathbf{E}}(\mathsf{D})$ be a Banach space of functions on $\mathsf{D}$.
For all parameter ${\varepsilon}\geq 0$ small consider a function
$u_{\varepsilon}\in{\mathbf{E}}(\mathsf{D})$ solving the variational problem
of the form

$a_{\varepsilon}[u_{\varepsilon}](\varphi)=f_{\varepsilon}(\varphi)\quad\text{
for all }\varphi\in{\mathbf{E}}(\mathsf{D}),$ (2.12)

where
$a_{\varepsilon}:{\mathbf{E}}(\mathsf{D})\rightarrow{\mathbf{E}}(\mathsf{D})^{\prime}$
is a nonlinear operator and $f_{\varepsilon}$ is a linear form on
${\mathbf{E}}(\mathsf{D})$, respectively. Throughout we assume that this
abstract state equation admits a unique solution and that
$u_{\varepsilon}-u_{0}\in{\mathbf{E}}(\mathsf{D})$ for all ${\varepsilon}$.
Consider now a cost function

$j({\varepsilon})=J_{\varepsilon}(u_{\varepsilon})\in{\mathbf{R}},$ (2.13)

where for all ${\varepsilon}\geq 0$ the functional
$J_{\varepsilon}:{\mathbf{E}}(\mathsf{D})\to{\mathbf{R}}$ is differentiable at
$u_{0}$. In the following sections we review methods how to obtain an
asymptotic expansion of $j({\varepsilon})$ at ${\varepsilon}=0$. For this
purpose we introduce the Lagrangian function

$\mathfrak{L}({\varepsilon},u,v)=J_{\varepsilon}(u)+a_{\varepsilon}[u](v)-f_{\varepsilon}(v),\quad
u\in{\mathbf{E}}(\mathsf{D}),\;v\in{\mathbf{E}}(\mathsf{D}).$

######  Proposition 1 ([10]).

Let $\ell:[0,\tau]\to{\mathbf{R}}$ be a function with $\ell({\varepsilon})>0$
for ${\varepsilon}>0$ and $\displaystyle\lim_{{\varepsilon}\searrow
0}\ell({\varepsilon})=0$. Furthermore, assume that the limits

$\displaystyle\mathfrak{R}_{1}(u_{0},p_{0}):=$
$\displaystyle\underset{{\varepsilon}\searrow
0}{\lim}\;\frac{1}{\ell({\varepsilon})}\int_{0}^{1}\left[\partial_{u}\mathfrak{L}({\varepsilon},su_{\varepsilon}+(1-s)u_{0},p_{0})-\partial_{u}\mathfrak{L}({\varepsilon},u_{0},p_{0})(u_{\varepsilon}-u_{0})\right]\;ds,$
(2.14) $\displaystyle\mathfrak{R}_{2}(u_{0},p_{0}):=$
$\displaystyle\underset{{\varepsilon}\searrow
0}{\lim}\;\frac{1}{\ell({\varepsilon})}(\partial_{u}\mathfrak{L}({\varepsilon},u_{0},p_{0})-\partial_{u}\mathfrak{L}(0,u_{0},p_{0}))(u_{\varepsilon}-u_{0}),$
(2.15) $\displaystyle\partial_{\ell}\mathfrak{L}(0,u_{0},p_{0}):=$
$\displaystyle\lim_{{\varepsilon}\searrow
0}\;\frac{1}{\ell({\varepsilon})}(\mathfrak{L}({\varepsilon},u_{0},p_{0})-\mathfrak{L}(0,u_{0},p_{0})),$
(2.16)

exist. Then the following expansion holds:

$\displaystyle
j({\varepsilon})=j(0)+\ell({\varepsilon})((\mathfrak{R}_{1}(u_{0},p_{0})+\mathfrak{R}_{2}(u_{0},p_{0})+\partial_{\ell}\mathfrak{L}(0,u_{0},p_{0}))+o(\ell({\varepsilon})).$
(2.17)

######  Remark 3.

Note that by the fundamental theorem of calculus $\mathfrak{R}_{1}$ can be
equivalently written as

$\mathfrak{R}_{1}(u_{0},p_{0})=\underset{{\varepsilon}\searrow
0}{\lim}\;\frac{1}{\ell({\varepsilon})}\left[\mathfrak{L}({\varepsilon},u_{\varepsilon},p_{0})-\mathfrak{L}({\varepsilon},u_{0},p_{0})-\partial_{u}\mathfrak{L}({\varepsilon},u_{0},p_{0})(u_{\varepsilon}-u_{0})\right].$

However, the form stated in the proposition is better suited for the formal
computation of the topological derivatives presented later on.

###  2.5. Formal computation of topological derivative

We are now applying formally Proposition 1 to compute the topological
derivative of problem (2.8). For this we let
$\mathfrak{L}({\varepsilon},\varphi,\psi):=G({\varepsilon},\varphi,\mathcal{D}\varphi,\psi,\mathcal{D}\psi)$.
For sake of better readability, from now on we drop the dependence of
$A_{i}^{({\varepsilon})}$, $i=1,2$ and $j^{({\varepsilon})}$ on the space
variable $x$. In what follows we compute the terms
$\mathfrak{R}_{1}(u_{0},p_{0})$, $\mathfrak{R}_{2}(u_{0},p_{0})$ and
$\partial_{\ell}\mathfrak{L}(0,u_{0},p_{0})$ of Proposition 1 for the
Lagrangian $\mathfrak{L}$ separately as follows. For this purpose we first
introduce the following abbreviations:

$\begin{aligned} a_{1}({\varepsilon},\varphi,\psi)&:=\langle
A_{1}^{({\varepsilon})}(\varphi,\mathcal{D}\varphi),\psi\rangle,\\\
f_{1}({\varepsilon},\psi)&:=\langle F_{1}^{({\varepsilon})},\psi\rangle,\\\
J_{1}({\varepsilon},\psi)&:=J^{\text{vol}}(\Omega_{\varepsilon},\psi,\mathcal{D}\psi),\end{aligned}\begin{aligned}
a_{2}({\varepsilon},\varphi,\psi)&:=\langle
A_{2}^{({\varepsilon})}(\varphi,\mathcal{D}\varphi),\mathcal{D}\psi\rangle,\\\
f_{2}({\varepsilon},\psi)&:=\langle
F_{2}^{({\varepsilon})},\mathcal{D}\psi\rangle,\\\
J_{2}(\psi)&:=J^{\text{bnd}}(\psi,\mathcal{D}\psi).\end{aligned}$

Then we split for instance $\mathfrak{R}_{1}(u_{0},p_{0})$ into the parts
coming from $a_{1},a_{2},f_{1},f_{2},J_{1},J_{2}$. The term coming from
$a_{1}$ contributing to $\mathfrak{R}_{1}(u_{0},p_{0})$ is for instance given
by

$\mathfrak{R}_{1}^{A_{1}}(u_{0},p_{0}):=\lim_{{\varepsilon}\searrow
0}\frac{1}{|\omega_{\varepsilon}|}\int_{0}^{1}\left[\partial_{\varphi}a_{1}({\varepsilon},u_{0}+s(u_{\varepsilon}-u_{0}),p_{0})-\partial_{\varphi}a_{1}({\varepsilon},u_{0},p_{0})\right](u_{\varepsilon}-u_{0})\;ds,$
(2.18)

where $\partial_{\varphi}a_{1}$ denotes the derivative with respect to the
second argument. Similarly we denote by
$\mathfrak{R}_{1}^{A_{2}}(u_{0},p_{0})$ and
$\mathfrak{R}_{1}^{J_{1}}(u_{0},p_{0})$ the contributions of $a_{2}$ and
$J_{1}$ to the term $\mathfrak{R}_{1}(u_{0},p_{0})$, respectively. We proceed
in the same fashion for each term. The detailed computations are outlined in
the following subsections. Additionally to Assumption 1 we make the following
assumption.

######  Assumption 2.

We assume that $p_{0}$ is continuously differentiable at $z$.

####  2.5.1. Terms coming from $A_{1}^{({\varepsilon})}$

For $A_{1}^{({\varepsilon})}:{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow{\mathbf{R}}^{m}$ let
$\partial_{y_{1}}A_{1}^{({\varepsilon})}:{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow\mathcal{L}({\mathbf{R}}^{m},{\mathbf{R}}^{m})$ denote the
derivative with respect to the first argument and
$\partial_{y_{2}}A_{1}^{({\varepsilon})}:{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow\mathcal{L}({\mathbf{R}}^{m\times d},{\mathbf{R}}^{m})$ the
derivative with respect to the second argument. We compute

$\displaystyle\mathfrak{R}_{1,{\varepsilon}}^{A_{1}}:=\int_{0}^{1}\left[\partial_{\varphi}a_{1}({\varepsilon},u_{0}+s(u_{\varepsilon}-u_{0}),p_{0})-\partial_{\varphi}a_{1}({\varepsilon},u_{0},p_{0})\right](u_{\varepsilon}-u_{0})\;ds$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\int_{\mathsf{D}}[\partial_{y_{1}}A_{1}^{({\varepsilon})}(u_{0}+s(u_{\varepsilon}-u_{0}),\mathcal{D}u_{0}+s\mathcal{D}(u_{\varepsilon}-u_{0}))-\partial_{y_{1}}A_{1}^{({\varepsilon})}(u_{0},\mathcal{D}u_{0})](u_{\varepsilon}-u_{0})\cdot
p_{0}\;dx\,ds$
$\displaystyle+\int_{0}^{1}\int_{\mathsf{D}}[\partial_{y_{2}}A_{1}^{({\varepsilon})}(u_{0}+s(u_{\varepsilon}-u_{0}),\mathcal{D}u_{0}+s\mathcal{D}(u_{\varepsilon}-u_{0}))$
$\displaystyle\hskip
142.26378pt-\partial_{y_{2}}A_{1}^{({\varepsilon})}(u_{0},\mathcal{D}u_{0})](\mathcal{D}(u_{\varepsilon}-u_{0}))\cdot
p_{0}\;dx\,ds.$

Thus changing variables and using $(\mathcal{D}\varphi)\circ
T_{\varepsilon}=\frac{1}{{\varepsilon}}\mathcal{D}(\varphi\circ
T_{\varepsilon})$ and $K_{\varepsilon}:=\frac{(u_{\varepsilon}-u_{0})\circ
T_{\varepsilon}}{{\varepsilon}}$ yields

$\displaystyle\mathfrak{R}_{1,{\varepsilon}}^{A_{1}}=$
$\displaystyle{\varepsilon}^{d}\int_{0}^{1}\int_{\mathsf{D}_{\varepsilon}}[\partial_{y_{1}}A_{1}^{\omega}(u_{0}\circ
T_{\varepsilon}+s{\varepsilon}K_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon}+s\mathcal{D}K_{\varepsilon})$ $\displaystyle\hskip
142.26378pt-\partial_{y_{1}}A_{1}^{\omega}(u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})]({\varepsilon}K_{\varepsilon})\cdot(p_{0}\circ
T_{\varepsilon})\;dx\,ds$
$\displaystyle+{\varepsilon}^{d}\int_{0}^{1}\int_{\mathsf{D}_{\varepsilon}}[\partial_{y_{2}}A_{1}^{\omega}(u_{0}\circ
T_{\varepsilon}+s{\varepsilon}K_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon}+s\mathcal{D}K_{\varepsilon})$ $\displaystyle\hskip
142.26378pt-\partial_{y_{2}}A_{1}^{\omega}(u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})](\mathcal{D}K_{\varepsilon})\cdot(p_{0}\circ
T_{\varepsilon})\;dx\,ds.$

Therefore setting
$\mathfrak{R}_{1}^{A_{1}}(u_{0},p_{0})=\underset{{\varepsilon}\searrow
0}{\lim}\frac{1}{|\omega_{\varepsilon}|}\mathfrak{R}_{1,{\varepsilon}}^{A_{1}}$,
we obtain

$\displaystyle\mathfrak{R}_{1}^{A_{1}}(u_{0},p_{0})=\frac{1}{|\omega|}\int_{{\mathbf{R}}^{d}}[A_{1}^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z)+\mathcal{D}K)$
$\displaystyle-A_{1}^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z))$
$\displaystyle-\partial_{y_{2}}A_{1}^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z))(\mathcal{D}K)]\cdot
p_{0}(z)\;dx.$

Next, we compute

$\displaystyle\mathfrak{R}_{2,{\varepsilon}}^{A_{1}}:=$
$\displaystyle\left[\partial_{\varphi}a_{1}({\varepsilon},u_{0},p_{0})-\partial_{\varphi}a_{1}(0,u_{0},p_{0})\right](u_{\varepsilon}-u_{0})$
$\displaystyle=$
$\displaystyle\int_{\mathsf{D}}[\partial_{y_{1}}A_{1}^{({\varepsilon})}(u_{0},\mathcal{D}u_{0})-\partial_{y_{1}}A_{1}^{(0)}(u_{0},\mathcal{D}u_{0})](u_{\varepsilon}-u_{0})\cdot
p_{0}\;dx$
$\displaystyle+\int_{\mathsf{D}}[\partial_{y_{2}}A_{1}^{({\varepsilon})}(u_{0},\mathcal{D}u_{0})-\partial_{y_{2}}A_{1}^{(0)}(u_{0},\mathcal{D}u_{0})](\mathcal{D}(u_{\varepsilon}-u_{0}))\cdot
p_{0}\;dx.$

Thus changing variables and using $(\mathcal{D}\varphi)\circ
T_{\varepsilon}=\frac{1}{{\varepsilon}}\mathcal{D}(\varphi\circ
T_{\varepsilon})$ and $K_{\varepsilon}:=\frac{(u_{\varepsilon}-u_{0})\circ
T_{\varepsilon}}{{\varepsilon}}$ yields

$\displaystyle\mathfrak{R}_{2,{\varepsilon}}^{A_{1}}=$
$\displaystyle{\varepsilon}^{d}\int_{\omega}[\partial_{y_{1}}A_{1}^{\text{in}}(u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})-\partial_{y_{1}}A_{1}^{\text{out}}(u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})]({\varepsilon}K_{\varepsilon})\cdot(p_{0}\circ
T_{\varepsilon})\;dx$
$\displaystyle+{\varepsilon}^{d}\int_{\omega}[\partial_{y_{2}}A_{1}^{\text{in}}(u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})-\partial_{y_{2}}A_{1}^{\text{out}}(u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})](\mathcal{D}K_{\varepsilon})\cdot(p_{0}\circ
T_{\varepsilon})\;dx.$

Therefore setting
$\mathfrak{R}_{2}^{A_{1}}(u_{0},p_{0})=\underset{{\varepsilon}\searrow
0}{\lim}\frac{1}{|\omega_{\varepsilon}|}\mathfrak{R}_{1,{\varepsilon}}^{A_{1}}$,
we obtain

$\displaystyle\mathfrak{R}_{2}^{A_{1}}(u_{0},p_{0})=\frac{1}{|\omega|}\int_{\omega}[\partial_{y_{2}}A_{1}^{\text{in}}(u_{0}(z),\mathcal{D}u_{0}(z))-\partial_{y_{2}}A_{1}^{\text{out}}(u_{0}(z),\mathcal{D}u_{0}(z))](\mathcal{D}K)\cdot
p_{0}(z)\;dx.$

Finally we compute the term contributing to
$\partial_{\ell}\mathfrak{L}(0,u_{0},p_{0})$ coming from
$A_{1}^{({\varepsilon})}$ by

$\displaystyle\partial_{\ell}\mathfrak{L}^{A_{1}}(0,u_{0},p_{0})=[A_{1}^{\text{in}}(u_{0}(z),\mathcal{D}u_{0}(z))-A_{1}^{\text{out}}(u_{0}(z),\mathcal{D}u_{0}(z))]\cdot
p_{0}(z).$

####  2.5.2. Terms coming from $A_{2}^{({\varepsilon})}$

For $A_{2}^{({\varepsilon})}:{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow{\mathbf{R}}^{m\times d}$ let
$\partial_{y_{1}}A_{2}^{({\varepsilon})}:{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow\mathcal{L}({\mathbf{R}}^{m},{\mathbf{R}}^{m\times d})$ denote
the derivative with respect to the first argument and
$\partial_{y_{2}}A_{2}^{({\varepsilon})}:{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow\mathcal{L}({\mathbf{R}}^{m\times d},{\mathbf{R}}^{m\times d})$
the derivative with respect to the second argument. We have

$\displaystyle\mathfrak{R}_{1,{\varepsilon}}^{A_{2}}=$
$\displaystyle\int_{0}^{1}\left(\partial_{\varphi}a_{2}({\varepsilon},u_{0}+s(u_{\varepsilon}-u_{0}),p_{0})-\partial_{\varphi}a_{2}({\varepsilon},u_{0},p_{0})\right)(u_{\varepsilon}-u_{0})\;ds$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\int_{\mathsf{D}}[\partial_{y_{1}}A_{2}^{({\varepsilon})}(u_{0}+s(u_{\varepsilon}-u_{0}),\mathcal{D}u_{0}+s\mathcal{D}(u_{\varepsilon}-u_{0}))$
$\displaystyle\hskip
170.71652pt-\partial_{y_{1}}A_{2}^{({\varepsilon})}(u_{0},\mathcal{D}u_{0})](u_{\varepsilon}-u_{0}):\mathcal{D}p_{0}\;dx\,ds$
$\displaystyle+\int_{0}^{1}\int_{\mathsf{D}}[\partial_{y_{2}}A_{2}^{({\varepsilon})}(u_{0}+s(u_{\varepsilon}-u_{0}),\mathcal{D}u_{0}+s\mathcal{D}(u_{\varepsilon}-u_{0}))$
$\displaystyle\hskip
170.71652pt-\partial_{y_{2}}A_{2}^{({\varepsilon})}(u_{0},\mathcal{D}u_{0})](\mathcal{D}(u_{\varepsilon}-u_{0})):\mathcal{D}p_{0}\;dx\,ds.$

Letting
$\mathfrak{R}_{1}^{A_{2}}(u_{0},p_{0}):=\underset{{\varepsilon}\searrow
0}{\lim}\frac{1}{|\omega_{\varepsilon}|}\mathfrak{R}_{1,{\varepsilon}}^{A_{2}}$
and following the same steps as in the computation of
$\mathfrak{R}_{1}^{A_{1}}(u_{0},p_{0})$ above leads to

$\begin{split}\mathfrak{R}_{1}^{A_{2}}(u_{0},p_{0})=\frac{1}{|\omega|}\int_{{\mathbf{R}}^{d}}[A_{2}^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z)+\mathcal{D}K)&-A_{2}^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z))\\\
&-\partial_{y_{2}}A_{2}^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z))(\mathcal{D}K)]:\mathcal{D}p_{0}(z)\;dx.\end{split}$

We further have

$\displaystyle\mathfrak{R}_{2,{\varepsilon}}^{A_{2}}:=$
$\displaystyle\left[\partial_{\varphi}a_{2}({\varepsilon},u_{0})-\partial_{\varphi}a_{2}(0,u_{0})\right](u_{\varepsilon}-u_{0})$
$\displaystyle=$
$\displaystyle\int_{\mathsf{D}}[\partial_{y_{1}}A_{2}^{({\varepsilon})}(u_{0},\mathcal{D}u_{0})-\partial_{y_{1}}A_{2}^{(0)}(u_{0},\mathcal{D}u_{0})](u_{\varepsilon}-u_{0}):\mathcal{D}p_{0}\;dx$
$\displaystyle+\int_{\mathsf{D}}[\partial_{y_{2}}A_{2}^{({\varepsilon})}(u_{0},\mathcal{D}u_{0})-\partial_{y_{2}}A_{2}^{(0)}(u_{0},\mathcal{D}u_{0})](\mathcal{D}(u_{\varepsilon}-u_{0})):\mathcal{D}p_{0}\;dx.$

Letting
$\mathfrak{R}_{2}^{A_{2}}(u_{0},p_{0}):=\underset{{\varepsilon}\searrow
0}{\lim}\frac{1}{|\omega_{\varepsilon}|}\mathfrak{R}_{2,{\varepsilon}}^{A_{2}}$
and following the same steps as in the computation of
$\mathfrak{R}_{2}^{A_{1}}(u_{0},p_{0})$ leads to

$\displaystyle\mathfrak{R}_{2}^{A_{2}}(u_{0},p_{0})$
$\displaystyle=\frac{1}{|\omega|}\int_{\omega}[\partial_{y_{2}}A_{2}^{\text{in}}(u_{0}(z),\mathcal{D}u_{0}(z))-\partial_{y_{2}}A_{2}^{\text{out}}(u_{0}(z),\mathcal{D}u_{0}(z))](\mathcal{D}K):\mathcal{D}p_{0}(z)\;dx.$
(2.19)

Finally, the part of $\partial_{\ell}\mathfrak{L}(0,u_{0},p_{0})$ coming from
$A_{2}^{({\varepsilon})}$ reads

$\displaystyle\partial_{\ell}\mathfrak{L}^{A_{2}}(0,u_{0},p_{0})=[A_{2}^{\text{in}}(u_{0}(z),\mathcal{D}u_{0}(z))-A_{2}^{\text{out}}(u_{0}(z),\mathcal{D}u_{0}(z))]:\mathcal{D}p_{0}(z).$

####  2.5.3. Terms coming from right hand side

Since the right hand side does not depend on the solution and therefore
$\partial_{u}F=0$, there is no contribution to $\mathfrak{R}_{1}$ and
$\mathfrak{R}_{2}$. It only remains

$\displaystyle\partial_{\ell}\mathfrak{L}^{F}(0,u_{0},p_{0})=-(F_{1}^{\text{in}}-F_{1}^{\text{out}})\cdot
p_{0}(z)-(F_{2}^{\text{in}}-F_{2}^{\text{out}}):\mathcal{D}p_{0}(z).$ (2.20)

####  2.5.4. Terms coming from $J^{\text{vol}}$

Let $j_{y_{1}}^{({\varepsilon})}:{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow\mathcal{L}({\mathbf{R}}^{m},{\mathbf{R}})$ denote the derivative
of $j^{({\varepsilon})}$ with respect to the first argument and
$j_{y_{2}}^{({\varepsilon})}:{\mathbf{R}}^{m}\times{\mathbf{R}}^{m\times
d}\rightarrow\mathcal{L}({\mathbf{R}}^{m\times d},{\mathbf{R}})$ the
derivative with respect to the second argument such that

$\displaystyle\partial_{u}J^{\text{vol}}(\Omega_{\varepsilon},u,\mathcal{D}u)(\hat{u})=\int_{\mathsf{D}}j_{y_{1}}^{({\varepsilon})}(u,\mathcal{D}u)(\hat{u})+j_{y_{2}}^{({\varepsilon})}(u,\mathcal{D}u)(\mathcal{D}\hat{u})\;dx.$
(2.21)

We have

$\displaystyle\mathfrak{R}_{1,{\varepsilon}}^{J^{\text{vol}}}:=$
$\displaystyle\int_{0}^{1}(\partial_{u}J_{1}({\varepsilon},u_{0}+s(u_{\varepsilon}-u_{0}))-\partial_{u}J_{1}({\varepsilon},u_{0}))(u_{\varepsilon}-u_{0})\;ds$
$\displaystyle=$
$\displaystyle\int_{0}^{1}[\partial_{u}J^{\text{vol}}(\Omega_{\varepsilon},u_{0}+s(u_{\varepsilon}-u_{0}),\mathcal{D}u_{0}+s\mathcal{D}(u_{\varepsilon}-u_{0}))$
$\displaystyle\hskip
142.26378pt-\partial_{u}J^{\text{vol}}(\Omega_{\varepsilon},u_{0},\mathcal{D}u_{0})](u_{\varepsilon}-u_{0})\;ds$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\Bigg{\\{}\int_{\mathsf{D}}[j^{({\varepsilon})}_{y_{1}}(u_{0}+s(u_{\varepsilon}-u_{0}),\mathcal{D}u_{0}+s\mathcal{D}(u_{\varepsilon}-u_{0}))-j^{({\varepsilon})}_{y_{1}}(u_{0},\mathcal{D}u_{0})](u_{\varepsilon}-u_{0})\;dx$
$\displaystyle+\int_{\mathsf{D}}[j^{({\varepsilon})}_{y_{2}}(u_{0}+s(u_{\varepsilon}-u_{0}),\mathcal{D}u_{0}+s\mathcal{D}(u_{\varepsilon}-u_{0}))-j^{({\varepsilon})}_{y_{2}}(u_{0},\mathcal{D}u_{0})](\mathcal{D}(u_{\varepsilon}-u_{0}))\;dx\Bigg{\\}}\;ds.$

Therefore changing variables and using $(\mathcal{D}\varphi)\circ
T_{\varepsilon}=\frac{1}{{\varepsilon}}\mathcal{D}(\varphi\circ
T_{\varepsilon})$ and $K_{\varepsilon}:=\frac{(u_{\varepsilon}-u_{0})\circ
T_{\varepsilon}}{{\varepsilon}}$ yields

$\displaystyle\mathfrak{R}_{1,{\varepsilon}}^{J^{\text{vol}}}=$
$\displaystyle{\varepsilon}^{d}\int_{0}^{1}\Bigg{\\{}\int_{\mathsf{D}_{\varepsilon}}[j^{({\varepsilon})}_{y_{1}}(u_{0}\circ
T_{\varepsilon}+s{\varepsilon}K_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon}+s\mathcal{D}K_{\varepsilon})$ $\displaystyle\hskip
142.26378pt-j^{({\varepsilon})}_{y_{1}}(u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})]({\varepsilon}K_{\varepsilon})\;dx$
$\displaystyle+{\varepsilon}^{d}\int_{\mathsf{D}_{\varepsilon}}[j^{({\varepsilon})}_{y_{2}}(u_{0}\circ
T_{\varepsilon}+s{\varepsilon}K_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon}+s\mathcal{D}K_{\varepsilon})$ $\displaystyle\hskip
142.26378pt-j^{({\varepsilon})}_{y_{2}}(u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})](\mathcal{D}K_{\varepsilon})\;dx\Bigg{\\}}\;ds.$

Letting
$\mathfrak{R}_{1}^{J^{\text{vol}}}(u_{0},p_{0}):=\lim_{{\varepsilon}\searrow
0}\frac{1}{|\omega_{\varepsilon}|}\mathfrak{R}_{1,{\varepsilon}}^{J^{\text{vol}}}$
and using ${\varepsilon}K_{\varepsilon}\to 0$ and
$\mathcal{D}K_{\varepsilon}\to\mathcal{D}K$, we obtain

$\displaystyle\mathfrak{R}_{1}^{J^{\text{vol}}}(u_{0},p_{0})$
$\displaystyle=\frac{1}{|\omega|}\int_{{\mathbf{R}}^{d}}[j^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z)+\mathcal{D}K)-j^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z))$
$\displaystyle\hskip
142.26378pt-j^{\omega}_{y_{2}}(u_{0}(z),\mathcal{D}u_{0}(z))(\mathcal{D}K)]\;dx.$

Here,
$j^{\omega}(y_{1},y_{2})=\chi_{\omega}\,j^{\text{in}}(y_{1},y_{2})+\chi_{{\mathbf{R}}^{d}\setminus\omega}\,j^{\text{out}}(y_{1},y_{2})$.
We proceed with the contribution of $J^{\text{vol}}$ to the term
$\mathfrak{R}_{2}(u_{0},p_{0})$. For this we compute

$\displaystyle\mathfrak{R}_{2,{\varepsilon}}^{J^{\text{vol}}}:=$
$\displaystyle[\partial_{u}J_{1}({\varepsilon},u_{0})-\partial_{u}J_{1}(0,u_{0})](u_{\varepsilon}-u_{0})$
$\displaystyle=$
$\displaystyle[\partial_{u}J^{\text{vol}}({\varepsilon},u_{0},\mathcal{D}u_{0})-\partial_{u}J^{\text{vol}}(0,u_{0},\mathcal{D}u_{0})](u_{\varepsilon}-u_{0})$
$\displaystyle=$
$\displaystyle\int_{\mathsf{D}}(j_{y_{1}}^{({\varepsilon})}(u_{0},\mathcal{D}u_{0})-j_{y_{1}}^{(0)}(u_{0},\mathcal{D}u_{0}))(u_{\varepsilon}-u_{0})$
$\displaystyle+\int_{\mathsf{D}}(j_{y_{2}}^{({\varepsilon})}(u_{0},\mathcal{D}u_{0})-j_{y_{2}}^{(0)}(u_{0},\mathcal{D}u_{0}))(\mathcal{D}(u_{\varepsilon}-u_{0})).$

Therefore setting
$\mathfrak{R}_{2}^{J^{\text{vol}}}(u_{0},p_{0}):=\lim_{{\varepsilon}\searrow
0}\frac{1}{|\omega_{\varepsilon}|}\mathfrak{R}_{2,{\varepsilon}}^{J^{\text{vol}}}$,
we obtain

$\displaystyle\mathfrak{R}_{2}^{J^{\text{vol}}}(u_{0},p_{0})=$
$\displaystyle\underset{{\varepsilon}\searrow
0}{\lim}\frac{1}{|\omega_{\varepsilon}|}[\partial_{u}J^{\text{vol}}({\varepsilon},u_{0},\mathcal{D}u_{0})-\partial_{u}J^{\text{vol}}(0,u_{0},\mathcal{D}u_{0})](u_{\varepsilon}-u_{0})$
$\displaystyle=$ $\displaystyle\underset{{\varepsilon}\searrow
0}{\lim}\frac{1}{|\omega_{\varepsilon}|}\Bigg{[}\int_{\mathsf{D}}(j_{y_{1}}^{({\varepsilon})}(u_{0},\mathcal{D}u_{0})-j_{y_{1}}^{(0)}(u_{0},\mathcal{D}u_{0}))(u_{\varepsilon}-u_{0})$
$\displaystyle\hskip
56.9055pt+(j_{y_{2}}^{({\varepsilon})}(u_{0},\mathcal{D}u_{0})-j_{y_{2}}^{(0)}(u_{0},\mathcal{D}u_{0}))(\mathcal{D}(u_{\varepsilon}-u_{0}))\;dx\Bigg{]}$
$\displaystyle=$ $\displaystyle\underset{{\varepsilon}\searrow
0}{\lim}\frac{1}{|\omega|}\Bigg{[}\int_{\omega}[j_{y_{1}}^{({\varepsilon})}(u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})-j_{y_{1}}^{(0)}(u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})]({\varepsilon}K_{\varepsilon})$ $\displaystyle\hskip
28.45274pt+[j_{y_{2}}^{({\varepsilon})}(u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})-j_{y_{2}}^{(0)}(u_{0}\circ
T_{\varepsilon},(\mathcal{D}u_{0})\circ
T_{\varepsilon})](\mathcal{D}K_{\varepsilon})\;dx\Bigg{]}$ $\displaystyle=$
$\displaystyle\frac{1}{|\omega|}\int_{\omega}[j_{y_{2}}^{\text{in}}(u_{0}(z),\mathcal{D}u_{0}(z))-j_{y_{2}}^{\text{out}}(u_{0}(z),\mathcal{D}u_{0}(z))](\mathcal{D}K)\;dx.$

Finally, the term of $J^{\text{vol}}$ contributing to
$\partial_{\ell}\mathfrak{L}$ reads

$\displaystyle\partial_{\ell}\mathfrak{L}^{J^{\text{vol}}}(0,u_{0},p_{0})=j^{\text{in}}(u_{0}(z),\mathcal{D}u_{0}(z))-j^{\text{out}}(u_{0}(z),\mathcal{D}u_{0}(z)).$
(2.22)

####  2.5.5. Terms coming from $J^{\text{bnd}}$

There is no contribution from the term $J^{\text{bnd}}$, since the rescaled
term $\partial\mathsf{D}_{\varepsilon}$ tends to ”infinity” and is therefore
not present.

####  2.5.6. Summary

In summary we showed that

$\displaystyle\mathfrak{R}_{1}(u_{0},p_{0})$
$\displaystyle=\mathfrak{R}_{1}^{A_{1}}(u_{0},p_{0})+\mathfrak{R}_{1}^{A_{2}}(u_{0},p_{0})+\mathfrak{R}_{1}^{J^{\text{vol}}}(u_{0},p_{0}),$
$\displaystyle\mathfrak{R}_{2}(u_{0},p_{0})$
$\displaystyle=\mathfrak{R}_{2}^{A_{1}}(u_{0},p_{0})+\mathfrak{R}_{2}^{A_{2}}(u_{0},p_{0})+\mathfrak{R}_{2}^{J^{\text{vol}}}(u_{0},p_{0}),$
$\displaystyle\partial_{\ell}\mathfrak{L}(u_{0},p_{0})$
$\displaystyle=\partial_{\ell}\mathfrak{L}^{A_{1}}(u_{0},p_{0})+\partial_{\ell}\mathfrak{L}^{A_{2}}(u_{0},p_{0})+\partial_{\ell}\mathfrak{L}^{F}(u_{0},p_{0})+\partial_{\ell}\mathfrak{L}^{J^{\text{vol}}}(u_{0},p_{0}).$

Therefore we proved the following theorem under the assumption that formally
$\nabla K_{\varepsilon}\to\nabla K$ in ${\mathbf{E}}({\mathbf{R}}^{d})$ and
${\varepsilon}K_{\varepsilon}\rightarrow 0$ and that $u_{0}$ is sufficiently
smooth at $z$.

######  Theorem 1.

The first order topological expansion of $\mathcal{J}$ at $\Omega)$ in the
points $z$ reads

$\displaystyle
d\mathcal{J}(\Omega,\omega)(z)=\mathfrak{R}_{1}(u_{0},p_{0})+\mathfrak{R}_{2}(u_{0},p_{0})+\partial_{\ell}\mathfrak{L}(0,u_{0},p_{0}),$
(2.23)

where

$\displaystyle\mathfrak{R}_{1}(u_{0},p_{0})=$
$\displaystyle\frac{1}{|\omega|}\int_{{\mathbf{R}}^{d}}[A_{1}^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z)+\mathcal{D}K)-A_{1}^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z))$
$\displaystyle\hskip
142.26378pt-\partial_{y_{2}}A_{1}^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z))(\mathcal{D}K)]\cdot
p_{0}(z)\;dx$
$\displaystyle+\frac{1}{|\omega|}\int_{{\mathbf{R}}^{d}}[A_{2}^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z)+\mathcal{D}K)-A_{2}^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z))$
$\displaystyle\hskip
142.26378pt-\partial_{y_{2}}A_{2}^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z))(\mathcal{D}K)]:\mathcal{D}p_{0}(z)\;dx$
$\displaystyle+\frac{1}{|\omega|}\int_{{\mathbf{R}}^{d}}[j^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z)+\mathcal{D}K(x))-j^{\omega}(u_{0}(z),\mathcal{D}u_{0}(z))$
$\displaystyle\hskip
170.71652pt-j^{\omega}_{y_{2}}(u_{0}(z),\mathcal{D}u_{0}(z))(\mathcal{D}K)]\;dx$
$\displaystyle\mathfrak{R}_{2}(u_{0},p_{0})=$
$\displaystyle\frac{1}{|\omega|}\int_{\omega}[\partial_{y_{2}}A_{1}^{\text{in}}(u_{0}(z),\mathcal{D}u_{0}(z))-\partial_{y_{2}}A_{1}^{\text{out}}(u_{0}(z),\mathcal{D}u_{0}(z))](\mathcal{D}K)\cdot
p_{0}(z)\;dx$
$\displaystyle+\frac{1}{|\omega|}\int_{\omega}[\partial_{y_{2}}A_{2}^{\text{in}}(u_{0}(z),\mathcal{D}u_{0}(z))-\partial_{y_{2}}A_{2}^{\text{out}}(u_{0}(z),\mathcal{D}u_{0}(z))](\mathcal{D}K):\mathcal{D}p_{0}(z)\;dx$
$\displaystyle+\frac{1}{|\omega|}\int_{\omega}[j_{y_{2}}^{\text{in}}(u_{0}(z),\mathcal{D}u_{0}(z))-j_{y_{2}}^{\text{out}}(u_{0}(z),\mathcal{D}u_{0}(z))](\mathcal{D}K)\;dx$
and $\displaystyle\partial_{\ell}\mathfrak{L}(0,u_{0},p_{0})=$
$\displaystyle[A_{1}^{\text{in}}(u_{0}(z),\mathcal{D}u_{0}(z))-A_{1}^{\text{out}}(u_{0}(z),\mathcal{D}u_{0}(z))]\cdot
p_{0}(z)$
$\displaystyle+[A_{2}^{\text{in}}(u_{0}(z),\mathcal{D}u_{0}(z))-A_{2}^{\text{out}}(u_{0}(z),\mathcal{D}u_{0}(z))]:\mathcal{D}p_{0}(z)$
$\displaystyle-(F_{1}^{\text{in}}-F_{1}^{\text{out}})\cdot
p_{0}(z)-(F_{2}^{\text{in}}-F_{2}^{\text{out}}):\mathcal{D}p_{0}(z)$
$\displaystyle+j^{\text{in}}(u_{0}(z),\mathcal{D}u_{0}(z))-j^{\text{out}}(u_{0}(z),\mathcal{D}u_{0}(z)).$

The function $K$ solves

$\displaystyle\begin{aligned}
\int_{{\mathbf{R}}^{d}}&[A_{2}^{\omega}(z,u_{0}(z),\mathcal{D}u_{0}(z)+\mathcal{D}K)-A_{2}^{\omega}(z,u_{0}(z),\mathcal{D}u_{0}(z))]:\mathcal{D}\psi\;dx\\\
&=\int_{\omega}(F_{2}^{\text{in}}-F_{2}^{\text{out}})(z):\mathcal{D}\psi\;dx-\int_{\omega}[A_{2}^{\text{in}}-A_{2}^{\text{out}}](z,u_{0}(z),\mathcal{D}u_{0}(z)):\mathcal{D}\psi\;dx\end{aligned}$
(2.24)

for all $\psi\in{\mathbf{E}}({\mathbf{R}}^{d})$.

##  3\. Numerical experiments

In this section, we consider several linear and nonlinear problems in two and
three space dimensions for which we numerically compute the topological
derivative by means of the procedure outlined in the previous section, see the
generic formulas of Theorem 1. We solve all underlying PDEs by means of
piecewise linear and globally continuous finite elements on triangular grids
using the finite element software package NGSolve [21]. We compute an
approximation of the corrector function $K$, which is defined as the solution
of a PDE on the unbounded domain ${\mathbf{R}}^{d}$, by solving (2.24) on a
large, but bounded domain $B_{R}:=B(\mathbf{0},R)$ with homogeneous Dirichlet
boundary conditions on $\partial B_{R}$. This procedure is motivated by the
fact that the solution $K$ often can be shown to exhibit a decay behaviour as
$|x|\rightarrow\infty$. We will show that the numerically computed topological
derivative matches the analytical formula well for the topological derivative
in cases where this formula is known. Moreover, we perform numerical tests to
verify the topological asymptotic expansion

$\displaystyle\mathcal{J}(\Omega_{\varepsilon})=\mathcal{J}(\Omega)+{\varepsilon}^{d}|\omega|d\mathcal{J}(\Omega)(z)+\mathcal{O}({\varepsilon}^{d+1}),$
(3.1)

or, in other words

$\displaystyle\delta\mathcal{J}:=|\mathcal{J}(\Omega_{\varepsilon})-\mathcal{J}(\Omega)+{\varepsilon}^{d}|\omega|d\mathcal{J}(\Omega)(z)|=\mathcal{O}({\varepsilon}^{d+1}).$
(3.2)

We consider several different inclusion shapes in two and three space
dimensions. In 2D, we consider

  1. (1)

$\omega=\omega_{2D}^{(1)}:=B((0,0),1)$ the unit disk in ${\mathbf{R}}^{2}$

  2. (2)

$\omega=\omega_{2D}^{(2)}:=B((0.5,0.5),1)$ a shifted unit disk in
${\mathbf{R}}^{2}$

  3. (3)

$\omega=\omega_{2D}^{(3)}:=Ell((0,0),\frac{3}{2},\frac{2}{3})$ an axis-aligned
ellipse centered at the origin with axes lengths $\frac{3}{2}$ and
$\frac{2}{3}$,

  4. (4)

$\omega=\omega_{2D}^{(4)}:=Ell((0.5,0.5),\frac{3}{2},\frac{2}{3})$ a shifted
axis-aligned ellipse centered at the point $(0.5,0.5)^{\top}$ with the same
axes lengths $\frac{3}{2}$ and $\frac{2}{3}$,

  5. (5)

$\omega=\omega_{2D}^{(5)}:=$ an L-shaped domain of area $\pi$,

and in 3D the shapes

  1. (1)

$\omega=\omega_{3D}^{(1)}:=B((0,0,0),1)$ the unit ball in ${\mathbf{R}}^{3}$

  2. (2)

$\omega=\omega_{3D}^{(2)}:=B((0.5,0.5,0.5),1)$ a shifted unit ball in
${\mathbf{R}}^{3}$

  3. (3)

$\omega=\omega_{3D}^{(3)}:=Ell((0,0,0),\frac{3}{2},\frac{2}{3},1)$ an axis-
aligned ellipsoid centered at the origin with axes lengths $\frac{3}{2}$,
$\frac{2}{3}$ and $1$,

  4. (4)

$\omega=\omega_{3D}^{(4)}:=Ell((0.5,0.5,0.5),\frac{3}{2},\frac{2}{3},1)$ a
shifted axis-aligned ellipsoid centered at the point $(0.5,0.5,0.5)^{\top}$
with the same axes lengths $\frac{3}{2}$, $\frac{2}{3}$ and $1$.

Note that $\omega_{2D}^{(1)}$, $\omega_{2D}^{(3)}$, $\omega_{3D}^{(1)}$ and
$\omega_{3D}^{(3)}$ are symmetric with respect to the $x_{1}$-, the $x_{2}$\-
(and the $x_{3}$-)axes while the other shapes are not.

The finite element software NGSolve allows to define PDEs in weak form in a
symbolic way and also supports automated differentiation of expressions, see
[16] for applications of these capabilities in the context of shape
derivatives. Our implementation is available from [15]. It consists of a main
file, which is completely independent of the concrete topology optimization
problem at hand, and four other files defining the geometry, the PDE, the cost
function as well as some algorithmic parameters. The main file implements the
topological derivative of a general problem of the form (2.8), i.e., it
implements the solution of the state and adjoint equation, of the
corresponding corrector equation (2.24) on the large domain $B_{R}$ as well as
the computation of the terms $\mathfrak{R}_{1}(u_{0},p_{0})$,
$\mathfrak{R}_{2}(u_{0},p_{0})$ and
$\partial_{\ell}\mathfrak{L}(0,u_{0},p_{0})$ as they are given in Section
2.5.6.

###  3.1. Diffusion, convection, reaction

Here, we consider the class of topology optimization problems with tracking-
type cost functionals and a scalar diffusion-convection-reaction equation as a
PDE constraint. Given the computational domain whose boundary is divided into
the disjoint Dirichlet and Neumann boundaries,
$\partial\mathsf{D}=\Gamma_{D}\cup\Gamma_{N}$, the PDE-constrained topology
optimization problem is to find $u\in H^{1}_{\Gamma_{D}}(\mathsf{D}):=\\{v\in
H^{1}(\mathsf{D}):v|_{\Gamma_{D}}=0\\}$ and $\Omega\in\mathcal{A}$ for some
set of admissible shapes $\mathcal{A}$ as a solution to

$\displaystyle\underset{\Omega}{\mbox{min
}}J(u,\Omega):=\int_{\mathsf{D}}\tilde{\alpha}_{\Omega}(x)|u-u_{d}|^{2}+\tilde{\beta}_{\Omega}(x)|\nabla(u-u_{d})|^{2}\;\mbox{d}x+$
$\displaystyle\tilde{\gamma}\int_{\Gamma_{N}}|u-u_{d}|^{2}\;\mbox{d}s$ (3.3)

such that

$\displaystyle\int_{\mathsf{D}}\beta_{\Omega}(x,|\nabla u|)\nabla u$
$\displaystyle\cdot\nabla\psi+(\mathbf{b}_{\Omega}(x)\cdot\nabla
u)\psi+\alpha_{\Omega}(x,u)\psi\;\mbox{d}x$ (3.4) $\displaystyle=$
$\displaystyle\int_{\mathsf{D}}f_{\Omega}(x)\psi+\mathbf{M}_{\Omega}(x)\cdot\nabla\psi\;\mbox{d}x+\int_{\Gamma_{N}}g_{N}\cdot\psi\;\mbox{d}S$
(3.5)

for all $\psi\in H^{1}_{\Gamma_{D}}(\mathsf{D})$. Here,
$c_{\Omega}(x):=\chi_{\Omega}(x)c_{1}+\chi_{\mathsf{D}\setminus\Omega}(x)c_{2}$
for some given constants $c_{1}$, $c_{2}$ for
$c\in\\{\tilde{\alpha},\tilde{\beta},\mathbf{b},f,\mathbf{M}\\}$. The
functions $\beta_{\Omega}(x,|\nabla u|)$ and $\alpha_{\Omega}(x,u)$ will be
defined piecewise in the subsequent subsections.

This problem fits into the framework considered in Section 2 with the choices

$\displaystyle j^{\Omega}(x,u,\nabla u)=$
$\displaystyle\tilde{\alpha}_{\Omega}(x)|u-u_{d}|^{2}+\tilde{\beta}_{\Omega}(x)|\nabla(u-u_{d})|^{2},$
$\displaystyle j^{bnd}(x,u,\nabla u)=$
$\displaystyle\tilde{\gamma}|u-u_{d}|^{2},$ $\displaystyle
A_{1}^{\Omega}(x,u,\nabla u)=$
$\displaystyle(\mathbf{b}_{\Omega}(x)\cdot\nabla u)+\alpha_{\Omega}(x,u),$
$\displaystyle A_{2}^{\Omega}(x,u,\nabla u)=$
$\displaystyle\beta_{\Omega}(x,|\nabla u|)\nabla u$ $\displaystyle
F_{1}^{\Omega}(x)=$ $\displaystyle f_{\Omega}(x),$ $\displaystyle
F_{2}^{\Omega}(x)=$ $\displaystyle\mathbf{M}_{\Omega}(x).$

####  3.1.1. Example 1: A linear diffusion-convection-reaction problem in 2D

We begin with a simple two-dimensional, linear version of problem (3.3)–(3.5)
where we set $u_{d}=0$, $\tilde{\alpha}_{1}=1$, $\tilde{\alpha}_{2}=2$,
$\tilde{\beta}_{1}=\tilde{\beta}_{2}=0$, $\tilde{\gamma}=0$,
$\mathbf{b}_{1}=(1,0)^{\top}$, $\mathbf{b}_{2}=(0,1)^{\top}$, $f_{1}=1$,
$f_{2}=2$, $\mathbf{M}_{1}=\mathbf{M}_{2}=(0,0)^{\top}$,
$g_{N}(x_{1},x_{2})=x_{1}x_{2}$. Moreover, we choose the piecewise constant
functions

$\displaystyle\beta_{\varepsilon}(x,|\nabla u|)=$
$\displaystyle\chi_{\Omega_{\varepsilon}}(x)\beta_{1}+\chi_{\mathsf{D}\setminus\Omega_{\varepsilon}}(x)\beta_{2},$
$\displaystyle\alpha_{\varepsilon}(x,u)=$
$\displaystyle\chi_{\Omega_{\varepsilon}}(x)\alpha_{1}+\chi_{\mathsf{D}\setminus\Omega_{\varepsilon}}(x)\alpha_{2},$

with the values $\beta_{1}=1$, $\beta_{2}=2$, $\alpha_{1}=1$, $\alpha_{2}=2$,
which thus do not depend on $u$ or $\nabla u$.

In this setting, it is well-known that, for a point
$z\in\mathsf{D}\setminus\Omega$ and inclusion shape $\omega=B_{1}(0)$ the unit
disk, the corrector term $K$ is linear inside $\omega$ and satisfies

$\displaystyle\nabla
K|_{\omega}=-\frac{\beta_{1}-\beta_{2}}{\beta_{1}+\beta_{2}}\nabla u(z),$
(3.6)

see e.g. [13]. Thus, it follows that

$\displaystyle\mathfrak{R}_{2}^{A_{2}}(u,p)=$
$\displaystyle\frac{1}{|\omega|}\int_{\omega}(\beta_{1}-\beta_{2})\nabla
K\cdot\nabla
p(z)\;dx=-(\beta_{1}-\beta_{2})\frac{\beta_{1}-\beta_{2}}{\beta_{1}+\beta_{2}}\nabla
u(z)\cdot\nabla p(z)$ (3.7)

which together with the term
$\partial_{\ell}\mathfrak{L}^{A_{2}}(u,p)=(\beta_{1}-\beta_{2})\nabla
u(z)\cdot\nabla p(z)$ sums up to

$\displaystyle\mathfrak{R}_{2}^{A_{2}}(u,p)+\partial_{\ell}\mathfrak{L}^{A_{2}}(u,p)=$
$\displaystyle 2\beta_{2}\frac{\beta_{1}-\beta_{2}}{\beta_{1}+\beta_{2}}\nabla
u(z)\cdot\nabla p(z).$ (3.8)

Similarly,
$\partial_{\ell}\mathfrak{L}^{A_{1}}(u,p)=(\mathbf{b}_{1}-\mathbf{b}_{2})\cdot\nabla
u(z)\,p(z)+(\alpha_{1}-\alpha_{2})u(z)p(z)$ together with

$\displaystyle\mathfrak{R}_{2}^{A_{1}}(u,p)=$
$\displaystyle\frac{1}{|\omega|}\int_{\omega}(\mathbf{b}_{1}-\mathbf{b}_{2})\cdot\nabla
K\,p(z)\;dx=-\frac{\beta_{1}-\beta_{2}}{\beta_{1}+\beta_{2}}(\mathbf{b}_{1}-\mathbf{b}_{2})\cdot\nabla
u(z)\,p(z)$ (3.9)

adds up to

$\displaystyle\mathfrak{R}_{2}^{A_{1}}(u,p)+\partial_{\ell}\mathfrak{L}^{A_{1}}(u,p)=2\frac{\beta_{2}}{\beta_{1}+\beta_{2}}(\mathbf{b}_{1}-\mathbf{b}_{2})\cdot\nabla
u(z)\,p(z)+(\alpha_{1}-\alpha_{2})u(z)p(z).$ (3.10)

Moreover, we have $\mathfrak{R}_{1}(u_{0},p_{0})=0$ since the PDE constraint
is linear and the cost function does not depend on $\nabla u$. Together with
the remaining terms $\partial_{\ell}\mathfrak{L}^{F}$ and
$\partial_{\ell}\mathfrak{L}^{J^{\text{vol}}}$, the topological derivative
reads in closed form,

$\displaystyle\begin{aligned}
d\mathcal{J}(\Omega,\omega)(z)=&2\beta_{2}\frac{\beta_{1}-\beta_{2}}{\beta_{1}+\beta_{2}}\nabla
u(z)\cdot\nabla
p(z)+2\frac{\beta_{2}}{\beta_{1}+\beta_{2}}(\mathbf{b}_{1}-\mathbf{b}_{2})\cdot\nabla
u(z)\,p(z)+(\alpha_{1}-\alpha_{2})u(z)p(z)\\\
&-(f_{1}-f_{2})p(z)+(\tilde{\alpha}_{1}-\tilde{\alpha}_{2})u(z)^{2},\end{aligned}$
(3.11)

see also, e.g., [3] for the case where
$\mathbf{b}_{1}=\mathbf{b}_{2}=(0,0)^{\top}$. A similar problem including a
convection term is covered by the analysis in [7], however, there the factor
$2\beta_{2}/(\beta_{1}+\beta_{2})$ (and thus the term
$\mathfrak{R}_{2}^{A_{1}}(u,p)$) is missing. We remark that our numerical
experiments indicate that this factor is important for having the correct
topological derivative formula.

We numerically compute the topological derivative by the procedure explained
in Section 2 for $\mathsf{D}=(-1,1)^{2}$
$\Gamma_{D}=\\{(x_{1},x_{2}):x_{1}=-1\mbox{ or }x_{2}=-1\\}$,
$\Gamma_{N}=\partial\mathsf{D}\setminus\Gamma_{D}$ and
$\Omega=B((0,-0.5)^{\top},0.3)$ at the point
$z=(0,0.5)^{\top}\in\mathsf{D}\setminus\Omega$ for the five inclusion shapes
$\omega_{2D}^{(i)}$, $i=1,\dots,5$ defined above. Figure 1 depicts the
computational domain $\mathsf{D}$ as well as the unperturbed state obtained on
a mesh with 528724 vertices.

|  |
---|---|---
(a) | (b) | (c)
Figure 1. (a) Computational domain $\mathsf{D}$ with subdomain $\Omega$ for
Examples 1 and 2. (b) Solution $u$ to unperturbed state equation for linear
Example 1. (c) Solution $u$ to unperturbed state equation for nonlinear
Example 2.

We solved the corresponding problem (2.24) for the corrector $K$ on the
bounded domain $B(\mathbf{0},R)$ with $R=1000$ using a finite element mesh
consisting of 371950 vertices. In order to numerically verify the computed
values, we performed a Taylor test for the topological asymptotic expansion
(3.1). We chose a set of decreasing values for the inclusion radius
${\varepsilon}$, namely

$\underline{{\varepsilon}}:=\\{{\varepsilon}_{0}\,\delta^{9},\dots,{\varepsilon}_{0}\,\delta^{1},{\varepsilon}_{0}\,\delta^{0}\\}$
(3.12)

with two constants ${\varepsilon}_{0}=0.005$ and $\delta=1.5$, and solved the
perturbed PDE constraint (2.9) on a mesh which is highly refined around the
point $z$ and evaluated the cost function to get perturbed cost function
values $\mathcal{J}(\Omega_{\varepsilon})$. For all
${\varepsilon}\in\underline{{\varepsilon}}$ and for each of the five inclusion
shapes $\omega_{2D}^{(i)}$, we computed the quantities $\delta\mathcal{J}$
defined in (3.2). In Figure 2, it can be observed that, for each inclusion
shape, $\delta\mathcal{J}$ behaves (at least) like ${\varepsilon}^{3}$. In the
case of the two inclusions $\omega_{2D}^{(1)}$, $\omega_{2D}^{(3)}$ which are
symmetric with respect to both the $x_{1}$\- and the $x_{2}$-axes, we can even
observe that $\delta\mathcal{J}=\mathcal{O}({\varepsilon}^{4})$. This is
expected since it is known that the second term in the expansion (3.1)
vanishes, $d^{2}\mathcal{J}(\Omega,\omega)(z)=0$ for this example if symmetric
inclusion shapes $\omega$ are considered, see e.g. [8]. In Figure 2 as well as
in the subsequent Taylor test graphs, the data for $\omega_{2D}^{(i)}$, $i>1$,
is scaled such that they coincide with the data for $\omega_{2D}^{(1)}$ for
the largest considered value of ${\varepsilon}$. This is done for better
comparison of the convergence rates.

Figure 2. Numerical verification of topological asymptotic expansion (3.1) for
Example 1 for five different inclusion shapes. Figure 3. Numerical
verification of topological asymptotic expansion (3.1) for Example 2 for five
different inclusion shapes.

We will further investigate the numerically computed topological derivative
for this example in Section 4 where we will also discuss the efficient
evaluation of the topological derivative in the full computational domain
$\mathsf{D}$.

####  3.1.2. Example 2: A nonlinear diffusion-convection-reaction problem in
2D

Next, we consider a more general, quasilinear setting of problem (3.3)–(3.5)
where we set $u_{d}=0$, $\tilde{\alpha}_{1}=1$, $\tilde{\alpha}_{2}=2$,
$\tilde{\beta}_{1}=1$, $\tilde{\beta}_{2}=2$, $\tilde{\gamma}=1$,
$\mathbf{b}_{1}=(1,0)^{\top}$, $\mathbf{b}_{2}=(0,1)^{\top}$, $f_{1}=1$,
$f_{2}=2$, $\mathbf{M}_{1}=(1,0)^{\top}$, $\mathbf{M}_{2}=(0,1)^{\top}$,
$g_{N}(x_{1},x_{2})=x_{1}x_{2}$ as well as the nonlinear functions

$\displaystyle\beta_{\varepsilon}(x,|\nabla
u(x)|)=\begin{cases}1&x\in\Omega_{\varepsilon},\\\ \hat{\beta}_{2}(|\nabla
u(x)|)&\mbox{else},\end{cases}\qquad\alpha_{\varepsilon}(x,u(x))=\begin{cases}1&x\in\Omega_{\varepsilon},\\\
u(x)^{3}&\mbox{else},\end{cases}$

where $\hat{\beta}_{2}(s)=\nu_{0}-(\nu_{0}-200)e^{-s^{6}/1000}$ with
$\nu_{0}=10^{7}/(4\pi)$. This function is sometimes used as a magentic
reluctivity function in electromagnetics and satisfies the monotonicity and
Lipschitz conditions which ensure existence of a unique solution to the PDE
constraint. Topological derivatives for problems involving quasilinear PDE
constraints are challenging from both the analytical and the numerical point
of view, see e.g. [5, 6, 12, 14]. The numerical experiments depicted in Figure
3 seem to exhibit the behavior $\mathcal{O}({\varepsilon}^{3})$ for all five
inclusion shapes which is in accordance with (3.2) and therefore confirms the
numerically computed value for the topological derivative.

####  3.1.3. Example 3: A linear diffusion-convection-reaction problem in 3D

Next, we consider a linear three-dimensional version of (3.3)–(3.5) on the
domain $\mathsf{D}=B(\mathbf{0},1)$ with the subdomain
$\Omega=B((0,-0.7,0)^{\top},0.2)$. We chose the parameters $u_{d}=0$,
$\tilde{\alpha}_{1}=1$, $\tilde{\alpha}_{2}=2$,
$\tilde{\beta}_{1}=\tilde{\beta}_{2}=0$, $\tilde{\gamma}=0$, $\beta_{1}=1$,
$\beta_{2}=2$, $\mathbf{b}_{1}=\mathbf{b}_{2}=(0,0)^{\top}$, $f_{1}=1$,
$f_{2}=2$, $\mathbf{M}_{1}=\mathbf{M}_{2}=(0,0)^{\top}$,
$g_{N}(x_{1},x_{2})=x_{1}x_{2}$ and evaluate the topological derivative at the
point $z=(0,0.1,0)^{\top}$. For each of the four inclusion shapes
$\omega_{3D}^{(i)}$, $i=1,2,3,4$, we solved the corresponding corrector
equation (2.24) on a three-dimensional ball $B_{R}$ of radius $R=1000$ using a
tetrahedral mesh with 301116 vertices. Subsequently, we computed the
topological derivative according to Section 2.5.6 and computed the quantities
$\delta\mathcal{J}$ (3.2) for the vector of radius values (3.12) with
${\varepsilon}_{0}=0.05$ and $\delta=1.25$. Figure 4 shows that the quantities
$\delta\mathcal{J}$ decay at least as fast as ${\varepsilon}^{4}$ for all four
inclusion shapes, which is in accordance with (3.1). The deterioration of the
rates in Figure 4 can be attributed to the discretization error. We used a
mesh consisting of 245177 vertices for the computational domain $\mathsf{D}$
which is highly refined around the point $z$.

Figure 4. Numerical verification of topological asymptotic expansion (3.1) for
Example 3 for four different inclusion shapes. Figure 5. Numerical
verification of topological asymptotic expansion (3.1) for Example 4 for four
different inclusion shapes.

###  3.2. Elasticity

In this section, we consider vector-valued partial differential equations
coming from elasticity as PDE constraints. We consider the computational
domain $\mathsf{D}=(0,2)\times(0,1)$ which is clamped at the top left and
bottom left, $\Gamma_{D}=\\{(x_{1},x_{2}):x_{1}=0\mbox{ and
}x_{2}\in[0,0.12]\cup[0.88,1]\\}$ and is subject to a downward directed force
at $\Gamma_{N}={1}\times(0.45,0.55)$. On the rest of the boundary,
$\Gamma:=\partial\mathsf{D}\setminus(\Gamma_{D}\cup\Gamma_{N})$, homogeneous
Neumann boundary conditions are set. The domain consists of a strong material
with Young’s modulus $E=1000$, and a very weak material which mimicks void
inside $\Omega=B((0.5,0.5)^{\top},0.3)\cup B((1.8,0.25)^{\top},0.1)\cup
B((1.8,0.75)^{\top},0.1)$ where $E=0.1$, see Figure 6(a). The Poisson ratio is
set to $\nu=\frac{1}{3}$ in both subdomains. We consider the problem to
minimize the compliance of the structure which is subject to self-weight as
well as an external force acting on the boundary $\Gamma_{N}$. The problem
reads

$\displaystyle\underset{\Omega}{\mbox{min
}}J(u,\Omega):=\frac{1}{2}\int_{\mathsf{D}}S_{\Omega}(\mathcal{D}u):\mathcal{D}u\;\mbox{d}x$
(3.13)

such that

$\displaystyle u\in
H^{1}_{\Gamma_{D}}(\mathsf{D})^{d}:\int_{\mathsf{D}}S_{\Omega}(\mathcal{D}u):\mathcal{D}\psi\;\mbox{d}x=$
$\displaystyle\int_{\mathsf{D}}f_{\Omega}\cdot\psi\;\mbox{d}x+\int_{\Gamma_{N}}g_{N}\cdot\psi\;\mbox{d}S$
(3.14)

for all $\psi\in H^{1}_{\Gamma_{D}}(\mathsf{D})^{d}$, where
$S_{\Omega}:{\mathbf{R}}^{d\times d}\rightarrow{\mathbf{R}}^{d\times d}$
represents a stress tensor and $f_{\Omega}$ is a piecewise constant vector
function. Also this problem fits into the framework considered in Section 2 by
setting

$\displaystyle j^{\Omega}(x,u(x),\mathcal{D}u(x))=$
$\displaystyle\frac{1}{2}S_{\Omega}(\mathcal{D}u(x)):\mathcal{D}u(x),$
$\displaystyle A_{1}^{\Omega}(x,u(x),\mathcal{D}u(x))=$ $\displaystyle 0,$
$\displaystyle A_{2}^{\Omega}(x,u(x),\mathcal{D}u(x))=$ $\displaystyle
S_{\Omega}(\mathcal{D}u(x)),$ $\displaystyle F_{1}^{\Omega}(x)=$
$\displaystyle f_{\Omega}(x),$ $\displaystyle F_{2}^{\Omega}(x)=$
$\displaystyle 0.$

|
---|---
(a) | (b)
|
(c) | (d)
Figure 6. (a) Computational domain $\mathsf{D}$ with subdomain $\Omega$ for
Examples 4 and 5. (b) Solution to linear state equation of Example 4
(undeformed). (c) Solution to nonlinear state equation of Example 5 in
deformed configuration. (d) Deformed domain of Example 5.

####  3.2.1. Example 4: Linear elasticity in 2D

We consider problem (3.13)–(3.14) for the two-dimensional cantilever example
introduced above under a linear stress-strain relation,

$\displaystyle
S_{\Omega}(\mathcal{D}u)=2\mu_{\Omega}e(u)+\lambda_{\Omega}Tr(e(u))I_{2}$
(3.15)

with the piecewise defined Lamé coefficients

$\displaystyle\mu_{\Omega}(x)=\chi_{\Omega}(x)\mu_{1}+\chi_{\mathsf{D}\setminus\Omega}(x)\mu_{2},\qquad\lambda_{\Omega}(x)=\chi_{\Omega}(x)\lambda_{1}+\chi_{\mathsf{D}\setminus\Omega}(x)\lambda_{2},$

corresponding to Young’s modulus $E_{2}=1000$, $E_{1}=0.1$ and Poisson ratio
$\nu_{2}=\nu_{1}=\frac{1}{3}$, i.e.,

$\displaystyle\mu_{i}=E_{i}\frac{1}{2(1+\nu_{i})}\qquad\lambda_{i}=E_{i}\frac{\nu_{i}}{(1+\nu_{i})(1-2\nu_{i})},\quad\mbox{for
}i=1,2.$ (3.16)

Here, $e(u)=\frac{1}{2}(\mathcal{D}u+\mathcal{D}u^{\top})$ denotes the
linearized strain tensor, $Tr(A)$ denotes the trace of a matrix $A$ and
$I_{2}$ denotes the two-dimensional identity matrix. Moreover, in this example
we neglect self-weight, $f_{\Omega}=0$, and consider the external load
$g_{N}=(0,-1)^{\top}$ acting on $\Gamma_{N}$.

We computed the unperturbed state and adjoint on a mesh with about 150 000
vertices (resulting in about 300 000 degrees of freedom) which is highly
refined around the point $z=(1.2,0.5)^{\top}$ for which we compute and verify
the topological derivative. We solve the corrector equation (2.24) with four
considered two-dimensional inclusion shapes $\omega_{2D}^{(i)}$, $i=1,\dots 4$
where we replace the unbounded domain ${\mathbf{R}}^{d}$ by the large, but
bounded domain $B_{R}=B(\mathbf{0},R)$ with $R=1000$ using a mesh with about
370 000 vertices. We again compute the quantities $\delta\mathcal{J}$ defined
in (3.2) for a range of values for the inclusion radius ${\varepsilon}$ (3.12)
with ${\varepsilon}_{0}=0.005$ and $\delta=1.5$. The behavior of
$\delta\mathcal{J}$ for these five inclusion shapes can be seen in Figure 5
where we can again observe the behavior
$\delta\mathcal{J}=\mathcal{O}({\varepsilon}^{3})$ for all five inclusion
shapes.

####  3.2.2. Example 5: Nonlinear elasticity in 2D

Finally, we consider problem (3.13)–(3.14) with the nonlinear St. Venant-
Kirchhoff material

$\displaystyle
S_{\Omega}(\mathcal{D}u):=(I_{2}+\mathcal{D}u)\left[\lambda_{\Omega}\mbox{Tr}(\frac{1}{2}(C(\mathcal{D}u)-I_{2}))I_{2}+\mu_{\Omega}(C(\mathcal{D}u)-I_{2})\right],$
(3.17)

with
$C(\mathcal{D}u):=(I_{2}+\mathcal{D}u)^{\top}(I_{2}+\mathcal{D}u)=I_{2}+\mathcal{D}u+\mathcal{D}u^{\top}+\mathcal{D}u^{\top}\mathcal{D}u$,
and
$\lambda_{\Omega}=\chi_{\Omega}\lambda_{1}+\chi_{\mathsf{D}\setminus\Omega}\lambda_{2}$,
$\mu_{\Omega}=\chi_{\Omega}\mu_{1}+\chi_{\mathsf{D}\setminus\Omega}\mu_{2}$
with $E_{2}=1000$, $E_{1}=0.1$ and $\nu_{2}=\nu_{1}=0.3$ and
$\lambda_{1},\lambda_{2},\mu_{1},\mu_{2}$ defined as in (3.16). Here we use
the larger external force $g_{N}=(0,-20)^{\top}$ acting on $\Gamma_{N}$ and
also mimic self-weight with $f_{2}=(0,-5)^{\top}$, $f_{1}=0$.

For solving the state, adjoint and corrector equations, we use the same meshes
as in Example 4. For solving the nonlinear state equation (2.4) we use a load
stepping scheme with $N_{ls}=20$ load steps where we gradually increase the
loads given by $f_{2}$ and $g_{N}$ until they reach their full values given
above, i.e., we solve (2.4) with the data given above and

$\displaystyle f_{2}=\frac{k}{N_{ls}}\left(\begin{array}[]{c}0\\\
-5\end{array}\right)\qquad g_{N}=\frac{k}{N_{ls}}\left(\begin{array}[]{c}0\\\
-20\end{array}\right)$

for $k=1,\dots N_{ls}$. For each load step, we use Newton’s method to solve
(2.4) for the given loads. The solution to the nonlinear state equation in
deformed configuration can be seen in Figure 6(c) and the deformed domain in
Figure 6(d). Moreover, we solve the corrector equation (2.24) by means of a
damped Newton method with a very conservative damping factor of $0.002$.

We performed a Taylor test by computing the quantities $\delta\mathcal{J}$ for
the four inclusion shapes $\omega_{2D}^{(i)}$, $i=1,2,3,4$, using the same
parameters $\underline{{\varepsilon}}$ as in Example 4. Figure 7 shows that
$\delta\mathcal{J}$ decays at least as fast as ${\varepsilon}^{3}$ which
confirms the topologcial asymptotic expansion (3.1).

Figure 7. Numerical verification of topological asymptotic expansion (3.1) for
Example 5 for four different inclusion shapes.

##  4\. On the efficient evaluation of topological derivatives in the full
domain

In Section 3, we applied the systematic procedure for the computation of
topological derivatives introduced in Section 2 to different model problems
involving linear and nonlinear PDE constraints in two and three space
dimensions. Moreover we verified the obtained values by Taylor tests. However,
looking at the formulas of Section 2.5.6, it can be seen that the evaluation
of the topological derivative at a spatial point $z$ requires the solution of
the corrector equation (2.24). On the other hand, in order to employ an
optimization algorithm, one usually is interested in the topological
derivative in the full design domain. Of course, solving (2.24) for every
spatial point in the design domain (or every element of a mesh) is
prohibitively expensive.

In the following, we show how, for linear or semilinear PDE constraints, the
corresponding solution $K$ can be obtained as a linear combination of some
functions which can be precomputed. Note that, since we assume the PDE
constraint to be linear or semilinear, both $A_{2}^{\text{in}}$ and
$A_{2}^{out}$ are linear with respect to the third argument $\mathcal{D}u$ and
constant with respect to the second argument $u$, i.e., there exist
$a_{2}^{\text{in}}(z),a_{2}^{out}(z)$ such that

$\displaystyle
A_{2}^{\text{in}}(z,u_{0}(z),\mathcal{D}u_{0}(z))=a_{2}^{\text{in}}(z)\mathcal{D}u_{0}(z)\quad\mbox{and}\quad
A_{2}^{out}(z,u_{0}(z),\mathcal{D}u_{0}(z))=a_{2}^{out}(z)\mathcal{D}u_{0}(z).$

We introduce the notation

$\displaystyle
a_{2}^{\omega}(z,\mathcal{D}u_{0}(z)):=\chi_{\omega}(x)a_{2}^{\text{in}}(z)\mathcal{D}u_{0}(z)+\chi_{{\mathbf{R}}^{d}\setminus\omega}(x)a_{2}^{out}(z)\mathcal{D}u_{0}(z).$

Now let $\hat{K}$ be the solution to

$\displaystyle\int_{{\mathbf{R}}^{d}}$ $\displaystyle
a_{2}^{\omega}(z,\mathcal{D}\hat{K}):\mathcal{D}\psi\;dx=\int_{\omega}(F_{2}^{\text{in}}-F_{2}^{out})(z):\mathcal{D}\psi\;dx$
(4.1)

for all $\psi$ and, for $i\in\\{1,\dots,m\\}$, $j\in\\{1,\dots,d\\}$, let
$\tilde{K}_{e_{ij}}$ be the solution to

$\displaystyle\begin{aligned}
\int_{{\mathbf{R}}^{d}}&a_{2}^{\omega}(z,\mathcal{D}\tilde{K}_{e_{ij}}):\mathcal{D}\psi\;dx&=-\int_{\omega}(a_{2}^{\text{in}}-a_{2}^{out})(z)e_{ij}:\mathcal{D}\psi\;dx\end{aligned}$
(4.2)

for all $\psi$. Here $e_{ij}\in{\mathbf{R}}^{m\times d}$ denotes the unit
basis element with the value $1$ at the position $(i,j)$ and the value $0$
else.

These considerations immediately yield the following lemma.

######  Lemma 1.

Assume that the PDE constraint (2.4) is linear or semilinear and let
$\mathcal{D}u_{0}(z)\in{\mathbf{R}}^{m\times d}$ be given. Then the solution
$K$ to (2.24) is given as the linear combination

$\displaystyle
K=\hat{K}+\sum_{i=1}^{m}\sum_{j=1}^{d}(a_{2}^{\text{in}}-a_{2}^{out})(z)\mathcal{D}u_{0}(z)[i,j]\,\tilde{K}_{e_{ij}}$
(4.3)

with $\hat{K}$ and $\tilde{K}_{e_{ij}}$ given in (4.1) and (4.2),
respectively.

Using Lemma 1 it is possible to compute the solution $K$ of (2.24) for any
value of $\mathcal{D}u_{0}(z)$ by linear combination once the $md+1$ functions
$\hat{K}$, $\tilde{K}_{e_{ij}}$ have been computed. By integration over
$\omega$ or ${\mathbf{R}}^{d}$ (which is numerically approximated by
$B(\mathbf{0},R)$ with $R=1000$), the topological derivative can then be
approximately evaluated at every point without having to solve an additional
boundary value problem.

In the case of Example 1 (Section 3.1.1), it follows from
$\mathbf{M}_{1}=\mathbf{M}_{2}=(0,0)^{\top}$ that $\hat{K}=0$. Since $m=1$,
$d=2$, we need to precompute the solutions to (4.2) for the two unit vectors
in ${\mathbf{R}}^{2}$ and can obtain the solution $K$ to (2.24) by Lemma 1.
The solutions $\tilde{K}_{e_{11}}$ and $\tilde{K}_{e_{12}}$ to (4.2) with
${\mathbf{R}}^{d}$ replaced by $B(\mathbf{0},1000)$ are depicted in Figure 8.

|
---|---
Figure 8. Solutions to (4.2) for $e_{11}=(1,0)^{\top}$ (left) and
$e_{12}=(0,1)^{\top}$ (right).

Next, we chose a (comparably coarse) mesh with only background material, i.e.
$\Omega=\emptyset$, and numerically computed the topological derivative for
the centroid of every triangular element by evaluating the formulas of Section
2.5.6. Also here, integrals over ${\mathbf{R}}^{d}$ were approximated by
integrals over the large, but bounded ball $B(\mathbf{0},1000)$. A comparison
of the obtained results with the analytical formula for the topological
derivative given in (3.11) showed good accordance, see Figure 9.

|  |
---|---|---
Figure 9. Left: Numerically computed topological derivative for Example 1.
Center: Analytical formula (3.11). Right: Difference between numerical and
analytical formulas. Observe the different orders of magnitude.

######  Remark 4.

The procedure described in this section is still computationally expensive
since it involves the numerical computation of several integrals over $\omega$
and $B(\mathbf{0},1000)$ for each topological derivative evaluation.

We remark that in many cases it is possible to exploit the (affine) linearity
of the topological derivative $d\mathcal{J}$ (2.23) with respect to the
adjoint state as well as some rotational symmetry property of $K$ with respect
to $\mathcal{D}u_{0}(z)$ to directly precompute the topological derivative
$d\mathcal{J}$ for some basis elements and to obtain the full topological
derivative as a linear combination of these values. We remark that such a
procedure – if applicable – can be computationally cheaper and allow for the
use of the obtained formulas in iterative topology optimization algorithms,
see e.g. [6, 14] for more details.

## Conclusion and Outlook

We have illustrated a systematic way of computing topological derivatives for
a large class of PDE-constrained topology optimization problems. Using
automated differentiation in NGSolve, it was possible to compute all potential
terms of the topological derivative for an abstract problem class. We showed
the effectivity of the method in several linear and nonlinear model problems
and verified the topological derivative formulas by computing the values of
the topological asymptotic expansion that should be satisfied by the
topological derivative.

The presented work can be extended in several directions. On the one hand, we
only considered problems posed in an $H^{1}$ setting. We remark that the same
procedure is possible for any differential operator $\mathcal{D}$ that scales
as $(\mathcal{D}\varphi)\circ
T_{\varepsilon}=\frac{1}{{\varepsilon}}\mathcal{D}(\varphi\circ
T_{\varepsilon})$. In particular, this also includes the $curl$ operator, thus
also allowing to treat $H(curl)$ problems as they arise in electromagnetics
[14]. Moreover, for some linear problems, the computational effort for solving
the corrector equation (2.24) could be reduced by employing a boundary element
method rather than a finite element method on an approximation of the
unbounded domain. Finally, a further topic of future research is the numerical
analysis of the topological derivative with respect to mesh sizes in the
bounded domain $\mathsf{D}$ and in the blown-up domain $B_{R}$ as well as with
respect to the radius $R$ of the large domain.

## References

  * [1] G. Allaire and F. Jouve.  Coupling the level set method and the topological gradient in structural optimization.  In Springer, editor, IUTAM Symposium on Topological Design Optimization of Structures, Machines and Materials, pages 3–12, 2006.
  * [2] H. Ammari and H. Kang.  Polarization and Moment Tensors.  Springer-Verlag New York, 2007.
  * [3] S. Amstutz.  Sensitivity analysis with respect to a local perturbation of the material property.  Asymptotic analysis, 49(1), 2006.
  * [4] S. Amstutz and H. Andrä.  A new algorithm for topology optimization using a level-set method.  Journal of Computational Physics, 216(2):573–588, 2006.
  * [5] S. Amstutz and A. Bonnafé.  Topological derivatives for a class of quasilinear elliptic equations.  Journal de Mathématiques Pures et Appliquées, 107(4):367–408, 2017.
  * [6] S. Amstutz and P. Gangl.  Topological derivative for the nonlinear magnetostatic problem.  Electron. Trans. Numer. Anal., 51:169–218, 2019.
  * [7] S. Amstutz, A.A. Novotny, and N. Van Goethem.  Topological sensitivity analysis for elliptic differential operators of order 2m.  Journal of Differential Equations, 256(4):1735–1770, 2014.
  * [8] P. Baumann and K. Sturm.  Adjoint based methods for the computation of higher order topological derivatives with an application to linear elasticity.  arXiv, 2021.
  * [9] M. Burger, B. Hackl, and W. Ring.  Incorporating topological derivatives into level set methods.  Journal of Computational Physics, 194(1):344–362, 2004.
  * [10] M.C. Delfour.  Control, shape, and topological derivatives via minimax differentiability of lagrangians.  In Springer INdAM Series, pages 137–164. Springer International Publishing, 2018.
  * [11] H. A. Eschenauer, V. V. Kobelev, and A. Schumacher.  Bubble method for topology and shape optimization of structures.  Structural optimization, 8(1):42–51, 1994.
  * [12] P. Gangl and K. Sturm.  A simplified derivation technique of topological derivatives for quasi-linear transmission problems.  ESAIM Control Optim. Calc. Var., 26:Paper No. 106, 20, 2020.
  * [13] P. Gangl and K. Sturm.  Topological derivative for PDEs on surfaces, 2020.
  * [14] P. Gangl and K. Sturm.  Asymptotic analysis and topological derivative for 3d quasi-linear magnetostatics.  ESAIM: Mathematical Modelling and Numerical Analysis, 55:S853–S875, 2021.
  * [15] P. Gangl and K. Sturm.  Supplementary code for article ”Automated computation of topological derivatives with application to nonlinear elasticity and reaction-diffusion problems”, 2021.  https://doi.org/10.5281/zenodo.5118379.
  * [16] P. Gangl, K. Sturm, M. Neunteufel, and J. Schöberl.  Fully and semi-automated shape differentiation in NGSolve.  Structural and Multidisciplinary Optimization, 63(3):1579–1607, November 2020.
  * [17] S. Garreau, P. Guillaume, and M. Masmoudi.  The Topological Asymptotic for PDE Systems: The Elasticity Case.  SIAM Journal on Control and Optimization, 39(6):1756–1778, 2001\.
  * [18] M. Hintermüller.  Fast level set based algorithms using shape and topological sensitivity information.  Control and Cybernetics, 34(1), 2005.
  * [19] M. Hintermüller, A. Laurain, and A. A. Novotny.  Second-order topological expansion for electrical impedance tomography.  Advances in Computational Mathematics, 36(2):235–265, 2012.
  * [20] A. A. Novotny and J. Sokołowski.  Topological derivatives in shape optimization.  Interaction of Mechanics and Mathematics. Springer, Heidelberg, 2013.
  * [21] J. Schöberl.  C++11 implementation of finite elements in NGSolve.  Technical Report 30, Institute for Analysis and Scientific Computing, Vienna University of Technology, 2014.
  * [22] J. Sokołowski and A. Zochowski.  On the topological derivative in shape optimization.  SIAM Journal on Control and Optimization, 37(4):1251–1272, 1999\.
  * [23] K. Sturm.  Topological sensitivities via a Lagrangian approach for semilinear problems.  Nonlinearity, 33(9):4310–4337, jul 2020.