Shape optimization for interface identification in nonlocal models

Shape optimization methods have been proven useful for identifying interfaces in models governed by partial differential equations. Here we consider a class of shape optimization problems constrained by nonlocal equations which involve interface-dependent kernels. We derive a novel shape derivative associated to the nonlocal system model and solve the problem by established numerical techniques.


Introduction
Many physical relations and data-based coherences cannot satisfactorily be described by classical differential equations.Often they inherently possess some features, which are not purely local.In this regard, mathematical models which are governed by nonlocal operators enrich our modeling spectrum and present useful alternative as well as supplemental approaches.That is why they appear in a large variety of applications including among others, anomalous or fractional diffusion [9,10,18], peridynamics [24,51,60,26], image processing [29,35,39], cardiology [13], machine learning [41], as well as finance and jump processes [34,6,5,55,25].Nonlocal operators are integral operators allowing for interactions between two distinct points in space.The nonlocal models investigated in this paper involve kernels that are not necessarily symmetric and which are assumed to have a finite range of nonlocal interactions; see, e.g, [22,56,23,25] and the references therein.
Not only the problem itself but also various optimization problems involving nonlocal models of this type are treated in literature.For example matching-type problems are treated in [19,17,20] to identify system parameters such as the forcing term or a scalar diffusion parameter.The control variable is typically modeled to be an element of a suitable function space.Moreover, nonlocal interface problems have become popular in recent years [16,40,12].However, shape optimization techniques applied to nonlocal models can hardly be found in literature.For instance, the articles [8,52,38] deal with minimizing (functions of) eigenvalues of the fractional Laplacian with respect to the domain of interest.Also, in [14,7] the energy functional related to fractional equations is minimized.In [11] a functional involving a more general kernel is considered.All of the aforementioned papers are of theoretical nature only.To the best of our knowledge, shape optimization problems involving nonlocal constraint equations with truncated kernels and numerical methods for solving such problems cannot yet be found in literature.
Instead, shape optimization problems which are constrained by partial differential equations appear in many fields of application [43,31,49,50] and particularly for inverse problems where the parameter to be estimated, e.g., the diffusivity in a heat equation model, is assumed to be defined piecewise on certain subdomains.Given a rough picture of the configuration, shape optimization techniques can be successfully applied to identify the detailed shape of these subdomains [47,45,46,58].
In this paper we transfer the problem of parameter identification into a nonlocal regime.Here, the parameter of interest is given by the kernel which describes the nonlocal model.We assume that this kernel is defined piecewise with respect to a given partition {Ω i } i of the domain of interest Ω.Thereby, the state of such a nonlocal model depends on the interfaces between the respective subdomains Ω i .Under the assumption that we know the rough setting but are lacking in details, we can apply the techniques developed in the aforementioned shape optimization papers to identify these interfaces from a given measured state.
For this purpose we formulate a shape optimization problem which is constrained by an interface-dependent nonlocal convection-diffusion model.Here, we do not aim at investigating conceptual improvements of existing shape optimization algorithms.On the contrary, we want to study the applicability of established methods for problems of this type.Thus this paper can be regarded as a feasibility study where we set a focus on the numerical implementation.
The realization of this plan basically requires two ingredients both of which are worked out here.First, we define a reasonable interface-dependent nonlocal model and provide a finite element code which discretizes a variational formulation thereof.Second, we need to derive the shape derivative of the corresponding nonlocal bilinear form which is then implemented into an overall shape optimization algorithm.
This leads to the following organization of the present paper.In Section 2 we formulate the shape optimization problem including an interface-dependent nonlocal model.Once established, we briefly recall basic concepts from the shape optimization regime in Section 3. Then Section 4 is devoted to the task of computing the shape derivative of the nonlocal bilinear form and the reduced objective functional.Finally we present numerical illustrations in Section 5 which corroborate theoretical findings.

Problem formulation
The system model to be considered is the homogeneous steady-state nonlocal Dirichlet problem with volume constraints, given by posed on a bounded domain Ω ⊂ R d ; see, e.g, [22,4,23,25,56] and the references therein.Here, we assume that this domain is partitioned into a simply connected interior subdomain Ω 1 ⊂ Ω with boundary Γ := ∂Ω 1 and a domain Ω 2 := Ω\Ω 1 .Thus we have Ω = Ω(Γ) = Ω 1 ∪Γ ∪Ω 2 , where ∪ denotes the disjoint union.In the following, the boundary Γ of the interior domain Ω 1 is called the interface and is assumed to be an element of an appropriate shape space; see also Section 3 for a related discussion.The governing operator L Γ is an interface-dependent, nonlocal convection-diffusion operator of the form which is determined by a nonnegative, interface-dependent (interaction) kernel γ Γ : R d ×R d → R.
The second equation in ( 1) is called Dirichlet volume constraint.It specifies the values of u on the interaction domain which consists of all points in the complement of Ω that interact with points in Ω.For ease of exposition, we set u = 0 on Ω I , but generally we can use the constrained u = g on Ω I , if g satisfies some appropriate assumptions.Furthermore, we assume that the kernel depends on the interface in the following way where χ Ωi×Ωj denotes the indicator of the set Ω i × Ω j .For instance, in [48] the authors refer to γ ij and γ iI as inter-and intra-material coefficients.Notice that we do not need kernels γ Ii , since u = 0 on Ω I .Furthermore, throughout this work we consider truncated interaction kernels, which can be written as for appropriate positive functions φ ij : R d × R d → R and φ iI : R d × R d → R, which we refer to as kernel functions.We assume for i ∈ {1, 2} that there exist two radii 0 denotes the Euclidean ball of radius ε k i .In this paper we differentiate between square integrable kernels and singular symmetric kernels.For square integrable kernels we require We do not assume that (3) is symmetric for this type of kernels.In the case of singular symmetric kernels we require the existence of constants 0 for x ∈ Ω and y ∈ S 1 (x) ∪ S 2 (x).Also, since the singular kernel is required to be symmetric, the condition γ(x, y) = γ(y, x), and, respectively, φ 12 (x, y) = φ 21 (y, x), φ ii (x, y) = φ ii (y, x) has to hold.Because we do not need to define γ Ii , as described above, there is no further symmetry condition for γ iI required.
Example 2.1.One example of such a singular symmetric kernel is given by where 0 < ε < ∞ and the functions σ ij , σ iI : R d × R d → R are bounded from below and above by some positive constants γ * and γ * .Additionally, the σ ii are assumed to be symmetric on Ω × Ω and σ 12 (x, y) = σ 21 (y, x) holds for x, y ∈ Ω ∪ Ω I .
For the forcing term f Γ in (1) we assume a dependency on the interface in the following way where we assume that f i ∈ H 1 (Ω), i = 1, 2, because we need that f is weakly differentiable in Section 4. Next, we introduce a variational formulation of problem (1).For this purpose we define the corresponding forms and for some functions u, v : Ω ∪ Ω I → R. By inserting the definitions of the nonlocal operator (2) with the kernel given in (4) and the definition of the forcing term (5), we obtain the nonlocal bilinear form and the linear functional In order to derive the second bilinear form (8) we used Fubini's theorem.We make use of both representations (7) and (8) of the nonlocal bilinear form in the proofs of Section 4 .For singular symmetric kernels we also use another equivalent representation of the nonlocal bilinear form given by where we again used Fubini's theorem and applied that u, v = 0 on Ω I .Next, we use the nonlocal bilinear form to define a seminorm With this seminorm, we further define the energy spaces We now formulate the variational formulation corresponding to problem (1) as follows For square integrable kernels one can show under appropriate assumptions on the kernel the equivalence between V (Ω ∪ Ω I ) and L 2 (Ω ∪ Ω I ) and, respectively, the equivalence of V c (Ω ∪ Ω I ) and L 2 c (Ω ∪ Ω I ), see related results in [23,27,57].Moreover, the well-posedness of problem (10) for symmetric (square integrable) kernels is proven in [23] and in [27] the well-posedness for some nonsymmetric cases is also covered (again under certain conditions on the kernel and the forcing term f ).For the singular symmetric kernels the well-posedness of problem (10), the equivalence between V (Ω ∪ Ω I ) and the fractional Sobolev space H s (Ω ∪ Ω I ) and between V c (Ω ∪ Ω I ) and H s c (Ω ∪ Ω I ) is shown in [23].Finally, let us suppose we are given certain measurements ū : Ω → R on the domain Ω, which we assume to follow the nonlocal model (10) with the interface-dependent kernel γ Γ and the forcing term f Γ defined in (3) and (5), respectively.In oder to formulate the shape derivative in Chapter 4 we need ū ∈ H 1 (Ω).Then, given the data ū we aim at identifying the interface Γ for which the corresponding nonlocal solution u(Γ) is the "best approximation" to these measurements.Mathematically spoken, we formulate an optimal control problem with a tracking-type objective functional where the interface Γ, modeled as a shape, represents the control variable.We now assume Ω := (0, 1) 2 and introduce the following nonlocally constrained shape optimization problem min The objective functional is given by The first term j(u, Γ) is a standard L 2 tracking-type functional "projecting" the data on the set of reachable solutions, whereas the second term j reg (Γ) is known as the perimeter regularization, which is commonly used in the related literature to overcome possible ill-posedness of optimization problems [3].

Basic concepts in shape optimization
For solving the constrained shape optimization problem (11) we want to use the same shape optimization algorithms as they are developed in [47,45,44] for problem classes that are comparable in structure.Thus, in this section we briefly introduce the basic concepts and ideas of the therein applied shape formalism.For a rigorous introduction to shape spaces, shape derivatives and shape calculus in general, we refer to the monographs [15,53,58].

Notations and definitions
Based on our perception of the interface, we now refer to the image of a simple closed and smooth curve as a shape, i.e., the spaces of interest are subsets of By the Jordan curve theorem [30] such a shape Γ ∈ A divides the plane into two (simply) connected components with common boundary Γ.One of them is the bounded interior, which in our situation can then be identified with Ω 1 .Functionals J : A → R which assign a real number to a shape are called shape functionals.Since this paper deals with minimizing such shape functionals, i.e., with so-called shape optimization problems, we need to introduce the notion of an appropriate shape derivative.To this end we consider a family of mappings F t : Ω → R d with F 0 = id, where t ∈ [0, T ] and T > 0, which transform a shape Γ into a family of perturbed shapes Here the family of mappings {F t } t∈[0,T ] is described by the perturbation of identity, which for a smooth vector field We note that for sufficiently small t ∈ [0, T ] the function F t is injective, and thus Γ t ∈ A. Then the Eulerian or directional derivative of a shape functional J at a shape Γ in direction of a vector field At this point, let us also define the material derivative of a family of functions For functions v, which do not explicitly depend on the shape, i.e., v t = v for all t ∈ [0, T ], we derive For more details on shape optimization we refer to the literature, e.g., [37].

Optimization approach: Averaged adjoint method
Let us assume, that for each admissible shape Γ, there exists a unique solution u(Γ) of the constraint equation, i.e., u(Γ) Then we can consider the reduced problem In order to employ derivative based minimization algorithms, we need to derive the shape derivative of the reduced objective functional J red .By formally applying the chain rule, we obtain where D u J and D Γ J denote the partial derivatives of the objective J with respect to the state variable u and the control Γ, respectively.In applications we typically do not have an explicit formula for the control-to-state mapping u(Γ), so that we cannot analytically quantify the sensitivity of the unique solution u(Γ) with respect to the interface Γ. Thus, a formula for the shape derivative D Γ u(Γ)[V] is unattainable.One possible approach to access this derivative is the averaged adjoint method (AAM) of [33,54], which is a Lagrangian method, where the so-called Langrangian functional is defined as The basic idea behind Lagrangian methods is the aspect, that we can express the reduced functional as Now let Γ be fixed and denote by Γ t := F t (Γ) and Ω t i := F t (Ω i ) the deformed interior boundary and respectively the deformed domains.Furthermore let Ω t := F t (Ω) be the domain which corresponds to the interior boundary Γ t .Then we consider the reduced objective functional regarding Γ t , i.e., where u(Γ t ) ∈ V c (Ω t ∪ Ω I ).If we now try to differentiate L with respect to t in order to derive the shape derivative, we would have to compute the derivative for u(Γ t ) , and therefore the space Moreover we define 3.2 Optimization approach: Averaged adjoint method Then we can reformulate (15) as where , which is one prerequisite of the AAM.Then, in order to use the AAM to compute the shape derivative, the following additional assumptions have to be met.
• Assumption (H0): • For every t ∈ [0, T ] there exists a unique solution v t ∈ L 2 (Ω), such that v t solves the average adjoint equation • Assumption (H1): Assume that the following equation holds Then the next theorem yields a practical formula for deriving the shape derivative: ).Let the assumptions (H0) and (H1) be satisfied and suppose there exists a unique solution v t to the average adjoint equation (17).Then for v ∈ V c (Ω ∪ Ω I ) we obtain Proof.See proof of [33, Theorem 3.1].
For t = 0 the average adjoint equation ( 17) can be written as Here we also call (19) adjoint equation and the solution v 0 is referred to as the adjoint solution.
Moreover the nonlocal problem ( 10) is also called state equation and the solution u 0 is named state solution.

Optimization algorithm
Let us assume for a moment that we have an explicit formula for the shape derivative of the reduced objective functional.We now briefly recall the techniques developed in [47] and describe how to exploit this derivative for implementing gradient based optimization methods or even Quasi-Newton methods, such as L-BFGS, to solve the constrained shape optimization problem (11).
In order to identify gradients we need to require the notion of an inner product, or more generally a Riemannian metric.Unfortunately, shape spaces typically do not admit the structure of a linear space.However, in particular situations it is possible to define appropriate quotient spaces, which can be equipped with a Riemannian structure.For instance consider the set A introduced in ( 12).Since we are only interested in the image of the defining embedding, a re-parametrization thereof does not lead to a different shape.Consequently, two curves that are equal modulo (diffeomorphic) re-parametrizations define the same shape.This conception naturally leads to the quotient space Emb(S 1 , R d )/ Diff(S 1 , S 1 ), which can be considered an infinite-dimensional Riemannian manifold [36].This example already intimates the difficulty of translating abstract shape derivatives into discrete optimization methods; see, e.g., the thesis [59] on this topic.A detailed discussion of these issues is not the intention of this work and we now outline Algorithm 1.
The basic idea can be intuitively explained in the following way.Starting with an initial guess Γ 0 , we aim to iterate in a steepest-descent fashion over interfaces Γ k until we reach a "stationary point" of the reduced objective functional J red .The interface Γ k is encoded in the finite element mesh and transformations thereof are realized by adding vector fields U : Ω → R d (which can be interpreted as tangent vectors at a fixed interface) to the finite element nodes which we denote by Ω k .Thus, the essential part is to update the finite element mesh after each iteration by adding an appropriate transformation vector field.For this purpose, we use the solution U(Γ) : Ω(Γ) → R d of the so-called deformation equation The right-hand side of this equation is given by the shape derivative of the reduced objective functional (18) and the left-hand side denotes an inner product on the vector field space H 1 0 (Ω, R 2 ).In the view of the manifold interpretation, we can consider a Γ as inner product on the tangent space at Γ, so that U(Γ) is interpretable as the gradient of the shape functional J red at Γ.The solution U(Γ) : Ω → R 2 of ( 20) is then added to the coordinates Ω k of the finite element nodes.A common choice for a Γ is the bilinear form associated to the linear elasticity equation given by , where and are the strain and stress tensors, respectively.Deformation vector fields V which do not change the interface do not have an impact on the reduced objective functional, so that Therefore, the right-hand side D Γ J red (Γ)[V] is only assembled for test vector fields whose support intersects with the interface Γ and set to zero for all other basis vector fields.This prevents wrong mesh deformations resulting from discretization errors as outlined and illustrated in [45].Furthermore, λ and µ in (21) denote the Lamé parameters which do not need to have a physical meaning here.It is more important to understand their effect on the mesh deformation.They enable us to control the stiffness of the material and thus can be interpreted as some sort of step size.In [44], it is observed that locally varying Lamé parameters have a stabilizing effect on the mesh.A good strategy is to choose λ = 0 and µ as solution of the following Laplace equation Therefore µ min , µ max ∈ R influence the step size of the optimization algorithm.A small step is achieved by the choice of a large µ max .Note that a Γ then depends on the interface Γ through the parameter µ = µ(Γ) : Ω(Γ) → R.
Algorithm 1: Shape optimization algorithm Interpolate ū onto the current finite element mesh Ω k 4 Assemble A Γ and solve state (10) and adjoint equation ( 19) Compute the mesh deformation 7 Assemble shape derivative Update mesh 24 How to perform the limited memory L-BFGS update in Line 14 of Algorithm 1 within the shape formalism is investigated in [46,Section 4].Here, we only mention that the therein examined vector transport is approximated with the identity operator, so that we finally treat the gradients U k : Ω k → R d as vectors in R d|Ω k | and implement the standard L-BFGS update [44,Section 5].

Shape derivative of the reduced objective functional
In Section 3 we have depicted the optimization methodology, that we follow in this work to numerically solve the constrained shape optimization problem (11).In order to proof the requirements of the AAM, we need some additional assumptions.
• Additionally assume that there exists a constant 0 < C 0 < ∞, such that where A is defined as in (16).

Assumption (P1):
• For singular kernels: • For square integrable kernels: Let the kernel functions satisfy the requirements Singular kernels already satisfy assumption (P0) since the first condition is fulfilled by the theory of [23,57] and the second requirement is shown in the following Lemma: Lemma 4.1.In the case of a singular kernel, there exists a constant 0 < C 0 < ∞, so that Proof.Let ε := min{ε Since T is chosen small enough, [0, T ] × Ω is a compact set and ξ t is continuous on [0, T ] × Ω, there exists ξ * > 0, s.t.ξ t (x) ≥ ξ * for every t ∈ [0, T ] and x ∈ Ω.Therefore, by using that F t (Ω) = Ω, we derive In the following we proof that assumption (P1) also holds for a standard example of a singular symmetric kernel.
where we used that V ∈ C k 0 (Ω) is Lipschitz continuous for some Lipschitz constant L > 0 and that there exists a We will see in the proof of the following Lemma 4.3, that in our case the average adjoint equation ( 17) is equivalent to equation (23).Now we can show, that the additional requirements of AAM are satisfied by problem (11): Lemma 4.3.Let G be defined as in (16) and let the assumptions (P0) and (P1) be fulfilled.Then the assumptions (H0) and (H1) are satisfied and for every t ∈ [0, T ] there exists a solution v t ∈ V c (Ω ∪ Ω I ) that solves the average adjoint equation (17).
Proof.Because of the length of the proof, we move it to Appendix A.
The missing piece to implement the respective algorithmic realization presented in Subsection 3.3 is the shape derivative of the reduced objective functional, which is used in Line 8 of Algorithm 1 and given by As a first step, we formulate the shape derivative of the objective functional J and the linear functional F , which can also be found in the standard literature.

Theorem 4.4 (Shape derivative of the reduced objective functional). Let the assumptions (P0
) and (P1) be satisfied.Further let Γ be a shape with corresponding state variable u 0 and adjoint variable v 0 .Then, for a vector field Proof.In order to proof this theorem, we just have to compute the shape derivative of the objective function J(u 0 , Γ) and of the linear functional F Γ (v 0 ).Therefore, let ξ t (x) := det DF t (x).Then, we have ξ 0 (x) = det DF 0 (x) = det(I) = 1 and d dt t=0 + ξ t = div V(see e.g.[42]), such that the shape derivative of the right-hand side F Γ can be derived as follows Moreover, the shape derivative of the objective functional can be written as Here the shape derivative of the regularization term is an immediate consequence of [58, Theorem 4.13] and is given by where n denotes the outer normal of Ω 1 .Additionally, we obtain for the shape derivative of the tracking-type functional Putting the above terms into equation ( 24) yields the formula of Theorem 4.4.
The last step to derive the shape derivative of the reduced objective functional (24) is to compute the shape derivative of the nonlocal bilinear form A Γ , which is shown in the next Lemma.

Lemma 4.5 (Shape derivative of the nonlocal bilinear form). Let the assumptions (P0)
and (P1) be satisfied.Further let Γ be a shape with corresponding state variable u 0 and adjoint variable v 0 .Then for a vector field V ∈ C k 0 (Ω, R d ) we find for a square integrable kernel γ that and for a singular kernel γ that Proof.Define ξ t (x) := det DF t (x) and γ t ij (x, y) := γ ij (F t (x), F t (y)).Case 1: Square integrable kernels Then, we can write by using representation (8)

of the nonlocal bilinear form
So we derive the shape derivative of the nonlocal bilinear form For the second equation, the following computations are used, which can be obtained by applying Fubini's theorem and by swapping x and y we obtain Case 2: Singular kernels Analogously we get for the singular symmetric kernel 1 2 We can now derive the shape derivative of the reduced objective functional.If we formally set u 0 = ū and ν = 0, we can conclude from ( 19) that v 0 = 0 and therefore So, if there is a shape Γ, such that u(Γ) = u 0 = ū, then Γ is a stationary point of the reduced objective functional (14).

Numerical experiments
In this section, we want to put the above derived formula (25) for the shape derivative of the reduced objective functional into numerical practice.In the following numerical examples we test one singular symmetric and one nonsymmetric square integrable kernel.Specifically, where with with scaling constants d δ := 2−2s πδ 2−2s and where with scaling constants c δ := 1 δ 4 .We truncate all kernel functions by • 2 -balls of radius δ = 0.1 so that Ω ∪ Ω I ⊂ [−δ, 1 + δ] 2 .As a right-hand side we choose a piecewise constant function i.e., f 1 = 100 and f 2 = −10.We note that the nonsymmetric kernel γ nonsym satisfies the conditions for the class of integrable kernels considered in [57], such that the corresponding nonlocal problem is well-posed and also assumption (P1) can easily be verified in this case.The symmetric kernel γ sym is a special case of Example 2.1 and therefore the assumptions (P0) and (P1) are met.The well-posedness of the nonlocal problem regarding the singular kernel is shown in [23].As a perimeter regularization we choose ν = 0.002 and, since we only utilize V with supp(V) ∩ Γ k = ∅, we additionally assume that the nonlocal boundary has no direct influence on the shape derivative of the nonlocal bilinear form D Γ A Γ , such that for all V ∈ C k 0 (Ω, R d ) with supp(V) ∩ Γ k = ∅ we have for the square integrable kernel and for the singular symmetric kernel In order to solve problem (11), we apply a finite element method, where we employ continuous piecewise linear basis functions on triangular grids for the discretization of the nonlocal constraint equation.In particular we use the free meshing software Gmsh [28] to construct the meshes and the Python package nlfem [32] to assemble the stiffness matrices of the nonlocal state and adjoint equation.Moreover, to compute the load vector of the state and adjoint equation and the shape derivatives D Γ J and D Γ F Γ , we employ the open-source finite element software FEniCS [2,1].For a detailed discussion on the assembly of the nonlocal stiffness matrix we refer to [21,32].Here we solely emphasize how to implement a subdomain-dependent kernel of type (3).
During the mesh generation each triangle is labeled according to its subdomain affiliation.Thus, whenever we integrate over a pair of two triangles, we can read out the labels (i, j) and choose the corresponding kernel γ ij .The data ū is generated as solution u(Γ) of the constraint equation associated to a target shape Γ.Thus the data is represented as a linear combination of basis functions from the finite element basis.For the interpolation task in Line 3 of Algorithm 1 we solely need to translate between (non-matching) finite element grids by using the project function of FEniCS.In all examples below the target shape Γ is chosen to be a circle of radius 0.25 centered at (0.5, 0.5).
We now present two different non-trivial examples which differ in the choice of the initial guess Γ 0 .They are presented and described in the Since the start shapes are smaller than the target shape, the shape needs to expand in the first few iterations.Thereby the nodes of the mesh are pushed towards the boundary, so that the mesh quality decreases and the algorithm stagnates, because nodes are prohibited to be pushed outside of Ω.Therefore, we apply a re-meshing technique, where we re-mesh after the fifth and tenth iteration.In order to re-mesh, we save the points of our current shape as a spline in a dummy .geofile, that also contains the information of the nonlocal boundary, and then compute a new mesh with Gmsh.In this new mesh the distance between the nodes and the boundary is sufficiently large enough to attain a better improvement regarding the objective function value by the new mesh deformations.
It is important to mention that computation times and the performance of Algorithm 1 in general are very sensitive to the choice of parameters and may strongly vary, which is why reporting exact computation times is not very meaningful at this stage.Particularly delicate Moreover, especially in the case of system parameters with high interface-sensitivity in combination with an inconveniently small µ max , mesh deformations may be large in the early phase of the algorithm.Thus, such mesh deformations Ũk of high magnitude lead to destroyed meshes so that an evaluation of the reduced objective functional J red ((id + α Ũk )(Ω k )), which requires the assembly of the nonlocal stiffness matrix, becomes a pointless computation.In order to avoid such computations we first perform a line search depending on one simple mesh quality criterion.More precisely, we downscale the step size, i.e., α = τ α, until all finite element nodes of the resulting mesh (id + α Ũk )(Ω k ) are a subset of Ω.After that, we continue with the backtracking line search in Line 19 of Algorithm 1.

Concluding remarks and future work
We have conducted a numerical investigation of shape optimization problems which are constrained by nonlocal system models.We have proven through numerical experiments the applicability of established shape optimization techniques for which the shape derivative of the nonlocal bilinear form represents the novel ingredient.All in all, this work is only a first step along the exploration of the interesting field of nonlocally constrained shape optimization problems. that We can compute where Furthermore the second criterion of (H0) is also satisfied: As we saw above in the proof of assumption (H0), the left hand side of the average adjoint equation can be written as follows Then, (17) can be reformulated as Since the right hand side of ( 27) is a linear and continuous operator with regards to ũ, equation ( 27) is a well-defined nonlocal problem, which has a unique solution v t due to the assumptions (P0).By further using assumptions (P0), we can conclude for u t , that there exists a C 2 > 0, such that where we used in the last step, that ||f

Assumption (H1):
Since u t and v t are bounded for t ∈ [0, T ], then for every sequence {t n } n∈N with t n → 0 there exist subsequences {t n k } k∈N and {t n l } l∈N , such that there exist functions In the following we use that for functions Case 1: Proof of (H1) for square integrable kernels Since φ ij is essentially bounded on Ω×Ω, φ iI is essentially bounded on Ω×Ω I and ξ t is continuous and therefore bounded on Ω, we can conclude that and by using the dominated convergence theorem we get Thus, by (28) we derive Due to the continuity of parameter integrals, we have Analogously, we can show

So we can conclude lim
Since the solution is unique we derive q 1 = u 0 and u t u 0 .Similarly, we have for q 2 and for all ũ ∈ L 2 c (Ω) So we conclude q 2 = v 0 and v t v 0 (t 0).By using the mean value theorem, there exist s t ∈ (0, t), s.t.s t → 0(t 0) and Therefore we now prove assumption (H1) by showing lim Computing the derivative regarding t yields First we can show By applying [54, Lemma 2.16], we obtain ∇f Moreover for every t ∈ [0, T ) the derivative d dr r=t + ξ r = d dr r=t + det(I + rDV) is continuous in r and d dr r=0 + ξ r = div V(see e.g.[42]), so we derive f t Again by using (28), we obtain the convergence in (29).Furthermore, we now employ representation (7) of the nonlocal bilinear form A to compute the partial derivative of A regarding t ∂ t A(t, u 0 , v s ) = ∂ t A(t, u 0 , v s ) = Ω Ω∪Ω I v s (x) u 0 (x)∇ x γ t Γ (x, y) − u 0 (y)∇ y γ t Γ (y, x) V(x)ξ t (x)ξ t (y) dydx.

Figure 5 . 5 :
Figure 5.5:In the first six or seven iterations the improvement regarding the objective function value is quite high.After that the objective function value decreases in a much slower fashion.Due to the regularization term the objective functional will not converge to zero.