# Numerical shape and topology optimization of regions supporting the boundary conditions of a physical problem

Eric Bonnetier, Carlos Brito-Pacheco, Charles Dapogny, Rafael Estevez

► To cite this version:

Eric Bonnetier, Carlos Brito-Pacheco, Charles Dapogny, Rafael Estevez. Numerical shape and topology optimization of regions supporting the boundary conditions of a physical problem. ESAIM: Control, Optimisation and Calculus of Variations, 2025, 31, pp.75. 10.1051/cocv/2025061 hal-05291243

HAL Id: hal-05291243

https://hal.science/hal-05291243v1

Submitted on 3 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
NUMERICAL SHAPE AND TOPOLOGY OPTIMIZATION
OF REGIONS SUPPORTING THE BOUNDARY CONDITIONS
OF A PHYSICAL PROBLEM

ERIC BONNETIER<fnref n="1" />, CARLOS BRITO-PACHECO<fnref n="2" />, CHARLES DAPOGNY<fnref n="2" />,*, AND RAFAEL ESTEVEZ<fnref n="3" />

**Abstract.** This article deals with a particular class of shape and topology optimization problems: the optimized design is a region *G* of the boundary ∂Ω of a given domain Ω, which supports a particular type of boundary conditions in the state problem characterizing the physical situation. In our analyses, we develop adapted versions of the notions of shape and topological derivatives, which are classically tailored to functions of a “bulk” domain. This leads to two complementary notions of derivatives for a quantity of interest *J*(*G*) depending on a region *G* ⊂ ∂Ω: on the one hand, we elaborate on the boundary variation method of Hadamard for evaluating the sensitivity of *J*(*G*) with respect to “small” perturbations of the boundary of *G* within ∂Ω. On the other hand, we use techniques from asymptotic analysis to appraise the sensitivity of *J*(*G*) with respect to the addition of a new connected component to the region *G*, shaped as a “small” surface disk. The calculation of both types of derivatives raises original difficulties, which are closely related to the weakly singular behavior of the solution to a boundary value problem at the points of ∂Ω where the boundary conditions change types. These aspects are carefully detailed in a simple mathematical setting based on the conductivity equation. We notably propose formal arguments to calculate our derivatives with a minimum amount of technicality, and we show how they can be generalized to handle more intricate problems, arising for instance in the physical contexts of acoustics and structural mechanics, respectively governed by the Helmholtz equation and the linear elasticity system. In numerical applications, our derivatives are incorporated into a recent algorithmic framework for tracking arbitrarily dramatic motions of a region *G* within a fixed ambient surface, which combines the level set method with remeshing techniques to offer a clear, body-fitted discretization of the evolving region. Finally, various 3d numerical examples are presented to illustrate the salient features of our analysis.

**Mathematics Subject Classification.** 49M41, 49Q10, 49J50, 65K10, 74P15, 74P20.

Received March 10, 2025. Accepted July 3, 2025.

<footnote>*Keywords and phrases:* Shape optimization, topology optimization, boundary conditions, partial differential equations, level set method, asymptotic analysis.</footnote>

<footnote n="1">¹ Institut Fourier, Université Grenoble-Alpes, BP 74, 38402 Saint-Martin-d'Hères Cedex, France.</footnote>

<footnote n="2">² Univ. Grenoble Alpes, CNRS, Grenoble INP, LJK, 38000 Grenoble, France.</footnote>

<footnote n="3">³ Univ. Grenoble-Alpes – CNRS UMR 5266, SIMaP, F-38000 Grenoble, France.</footnote>

<footnote>* Corresponding author: charles.dapogny@univ-grenoble-alpes.fr</footnote>
# 1. INTRODUCTION

The need for energy savings and the looming concerns about material scarcity have raised a considerable interest in shape and topology optimization in the academic and industrial communities. These techniques find applications in a large variety of areas, such as structural mechanics [1, 2], civil engineering and architecture [3, 4], fluid mechanics [5, 6], electromagnetism [7–10], and biomedical engineering [11–13].

In the classical instances of such problems, the design under scrutiny is a “bulk” domain $\Omega$ in $\mathbb{R}^d$ ($d = 2, 3$), which is optimized with respect to a performance criterion $J(\Omega)$, under constraints about e.g. its volume or perimeter. In most applications, $J(\Omega)$ depends on a physical “state” function $u$, characterized as the solution to a partial differential equation posed on $\Omega$. Most often, the regions of $\partial\Omega$ supporting specific boundary conditions attached to the latter are imposed by the context, and are not subject to optimization.

The present article investigates optimal design problems where the design variable is precisely one of those regions of $\partial\Omega$ supporting a particular type of boundary conditions in the formulation of the underlying physical problem. Here are a few examples:

* In electrostatics, $\Omega$ represents a conductor and the state is the voltage potential $u: \Omega \to \mathbb{R}$, solution to the conductivity equation. It is grounded on a subset $\Gamma_D$ of $\partial\Omega$ and a flux $g: \Gamma_N \to \mathbb{R}$ is imposed on a disjoint region $\Gamma_N \subset \partial\Omega$: these effects are modeled by a homogeneous Dirichlet condition on $\Gamma_D$ and an inhomogeneous Neumann condition on $\Gamma_N$. The remaining part $\Gamma$ of $\partial\Omega$, which is insulated from the outside, is subject to a homogeneous Neumann boundary condition. Although $\Gamma_D$ and $\Gamma_N$ are usually fixed, one may wish to minimize (or maximize) the amplitude of the electric field in $\Omega$ with respect to their placement on $\partial\Omega$.

* In acoustics, $u: \Omega \to \mathbb{R}$ is the sound pressure in a room $\Omega$, solution to the Helmholtz equation. The boundary $\partial\Omega$ is decomposed into two regions $\Gamma_N$ and $\Gamma_R$: Neumann boundary conditions are imposed on $\Gamma_N$, where an incoming wave undergoes perfect reflection, while $\Gamma_R$ bears Robin boundary conditions, which account for a partial absorption. One may then seek to arrange $\Gamma_N$ and $\Gamma_R$ within $\partial\Omega$ so as to minimize the sound pressure in $\Omega$.

* In structural mechanics, $\Omega$ is a mechanical part, attached on a subset $\Gamma_D$ of its boundary $\partial\Omega$, and submitted to surface loads $g: \Gamma_N \to \mathbb{R}^d$ on a disjoint region $\Gamma_N \subset \partial\Omega$. The vector field $u: \Omega \to \mathbb{R}^d$, represents the displacement of the structure. It is the solution to the linear elasticity system. Usually, $\Gamma_D$ and $\Gamma_N$ are given by the context, and only the remaining, traction-free boundary $\Gamma$ is optimized. However, it may be relevant to optimize the placement of the fixation region $\Gamma_D$ to minimize the displacement of the structure.

These questions fit in the general shape optimization framework of a subset $G$ of a fixed ambient surface $S \subset \mathbb{R}^d$. Early studies in this context are devoted to the simulation of geometric flows within $S$, notably the mean curvature flow. In [14], this task is investigated thanks to the level set method, in a situation where $S$ is equipped with a triangular mesh; it is also considered in [15], where $S$ itself is defined in an implicit way. The perhaps most natural instance of a physical optimal design problem posed on a surface $S \subset \mathbb{R}^d$ concerns the optimal reinforcement of a shell structure: $S$ then plays the role of the midsurface, in which the optimized region $G$ is that made of a stiffer material. Popular numerical strategies feature a fixed mesh of $S$ which serves as the support of density-based topology optimization techniques [16, 17]. The level set method is also employed in such setting in [18], and it is coupled with a geometric optimization procedure for the midsurface $S$ itself in [19]. On a different note, in [20], the level set method is combined with a conformal mapping strategy, reducing the surface $S$, and thereby the whole shape optimization problem, to a more classical planar situation.

In spite of their natural character and ubiquity in concrete applications, optimization problems of regions supporting the boundary conditions of a physical problem have been rarely considered in the literature. Without anticipating too much on the more complete overview of related contributions in the particular application contexts of Section 7, let us mention that density-based topology optimization methods are prevailing in this context also, see e.g. [21, 22] in the context of the optimal design of a fixture system. In [23, 24], the level set method is used on a fixed mesh of a computational domain for the concurrent optimization of the 2d shape and of the region of the boundary where Dirichlet boundary conditions are imposed; see also [25] for similar ideas.
Closer to the framework of the present article, the article [26] leverages the level set method on a fixed mesh of a box-shaped room $\Omega$ to optimize the distribution of absorbent and sound-hard materials on its boundary.

Optimizing a function $J(G)$ depending on a region $G$ supporting boundary conditions raises challenging issues from various perspectives. From the theoretical viewpoint, the realization of this task requires the calculation of the derivative of $J(G)$. This information is indeed the basic expression of the optimality conditions of the optimization problem. It is also the main ingredient of iterative optimization algorithms, starting from the simple, unconstrained gradient descent method to more advanced constrained optimization algorithms, such as those proposed in [27–29]. Interestingly, these derivatives can also be used to enforce the robustness of the optimization problem with respect to small perturbations of the geometry of $G$, in the spirit of e.g. [30–32]. From the numerical viewpoint, optimal design problems of regions supporting boundary conditions raise in particular the need to track the possibly dramatic evolution of a region within an ambient surface – a task which is already notoriously difficult when the ambient medium is (a bounded domain of) the Euclidean space $\mathbb{R}^d$.

The present article is the natural continuation of our previous contributions [33–35]. The article [35] deals with the shape sensitivity of a function $J(G)$ depending on a region $G$ of the boundary of a domain $\Omega$ bearing the boundary conditions of a physical problem, in the spirit of the method of Hadamard: the derivative of $J(G)$ with respect to “small”, diffeomorphic perturbations of $G$ is considered. The situation where $G$ is the support of Dirichlet conditions and the complement $\partial\Omega \setminus \bar{G}$ is equipped with Neumann conditions is of particular interest in that work. Indeed, the weakly singular behavior of the state function $u$ in that case renders the treatment of shape derivatives particularly difficult – a fact which was previously acknowledged in [36]. To alleviate this issue, an approximation of the state problem is proposed, which lends itself to simpler calculations and numerical treatment. These developments pave the way to a numerical algorithm for the shape optimization of the region $G$, see Section 3.3 below for a brief presentation. The article [33] deals with singular perturbations of $G$, at the theoretical level: asymptotic formulas are derived for the solution $u$ to the conductivity equation in the case where homogeneous Neumann boundary conditions are replaced by homogeneous Dirichlet equations (and vice-versa) in a “vanishing” zone $\omega_\varepsilon \subset \partial\Omega$. A preliminary application of these results to the device of a notion of topological derivative for functions $J(G)$ depending on regions $G$ bearing boundary conditions was described in the recent article [34]. The latter stands at the numerical level; elaborating on the ideas of [37], it introduces a body-fitted mesh evolution method to track the (possibly dramatic) motion of a region within a fixed ambient surface, which efficiently combines the level set method with remeshing algorithms. Note that this algorithm was used in the very recent work [38] in view of optimizing Neumann eigenvalues on the unit 3d sphere.

In the present article, we leverage these ideas to propose a general shape and topology optimization framework for a region $G \subset \partial\Omega$ supporting the boundary conditions attached to a partial differential equation set on the fixed ambient domain $\Omega$. Our strategy combines (adapted versions of) the notions of shape and topological derivatives to appraise the sensitivity of a function $J(G)$ with respect to small, diffeomorphic perturbations of $\partial G$ and to singular perturbations, via the addition of a small surface disk $\omega_{x_0,\varepsilon}$ around a point $x_0 \in \partial\Omega$, respectively. From the theoretical vantage, one of our contributions is to provide formal calculation methods for both types of derivatives. These are detailed in the simple setting of electrostatics, and they can be adapted to treat several novel situations, in the more intricate contexts of acoustics and structural mechanics. Notably, the calculation of a topological derivative for a function $J(G)$ depending on a region $G$ bearing boundary conditions is rather subtle, and we describe how our methods can be adapted to achieve this purpose in each situation. Interestingly, this article illustrates two different forays of asymptotic analysis in the realm of shape and topology optimization: on the one hand, these concepts are used to smoothen a singular transition between two zones bearing different boundary conditions, leading to a simplified calculation of shape derivatives. On the other hand, they allow to investigate singular perturbations of a smooth background problem, via the insertion of a “small” zone where boundary conditions are altered.

This article is organized as follows. In Section 2, we define the versions of shape and topological derivatives adapted to our purpose of evaluating the sensitivity of a function $J(G)$ depending on a boundary region $G$, and we sketch a generic optimal design algorithm based on these ingredients. The next Sections 3 and 4 concern the model physical context of electrostatics, which conveniently allows us to expose the salient features of the proposed methodology with a minimum amount of technicality. In Section 3, we describe the calculation of
the shape derivative of a functional $J(G)$ which depends on the region supporting the homogeneous Dirichlet boundary conditions attached to the conductivity equation. We also outline the treatment of the simpler problem where the optimized region $G$ is that supporting inhomogeneous Neumann boundary conditions, and we address other easy extensions of this analysis. In the same physical situation, Section 4 deals with the delicate calculation of the topological derivative of $J(G)$. In Section 5 and Section 6, we generalize this material to acoustics and structural mechanics, accounted for by the Helmholtz and linear elasticity equations, respectively: the calculations of approximate or exact shape derivatives for the typical functionals $J(G)$ of interest follow exactly the same trail as in the electrostatics setting. We explicitly compute the topological derivatives in these cases. We then turn to numerical applications in Section 7, where several examples are analyzed in the settings of electrostatics, acoustics and structural mechanics, in order to appraise the main features of our methods. Finally, a few conclusive remarks and perspectives for future research are drawn in Section 8.

## 2. SENSITIVITY OF A FUNCTION DEPENDING ON A SUBSET OF THE BOUNDARY OF A DOMAIN

Let $\Omega$ be a smooth bounded domain in $\mathbb{R}^d$ ($d = 2, 3$). We consider a model shape and topology optimization problem of the form:

$$ \min_{G \subset \partial \Omega} J(G), \qquad (P) $$

where $J(G)$ is a given objective function depending on the region $G \subset \partial \Omega$ to be optimized. This formulation omits constraints for simplicity of the presentation, but their treatment with the material developed in this article does not entail any additional difficulty from the conceptual viewpoint.

The analysis of the problem (P) hinges on the sensitivity of $J(G)$ with respect to “small” variations of $G$. In this section, we introduce two complementary means to appraise the notion of derivative with respect to a boundary region:

* In Section 2.1, we adapt the “classical” boundary variation method of Hadamard to the context of a subset $G$ of the boundary $\partial \Omega$. This paves the way to a notion of shape derivative, that accounts for the sensitivity of $J(G)$ with respect to small perturbations of the boundary of $G$.

* In Section 2.2, we define a suitable version of the concept of topological derivative, which appraises the sensitivity of $J(G)$ with respect to the addition of a “small” surface disk to $G$.

A sketch of our numerical solution strategy for (P) based on these ingredients is provided in Section 2.3.

**Notation.** Throughout this article, $\Omega \subset \mathbb{R}^d$ is a fixed, smooth bounded domain, and the optimized design $G$ is a smooth open subset of its boundary $\partial \Omega$. Moreover,

* At any point $x \in \partial \Omega$, we denote by $n(x)$ the unit normal vector to $\partial \Omega$ at $x$, pointing outward $\Omega$.

* The tangent plane $T_x \partial \Omega$ at $x \in \partial \Omega$ is the vector hyperplane defined by

$$ T_x \partial \Omega := \{ \tau = (\tau_1, \dots, \tau_d) \in \mathbb{R}^d \text{ s.t. } n(x) \cdot \tau = 0 \}. $$

* The tangential gradient of a smooth function $f: \partial \Omega \to \mathbb{R}$ is denoted by $\nabla_{\partial \Omega} f$. By definition, $\nabla_{\partial \Omega} f = \nabla \tilde{f} - (\nabla \tilde{f} \cdot n)n$, where $\tilde{f}: \mathbb{R}^d \to \mathbb{R}$ is an arbitrary smooth extension of $f$ to $\mathbb{R}^d$.

* We denote by $\text{div}_{\partial \Omega}(V)$ the tangential divergence of a smooth vector field $V: \partial \Omega \to \mathbb{R}^d$. By definition, $\text{div}_{\partial \Omega}(V) = \text{div}(\tilde{V}) - \nabla \tilde{V} \cdot n$, for any smooth extension $\tilde{V}: \mathbb{R}^d \to \mathbb{R}^d$ of $V$ to $\mathbb{R}^d$.

* We denote by $\Sigma = \partial G$ the boundary of $G$, and by $d\ell$ the integration on $\Sigma$, i.e. the restriction of the $(d-2)$-dimensional Hausdorff measure on this set.

* For any point $x \in \Sigma$, $n_{\Sigma}(x) \in T_x \partial \Omega$ is the unit conormal vector to $\Sigma$ at $x$, pointing outward $G$.

* Obviously, the functionals $J(G)$ of interest depend not only on the region $G$ of $\partial \Omega$, but also on the whole domain $\Omega$. As the latter is fixed, this dependence is omitted.

We refer to Figure 1 for an illustration of some of these definitions.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Perturbation $G_{\theta}$ of a region $G \subset \partial\Omega$ induced by a “small” vector field $\theta$ in Section 2.1. When $\theta \cdot n \equiv 0$, the boundary $\partial\Omega_{\theta}$ of the perturbed domain $\Omega_{\theta}$ coincides with $\partial\Omega$ at first order, and the deformation $\theta$ causes G to “slide” within $\partial\Omega$.</figcaption></figure>

## 2.1. Hadamard's boundary variation method

In this section, we adapt the classical method of Hadamard to estimate the sensitivity of a functional $J(G)$ with respect to perturbations of the boundary $\Sigma$ of G. Following [39–42], we consider variations of the reference domain $\Omega$ of the form:

$$ \Omega_{\theta} := (\mathrm{Id} + \theta)(\Omega), \quad \text{where} \quad \theta \in W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d), \quad \|\theta\|_{W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)} < 1. \qquad (2.1) $$

The variations $G_{\theta}$ of $G \subset \partial\Omega$ induced by this operation read:

$$ G_{\theta} := (\mathrm{Id} + \theta)(G), \qquad (2.2) $$

see Figure 1. This leads to the following definition of the shape derivative of a functional J depending on a region $G \subset \partial\Omega$.

**Definition 2.1.** The functional $J(G)$ is shape differentiable at $G \subset \partial\Omega$ if the mapping $\theta \mapsto J(G_{\theta})$, defined from a neighborhood of 0 in $W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)$ into $\mathbb{R}$, is Fréchet differentiable at $\theta = 0$. The corresponding shape derivative $J'(G)(\theta)$ then satisfies the following expansion:

$$ J(G_{\theta}) = J(G) + J'(G)(\theta) + o(\theta), \quad \text{where} \quad \frac{o(\theta)}{\|\theta\|_{W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)}} \xrightarrow{\theta \to 0} 0. \qquad (2.3) $$

The shape derivative $J'(G)(\theta)$ of a functional $J(G)$ allows to identify descent directions for $J(G)$, i.e. vector fields $\theta$ such that $J'(G)(\theta) < 0$. Loosely speaking, the region $G_{\tau\theta}$ resulting from the deformation of G via such a vector field $\theta$, for a small enough (pseudo-)time step $\tau > 0$, performs “better” than G, i.e.:

$$ J(G_{\tau\theta}) \approx J(G) + \tau J'(G)(\theta) < 0. $$
**Remark 2.2.** In practice, we shall restrict the set of deformations $\theta$ featured in our implementation of Hadamard's method so as to accommodate some specificities of our applications:

(i) We are especially interested in deformations $\theta$ which leave the domain $\Omega$ invariant, thus giving rise to perturbations $G_\theta$ resulting from a “sliding” of $G$ within $\partial\Omega$. Hence, we shall often restrict ourselves to tangential vector fields $\theta$, satisfying $\theta \cdot n \equiv 0$ on $\partial\Omega$. Note that, strictly speaking, these do not leave $\Omega$ invariant: they only do so “at first order”. Actually, if variations of $\Omega$ were defined within the alternative framework of the so-called velocity method [43, 44] instead of that of the method of Hadamard, i.e. if we were to set

$\Omega_\theta = \chi_\theta(1, \Omega)$, where $(t, x) \mapsto \chi_\theta(t, x)$ is the flow of $\theta$, defined by

$$ \left\{ \begin{array}{ll} \frac{\partial \chi_\theta}{\partial t}(t, x) = \theta(\chi_\theta(t, x)) & \text{for } t > 0, \\ \chi_\theta(0, x) = x, & \end{array} \right. $$

then the action of a tangential vector field $\theta$ on $\Omega$ would ensure that $\Omega_\theta$ coincide with $\Omega$. For simplicity, and since the velocity method and the method of Hadamard induce the same notion of first-order shape derivative, we ignore this technicality in the following.

(ii) The boundary $\partial\Omega$ is often composed of several regions, $G$ being the only one subjected to optimization. Accordingly, the deformations $\theta$ considered in the method of Hadamard should vanish on those other regions of $\partial\Omega$ excluded from the optimization.

(iii) The handling of some geometric quantities attached to $G$ requires the latter to be “smooth enough”, as well as its variations $G_\theta$, see e.g. the results summarized in Section 3.2 about the geodesic signed distance function.

As a result, we shall restrict the set of deformations $\theta$ to a subset $\Theta_G \subset W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)$ made of smooth tangential vector fields, vanishing “far” from the optimized region $G$.

The shape derivatives of the functionals $J(G)$ considered in this article often turn out to be of the form

$$ \forall \theta \in \Theta_G, \quad J'(G)(\theta) = \int_{\Sigma} v_G \theta \cdot n_{\Sigma} \, d\ell, \qquad (2.4) $$

where the scalar field $v_G : \Sigma \to \mathbb{R}$ depends on the region $G$ and the considered objective function $J(G)$. It may involve the solution to one or several boundary value problems attached to $\partial\Omega$ and $G$, see for instance Section 3.3 below. The structure (2.4) is very reminiscent of the “classical” structure theorem for the derivative of functions of a bulk domain $\Omega$, see e.g. Section 5.9 in [41]. It conveniently reveals a descent direction for $J(G)$; indeed, letting $\theta = -v_G n_\Sigma$ leads to:

$$ J'(G)(\theta) = - \int_{\Sigma} v_G^2 d\ell < 0. $$

Before closing this section, let us introduce two simple and ubiquitous examples of functionals depending on the boundary region $G$, namely the area $\text{Area}(G)$ and contour $\text{Cont}(G)$, which are respectively defined by:

$$ \text{Area}(G) = \int_G ds, \quad \text{Cont}(G) = \int_{\Sigma} d\ell. \qquad (2.5) $$

The next proposition supplies the shape derivatives of slightly more general versions of these quantities. Its proof is an elementary adaptation of classical arguments, found in e.g. [41] and Chap. 17 of [45], and it is omitted for brevity.

**Proposition 2.3.** Let $G$ be smooth region of the boundary $\partial\Omega$. Then,
(i) For any smooth function $f : \mathbb{R}^d \to \mathbb{R}$, the functional $J(G)$ defined by:

$$J(G) := \int_G f \, ds$$

is shape differentiable; its shape derivative reads, for any tangential deformation $\theta$ (i.e. $\theta \cdot n = 0$):

$$J'(G)(\theta) = \int_{\Sigma} f \theta \cdot n_{\Sigma} \, d\ell.$$

(ii) For any smooth function $g : \mathbb{R}^d \to \mathbb{R}$, the functional $K(G)$ defined by:

$$K(G) := \int_{\Sigma} g \, d\ell$$

is shape differentiable; its shape derivative reads, for any tangential deformation $\theta$:

$$K'(G)(\theta) = \int_{\Sigma} (\nabla_{\partial\Omega} g \cdot n_{\Sigma} + \kappa g)\theta \cdot n_{\Sigma} \, d\ell,$$

where $\kappa := \text{div}_{\partial\Omega}(n_{\Sigma})$ is the mean curvature of $\Sigma$.

## 2.2. Topological perturbations of subsets of the boundary of a domain

The variations $G_{\theta}$ of the region $G \subset \partial\Omega$ featured in the method of Hadamard of Section 2.1 share the same topology as $G$, i.e. they have the same number of holes, and the same number of connected components [46]. Hence, strictly speaking, the notion of shape derivative in Definition 2.1 does not account for the sensitivity of $J(G)$ to topological changes in $G$. Admittedly, in practice, certain topological changes may occur as $G$ evolves according to (2.2): separate parts of the boundary $\Sigma$ of $G$ may collide and merge by a slight abuse of the theoretical framework (i.e. when the norm constraint $||\theta||_{W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)} < 1$ is omitted in (2.2)), but in any event, no hole can emerge inside $G$, and no new connected component can be added to $G$, see [39] for related discussions.

We introduce in this section a complementary notion of derivative for functions $J$ of a region $G$ of $\partial\Omega$, which proceeds via the addition of a small surface disk to $G$. This notion is an adaptation of that of topological derivative for functions depending on a “bulk” domain, which evaluates the removal of an infinitesimally small hole from the latter; we refer to [47–49] for the seminal works about topological derivatives, see also [50, 51] for overviews and the contribution [52] dealing with holes located close to the boundary of the domain. Let us denote by

$$H := \{x = (x_1, \dots, x_d) \text{ s.t. } x_d < 0\}, \quad \text{and} \quad D_{\varepsilon} := \{x = (x_1, \dots, x_{d-1}, 0) \in \partial H \text{ s.t. } |x| < \varepsilon\} \quad (2.6)$$

the lower half-space in $\mathbb{R}^d$ and the planar disk with center 0 and radius $\varepsilon$, respectively. Let $x \in \partial\Omega \setminus \bar{G}$; since the domain $\Omega$ is smooth, there exists an open neighborhood $O$ of 0 in $\mathbb{R}^d$ and a smooth diffeomorphism $T: O \to T(O)$ such that:

$$T(0) = x \quad \text{and} \quad T(H \cap O) = \Omega \cap T(O). \qquad (2.7)$$

We then define the surface disk $\omega_{x,\varepsilon} \subset \partial\Omega$ with center $x$ and radius $\varepsilon$ by:

$$\omega_{x,\varepsilon} = T(D_{\varepsilon}), \qquad (2.8)$$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Variation $G_{x,\epsilon}$ in (2.9) of a surface region $G \subset \partial\Omega$ obtained by addition of a “small” surface disk $\omega_{x,\epsilon}$.</figcaption></figure>

and we consider variations $G_{x,\epsilon}$ of $G$ of the form:

$$G_{x,\epsilon} := G \cup \omega_{x,\epsilon}, \qquad (2.9)$$

see Figure 2 for an illustration.

**Definition 2.4.** A function $J(G)$ has a topological derivative $d_T J(G)(x)$ at $x \in \partial\Omega \setminus \bar{G}$ if there exists a function $\rho: \mathbb{R}_+ \to \mathbb{R}_+$ such that $\rho(\epsilon) \to 0$ as $\epsilon \to 0$, and the following asymptotic expansion holds:

$$J(G_{x,\epsilon}) = J(G) + \rho(\epsilon)d_T J(G)(x) + o(\rho(\epsilon)).$$

Loosely speaking, the negative values of $d_T J(G)(x)$ indicate that it is beneficial to attach a “small” surface disk $\omega_{x,\epsilon}$ to $G$ to decrease the value of the criterion $J(G)$.

**Remark 2.5.** At first sight, this definition might seem to depend on the choice of the local diffeomorphism $T$ in (2.7); the analysis in the forthcoming sections reveals that this is not the case.

**Remark 2.6.** One could alternatively think of topological perturbations of $G$ of the form $G \setminus \bar{\omega}_{x,\epsilon}$, $x \in G$, i.e. perturbations resulting from the removal of a “small” surface disk from $G$. In our applications, where $G$ (resp. $\partial\Omega \setminus \bar{G}$) typically bears the homogeneous Dirichlet (resp. homogeneous Neumann) conditions of a boundary value problem, this would correspond to replacing a homogeneous Dirichlet boundary condition by a homogeneous Neumann condition in a “small” surface region $\omega_{x,\epsilon} \subset G$, instead of replacing a homogeneous Neumann condition by a homogeneous Dirichlet condition, as we do by considering perturbations of $G$ of the form (2.9). As explained in Remark 4.7, such issues can be handled with a similar analysis as that discussed in this article, so that Definition 2.4 is sufficient for our purpose.

## 2.3. Description of the general optimization strategy

The notions of shape and topological derivatives of functions of a region $G \subset \partial\Omega$, introduced in Sections 2.1 and 2.2, pave the way to an optimization strategy for the solution of the problem (P), a generic version of
which is sketched in Algorithm 1. This method is reminiscent of the algorithm coupling shape and topological derivatives proposed in [53] to deal with the optimization of “bulk” shapes.

Starting from an initial configuration $G^0$, the algorithm proceeds iteratively, producing shapes, states, descent directions, etc., indexed by the integer $n = 0, \dots$ and all the corresponding instances of the various objects at play are denoted with an $n$ superscript. At the $n^{th}$ stage of the process, the shape derivative $J'(G^n)(\theta)$ of the objective function is calculated at $G^n$, and a descent direction $\theta^n$ is inferred from the structure (2.4). The new shape $G^{n+1} := (\text{Id} + \tau^n \theta^n)(G^n)$ is then obtained by displacing $G^n$ in this direction for a suitably small time step $\tau^n$. Every $n_{\text{top}}$ iteration, this procedure is replaced with the calculation of the topological derivative $d_T J(G^n)(x)$, and the addition of a small surface disk $\omega_{x,\varepsilon}$ to $G^n$ around the point $x \in \partial\Omega \setminus G^n$ where $d_T J(G^n)(x)$ takes the largest negative value.

<caption>**Algorithm 1:** Shape and topology optimization of a region $G \subset \partial\Omega$.</caption>

**Initialization:**
* Fixed domain $\Omega \subset \mathbb{R}^d$,
* Initial region $G^0 \subset \partial\Omega$.

**for** $n = 0, \dots$, until convergence **do**
  **if** $n \bmod. n_{\text{top}} = 0$ **then**
    (1) Calculate the topological derivative $d_T J(G^n)(x)$;
    (2) Find the point $x \in \partial\Omega$ where $d_T J(G^n)(x)$ is minimum (negative);
    (3) Update $G^n$ as $G^{n+1} = G^n \cup \omega_{x,\varepsilon}$ for $\varepsilon$ small enough.
  **else**
    (1) Calculate the shape derivative $J'(G^n)(\theta)$;
    (2) Find a descent direction $\theta^n$;
    (3) Update $G^n$ as $G^{n+1} := (\text{Id} + \tau^n \theta^n)(G^n)$, where $\tau^n$ is a suitably small time step.
  **end if**
**end for**
**return** Optimized region $G^n \subset \partial\Omega$.

The key ingredients of this workflow are the analytical expressions of the shape and topological derivatives of $J(G)$, which are extensively discussed in the sequel. In the next Sections 3 and 4, we detail the case of the conductivity equation: we notably give a sketch of the proofs of the results contained in our former articles [33, 35], which can be easily reproduced and adapted to different situations. We discuss the more intricate physical contexts of the Helmholtz equation and the linear elasticity system, where such results are new to the best of our knowledge, in the next Sections 5 and 6, mainly highlighting the differences. A more practical version of Algorithm 1, tailored to our numerical setting is provided in Section 7, see Algorithm 2.

**Remark 2.7.** The framework introduced in Sections 2.1 and 2.2 can be adapted to handle the joint optimization of the ambient domain $\Omega$ and the region of interest $G$ on the boundary of the latter. Beyond technical implementation aspects, the main needed modification to achieve this goal is to allow for deformations $\theta$ that are not necessarily tangential to $\partial\Omega$ when computing the Hadamard derivative, see (2.1) and (2.2). This would lead to a more complex expression of the shape derivative of $J(G)$, that includes contributions from the normal component $\theta \cdot n$, thus encoding how to optimize the shape of $\Omega$ independently of that of the region $G$ on its boundary. This possibility, which was considered in our earlier work [35], is ignored in the present study, which focuses on the optimization of a region on the boundary of a fixed domain.
### 3. SHAPE DERIVATIVE OF A FUNCTIONAL OF A REGION SUPPORTING THE BOUNDARY CONDITIONS OF THE CONDUCTIVITY EQUATION

This section takes place in the physical setting of electrostatics, governed by the conductivity equation, where the salient points of our investigations can be presented with a minimum amount of technicality. After introducing the physical problem at stake in Section 3.1 and summarizing a few facts about geodesic signed distance functions in Section 3.2, we focus in Section 3.3 on the case of a functional depending on the region $G \subset \partial\Omega$ supporting a homogeneous Dirichlet boundary condition. We notably explain why “exact” shape derivatives with respect to a transition zone between homogeneous Dirichlet and Neumann boundary conditions are difficult to handle, and we propose a suitable approximate version of these derivatives. In Section 3.4, we provide practical, implementation friendly versions of these formulas. A few extensions of this material are discussed in Section 3.5, including the simpler treatment of a functional of the region of $\partial\Omega$ bearing inhomogeneous Neumann boundary conditions.

#### 3.1. A brief presentation of the conductivity equation

Let $\Omega \subset \mathbb{R}^d$ be a smooth bounded domain, whose boundary $\partial\Omega$ is made of three disjoint parts:

$$ \partial\Omega = \overline{\Gamma_D} \cup \overline{\Gamma_N} \cup \Gamma \text{ such that } \overline{\Gamma_D} \cap \overline{\Gamma_N} = \emptyset. $$

We assume that:

* The voltage potential $u$ is set to 0 on $\Gamma_D$;

* A smooth electric current $g: \mathbb{R}^d \to \mathbb{R}$ is applied on $\Gamma_N$;

* The region $\Gamma$ is insulated from the outside.

The domain $\Omega$ is filled by a material with smooth conductivity $\gamma \in C^\infty(\bar{\Omega})$, which is uniformly bounded away from 0 and $\infty$, i.e.

$$ \exists 0 < \gamma_{-} \le \gamma_{+} < \infty \text{ s.t. } \forall x \in \Omega, \quad \gamma_{-} \le \gamma(x) \le \gamma_{+}. \quad (3.1) $$

In this situation, the voltage potential $u$ induced by a smooth source $f: \mathbb{R}^d \to \mathbb{R}$ belongs to the space

$$ H_{\Gamma_D}^1(\Omega) := \{u \in H^1(\Omega), u = 0 \text{ on } \Gamma_D\}, $$

and it is the unique solution in the latter to the following boundary value problem:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(\gamma \nabla u) = f & \text{in } \Omega, \\ u = 0 & \text{on } \Gamma_D, \\ \gamma \frac{\partial u}{\partial n} = g & \text{on } \Gamma_N, \\ \gamma \frac{\partial u}{\partial n} = 0 & \text{on } \Gamma. \end{array} \right. \quad (3.2) $$

The performance of the configuration is evaluated in terms of the quantity

$$ \int_{\Omega} j(u) \, dx, \quad (3.3) $$

where $j \in C^2(\mathbb{R}^d)$ is a given function, satisfying the following growth conditions:

$$ \exists C > 0 \text{ s.t. } \forall t \in \mathbb{R}^d, \quad |j(t)| \le C(1 + |t|^2), \quad |j'(t)| \le C(1 + |t|), \quad \text{and} \quad |j''(t)| \le C. \quad (3.4) $$
## 3.2. Preliminaries about the geodesic signed distance function in $∂Ω$

In this section we briefly recall some material related to the signed distance function to a subset $G$ of the boundary of the smooth ambient domain $Ω$. We refer to [35] and the contributions mentioned therein for more details.

We first recall the notion of geodesic distance between points on $∂Ω$, see e.g. [54].

**Definition 3.1.** Let $x, y$ be two points in $∂Ω$; the geodesic distance $d^{∂Ω}(x, y)$ between $x$ and $y$ is the shortest length of a curve on $∂Ω$ with endpoints $x, y$:

$$d^{∂Ω}(x, y) = \inf \left\{ \int_0^1 |\gamma'(t)| dt : \gamma \in \mathcal{C}^1([0, 1], ∂Ω), \gamma(0) = x, \gamma(1) = y \right\}.$$

**Definition 3.2.** Let $G$ be a Lipschitz open subset of $∂Ω$ with boundary $Σ = ∂G$. The geodesic signed distance function $d_G^{∂Ω}$ to $G$ is:

$$\forall x \in ∂Ω, \quad d_G^{∂Ω}(x) = \begin{cases} -d^{∂Ω}(x, Σ) & \text{if } x \in G, \\ 0 & \text{if } x \in Σ, \\ d^{∂Ω}(x, Σ) & \text{if } x \in ∂Ω \setminus \overline{G}, \end{cases}$$

where $d^{∂Ω}(x, Σ)$ is the geodesic distance from $x$ to the compact subset $Σ \subset ∂Ω$:

$$d^{∂Ω}(x, Σ) = \min_{y \in Σ} d^{∂Ω}(x, y). \qquad (3.5)$$

**Definition 3.3.** Let $y$ be a point on $∂Ω$;

(1) The set $Π_{Σ}^{∂Ω}(y)$ of projections of $y$ onto $Σ$ is the set of minimizers in (3.5).

(2) When $Π_{Σ}^{∂Ω}(y)$ is a singleton, its unique element $π_{Σ}(y)$ is called the projection of $y$ onto $Σ$.

Let us next recall the definitions of the exponential and logarithm mappings on the hypersurface $∂Ω$.

**Definition 3.4.** Let $x$ be a point on $∂Ω$. The exponential at $x$ is the mapping $exp_x$, defined from a neighborhood $U$ of 0 in $T_x∂Ω$ into a neighborhood $V$ of $x$ in $∂Ω$ by:

$$\forall v \in U, \quad \exp_x(v) = \gamma_v(1),$$

where $t \mapsto \gamma_v(t)$ is the unique geodesic curve passing through $x$ at $t = 0$ with velocity $\gamma_v'(0) = v$.

Up to shrinking $U$ and $V$, $exp_x$ is a smooth diffeomorphism from $U$ onto $V$, whose inverse mapping is denoted by $log_x : V \to U$:

$$\forall y \in V, \quad \log_x(y) \text{ is the unique vector } v \in U \subset T_x∂Ω \text{ s.t. } \exp_x(v) = y.$$

The following result states that when the region $G$ is sufficiently smooth, there exists a tubular neighborhood of the boundary $Σ$ with small width $2δ$ which is diffeomorphic to the product set $Σ \times (-δ, δ)$. It follows from a simple adaptation of the proof of Theorem 3.1 in [55], which is based on the Implicit Function Theorem, see Figure 3 for an illustration.

**Theorem 3.5.** Let $G$ be a smooth open subset of $∂Ω$, with boundary $Σ = ∂G$. Then, there exists a number $δ > 0$ such that the mapping $F : Σ \times (-δ, δ) \to ∂Ω$ defined by:

$$\forall x \in Σ, \ t \in (-δ, δ), \quad F(x, t) = \exp_x(tn_Σ(x))$$
is a smooth diffeomorphism onto its image $U_{\delta} := F(\Sigma \times (-\delta, \delta)) \subset \partial\Omega$. Its inverse reads:

$$ \forall y \in U_{\delta}, \quad F^{-1}(y) = (\pi_{\Sigma}(y), d_{G}^{\partial\Omega}(y)). $$

Let us now provide a handful Fubini-like formula for rewriting integrals posed on the tubular neighborhood $U_{\delta}$ introduced in Theorem 3.5 as nested integrals over $\Sigma$ and the interval $(-\delta, \delta) \subset \mathbb{R}$. This result follows from the change of variables formula on a manifold, see e.g. Chap. XVI in [56]:

**Proposition 3.6.** Let $G$ be a smooth open subset of $\partial\Omega$, let $\delta > 0$ be small enough so that the conclusion of Theorem 3.5 holds true, and let $U_{\delta}$ be the tubular neighborhood of $\Sigma$ introduced in there. Then, a function $\varphi: U_{\delta} \to \mathbb{R}$ is integrable if and only if the composite function $\varphi \circ F$ is integrable on $\Sigma \times (-\delta, \delta)$, and:

$$ \int_{U_{\delta}} \varphi \, ds = \int_{\Sigma} \int_{-\delta}^{\delta} |\det\nabla F(x, t)| \varphi \circ F(x, t) \, dt \, d\ell(x). $$

Here, the definition of the $(d-1) \times (d-1)$ matrix $\nabla F(x, t)$ involves an arbitrary local orthonormal basis $(\tau_1, \dots, \tau_{d-2})$ of tangent vectors to the codimension 1 submanifold $\Sigma$ of $\partial\Omega$. This matrix gathers the components of the tangent vectors

$$ \frac{\partial F}{\partial \tau_1}(x, t), \dots, \frac{\partial F}{\partial \tau_{d-2}}(x, t), \frac{\partial F}{\partial t}(x, t) \in T_{F(x,t)}\partial\Omega, $$

in an arbitrary orthonormal basis of $T_{F(x,t)}\partial\Omega$.

**Remark 3.7.** For further reference, we note the following expression of the determinant $\det\nabla F(x, t)$ when $t=0$:

$$ \forall x \in \Sigma, \quad \det\nabla F(x, 0) = 1. $$

The next result deals with the derivative of the signed distance function.

**Proposition 3.8.** Let $G \subset \partial\Omega$ be a smooth region with boundary $\Sigma$ and let $\delta > 0$ and $U_{\delta}$ be as in the statement of Theorem 3.5. Then the signed distance function $d_{G}^{\partial\Omega}$ is differentiable at any point $x \in U_{\delta}$, and its tangential gradient reads:

$$ \nabla_{\partial\Omega} d_{G}^{\partial\Omega}(x) = -\frac{1}{d_{G}^{\partial\Omega}(x)} \log_{x}(\pi_{\Sigma}(x)). \quad (3.6) $$

Eventually, we shall require the expressions of the derivative of the signed distance function $d_{G}^{\partial\Omega}$ with respect to variations of the region $G$. The following result deals with the “Lagrangian” version of this derivative, that is, the derivative of the mapping

$$ \theta \mapsto d_{G_{\theta}}^{\partial\Omega}(y + \theta(y)), \quad y \in \partial\Omega, \quad (3.7) $$

evaluating the signed distance function to the perturbed region $G_{\theta}$ at the perturbed location $y + \theta(y)$ of a given point $y \in \partial\Omega$.

**Proposition 3.9.** Let $G$ be a smooth open subset of $\partial\Omega$ with boundary $\Sigma$, and let $U_{\delta}$ be a tubular neighborhood of $\Sigma$ supplied by Theorem 3.5. Let $y$ be a fixed point in $U_{\delta} \setminus \bar{\Sigma}$ and let $x := \pi_{\Sigma}(y) \in \Sigma$. For any smooth vector field $\theta$, we define:

$$ D(\theta) := d_{G_{\theta}}^{\partial\Omega}(y + \theta(y)). $$
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Illustration of some of the definitions related to the geodesic signed distance function introduced in Section 3.2.</figcaption></figure>

Then $D(\theta)$ is Fréchet differentiable at $\theta = 0$ and its derivative reads, for a smooth tangential deformation $\theta$:

$$D'(0)(\theta) = -\theta(y) \cdot \frac{\log_y(x)}{d_G^{\partial\Omega}(y)} - \theta(x) \cdot n_{\Sigma}(x).$$

**Remark 3.10.** The combination of Proposition 3.9 with Proposition 3.8 allows to compute the “Eulerian” version of the shape derivative of $G \mapsto d_G^{\partial\Omega}$, i.e. the derivative of the mapping $\theta \mapsto d_{G_\theta}^{\partial\Omega_\theta}(y)$ for a fixed point $y \in \partial\Omega$. Note that this mapping is not well-defined since the point $y$ does not a priori belong to the perturbed boundary $\partial\Omega_\theta$, even when $\theta$ is “small”. However, it is customary to define this “Eulerian” derivative by the following formula, for any tangential deformation $\theta$:

$$ \begin{aligned} (d_G^{\partial\Omega})'(\theta)(y) &:= D'(0)(\theta) - \nabla_{\partial\Omega} d_G^{\partial\Omega}(y) \cdot \theta(y) \\ &= -\theta(x) \cdot n_{\Sigma}(x), \end{aligned} $$

suggested by the (formal) application of the chain rule to (3.7). We refer to [39, 40, 57] about the concepts of “Lagrangian” and “Eulerian” derivatives for functions defined on the domain.

### 3.3. A smoothed approximation method for the optimization of the Dirichlet–Neumann transition

In this section, we focus on the situation where the performance criterion (3.3) is optimized with respect to the shape of the region $\Gamma_D$ supporting the homogeneous Dirichlet conditions in the boundary value problem (3.2). In the language of Section 2, this amounts to setting $G = \Gamma_D$, $\Sigma = \Sigma_D := \partial\Gamma_D$ and to considering the objective functional

$$J(\Gamma_D) = \int_{\Omega} j(u_{\Gamma_D}) dx, \qquad (3.8)$$

where $j$ satisfies (3.4) and we have denoted by $u_{\Gamma_D}$ the solution to (3.2) in this situation. Since in the present context the region $\Gamma_N \subset \partial\Omega$ is not subject to optimization, the set $\Theta_{\Gamma_D}$ of deformations considered in the practice of the method of Hadamard is:

$$\Theta_{\Gamma_D} := \{ \theta : \mathbb{R}^d \to \mathbb{R}^d \text{ is smooth, with } \theta \cdot n = 0 \text{ on } \partial\Omega \text{ and } \theta = 0 \text{ on } \Gamma_N \}, \qquad (3.9)$$

see Remark 2.2.
Computing the shape derivative of the functional $J(\Gamma_D)$ in (3.8) raises multiple difficulties, which are essentially due to the weakly singular behavior of $u_{\Gamma_D}$ at the points $x \in \Sigma_D$ marking the transition between homogeneous Dirichlet and homogeneous Neumann boundary conditions in (3.2). Without entering into details, the variational Lax-Milgram theory guarantees that $u_{\Gamma_D}$ belongs to $H^1(\Omega)$. Under our regularity assumptions on $\Omega$ and $\gamma$, this function is actually smooth, except in the vicinity of $\Sigma_D$, where it does not enjoy $H^2$ regularity. Loosely speaking, this phenomenon is caused by the existence of non trivial “weakly singular” solutions to the homogeneous version of (3.2), that belong to $H^1(\Omega)$, but not to $H^2(\Omega)$, see e.g. [58–60] for a comprehensive treatment of these questions.

From the theoretical viewpoint, this lack of regularity of $u_{\Gamma_D}$ implies that the calculation of the shape derivative $J'(\Gamma_D)(\theta)$ under the convenient form (2.4) cannot be conducted by standard techniques, as those exposed in [39, 40]. In particular, Céa’s formal method for the fast computation of shape derivatives (see [61]) leads to an erroneous result as it tacitly rests on the assumption that $u_{\Gamma_D}$ is “smooth”, see [35, 62] about this point. In addition to being mathematically difficult to calculate, the expression of $J'(\Gamma_D)(\theta)$ is awkward from the numerical viewpoint, as it brings into play the coefficients of a decomposition of $u_{\Gamma_D}$ along the above mentioned weakly singular solutions to the homogeneous version of (3.2), which cannot be computed in closed form in general.

In order to overcome both issues, we follow the idea presented in [35] (see also [63]): we trade the weakly singular solution $u_{\Gamma_D}$ to the “exact” problem (3.2), featuring a sharp transition between Dirichlet and Neumann boundary conditions, for the regularized version $u_{\Gamma_D, \varepsilon}$ characterized as the $H^1(\Omega)$ solution to the following boundary value problem:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(\gamma \nabla u_{\Gamma_D, \varepsilon}) = f & \text{in } \Omega, \\ \gamma \frac{\partial u_{\Gamma_D, \varepsilon}}{\partial n} + h_{\Gamma_D, \varepsilon} u_{\Gamma_D, \varepsilon} = 0 & \text{on } \Gamma \cup \Gamma_D, \\ \gamma \frac{\partial u_{\Gamma_D, \varepsilon}}{\partial n} = g & \text{on } \Gamma_N. \end{array} \right. \qquad (3.10) $$

Here, $h_{\Gamma_D, \varepsilon} : \partial\Omega \to \mathbb{R}$ is the function defined by:

$$ h_{\Gamma_D, \varepsilon}(x) = \frac{1}{\varepsilon} h \left( \frac{d_{\Gamma_D}^{\partial \Omega}(x)}{\varepsilon} \right), \qquad (3.11) $$

where $d_{\Gamma_D}^{\partial \Omega}$ is the geodesic signed distance function to $\Gamma_D$ on $\partial\Omega$ (see Sect. 3.2), and $h \in C^\infty(\mathbb{R})$ satisfies:

$$ 0 \le h \le 1, \quad h \equiv 1 \text{ on } (-\infty, -1), \quad h(0) > 0, \quad \text{and} \quad h \equiv 0 \text{ on } [1, \infty). \qquad (3.12) $$

Intuitively, as $\varepsilon$ gets “small”, $h_{\Gamma_D, \varepsilon}$ takes very large values on $\bar{\Gamma}_D$, where the Robin boundary condition in (3.10) then mimicks a homogeneous Dirichlet boundary condition; on the contrary, $h_{\Gamma_D, \varepsilon}$ vanishes “well inside” $\Gamma$, and so this condition boils down to a homogeneous Neumann boundary condition, see Figure 4. The key feature of this approximation procedure is that, since $\Omega$ and $h_{\Gamma_D, \varepsilon}$ are smooth, the classical elliptic regularity theory implies that the solution $u_{\Gamma_D, \varepsilon}$ to (3.10) is smooth, contrary to $u_{\Gamma_D}$, see e.g. Chap. 9 in [64] or [65].

We then replace the exact shape functional $J(\Gamma_D)$ in (3.8) by the following approximate counterpart:

$$ J_{\varepsilon}(\Gamma_D) = \int_{\Omega} j(u_{\Gamma_D, \varepsilon}) \, dx. \qquad (3.13) $$

The shape derivative of this regularized functional can now be calculated by standard means.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Approximation of a “sharp” transition between homogeneous Dirichlet and homogeneous Neumann conditions by a smooth Robin boundary condition in Section 3.3.</figcaption></figure>

**Proposition 3.11.** The functional $J_{\varepsilon}(\Gamma_D)$ in (3.13) is shape differentiable at $\Gamma_D$ when deformations $\theta$ are taken in the set $\Theta_{\Gamma_D}$ in (3.9), and its shape derivative reads:

$$ \forall \theta \in \Theta_{\Gamma_D}, \quad J'_{\varepsilon}(\Gamma_D)(\theta) = -\frac{1}{\varepsilon^2} \int_{\Gamma \cup \Gamma_D} h'\left(\frac{d_{\Gamma_D}^{\partial\Omega}(x)}{\varepsilon}\right) \theta(\pi_{\Sigma_D}(x)) \cdot n_{\Sigma_D}(\pi_{\Sigma_D}(x)) u_{\Gamma_D, \varepsilon}(x) p_{\Gamma_D, \varepsilon}(x) \,ds(x), \quad (3.14) $$

where the adjoint state $p_{\Gamma_D, \varepsilon}$ is the unique solution in $H^1(\Omega)$ to the following boundary value problem:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(\gamma \nabla p_{\Gamma_D, \varepsilon}) = -j'(u_{\Gamma_D, \varepsilon}) & \text{in } \Omega, \\ \gamma \frac{\partial p_{\Gamma_D, \varepsilon}}{\partial n} + h_{\Gamma_D, \varepsilon} p_{\Gamma_D, \varepsilon} = 0 & \text{on } \Gamma \cup \Gamma_D, \\ \gamma \frac{\partial p_{\Gamma_D, \varepsilon}}{\partial n} = 0 & \text{on } \Gamma_N. \end{array} \right. \quad (3.15) $$

*Sketch of the proof.* This result is proved rigorously in [35]. Here, we outline a formal and relatively simple calculation which can be readily adapted to various situations. The latter is based on the formal method of Céa, whose practice is legitimate in the treatment of the approximate functional $J_{\varepsilon}(\Gamma_D)$, since all the objects featured in its expression (the region $\Gamma_D$, the state function $u_{\Gamma_D, \varepsilon}$, etc.) are smooth, see again [39, 40, 61].

For a given region $\Gamma_D \subset \partial\Omega$, we introduce the Lagrangian $\mathcal{L}(\Gamma_D, \cdot, \cdot) : H^1(\mathbb{R}^d) \times H^1(\mathbb{R}^d) \to \mathbb{R}$ defined by:

$$ \mathcal{L}(\Gamma_D, u, p) = \int_{\Omega} j(u) \,dx + \int_{\Omega} \gamma \nabla u \cdot \nabla p \,dx + \int_{\partial\Omega} h_{\Gamma_D, \varepsilon} u p \,ds - \int_{\Omega} f p \,dx - \int_{\Gamma_N} g p \,ds. $$

Intuitively, the latter is obtained by replacing the state function $u_{\Gamma_D, \varepsilon}$ by a dummy function $u$ in the definition (3.13) of $J_{\varepsilon}(\Gamma_D)$, and then enforcing the constraint that $u = u_{\Gamma_D, \varepsilon}$ with a Lagrange multiplier $p$. By construction, it holds:

$$ \forall p \in H^1(\mathbb{R}^d), \quad J_{\varepsilon}(\Gamma_D) = \mathcal{L}(\Gamma_D, u_{\Gamma_D, \varepsilon}, p). \quad (3.16) $$
For a fixed region $\Gamma_D$, we now search for the saddle points $(u, p) \in H^1(\mathbb{R}^d) \times H^1(\mathbb{R}^d)$ of the function $\mathcal{L}(\Gamma_D, \cdot, \cdot)$, and to this end, we calculate its partial derivatives:

* The partial derivative $\frac{\partial \mathcal{L}}{\partial p}(\Gamma_D, u, p)$ equals:

$$ \forall \hat{p} \in H^1(\mathbb{R}^d), \quad \frac{\partial \mathcal{L}}{\partial p}(\Gamma_D, u, p)(\hat{p}) = \int_{\Omega} \gamma \nabla u \cdot \nabla \hat{p} \, dx + \int_{\partial \Omega} h_{\Gamma_D, \varepsilon} u \hat{p} \, ds - \int_{\Omega} f \hat{p} \, dx - \int_{\Gamma_N} g \hat{p} \, ds. $$

* The partial derivative $\frac{\partial \mathcal{L}}{\partial u}(\Gamma_D, u, p)$ reads:

$$ \forall \hat{u} \in H^1(\mathbb{R}^d), \quad \frac{\partial \mathcal{L}}{\partial u}(\Gamma_D, u, p)(\hat{u}) = \int_{\Omega} j'(u) \hat{u} \, dx + \int_{\Omega} \gamma \nabla \hat{u} \cdot \nabla p \, dx + \int_{\partial \Omega} h_{\Gamma_D, \varepsilon} \hat{u} p \, ds. $$

The condition that the derivative $\frac{\partial \mathcal{L}}{\partial p}(\Gamma_D, u, p)(\hat{p})$ vanish for arbitrary variations $\hat{p} \in H^1(\mathbb{R}^d)$ implies that $u = u_{\Gamma_D, \varepsilon}$, the solution to (3.10). Likewise, since $\frac{\partial \mathcal{L}}{\partial u}(\Gamma_D, u, p)(\hat{u})$ should vanish for arbitrary variations $\hat{u} \in H^1(\mathbb{R}^d)$, we recognize that the function $p$ is the unique solution $p_{\Gamma_D, \varepsilon}$ to the adjoint problem (3.15).

Now returning to (3.16) and taking derivatives with respect to $\Gamma_D$ in a direction $\theta \in \Theta_{\Gamma_D}$, we obtain from the chain rule that, for an arbitrary function $p \in H^1(\mathbb{R}^d)$:

$$ J'_{\varepsilon}(\Gamma_D)(\theta) = \frac{\partial \mathcal{L}}{\partial \Gamma_D}(\Gamma_D, u_{\Gamma_D, \varepsilon}, p)(\theta) + \frac{\partial \mathcal{L}}{\partial u}(\Gamma_D, u_{\Gamma_D, \varepsilon}, p)(u'_{\Gamma_D, \varepsilon}(\theta)), \quad (3.17) $$

where $u'_{\Gamma_D, \varepsilon}(\theta)$ is the “Eulerian” derivative of $\Gamma_D \mapsto u_{\Gamma_D, \varepsilon}$, i.e. for each $x \in \Omega$, $u'_{\Gamma_D, \varepsilon}(\theta)(x)$ is the derivative of the mapping $\theta \mapsto u_{(\Gamma_D)\theta, \varepsilon}(x)$, defined from a neighborhood of 0 in $\Theta_{\Gamma_D}$ into $\mathbb{R}$, see Remark 3.10. In order to eliminate the second term from the above formula, involving this “difficult” quantity, we make the particular choice $p = p_{\Gamma_D, \varepsilon}$ in (3.17). This yields:

$$ J'_{\varepsilon}(\Gamma_D)(\theta) = \frac{\partial \mathcal{L}}{\partial \Gamma_D}(\Gamma_D, u_{\Gamma_D, \varepsilon}, p_{\Gamma_D, \varepsilon})(\theta), $$

featuring the partial derivative of $\mathcal{L}$ with respect to its explicit dependence on the region $\Gamma_D$. The latter is now easily calculated thanks to Proposition 2.3 and 3.9 and Remark 3.10, thus leading to the desired result. $\square$

The above replacement procedure of the quantities $J(\Gamma_D)$, $u_{\Gamma_D}$, attached to a sharp transition between homogeneous Dirichlet and Neumann boundary conditions by the smoothed counterparts $J_{\varepsilon}(\Gamma_D)$, $u_{\Gamma_D, \varepsilon}$ has been proved to be consistent in our previous work [35]. Loosely speaking, under particular assumptions about the geometry of the region $\Sigma_D$, the approximate state $u_{\Gamma_D, \varepsilon}$, the approximate objective function $J_{\varepsilon}(\Gamma_D)$ and its shape derivative $J'_{\varepsilon}(\Gamma_D)(\theta)$ converge to their exact counterparts $u_{\Gamma_D}$, $J(\Gamma_D)$ and $J'(\Gamma_D)(\theta)$ as the regularization parameter $\varepsilon$ vanishes.

### 3.4. A practical approximate formula for the regularized shape derivative

However tractable, the numerical evaluation of the approximate shape derivative $J'_{\varepsilon}(\Gamma_D)(\theta)$ in (3.14) still proves difficult. Firstly, it requires the calculation of the projection $\pi_{\Sigma_D}(y)$ onto the boundary $\Sigma_D$ of $\Gamma_D$ for each point $y$ in a tubular neighborhood of $\Sigma_D$ – an issue which is already quite difficult to address with satisfying accuracy in the classical setting of a domain in the Euclidean space, see [66]. Secondly, the formula (3.14) does not conform to the structure (2.4), and thus it does not directly provide a descent direction for $J_{\varepsilon}(\Gamma_D)$. Both issues can be alleviated thanks to a subsequent approximation of (3.14) inspired by our previous work [67], that we now describe.
Since the derivative $h'$ of the smoothing profile $h$ has compact support inside $[-1, 1]$, the integrand of the formula (3.14) has compact support inside the tubular neighborhood $U_{\epsilon} = \{x \in \partial\Omega \text{ s.t. } d^{\partial\Omega}(x, \Sigma) < \epsilon\}$, and Proposition 3.6 yields:

$$ 
\begin{aligned}
J'_{\epsilon}(\Gamma_D)(\theta) &= -\frac{1}{\epsilon^2} \int_{U_{\epsilon}} h'\left(\frac{d^{\partial\Omega}(y)}{\epsilon}\right) \theta(\pi_{\Sigma_D}(y)) \cdot n_{\Sigma_D}(\pi_{\Sigma_D}(y)) u_{\Gamma_D, \epsilon}(y) p_{\Gamma_D, \epsilon}(y) \mathrm{d}s(y) \\
&= -\frac{1}{\epsilon^2} \int_{-\epsilon}^{\epsilon} \int_{\Sigma_D} h'\left(\frac{t}{\epsilon}\right) |\det\nabla F(x, t)| \theta(x) \cdot n_{\Sigma_D}(x) u_{\Gamma_D, \epsilon}(F(x, t)) p_{\Gamma_D, \epsilon}(F(x, t)) \mathrm{d}\ell(x) \mathrm{d}t.
\end{aligned}
$$

Yet another change of variables in the outermost integral yields:

$$ J'_{\epsilon}(\Gamma_D)(\theta) = -\frac{1}{\epsilon} \int_{-1}^{1} \int_{\Sigma_D} h'(t) |\det\nabla F(x, \epsilon t)| \theta(x) \cdot n_{\Sigma_D}(x) u_{\Gamma_D, \epsilon}(F(x, \epsilon t)) p_{\Gamma_D, \epsilon}(F(x, \epsilon t)) \mathrm{d}\ell(x) \mathrm{d}t. $$

Since $u_{\Gamma_D, \epsilon}$ and $p_{\Gamma_D, \epsilon}$ are smooth, it holds, for $\epsilon > 0$ small enough:

$$ \forall x \in \Sigma_D, \ t \in (-1, 1), \quad u_{\Gamma_D, \epsilon}(F(x, \epsilon t)) \approx u_{\Gamma_D, \epsilon}(F(x, 0)) = u_{\Gamma_D, \epsilon}(x), $$

and likewise for $p_{\Gamma_D, \epsilon}$. Hence, we obtain:

$$ 
\begin{aligned}
J'_{\epsilon}(\Gamma_D)(\theta) &\approx -\frac{1}{\epsilon} \int_{-1}^{1} \int_{\Sigma_D} h'(t) |\det\nabla F(x, 0)| \theta(x) \cdot n_{\Sigma_D}(x) u_{\Gamma_D, \epsilon}(x) p_{\Gamma_D, \epsilon}(x) \mathrm{d}\ell(x) \mathrm{d}t \\
&= -\frac{1}{\epsilon} \int_{-1}^{1} \int_{\Sigma_D} h'(t) \theta(x) \cdot n_{\Sigma_D}(x) u_{\Gamma_D, \epsilon}(x) p_{\Gamma_D, \epsilon}(x) \mathrm{d}\ell(x) \mathrm{d}t,
\end{aligned}
$$

where we have used Remark 3.7 to pass from the first line to the second. Exploiting once more the properties (3.12) of the function $h$, we arrive at:

$$ J'_{\epsilon}(\Gamma_D)(\theta) \approx -\frac{1}{\epsilon} \int_{\Sigma_D} \theta(x) \cdot n_{\Sigma_D}(x) u_{\Gamma_D, \epsilon}(x) p_{\Gamma_D, \epsilon}(x) \mathrm{d}\ell(x), $$

an approximate formula which is simple to evaluate in practice and fulfills the desirable structure (2.4).

### 3.5. Shape derivative of a function depending on the inhomogeneous Neumann region

In this section, we briefly discuss the situation where the region $G \subset \partial\Omega$ to be optimized is that $\Gamma_N$ supporting the inhomogeneous Neumann boundary condition in the problem (3.2). Accordingly, we set $\Sigma_N = \partial\Gamma_N$, and we now denote by $u_{\Gamma_N}$ the solution to (3.2); we then show how to calculate the shape derivative of the model objective function

$$ J(\Gamma_N) = \int_{\Omega} j(u_{\Gamma_N}) \mathrm{d}x, $$

where $j: \mathbb{R} \to \mathbb{R}$ again satisfies (3.4). In this context, we consider the following subset of $W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)$ for deformations $\theta$ in the practice (2.1) and (2.2) of the method of Hadamard:

$$ \Theta_{\Gamma_N} := \{ \theta : \mathbb{R}^d \to \mathbb{R}^d \text{ is smooth, } \theta \cdot n = 0 \text{ on } \partial\Omega \text{ and } \theta = 0 \text{ on } \Gamma_D \}. \quad (3.18) $$
This situation turns out to be much easier to handle than that considered in Sections 3.3 and 3.4. The potential $u$ in (3.2) is indeed less singular near the optimized boundary $\Sigma_N$, now marking the transition between homogeneous and inhomogeneous Neumann boundary conditions, than near the boundary $\Sigma_D$ between homogeneous Dirichlet and Neumann conditions, that was optimized in the previous developments. In particular, the calculation of the shape derivative $J'(\Gamma_N)(\theta)$ does not require an approximation procedure in the spirit of that developed in Sections 3.3 and 3.4, and its “exact” shape derivative can be used as is in practice.

The main result in the present situation is the following; we refer to our previous work [35] for a rigorous proof. Interestingly, the correct expression for the derivative can be obtained by an application of the formal Céa's method exemplified in the proof of Proposition 3.11.

**Theorem 3.12.** The shape functional $J(\Gamma_N)$ is shape differentiable when deformations are considered in the set $\Theta_{\Gamma_N}$ in (3.18), and its shape derivative reads:

$$ \forall \theta \in \Theta_{\Gamma_N}, \quad J'(\Gamma_N)(\theta) = - \int_{\Sigma_N} g p_{\Gamma_N} \theta \cdot n_{\Sigma_N} \, d\ell, $$

where the adjoint state $p_{\Gamma_N}$ is the unique solution in $H^1(\Omega)$ to the following boundary value problem:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(\gamma \nabla p_{\Gamma_N}) = -j'(u_{\Gamma_N}) & \text{in } \Omega, \\ p_{\Gamma_N} = 0 & \text{on } \Gamma_D, \\ \gamma \frac{\partial p_{\Gamma_N}}{\partial n} = 0 & \text{on } \Gamma \cup \Gamma_N. \end{array} \right. $$

**Remark 3.13.** Multiple variations of the above configurations can be treated with the same techniques. For instance, one could consider the case where inhomogeneous Dirichlet boundary conditions are imposed on the region $G = \Gamma_D$ to be optimized, of the form $u = u_{in}$, where $u_{in}$ is a given, smooth function, see Section 7.3 for a related numerical example. On a different note, one could deal with a function $J(G)$ depending on a region $G$ bearing Robin boundary conditions, see Section 5 for a related study in the physical context of acoustics, and Section 7.4 for a numerical example.

## 4. TOPOLOGICAL PERTURBATIONS OF REGIONS SUPPORTING BOUNDARY CONDITIONS IN THE SETTING OF THE CONDUCTIVITY EQUATION

As in Section 3.1, we consider the state equation (3.2) related to the conductivity equation. Again, in the main part of this section, the region $G \subset \partial\Omega$ to be optimized is $\Gamma_D$, that which supports Dirichlet boundary conditions. We then consider the calculation of the topological derivative of the model functional (3.3), in the sense of Definition 2.1.

As we have mentioned, this task leverages techniques from asymptotic analysis. For short, and out of consistency with the classical notation and jargon of that field, the unperturbed, “background” potential is denoted by $u_0$ throughout this section. It is characterized by the boundary value problem (3.2), that we recall for convenience:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(\gamma \nabla u_0) = f & \text{in } \Omega, \\ u_0 = 0 & \text{on } \Gamma_D, \\ \gamma \frac{\partial u_0}{\partial n} = g & \text{on } \Gamma_N, \\ \gamma \frac{\partial u_0}{\partial n} = 0 & \text{on } \Gamma. \end{array} \right. \qquad (4.1) $$
We fix $x_0 \in \Gamma$, and consider the “perturbed” potential $u_\varepsilon$, which is the unique $H^1(\Omega)$ solution to:

$$ egin{cases} 
-\operatorname{div}(\gamma \nabla u_\varepsilon) = f & \text{in } \Omega, \\
u_\varepsilon = 0 & \text{on } \Gamma_D \cup \omega_{x_0, \varepsilon}, \\
\gamma \frac{\partial u_\varepsilon}{\partial n} = g & \text{on } \Gamma_N, \\
\gamma \frac{\partial u_\varepsilon}{\partial n} = 0 & \text{on } \Gamma \setminus \overline{\omega_{x_0, \varepsilon}}, 
\end{cases} 
\quad (4.2) $$

where $\omega_{x_0, \varepsilon}$ is the surface disk defined in (2.8). Our goal is to determine the asymptotic behavior of the potential $u_\varepsilon$ as $\varepsilon \to 0$, as well as that of the quantity of interest

$$ J((\Gamma_D)_{x_0, \varepsilon}) = \int_{\Omega} j(u_{\varepsilon}) \, dx. \quad (4.3) $$

After a brief reminder about the suitable functional spaces involved in the treatment of this problem and the Green's function of the background problem (4.1) in Section 4.1, we recall in Section 4.2 the result from [33] giving the asymptotic expansion of $u_\varepsilon$. We provide a formal sketch of the (difficult) proof in a representative particular case, which will be adapted to address other situations, as well as more intricate physical contexts. We then take advantage of these results in Section 4.3 to calculate the asymptotic expansion of $J((\Gamma_D)_{x_0, \varepsilon})$ and therefore identify the topological derivative of the function $J(\Gamma_D)$. We eventually discuss several extensions of these results to different types of changes of boundary conditions in Section 4.4.

## 4.1. Preliminaries

In this section, we introduce the necessary background material from functional analysis and potential theory. Throughout, $\Omega$ stands for a smooth bounded domain in $\mathbb{R}^d$.

### 4.1.1. Sobolev spaces on a subregion of the boundary of a domain

To a large extent, the calculations of the topological derivatives at stake in this article bring into play functions defined on (a region of) the boundary $\partial\Omega$ of the ambient domain $\Omega$, as the traces of solutions to boundary value problems posed on $\Omega$. The definitions of the suitable energy spaces for these are recalled in the present section; we refer for instance to [59, 68] about these issues.

Let us first consider functions defined on the whole boundary $\partial\Omega$. For any real number $0 < s < 1$, the fractional Sobolev space $H^s(\partial\Omega)$ is defined by:

$$ H^s(\partial\Omega) = \{u \in L^2(\partial\Omega) \text{ s.t. } \|u\|_{H^s(\partial\Omega)} < \infty\}, $$ 
where $\|u\|_{H^s(\partial\Omega)} := \|u\|_{L^2(\partial\Omega)}^2 + \int_{\partial\Omega} \int_{\partial\Omega} \frac{|u(x) - u(y)|^2}{|x - y|^{d-1+2s}} \, ds(x)ds(y). $$

By convention, $H^0(\partial\Omega) = L^2(\partial\Omega)$, and for $0 < -s < 1$, $H^{-s}(\partial\Omega)$ is the topological dual of $H^s(\partial\Omega)$.

Let now $G$ be a Lipschitz open subset of $\partial\Omega$; we shall use two types of fractional Sobolev spaces of functions on $G$:

* For any $0 < s < 1$, we denote by $H^s(G)$ the space of restrictions to $G$ of functions in $H^s(\partial\Omega)$:

$$ H^s(G) = \{U|_G \text{ for some } U \in H^s(\partial\Omega)\}, \text{ equipped with the quotient norm} $$ 
$$ \|u\|_{H^s(G)} := \inf \{ \|U\|_{H^s(\partial\Omega)}, U \in H^s(\partial\Omega), U|_G = u \}. $$
* For any $0 < s < 1$, we denote by $\tilde{H}^s(G)$ the subspace of $H^s(\partial\Omega)$ made of functions with compact support in $\bar{G}$. It is equipped with the norm $||u||_{\tilde{H}^s(G)} = ||u||_{H^s(\partial\Omega)}$ induced by that of $H^s(\partial\Omega)$. Equivalently, $H^s(G)$ is the space of functions $u \in L^2(G)$ whose extension $\tilde{u}$ by 0 to $\partial\Omega$ belongs to $H^s(\partial\Omega)$.

As for the negative versions of these spaces,

* For $0 < s < 1$, $H^{-s}(G)$ is still defined as the space of distributions on $G$ obtained by restriction of a distribution in $H^{-s}(\partial\Omega)$. The space $H^{-s}(G)$ is naturally identified with the dual space of $\tilde{H}^s(G)$ via the following duality:

$$ \forall u \in H^{-s}(G), \quad v \in \tilde{H}^s(G), \quad \langle u, v \rangle_{H^{-s}(G), \tilde{H}^s(G)} := \langle U, \tilde{v} \rangle_{H^{-s}(\partial\Omega), H^s(\partial\Omega)}, $$

where $U$ is any element in $H^{-s}(\partial\Omega)$ such that $U|_G = u$ and $\tilde{v}$ is the extension of $v$ by 0 to the whole hypersurface $\partial\Omega$.

* For $0 < s < 1$, $\tilde{H}^{-s}(G)$ is again the subspace of $H^{-s}(\partial\Omega)$ made of elements with compact support in $\bar{G}$. This space is naturally identified with the dual of $H^s(G)$ via the following duality:

$$ \forall u \in \tilde{H}^{-s}(G), \quad v \in H^s(G), \quad \langle u, v \rangle_{\tilde{H}^{-s}(G), H^s(G)} := \langle \tilde{u}, V \rangle_{H^{-s}(\partial\Omega), H^s(\partial\Omega)}, $$

where $\tilde{u}$ is the extension of $u$ by 0 to $\partial\Omega$ and $V$ is any element in $H^s(\partial\Omega)$ such that $V|_G = v$.

### 4.1.2. The Green's function for the background problem

Let us recall the expression of the fundamental solution $\Gamma(x, y)$ of the operator $-\Delta$ in free space:

$$ \forall x, y \in \mathbb{R}^d, x \neq y, \quad \Gamma(x, y) = \begin{cases} -\frac{1}{2\pi} \log|x-y| & \text{if } d=2, \\ \frac{1}{(d-2)\alpha_d} |x-y|^{2-d} & \text{if } d \ge 3, \end{cases} \quad (4.4) $$

where $\alpha_d$ is the area of the unit sphere $S^{d-1} \subset \mathbb{R}^d$. For each $x \in \mathbb{R}^d$, the function $\Gamma(x, \cdot)$ satisfies:

$$ -\Delta_y \Gamma(x, y) = \delta_{y=x} \quad \text{in the sense of distributions on } \mathbb{R}^d, $$

where $\delta_{y,x}$ is the Dirac distribution at $x$. Here and throughout the article, the subscript $y$ indicates that the differential operator applies to the $y$ variable only. We set:

$$ \forall x, y \in \mathbb{R}^d, x \neq y, \quad \Gamma(x, y) := \Gamma(x - y). $$

The Green's function $N(x, y)$ for the background equation (4.1) is next defined by the following relation:

$$ \text{For all } x \in \Omega, y \mapsto N(x, y) \text{ satisfies} \quad \left\{ \begin{array}{ll} -\text{div}_y(\gamma(y)\nabla_y N(x, y)) = \delta_{y=x} & \text{in } \Omega, \\ N(x, y) = 0 & \text{for } y \in \Gamma_D, \\ \gamma(y) \frac{\partial N}{\partial n_y}(x, y) = 0 & \text{for } y \in \Gamma \cup \Gamma_N. \end{array} \right. \quad (4.5) $$

**Remark 4.1.** The function $N(x, y)$ is symmetric in its arguments (see [69] for a proof) and it is related to the fundamental solution $\Gamma(x, y)$ of the operator $-\Delta$ as follows:

$$ N(x, y) = \frac{1}{\gamma(x)} \Gamma(x, y) + R(x, y), \quad (4.6) $$
where the correction term $R(x, y)$ is the solution to:

$$ egin{cases} -\operatorname{div}_y(\gamma(y)\nabla_y R(x, y)) = \frac{1}{\gamma(y)}\nabla\gamma(y) \cdot \nabla_y\Gamma(x, y) & \text{in } \Omega, \\ R(x, y) = -\frac{1}{\gamma(y)}\Gamma(x, y) & \text{for } y \in \Gamma_D, \\ \gamma(y)\frac{\partial R}{\partial n_y}(x, y) = \frac{\gamma(y)}{\gamma(x)}\frac{\partial \Gamma}{\partial n_y}(x, y) & \text{for } y \in \Gamma \cup \Gamma_N. \end{cases} \quad (4.7) $$

Since $\gamma$ and $\Omega$ are smooth, for every open subset $U \in \mathbb{R}^d \setminus (\Sigma \cup \{x\})$, the function $R(x, \cdot)$ is of class $C^\infty$ on $\bar{\Omega} \cap U$, see [64, 65]. We refer to e.g. [70, 71] for further details about the construction of $N(x, y)$.

The key property of the Green's function $N(x, y)$ is the following relation, which holds (at least) for functions $\varphi \in C^1(\bar{\Omega})$ such that $\varphi = 0$ on $\Gamma_D$:

$$ \varphi(x) = \int_{\Omega} \gamma(y) \nabla_y N(x, y) \cdot \nabla \varphi(y) \, dy, \quad x \in \Omega. \quad (4.8) $$

In particular, one may integrate by parts to express the solution $u_0$ to (4.1) in terms of $N(x, y)$ as:

$$ u_0(x) = \int_{\Omega} f(y) N(x, y) \, dy. \quad (4.9) $$

### 4.1.3. The Neumann function of the lower half-space

In the following, we shall need another Green's function $L_\gamma(x, y)$, associated to the version of (4.1) featuring a constant conductivity $\gamma > 0$, posed on the lower half-space $H$ (see (2.6)), and equipped with homogeneous Neumann boundary conditions on $\partial H$. For $x \in H$, $y \mapsto L_\gamma(x, y)$ satisfies:

$$ \left\{ \begin{array}{ll} -\operatorname{div}_y(\gamma \nabla_y L_\gamma(x, y)) = \delta_{y=x} & \text{in } H, \\ \gamma \frac{\partial L_\gamma}{\partial n_y}(x, y) = 0 & \text{for } y \in \partial H. \end{array} \right. \quad (4.10) $$

Such a function can be constructed thanks to the so-called method of images, see e.g. [72]; precisely:

$$ L_{\gamma}(x, y) = \frac{1}{\gamma} (\Gamma(x - y) + \Gamma(x + y)). \quad (4.11) $$

It is indeed straightforward to see that the above function satisfies (4.10). Similarly to the formula (4.8) involving $N(x, y)$, for $\varphi \in C_c^1(\bar{H})$, the following relation holds:

$$ \varphi(x) = \int_H \gamma \nabla_y L(x, y) \cdot \nabla \varphi(y) \, dy, \quad x \in H. \quad (4.12) $$

### 4.1.4. The single layer potential operator

Let us recall the definition of the single layer potential operator $S_\Omega$, which represents the potential generated in $\mathbb{R}^d$ by a density of charges $\varphi : \partial\Omega \to \mathbb{R}$ on $\partial\Omega$.

**Definition 4.2.** The single layer potential associated to a smooth density $\varphi \in C^\infty(\partial\Omega)$ is the function defined by:

$$ S_{\Omega}\varphi(x) = \int_{\partial\Omega} \Gamma(x, y)\varphi(y) \, ds(y), \quad x \in \mathbb{R}^d \setminus \partial\Omega, $$

where $\Gamma(x, y)$ is the fundamental solution of the operator $-\Delta$ in free space, see (4.4).
It is well-known that the single layer potential is continuous across $∂Ω$, and that it induces an operator $S_{∂Ω} : C^∞(∂Ω) → C^∞(∂Ω)$ defined by:

$$S_{∂Ω}φ(x) = \int_{∂Ω} Γ(x, y)φ(y) \, ds(y), \quad x ∈ ∂Ω.$$

The next proposition gathers a few properties of this mapping.

**Proposition 4.3.** The following facts hold true:

(i) The mapping $S_{∂Ω}$ has an extension as a bounded mapping from $H^{-1/2}(∂Ω)$ into $H^{1/2}(∂Ω)$.

(ii) Let $G$ be a Lipschitz open subset of $∂Ω$; then $S_{∂Ω}$ induces a bounded operator $S_G : \tilde{H}^{-1/2}(G) → H^{1/2}(G)$ via the formula:

$$∀φ ∈ \tilde{H}^{-1/2}(G), \quad S_Gφ = (S_{∂Ω}φ)|_G,$$

where $\tilde{φ} ∈ H^{-1/2}(∂Ω)$ is the extension by 0 to $∂Ω$ of an element $φ ∈ \tilde{H}^{-1/2}(G)$.

(iii) If $d ≥ 3$, the mapping $S_G : \tilde{H}^{-1/2}(G) → H^{1/2}(G)$ is invertible.

The proof of (i) is found in [68], while (ii) follows almost immediately from the definitions of the functional spaces $\tilde{H}^{-1/2}(G)$ and $H^{1/2}(G)$. The last point (iii) is more subtle, and it is proved in [33]. Note that when $G$ is the unit disk $D_1$, this result also holds when $d = 2$. This fact is needed in the rigorous proof of the forthcoming Theorem 4.4, which is conducted in [33], but not in the formal calculation method proposed in the next Section 4.2.

## 4.2. Asymptotic expansion of the voltage potential $u_ε$

This section is devoted to the asymptotic expansion of the perturbed potential $u_ε$, solution to (4.2), in the limit $ε → 0$, where the support $ω_{x_0,ε}$ of the replacement of homogeneous Neumann boundary conditions by homogeneous Dirichlet boundary conditions vanishes. The main result in this setting is the following asymptotic expansion, whose proof is rigorously detailed in [33]. We provide a formal sketch of the latter, under additional technical assumptions, which lends itself to generalizations.

**Theorem 4.4.** Let $x_0$ be a given point in $Γ$; then for any point $x ∈ \bar{Ω} \setminus (Σ_D ∪ \{x_0\})$, the following asymptotic expansion holds:

$$u_{\varepsilon}(x) = u_{0}(x) - \frac{\pi}{|\log \varepsilon|} \gamma(x_{0}) u_{0}(x_{0}) N(x, x_{0}) + o\left(\frac{1}{|\log \varepsilon|}\right) \quad \text{if } d = 2, \qquad (4.13)$$

and

$$u_{\varepsilon}(x) = u_{0}(x) - 4\varepsilon \gamma(x_{0}) u_{0}(x_{0}) N(x, x_{0}) + o(\varepsilon) \quad \text{if } d = 3. \qquad (4.14)$$

*Sketch of proof.* Without loss of generality, we assume that the center $x_0$ of the surface disk where the substitution of boundary conditions occurs is the origin 0 and that the unit normal vector $n(0)$ to $∂Ω$ at 0 coincides with the $d^{th}$ coordinate vector $e_d$; we denote $ω_ε := ω_{x_0,ε}$ throughout the proof. Our formal argument hinges on the simplifying assumptions that $γ$ is constant, and that the boundary $∂Ω$ is completely flat around 0, that is:

$$∂Ω coincides with ∂H in a neighborhood O of 0 in R^d. \qquad (4.15)$$

This implies in particular that $ω_ε$ coincides with $D_ε$, see Figure 5.
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Illustration of the flatness (4.15) of $\partial\Omega$ near $x_0 = 0$ assumed in the proof of Theorem 4.4.</figcaption></figure>

Let us introduce the error $r_{\varepsilon} = u_{\varepsilon} - u_0 \in H^1(\Omega)$ between the perturbed and background potentials; this function satisfies the boundary value problem:

$$ egin{cases} -\operatorname{div}(\gamma \nabla r_{\varepsilon}) = 0 & \text{in } \Omega, \\ r_{\varepsilon} = 0 & \text{on } \Gamma_D, \\ r_{\varepsilon} = -u_0 & \text{on } \omega_{\varepsilon}, \\ \gamma \frac{\partial r_{\varepsilon}}{\partial n} = 0 & \text{on } \Gamma_N \cup (\Gamma \setminus \bar{\omega}_{\varepsilon}). \end{cases} \quad (4.16) $$

Standard a priori estimates show this error tends to 0:

$$ r_{\varepsilon} \xrightarrow{\varepsilon \to 0} 0 \text{ in } H^1(\Omega), \quad (4.17) $$

see again [33]. We now proceed in four steps to calculate an expansion of $r_{\varepsilon}(x)$ as $\varepsilon \to 0$.

**Step 1:** We construct a representation formula for the values of $r_{\varepsilon}$ “far” from 0 in terms of its values inside the vanishing region $\omega_{\varepsilon}$ and the Green’s function $N(x, y)$.

This task starts from the integral expression of $r_{\varepsilon}$ based on the Green’s function $N(x, y)$ in (4.5). For any point $x \in \Omega$, it holds:

$$ r_{\varepsilon}(x) = - \int_{\Omega} \operatorname{div}_{y}(\gamma \nabla_{y} N(x, y)) r_{\varepsilon}(y) \, dy. $$

Two successive integrations by parts in the above expression yield:

$$ r_{\varepsilon}(x) = - \int_{\partial\Omega} \gamma \frac{\partial N}{\partial n_y}(x, y) r_{\varepsilon}(y) \, ds(y) + \int_{\Omega} \gamma \nabla_y N(x, y) \cdot \nabla r_{\varepsilon}(y) \, dy \quad (4.18) $$

$$ = - \int_{\partial\Omega} \gamma \frac{\partial N}{\partial n_y}(x, y) r_{\varepsilon}(y) \, ds(y) + \int_{\partial\Omega} \gamma \frac{\partial r_{\varepsilon}}{\partial n}(y) N(x, y) \, ds(y), \quad (4.19) $$

where we have used the equation (4.16) satisfied by $r_{\varepsilon}$ for passing from the first line to the second one. Now, taking into account the boundary conditions satisfied by $r_{\varepsilon}$ and $y \mapsto N(x, y)$, the first integral in the above right-hand side vanishes. Likewise, the boundary conditions in (4.5) and (4.16) imply that the integrand of the
second term vanishes on $\Gamma_D$, $\Gamma_N$ and $\Gamma \setminus \bar{\omega}_\varepsilon$, which leaves us with the expression:

$$r_{\varepsilon}(x) = \int_{\omega_{\varepsilon}} \gamma \frac{\partial r_{\varepsilon}}{\partial n}(y) N(x, y) \, ds(y).$$

A change of variables in the above integral now yields:

$$r_{\varepsilon}(x) = \int_{\mathbb{D}_1} \varphi_{\varepsilon}(z) N(x, \varepsilon z) \, ds(z), \qquad (4.20)$$

where we have introduced the function $\varphi_\varepsilon(z) := \varepsilon^{d-1} (\gamma \frac{\partial r_\varepsilon}{\partial n})(\varepsilon z) \in \tilde{H}^{-1/2}(\mathbb{D}_1)$. This is the desired representation formula for $r_\varepsilon(x)$.

**Step 2: We characterize $\varphi_\varepsilon$ as the solution to an integral equation.**

To achieve this, we repeat the derivation of the first step, except that we replace the true (non-explicit) Green's function $N(x, y)$ of the problem (4.1) by the approximate, but explicit, function $L_\gamma(x, y)$ introduced in Section 4.1.3, which captures the behavior of (4.4) near 0. Again, for any point $x \in \Omega$, two consecutive integrations by parts yield:

$$r_{\varepsilon}(x) = - \int_{\Omega} \mathrm{div}_y (\gamma \nabla_y L_{\gamma}(x, y)) r_{\varepsilon}(y) \, dy \qquad (4.21)$$

$$= - \int_{\partial \Omega} \gamma \frac{\partial L_{\gamma}}{\partial n_y}(x, y) r_{\varepsilon}(y) \, ds(y) + \int_{\Omega} \gamma \nabla_y L_{\gamma}(x, y) \cdot \nabla r_{\varepsilon}(y) \, dy \qquad (4.22)$$

$$= - \int_{\partial \Omega} \gamma \frac{\partial L_{\gamma}}{\partial n_y}(x, y) r_{\varepsilon}(y) \, ds(y) + \int_{\partial \Omega} \gamma \frac{\partial r_{\varepsilon}}{\partial n}(y) L_{\gamma}(x, y) \, ds(y). \qquad (4.23)$$

Now using the fact that $\partial\Omega$ coincides with $\partial H$ in a neighborhood $\mathcal{O}$ of 0 in $\mathbb{R}^d$, and considering the boundary conditions satisfied by $r_\varepsilon$, we obtain that:

$$r_{\varepsilon}(x) = \int_{\omega_{\varepsilon}} \gamma \frac{\partial r_{\varepsilon}}{\partial n}(y) L_{\gamma}(x, y) \, ds(y) + K r_{\varepsilon}(x), \text{ where}$$

$$Kr_{\varepsilon}(x) := - \int_{(\Gamma_N \cup \Gamma) \setminus \bar{\Omega}} \gamma \frac{\partial L_{\gamma}}{\partial n_y}(x, y) r_{\varepsilon}(y) \, ds(y) + \int_{\Gamma_D} \gamma \frac{\partial r_{\varepsilon}}{\partial n}(y) L_{\gamma}(x, y) \, ds(y).$$

We now change variables in the above integrals and let $x$ approach $\omega_\varepsilon$. Since the single layer potential is continuous across $\partial\Omega$, as expressed in Proposition 4.3, we obtain that:

$$\forall x \in \mathbb{D}_1, \quad -u_0(\varepsilon x) = \int_{\mathbb{D}_1} \varphi_\varepsilon(z) L_\gamma(\varepsilon x, \varepsilon z) \, ds(z) + Kr_\varepsilon(\varepsilon x).$$

Now, since $\varepsilon x$ lies “far” from the support of the integrals featured in $Kr_\varepsilon$, the convergence (4.17) of the error $r_\varepsilon$ implies that the last term in the above right-hand side tends to 0 as $\varepsilon \to 0$. Taking advantage of the explicit form (4.11) of the Green's function $L_\gamma(x, y)$, this eventually leads to the following integral equation for the function $\varphi_\varepsilon \in \tilde{H}^{-1/2}(\mathbb{D}_1)$:

$$\frac{2}{\gamma} \int_{\mathbb{D}_1} \varphi_{\varepsilon}(z) \Gamma(\varepsilon x, \varepsilon z) \, ds(z) = -u_0(0) + o(1), \qquad (4.24)$$

where $o(1)$ is a remainder that converges to 0 in $H^{1/2}(\mathbb{D}_1)$.
Step 3: We analyze the integral equation (4.24) to glean information about $\varphi_{\varepsilon}$.

Let us denote by $S_1 := S_{\mathbb{D}_1}$ the single layer potential attached to $\mathbb{D}_1$, see Proposition 4.3. Our analysis of (4.24) relies on the explicit expression (4.4) of $\Gamma(x, y)$, which depends on the space dimension $d$.

• When $d = 2$, (4.24) rewrites:

$$-\int_{\mathbb{D}_1} \log|\varepsilon x - \varepsilon z|\varphi_{\varepsilon}(z) \, ds(z) = -\pi\gamma u_0(0) + o(1),$$

and so:

$$|\log \varepsilon| \left( \int_{\mathbb{D}_1} \varphi_{\varepsilon}(z) \, ds(z) \right) - S_1 \varphi_{\varepsilon}(z) = -\pi \gamma u_0(0) + o(1).$$

This relation immediately yields:

$$\int_{\mathbb{D}_1} \varphi_{\varepsilon}(z) \, ds(z) = -\frac{1}{|\log \varepsilon|} \pi \gamma u_0(0) + o\left(\frac{1}{|\log \varepsilon|}\right). \quad (4.25)$$

• In the case $d = 3$, (4.24) reads:

$$S_1 \varphi_{\varepsilon}(z) = -\frac{\varepsilon}{2} \gamma u_0(0).$$

Using the invertibility of $S_1: \tilde{H}^{-1/2}(\mathbb{D}_1) \to H^{1/2}(\mathbb{D}_1)$ asserted in Proposition 4.3, we obtain:

$$\varphi_{\varepsilon}(z) = -\frac{\varepsilon}{2} \gamma u_0(0) S_1^{-1} 1(z) + o(\varepsilon),$$

where the function $S_1^{-1} 1 \in \tilde{H}^{-1/2}(\mathbb{D}_1)$ is the so-called equilibrium distribution of the single layer potential $S_1$. From the physical viewpoint, it represents the density of charges on $\mathbb{D}_1$ that induces a constant, unit voltage potential on $\mathbb{D}_1$. Its explicit expression is derived in e.g. [72, 73]:

$$S_1^{-1} 1(z) = \frac{4}{\pi \sqrt{1 - |z|^2}}, \quad \text{and in particular:} \quad \int_{\mathbb{D}_1} S_1^{-1} 1(z) \, ds(z) = 8. \quad (4.26)$$

It thus follows that:

$$\int_{\mathbb{D}_1} \varphi_{\varepsilon}(z) \, ds(z) = -4\varepsilon\gamma u_0(0) + o(\varepsilon). \quad (4.27)$$

Step 4: We pass to the limit in the representation formula (4.20).

We return to this end to the representation formula (4.20) for the value $r_{\varepsilon}(x)$ of the error at a point $x \notin \Sigma_D \cup \{0\}$. A Taylor expansion of the smooth function $y \mapsto N(x, y)$ around $y = 0$ together with the Lebesgue dominated convergence theorem yield:

$$\begin{aligned} r_{\varepsilon}(x) &= \int_{\mathbb{D}_1} \varphi_{\varepsilon}(z) N(x, \varepsilon z) \, ds(z) \\ &= \left( \int_{\mathbb{D}_1} \varphi_{\varepsilon}(z) \, ds(z) \right) (N(x, 0) + o(1)), \end{aligned}$$

and the desired formulas (4.13) and (4.14) follow from the combination of this result with (4.25) and (4.27). $\square$
### 4.3. Asymptotic expansion of a quantity of interest

The asymptotic expansion of the voltage potential $u_ε$ derived in the previous section makes it possible to calculate the sensitivity of $J(\Gamma_D)$ with respect to the addition of a surface disk $\omega_{x_0,ε}$ to the Dirichlet zone $\Gamma_D$.

**Corollary 4.5.** Let $x_0 \in \Gamma$ be given. The perturbed value $J((\Gamma_D)_{x_0,ε})$ in (4.3), accounting for the replacement of the homogeneous Neumann boundary condition on $\omega_{x_0,ε} \subset \Gamma$ by a homogeneous Dirichlet boundary condition, has the following asymptotic expansion:

$$J((\Gamma_D)_{x_0,ε}) = J(\Gamma_D) + \frac{\pi}{|\log \varepsilon|} \gamma(x_0) u_0(x_0) p_0(x_0) + o\left(\frac{1}{|\log \varepsilon|}\right) \text{ if } d=2,$$

and

$$J((\Gamma_D)_{x_0,ε}) = J(\Gamma_D) + 4\varepsilon \gamma(x_0) u_0(x_0) p_0(x_0) + o(\varepsilon) \text{ if } d=3.$$

where $p_0$ is the unique $H^1(\Omega)$ solution to the boundary value problem:

$$\left\{
\begin{array}{ll}
-\mathrm{div}(\gamma \nabla p_0) = -j'(u_0) & \text{in } \Omega, \\
p_0 = 0 & \text{on } \Gamma_D, \\
\gamma \frac{\partial p_0}{\partial n} = 0 & \text{on } \Gamma \cup \Gamma_N.
\end{array}
\right.
\qquad (4.28)$$

*Sketch of the proof.* We assume for simplicity that $d=3$, the case $d=2$ being similar. At first, using the asymptotic expansion of Theorem 4.4, we have:

$$J((\Gamma_D)_{x_0,\varepsilon}) = J(\Gamma_D) + \int_{\Omega} j'(u_0(x)) (-4\varepsilon\gamma(x_0) u_0(x_0) N(x, x_0)) dx + o(\varepsilon).$$

This follows from an application of the Lebesgue dominated convergence theorem whose rigorous justification is detailed in [74]. Now using the representation formula (4.9) for the adjoint state $p_0$ in (4.28), we immediately see that:

$$\begin{align*}
J((\Gamma_D)_{x_0,\varepsilon}) &= J(\Gamma_D) - 4\varepsilon\gamma(x_0) u_0(x_0) \int_{\Omega} j'(u_0(x)) N(x, x_0) dx + o(\varepsilon) \\
&= J(\Gamma_D) + 4\varepsilon\gamma(x_0) u_0(x_0) p_0(x_0) + o(\varepsilon),
\end{align*}$$

which is the desired result.

<figure><img src="image_7.png" /></figure>

**Remark 4.6.** It is possible to replace the homogeneous Dirichlet boundary condition on $\Gamma_D$ in (4.1) by an inhomogeneous Dirichlet boundary condition, i.e. $u_ε = u_{in}$, for some given, smooth function $u_{in} \in C^\infty(\mathbb{R}^d)$. The above calculations can be straightforwardly adapted to this case, and the asymptotic expansion of $u_ε$ becomes:

$$u_{\varepsilon}(x) = u_0(x) + \frac{\pi}{|\log \varepsilon|} \gamma(x_0)(u_{in}(x_0) - u_0(x_0))N(x, x_0) + o\left(\frac{1}{|\log \varepsilon|}\right) \text{ if } d=2,$$

and

$$u_{\varepsilon}(x) = u_0(x) + 4\varepsilon\gamma(x_0)(u_{in}(x_0) - u_0(x_0))N(x, x_0) + o(\varepsilon) \text{ if } d=3.$$
In this case, the sensitivity of replacing the homogeneous Neumann boundary condition on $\omega_{x_0,\varepsilon} \subset \Gamma$ by the inhomogeneous Dirichlet conditions $u = u_{\text{in}}$ equals:

$$J((\Gamma_D)_{x_0,\varepsilon}) = J(\Gamma_D) + \frac{\pi}{|\log \varepsilon|} \gamma(x_0) (u_{\text{in}}(x_0) - u_0(x_0)) p_0(x_0) + o\left(\frac{1}{|\log \varepsilon|}\right) \text{ if } d=2,$$

and

$$J((\Gamma_D)_{x_0,\varepsilon}) = J(\Gamma_D) + 4\varepsilon \gamma(x_0) (u_{\text{in}}(x_0) - u_0(x_0)) p_0(x_0) + o(\varepsilon) \text{ if } d=3,$$

where the adjoint state $p_0$ is again the solution to (4.28).

**Remark 4.7.** A counterpart of the arguments and results of this section holds in the case where homogeneous Dirichlet boundary conditions are replaced by homogeneous Neumann boundary conditions on a small surface disk, i.e. in the case where $G = \Gamma$ and $x_0 \in \Gamma_D$, see [33].

## 4.4. Extension: replacement of the homogeneous Neumann boundary condition by an inhomogeneous Neumann condition on a small subset

This section considers the situation where the homogeneous Neumann boundary condition is replaced by an inhomogeneous Neumann condition on a small surface disk $\omega_{x_0,\varepsilon} \subset \Gamma$. The voltage potential $u_\varepsilon$ in this perturbed situation is the solution to the boundary value problem:

$$\left\{
\begin{array}{ll}
-\text{div}(\gamma \nabla u_{\varepsilon}) = f & \text{in } \Omega, \\
u_{\varepsilon} = 0 & \text{on } \Gamma_D, \\
\gamma \frac{\partial u_{\varepsilon}}{\partial n} = 0 & \text{on } \Gamma \setminus \overline{\omega_{x_0,\varepsilon}}, \\
\gamma \frac{\partial u_{\varepsilon}}{\partial n} = g & \text{on } \Gamma_N \cup \omega_{x_0,\varepsilon},
\end{array}
\right.
\qquad (4.29)$$

where we recall that $g: \mathbb{R}^d \to \mathbb{R}$ is a given smooth function.

**Theorem 4.8.** The following asymptotic expansions hold, at any point $x \in \bar{\Omega} \setminus (\Sigma_D \cup \{x_0\})$:

$$u_{\varepsilon}(x) = u_0(x) + 2\varepsilon g(x_0) N(x, x_0) + o(\varepsilon) \text{ if } d=2,$$

and

$$u_{\varepsilon}(x) = u_0(x) + \varepsilon^2 \pi g(x_0) N(x, x_0) + o(\varepsilon^2) \text{ if } d=3.$$

*Sketch of the proof.* The derivation of these formulas essentially follows the trail of the proof of Theorem 4.4, in a much simpler version. Again, for simplicity, we assume that $x_0 = 0$, $n(0) = e_d$, and that $\partial\Omega$ is flat around 0, see (4.15). Let $r_\varepsilon := u_\varepsilon - u_0$ be the error between the perturbed and background potentials, defined by (4.29) and (3.2), respectively. This function satisfies the following boundary value problem:

$$\left\{
\begin{array}{ll}
-\text{div}(\gamma \nabla r_{\varepsilon}) = 0 & \text{in } \Omega, \\
r_{\varepsilon} = 0 & \text{on } \Gamma_D, \\
\gamma \frac{\partial r_{\varepsilon}}{\partial n} = 0 & \text{on } (\Gamma_N \cup \Gamma) \setminus \bar{\omega}_{\varepsilon}, \\
\gamma \frac{\partial r_{\varepsilon}}{\partial n} = g & \text{on } \omega_{\varepsilon}.
\end{array}
\right.
\qquad (4.30)$$
From the definition of the Green's function $N(x, y)$ in (4.5), it holds:

$$r_{\varepsilon}(x) = - \int_{\Omega} \mathrm{div}_{y}(\gamma(y)\nabla_{y}N(x, y))r_{\varepsilon}(y) \, dy,$$

and so, after integration by parts:

$$r_{\varepsilon}(x) = - \int_{\partial\Omega} \gamma(y) \frac{\partial N}{\partial n_y}(x, y) r_{\varepsilon}(y) \, ds(y) + \int_{\partial\Omega} \gamma(y) \frac{\partial r_{\varepsilon}}{\partial n}(y) N(x, y) \, ds(y).$$

Now using the boundary conditions satisfied by the functions $r_{\varepsilon}$ and $y \mapsto N(x, y)$, it follows:

$$\begin{aligned} r_{\varepsilon}(x) &= \int_{\omega_{\varepsilon}} g(y) N(x, y) \, ds(y) \\ &= \varepsilon^{d-1} \int_{\mathbb{D}_1} g(\varepsilon z) N(x, \varepsilon z) \, ds(z) \end{aligned}$$

and so

$$r_{\varepsilon}(x) = \varepsilon^{d-1} g(0) N(x, 0) \left( \int_{\mathbb{D}_1} ds(z) \right) + o(\varepsilon^{d-1}),$$

which yields the desired result.

We now turn to the computation of the sensitivity of the objective function $J(\Gamma_N)$ in (3.3) in the present context, i.e. we consider the asymptotic expansion of the quantity of interest:

$$J((\Gamma_N)_{x_0, \varepsilon}) := \int_{\Omega} j(u_{\varepsilon}) \, dx.$$

The proof of the following result is omitted, as it is completely similar to that of Corollary 4.5.

**Corollary 4.9.** Let $x_0$ be a given point in $\Gamma$; the perturbed value $J((\Gamma_N)_{x_0, \varepsilon})$ of the objective function $J(\Gamma_N)$, accounting for the replacement of the homogeneous Neumann boundary condition by the inhomogeneous Neumann boundary condition $\gamma \frac{\partial u}{\partial n} = g$ on a “small” surface disk $\omega_{x_0, \varepsilon}$, has the following expansion:

$$J((\Gamma_N)_{x_0, \varepsilon}) = J(\Gamma_N) - 2\varepsilon g(x_0) p_0(x_0) + o(\varepsilon) \quad \text{if } d = 2,$$

and

$$J((\Gamma_N)_{x_0, \varepsilon}) = J(\Gamma_N) - \pi \varepsilon^2 g(x_0) p_0(x_0) + o(\varepsilon^2) \quad \text{if } d = 3,$$

where $p_0 \in H^1(\Omega)$ is the solution to (4.28).

## 5. OPTIMIZATION OF REGIONS BEARING BOUNDARY CONDITIONS OF THE HELMHOLTZ EQUATION

In this section, we slip into the physical context of acoustics, where the state function $u$ is the solution to the Helmholtz equation, and we adapt the previous material to this setting. As we have mentioned in the introduction, exact or approximate shape derivatives of functionals $J(G)$ depending on a region $G \subset \partial\Omega$ bearing boundary conditions of this physical problem can be calculated in the same spirit as in the context of electrostatics, tackled in Section 3; we refer to Section 7.4 for an application example. We therefore focus on the
calculation of the topological derivative of such a function $J(G)$, mostly presenting the differences with the analysis conducted in Section 4. To emphasize the parallel between both situations, we retain the same notations as in there, insofar as possible.

## 5.1. Presentation of the model

For simplicity, we focus on a model interior Helmholtz problem, excerpted from [26, 75, 76]. The arguments exposed in here can be readily adapted to the situation of an infinite propagation medium, see Section 7.4 for an example in this context. Let $\Omega$ be a smooth bounded domain in $\mathbb{R}^d$, whose boundary $\partial\Omega$ is composed of two disjoint parts, namely:

$$ \partial\Omega = \overline{\Gamma_N} \cup \overline{\Gamma_R}, $$

where:

* The region $\Gamma_N$ is the support of homogeneous Neumann boundary conditions,

* The region $\Gamma_R$ is the support of a homogeneous impedance (i.e. Robin) boundary condition.

In this situation, the state function describes the pressure of an acoustic wave. In the “background” situation, it is the solution $u_0$ to the following equation:

$$ \left\{ \begin{array}{ll} -\operatorname{div}(\gamma \nabla u_0) - k^2 u_0 = f & \text{in } \Omega, \\ \gamma \frac{\partial u_0}{\partial n} = 0 & \text{on } \Gamma_N, \\ \gamma \frac{\partial u_0}{\partial n} + iku_0 = 0 & \text{on } \Gamma_R. \end{array} \right. \qquad (5.1) $$

Here, we have denoted the wave number by $k > 0$; the coefficient $\gamma$ encodes the physical properties of the medium (it is the inverse of its admittance), and it is uniformly bounded away from 0 and $\infty$, see (3.1). The smooth right-hand side $f: \mathbb{R}^d \to \mathbb{R}$ represents a source acting in the medium. From the physical viewpoint, the homogeneous Neumann boundary condition on $\Gamma_N$ accounts for a “hard wall”, where perfect reflection of the wave occurs, while the Robin condition encodes a partial absorption of the energy.

**Remark 5.1.** The well-posedness of (5.1) is not a consequence of the usual Lax-Milgram theory, as the involved operator is not coercive. This property is a consequence of the fact that the coefficient of the impedance boundary condition has a non zero imaginary part; mathematically, it follows from the combination of a unique continuation principle with the so-called Banach–Necas–Babuška theorem, see for instance Chap. 25 and 35 in [77].

Let us recall that the fundamental solution $\Gamma(x, y)$ of the Helmholtz operator $u \mapsto -\Delta u - k^2 u$ in the free space $\mathbb{R}^d$ ($d = 2, 3$) is given by the following formulas:

$$ \Gamma(x, y) = \begin{cases} \frac{-1}{4i} H_0^{(1)}(k|x-y|) & \text{if } d = 2, \\ \frac{e^{ik|x-y|}}{4\pi|x-y|} & \text{if } d = 3, \end{cases} \qquad (5.2) $$

where $H_0^{(1)}$ is the Hankel function of the first kind and of order 0, which is the solution to the ordinary differential equation:

$$ \frac{1}{r} \frac{d}{dr} \left( r \frac{dH}{dr}(r) \right) + k^2 H(r) = 0, \quad \text{for } r > 0, $$

see for instance [78].
The Green's function $N(x, y)$ for the background problem (5.1) can be constructed from this datum by standard means, as the solution to:

$$ egin{cases} -\operatorname{div}_y(\gamma \nabla_y N(x, y)) - k^2 N(x, y) = \delta_{y=x} & \text{in } \Omega, \\ \gamma \frac{\partial N}{\partial n_y}(x, y) = 0 & \text{for } y \in \Gamma_N, \\ \gamma \frac{\partial N}{\partial n_y}(x, y) + ikN(x, y) = 0 & \text{for } y \in \Gamma_R. \end{cases} \tag{5.3} $$

**Remark 5.2.** In this context also, a Green's function $L_\gamma(x, y)$ for the version of the boundary value problem (5.1) posed on the lower half-space $H$, featuring a constant coefficient $\gamma > 0$ and equipped with homogeneous Neumann boundary conditions on $\partial H$ can be constructed by the method of images:

$$ L_{\gamma}(x, y) = \frac{1}{\gamma}(\Gamma(x, y) + \Gamma(x, -y)). $$

## 5.2. Sensitivity of the solution to the Helmholtz equation with respect to a singular perturbation of the Neumann boundary condition

To set ideas, we consider a perturbed version of the problem (5.1) where the homogeneous Neumann boundary condition is replaced by an impedance boundary condition on a “small” surface disk $\omega_{x_0, \varepsilon} \subset \Gamma_N$ around the point $x_0 \in \Gamma_N$, that is:

$$ egin{cases} -\operatorname{div}(\gamma \nabla u_{\varepsilon}) - k^2 u_{\varepsilon} = f & \text{in } \Omega, \\ \gamma \frac{\partial u_{\varepsilon}}{\partial n} = 0 & \text{on } \Gamma_N \setminus \overline{\omega_{x_0, \varepsilon}}, \\ \gamma \frac{\partial u_{\varepsilon}}{\partial n} + iku_{\varepsilon} = 0 & \text{on } \Gamma_R \cup \omega_{x_0, \varepsilon}, \end{cases} \tag{5.4} $$

Our main result concerning the asymptotic behavior of $u_\varepsilon$ when $\varepsilon \to 0$ is the following:

**Theorem 5.3.** Let $x_0 \in \Gamma_N$ be given. The following asymptotic expansions hold true, at any point $x \in \overline{\Omega} \setminus ((\overline{\Gamma_N} \cap \overline{\Gamma_R}) \cup \{x_0\})$:

$$ u_{\varepsilon}(x) = u_0(x) - 2\varepsilon iku_0(x_0)N(x, x_0) + o(\varepsilon) \text{ if } d=2, $$

and

$$ u_{\varepsilon}(x) = u_0(x) - \pi\varepsilon^2 iku_0(x_0)N(x, x_0) + o(\varepsilon^2) \text{ if } d=3. $$

*Sketch of the proof.* We rely on the formal argument employed in our treatment of Theorem 4.4, and we only sketch the main frame for brevity. Without loss of generality, we set $x_0 = 0$ and $\omega_\varepsilon := \omega_{x_0, \varepsilon}$; we furthermore rely on the simplifying assumption that $\gamma$ is constant, and $\partial\Omega$ is completely flat in a neighborhood of 0, see (4.15). The error $r_\varepsilon = u_\varepsilon - u_0$ is the unique solution in the space $H^1(\Omega; \mathbb{C})$ of complex-valued $H^1$ functions to the following boundary value problem:

$$ egin{cases} -\operatorname{div}(\gamma \nabla r_{\varepsilon}) - k^2 r_{\varepsilon} = 0 & \text{in } \Omega, \\ \gamma \frac{\partial r_{\varepsilon}}{\partial n} = 0 & \text{on } \Gamma_N \setminus \overline{\omega_{\varepsilon}}, \\ \gamma \frac{\partial r_{\varepsilon}}{\partial n} + ikr_{\varepsilon} = 0 & \text{on } \Gamma_R, \\ \gamma \frac{\partial r_{\varepsilon}}{\partial n} + ikr_{\varepsilon} = -iku_0 & \text{on } \omega_{\varepsilon}. \end{cases} \tag{5.5} $$
Like in the proof of Theorem 4.4, a priori estimates for (5.5) allow to prove that

$$r_{\varepsilon} \xrightarrow{\varepsilon \to 0} 0 \text{ in } H^1(\Omega; \mathbb{C}). \qquad (5.6)$$

We now proceed in four steps.

**Step 1:** We construct a representation formula for the values of $r_{\varepsilon}$ “far” from 0 in terms of its values inside the region $\omega_{\varepsilon}$.

To this end, we rely on the Green’s function $N(x, y)$ of the background problem (5.1); it holds:

$$ \begin{aligned} r_{\varepsilon}(x) &= -\int_{\Omega} (\operatorname{div}_y(\gamma \nabla_y N(x, y)) + k^2 N(x, y)) r_{\varepsilon}(y) \, dy \\ &= -\int_{\partial \Omega} \gamma \frac{\partial N}{\partial n_y}(x, y) r_{\varepsilon}(y) \, ds(y) + \int_{\Omega} (\gamma \nabla_y N(x, y) \cdot \nabla r_{\varepsilon}(y) - k^2 N(x, y) r_{\varepsilon}(y)) \, dy \\ &= -\int_{\partial \Omega} \gamma \frac{\partial N}{\partial n_y}(x, y) r_{\varepsilon}(y) \, ds(y) + \int_{\partial \Omega} \gamma \frac{\partial r_{\varepsilon}}{\partial n_y}(y) N(x, y) \, ds(y), \end{aligned} $$

where the second and third lines follow from integration by parts. Recalling the boundary conditions satisfied by $r_{\varepsilon}$ and $y \mapsto N(x, y)$ in (5.1) and (5.3), the above expression simplifies to:

$$r_{\varepsilon}(x) = \int_{\omega_{\varepsilon}} \gamma \frac{\partial r_{\varepsilon}}{\partial n_y}(y) N(x, y) \, ds(y),$$

and so, after rescaling:

$$r_{\varepsilon}(x) = \int_{\mathbb{D}_1} \varphi_{\varepsilon}(z) N(x, \varepsilon z) \, ds(z), \quad \text{where } \varphi_{\varepsilon}(z) := \varepsilon^{d-1} \left( \gamma \frac{\partial r_{\varepsilon}}{\partial n} \right) (\varepsilon z) \in \tilde{H}^{-1/2}(\mathbb{D}_1). \qquad (5.7)$$

**Step 2:** We characterize the function $\varphi_{\varepsilon}$ by an integral equation.

To achieve this, we essentially repeat the calculations from the first step, except that we use the explicit Green’s function for the half-space $L_{\gamma}(x, y)$ discussed in Remark 5.2 in place of the more abstract Green’s function $N(x, y)$ for the background equation (5.1). We first obtain:

$$r_{\varepsilon}(x) = -\int_{\partial \Omega} \gamma \frac{\partial L_{\gamma}}{\partial n_y}(x, y) r_{\varepsilon}(y) \, ds(y) + \int_{\partial \Omega} \gamma \frac{\partial r_{\varepsilon}}{\partial n_y}(y) L_{\gamma}(x, y) \, ds(y).$$

Invoking the boundary conditions satisfied by the functions $r_{\varepsilon}$ and $y \mapsto L_{\gamma}(x, y)$, it follows that:

$$r_{\varepsilon}(x) = \int_{\omega_{\varepsilon}} \gamma \frac{\partial r_{\varepsilon}}{\partial n_y} L_{\gamma}(x, y) \, ds(y) + K r_{\varepsilon}(x),$$

where the term $Kr_{\varepsilon}(x)$ gathers integrals of $r_{\varepsilon}$ and its derivatives whose supports are “far” from 0.

We now let $x$ tend to $\omega_{\varepsilon}$ in this formula, and insert the resulting expression in the impedance boundary condition in (5.5); this yields:

$$\forall x \in \omega_{\varepsilon}, \quad \gamma \frac{\partial r_{\varepsilon}}{\partial n}(x) + ik \int_{\omega_{\varepsilon}} \gamma \frac{\partial r_{\varepsilon}}{\partial n_y} L_{\gamma}(x, y) \, ds(y) + K r_{\varepsilon}(x) = -iku_0(x).$$
Finally, the convergence (5.6) and a change of variables in the above integral yield the desired integral equation for the function $\varphi_{\varepsilon}$ in (5.7):

$$ \forall x \in \mathbb{D}_1, \quad \frac{1}{\varepsilon^{d-1}} \varphi_{\varepsilon}(x) + ik \int_{\mathbb{D}_1} \varphi_{\varepsilon}(z)L(\varepsilon x, \varepsilon z) \, ds(z) = -iku_0(0) + o(1). \quad (5.8) $$

**Step 3: We use the integral equation (5.8) to glean information about the asymptotic behavior of $\varphi_{\varepsilon}$**.

To this end, we observe that, because of the homogeneity of $L(\cdot, \cdot)$, the integral term on the left-hand side of (5.8) is of order $|\log \varepsilon|$ if $d=2$, and of order $\varepsilon^{-(d-2)}$ if $d \ge 3$; it is thus negligible with respect to the first term in the left-hand side of this equation. Hence, by taking the mean value in (5.8), we immediately obtain the following relations:

$$ \int_{\mathbb{D}_1} \varphi_{\varepsilon}(z) \, ds(z) = \begin{cases} -2ik\varepsilon u_0(0) + o(\varepsilon) & \text{if } d=2, \\ -\pi\varepsilon^2 ik u_0(0)\pi + o(\varepsilon^2) & \text{if } d=3, \end{cases} $$

which encompass the needed information for our purpose.

**Step 4: We pass to the limit in the representation formula (5.7).**

The application of the Lebesgue dominated convergence theorem to the representation formula (5.7) yields:

$$ r_{\varepsilon}(x) = \left( \int_{\mathbb{D}_1} \varphi_{\varepsilon}(z) \, ds(z) \right) \left( N(x, 0) + o(1) \right), $$

and the desired result follows immediately from the formulas derived in Step 3. $\square$

### 5.3. Calculation of the topological derivative of a functional depending on the region $\Gamma_R$

In this section, we cast the previous analysis in the context of shape and topology optimization, along the lines of Section 2.2. We consider a function $J(\Gamma_R)$ of the region $\Gamma_R \subset \partial\Omega$ bearing the impedance boundary condition of the Helmholtz equation (5.1), of the form:

$$ J(\Gamma_R) = \int_{\Omega} j(u_{\Gamma_R}) \, dx, $$

where we have denoted by $u_{\Gamma_R}$ the solution to (5.1). Here, $j: \mathbb{C} \to \mathbb{R}$ is smooth when it is seen as a function defined on the real vector space $\mathbb{R}^2$, and it satisfies the growth conditions (3.4); with a small abuse of notation, we let

$$ \forall u = (u_1, u_2) \in \mathbb{C} \approx \mathbb{R}^2, \quad j'(u) := \frac{\partial j}{\partial u_1}(u) + i \frac{\partial j}{\partial u_2}(u) \in \mathbb{C}. $$

The sensitivity of the function $J(\Gamma_R)$ with respect to the addition of a small surface disk $\omega_{x_0, \varepsilon}$ to $\Gamma_R$ is then given by the next result.

**Corollary 5.4.** Let $x_0$ be a given point on $\Gamma_N$. The mapping $\varepsilon \mapsto J((\Gamma_R)_{x_0, \varepsilon})$ has the following asymptotic expansion:

$$ J((\Gamma_R)_{x_0, \varepsilon}) = J(\Gamma_R) + 2\varepsilon k \operatorname{Im} \left( \overline{u_0(x_0)} p_0(x_0) \right) + o(\varepsilon) \quad \text{if } d=2, $$

and

$$ J((\Gamma_R)_{x_0, \varepsilon}) = J(\Gamma_R) + \pi\varepsilon^2 k \operatorname{Im} \left( \overline{u_0(x_0)} p_0(x_0) \right) + o(\varepsilon^2) \quad \text{if } d=3, $$
where Im($u$) is the imaginary part $u_2$ of a complex number $u = (u_1, u_2) \in \mathbb{C}$ and the adjoint state $p_0 \in H^1(\Omega; \mathbb{C})$ is the unique solution to the following boundary value problem:

$$ egin{cases} -\operatorname{div}(\gamma \nabla p_0) - k^2 p_0 = -j'(u_0) & \text{in } \Omega, \\ \gamma \frac{\partial p_0}{\partial n} = 0 & \text{on } \Gamma_N, \\ \gamma \frac{\partial p_0}{\partial n} - ikp_0 = 0 & \text{on } \Gamma_R. \end{cases} \qquad (5.9) $$

*Hint of the proof.* The proof is similar to that of Corollary 4.5. Considering the case $d = 3$ to set ideas, the Lebesgue dominated convergence theorem shows that:

$$ J((\Gamma_R)_{x_0, \varepsilon}) = J(\Gamma_R) - \pi \varepsilon^2 k \operatorname{Re} \left( \int_{\Omega} j'(u_0(x)) \overline{(iu_0(x)N(x, x_0))} dx \right) + o(\varepsilon^2). $$

To obtain this formula, we have used the basic calculus rule:

$$ \forall u, h \in \mathbb{C}, \quad j(u+h) = j(u) + \operatorname{Re}(j'(u)\bar{h}) + o(h), \quad \text{where} \quad \frac{o(h)}{|h|} \xrightarrow{h \to 0} 0, $$

and $\operatorname{Re}(u) = u_1$ is the real part of a complex number $u = (u_1, u_2) \in \mathbb{C}$. Since the complex conjugate $\overline{N(x, y)}$ is the Green's function for the boundary value problem (5.9) (where the sign in the Robin boundary condition is changed compared to (5.1)), the representation formula (4.9) for $p_0$ reads:

$$ p_0(x_0) = - \int_{\Omega} j'(u_0(x)) \overline{N(x, x_0)} dx, $$

and so:

$$ \begin{aligned} J((\Gamma_R)_{x_0, \varepsilon}) &= J(\Gamma_R) - \pi \varepsilon^2 k \operatorname{Re} \left( \overline{iu_0(x_0)} \int_{\Omega} j'(u_0(x)) \overline{N(x, x_0)} dx \right) + o(\varepsilon^2) \\ &= J(\Gamma_R) + \pi \varepsilon^2 k \operatorname{Im} \left( \overline{u_0(x_0)} p_0(x_0) \right) + o(\varepsilon^2), \end{aligned} $$

which is the desired formula.

<figure><img src="image_8.png" /></figure>

# 6. OPTIMIZATION OF THE SUPPORT OF BOUNDARY CONDITIONS FOR THE LINEAR ELASTICITY SYSTEM

In this section, we adapt the findings of the previous Sections 3 to 5 to the realm of structural mechanics, governed by the system of linear elasticity. The calculation of the (exact or approximate) shape derivative of a functional depending on a region supporting the boundary conditions of this system follows exactly the trail conducted in Section 3, up to an increased level of tediousness. In contrast, the asymptotic analysis of the effect of singular changes in such regions – and thereby, the calculation of topological derivatives – involves non trivial specificities. For this reason, in this section, we focus on the calculation of topological derivatives of functionals of boundary regions featured in the linear elasticity equations, and we refer to the applications of Sections 7.5 and 7.6 for typical shape derivative formulas in this context and illustrations of their practical use. After presenting the mathematical model at stake in Section 6.1 and a few technical preliminaries in Section 6.2, we discuss in Sections 6.3 and 6.4 the calculation of the sensitivity of the elastic displacement and that of a related quantity of interest with respect to the addition of a vanishingly small Dirichlet region. Multiple variations of this model are available, which can be treated by simple adaptations of this material, see Remark 6.5. In order to emphasize the parallel between the present discussions and those in Sections 4 and 5, we retain our previous notations whenever possible.
## 6.1. Presentation of the linear elasticity setting

In this section, $\Omega$ stands for a mechanical structure in $\mathbb{R}^d$ ($d=2,3$), whose boundary $\partial\Omega$ is composed of three disjoint pieces:

$$ \partial\Omega = \overline{\Gamma_D} \cup \overline{\Gamma_N} \cup \overline{\Gamma}, $$

where:

* The region $\Gamma_D$ is held fixed;

* The region $\Gamma_N$ is subjected to surface loads $g: \mathbb{R}^d \to \mathbb{R}^d$;

* The remaining part $\Gamma$ is traction-free.

Assuming smooth body forces $f: \mathbb{R}^d \to \mathbb{R}^d$, the displacement $u$ of $\Omega$ is the unique solution in the space $H^1(\Omega)^d$ to the system of linear elasticity:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(Ae(u)) = f & \text{in } \Omega, \\ Ae(u)n = g & \text{on } \Gamma_N, \\ Ae(u)n = 0 & \text{on } \Gamma, \\ u = 0 & \text{on } \Gamma_D. \end{array} \right. \qquad (6.1) $$

Here, $e(u) := \frac{1}{2}(\nabla u + \nabla u^T)$ is the strain tensor associated to a displacement field $u: \Omega \to \mathbb{R}^d$, and $A$ is the Hooke's tensor, defined by:

$$ \text{For all symmetric } d \times d \text{ matrix } e, \quad Ae = 2\mu e + \lambda \mathrm{tr}(e)I, $$

where $\lambda$ and $\mu$ are the Lamé coefficients of the elastic material. Note that these coefficients are often better expressed in terms of the more physical Young's modulus $E$ and Poisson's ratio $\nu$:

$$ \mu = \frac{E}{2(1 + \nu)}, \quad \text{and} \quad \lambda = \begin{cases} \frac{E\nu}{(1+\nu)(1-2\nu)} & \text{if } d=3, \\ \frac{E\nu}{1-\nu^2} & \text{if } d=2 \text{ (plane stress).} \end{cases} $$

We refer to classical treaties such as [79] for more exhaustive introductions to the physical context of linear elasticity.

## 6.2. The fundamental solution to the linear elasticity system

The fundamental solution to the linear elasticity operator $u \mapsto -\mathrm{div}(Ae(u))$ in the free space $\mathbb{R}^d$ is the Kelvin matrix, defined by

$$ \forall i, j = 1, \dots, d, \quad \Gamma_{ij}(x, y) = \begin{cases} \frac{\alpha}{4\pi} \frac{\delta_{ij}}{|x-y|} + \frac{\beta}{4\pi} \frac{(x_i - y_i)(x_j - y_j)}{|x-y|^3} & \text{if } d=3, \\ -\frac{\alpha}{2\pi} \delta_{ij} \log|x-y| + \frac{\beta}{2\pi} \frac{(x_i - y_i)(x_j - y_j)}{|x-y|^2} & \text{if } d=2, \end{cases} \qquad (6.2) $$

where the constants $\alpha$ and $\beta$ are given by:

$$ \alpha = \frac{1}{2} \left( \frac{1}{\mu} + \frac{1}{2\mu + \lambda} \right), \quad \text{and} \quad \beta = \frac{1}{2} \left( \frac{1}{\mu} - \frac{1}{2\mu + \lambda} \right), $$
see e.g. [80, 81]. This means that, for each point $x \in \mathbb{R}^d$ and each index $j = 1, \dots, d$, the $j$th column vector $\Gamma_j(x, \cdot)$ of $\Gamma(x, \cdot)$ is the solution to the equation:

$$-\mathrm{div}_y(Ae_y(\Gamma_j(x, y))) = \delta_{y=x}e_j \quad \text{in the sense of distributions in } \mathbb{R}^d,$$

where $e_j$ is the $j$th vector of the canonical basis of $\mathbb{R}^d$. From the physical point of view, if $a \in \mathbb{R}^d$ is a vector, $y \mapsto \Gamma(x, y)a$ is the displacement resulting from a point load $a$ applied at $x$.

The variational counterpart of (6.2) reads as follows: for a sufficiently smooth function $\varphi \in C_c^\infty(\mathbb{R}^d)^d$, it holds:

$$\forall x \in \mathbb{R}^d, \quad \varphi_j(x) = \int_{\mathbb{R}^d} Ae_y(\Gamma_j(x, y)) : e(\varphi)(y) \, dy.$$

Using the fundamental solution (6.2), it is possible to construct the Green's function $N(x, y)$ for the problem (6.1): for $j = 1, \dots, d$, the $j$th column $N_j$ of $N$ satisfies the equation

$$\left\{
\begin{array}{ll}
-\mathrm{div}_y(Ae_y(N_j(x, y))) = \delta_{y=x}e_j & \text{in } \Omega, \\
N_j(x, y) = 0 & \text{for } y \in \Gamma_D, \\
Ae_y(N_j(x, y)) = 0 & \text{for } y \in \Gamma \cup \Gamma_N.
\end{array}
\right.
\qquad (6.3)$$

This system has the following variational characterization:

$$\text{For all } \varphi \in C_c^\infty(\mathbb{R}^d)^d \text{ s.t. } \varphi = 0 \text{ on } \Gamma_D, \quad \varphi_j(x) = \int_{\Omega} Ae_y(N_j(x, y)) : e(\varphi)(y) \, dy. \qquad (6.4)$$

**Remark 6.1.** As in Sections 4 and 5, we shall need another Green's function $L(x, y)$ for the linear elasticity system (6.1), tuned to the lower half-space $H$ (see (2.6)), satisfying homogeneous Neumann conditions on $\partial H$. Unfortunately, due to the vector-valued nature of the elasticity system, the expression of $L(x, y)$ cannot be obtained from that of $\Gamma(x, y)$ by the method of images. This function – called the Mindlin function in the literature – has the following expression in 3d, for two points $x \neq y$ on $\partial H$:

$$L_{ij}(x, y) = \frac{1 - \nu}{2\pi\mu|x - y|} \delta_{ij} + \frac{\nu}{2\pi\mu} \frac{(x_i - y_i)(x_j - y_j)}{|x - y|^3}, \quad i, j = 1, 2, \qquad (6.5)$$

$$L_{3j}(x, y) = -\frac{1 - 2\nu}{4\pi\mu} \frac{x_j - y_j}{|x - y|^2}, \quad j = 1, 2, \text{ and } L_{33}(x, y) = \frac{1 - \nu}{2\pi\mu|x - y|}. \qquad (6.6)$$

We refer to [82] for the original derivation of this formula, see also [83, 84]. Note that, in these references, the formulas look different, as the Green's function is provided for the upper half-space.

In 2d, setting $\bar{\nu} = \frac{\nu}{1+\nu}$, the components of this Green's function read, according to [83] Section 2.8.2:

$$L_{11}(x, y) = -\frac{1 - \bar{\nu}}{\pi\mu} \log|x - y| + \frac{3 - 4\bar{\nu}}{8\pi\mu(1 - \bar{\nu})}, \qquad (6.7)$$

$$L_{12}(x, y) = -\frac{1}{2\pi\mu}(1 - 2\bar{\nu})\Theta, \text{ where } \Theta = \begin{cases} \frac{\pi}{2} & \text{if } x_1 > y_1, \\ -\frac{\pi}{2} & \text{otherwise,} \end{cases} \text{ and } L_{22}(x, y) = -\frac{1 - \bar{\nu}}{\pi\mu}\log|x - y|. \qquad (6.8)$$
### 6.3. Asymptotic expansion of the elastic displacement under singular perturbations of the Dirichlet region

In this section, we consider the situation where the reference “background” problem (6.1), whose solution is as usual denoted by $u_0$, is perturbed by the replacement of the homogeneous Neumann boundary conditions by homogeneous Dirichlet boundary conditions on a small surface disk $\omega_{x_0,\varepsilon} \subset \Gamma_N$ centered at a given point $x_0 \in \Gamma$. The displacement $u_\varepsilon$ of $\Omega$ in this perturbed situation is the unique solution to the following boundary value problem:

$$ egin{cases} -\operatorname{div}(Ae(u_\varepsilon)) = f & \text{in } \Omega, \\ u_\varepsilon = 0 & \text{on } \Gamma_D \cup \omega_{x_0,\varepsilon}, \\ Ae(u_\varepsilon)n = g & \text{on } \Gamma_N, \\ Ae(u_\varepsilon)n = 0 & \text{on } \Gamma \setminus \overline{\omega_{x_0,\varepsilon}}. \end{cases} \qquad (6.9) $$

The main result is the following.

**Theorem 6.2.** Let $x_0$ be a given point in $\Gamma$. Then the following asymptotic expansions hold true about the perturbed displacement $u_\varepsilon$ in (6.9):

$$ u_{\varepsilon,j}(x) = u_{0,j}(x) - \frac{1}{|\log \varepsilon|} \frac{\pi \mu}{1-\bar{\nu}} u_0(x_0) \cdot N_j(x, x_0) + o\left(\frac{1}{|\log \varepsilon|}\right), \quad j=1,2, \text{ if } d=2, $$

and

$$ u_{\varepsilon,j}(x) = u_{0,j}(x) - \varepsilon M u_0(x_0) \cdot N_j(x, x_0) + o(\varepsilon), \quad j=1,2,3, \text{ if } d=3. $$

In the above formula, $M \in \mathbb{R}^{3 \times 3}$ is a polarization tensor, whose entries read:

$$ M_{ij} = \int_{\mathbb{D}_1} T_L^{-1} e_j \cdot e_i \, ds, \quad i, j = 1, 2, 3. \qquad (6.10) $$

This definition involves the integral operator

$$ T_L : \tilde{H}^{-1/2}(\mathbb{D}_1)^d \to H^{1/2}(\mathbb{D}_1)^d, \quad T_L\varphi(x) = \int_{\mathbb{D}_1} L(x, z)\varphi(z) \, ds(z), \quad x \in \mathbb{D}_1, \qquad (6.11) $$

whose kernel is the Mindlin function $L(x, y)$ given in (6.5) and (6.6).

### Remark 6.3.

* The above definition of the polarization tensor $M$ in (6.10) hinges on the invertibility of the operator $T_L$ in (6.11). To the best of our knowledge, this fact is not known in the literature, and it is not a straightforward adaptation of the counterpart result in the setting of the conductivity equation (see Prop. 4.3). The derivation presented below proceeds formally under the assumption that this fact holds true.

* The operator $T_L$ has an enlightening physical interpretation: $T_L\varphi : \mathbb{D}_1 \to \mathbb{R}^d$ is the displacement induced by a force $\varphi : \mathbb{D}_1 \to \mathbb{R}^d$ applied to the unit disk $\mathbb{D}_1$. Hence, for $j = 1, 2, 3$, the $j^{th}$ column of the polarization tensor $M$ is the total magnitude of the force that should be applied on $\mathbb{D}_1$ to realize a uniform, unit displacement in the direction $e_j$. Interestingly, “classical” calculations in contact mechanics allow to identify some of the entries of this tensor, see for instance the so-called “flat punch” or “indentation” problems in [85–87]. However, to the best of our knowledge, the complete analytical calculation of $M$ is not available in the literature, and this task has to be realized numerically: see Section 7.2 about the numerical method employed to solve the integral equation associated to $T_L$.
Sketch of the proof. Again, we rely on a formal calculation along the trail of Theorem 4.4, under the assumption (4.15) that $x_0 = 0$, $n(0) = e_d$, and that the boundary $\partial\Omega$ is flat in a neighborhood of 0. Let us introduce the error $r_\varepsilon = u_\varepsilon - u_0 \in H^1(\Omega)^d$, which is the unique solution to the following boundary value problem:

$$ egin{cases} -\operatorname{div}(Ae(r_\varepsilon)) = 0 & \text{in } \Omega, \\ r_\varepsilon = 0 & \text{on } \Gamma_D, \\ r_\varepsilon = -u_0 & \text{on } \omega_\varepsilon, \\ Ae(r_\varepsilon)n = 0 & \text{on } (\Gamma_N \cup \Gamma) \setminus \bar{\omega}_\varepsilon. \end{cases} \tag{6.12} $$

Step 1: We construct a representation formula for the value $r_\varepsilon(x)$ of the error at an arbitrary point $x \in \bar{\Omega} \setminus \{0\}$ in terms of its values inside $\omega_\varepsilon$.

To this end, we use the definition (6.3) of the Green's function $N(x, y)$ for the background elasticity problem (6.1), which yields, for each component $j = 1, \dots, d$:

$$ egin{aligned} r_{\varepsilon,j}(x) &= -\int_{\Omega} \operatorname{div}_y (Ae_y(N_j(x, y))) \cdot r_\varepsilon(y) \, dy \\ &= -\int_{\partial\Omega} Ae_y(N_j(x, y))n(y) \cdot r_\varepsilon(y) \, ds(y) + \int_{\Omega} Ae_y(N_j(x, y)) \cdot e(r_\varepsilon)(y) \, dy. \end{aligned} $$

Using the boundary conditions in (6.3) and (6.12) for the functions $r_\varepsilon$ and $y \mapsto N(x, y)$, as well as another integration by parts, we obtain:

$$ r_{\varepsilon,j}(x) = - \int_{\Omega} \operatorname{div}(Ae(r_{\varepsilon}))(y) \cdot N_j(x, y) \, dy + \int_{\partial\Omega} Ae(r_{\varepsilon})(y)n(y) \cdot N_j(x, y) \, ds(y). $$

Invoking again the problem satisfied by $r_\varepsilon$ in (6.12), we end up with

$$ r_{\varepsilon,j}(x) = \int_{\omega_\varepsilon} Ae(r_\varepsilon)(y)n(y) \cdot N_j(x, y) \, ds(y). $$

By a change of variables in the above integral, we obtain the desired representation formula:

$$ r_{\varepsilon,j}(x) = \int_{\mathbb{D}_1} \varphi_\varepsilon(z) \cdot N_j(x, \varepsilon z) \, ds(z), \tag{6.13} $$

where we have introduced the rescaled quantity:

$$ \varphi_\varepsilon(z) = \varepsilon^{d-1} (Ae(r_\varepsilon)n)(\varepsilon z) \in \tilde{H}^{-1/2}(\mathbb{D}_1)^d. $$

Step 2: We construct an integral equation characterizing the function $\varphi_\varepsilon(z)$.

As in the proof of Theorem 4.4, to achieve this, we repeat the above calculation, this time with the explicit Green's function $L(x, y)$ for the lower half space $H$, equipped with homogeneous Neumann boundary conditions on $\partial H$ (see Rem. 6.1), in place of the “difficult” Green's function $N(x, y)$ for the background equation (6.1). This yields the integral equation:

$$ \forall x \in \mathbb{D}_1, \quad \int_{\mathbb{D}_1} L(\varepsilon x, \varepsilon z) \varphi_\varepsilon(z) \, ds(z) = -u_0(\varepsilon x) + o(1). \tag{6.14} $$
Step 3: We use this integral equation to glean information about the function $\varphi_{\varepsilon}(z)$.

This task proceeds differently, depending on the space dimension:

* In 2d, using the expression (6.7) and (6.8) of $L(x, y)$, (6.14) rewrites:

$$ \frac{1 - \bar{\nu}}{\pi \mu} (\log \varepsilon) \int_{\mathbb{D}_1} \varphi_{\varepsilon}(y) \, ds(y) + K \varphi_{\varepsilon}(y) = -u_0(\varepsilon x) + o(1), $$

where $K: \tilde{H}^{-1/2}(\mathbb{D}_1)^d \to H^{1/2}(\mathbb{D}_1)^d$ is a bounded operator. We directly obtain from this equation that:

$$ \int_{\mathbb{D}_1} \varphi_{\varepsilon}(y) \, ds(y) = \frac{1}{|\log \varepsilon|} \frac{\pi \mu}{1 - \bar{\nu}} u_0(0) + o\left(\frac{1}{|\log \varepsilon|}\right), \quad (6.15) $$

which is the needed information for our purpose.

* In 3d, using the homogeneity of the kernel $L(x, y)$ in (6.5) and (6.6), (6.14) becomes:

$$ T_L \varphi_{\varepsilon} = -\varepsilon u_0(0) + o(\varepsilon), $$

where the integral operator $T_L$ is defined in (6.11). This immediately yields:

$$ \varphi_{\varepsilon} = \sum_{j=1}^{3} u_{0,j}(0) T_L^{-1} e_j, $$

and so:

$$ \int_{\mathbb{D}_1} \varphi_{\varepsilon} \, ds = M u_0(0), \quad (6.16) $$

involving the polarization tensor $M$ defined in (6.10).

Step 4: We pass to the limit in the representation formula (6.13).

Applying the Lebesgue dominated convergence theorem to the representation formula (6.13), we obtain:

$$ r_{\varepsilon,j}(x) = \left( \int_{\mathbb{D}_1} \varphi_{\varepsilon}(z) \, ds(z) \right) \cdot \left( N_j(x, 0) + o(1) \right), $$

and the desired result follows from the combination of this identity with (6.15) and (6.16). $\square$

## 6.4. Sensitivity of a quantity of interest with respect to the addition of a Dirichlet region

Let us consider the following quantity of interest, depending on the region $\Gamma_D$ of $\partial\Omega$ supporting the homogeneous Dirichlet boundary conditions of the problem (6.1):

$$ J(\Gamma_D) = \int_{\Omega} j(u_{\Gamma_D}) \, dx, $$

for a smooth function $j: \mathbb{R}^d \to \mathbb{R}$, satisfying suitable growth conditions, see (3.4). Here, $u_{\Gamma_D}$ denotes the solution to the boundary value problem (6.1).
**Corollary 6.4.** The perturbed value $J((\Gamma_D)_{x_0,\varepsilon})$ of $J(\Gamma_D)$, accounting for the replacement of the homogeneous Neumann boundary condition on $\omega_{x_0,\varepsilon} \subset \Gamma$ by a homogeneous Dirichlet boundary condition, has the following asymptotic expansion:

$$J((\Gamma_D)_{x_0,\varepsilon}) = J(\Gamma_D) + \frac{1}{|\log \varepsilon|} \frac{\pi \mu}{1-\bar{\nu}} u_0(x_0) \cdot p_0(x_0) + o\left(\frac{1}{|\log \varepsilon|}\right) \text{ if } d=2,$$

and

$$J((\Gamma_D)_{x_0,\varepsilon}) = J(\Gamma_D) + \varepsilon M u_0(x_0) \cdot p_0(x_0) + o(\varepsilon) \text{ if } d=3.$$

Here, the polarization tensor $M$ is defined by (6.10) and the adjoint state $p_0$ is the unique solution in $H^1(\Omega)^d$ to the boundary value problem:

$$\left\{
\begin{array}{ll}
-\mathrm{div}(Ae(p_0)) = -j'(u_0) & \text{in } \Omega, \\
p_0 = 0 & \text{on } \Gamma_D, \\
Ae(p_0)n = 0 & \text{on } \Gamma_N \cup \Gamma.
\end{array}
\right.
\qquad (6.17)$$

**Remark 6.5.** As in the case of electrostatics addressed in Section 4, multiple variations of the present study could be considered. For instance, one may be interested in accounting for the effect of the replacement of homogeneous Neumann boundary conditions by inhomogeneous Dirichlet or Neumann conditions, etc. Numerical examples associated to such variations are presented in Section 7.6.

## 7. NUMERICAL ILLUSTRATIONS

In this section, we apply our optimization method to several examples of shape and topology optimization problems related to regions bearing the boundary conditions of a physical problem. After providing a few general details about our numerical implementation in Section 7.1, we describe our treatment of integral equations in Section 7.2. Turning to physical applications, we first address in Section 7.3 the optimal design of an electroosmotic mixer, whose behavior is governed by the conductivity equation. An instance of such optimal design problems is then considered in the realm of acoustics in Section 7.4. Eventually, two instances of our boundary optimization framework are tackled in the physical setting of mechanical structures in Sections 7.5 and 7.6.

### 7.1. Description of the numerical framework

The practical realization of the general shape and topology optimization Algorithm 1 for the resolution of (P) deserves a few comments, see Algorithm 2 below for a more specific sketch of our numerical workflow.

Optimizing a region $G \subset \partial\Omega$ bearing boundary conditions raises antagonistic needs about the numerical representation of $G$. On the one hand, at each stage of the optimization process, one needs to solve one or several boundary value problems on $\Omega$, whose boundary conditions are supported by the region $G$. Ideally, this operation ought to be conducted on a mesh $\mathcal{T}$ of $\Omega$ whose associated surface mesh of $\partial\Omega$ encloses an explicit discretization of $G$, see Figure 6(a). On the other hand, the numerical representation of $G$ should allow to account for possibly dramatic deformations of this region from one optimization iteration to the next in a robust manner, including changes in its topology. This task is notoriously difficult to carry out under a meshed representation.

To reconcile these requirements, we use the level set based mesh evolution strategy introduced in [37, 88, 89], recently adapted to the case of regions on surfaces in [34]. Briefly, two complementary representations of $G$ are available at each optimization iteration:
<figure><img src="image_9.png" /><figcaption>FIGURE 6. (a) Meshed representation of a region G (in blue) of the boundary of a 3d domain Ω; (b) Corresponding representation of G as the negative subdomain of a level set function defined at the vertices of the surface mesh S of ∂Ω induced by the mesh T of Ω.</figcaption></figure>

• (Meshed representation) The domain Ω is equipped with a high-quality simplicial mesh T, i.e. made of triangles in 2d, tetrahedra in 3d. Its boundary S, (a line mesh if d = 2, a surface triangulation if d = 3), encloses an explicit discretization of the region G as a submesh S$_{int}$, see Figure 6(a).

• (Level set representation) The region G is described as the negative subset of a “level set” function φ : ∂Ω → R, that is:

$$
\forall x \in \partial\Omega, \quad \begin{cases} \phi(x) < 0 & \text{if } x \in G, \\ \phi(x) = 0 & \text{if } x \in \Sigma, \\ \phi(x) > 0 & \text{otherwise.} \end{cases}
$$

In practice φ is discretized at the vertices of a surface mesh S of ∂Ω, see [90] for the seminal reference
about the level set method, the articles [91–94] about its introduction in the field of shape and topol-
ogy optimization, and the books [95, 96] for further details, and Figure 6(b) for an illustration of this
representation.

Every numerical operation in the process is performed by using the most suitable of these two representations of G. Notably, the resolution of boundary value problems on Ω whose boundary conditions are supported on G and ∂Ω \ G takes advantage of the meshed representation of G: the Finite Element Method can be conveniently applied on the mesh T of Ω, and in principle, any numerical solver could be used to this end, in a black-box fashion. On the other hand, the deformation of G according to a tangential velocity field V : ∂Ω → R d (which is, in our particular context, a descent direction for the considered optimization problem (P)) for a (pseudo-) time period (0, τ) is better accounted for under a level set representation φ of G: it is captured by solving the following advection-like equation on the surface mesh S of ∂Ω [95, 96]:

$$
\left\{
\begin{array}{ll}
\frac{\partial \psi}{\partial t}(t, x) + V(t, x) \cdot \nabla_{\partial\Omega} \psi(t, x) = 0 & \text{for } t \in (0, \tau), \ x \in \partial\Omega, \\
\psi(0, x) = \phi(x) & \text{for } x \in \partial\Omega.
\end{array}
\right.
\qquad (7.1)
$$
This representation makes it even simpler to account for the addition of a tiny surface disk $\omega_{x_0, \varepsilon} \subset \partial\Omega$ to $G$: one replaces $\phi$ with the level set function

$$ \min(\phi(x), |x - x_0| - \varepsilon), \quad x \in \partial\Omega. \qquad (7.2) $$

Our level set based mesh evolution strategy crucially hinges on a set of efficient numerical algorithms allowing to switch from one of these representations to the other. Notably, a level set function $\phi$ for a region $G \subset \partial\Omega$ available under meshed representation is generated by calculating the signed distance function $d_G^{\partial\Omega}$ to $G$ (see Def. 3.2). Conversely, the creation of a meshed representation of $G$ from the datum of a level set function $\phi: \partial\Omega \to \mathbb{R}$ relies on suitable mesh modification techniques, described in e.g. [97–99].

<table><caption>**Algorithm 2:** Mesh evolution algorithm for the optimal design of a region $G \subset \partial\Omega$.</caption><tr><td>Initialization:</td><td>Simplicial mesh T<sup>0</sup> of Ω whose associated surface triangulation S<sup>0</sup> features a submesh of the initial region G<sup>0</sup> ⊂ ∂Ω.</td></tr><tr><td>for n = 0, ..., until convergence do</td><td>(1) Calculate the state function u<sub>G<sup>n</sup></sub> on the mesh T<sup>n</sup>.</td></tr><tr><td></td><td>(2) Calculate the adjoint state p<sub>G<sup>n</sup></sub> on T<sup>n</sup>.</td></tr><tr><td></td><td>(3) Calculate the signed distance function φ<sup>n</sup> to G<sup>n</sup> on the surface mesh S<sup>n</sup>.</td></tr><tr><td></td><td>if n mod. n<sub>top</sub> = 0 then</td></tr><tr><td></td><td>(4) Calculate the topological derivative d<sub>T</sub>J(G<sup>n</sup>)(x) at the vertices x of S<sup>n</sup>.</td></tr><tr><td></td><td>(5) Find the point x<sub>0</sub> ∈ ∂Ω \ G<sup>n</sup> where d<sub>T</sub>J(G<sup>n</sup>) takes the largest negative value.</td></tr><tr><td></td><td>(6) Obtain the new level set function φ<sup>n+1</sup> from φ<sup>n</sup> by (7.2), for a “small” parameter ε > 0.</td></tr><tr><td></td><td>else</td></tr><tr><td></td><td>(4) Calculate the shape derivative J'(G<sup>n</sup>)(θ) of J, or that J'<sub>ε</sub>(G<sup>n</sup>)(θ) of an approximate version, in the spirit of Section 3.3.</td></tr><tr><td></td><td>(5) Find a descent direction θ<sup>n</sup>.</td></tr><tr><td></td><td>(6) Select a time step τ<sup>n</sup> > 0 and calculate the solution ψ(t, x) to the advection equation (7.1) with velocity V = θ<sup>n</sup> on the surface mesh S<sup>n</sup> over the interval (0, τ<sup>n</sup>); set φ<sup>n+1</sup> = ψ(τ<sup>n</sup>, ·).</td></tr><tr><td></td><td>end if</td></tr><tr><td></td><td>(7) Modify T<sup>n</sup> into a mesh T<sup>n+1</sup> of Ω<sup>n+1</sup> whose surface S<sup>n+1</sup> contains a submesh S<sup>n+1</sup><sub>int</sub> of G<sup>n+1</sup>.</td></tr><tr><td></td><td>end for</td></tr><tr><td></td><td>return Mesh T<sup>n</sup> of Ω, whose surface part S<sup>n</sup> contains a submesh S<sup>n</sup><sub>int</sub> of the optimized region G<sup>n</sup>.</td></tr></table>

The operations involved in this framework are carried out by black-box uses of several open-source libraries developed in our previous works. The boundary value problems at stake are assembled and solved with the Rodin open-source library [100]. The signed distance function $d_G$ to a region $G \subset \partial\Omega$ is computed using the Fast Marching Method [101, 102], implemented in the ISCD Mshdist software [103]; the solution of the level set evolution equation (7.1) on the mesh $S$ of $\partial\Omega$ relies on the method of characteristics, implemented in the ISCD Advection software [104]. Eventually, the mesh modification operations of our framework leverage the mmg library [97, 98]. We refer to [105] for an educational implementation of this framework, in the context of the optimization of a “bulk” shape $\Omega \subset \mathbb{R}^d$. All our numerical experiments are conducted on a regular laptop Apple MacBookPro 18,3 (M1 Pro chip) with 10 cores and 16 GB of memory.

## 7.2. Numerical resolution of boundary integral equations

According to Sections 4, 5 and 6, the determination of the asymptotic sensitivity of the solution to a boundary value problem with respect to the addition of a “small” disk $\omega_{x_0, \varepsilon}$ to a region $G \subset \partial\Omega$ supporting its boundary
conditions raises an integral equation, of the form:

$$ \text{Search for } \varphi \in \tilde{H}^{-1/2}(\mathbb{D}_1) \text{ s.t. } T_L\varphi(x) = f(x), \quad x \in \mathbb{D}_1, \qquad (\mathcal{B}) $$

where the unknown $\varphi$ is scalar-valued for simplicity of the presentation, and the integral operator $T_L$ : $\tilde{H}^{-1/2}(\mathbb{D}_1) \to H^{1/2}(\mathbb{D}_1)$ is defined by:

$$ T_L\varphi(x) = \int_{\mathbb{D}_1} L(x, z)\varphi(z) \, ds(z). \qquad (7.3) $$

It involves a homogeneous kernel $L$ of the form of those in (4.4), (5.2), (6.2), (6.5) and (6.6), and a known source term $f \in H^{1/2}(\mathbb{D}_1)$. In fortunate cases, as in electrostatics, the solution to this equation can be computed analytically, leading to a fully explicit asymptotic formula for the state $u_\epsilon$, as in electrostatics, see Theorem 4.4. However, in several realistic situations, no closed form expression of the solution is available and the latter has to be computed numerically, see notably the case of the linear elasticity system in Theorem 6.2.

The numerical resolution of integral equations is extensively discussed in the numerical analysis literature [106–109]. However, the particular situation of an equation set on a “screen”, i.e. an open subregion of the total boundary of a domain, is much less classical, see Section 4.1.11 in [110] for a related discussion. For completeness, we briefly discuss the use of the Boundary Element Method implemented in our work to solve such equations in the case where the space dimension $d$ equals 3, referring to [110, 111] for details.

### 7.2.1. The discrete Galerkin approximation of (B)

Like the “classical” Finite Element Method [112, 113], the Boundary Element Method is based on a variational formulation of the considered equation (B). The latter is obtained by multiplying (B) by an arbitrary test function $\psi \in \tilde{H}^{-1/2}(\mathbb{D}_1)$ and integrating:

$$ \text{Search for } \varphi \in \tilde{H}^{-1/2}(\mathbb{D}_1) \text{ s.t. } \forall \psi \in \tilde{H}^{-1/2}(\mathbb{D}_1), \quad a(\varphi, \psi) = \ell(\psi), \\ \text{where } a(\varphi, \psi) := \int_{\mathbb{D}_1} T_L\varphi(x)\psi(x) \, ds(x), \quad \text{and } \ell(\psi) := \int_{\mathbb{D}_1} f(x)\psi(x) \, ds(x). \qquad (\mathcal{G}) $$

A discrete approximation of the solution $\varphi$ to this problem is then sought in a finite-dimensional subspace of $\tilde{H}^{-1/2}(\mathbb{D}_1)$. As regards the latter, we choose the space $V_S$ of $\mathbb{P}_1$ Lagrange elements on a given surface mesh $S$ of $\mathbb{D}_1$ made of $K$ non overlapping and conforming triangles $T_1, \dots, T_K$, and $N$ vertices $a_1, \dots, a_N$:

$$ V_S = \{ v : \mathbb{D}_1 \to \mathbb{R}, \text{ is continuous and } v|_{T_k} \text{ is affine}, k = 1, \dots, K \}. $$

As usual, a basis of $V_S$ is given by the collection $\{\varphi_n\}_{n=1,\dots,N}$ of functions defined by the relations:

$$ \text{For all } n = 1, \dots, N, \text{ } \varphi_n \text{ is the unique function in } V_S \text{ s.t. } \varphi_n(a_m) = \begin{cases} 1 & \text{if } n = m, \\ 0 & \text{otherwise,} \end{cases} \quad m = 1, \dots, N. $$

For any $\varphi, \psi \in V_S$, we decompose the integral $a(\varphi, \psi)$ as:

$$ a(\varphi, \psi) = \sum_{k=1}^{K} \sum_{l=1}^{K} a_{kl}(\varphi, \psi), \text{ where } a_{kl}(\varphi, \psi) := \int_{T_k} \varphi(x) \left( \text{p.v.} \int_{T_l} L(x, z) \psi(z) \, ds(z) \right) ds(x), \qquad (7.4) $$

where p.v. stands for the Cauchy principal value, and all the integrals in the above sum can be proved to be well-defined.
### 7.2.2. Decomposition of singular integrals through relative coordinates

The evaluation of $a(\varphi, \psi)$ for given functions $\varphi, \psi \in V_S$ via the decomposition (7.4) leaves us with the task of approximating the integrals $a_{kl}(\varphi, \psi)$, $k, l = 1, \dots, K$. To achieve this, we express each triangle $T_k \in S$ as the image $T_k = A_k(\hat{T})$ of the reference simplex

$$ \hat{T} := \{(x_1, x_2) \in \mathbb{R}^2 \text{ s.t. } 0 \le x_1 \le 1 \text{ and } 0 \le x_2 \le 1 - x_1\} \subset \mathbb{R}^2 $$

via a suitable affine mapping $A_k : \mathbb{R}^2 \to \mathbb{R}^3$. For each triangle, a local basis $\{b_{k,i}\}_{1 \le i \le 3}$ of affine functions on $T_k$ is given by:

$$ b_{k,i} = \hat{b}_i \circ A_k, $$

where the functions $\hat{b}_i$ constitute the basis of the first-order polynomial functions on $\hat{T}$, namely:

$$ \forall x = (x_1, x_2) \in \hat{T}, \quad \hat{b}_1(x) = x_1, \quad \hat{b}_2(x) = x_2, \quad \hat{b}_3(x) = 1 - x_1 - x_2. $$

Let $\varphi, \psi$ be given functions in $V_S$, and let us decompose $\varphi$ and $\psi$ on the basis $\{\varphi_n\}_{n=1,\dots,N}$ as:

$$ \varphi(x) = \sum_{n=1}^{N} \alpha_n \varphi_n(x), \quad \psi(x) = \sum_{n=1}^{N} \beta_n \varphi_n(x) \quad x \in \mathbb{D}_1. \quad (7.5) $$

Pulling back the integrals in $a_{kl}(\varphi, \psi)$ onto the reference triangle $\hat{T}$, we obtain, for $k, l = 1, \dots, K$:

$$ a_{kl}(\varphi, \psi) = \sum_{i,j=1}^{3} \alpha_{n(k,i)} \beta_{n(l,j)} \hat{I}_{klij}, \quad \text{where } \hat{I}_{klij} := \int_{\hat{T}} \text{p.v.} \int_{\hat{T}} \hat{L}_{klij}(x,y) \, dy \, dx \quad \text{and} \\ \hat{L}_{klij}(x,y) := \hat{b}_j(x) \hat{b}_i(y) \text{ Jac}(A_k) \text{Jac}(A_l) L(A_k(x), A_l(y)). \quad (7.6) $$

Here, for $k = 1, \dots, K$ and $i = 1, 2, 3$, we have denoted by $n(k, i) \in \{1, \dots, N\}$ the global index of the $i^{th}$ vertex in the $k^{th}$ triangle $T_k$. The Jacobian determinants of the mappings $A_k : \hat{T} \subset \mathbb{R}^2 \to \mathbb{R}^3$ are defined by $Jac(A_k) := \det(\nabla A_k^T \nabla A_k)^{1/2}, k = 1, \dots, K$.

Suitable quadrature formulas are then used to approximate these integrals while taking into account the weak singularity of the kernel $L(x, y)$. More precisely, we rely on the regularized coordinate versions of the quantities $\hat{I}_{klij}$ proposed in Section 5.2.4 of [110] which can be approximated via standard Gaussian quadrature methods, see also [114–117] about this issue.

### 7.2.3. Regularization of the variational problem

Unfortunately, the direct resolution of (G) with the above strategy is plagued with numerical artifacts. These are mainly due to the nature of the considered surface $\mathbb{D}_1$, which is open, contrary to the common setting of integral equations where it is a closed boundary. In such a situation, it is indeed well-known that the solution $\varphi$ blows up near the boundary $\partial \mathbb{D}_1$, see [118–120] about this general behavior, and the numerical experiments in [121] about the induced numerical instabilities.

To remedy this issue, we resort to a simple, heuristic procedure. We approximate (G) by the following regularized problem:

$$ \text{Search for } \varphi_{\eta} \in H^1(\mathbb{D}_1) \text{ s.t. } \forall \psi \in H^1(\mathbb{D}_1), \\ \eta \int_{\mathbb{D}_1} \nabla_{\mathbb{D}_1} \varphi_{\eta}(x) \cdot \nabla_{\mathbb{D}_1} \psi(x) \, ds(x) + \int_{\mathbb{D}_1} T_L \varphi_{\eta}(x) \psi(x) \, ds(x) = \int_{\mathbb{D}_1} f(x) \psi(x) \, ds(x). \quad (7.7) $$
<figure><img src="image_10.png" /><figcaption>FIGURE 7. (a) Exact, and (b) Approximate solution (using $\eta = 1e-5$, $h = 0.039$) to the boundary integral equation (B) in the case of the conductivity equation considered in Section 7.2.4; a threshold is applied to both functions near $\partial D_1$ for vizualization, because of the blow up of the exact solution there.</figcaption></figure>

where we have slightly perturbed the main bilinear form $a(\cdot, \cdot)$ in (G) with the addition of a regularizing elliptic term, penalized by a “small” parameter $\eta > 0$.

In the next section, we test this regularization method on an example.

### 7.2.4. An example in the setting of the conductivity equation

Let us consider the resolution of (B) in the context of three-dimensional electrostatics discussed in Section 4.2: $L(x, y)$ is the fundamental solution (4.4) of the Laplace equation in free space. We set the right-hand side to $f = 1$, in which case the solution $\varphi$ has the analytical expression (4.26), see Figure 7(a).

We solve the approximate problem (7.7) for several values of the mesh size $h$ (i.e. the average length of an edge in the mesh) and of the regularization parameter $\eta$, see Figure 7(b) for an illustration. The accuracy of the resolution is evaluated in terms of the following three measures of error:

(1) The squared residual

$$ \mathcal{R}(h, \eta) := \int_{\mathcal{T}_h} \left( \frac{1}{4\pi} \int_{\mathcal{T}_h} \frac{1}{|x-y|} \varphi_{\eta}(y) \, ds(y) - 1 \right)^2 ds(x) \quad (7.8) $$

evaluates how far the approximate solution $\varphi_{\eta}$ is from satisfying the original equation (B).

(2) The squared error over the mean value

$$ \mathcal{A}(h, \eta) := \left( \int_{\mathcal{T}_h} \varphi_{\eta}(x) \, ds(x) - 8 \right)^2 $$

measures the discrepancy between the mean values of the approximate and exact solutions.
<figure><img src="image_11.png" /><figcaption>FIGURE 8. (a) Plot of the numerical resolution error $\mathcal{R}(h, \eta)$; (b) Plot of the numerical error of the average $\mathcal{A}(h, \eta)$; (c) Plot of the numerical squared approximation error $\mathcal{E}(h, \eta)$.</figcaption></figure>

(3) The squared error

$$ \mathcal{E}(h, \eta) := \int_{\{x \in \mathbb{D}_1, |x| < 0.9\}} |\varphi_{\eta}(x) - \varphi(x)|^2 \mathrm{d}s(x) $$

quantifies the error between the approximate and exact solutions. Note that this error is measured on the subset $\{x \in \mathbb{D}_1, |x| < 0.9\}$ lying at a distance from the boundary $\partial \mathbb{D}_1$ where the exact solution $\varphi$ blows up.

We solve (7.7) and evaluate these three quantities on 50 uniformly spaced values in the intervals $0.00001 < \eta < 0.1$ and $0.1 < h < 1$, respectively. The results are depicted in terms of both parameters in Figure 8. As expected, these measures of the error tend to 0 as $\eta$ and $h$ go to zero. We also observe that the regularization parameter $\eta$ should not be taken too small when compared to the mesh size $h$, see notably Figures 8(b,c). This effect is quite understandable, as the regularizing term in (7.7) expresses a diffusion over a zone with radius $\eta^{1/2}$ – an operation which is numerically unstable when the mesh is too coarse with respect to this length scale.

### 7.2.5. Extension to vector-valued boundary integral equations

The technique presented in the previous sections can be extended to deal with vector-valued integral equations. To set ideas, let us now consider an integral equation of the form (B), involving the operator $T_L: \tilde{H}^{-1/2}(\mathbb{D}_1)^3 \to H^{1/2}(\mathbb{D}_1)^3$ defined by:

$$ \forall \varphi \in \tilde{H}^{-1/2}(\mathbb{D}_1)^3, \quad T_L\varphi(x) = \int_{\mathbb{D}_1} L(x, z)\varphi(z) \mathrm{d}s(z), \quad x \in \mathbb{D}_1, \qquad (7.9) $$

where the kernel $L(x, z)$ is now a $3 \times 3$ matrix-valued function and the right-hand side $f$ belongs to $H^{1/2}(\mathbb{D}_1)^3$.
An identical strategy to that outlined in Sections 7.2.1, 7.2.2 and 7.2.3 allows to treat the present case: a variational formulation is introduced for (B), and the left-hand side of the latter is decomposed onto the triangles $T_k$ of the surface triangulation $\mathcal{S}$ of $\mathbb{D}_1$, see (7.6). In this case the function $\hat{L}_{klij} : \mathbb{D}_1 \times \mathbb{D}_1 \to \mathbb{R}^d$ defined in (7.6) is written in terms of the matrix product:

$$ \hat{L}_{klij}(x, y) = \mathrm{Jac}(A_k) \mathrm{Jac}(A_l) \left[ L(A_k(x), A_l(y)) \hat{b}_i(y) \right] \cdot \hat{b}_j(x), \quad (7.10) $$

and the same quadrature rules as in the scalar case can be used for the calculation of the integrals $\hat{I}_{klij}$. In this situation also, the numerical resolution is riddled with numerical artifacts, caused by the blow up of the solution near the boundary $\partial\mathbb{D}_1$: a similar regularization strategy as that proposed in Section 7.2.3 is used to guarantee its stability.

Let us discuss an example of this methodology in the linear elasticity context of Section 6: we aim to calculate the entries of the polarization tensor $M$ in (6.10), featured in the asymptotic expansion of the perturbed displacement by the addition of a “small” surface disk to the region of $\partial\Omega$ supporting Dirichlet boundary conditions, see Theorem 6.2. We consider the equation (B) where the operator $T_L$ is associated to the 3d Mindlin kernel $L(x, y)$ with physical parameters $\mu = 67.5676$ and $\nu = 0.48$, see (6.5) and (6.6). The right-hand side $f$ is one of the three coordinate vectors $e_i \in \mathbb{R}^3$, $i = 1, 2, 3$. The regularizing bilinear form chosen in this situation is that associated to the linearized elasticity system, i.e. we solve the following variational problem:

Search for $\varphi_\eta \in H^1(\mathbb{D}_1)^3$ s.t. $\forall \psi \in H^1(\mathbb{D}_1)^3$,

$$ \eta \int_{\mathbb{D}_1} A e_{\mathbb{D}_1}(\varphi_\eta) : e_{\mathbb{D}_1}(\psi) \, ds(x) + \int_{\mathbb{D}_1} T_L \varphi_\eta(x) \cdot \psi(x) \, ds(x) = \int_{\mathbb{D}_1} f(x) \cdot \psi(x) \, ds(x). \quad (7.11) $$

We solve this equation for $i = 1, 2, 3$ and various values of the mesh size $h$ and regularization coefficient $\eta$. For each of these parameters, we take 50 uniformly spaced values within the intervals $0.5 < h < 1$ and $0.00001 < \eta < 0.1$, and we evaluate the accuracy of the solution by computing the resolution error $\mathcal{R}(h, \eta)$ in (7.8). The error plots are presented in Figure 9, see also Figure 10 for an illustration of the solution.

The error behaves as in the previous Section 7.2.4. As both parameters $h$ and $\eta$ decrease, so does the quantity $\mathcal{R}(h, \eta)$. Also, as expected due to the symmetry of the physical configuration, the solutions associated to the right-hand sides $e_1$ and $e_2$ are identical up to a $90^\circ$ rotation. On the other hand, the solution associated to the right-hand side $e_3$ differs from the previous two: it is totally symmetric, see Figure 7.

## 7.3. Optimization of the repartition of the cathode and anode regions on the boundary of a direct current electroosmotic mixer

Our first example of optimal design problems featuring regions supporting boundary condition regions arises in the field of microfluidics, which deals with the handling of very small volumes of fluid, ranging between $10^{-18}$ and $10^{-9}$ L. This field has been notably active over the past decade. It notably heralds decisive advances in analytical chemistry, molecular biology and biomedical engineering (with the design of biochips and DNA micro-arrays and the possibility to realize electrophoresis and liquid chromatography for proteins and DNA), in optics, where it opens the way to the design of microlens arrays, etc. We refer to e.g. [122, 123] for comprehensive introductions to this subject and its challenges.

One basic operation of particular interest in microfluidics consists in mixing two fluids by electroosmosis [124–128]. This process typically takes place in a so-called electroosmotic mixer: in broad outline, two fluids are injected through an entrance channel, and an electric field is imposed inside the ring-shaped mixing chamber thanks to a suitable placement of anodes (positively charged electrodes) and cathodes (negatively charged electrodes) on its boundary. The induced Coulomb force electrically actuates the charged particles and ions and triggers the pumping of electrolytic fluid through drag forces. The more intense the electric field inside the
<figure><img src="image_12.png" /><figcaption>FIGURE 9. Plot of the error $\mathcal{R}(h, \eta)$ in (7.8) committed in the solution of the integral equation (7.11) when the right-hand side $f$ equals (a) $e_1$; (b) $e_2$; (c) $e_3$; the coordinates axes for $h$, $\eta$ and $\mathcal{R}(h, \eta)$ are expressed in logarithmic scale.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 10. Plot of the solution to the integral equation (7.11) in the linear elasticity setting of Section 7.2.5 when the right-hand side $f$ equals (a) $e_1$; (b) $e_2$; (c) $e_3$; the numerical parameters of the computation are $h = 0.039$, $\eta = 0.00001$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIGURE 11. Illustration of the electrosmotic mixer considered in Section 7.3: two fluids enter via the inflow channel, get mixed in the central chamber under the action of the difference between the voltage potentials applied on the anode and cathode regions; the resulting mixture exits through the outflow channel.</figcaption></figure>

chamber, the more efficient the mixing process. The resulting mixture finally leaves the device through an exit channel, see Figure 11.

In this section, we seek to optimize the design of such an electroosmotic mixer. Our study relies on a simplified version of the physical model presented in [129]. We aim to optimize the distribution of anodes and cathodes on the boundary of the mixing chamber so as to generate an electric field with maximum amplitude inside the device, and thereby improve the efficiency of the mixer.

### 7.3.1. Presentation of the optimization problem

The optimal design problem of this section arises in the mathematical setting of the conductivity equation, see Section 3.1. The electroosmotic mixer is represented by a bounded domain $\Omega \subset \mathbb{R}^3$, whose boundary is decomposed into three disjoint parts:

$$ \partial\Omega = \overline{\Gamma_A} \cup \overline{\Gamma_C} \cup \overline{\Gamma}. $$

In this representation,

* The region $\Gamma_C$ is the cathode, where the voltage potential is set to 0;

* The region $\Gamma_A$ is the anode, where a value $u_{in} > 0$ is imposed for the potential;

* The device is perfectly insulated from the outside in the remaining region $\Gamma$.

We also let $\Sigma_C := \partial\Gamma_C$, $\Sigma_A := \partial\Gamma_A$, and we assume that $\Gamma_C$ and $\Gamma_A$ are “well-separated” from each other, in the sense that:

$$ \mathrm{dist}(\Gamma_C, \Gamma_A) > \delta \text{ for some fixed value } \delta > 0. $$
In this setting, we consider the following optimal design problem:

$$
\begin{equation}
\begin{array}{ll}
\displaystyle \min_{\Gamma_C, \Gamma_A \subset \partial\Omega} & J(\Gamma_C, \Gamma_A) + \ell_A \operatorname{Area}(\Gamma_A) + \ell_C \operatorname{Area}(\Gamma_C) + m_A \operatorname{Cont}(\Gamma_A) + m_C \operatorname{Cont}(\Gamma_C), \\
\text{where } & J(\Gamma_C, \Gamma_A) = \int_{\Omega} j(\nabla u_{\Gamma_C, \Gamma_A}) \, dx.
\end{array}
\tag{7.12}
\end{equation}
$$

Here, $\ell_C$, $\ell_A > 0$ are penalization parameters for the surface areas of the cathode $\Gamma_C$ and anode $\Gamma_A$, and $m_C > 0$ and $m_A > 0$ are parameters that penalize their contours, see the definitions in (2.5). The function $j: \mathbb{R}^3 \to \mathbb{R}$ satisfies the growth conditions (3.4), and $u_{\Gamma_C, \Gamma_A} \in H^1(\Omega)$ is the solution to:

$$
\left\{
\begin{array}{ll}
-\operatorname{div}(\gamma \nabla u_{\Gamma_C, \Gamma_A}) = 0 & \text{in } \Omega, \\
u_{\Gamma_C, \Gamma_A} = 0 & \text{on } \Gamma_C, \\
u_{\Gamma_C, \Gamma_A} = u_{\text{in}} & \text{on } \Gamma_A, \\
\gamma \frac{\partial u_{\Gamma_C, \Gamma_A}}{\partial n} = 0 & \text{on } \Gamma.
\end{array}
\right.
\qquad (7.13)
$$

7.3.2. Approximation of the functional J(Γ<sub>C</sub>, Γ<sub>A</sub>) and of its shape derivative

As noted in Section 3, the solution $u_{\Gamma_C, \Gamma_A}$ to the boundary value problem (7.13) is weakly singular near the contours $\Sigma_C$ and $\Sigma_A$, where the boundary conditions change types. To overcome the difficulties induced by this behavior in the calculation of the functional $J(\Gamma_C, \Gamma_A)$ and its shape derivative, we replace it with the following counterpart:

$$
J_{\varepsilon}(\Gamma_C, \Gamma_A) = \int_{\Omega} j(\nabla u_{\Gamma_C, \Gamma_A, \varepsilon}) \, dx, \quad (7.14)
$$

defined in terms of the solution $u_{\Gamma_C, \Gamma_A, \varepsilon} \in H^1(\Omega)$ to an approximate version of (7.13) where the sharp transitions between Dirichlet and Neumann boundary conditions are smoothed into a Robin boundary condition:

$$
\left\{
\begin{array}{ll}
-\operatorname{div}(\gamma \nabla u_{\Gamma_C, \Gamma_A, \varepsilon}) = 0 & \text{in } \Omega, \\
\gamma \frac{\partial u_{\Gamma_C, \Gamma_A, \varepsilon}}{\partial n} + (h_{\Gamma_C, \varepsilon} + h_{\Gamma_A, \varepsilon}) u_{\Gamma_C, \Gamma_A, \varepsilon} - h_{\Gamma_A, \varepsilon} u_{\text{in}} = 0 & \text{on } \partial\Omega.
\end{array}
\right.
\qquad (7.15)
$$

In this formulation, the coefficient $h_{\Gamma_C, \varepsilon} : \partial\Omega \to \mathbb{R}$ is defined by:

$$
\forall x \in \partial\Omega, \quad h_{\Gamma_C, \varepsilon}(x) = \frac{1}{\varepsilon} h\left(\frac{d_{\Gamma_C}^{\partial\Omega}(x)}{\varepsilon}\right), \qquad (7.16)
$$

where $d_{\Gamma_C}^{\partial\Omega}$ is the geodesic signed distance function to $\Gamma_C$ (see Sect. 3.2) and the transition profile $h \in C^\infty(\mathbb{R})$ satisfies (3.12); the function $h_{\Gamma_A, \varepsilon}$ is defined analogously.

The shape derivative of the functional $J_{\varepsilon}(\Gamma_C, \Gamma_A)$ is the subject of the next proposition. The proof, which is fairly similar to that of Proposition 3.11, is omitted for brevity. We refer to [121] for a detailed calculation.

**Proposition 7.1.** The criterion $J_{\varepsilon}(\Gamma_C, \Gamma_A)$ is shape differentiable, and its shape derivative reads, for any tangential deformation $\theta$ (i.e. $\theta \cdot n \equiv 0$ on $\partial\Omega$):

$$
\begin{align*}
J'_{\varepsilon}(\Gamma_C, \Gamma_A)(\theta) &= \frac{1}{\varepsilon^2} \int_{\partial\Omega} h'\left(\frac{d_{\Gamma_A}^{\partial\Omega}(x)}{\varepsilon}\right) \theta(\pi_{\Sigma_A}(x)) \cdot n_{\Sigma_A}(\pi_{\Sigma_A}(x))(u_{\text{in}} - u_{\Gamma_C, \Gamma_A, \varepsilon}(x)) p_{\Gamma_C, \Gamma_A, \varepsilon}(x) \,ds(x) \\
&\quad - \frac{1}{\varepsilon^2} \int_{\partial\Omega} h'\left(\frac{d_{\Gamma_C}^{\partial\Omega}(x)}{\varepsilon}\right) \theta(\pi_{\Sigma_C}(x)) \cdot n_{\Sigma_C}(\pi_{\Sigma_C}(x)) u_{\Gamma_C, \Gamma_A, \varepsilon}(x) p_{\Gamma_C, \Gamma_A, \varepsilon}(x) \,ds(x),
\end{align*}
$$
where the adjoint state $p_{\Gamma_C, \Gamma_A, \varepsilon} \in H^1(\Omega)$ is the solution to:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(\gamma \nabla p_{\Gamma_C, \Gamma_A, \varepsilon}) = \mathrm{div}(j'(\nabla u_{\Gamma_C, \Gamma_A, \varepsilon})) & \text{in } \Omega, \\ \gamma \frac{\partial p_{\Gamma_C, \Gamma_A, \varepsilon}}{\partial n} + (h_{\Gamma_C, \varepsilon} + h_{\Gamma_A, \varepsilon}) p_{\Gamma_C, \Gamma_A, \varepsilon} = -j'(\nabla u_{\Gamma_C, \Gamma_A, \varepsilon}) \cdot n & \text{on } \partial\Omega. \end{array} \right. $$

An approximate counterpart for this formula which conveniently simplifies the numerical implementation is obtained by arguing along the lines of Section 3.4:

$$ J'_{\varepsilon}(\Gamma_C, \Gamma_A)(\theta) = \frac{1}{\varepsilon} \int_{\Sigma_A} (u_{\mathrm{in}} - u_{\Gamma_C, \Gamma_A, \varepsilon}) p_{\Gamma_C, \Gamma_A, \varepsilon} \theta \cdot n_{\Sigma_A} d\ell - \frac{1}{\varepsilon} \int_{\Sigma_C} u_{\Gamma_C, \Gamma_A, \varepsilon} p_{\Gamma_C, \Gamma_A, \varepsilon} \theta \cdot n_{\Sigma_C} d\ell. \quad (7.17) $$

### 7.3.3. The topological derivative of $J(\Gamma_C, \Gamma_A)$

Let us now deal with the topological sensitivity of the criterion $J(\Gamma_C, \Gamma_A)$ when the homogeneous Neumann boundary condition is replaced by a (homogeneous or inhomogeneous) Dirichlet boundary condition on a “small” surface disk $\omega_{x_0, \varepsilon} \subset \Gamma$. An elementary adaptation of the proofs of Theorem 4.4 and Corollary 4.5 (see also Rem. 4.6) yields the following result.

**Proposition 7.2.** Let $\Gamma_C, \Gamma_A$ be disjoint regions of the smooth boundary $\partial\Omega$ as in Section 7.3.1, and let $x_0 \in \Gamma$. Then,

(i) The perturbed criterion $J((\Gamma_C)_{x_0, \varepsilon}, \Gamma_A)$, that accounts for the addition of the surface disk $\omega_{x_0, \varepsilon} \subset \Gamma$ to $\Gamma_C$, has the following asymptotic expansion:

$$ J((\Gamma_C)_{x_0, \varepsilon}, \Gamma_A) = J(\Gamma_C, \Gamma_A) - \frac{\pi}{|\log \varepsilon|} \gamma(x_0) u_{\Gamma_C, \Gamma_A}(x_0) p_{\Gamma_C, \Gamma_A}(x_0) + o\left(\frac{1}{|\log \varepsilon|}\right) \quad \text{if } d=2, $$

and

$$ J((\Gamma_C)_{x_0, \varepsilon}, \Gamma_A) = J(\Gamma_C, \Gamma_A) - 4\varepsilon \gamma(x_0) u_{\Gamma_C, \Gamma_A}(x_0) p_{\Gamma_C, \Gamma_A}(x_0) + o(\varepsilon) \quad \text{if } d=3. $$

(ii) The perturbed criterion $J(\Gamma_C, (\Gamma_A)_{x_0, \varepsilon})$, that accounts for the addition of $\omega_{x_0, \varepsilon} \subset \Gamma$ to $\Gamma_A$, has the following asymptotic expansion:

$$ J(\Gamma_C, (\Gamma_A)_{x_0, \varepsilon}) = J(\Gamma_C, \Gamma_A) + \frac{\pi}{|\log \varepsilon|} \gamma(x_0) (u_{\mathrm{in}} - u_{\Gamma_C, \Gamma_A}(x_0)) p_{\Gamma_C, \Gamma_A}(x_0) + o\left(\frac{1}{|\log \varepsilon|}\right) \quad \text{if } d=2, $$

and

$$ J(\Gamma_C, (\Gamma_A)_{x_0, \varepsilon}) = J(\Gamma_C, \Gamma_A) + 4\varepsilon \gamma(x_0) (u_{\mathrm{in}} - u_{\Gamma_C, \Gamma_A}(x_0)) p_{\Gamma_C, \Gamma_A}(x_0) + o(\varepsilon) \quad \text{if } d=3. $$

In all the above formulas, the adjoint state $p_{\Gamma_C, \Gamma_A} \in H^1(\Omega)$ is the solution to:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(\gamma \nabla p_{\Gamma_C, \Gamma_A}) = \mathrm{div}(j'(\nabla u_{\Gamma_C, \Gamma_A})) & \text{in } \Omega, \\ \gamma \frac{\partial p_{\Gamma_C, \Gamma_A}}{\partial n} = -j'(\nabla u_{\Gamma_C, \Gamma_A}) \cdot n & \text{on } \Gamma, \\ p_{\Gamma_C, \Gamma_A} = 0 & \text{on } \Gamma_C \cup \Gamma_A. \end{array} \right. $$
<figure><img src="image_15.png" /><figcaption>FIGURE 12: Different views of the initial mesh $\mathcal{T}^0$ of the electroosmotic mixer considered in Section 7.3. The red and pink regions correspond to the non optimizable components of the anode and cathode regions, and the green region represents the inflow region.</figcaption></figure>

### 7.3.4. Numerical results

We optimize (7.12) when $\Omega$ is the ring-shaped electroosmotic mixer depicted on Figure 12 [124, 130, 131]. In the perspective of maximizing the power of the electric field inside $\Omega$, the function $j: \mathbb{R}^3 \to \mathbb{R}$ featured in the objective function (7.14) is defined by $j(V) = -\gamma|V|^2$.

The initial mesh $\mathcal{T}^0$ of the domain $\Omega$ is composed of approximately 29,000 vertices and 135,000 tetrahedra. Throughout the optimization process, we impose that the top (resp. bottom) of the inflow and outflow channels
<table><caption>TABLE 1. Values of the parameters used in the (left) first, (center) second, and (right) third experiments of optimal design of the anode and cathode of an electroosmotic mixer considered in Section 7.3.</caption><thead><tr><th></th><th>l<sub>C</sub></th><th>l<sub>A</sub></th><th>m<sub>C</sub></th><th>m<sub>A</sub></th></tr></thead><tbody><tr><td>Experiment 1</td><td>0.0001</td><td>0.0001</td><td>0</td><td>0</td></tr><tr><td>Experiment 2</td><td>0.001</td><td>0.001</td><td>0</td><td>0</td></tr><tr><td>Experiment 3</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td></tr></tbody></table>

are part of the anode (resp. cathode). At each iteration of Algorithm 1, we alternately optimize either $\Gamma_C$ or $\Gamma_A$: we rely on the use of topological derivatives to attach new connected components to these regions every 10 iterations (i.e. $n_{top} = 10$ in Algorithm 2) until 100 iterations are reached, after which only geometric optimization steps are performed, based on the use of shape derivatives.

Three numerical experiments are conducted. In all cases, the minimum (resp. maximum) size of an element in the mesh is $h_{min} = 0.05$ (resp. $h_{max} = 0.5$), and the parameter $\varepsilon$ in the regularization process of shape derivatives equals $\varepsilon = 0.001$. Different values are used for the weights $l_C, l_A, m_C, m_A$ in Algorithm 2, as reported in Table 1.

Our first experiment deals with the maximization of the electric field inside $\Omega$ under a penalization of the areas of $\Gamma_C$ and $\Gamma_A$, without penalization of their contours. We perform 250 iterations of our optimization Algorithm 2 and a few snapshots of the optimization process are given in Figure 13, see Figure 14(e) for the convergence history. The method tends to assign the cathode and anode to opposite sides of the mixer $\Omega$. Initially, an anode is formed at the middle of the upper side and a corresponding cathode appears symmetrically on its lower side. The geometric optimization process then results in the expansion of the existing regions, which gradually attempt to cover the entire upper and lower sides of the mixer. Finally, the areas of both regions decrease because of the penalization of their values in the objective (7.12), revealing a pattern which is strongly reminiscent of the homogenization effect whereby shapes get close to optimal by developing very thin features, at the microscopic level, see e.g. [41, 132, 133]. The final mesh $\mathcal{T}$ of $\Omega$ contains approximately 5 million tetrahedrons, and the total simulation lasted around 3 hours.

In the second experiment, we use larger penalization parameters $l_A, l_C$ for the areas of $\Gamma_A$ and $\Gamma_C$. The results are presented in Figure 15. As in the previous experiment, both regions independently develop multiple branches, while occupying a lower area of the boundary $\partial\Omega$ of the mixer. The final mesh consists of 1.8 million tetrahedra. The total computational time is approximately 2 hours.

In our third experiment, we add a penalization on the lengths of the contours of $\Gamma_A$ and $\Gamma_C$, while keeping the same penalization parameters for their areas as in the second experiment. The results are reported on Figure 15 (b,d). As expected, the optimized regions have simpler shapes, containing fewer branches. The final mesh also consists of 1.8 million tetrahedra, for a total simulation time of about 2 hours.

## 7.4. Optimization of absorbent regions on the surface of a sound-hard obstacle for acoustic cloaking

In acoustics, scattering occurs when an incident wave hits an obstacle $\Omega$ or a discontinuity of the material parameters in the propagation medium. This interaction induces a scattered wave by deflection, absorption, or transmission in various directions, which is testament to the presence of the obstacle. In medical applications, one aims to exacerbate this effect in order to probe a patient's body, see e.g. [134] about tomography imaging methods used for medical diagnosis. On the contrary, scattering is highly undesirable in military applications, where one would rather try to hide the obstacle $\Omega$, which represents for instance a stealth submarine, or an aircraft.

In the latter type of situations, acoustic cloaking techniques have been developed to make an obstacle invisible to an incident wave, i.e. to attenuate or suppress the wave scattered by the obstacle. Among them, passive techniques consist in surrounding the obstacle with a shell made of a metamaterial with very specific properties, see e.g. [135–141], or [142] for a review. Recent contributions such as [143–145] have proposed density-based topology optimization strategies for the distribution of such a material around $\Omega$. In this section, we revisit
<figure><img src="image_16.png" /><figcaption>FIGURE 13. Snapshots of the optimization process of the distribution of the anode (in orange) and cathode (in blue) on the boundary of the electroosmotic mixer, in the first experiment of Section 7.3.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIGURE 14. (a,b) Optimized designs of the anode and cathode obtained in the first experiment of Section 7.3; (c,d) perspective views of the top and bottom parts of the mixer; (e) Convergence history in this first experiment.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>FIGURE 15. (a) Optimized design obtained in the second experiment of Section 7.3; (b) Optimized design obtained in the third experiment; (c) Convergence history for the second experiment; (d) Convergence history for the third experiment.</figcaption></figure>
the design of a cloaking mechanism for an obstacle $\Omega$ from a slightly different perspective: we aim to make $\Omega$ invisible to detection by optimizing the arrangement of constituent materials on its boundary rather than the properties of the surrounding region.

### 7.4.1. Presentation of the scattering problem and of the optimization problem

Our physical model is based on that proposed in [146, 147], see also Figure 16 for an illustration. Let $\Omega \subset \mathbb{R}^3$ be a smooth obstacle immersed in a fluid with specific volume $\gamma \in C^\infty(\mathbb{R}^3)$; this function is uniformly bounded away from 0 and $\infty$, as in (3.1). In the time-harmonic regime, we assume that an incident wave with wave number $k$ and complex-valued amplitude $u_{\text{in}}(x)$ passes through the medium, satisfying the Helmholtz equation (postulating a time dependency of the form $e^{-ikt}$):

$$-\operatorname{div}(\gamma(x)\nabla u_{\text{in}}(x)) - k^2 u_{\text{in}}(x) = 0, \quad x \in \mathbb{R}^3.$$

Denoting by $u_{\Gamma_R}(x)$ and $u_{\text{tot}}(x) = u_{\text{in}}(x) + u_{\Gamma_R}(x)$ the scattered and total amplitudes of the sound pressure field, respectively, it holds:

$$-\operatorname{div}(\gamma(x)\nabla u_{\text{tot}}(x)) - k^2 u_{\text{tot}}(x) = 0, \quad x \in \mathbb{R}^3 \setminus \bar{\Omega},$$

whence the following equation for the scattered wave:

$$-\operatorname{div}(\gamma(x)\nabla u_{\Gamma_R}(x)) - k^2 u_{\Gamma_R}(x) = 0, \quad x \in \mathbb{R}^3 \setminus \bar{\Omega}.$$

The latter is complemented with the Sommerfeld radiation condition at infinity, expressing that the scattered field is outgoing:

$$\lim_{r \to \infty} r \left( \frac{\partial u_{\Gamma_R}}{\partial r}(x) - iku_{\Gamma_R}(x) \right) = 0, \qquad (7.18)$$

where $r = |x|$ is the radial spherical coordinate.

In practice, it is necessary to replace the infinite, free space where the wave is propagating by a bounded computational box $D$, and to somehow impose suitable conditions near its boundary $\partial D$ which mimick the outgoing behavior (7.18). To achieve this, we rely on the simple method proposed in [148]: we consider an artificial boundary $\Gamma_E \subset D$ around the obstacle $\Omega$; an approximation for (7.18) is then obtained from a first-order approximation of the scattered field $u(x)$ “far away” from $\Omega$:

$$\left[ \gamma \frac{\partial u_{\Gamma_R}}{\partial n} \right] - \left( ik - \frac{1}{R} \right) u_{\Gamma_R} = 0 \quad \text{on } \Gamma_E, \qquad (7.19)$$

where

$$[\alpha](x) := \lim_{\substack{t \to 0 \\ t > 0}} \alpha(x + tn(x)) - \lim_{\substack{t \to 0 \\ t > 0}} \alpha(x - tn(x)), \quad x \in \Gamma_E,$$

denotes the jump of a quantity $\alpha$ which is discontinuous across $\Gamma_E$, and $R$ is the distance between $\Gamma_E$ and $\partial D$. Admittedly, more advanced numerical methods are available to impose the outgoing behavior (7.18) to $u_{\Gamma_R}$, such as the well-known Perfectly Matched Layer method [149].

The boundary of the obstacle $\partial \Omega$ is decomposed into two disjoint parts,

$$\partial \Omega = \overline{\Gamma_R} \cup \overline{\Gamma_N},$$
<figure><img src="image_19.png" /><figcaption>FIGURE 16. Illustration of the acoustic cloaking setting considered in Section 7.4.</figcaption></figure>

where:

* The region $\Gamma_N$ is covered with a “sound-hard” material, i.e. the incident wave $u_{in}$ is entirely reflected on this region, which translates as an inhomogeneous Neumann condition for the scattered field $u_{\Gamma_R}$:

$$ \gamma \frac{\partial u_{\Gamma_R}}{\partial n} = -\gamma \frac{\partial u_{in}}{\partial n}, \quad x \in \Gamma_N, $$

where the unit normal vector $n$ to $\partial\Omega$ is pointing outward the obstacle $\Omega$.

* The region $\Gamma_R$ is made of an absorbent material. The incoming wave is partially absorbed in there, and the following impedance boundary condition is satisfied:

$$ \gamma \frac{\partial u_{\Gamma_R}}{\partial n} + \frac{ik}{Z} u_{\Gamma_R} = -\gamma \frac{\partial u_{in}}{\partial n} - \frac{ik}{Z} u_{in}, \quad x \in \Gamma_R, $$

where $Z > 0$ is the acoustic impedance of the material.

Gathering the above information, the scattered field $u_{\Gamma_R} \in H^1(\Omega; \mathbb{C})$ is the unique complex-valued solution to the following boundary value problem:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(\gamma \nabla u_{\Gamma_R}) - k^2 u_{\Gamma_R} = 0 & \text{in } D \setminus \bar{\Omega} \\ \gamma \frac{\partial u_{\Gamma_R}}{\partial n} = -\gamma \frac{\partial u_{in}}{\partial n} & \text{on } \Gamma_N, \\ \gamma \frac{\partial u_{\Gamma_R}}{\partial n} + \frac{ik}{Z} u_{\Gamma_R} = -\gamma \frac{\partial u_{in}}{\partial n} - \frac{ik}{Z} u_{in} & \text{on } \Gamma_R, \\ \left[ \gamma \frac{\partial u_{\Gamma_R}}{\partial n} \right] - \left( ik - \frac{1}{R} \right) u_{\Gamma_R} = 0 & \text{on } \Gamma_E, \\ \gamma \frac{\partial u_{\Gamma_R}}{\partial n} = 0 & \text{on } \partial D. \end{array} \right. \quad (7.20) $$
The attached variational formulation reads:

$$
\text{Search for } u_{\Gamma_R} \in H^1(\Omega; \mathbb{C}) \text{ s.t. } \forall v \in H^1(\Omega; \mathbb{C}), \\
\int_{D \setminus \bar{\Omega}} \gamma \nabla u_{\Gamma_R} \cdot \bar{\nabla} v \, dx - k^2 \int_{D \setminus \bar{\Omega}} u_{\Gamma_R} \bar{v} \, dx - \left(ik - \frac{1}{R}\right) \int_{\Gamma_E} u_{\Gamma_R} \bar{v} \, ds - \frac{ik}{Z} \int_{\Gamma_R} u_{\Gamma_R} \bar{v} \, ds = \frac{ik}{Z} \int_{\Gamma_R} u_{\text{in}} \bar{v} \, ds + \int_{\partial \Omega} \gamma \frac{\partial u_{\text{in}}}{\partial n} \bar{v} \, ds,
$$

see Remark 5.1 about the well-posedness of this problem.

In this setting, we aim to optimize the repartition of the absorbent and sound-hard regions $\Gamma_R$, $\Gamma_N$ within $\partial\Omega$ so as to minimize the amplitude of the scattered wave $u_{\Gamma_R}$ while keeping the amount of used absorbent material reasonably low; we thus consider the following shape and topology optimization problem:

$$
\min_{\Gamma_R \subset \partial \Omega} J(\Gamma_R) + \ell \operatorname{Area}(\Gamma_R), \quad \text{where } J(\Gamma_R) := \frac{1}{2\operatorname{Vol}(D \setminus \bar{\Omega})} \int_{D \setminus \bar{\Omega}} |u_{\Gamma_R}|^2 dx, \qquad (7.21)
$$

and $\ell > 0$ is a penalization parameter.

7.4.2. Shape derivative of the objective function J(Γ<sub>R</sub>)

The solution $u_{\Gamma_R}$ to the present boundary value problem (7.20) does not show such a weakly singular behavior as the function $u_{\Gamma_C, \Gamma_A}$ considered in Section 7.3, where a transition between homogeneous Dirichlet and Neumann boundary conditions is at play. The shape derivative of the objective function $J(\Gamma_R)$ in (7.21) can be calculated by standard means, and we omit the details for brevity, see [121] about this point.

**Proposition 7.3.** The criterion $J(\Gamma_R)$ is shape differentiable at $\theta = 0$ and its shape derivative reads, for any tangential deformation $\theta$:

$$
J'(\Gamma_R)(\theta) = -\operatorname{Im} \left( \int_{\Sigma_R} \frac{k}{Z} \overline{(u_{\Gamma_R} + u_{\text{in}})} p_{\Gamma_R} \theta \cdot n_{\Sigma_R} \, d\ell \right),
$$

where the adjoint state $p_{\Gamma_R} \in H^1(\Omega; \mathbb{C})$ is the unique (complex-valued) solution to the following boundary value problem:

$$
\left\{
\begin{array}{ll}
-\mathrm{div}(\gamma \nabla p_{\Gamma_R}) - k^2 p_{\Gamma_R} = - \frac{u_{\Gamma_R}}{\mathrm{Vol}(D \setminus \bar{\Omega})} & \text{in } D \setminus \bar{\Omega} \\
\gamma \frac{\partial p_{\Gamma_R}}{\partial n} = 0 & \text{on } \Gamma_N, \\
\gamma \frac{\partial p_{\Gamma_R}}{\partial n} - \frac{ik}{Z} p_{\Gamma_R} = 0 & \text{on } \Gamma_R, \\
\left[ \gamma \frac{\partial p_{\Gamma_R}}{\partial n} \right] - \left(-ik - \frac{1}{R}\right) p_{\Gamma_R} = 0 & \text{on } \Gamma_E, \\
\gamma \frac{\partial p_{\Gamma_R}}{\partial n} = 0 & \text{on } \partial D.
\end{array}
\right.
\qquad (7.22)
$$

7.4.3. Topological derivative of the functional J(Γ<sub>R</sub>)

The setting considered here is slightly different from the model of Section 5, insofar as the boundary-value problem (7.20) is posed on the exterior of $\Omega$, and does not involve a homogeneous Dirichlet region. Nevertheless, a simple variation of the proof of Section 5.3 allows to calculate the topological derivative of $J(\Gamma_R)$.

**Proposition 7.4.** Let $x_0 \in \Gamma_N$; the perturbed criterion $J((\Gamma_R)_{x_0,\varepsilon})$ has the following asymptotic expansion:

$$
J((\Gamma_R)_{x_0, \varepsilon}) = J(\Gamma_R) + \pi \varepsilon^2 k \operatorname{Im} \left( \overline{(u_{\Gamma_R}(x_0) + u_{\text{in}}(x_0))} p_{\Gamma_R}(x_0) \right) + o(\varepsilon^2) \quad \text{if } d=3,
$$

where the adjoint state $p_{\Gamma_R} \in H^1(\Omega; \mathbb{C})$ is the complex-valued solution to (7.22).
<figure><img src="image_20.png" /><figcaption>FIGURE 17. (a) Mesh K of the initial design of the aircraft in Section 7.4, whose boundary is made only of sound-hard material; (b) Computational mesh T of the computational domain D, enclosing the fictitious surface ΓE bearing the approximate radiation conditions (7.19) and the mesh K of the aircraft.</figcaption></figure>

### 7.4.4. Setup and analysis of the numerical results

The obstacle Ω is an aircraft embedded in a large computational box D. The specific volume of the surrounding air is γ ≡ 1; the wave number is k = 2π/λ, where the wavelength λ equals 20. The considered incoming wave u_in is a plane wave, travelling vertically in the direction ξ = (0, 0, 1) ∈ R³, that is:

$$ \forall x \in \mathbb{R}^3, \quad u_{\text{in}}(x) = e^{ik\xi \cdot x}. $$

The domain D is equipped with a tetrahedral mesh T, which consistently encloses a submesh K of the obstacle, see Figure 17. At every iteration of the optimization process, the lengths of the edges of K range from h_max = λ/3 to h_min = λ/32, thus ensuring a sufficient resolution of the acoustic phenomenon at play [150]. For instance, the initial mesh T⁰ consists of about 334,000 vertices, with 1,970,000 tetrahedra, and the submesh K⁰ of the aircraft consists of 93,000 vertices and 464,000 tetrahedra, see Table 2 for a summary of the parameters of the experiment.

We solve the problem (7.21) with Algorithm 2, starting from a situation where the boundary ∂Ω is only made of sound-hard material: ΓR⁰ = ∅. The optimization process starts with 6 iterations during which new connected components are added to ΓR according to the information contained in the topological derivative of Proposition 7.4. These are followed by 4 iterations of motion via boundary variation, based on the result of Proposition 7.3. After the 10th iteration, we carry out one topological update of ΓR every 10 iterations,
<table><caption>TABLE 2. *Numerical parameters used in the optimization process of Section 7.4.*</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>ℓ</td><td>1e-7</td></tr><tr><td>λ</td><td>20</td></tr><tr><td>k</td><td>2π/λ</td></tr><tr><td>h<sub>max</sub></td><td>λ/3</td></tr><tr><td>h<sub>min</sub></td><td>λ/32</td></tr></tbody></table>

intertwined with updates of its geometry, up to the stage $n = 100$, after which only geometrical updates are considered.

A few snapshots of the evolution process are shown in Figure 18, and the convergence history is reported in Figure 19(d); the total computational time is about 19 hours. The first snapshot, associated to the iteration $n = 6$, displays the design after introducing 6 absorbent regions, which are located on the wings and the center of the aircraft. By iteration $n = 10$, the two central zones have vanished. As the optimization proceeds, the components added during the topological updates are systematically removed during the geometric optimization updates. The design gradually develops homogenization patterns, similar to those observed in Section 7.3; the absorbent region develops new connected components by “stretching and cutting off” some of its regions, rather than by the addition of new surface disks during the topological update steps.

The optimized design, attained at iteration $n = 142$, is illustrated in Figure 19(a,b,c): it suggests that the absorbent material should preferably be concentrated near the wings of the aircraft. This trend is quite intuitive, since these regions are perpendicular to the direction $\xi = (0,0,1)$ of the incident wave, which causes them to reflect the latter strongly. This fact is confirmed by examination of the magnitude of the scattered pressure field, which is displayed in Figure 20. Obviously, in the optimized design, the pressure field $u$ has been significantly reduced near the wings, making them less visible. The field has now the largest magnitude near the aircraft’s rear wings, its value being twice smaller compared to that in Figure 20(b). However, we believe that no absorbent material was added there due to the imposed penalization on the area of $\Gamma_R$, which urges the algorithm to concentrate the addition of absorbent material on the wings.

## 7.5. Optimization of the distribution of the structural supports of a water tank

In this section, we slip into structural mechanics, where we implement the techniques developed in Sections 3, 4 and 6 to optimize the supporting regions of a water tank. As illustrated in Figure 21, the structure features a cylindrical body and a domed upper part equipped with an access hatch. In the context of use, it is elevated on supporting legs which are anchored to a solid base plate, thus supporting the structure against gravity loads. Our goal is to optimize the distribution of these supporting regions on the boundary of the tank so as to minimize its total displacement.

### 7.5.1. Description of the optimization problem

The water tank is represented by a bounded domain $\Omega \subset \mathbb{R}^3$, where the boundary $\partial\Omega$ is decomposed into three disjoint regions, as:

$$ \partial\Omega = \overline{\Gamma_D} \cup \overline{\Gamma_F} \cup \overline{\Gamma}, $$

where:

* The region $\Gamma_D$, with boundary $\Sigma_D := \partial\Gamma_D$, is that where the water tank is fixed;

* The region $\Gamma_F$ is the upper part of the tank; it is traction-free and non optimizable (it cannot be riveted, as it serves as the lid of the tank);

* The remaining region $\Gamma$ is traction-free.
<figure><img src="image_21.png" /><figcaption>FIGURE 18. Snapshots of the optimization process of the distribution of absorbent (in pink) and “sound-hard” (in white) materials on the boundary of an aircraft (Sect. 7.4).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>FIGURE 18. (cont.) Snapshots of the optimization process of the distribution of “sound-soft” (in pink) and “sound-hard” (in white) materials on the boundary of an aircraft (Sect. 7.4).</figcaption></figure>

The water tank $\Omega$ is subject to body forces $f: \mathbb{R}^3 \to \mathbb{R}^3$ representing gravity. Its displacement $u_{\Gamma_D}$ in these circumstances is the unique solution in $H^1(\Omega)^3$ to the linearized elasticity system:

$$ \left\{ \begin{array}{ll} -\operatorname{div}(Ae(u_{\Gamma_D})) = f & \text{in } \Omega, \\ Ae(u_{\Gamma_D})n = 0 & \text{on } \Gamma \cup \Gamma_F, \\ u_{\Gamma_D} = 0 & \text{on } \Gamma_D. \end{array} \right. \qquad (7.23) $$
<figure><img src="image_23.png" /><figcaption>FIGURE 19. (a,b,c) Several views of the optimized distribution of absorbent and “sound-hard” materials on the boundary of the aircraft considered in Section 7.4; (d) History of the values of $J(\Gamma_R)$ during the computation.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIGURE 19. (cont.) (a,b,c) Several views of the optimized distribution of absorbent and “sound-hard” materials on the boundary of the aircraft considered in Section 7.4; (d) History of the values of $J(\Gamma_R)$ during the computation.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>FIGURE 20. Magnitude of the scattered pressure field, with (a,b) Only a “sound-hard” boundary; (c,d) The optimized repartition of absorbent and “sound-hard” materials on the boundary of the aircraft in the example of Section 7.4.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>FIGURE 21. (a) Sketch of the physical setting and (b) Tetrahedral mesh $\mathcal{T}$ of the initial guess in the optimal design example of supports on the surface of a water tank considered in Section 7.5.</figcaption></figure>
In this setting, we wish to minimize the mean displacement of the structure under the effect of gravity, while utilizing a reasonably low amount of supports, i.e. we consider the problem:

$$ \min_{\Gamma_D \subset \partial\Omega} J(\Gamma_D) + \ell \text{Area}(\Gamma_D) + m \text{Cont}(\Gamma_D), \quad \text{where } J(\Gamma_D) = \frac{1}{2\text{Vol}(\Omega)} \int_{\Omega} |u_{\Gamma_D}|^2, \qquad (7.24) $$

and $\ell$ and $m$ are fixed penalization parameters for the area and contour length of $\Gamma_D$, respectively.

### 7.5.2. Shape derivative of the functional $J(\Gamma_D)$

Like in the situation addressed in the previous Section 7.3, the weak singularity of the state function $u_{\Gamma_D}$ in (7.23) near the transition region $\Sigma_D$ between homogeneous Dirichlet and homogeneous Neumann boundary conditions makes it difficult to compute the shape derivative of the functional $J(\Gamma_D)$. We therefore replace the function $J(\Gamma_D)$ in (7.24) by that $J_{\varepsilon}(\Gamma_D)$ defined by:

$$ J_{\varepsilon}(\Gamma_D) = \frac{1}{2\text{Vol}(\Omega)} \int_{\Omega} |u_{\Gamma_D, \varepsilon}|^2 \, dx, \qquad (7.25) $$

where $u_{\Gamma_D, \varepsilon} \in H^1(\Omega)^3$ is the unique solution to the following approximate boundary value problem:

$$ \left\{ \begin{array}{ll} -\text{div}(Ae(u_{\Gamma_D, \varepsilon})) = f & \text{in } \Omega, \\ Ae(u_{\Gamma_D, \varepsilon})n + h_{\Gamma_D, \varepsilon}u_{\Gamma_D, \varepsilon} = 0 & \text{on } \partial\Omega. \end{array} \right. $$

Like in the previous investigations, the function $h_{\Gamma_D, \varepsilon} : \partial\Omega \to \mathbb{R}$ is defined by:

$$ \forall x \in \partial\Omega, \quad h_{\Gamma_D, \varepsilon}(x) = h\left(\frac{d_{\Gamma_D}^{\partial\Omega}(x)}{\varepsilon}\right), $$

where $d_{\Gamma_D}^{\partial\Omega}$ is the geodesic signed distance function to $\Gamma_D$ (see Sect. 3.2) and the transition profile $h \in C^{\infty}(\mathbb{R})$ satisfies (3.12).

The shape derivative of $J_{\varepsilon}(\Gamma_D)$ is supplied by the next proposition, whose proof is omitted for brevity.

**Proposition 7.5.** The functional $J_{\varepsilon}(\Gamma_D)$ in (7.25) is shape differentiable and its shape derivative reads, for an arbitrary tangential deformation $\theta$ (i.e. such that $\theta \cdot n = 0$):

$$ J'_{\varepsilon}(\Gamma_D)(\theta) = -\frac{1}{\varepsilon^2} \int_{\partial\Omega} h' \left( \frac{d_{\Gamma_D}^{\partial\Omega}(x)}{\varepsilon} \right) \theta(\pi_{\Sigma_D}(x)) \cdot n_{\Sigma_D}(\pi_{\Sigma_D}(x)) u_{\Gamma_D, \varepsilon}(x) \cdot p_{\Gamma_D, \varepsilon}(x) \, ds(x), $$

where the adjoint state $p_{\Gamma_D, \varepsilon} \in H^1(\Omega)^3$ is characterized by the following boundary value problem:

$$ \left\{ \begin{array}{ll} -\text{div}(Ae(p_{\Gamma_D, \varepsilon})) = -\frac{u_{\Gamma_D, \varepsilon}}{\text{Vol}(\Omega)} & \text{in } \Omega, \\ Ae(p_{\Gamma_D, \varepsilon})n + h_{\Gamma_D, \varepsilon}p_{\Gamma_D, \varepsilon} = 0 & \text{on } \partial\Omega. \end{array} \right. $$

### 7.5.3. Topological derivative of the functional $J(\Gamma_D)$

The sensitivity of the objective function $J(\Gamma_D)$ in (7.25) with respect to the addition to $\Gamma_D$ of a small surface disk $\omega_{x_0, \varepsilon}$ centered at $x_0 \in \Gamma$ is given by the following result, which directly follows from the study conducted in Section 6.4.
**Proposition 7.6.** Let $\Gamma_D$ be a region of the smooth boundary $\partial\Omega$ and let $x_0 \in \Gamma$ be given. Then, the perturbed criterion $J((\Gamma_D)_{x_0,\varepsilon})$, accounting for the addition of the surface disk $\omega_{x_0,\varepsilon} \subset \Gamma$ to $\Gamma_D$, has the following asymptotic expansion:

$$J((\Gamma_D)_{x_0,\varepsilon}) = J(\Gamma_D) + \frac{1}{|\log \varepsilon|} \frac{\pi \mu}{1 - \bar{\nu}} u_{\Gamma_D}(x_0) \cdot p_{\Gamma_D}(x_0) + o\left(\frac{1}{|\log \varepsilon|}\right) \text{ if } d = 2,$$

and

$$J((\Gamma_D)_{x_0,\varepsilon}) = J(\Gamma_D) + \varepsilon M u_{\Gamma_D}(x_0) \cdot p_{\Gamma_D}(x_0) + o(\varepsilon) \text{ if } d = 3.$$

Here, the polarization tensor $M$ is defined by (6.10) and the adjoint state $p_{\Gamma_D} \in H^1(\Omega)^3$ is characterized by the following boundary value problem:

$$\left\{ \begin{array}{ll} -\mathrm{div}(Ae(p_{\Gamma_D})) = -\frac{u_{\Gamma_D}}{\mathrm{Vol}(\Omega)} & \text{in } \Omega, \\ Ae(p_{\Gamma_D})n = 0 & \text{on } \Gamma \cup \Gamma_F, \\ p_{\Gamma_D} = 0 & \text{on } \Gamma_D. \end{array} \right.$$

### 7.5.4. Description of the numerical experiments and results

The initial configuration is supplied by a tetrahedral mesh $\mathcal{T}^0$ of the water tank $\Omega$, consisting of 17,000 vertices and 81,000 tetrahedra, which is displayed in Figure 21(b). The initial region $\Gamma_D^0$ is a disk with unit radius, located on the bottom of the water tank.

The parameters of the constituent elastic material of $\Omega$ are $\lambda = 0.5769$, $\mu = 0.3846$ (corresponding to a Young's modulus $E = 1$ and Poisson's ratio $\nu = 0.3$) and the gravity force $f$ equals $f = (0, 0, -0.1)$. The minimum (resp. maximum) size of an element in the mesh is $h_{\min} = 0.05$ (resp. $h_{\max} = 0.5$) and the regularization parameter $\varepsilon$ featured in the approximate functional (7.25) equals $1e-6$. The penalization parameter $\ell$ for the area of $\Gamma_D$ equals $1e-5$, and we conduct two experiments in this context, associated to different penalization parameters $m$ for the length of the contour $\Sigma_D$. Each resolution starts with 5 iterations, during which the topological derivative of $J(\Gamma_D)$ is used to add new small regions to $\Gamma_D$, followed by 45 iterations of geometric optimization, during which the newly created regions evolve by deformation of their boundaries. We then carry out one application of topological derivative every 10 iterations until reaching $n = 100$, at which point only geometric optimization updates are used.

In the first experiment, no penalization on the contour length of $\Gamma_D$ is imposed: $m = 0$. The results are depicted on Figures 22 and 23 and the total simulation takes approximately 3 hours. The algorithm clearly prioritizes the “ribs” around the flat indentations for the placement of supporting regions. During the initial iterations, based on the topological derivative of $J(\Gamma_D)$, supports are immediately positioned near these ribs, and the subsequent geometric updates focus on covering most of these areas. The bottom region of the tank remains essentially unchanged, highlighting its essential role in minimizing the displacement of the water tank.

Our second experiment features a penalization of the contour $\mathrm{Cont}(\Gamma_D)$ of the support region with the parameter $m = 1e-5$. The results are depicted on Figures 24 and 25; the complete simulation requires about 2 hours. Similar trends are observed, while the optimized design is more regular. In this final design, all the supporting areas are concentrated near the ribs of the water tank or its bottom, in good agreement with observations made in previous contributions to the analysis of this problem [151, 152]. This demonstrates that corrugated structures exhibit high resistance to loads applied perpendicular to the direction of the ribs. The smooth minimization of the objective by the algorithm further reinforces the validity of our results, confirming the effectiveness and reliability of our approach.
<figure><img src="image_27.png" /><figcaption>FIGURE 22. A few intermediate shapes produced during the first experiment of the optimization of the support regions (in pink) of the water tank considered in Section 7.5; the blue region is not subject to optimization.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>FIGURE 23. (a) Final design in the first experiment of the optimization of the supporting region of the water tank considered in Section 7.5; (b) Associated convergence history.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>FIGURE 24. *A few intermediate shapes produced during the second experiment of the optimization of the support regions (in pink) of the water tank considered in Section 7.5; the blue region is not subject to optimization.*</figcaption></figure>

## 7.6. Optimization of a clamping-locator fixture system

So-called “clamping-locator” fixture systems are crucial ingredients in various engineering processes, where they allow to precisely position and hold objects during operations. For instance, such devices are widespread components of biological experiments, where delicate samples are observed by microscopes and must be handled without damage [153], or in medical operations, during which surgical tools need to be positioned tightly in place. They are also ubiquitous in industrial manufacturing processes, such as CNC milling [154] or other workflows in robotics.

As their name suggests, clamping-locator systems are made of two parts: locators hamper the motion of the target component in one or several directions, while clamps apply a very strong pressing force against locators, see Figure 26. Clamping-locator systems have numerous advantages in machining; however, the application of an excessive clamping force can deform delicate workpieces, leading to higher scrap rates and higher material costs. Therefore, it is essential to optimize the geometry and placement of the clamping and locator regions to minimize the deformation of the mechanical part incurred by the system. In this section, we address this issue by using our boundary optimization technology.
<figure><img src="image_30.png" /><figcaption>FIGURE 25. (a) Final design in the second experiment of the optimization of the support regions of the water tank considered in Section 7.5; (b) Bottom view of the optimized shape; (c) Associated convergence history.</figcaption></figure>

### 7.6.1. The optimization problem

The considered mechanical structure is represented by a bounded domain $\Omega \subset \mathbb{R}^3$, whose boundary $\partial\Omega$ is made of five disjoint pieces:

$$ \partial\Omega = \overline{\Gamma_D} \cup \overline{\Gamma_N} \cup \overline{\Gamma_T} \cup \overline{\Gamma_F} \cup \overline{\Gamma}. \qquad (7.26) $$

In this decomposition:

* The region $\Gamma_D$ is that where locators are operating, i.e. the displacement of the piece is prevented.

* A strong force $f$, where $f: \Gamma_N \to \mathbb{R}$ is applied in the normal direction on the clamping region $\Gamma_N$.

* The tool applies a force $g: \Gamma_T \to \mathbb{R}^3$ on the region $\Gamma_T$, which is fixed and not subject to optimization.

* No efforts are applied on the region $\Gamma_F$, which is not subject to optimization.

* The remaining region $\Gamma$ is also traction-free, but it is subject to optimization.

We also set $\Sigma_D := \partial\Gamma_D$, $\Sigma_N := \partial\Gamma_N$, and we assume that the sets $\Gamma_D$ and $\Gamma_N$ always stay well-separated by $\Gamma$, i.e. $\overline{\Gamma_D} \cap \overline{\Gamma_N} = \emptyset$, see Figure 26 for an illustration.
<figure><img src="image_31.png" /><figcaption>FIGURE 26. Schematics of a clamping-locator system as considered in Section 7.6.</figcaption></figure>

Neglecting body forces, the displacement of the structure in this situation is the unique solution $u \in H^1(\Omega)^3$ to the following version of the linearized elasticity system:

$$ egin{cases} -\operatorname{div}(Ae(u_{\Gamma_D, \Gamma_N})) = 0 & \text{in } \Omega, \\ Ae(u_{\Gamma_D, \Gamma_N})n = fn & \text{on } \Gamma_N, \\ Ae(u_{\Gamma_D, \Gamma_N})n = g & \text{on } \Gamma_T, \\ Ae(u_{\Gamma_D, \Gamma_N})n = 0 & \text{on } \Gamma_F \cup \Gamma, \\ u_{\Gamma_D, \Gamma_N} = 0 & \text{on } \Gamma_D, \end{cases} $$

for given functions $f: \mathbb{R}^3 \to \mathbb{R}$ and $g: \mathbb{R}^3 \to \mathbb{R}^3$. We seek to minimize the total displacement of the structure, while utilizing a reasonably low amount of clamps and locators. This problem is formulated as:

$$ \min_{\Gamma_D, \Gamma_N \subset \partial \Omega} J(\Gamma_D, \Gamma_N) + \ell_N \operatorname{Area}(\Gamma_N) + \ell_D \operatorname{Area}(\Gamma_D), \quad \text{where } J(\Gamma_D, \Gamma_N) := \frac{1}{2\operatorname{Vol}(\Omega)} \int_{\Omega} |u_{\Gamma_D, \Gamma_N}|^2 dx \quad (7.27) $$

where $\ell_D, \ell_N > 0$ are penalization parameters for the areas of the clamping and locator regions, respectively.

### 7.6.2. Approximate version of the functional $J(\Gamma_D, \Gamma_N)$ and its shape derivative

We follow once again the procedure of Section 3.3 for deriving a smooth approximation of the shape functional $J(\Gamma_D, \Gamma_N)$. More precisely, we trade the exact shape and topology optimization problem (7.27) for an approximate counterpart, where the transition between homogeneous Dirichlet and Neumann boundary conditions around $\Sigma_D$ is smoothed by a Robin boundary condition; the transition between homogeneous and inhomogeneous Neumann conditions near $\Sigma_N$ is unaltered, as it does not pose any particular difficulty, see Section 3.5.