Approximate Riesz Representatives of Shape Gradients

. We study ﬁnite element approximations of Riesz representatives of shape gradients. First, we provide a general perspective on its error analysis. Then, we focus on shape functionals constrained by elliptic boundary value problems and H 1 -representatives of shape gradients. We prove linear convergence in the energy norm for linear Lagrangian ﬁnite element approximations. This theoretical result is conﬁrmed by several numerical experiments.


Introduction
A shape functional is a real map defined on a set of admissible shapes.The goal of shape optimization is to modify an initial shape so that a shape functional attains an extremal value.A common approach is to employ steepest descent algorithms [8,Ch. 3.4].Shapes may be parameterized by C 1 -mappings acting on reference configurations.Then the shape gradient is a linear continuous operator on the non-reflexive Banach space C 1 , and the concept of steepest descent may not be well-defined; see [7,Page 103].A compromise is to replace "steepest descents" with Riesz representatives of shape gradients with respect to a Hilbert space X.Henceforth, we refer to these representatives as X-representatives.
After recalling basic definitions of shape calculus, we provide a general perspective on error analysis in the energy norm for finite element approximations of X-representatives of shape gradients.Then, we zero in on shape functionals constrained to elliptic boundary value problems.For this case, insight into shape Hessians [12,14] suggests to select representatives of shape gradients with respect to X = H 1 0 (D), where D is a hold-all domain that encloses the initial guess Ω; see Figure 1.For the choice X = H 1 0 (D), it is natural to consider discretization by means of linear Lagrangian finite elements [2,14,15].We show that linear Lagrangian finite element approximations of H 1 -representatives of shape gradients converge linearly with respect to the mesh width.Additionally, this convergence rate does not deteriorate when state and adjoint variables are replaced by linear Lagrangian finite elements solutions.This is an improvement on the result presented in [2], which involves approximations of state and adjoint variables with quadratic finite elements.Finally, we provide numerical evidence of the linear convergence rate predicted.Ω D Fig. 1.The hold-all domain D encloses the domain Ω.

Shape Functionals and Shape Gradients
Let Ω ⊂ R d , d = 2, 3, be an open bounded domain with piecewise smooth boundary ∂Ω, and let J (Ω) ∈ R be a real-valued quantity of interest associated to it.One is often interested in shape sensitivity, which quantifies the impact of small perturbations of ∂Ω on the value J (Ω).
We model perturbations of the domain Ω through maps of the form where The value J (Ω) is interpreted as the realization of a shape functional, a real map The sensitivity of J (Ω) with respect to the perturbation direction V is given by the Eulerian derivative of the shape functional J in the direction V, that is, We say that the shape functional is shape differentiable if Formula (2) defines a linear and bounded operator V → dJ (Ω; V).In literature, this operator is called shape gradient [9, Ch. 9, Sect.3.4].As mentioned in the introduction, Xrepresentatives of shape gradients can be employed to solve shape optimization problems, that is, to find where U ad denotes a set of admissible shapes.
Often, the quantity of interest takes the form where the state function u is the solution of a boundary value problem stated on Ω, B ⊂ Ω, α and β are two real constants and g is a sufficiently smooth target function.In this work, u ∈ H 1 0 (Ω) is the (weak) solution of the elliptic boundary value problem with homogeneous Dirichlet boundary conditions that is, where f ∈ H 1 (Ω).For the sake of brevity, we set g = 0.Then, the shape gradient of the shape functional associated to (3) and constrained to (5) reads [4, Formula (2.9)] where the adjoint function p ∈ H 1 0 (Ω) is the solution of Formula ( 3) is a prototypical PDE-constrained shape functional.In this work, Formula ( 6) is used as test case for proving convergence estimates and performing numerical experiments.
Remark 1. Formula (6) holds even if homogeneous Dirichlet boundary conditions in (4) are replaced by homogeneous Neumann boundary conditions, in which case the test and the trial spaces in ( 5) and ( 7) are replaced with H 1 (Ω).
Remark 2. For the sake of simplicity, we restrict our considerations to homogeneous boundary conditions.However, we expect that the results of this work hold true for (sufficiently regular) inhomogeneous boundary conditions, too.Note that Formula (6) should be adjusted accordingly; see [4, Section 2].

The General Case
Let (•, •) X denote the inner product of a Hilbert space X, and let us assume that the shape gradient dJ is well-defined on X.The X-representative V X of dJ can be computed by solving Next, for an index set N , we introduce a family {X n } n∈N of finite-dimensional subspaces of X.Let V Xn n∈N be a sequence of approximate X-representatives of dJ defined by By and large, the shape gradient of a PDE-constrained shape functional depends also on the state and the adjoint variables u and p.These functions are solutions of boundary value problem.Usually, only numerical approximations u h and p h are available.In that case, the approximate X-representative V Xn has to be replaced with the solution where dJ h is an approximation of the operator dJ obtained by replacing the functions u and p with their numerical approximations u h and p h .By Strang Lemma [11, Thm 4.1.1],the estimate (8) should be corrected by adding a consistency term, that is, for a constant C > 0 independent of n and h.

H 1 -representatives and Linear Lagrangian Finite Elements
A popular approach in shape optimization consists of replacing the initial domain Ω with a polygon/polyhedron equipped with a finite element mesh Ω h .This mesh is used to compute linear Lagrangian finite element approximations of the functions u and p.Then, the coordinates of the mesh nodes are (iteratively) updated according to the shape gradient [8,Ch. 6.5].This is equivalent to extending Ω h to a mesh D h that covers a hold-all domain D and choosing linear Lagrangian finite elements to construct the finite-dimensional subspace X n .Formula (10), standard finite element estimates, and Proposition 1 readily imply that, for this discretization, the approximate H 1 -representative of (6) satisfies which is the main result of this work.
Proposition 1.Let Ω ⊂ R d be a polyhedral domain, let f ∈ W 1,4 (Ω) in (5), and let us assume that the solution u of (5) satisfies Let (V h ) h∈(0,1] be a family of H 1 0 (Ω)-conforming piecewise linear Lagrangian finite element spaces built on a quasi-uniform family of simplicial meshes (T h ) h∈(0,1] , that is, a family of meshes such that max{diam(T ) : for a ρ > 0, where B T is the largest ball contained in the simplex T [10, Def.4.4.13].Let u h , p h ∈ V h be solutions of where α, β ∈ R, B ⊂ Ω, and α = 0 or B = Ω, if d = 3.Let dJ h (Ω; W n ) denote the operator defined by Formula (6) with u and p replaced by u h and p h , respectively.Then, for a constant C(Ω, f, u, p) > 0 independent of n and h.
Proof.First of all, note that We recall that, for generic functions q 0 ∈ L 2 (Ω) and q 1 , q 2 ∈ L 4 (Ω), the Cauchy-Schwarz inequality implies Thus, the first integral in ( 16) may be estimated as follows 1| Ω (∇f The second integral in ( 16) may be estimated as follows The third and the fourth integral in ( 16) may be estimated similarly.
Stability of the Ritz projection with respect to which in turn implies p − p h W 1,4 (Ω) = O(h), it is necessary to repeat the proof of (18) given in [3] tracking the consistency term The discrete Green's function g z h ∈ V h is given in [3] and satisfies ) .By the Cauchy-Schwarz inequality and standard finite element estimates, The stability result (19) holds if ( 21) is bounded independently of h.For this reason, we need to set α = 0 when d = 3, unless B = Ω.In this latter case, by Galerkin orthogonality, (20) is bounded by Remark 3. In Proposition 1, we assume W2,4 -regularity of the solution u of ( 5).
This assumption is made to achieve linear convergence with respect to h in the estimate (15).However, a three-dimensional polyhedral domain must satisfy tight geometric conditions for u to be in W 2,4 [6, Thm 7.1].Nevertheless, in [5] the authors show W 1,∞ -stability of the Ritz projection for general convex polyhedral domains.Therefore, we expect that (in the latter case) the righthand side of ( 15) can be replaced with a term of order O(h α ), where the rate α depends on the regularity of u and satisfies 0 < α ≤ 1.
Remark 4. In [4,13], the authors show that one can expect superconvergence in the approximation of the shape gradient dJ .In particular, they show that However, in the right-hand side of ( 22) appears the W 2,4 (Ω)-norm of W. But to prove convergence in the approximation of a H 1 -representative of dJ , the upper bound of |dJ (Ω; W) − dJ h (Ω; W)| cannot involve a norm stronger than the H 1 -norm; see Equation (10).
Remark 5.By the Hadamard structure theorem [9, Ch. 9, Thm 3.6], most shape gradients admit representatives g(Ω) in the space of distributions D k (∂Ω), that is, where γ ∂Ω V • n is the normal component of V on the boundary ∂Ω.For instance, if u, p ∈ H 2 (Ω), Formula ( 6) is equivalent to [4, Formula (2.10)] We advise against the use of g(Ω) (which corresponds to the L 2 (∂Ω)-representative of dJ ) to define descent directions because L 2 -representatives might bristle with undesirable oscillations [1].

Numerical Experiments
We provide numerical evidence of the estimate (11).We employ linear Lagrangian finite elements on quasi-uniform triangular meshes.The experiments are performed in MATLAB and are partly based on the library LehrFEM developed at ETHZ.Mesh generation and uniform refinement are performed with the functions initmesh and refinemesh of the MATLAB PDE Toolbox [16].The boundary of computational domains is approximated by a polygon, which is generally believed not to affect the convergence of linear finite elements [10, Sect.10.2].For domains with curved boundaries, the refined mesh is always adjusted to fit the boundary.Integrals in the domain are computed with a 3-point quadrature rule of order 3 in each triangle and line integrals with a 6-point Gauss quadrature on each segment.We consider three different geometries for the domain Ω (see Figure 2): 1.A disc of radius 6/5 centered in (0.01,0.02).
The hold-all domain D is a square with edges of length 3 centered in the origin.The region of interest B is the whole domain Ω.We set α = 0 and β = 1 in (3).
Fig. 2. The domain Ω is chosen to be either a disc or a triangle or a sector.
The reference value V X is approximated by computing V Xn h on a mesh with an extra level of refinement.In light or Remark 5, we employ both Formula (6) and Formula (24) to evaluate the right-hand side dJ h in (9).To avoid biased results, we display convergence history of V X − V Xn h H 1 (D) both with self-and cross-comparison.
In Figure 3, we plot the convergence history when the domain Ω is either a disc (first row) or a triangle (second row).As predicted by (11), we observe linear convergence when the right-hand side in ( 9) is evaluated according to (6).Interestingly, using Formula (24) seems not to affect the convergence rate.The same behavior is observed when homogeneous Dirichlet boundary conditions are replaced by homogeneous Neumann boundary conditions.Note that, in this latter case, the boundary-integral counterpart of Formula (6) reads [4] For the sake of brevity, we omit these plots.
In Figure 4 (first row), we plot the convergence history when the domain Ω is a sector.This domain does not guarantee that u and p are in H 2 (Ω) because it has a re-entrant corner.We observe that the convergence rates decrease to fractional values.This is a consequence of the lower regularity of the functions u and p.Additionally, the convergence rates depend on the formula used to evaluate dJ h .In particular, in the cross-comparison, the convergence line saturates when Formula ( 6) is used.This may be due to a poor accuracy of the reference solution.However, we point out that Formulas ( 6) and (24) may not be equivalent due to the lack of regularity of the functions u and p; cf.Remark 5. Curiously, for homogeneous Neumann boundary conditions, the presence of the re-entrant corner seems to have a milder impact on convergence rates; see Figure 4 (second row).However, note that the approximate algebraic convergence rates of u − u h H 1 (Ω) and p − p h H 1 (Ω) Fig. 4. Convergence history when Ω is a sector.Line refers to evaluation of dJ h according to Formula (6); line to Formula (24) (in the first row ) and to Formula (25) (in the second row ).For Dirichlet boundary conditions, convergence rates decay to fractional values.

Conclusion
Most shape optimization algorithms rely on Riesz representatives of shape gradients with respect to a chosen Hilbert space.Numerical discretization is inevitable when the shape functional is constrained to a boundary value problem.Formula (10) indicates how to estimate the discretization error when the Riesz representative is computed on a finite-dimensional trial space and the shape gradient can be evaluated only approximately.
For linear Lagrangian approximations of H 1 -representatives, Proposition 1 implies that the discretization error decays linearly with respect to the mesh width h.This convergence behavior is observed in several numerical experiments.
As a consequence of the Hadamard structure theorem, most shape gradients can be equivalently formulated as boundary or volume integrals.Although Proposition 1 relies on the volume formulation of the shape gradient, we have observed linear convergence independently of the formula employed to evaluate dJ .However, we advise to rely on the volume-based formula because it imposes lower regularity assumptions on the state and the adjoint variables [4,9,15].