A Cut Finite Element Method for the Bernoulli Free Boundary Value Problem
Erik Burman, Daniel Elfverson, Peter Hansbo, Mats G. Larson, Karl Larsson
AA Cut Finite Element Method for theBernoulli Free Boundary Value Problem ∗ Erik Burman † Daniel Elfverson ‡ Peter Hansbo § Mats G. Larson ¶ Karl Larsson (cid:107)
October 8, 2018
Abstract
We develop a cut finite element method for the Bernoulli free boundary problem.The free boundary, represented by an approximate signed distance function on afixed background mesh, is allowed to intersect elements in an arbitrary fashion. Thisleads to so called cut elements in the vicinity of the boundary. To obtain a stablemethod, stabilization terms is added in the vicinity of the cut elements penalizing thegradient jumps across element sides. The stabilization also ensures good conditioningof the resulting discrete system. We develop a method for shape optimization basedon moving the distance function along a velocity field which is computed as the H Riesz representation of the shape derivative. We show that the velocity field is thesolution to an interface problem and we prove an a priori error estimate of optimalorder, given the limited regularity of the velocity field across the interface, for thethe velocity field in the H norm. Finally, we present illustrating numerical results. Keywords.
Free boundary value problem; CutFEM; Shape optimization; Level set; Fic-titious domain method ∗ This research was supported in part by the Swedish Foundation for Strategic Research Grant No.AM13-0029, the Swedish Research Council Grants Nos. 2011-4992, 2013-4708, the Swedish Research Pro-gramme Essence, and EPSRC, UK, Grant Nr. EP/J002313/1. † Department of Mathematics, University College London, Gower Street, London WC1E 6BT, UK,[email protected] ‡ Department of Mathematics and Mathematical Statistics, Ume˚a University, SE–901 87 Ume˚a, Swe-den,[email protected] § Department of Mechanical Engineering, J¨onk¨oping University, SE-551 11 J¨onk¨oping, Sweden, [email protected] ¶ Department of Mathematics and Mathematical Statistics, Ume˚a University, SE–901 87 Ume˚a, Sweden,[email protected] (cid:107)
Department of Mathematics and Mathematical Statistics, Ume˚a University, SE–901 87 Ume˚a, Sweden,[email protected] a r X i v : . [ m a t h . NA ] S e p Introduction
In this paper we consider the application of the recently developed cut finite elementmethod (CutFEM) [4, 6] to the Bernoulli free boundary problem. This problem appearsin a variety of applications such as stationary water waves (Stokes waves) and the optimalinsulation problem. The Bernoulli free boundary problem is very well understood fromthe mathematical point of view, see [3, 14, 26] and the references therein, and it alsoserves as a standard test problem for different optimization and numerical methods, seee.g. [18, 21, 22] among others. When numerically solving free boundary problems it ishighly beneficial to avoid updating the computational mesh when updating the boundary,since large motions of the boundary may require complete remeshing. This can be achievedby the use of fictitious domain methods, which, however, are known not to perform verywell for shape optimization or free surface problems due to their lack of accuracy close tothe boundary [17]. An exception is the least squares formulation suggested in [13], wheresimilarly as in [6] the formulation is restricted to the physical domain.A fictitious domain method which does not lose accuracy close to the boundary isthe recently developed CutFEM method, see [4]. CutFEM uses weak enforcement of theboundary conditions and a sufficiently accurate representation of the domain together withcertain consistent stabilization terms to guarantee stability, optimal accuracy, and condi-tioning independent of the position of the boundary in the background mesh. Furthermore,no expensive and complicated mesh operations (edge split, edge collapse, remeshing, etc.)need to be performed when updating the boundary. This is a significant gain especially forcomplicated boundaries and for 3D applications. CutFEM has successfully been appliedfor problems with unknown or moving boundaries in, e.g., [8, 16].As in [1, 2] we consider a shape optimization approach to solve the Bernoulli freeboundary problem using sensitivity analysis and a level set representation [24] to trackthe evolution of the free boundary. In the sensitivity analysis we do not use the standard
Hadamard structure of the shape functional, i.e., we do not express the shape functionalas a normal perturbation of the boundary. Instead we use a volume representation whichrequires less smoothness and has proved to possess certain superconvergence propertiescompared to the boundary formulation [20], see also [19, 23]. To obtain a velocity fieldfrom the shape derivate, we use the Hilbertian regularization suggested in [12], where weessentially let the velocity field on the domain be defined as the solution to the weak ellipticproblem associated with the H inner product with right hand side given by the shapederivative functional. We may thus view the velocity field as the H Riesz representationof the shape derivative functional acting on H . This procedure leads to an elliptic interfaceproblem for the velocity field. The free boundary is then updated by moving the level setalong the velocity field.We derive a priori error estimates for the CutFEM approximation of the primal problemand the dual problem, involved in the computation of the shape derivative in the W p ,2 ≤ p < ∞ , norm, and then we use these estimates to prove an a priori error estimatefor the discrete approximation of the velocity field in the H norm. In the error estimatesfor the primal and dual problems we use inverse estimates, for 2 < p < ∞ , which leads to2uboptimal convergence rates, but it turns out these bounds are indeed sharp enough toprove optimal order estimates for the discrete velocity field.An outline of the paper is as follows: in Section 2 we present the model problemand CutFEM discretization, in Section 3 we use sensitivity analysis to derive the shapederivative, in Section 4 we discuss how to compute a regularized descent direction from theshape derivative, in Section 5 we present a level set representation of the free boundary andmethod for computing its evolution, in Section 6 we present an optimization algorithm, inSection 7 we present a priori error estimates of the CutFEM approximation of the primaland dual problems as well as for the discrete approximation of the velocity field, and finally,in Section 8 we present numerical experiments to verify the convergence rates and overallbehavior of the optimization algorithm. We consider the Bernoulli free boundary value problem: − ∆ u = f in Ω (2.1) u = g D on ∂ Ω (2.2) n · ∇ u = g N on Γ (2.3)where Γ ⊂ ∂ Ω is the free and ∂ Ω \ Γ is the fixed part of the boundary, g D = 0 and g N is constant on Γ. Note the double boundary conditions on Γ. We seek to determine thedomain Ω such that there exist a solution u to (2.1)-(2.3). For f = 0 we refer to [3, 14, 26]and the references therein for theoretical background of the Bernoulli free boundary valueproblem and for f = 0 we assume that f is such that there exist a unique solution.In order to obtain a formulation which is suitable as a starting point for a numericalalgorithm we recast the overdetermined boundary value problem as a constrained mini-mization problem as follows. We seek to minimize the functional J (Ω) := J (Ω; u (Ω)) = min Ω (cid:90) Γ u dΓ (2.4)where the function u solves the boundary value problem − ∆ u = f in Ω (2.5) u = g D on Γ fix := ∂ Ω \ Γ (2.6) n · ∇ u = g N on Γ (2.7)Note that we keep the Neumann condition on the free boundary and enforce the Dirichletcondition through the minimization of the functional J (Ω).3he weak formulation of (2.5-2.7) reads: find u ∈ V g D (Ω) := { v ∈ H (Ω) : v | Γ fix = g D } such that a (Ω; u, v ) := (cid:90) Ω ∇ u · ∇ v dΩ = (cid:90) Ω f dΩ =: F (Ω; v ) ∀ v ∈ V (Ω) (2.8)Let O be the set of admissible domains; then the constrained minimization problem reads:find Ω ∈ O such that J (Ω) = min Ω ∈O J (Ω; v ) (2.9)for u ∈ V g D (Ω) s.t. a (Ω; u, v ) = F (Ω; v ) ∀ v ∈ V (Ω) (2.10)To solve the minimization problem (2.9)-(2.10), we define the corresponding Lagrangianas L ( ω, v, q ) := J ( ω ; u ) − a ( ω ; v, q ) + F ( ω ; q ) (2.11)That is, we seek the domain Ω such thatΩ = arg min { ω ∈ O | min v ∈ V gD ( ω ) max q ∈ V ( ω ) L ( ω, v, q ) } (2.12) We will use a cut finite element method to discretize the boundary value problem (2.8).Before formulating the method we introduce some notation.
The Mesh and Finite Element Space.
Let Ω be a polygonal domain such that alladmissible domains Ω ∈ O are subsets of Ω , i.e. Ω ⊂ Ω . Let T h, denote a family of quasi-uniform triangulations of Ω with mesh parameter h ∈ (0 , h ] and define the correspondingspace of continuous piecewise linear polynomials V h (Ω ) = { v ∈ H (Ω ) : v | T ∈ P ( T ) , ∀ T ∈ T h, } (2.13)Given Ω ∈ O we define the active mesh T h = { T ∈ T h, : T ∩ Ω (cid:54) = ∅} (2.14)the union of the active elements Ω h = ∪ T ∈T h T (2.15)and the finite element space on the active mesh V h (Ω) = V h (Ω ) | Ω h (2.16)Let also F h denote the set of interior faces in T h such that at least one of its neighboringelements intersect the boundary ∂ Ω, F h = { F : T + F ∩ ∂ Ω (cid:54) = ∅ or T − F ∩ ∂ Ω (cid:54) = ∅} (2.17)where, T + F and T − F are the two elements sharing the face F . On a face F we define thejump (cid:74) v (cid:75) = v | T + F − v | T − F (2.18)where T + F is the element with the higher index.4 he Method. We define the forms A h (Ω; v, w ) = a h (Ω; v, w ) + s h ( v, w ) (2.19) a h (Ω; v, w ) = ( ∇ v, ∇ w ) Ω − ( ∂ n v, w ) Γ fix − ( ∂ n w, v ) Γ fix + ( γ D h − v, w ) Γ fix s h (Ω; v, w ) = (cid:88) F ∈F h ( γ h (cid:74) ∂ n v (cid:75) , (cid:74) ∂ n w (cid:75) ) F (2.20) F h (Ω; w ) = ( f, w ) Ω + ( g D , γ D h − w − ∂ n w ) Γ fix + ( g N , w ) Γ (2.21)were ( u, v ) ω := (cid:82) ω u · v d ω is the L inner product over the set ω equipped with theappropriate measure. Our method for the approximation of (2.5)–(2.7) takes the form:find u h ∈ V h (Ω) such that A h (Ω; u h , v ) = F h (Ω; v ) ∀ v ∈ V h (Ω) (2.22)We recognize the weak enforcement of Dirichlet boundary conditions by Nitsche’s method,cf. [15]. Furthermore, the term s h , first suggested in this context in [6], is added to stabilizethe method in the vicinity of the boundary. For O ∈ O we let W (Ω , R d ) denote the space of sufficiently smooth vector fields and for avector field ϑ ∈ W we define the map M ϑ : Ω × I (cid:51) ( x, t ) (cid:55)→ x + tϑ ( x ) ∈ M ϑ (Ω , t ) ⊂ R d (3.1)where I = ( − δ, δ ), δ >
0. For small enough δ , the mapping Ω (cid:55)→ M ϑ (Ω , t ) is a bijectionand M ϑ (Ω ,
0) = Ω. We also assume that the vector field ϑ is such that M ϑ (Ω , t ) ∈ O for t ∈ I with δ small enough.Let J (Ω) be a shape functional, i.e., a mapping J : O (cid:51) Ω (cid:55)→ J (Ω) ∈ R . We then havethe composition I (cid:51) t (cid:55)→ J ◦ M (Ω , t ) ∈ R and we define the shape derivative D Ω ,ϑ of J inthe direction ϑ by D Ω ,ϑ J (Ω) = ddt J ◦ M ϑ (Ω , t ) | t =0 = lim t → J ( M ϑ (Ω , t )) − J (Ω) t (3.2)Note that if M ϑ (Ω , t ) = Ω we have D Ω ,ϑ J = 0, even if M ϑ change points in the interior ofthe domain.We finally define the shape derivative D Ω J | Ω : W (Ω , R d ) → R by D Ω J | Ω ( ϑ ) = D Ω ,ϑ J (Ω) (3.3)In cases when the functional J depend on other arguments we use ∂ Ω to denote the partialderivative with respect to Ω and ∂ Ω ,ϑ to denote the partial derivative with respect to Ω inthe direction ϑ . 5 .2 Leibniz Formulas For v : Ω × I → R d we define the material time derivative in the direction ϑ by D t,ϑ v = lim t → v ( M ϑ ( x, t ) , t ) − v ( x, t (3.4)and the partial time derivative by ∂ t v = lim t → v ( x, t ) − v ( x, t (3.5)From the chain rule it follows that D t,ϑ v = ∂ t v + ϑ · ∇ v (3.6)The material derivative does not commute with the gradient and we have the commutator[ D t,ϑ , ∇ ] v = D t,ϑ ( ∇ v ) − ∇ ( D t,ϑ v ) = − ( Dϑ ) T ∇ v (3.7)where Dϑ = V ⊗ ∇ is the derivative (or Jacobian) of the vector field ϑ , and the usualproduct rule D t,ϑ ( vw ) = ( D t,ϑ v ) w + v ( D t,ϑ w ) (3.8)holds. To derive a expression of the shape derivative, the following lemma will be usedfrequently. Lemma 3.1.
Let f, g : R d → R be functions smooth enough for the following expressionsto be well defined. Then the following relationships hold D Ω ,ϑ (cid:90) Ω f dΩ = (cid:90) Ω ( D t,ϑ f + ( ∇ · ϑ ) f ) dΩ (3.9) D Ω ,ϑ (cid:90) Γ g dΓ = (cid:90) Γ ( D t,ϑ g + ( ∇ Γ · ϑ ) g ) dΓ (3.10) where ∇ Γ · v = ∇ · v − n · Dϑ · n and n is the unit normal to Γ .Proof. The proof can e.g. be found in [27].
Remark 3.1.
An alternative definition of the surface divergence is ∇ Γ · ϑ = tr( ϑ ⊗ ∇ Γ )where ∇ Γ v = (1 − n ⊗ n ) ∇ v is the tangent gradient. Recall the Lagrangian L ( ω, v, q ) = J ( ω ; v ) − a ( ω ; v, q ) + F ( ω ; q ) (3.11)6or a fixed Ω ∈ O we take the Fr´echet derivative D L (Ω , · , · ) : V (Ω) × V (Ω) → R of L , D L| (Ω ,u,p ) ( δu, δp ) = (cid:10) D L| (Ω ,u,p ) , ( δu, δp ) (cid:11) = (cid:10) ∂ v L| (Ω ,u,p ) , δu (cid:11)(cid:124) (cid:123)(cid:122) (cid:125) =0 + (cid:10) ∂ q L| (Ω ,u,p ) , δp (cid:11)(cid:124) (cid:123)(cid:122) (cid:125) =0 = 0 (3.12)for any direction δu, δp ∈ V (Ω) and use the following identities (cid:104) ∂ v L| (Ω ,u,p ) , δu (cid:105) = (cid:104) ∂ v J (Ω; u ) , δu (cid:105) − a (Ω , δu, p ) = 0 (3.13)which hold if p solves the dual problem a (Ω , v, p ) = m ( v ) := (cid:104) ∂ v J (Ω; u ) , v (cid:105) ∀ v ∈ V (Ω) (3.14)where (cid:104) ∂ v J (Ω; u ) , v (cid:105) = ( u, v ) Γ and (cid:104) ∂ q L| (Ω ,u,p ) , δp (cid:105) = F (Ω , δp ) − a (Ω , u, δp ) = 0 (3.15)which holds since u is the solution to the primal problem (2.8). The Correa-Seeger theorem[10] states that D Ω ,ϑ (cid:18) min v ∈ V gD (Ω) max q ∈ V (Ω) L (Ω , v, q ) (cid:19) = ∂ Ω ,ϑ L (Ω , u, p ) (3.16)and we obtain the shape derivative D Ω ,ϑ J (Ω) = ∂ Ω ,ϑ L (Ω , u, p ) (3.17) Lemma 3.2.
For u and p solving (2.8) and (3.14) , respectively, the shape derivative of L ( ω, v, q ) in the point (Ω , u, p ) is given by ∂ Ω ,ϑ L| (Ω ,u,p ) = − (cid:90) Ω ∇ u · ( Dϑ + ( Dϑ ) T ) · ∇ p dΩ (3.18)+ (cid:90) Ω ( ϑ · ∇ f ) p dΩ+ (cid:90) Ω ( ∇ · ϑ ) ( f − ∇ u · ∇ p ) dΩ+ (cid:90) Γ ( ∇ Γ · ϑ ) (cid:0) − u + g N p (cid:1) dΓ Proof.
This is a well known result but we include the proof for the convenience of thereader. We have ∂ Ω ,ϑ L (Ω , u, p ) = ∂ Ω ,ϑ J (Ω; u ) − ∂ Ω ,ϑ a (Ω , u, p ) + ∂ Ω ,ϑ F (Ω , p ) (3.19)7sing Lemma 3.1 we obtain the identities ∂ Ω ,ϑ J (Ω; u ) = (cid:90) Γ
12 ( ∇ Γ · ϑ ) u dΓ ∂ Ω ,ϑ a (Ω; u, p ) = (cid:90) Ω ( D t,ϑ ( ∇ u · ∇ p ) + ( ∇ · ϑ ) ∇ u · ∇ p ) dΩ= (cid:90) Ω (cid:0) ∇ u · ( Dϑ + ( Dϑ ) T ) · ∇ p + ( ∇ · ϑ ) ∇ u · ∇ p (cid:1) dΩ ∂ Ω ,ϑ F (Ω; p ) = (cid:90) Γ ( D t,ϑ ( g N p ) + ( ∇ Γ · ϑ )( g N p )) dΓ + (cid:90) Ω ( D t,ϑ ( f p ) + ( ∇ · ϑ )( f p )) dΩ= (cid:90) Γ ( ∇ Γ · ϑ )( g N p ) dΓ + (cid:90) Ω ( ∇ · ϑ )( f p ) dΩInserting these expressions into (3.19) we arrive at (3.18). In order to compute an approximation of the shape derivatives we need pproximationsof the solutions to the primal equation (2.8) and the dual equation (3.14). We employCutFEM formulations: find u h ∈ V h (Ω) such that A h (Ω; u h , w ) = F h (Ω; w ) ∀ w ∈ V h (Ω) (3.20)and p h ∈ V h (Ω) such that A h (Ω; p h , w ) = m h ( w ) := ( u h , w ) Γ ∀ w ∈ V h (Ω) (3.21)The discrete approximation of the shape derivative is obtained by inserting the discretequantities u h , p h into (3.18), i.e., D Ω ,ϑ L| (Ω ,u h ,p h ) = − (cid:90) Ω ∇ u h · ( Dϑ + ( Dϑ ) T ) · ∇ p h dΩ (3.22)+ (cid:90) Ω ( ϑ · ∇ f ) p h dΩ+ (cid:90) Ω ( ∇ · ϑ ) ( f − ∇ u h · ∇ p h ) dΩ+ (cid:90) Γ ( ∇ Γ · ϑ ) (cid:0) − u h + g N p h (cid:1) dΓ We follow [12] to define a velocity field β given the shape derivative. We seek the velocityfield β such that we obtain the largest decreasing direction of D Ω ,β L| (Ω ,u,p ) under some8egularity constraint, for instance assume that the velocity field is in [ H (Ω )] d we obtain β = arg min (cid:107) ϑ (cid:107) [ H d =1 D Ω ,ϑ L| (Ω ,u,p ) (4.1)Let b be the [ H (Ω )] d inner product b ( v, w ) := (cid:90) Ω ( Dv : Dw + v · w ) dΩ (4.2)An equivalent formulation of the minimization problem (4.1) is: find β (cid:48) ∈ [ H (Ω )] d suchthat b ( β (cid:48) , ϑ ) = − D Ω ,ϑ L| (Ω ,u,p ) ∀ ϑ ∈ [ H (Ω )] d (4.3)and set β = β (cid:48) (cid:107) β (cid:48) (cid:107) [ H (Ω )] d (4.4)It is then clear that β is a descent direction since D Ω ,β L| (Ω ,u,p ) = − b ( β (cid:48) , β ) = −(cid:107) β (cid:48) (cid:107) [ H (Ω )] d ≤ Remark 4.1.
To prove the equivalence between the minimization problem (4.1) and (4.3)–(4.4) we compute the saddle point to the Lagrangian corresponding to (4.1). We obtainthe Lagrangian K ( τ, λ ) = D Ω ,τ L| (Ω ,u,p ) + λ (cid:16) ( τ, τ ) [ H (Ω )] d − (cid:17) (4.6)In the saddle point ( τ, λ ), we have that0 = (cid:68) ∂∂τ K ( τ, λ ) , φ (cid:69) = D Ω ,φ L| (Ω ,u,p ) + 2 λ ( τ, φ ) [ H (Ω )] d (4.7)and 0 = ∂∂λ K ( τ, λ ) = ( τ, τ ) [ H (Ω )] d − β (cid:48) = 2 λτ where β = β (cid:48) (cid:107) β (cid:48) (cid:107) [ H (Ω )] d = 2 λτ (cid:107) λτ (cid:107) [ H (Ω )] d = τ (4.9)and λ = (cid:107) β (cid:48) (cid:107) [ H (Ω )] d /
2. Hence the two formulations are equivalent.9 .2 Regularity of the Velocity Field
Next we investigate the regularity of the velocity field β . For smooth domains and understronger regularity requirements the shape derivative can be formulated using Hadamard’sstructure theorem as an integral over the boundary, D Ω ,ϑ L (Ω , u, p ) = ( G , n · ϑ ) L (Γ) (4.10)where G is a function of the primal and dual solutions u and p , the right hand side, theboundary condition, and the mean curvature. We thus note that the velocity field β is asolution to the problem: find β ∈ [ H (Ω )] d such that b ( β, ϑ ) = − ( G , n · ϑ ) L (Γ) ∀ ϑ ∈ [ H (Ω )] d (4.11)The corresponding strong problem for each of the components β i , 1 = 1 , . . . , d of β is − ∆ β i = 0 , in Ω \ Γ (4.12) β i = 0 , on ∂ Ω (4.13)[ β i ] = 0 , on Γ (4.14)[ n · ∇ β i ] = G n i , on Γ (4.15)which is an interface problem. Given that Γ is smooth and G ∈ H / (Γ), we have theregularity estimate (cid:107) β i (cid:107) H (Ω ) + (cid:107) β i (cid:107) H (Ω \ Γ) (cid:46) (cid:107)G(cid:107) H / (Γ) (4.16)see [9], and hence β ∈ [ H (Ω )] d ∩ [ H (Ω \ Γ)] d . We define a discrete velocity field using a standard finite element discretization of (4.3) withpiecewise linear continuous trial and test functions V h (Ω ) on Ω . The discrete problemtakes the form: find β (cid:48) h ∈ [ V h (Ω )] d such that b ( β (cid:48) h , ϑ ) = − D Ω ,ϑ L| (Ω ,u h ,p h ) ∀ ϑ ∈ [ V h (Ω )] d (4.17)and set β h = β (cid:48) h (cid:107) β (cid:48) h (cid:107) [ H (Ω )] d (4.18) A level set function describing an interface needs to be evolved in order to find a minimumto (2.9)-(2.10). Let ρ ( x, Γ) be a distance function defined as the minimal Euclidean distance10etween x and Γ. The level set function is the signed distance function φ ( x ) = ρ ( x, Γ) x ∈ Ω \ Ω0 x ∈ Γ − ρ ( x, Γ) x ∈ Ω (5.1)This function is moved by solving a Hamilton-Jacobi equation of the form ∂ t φ + β · ∇ φ = 0 (5.2)After some time φ no longer resembles a discrete signed distance function and so calledreinitialization needs to be performed to restore the distance properties. Reinitializationcan be done by solving the Eikonal equation (cid:40) ∂ t ϕ + sign( φ )( |∇ ϕ | −
1) = 0 t ∈ (0 , T ] ϕ = φ t = 0 (5.3)for the unknown ϕ . Setting φ = lim T →∞ ϕ ( · , T ) yields a signed distance function on Ω. Inthe present paper we use a fast sweeping method to approximate (5.3) as suggested in [11]. To evolve the interface we use a standard finite element discretization of (5.2), using thespace V h (Ω ) of continuous piecewise linear elements on Ω , with symmetric interior penaltystabilization, see [5], in space and a Crank-Nicolson scheme in time. Given a time t wefirst determine a suitable time T such that we may use β h ( t ) as an approximation of β h ( t )on the interval [ t , t + T ), then we divide [ t , t + T ) into N Crank-Nicolson steps of equallength. This procedure is repeated until a stopping criteria is satisfied.The time T may be estimated using Taylor’s formula L ( M β h (Ω , ∆ t ) , u h ◦ M − β h (Ω , ∆ t ) , p h ◦ M − β h (Ω , ∆ t )) ≈ L (Ω , u h , p h ) | t = t + D t,β h L (Ω , u h , p h ) | t = t ∆ t (5.4)Given a damping parameter α ∈ [0 ,
1) we set L (Ω t,β h , u h ◦ M − β h (Ω , t ) , p h ◦ M − β h (Ω , t )) = α L (Ω , u h , p h ) which yields the estimate T = ( α − L (Ω , u h , p h ) D t,β L (Ω , u h , p h ) (5.5)To formulate the finite element method we divide [ t , t + T ) into N time steps [ t n − , t n ),of equal length k = T /N and we use the notation φ nh = φ h ( t n ) for the solution at time t n .Given φ h ∈ V h (Ω ), find φ nh ∈ V h (Ω ) for n = 1 , . . . , N, such that (cid:18) φ nh − φ n − h k , w (cid:19) Ω + (cid:18) β h ( t ) · ∇ φ nh + φ n − h , w (cid:19) Ω + r h (cid:18) φ nh + φ n − h , w (cid:19) = 0 ∀ w ∈ V h (Ω ) (5.6)11here r h is the stabilization term r h ( v, w ) = (cid:88) F ∈F h, ( γ h [ ∂ n v ] , [ ∂ n w ]) L ( F ) (5.7)where γ > F h, is the set of interior faces in the background mesh T h, . In this section we summarize the optimization procedure and propose an algorithm to solve(2.9)-(2.10). During the optimization procedure we use sensitivity analysis to compute thediscrete shape derivate (3.22), see Section 3. From the discrete shape derivate we computea velocity field β h (4.17) using a H (Ω ) regularization, which corresponds to the thegreatest descent direction of the shape derivative in H (Ω ), see Section 4. The velocityfield is then used to move the level set and update the free boundary, see Section 5. Thesesteps are presented in Algorithm 1. As a stopping criterion we require that the residualindicator R Γ ( u h ) = (cid:107) u h (cid:107) Γ ≤ TOL (6.1)for some tolerance 0 < TOL.
Algorithm 1
Bernoulli free boundary value problemInput: A initial level set φ h , a damping parameter α , and a tolerance TOL.Compute primal solution u h (3.20) and the residual indicator R Γ ( u h ) (6.1) while R Γ ( u h ) > TOL do Compute the dual solution p h (3.21)Compute the discrete shape derivative D Ω ,ϑ L (Ω h ,u h ,p h ) (3.22)Compute the velocity field β h (4.17)Move the interface (5.6)Compute the primal solution u h (3.20)Compute the residual indicator R Γ ( u h ) (6.1) end while In this section we derive an a priori error estimate for the velocity field in the H (Ω ) norm.Recall that the regularity of the velocity field is given by (4.16) and thus the best possibleorder of convergence is O ( h / ) in the H (Ω ) norm and O ( h / ) in the L (Ω ) norm, sincewe use a standard finite element method to approximate the velocity field. To prove theerror estimate for the velocity field we will need bounds for the discretization error ofthe primal and dual solutions in L norms since the right hand side of the problem (4.3)12efining the velocity field is the shape derivative functional (3.18), which is a trilinear form,depending on the primal and dual solutions as well as the test function. For simplicity, wederive error estimates in L p norms for the primal and dual solutions using inverse boundsin combination with L error estimates. These bounds are of course not of optimal orderbut, in the relevant case d ≤
3, they are sharp enough to establish optimal order boundsfor the velocity field, given the restricted regularity of the velocity field. We employ thenotation a (cid:46) b to abbreviate the inequality a ≤ Cb where the constant 0 ≤ C is genericconstant independent of the mesh size. Definition of the Energy Norm.
For 2 ≤ p < ∞ we define the energy norm ||| v ||| pp,h = (cid:107)∇ v (cid:107) pL p (Ω) + h (cid:107) n · ∇ v (cid:107) pL p (Γ fix ) + h − p (cid:107) v (cid:107) pL p (Γ fix ) + h ||| v ||| pL p ( F h ) (7.1)where ||| v ||| pL p ( F h ) = (cid:88) F ∈F h (cid:107) [ n · ∇ v ] (cid:107) pL p ( F ) (7.2) An Inverse Estimate.
We have the inverse estimate: for all v ∈ V h (Ω) it holds ||| v ||| h,p (cid:46) h d (1 /p − / ||| v ||| h, (7.3)To verify (7.3) we first note that using the inverse estimate h /p (cid:107) n · ∇ v (cid:107) L p ( F ) (cid:46) (cid:107)∇ v (cid:107) L p ( T ) (7.4)where F is a face on the boundary of T , and the fact that the mesh T h (which we recallconsists of full elements) covers Ω we have ||| v ||| h,p (cid:46) (cid:107)∇ v (cid:107) L p ( T h ) + h /p − (cid:107) v (cid:107) L p (Γ fix ) (7.5)Next using the inverse estimates (cid:107)∇ v (cid:107) L p ( T ) (cid:46) h d (1 /p − / (cid:107)∇ v (cid:107) T (7.6) h /p − (cid:107) v (cid:107) L p ( F ) (cid:46) h /p − h ( d − /p − / (cid:107) v (cid:107) F (cid:46) h d (1 /p − / h − / (cid:107) v (cid:107) F (7.7)to pass from L p to L norms we obtain ||| v ||| h,p (cid:46) (cid:107)∇ v (cid:107) L p ( T h ) + h /p − (cid:107) v (cid:107) Γ fix (7.8) (cid:46) h d (1 /p − / (cid:16) (cid:107)∇ v (cid:107) T h + h − / (cid:107) v (cid:107) Γ fix (cid:17) (7.9) (cid:46) h d (1 /p − / ||| v ||| h, (7.10)where in the last step we used the estimate (cid:107)∇ v (cid:107) T h (cid:46) (cid:107)∇ v (cid:107) Ω + h / ||| v ||| F h (7.11)see [4]. 13 .2 Interpolation Definition of the Interpolation Operator.
We recall that there is an extension op-erator E : W sp (Ω) → W sp (Ω δ ), for 0 ≤ s and 1 ≤ p ≤ ∞ , where Ω δ = Ω ∪ U δ (Γ) with U δ (Γ)the tubular neighborhood { x ∈ R d : ρ ( x, Γ) < δ } . For h ∈ (0 , h ], with h small enough,we have Ω ⊂ Ω h ⊂ Ω δ . Let π h : L (Ω h ) → V h (Ω) be a Scott-Zhang type interpolationoperator, see [25], and for u ∈ L (Ω) we define π h v = π h ( Ev ). For convenience we will usethe simplified notation v = Ev on Ω δ . Interpolation Error Estimates.
We have the elementwise interpolation estimate h − (cid:107) v − π h v (cid:107) L p ( T ) + (cid:107)∇ ( v − π h v ) (cid:107) L p ( T ) (cid:46) h (cid:107) v (cid:107) W p ( N ( T )) (7.12)where N ( T ) is the the set of neighboring elements in T h to element T . Summing over theelements and using the stability of the extension operator we obtain the interpolation errorestimate h − (cid:107) v − π h v (cid:107) L p (Ω h ) + ||| v − π h v ||| p,h (cid:46) h (cid:107) v (cid:107) W p (Ω h ) (cid:46) h (cid:107) v (cid:107) W p (Ω) (7.13)We also have the following interpolation error estimate in the energy norm ||| v − π h v ||| h,p (cid:46) h (cid:107) v (cid:107) W p (Ω) (7.14) Verification of (7.14).
Using the element wise trace inequality (cid:107) w (cid:107) pL p ( F ) (cid:46) h − (cid:107) w (cid:107) pL p ( T ) + h p − (cid:107)∇ w (cid:107) pL p ( T ) (7.15)where F is a face on ∂T , to estimate the terms on Γ fix , h (cid:107) n · ∇ w (cid:107) pL p (Γ fix ) (cid:46) (cid:107)∇ w (cid:107) pL p ( T h (Γ fix )) + h p (cid:107)∇ ⊗ ∇ w (cid:107) pL p ( T h (Γ fix )) (7.16) h − p (cid:107) w (cid:107) pL p (Γ fix ) (cid:46) h − p (cid:107) w (cid:107) pL p ( T h (Γ fix )) + (cid:107)∇ w (cid:107) pL p ( T h (Γ fix )) (7.17)and the face stabilization term h ||| w ||| p F h (cid:46) (cid:107)∇ w (cid:107) p T h ( F h ) + h p (cid:107)∇ ⊗ ∇ w (cid:107) p T h ( F h ) (7.18)We thus conclude that ||| w ||| ph,p (cid:46) (cid:107)∇ w (cid:107) p T h + h − p (cid:107) w (cid:107) pL p ( T h (Γ fix )) + (cid:107)∇ w (cid:107) pL p ( T h (Γ fix )) + h p (cid:107)∇ ⊗ ∇ w (cid:107) pL p ( T h (Γ fix )) (7.19)+ (cid:107)∇ w (cid:107) p T h ( F h ) + h p (cid:107)∇ ⊗ ∇ w (cid:107) p T h ( F h ) Setting w = v − π h v and using the interpolation error estimate (7.12) and the identity h p (cid:107)∇ ⊗ ∇ ( v − π h v ) (cid:107) pL p ( T ) = h p (cid:107)∇ ⊗ ∇ v (cid:107) pL p ( T ) , which holds since we consider piecewiselinear elements, we conclude that ||| v − π h v ||| h,p (cid:46) h (cid:107) v (cid:107) W p (Ω) (7.20)14 .3 Error Estimates for the Primal and Dual Solutions Lemma 7.1.
The finite element approximation u h defined by (2.22) of the solution u tothe primal problem (2.8) satisfies the a priori error estimate h − (cid:107) u − u h (cid:107) L p (Ω) + ||| u − u h ||| p,h ≤ h d (1 /p − / (cid:107) u (cid:107) W p (Ω) (7.21) for ≤ p < ∞ .Proof. Using the triangle inequality we obtain ||| u − u h ||| p,h ≤ ||| u − π h u ||| p,h + ||| π h u − u h ||| p,h (7.22) (cid:46) h (cid:107) u (cid:107) W p (Ω) + ||| π h u − u h ||| p,h (7.23)where we employed the energy norm interpolation estimate (7.13). For the second term onthe right hand side of (7.23) we employ the inverse inequality (7.3) with v = π h u − u h , ||| π h u − u h ||| p,h (cid:46) h d (1 /p − / ||| π h u − u h ||| ,h (7.24) (cid:46) h d (1 /p − / ( ||| u − π h u ||| ,h + ||| u − u h ||| ,h ) (7.25) (cid:46) h d (1 /p − / h (cid:107) u (cid:107) H (Ω) (7.26) (cid:46) h d (1 /p − / h (cid:107) u (cid:107) W p (Ω) (7.27)where in (7.25) we added and subtracted u and used the triangle inequality, in (7.26) weused the interpolation error estimate (7.14) with p = 2 together with the standard errorestimate ||| u − u h ||| ,h (cid:46) h (cid:107) u (cid:107) H (Ω) (7.28)see [6], and in (7.27) we used the fact that p > h − (cid:107) u − u h (cid:107) L p (Ω) using a standard duality argument. Let φ ∈ V (Ω)be the solution to the dual problem a (Ω; v, φ ) = ( ψ, v ) ∀ v ∈ V (Ω) (7.29)with ψ ∈ L q (Ω) and 1 /p +1 /q = 1. Then we have the elliptic regularity estimate (cid:107) φ (cid:107) W q (Ω) (cid:46) (cid:107) ψ (cid:107) L q (Ω) , and using concistency we conclude that A h (Ω; v, φ ) = ( ψ, v ) Ω ∀ v ∈ V h (Ω) + V (Ω) (7.30)Setting ψ = ( u − u h ) | u − u h | p − and v = u − u h we obtain (cid:107) u − u h (cid:107) pL p (Ω) = a ( u − u h , φ ) (7.31)= A h (Ω; u − u h , φ ) (7.32)= A h (Ω; u − u h , φ − π h φ ) (7.33) ≤ ||| u − u h ||| p,h ||| φ − π h φ ||| q,h (7.34) (cid:46) h ||| u − u h ||| p,h (cid:107) φ (cid:107) W q (Ω) (7.35) (cid:46) h ||| u − u h ||| p,h (cid:107) u − u h (cid:107) p/qL p (Ω) (7.36)15here we used the identity (cid:107) ψ (cid:107) L q (Ω) = (cid:107) u − u h (cid:107) p/qL q (Ω) , and thus we conclude that (cid:107) u − u h (cid:107) L p (Ω) (cid:46) h ||| u − u h ||| p,h (7.37)since p − p/q = 1. Lemma 7.2.
The finite element approximation p h defined by (3.21) of the solution p tothe dual problem (3.14) satisfies the a priori error estimate h − (cid:107) p − p h (cid:107) L p (Ω) + ||| p − p h ||| p,h (cid:46) h d (1 /p − / (cid:16) (cid:107) u (cid:107) W p (Ω) + (cid:107) p (cid:107) W p (Ω) (cid:17) (7.38) for ≤ p < ∞ .Proof. We proceed as in the proof of Lemma 7.1, with the difference that we need toaccount for the error in the right hand side. We obtain ||| π h p − p h ||| p,h (cid:46) h d (1 /p − / ||| π h p − p h ||| ,h (7.39) (cid:46) h d (1 /p − / A h (Ω; π h p − p h , π h p − p h ) (7.40) (cid:46) h d (1 /p − / ( A h (Ω; π h p − p, π h p − p h ) + a (Ω , p − p h , π h p − p h )) (7.41) (cid:46) h d (1 /p − / (cid:16) ||| π h p − p ||| ,h ||| π h p − p h ||| ,h (7.42)+ | m ( π h p − p h ) − m h ( π h p − p h ) | (cid:17) (cid:46) h d (1 /p − / (cid:16) h | p | H (Ω) ||| π h p − p h ||| ,h (7.43)+ h (cid:107) u (cid:107) H (Ω) ||| π h p − p h ||| ,h (cid:17) (cid:46) h d (1 /p − / (cid:16) | p | H (Ω) + | u | H (Ω) (cid:17) ||| π h p − p h ||| p,h (7.44)where we used a trace inequality and (7.21) to conclude that m ( v ) − m h ( v ) = ( u − u h , v ) Γ = (cid:107) u − u h (cid:107) H (Ω) (cid:107) v (cid:107) H (Ω) (cid:46) h (cid:107) u (cid:107) H (Ω) ||| v ||| ,h (7.45)To bound (cid:107) p − p h (cid:107) L p (Ω) we use a duality argument as in the proof of Lemma 7.1. Remark 7.1.
Lemma 7.1 and 7.2 are suboptimal for p >
2. Numerical test shows thatthe optimal error estimates, obtained by setting d = 0 in the bounds (7.21) and (7.45),hold for sufficiently smooth u and p . Remark 7.2.
In the analysis we have for simplicity assumed that the boundary is exact.The discrete approximation of the boundary may, however, be taken into account in theanalysis using the techniques in [7], under the assumption that the piecewise linear level setrepresentation of the boundary is second order accurate and that the associated discretenormal is first order accurate. Such an analysis shows that the geometric error is of order O ( h ) and thus of optimal order. 16 .4 Error Estimate for the Velocity Field Theorem 7.3.
Let d ≤ , β be the solution to (4.3), and β h be the solution to (4.17), then (cid:107) β − β h (cid:107) H (Ω ) ≤ M / h / (7.46) where M = (cid:107) β (cid:107) H (Ω \ Γ) + (cid:107) u (cid:107) W (Ω) + (cid:107) p (cid:107) W (Ω) + (cid:107) f (cid:107) L (Ω) + (cid:107) g N (cid:107) L (Γ) (7.47) Proof.
Adding and subtracting a Scott-Zhang interpolant π h β and using the weak formu-lations (4.3) and (4.17) we obtain (cid:107) β − β h (cid:107) H (Ω ) = ( β − β h , β − π h β ) H (Ω ) + ( β − β h , π h β − β h ) H (Ω ) (7.48)= ( β − β h , β − π h β ) H (Ω ) + D Ω ,e h L (Ω , u, p ) − D Ω ,e h L (Ω , u h , p h ) (7.49)where e h = π h β − β h . Estimating the right hand side we arrive at the bound (cid:107) β − β h (cid:107) H (Ω ) (cid:46) (cid:107) β − π h β (cid:107) H (Ω ) (cid:124) (cid:123)(cid:122) (cid:125) I + | D Ω ,e h L (Ω , u, p ) − D Ω ,e h L (Ω , u h , p h ) (cid:124) (cid:123)(cid:122) (cid:125) II | (7.50)Here Term I is an interpolation error term which needs special treatment due to thelimited regularity (4.16) of β across the interface Γ and Term II accounts for the error inthe velocity field that emanates from the approximation of the primal and dual solutionsin the discrete problem (4.17). Term I . Let T h, (Γ) be the set of all elements T ∈ T h, such that N ( T ) ∩ Γ (cid:54) = ∅ , where N ( T ) is the set of all elements that are neighbors to T . Then we have the estimates (cid:107) β − π h β (cid:107) H ( T ) (cid:46) h (cid:107) β (cid:107) H ( N ( T )) T ∈ T h \ T h (Γ) (7.51)and (cid:107) β − π h β (cid:107) H ( T ) (cid:46) (cid:107) β (cid:107) H ( N ( T )) T ∈ T h (7.52)see [25]. Summing over all elements we obtain I = (cid:107) β − π h β (cid:107) H (Ω ) (7.53) (cid:46) (cid:88) T ∈T h, \T h, (Γ) h (cid:107) β (cid:107) H ( N ( T )) + (cid:88) T ∈T h, \T h, (Γ) (cid:107) β (cid:107) H ( N ( T )) (7.54) (cid:46) h (cid:16) (cid:107) β (cid:107) H (Ω ) + (cid:107) β (cid:107) H (Ω ) (cid:124) (cid:123)(cid:122) (cid:125) = (cid:107) β (cid:107) H \ Γ) =: M (cid:17) + (cid:107) β (cid:107) H ( U δ (Γ)) (cid:124) (cid:123)(cid:122) (cid:125) (cid:70) (7.55)where Ω = φ − (( −∞ , = φ − ([0 , ∞ )), and U δ (Γ) = ∪ x ∈ Γ B δ ( x ) is the tubularneighborhood of Γ of thickness δ . 17bserving that δ ∼ h we may estimate (cid:70) by taking the L ∞ norm in the directionorthogonal to Γ, in the following way (cid:70) = (cid:107) β (cid:107) H ( U δ (Γ)) (cid:46) h sup t ∈ ( − δ,δ ) (cid:107) β (cid:107) H (Γ t ) (7.56)where Γ t = φ − ( t ). Next, we note that defining the domainsΩ ,t = φ − (( −∞ , t ]) , Ω ,t = φ − ([ t, ∞ )) (7.57)we have Ω ,t ⊆ Ω , = Ω t ∈ ( − δ, , Ω ,t ⊆ Ω , = Ω t ∈ [0 , δ ) (7.58)Therefore we have the trace inequalities (cid:107) v (cid:107) H (Γ t ) (cid:46) (cid:107) v (cid:107) H (Ω ,t ) (cid:46) (cid:107) v (cid:107) H (Ω ) t ∈ ( − δ, v ∈ H (Ω ) (7.59) (cid:107) v (cid:107) H (Γ t ) (cid:46) (cid:107) v (cid:107) H (Ω ,t ) (cid:46) (cid:107) v (cid:107) H (Ω ) t ∈ [0 , δ ) v ∈ H (Ω ) (7.60)which we may use to conclude thatsup t ∈ ( − δ, (cid:107) β (cid:107) H (Γ t ) (cid:46) (cid:107) β (cid:107) H (Ω ) , sup t ∈ [0 ,δ ) (cid:107) β (cid:107) H (Γ t ) (cid:46) (cid:107) β (cid:107) H (Ω ) (7.61)Thus we obtain the estimatesup t ∈ ( − δ,δ ) (cid:107) β (cid:107) H (Γ t ) (cid:46) (cid:107) β (cid:107) H (Ω ) + (cid:107) β (cid:107) H (Ω ) = M (7.62)which together with (7.56) gives (cid:70) (cid:46) M h (7.63)Finally, combining estimates (7.55) and (7.63) we obtain I = (cid:107) β − π h β (cid:107) H (Ω ) (cid:46) M ( h + h ) (cid:46) M h (7.64)for 0 < h ≤ h . Similar bounds are used in [7]. Term II . Using the notation e h = π h β − β h and B = De h + ( De h ) T , we decomposeTerm II as follows II = D Ω ,e h L (Ω , u, p ) − D Ω ,e h L (Ω , u h , p h ) (7.65)= (cid:90) Ω ( B ∇ u · ∇ p − B ∇ u h · ∇ p h ) dΩ (7.66)+ (cid:90) Ω ( ∇ · e h )( ∇ u · ∇ p − ∇ u h · ∇ p h ) dΩ+ (cid:90) Ω ( ∇ · e h ) f ( p − p h ) dΩ+ (cid:90) Γ ( ∇ Γ · e h ) 12 ( u − u h ) dΓ+ (cid:90) Γ ( ∇ Γ · e h ) g N ( p − p h ) dΓ= II + II + II + II + II (7.67)18tilizing the a priori error estimates (7.21) and (7.38) for the primal and dual problemswe obtain the following bounds II (cid:46) δ (cid:107) e h (cid:107) H (Ω ) + δ − ( h + h − d ) (cid:16) (cid:107) u (cid:107) W (Ω) + (cid:107) p (cid:107) W (Ω) (cid:17) (7.68) II (cid:46) δ (cid:107) e h (cid:107) H (Ω ) + δ − ( h + h − d ) (cid:16) (cid:107) u (cid:107) W (Ω) + (cid:107) p (cid:107) W (Ω) (cid:17) (7.69) II (cid:46) δ (cid:107) e h (cid:107) H (Ω ) + δ − h (8 − d ) / (cid:16) (cid:107) f (cid:107) L (Ω) + (cid:107) p (cid:107) W (Ω) (cid:17) (7.70) II (cid:46) δ (cid:107) e h (cid:107) H (Ω ) + δ − h (5 − d ) / (cid:107) u (cid:107) W (Ω) (7.71) II (cid:46) δ (cid:107) e h (cid:107) H (Ω ) + δ − h (5 − d ) / (cid:16) (cid:107) g N (cid:107) L (Γ) + (cid:107) u (cid:107) W (Ω) (cid:17) (7.72)Detailed derivations of estimates (7.68) -(7.72) are included in Appendix A. Collecting theestimates (7.68)-(7.72), using the fact d ≤
3, and defining M = (cid:107) u (cid:107) W (Ω) + (cid:107) p (cid:107) W (Ω) + (cid:107) f (cid:107) L (Ω) + (cid:107) g N (cid:107) L (Γ) (7.73)we obtain the bound II (cid:46) δ (cid:107) e h (cid:107) H (Ω ) + δ − M h (7.74) (cid:46) δ (cid:107) β − β h (cid:107) H (Ω ) + δ (cid:107) β − π h β (cid:107) H (Ω ) (cid:124) (cid:123)(cid:122) (cid:125) = I (cid:46) M h (7.64) + δ − M h (7.75)where we added and subtracted β in the first term. Conclusion of the Proof.
Starting from (7.50) and using the estimates (7.64) and(7.75) of I and II we obtain (cid:107) β − β h (cid:107) H (Ω ) (cid:46) δ (cid:107) β − β h (cid:107) H (Ω ) + (1 + δ ) M h + M δ − h (7.76)and thus, taking δ = 1 /
2, we obtain (cid:107) β − β h (cid:107) H (Ω ) (cid:46) ( M + M ) h = M h (7.77)where M = M + M , which completes the proof. We use the following settings in the numerical examples • Optimization algorithm – α = 0 .
5: Damping parameter in (5.5)19 N = 3: Number of time steps in (5.6) – TOL = 10 − : Tolerance in (6.1) • Finite element methods – γ = 1: Penalty parameter for the gradient jump (2.20) – γ D = 10: Penalty parameter for the Dirichlet boundary condition (2.21) – γ = 1: Penalty parameter for the gradient jump (5.7) Model Problem 1.
To define the domain Ω we let Ω = [0 , be the unit square andΩ ⊂ Ω be a domain in the interior of Ω with boundary Γ, finally, let Ω = Ω \ Ω . Wenote that ∂ Ω = Γ ∪ ∂ Ω and that Γ ∩ ∂ Ω = ∅ .With this set up we consider a Bernoulli free boundary value problem where the exactposition of the free boundary Γ is a circle of radius r = 0 .
25 centered in (0 . , .
5) and theexact solution is u ref = 4(( x − . + ( y − . ) / −
1. The corresponding Bernoulli freeboundary problem takes the form − ∆ u = − ∆ u ref in Ω (8.1) u = u ref on ∂ Ω (8.2) n · ∇ u = − u = 0 on Γ (8.4)We will use a level set function corresponding to the domain displayed in Figure 1 (rightsub-figure) as an initial guess.Figure 1: The final domain (left) and initial guess (right) for Model Problem 1. The grayarea is the computational domain Ω, the outer square boundary is fixed, and Γ is the freeboundary. 20 odel Problem 2. Let Ω = [0 , be the unit square as before and Ω ⊂ Ω be asubdomain in the interior of Ω with boundary Γ. Next let Ω and Ω be the balls of radius R = 1 /
12 centered in the points (1 / , /
3) and (2 / , / \ (Ω ∪ Ω )and consider the boundary conditions u = 1 on ∂ Ω \ Γ (8.5) n · ∇ u = − u = 0 on Γ (8.7)In this example there is no known exact position of the free boundary. We will use a levelset function corresponding to the domain Figure 2 (right sub-figure) as an initial guess.Figure 2: The final domain (left) and initial guess (right) for Model Problem 2. The grayarea is the computational domain Ω, the outer square boundary and the two inner mostcircles are fixed, and Γ is the free bondary. We investigate the convergence rate of the discrete velocity field for Model Problem 1. InFigure 3 we display the error in the discrete velocity field in the H -norm and L -normwhere the reference solution β ref is computed on a quasi-uniform mesh with 526338 degreesof freedom. We obtain slightly better convergence rates than O ( h / ) and O ( h / ) in H and L -norm, respectively, which is in agreement with Theorem 7.3. In Figure 4 and Figure 5 we present the convergence history of R Γ , see (6.1), for ModelProblem 1 and 2. In Figure 6 we show the approximation of Ω obtained after 0 , , , and 46 iterations, where iteration 46 is the final domain. In Figure 6 we note that werapidly obtain a domain which resembles the final domain, but to straighten the kinks inthe boundary takes some extra effort. 21 Degrees of freedom (dofs) -3 -2 -1 E rr o r H -errorH -errordofs dofs Figure 3: Convergence of the error in velocity field in Model Problem 1 in the H and L -norm. Iterations -6 -5 -4 -3 -2 -1 R ! Model 1
Figure 4: Convergence of R Γ for Model Problem 1.22
10 20 30 40 50
Iterations -6 -5 -4 -3 -2 -1 R ! Model 2
Figure 5: Convergence of R Γ for Model Problem 2. (a) Iteration 0. (b) Iteration 5.(c) Iteration 15. (d) Iteration 46. Figure 6: The computational domain for Model Problem 2 after 0, 5, 15, and 46 iterations.
A Bounds for
I I − I I in the Proof of Theorem 7.3 Term II . Dividing II into three suitable terms and then using H¨older’s inequality weobtain 23 I = (cid:90) Ω ( B ∇ u · ∇ p − B ∇ u h · ∇ p h ) dΩ (A.1)= ( B ∇ ( u − u h ) , ∇ p ) L (Ω) (A.2)+ ( B ∇ u, ∇ ( p − p h )) L (Ω) − ( B ∇ ( u − u h ) , ∇ ( p − p h )) L (Ω) ≤ (cid:107) B (cid:107) L (Ω) (cid:107)∇ ( u − u h ) (cid:107) L (Ω) (cid:107)∇ p (cid:107) L ∞ (Ω) (A.3)+ (cid:107) B (cid:107) L (Ω) (cid:107)∇ u (cid:107) L ∞ (Ω) (cid:107)∇ ( p − p h ) (cid:107) L (Ω) + (cid:107) B (cid:107) L (Ω) (cid:107)∇ ( u − u h ) (cid:107) L (Ω) (cid:107)∇ ( p − p h ) (cid:107) L (Ω) ≤ (cid:107) e h (cid:107) H (Ω) h (cid:107) u (cid:107) W (Ω) (cid:107) p (cid:107) W (Ω) (A.4)+ (cid:107) e h (cid:107) H (Ω) (cid:107) u (cid:107) W (Ω) h (cid:107) p (cid:107) W (Ω) + (cid:107) e h (cid:107) H (Ω) h (4 − d ) / (cid:107) u (cid:107) W (Ω) (cid:107) p (cid:107) W (Ω) (cid:46) δ (cid:107) e h (cid:107) H (Ω) + δ − (cid:16) h + h (4 − d ) (cid:17) (cid:107) u (cid:107) W (Ω) (cid:107) p (cid:107) W (Ω) (A.5)where in (A.4) we used the a priori error estimates (7.21) and (7.38), with p = 2 for thefirst two terms and with p = 4 for the third term, and the Sobolev embedding theorem toconclude that (cid:107)∇ v (cid:107) L ∞ (Ω) (cid:46) (cid:107) v (cid:107) W (Ω) (A.6)for v = p and v = u since d ≤
3, and finally in (A.5) we used the basic bound (cid:107) v (cid:107) W (Ω) (cid:46) (cid:107) v (cid:107) W (Ω) for v = p and v = u . Term II . Using the same approach as for Term II (with B replaced by ∇ · e h ) weobtain II (cid:46) δ (cid:107) e h (cid:107) H (Ω) + δ − h (cid:107) u (cid:107) W (Ω) (cid:107) p (cid:107) W (Ω) (A.7) Term II . Using H¨older’s inequality II = (cid:90) Ω ( ∇ · e h ) f ( p − p h ) dΩ (A.8) ≤ (cid:107)∇ · e h (cid:107) L (Ω ) (cid:107) f (cid:107) L (Ω) (cid:107) p − p h (cid:107) L (Ω) (A.9) (cid:46) (cid:107) e h (cid:107) H (Ω ) (cid:107) f (cid:107) L (Ω) h (8 − d ) / (cid:107) p (cid:107) W (Ω) (A.10) (cid:46) δ (cid:107) e h (cid:107) H (Ω ) + δ − h (8 − d ) / (cid:107) f (cid:107) L (Ω) (cid:107) p (cid:107) W (Ω) (A.11)24 erm II . Using the conjugate rule followed by H¨older’s inequality, (cid:90) Γ ( ∇ Γ · e h )( u − u h )) dΓ = (cid:90) Γ ( ∇ Γ · e h )( u + u h )( u − u h )) dΓ (A.12) ≤ (cid:107)∇ Γ · e h (cid:107) L (Γ) (cid:107) u + u h (cid:107) L (Γ) (cid:107) u − u h (cid:107) L (Γ) (A.13) ≤ δh (cid:107)∇ Γ · e h (cid:107) L (Γ) + δ − h − (cid:107) u + u h (cid:107) L (Γ) (cid:107) u − u h (cid:107) L (Γ) (A.14) ≤ δ (cid:107) e h (cid:107) H (Ω ) + δ − h (5 − d ) / (cid:107) u (cid:107) W (Ω) (A.15)Here we used the trace inequality (cid:107) v (cid:107) L p (Γ) ≤ (cid:107) v (cid:107) L p ( ∂ Ω) (cid:46) (cid:107) v (cid:107) − /pL p (Ω) (cid:107) v (cid:107) /pW p (Ω) v ∈ W p (Ω) (A.16)with p = 4 and v = u − u h followed by the a priori error estimate (7.21) to conclude that (cid:107) u − u h (cid:107) L (Γ) (cid:46) (cid:107) u − u h (cid:107) / L (Ω) (cid:107) u − u h (cid:107) / W (Ω) (A.17) (cid:46) ( h (8 − d ) / ) / ( h (4 − d ) / ) / (cid:107) u (cid:107) W (Ω) (A.18) (cid:46) h (7 − d ) / (cid:107) u (cid:107) W (Ω) (A.19)and the following estimate (cid:107) u + u h (cid:107) L (Γ) (cid:46) (cid:107) u (cid:107) L (Γ) + (cid:107) u − u h (cid:107) L (Γ) (A.20) (cid:46) (cid:107) u (cid:107) L (Γ) + (cid:107) u − u h (cid:107) W (Ω) (A.21) (cid:46) (cid:107) u (cid:107) L (Γ) + h (4 − d ) / (cid:107) u (cid:107) W (Ω) (A.22) (cid:46) (cid:107) u (cid:107) W (Ω) (A.23)which holds since h ∈ (0 , h ]. Term II . Using H¨older’s inequality II = (cid:90) Γ ( ∇ Γ · e h ) g N ( p − p h ) dΓ (A.24) ≤ (cid:107)∇ Γ · e h (cid:107) L (Γ) (cid:107) g N (cid:107) L (Γ) (cid:107) p − p h (cid:107) L (Γ) (A.25) ≤ δh (cid:107)∇ Γ · e h (cid:107) L (Γ) + δh − (cid:107) g N (cid:107) L (Γ) (cid:107) p − p h (cid:107) L (Γ) (A.26) ≤ δ (cid:107) e h (cid:107) H (Ω ) + h (5 − d ) / (cid:107) g N (cid:107) L (Γ) (cid:107) p (cid:107) W (Ω) (A.27)where in (A.27) the first term was estimated using the inverse estimate h (cid:107) v (cid:107) L (Γ ∩ T ) (cid:46) (cid:107) v (cid:107) H ( T ) v ∈ P ( T ) (A.28)for T ∈ T h such that Γ ∩ T (cid:54) = ∅ with v = e h , and for the second term we used the estimate (cid:107) p − p h (cid:107) L (Γ) (cid:46) h (7 − d ) / (cid:107) p (cid:107) W (Ω) (A.29)which follows in the same way as in (A.17-A.19).25 eferences [1] G. Allaire, C. Dapogny, and P. Frey. Shape optimization with a level set based meshevolution method. Comput. Methods Appl. Mech. Engrg. , 282:22–53, 2014.[2] G. Allaire, F. Jouve, and A.-M. Toader. Structural optimization by the level-setmethod. In
Free boundary problems (Trento, 2002) , volume 147 of
Internat. Ser.Numer. Math. , pages 1–15. Birkh¨auser, Basel, 2004.[3] A. Beurling.
On free boundary problems for the Laplace equation . Seminars on analyticfunctions I. Institute Advanced Studies Seminars, Princeton, 1957.[4] E. Burman, S. Claus, P. Hansbo, M. G. Larson, and A. Massing. CutFEM: discretizinggeometry and partial differential equations.
Internat. J. Numer. Methods Engrg. ,104(7):472–501, 2015.[5] E. Burman and M. A. Fern´andez. Finite element methods with symmetric stabilizationfor the transient convection-diffusion-reaction equation.
Comput. Methods Appl. Mech.Engrg. , 198(33-36):2508–2519, 2009.[6] E. Burman and P. Hansbo. Fictitious domain finite element methods using cut ele-ments: II. A stabilized Nitsche method.
Appl. Numer. Math. , 62(4):328–341, 2012.[7] E. Burman, P. Hansbo, M. G. Larson, and S. Zahedi. Cut finite element methods forcoupled bulk–surface problems.
Numer. Math. , 133(2):203–231, 2016.[8] M. Cenanovic, P. Hansbo, and M. G. Larson. Minimal surface computation using afinite element method on an embedded surface.
Internat. J. Numer. Methods Engrg. ,104(7):502–512, 2015.[9] Z. Chen and J. Zou. Finite element methods and their convergence for elliptic andparabolic interface problems.
Numer. Math. , 79(2):175–202, 1998.[10] R. Correa and A. Seeger. Directional derivates in minimax problems.
Numer. Funct.Anal. Optim. , 7(2-3):145–156, 1984/85.[11] C. Dapogny and P. Frey. Computation of the signed distance function to a discretecontour on adapted triangulation.
Calcolo , 49(3):193–219, 2012.[12] F. de Gournay. Velocity extension for the level-set method and multiple eigenvaluesin shape optimization.
SIAM J. Control Optim. , 45(1):343–367, 2006.[13] K. Eppler, H. Harbrecht, and M. S. Mommer. A new fictitious domain method inshape optimization.
Comput. Optim. Appl. , 40(2):281–298, 2008.[14] M. Flucher and M. Rumpf. Bernoulli’s free-boundary problem, qualitative theory andnumerical approximation.
J. Reine Angew. Math. , 486:165–204, 1997.2615] P. Hansbo. Nitsche’s method for interface problems in computational mechanics.
GAMM-Mitt. , 28(2):183–206, 2005.[16] P. Hansbo, M. G. Larson, and S. Zahedi. A cut finite element method for coupledbulk-surface problems on time-dependent domains.
Comput. Methods Appl. Mech.Engrg. , 307:96–116, 2016.[17] H. Harbrecht. Analytical and numerical methods in shape optimization.
Math. Meth-ods Appl. Sci. , 31(18):2095–2114, 2008.[18] H. Harbrecht. A Newton method for Bernoulli’s free boundary problem in threedimensions.
Computing , 82(1):11–30, 2008.[19] R. Hiptmair and A. Paganini. Shape optimization by pursuing diffeomorphisms.
Com-put. Methods Appl. Math. , 15(3):291–305, 2015.[20] R. Hiptmair, A. Paganini, and S. Sargheini. Comparison of approximate shape gradi-ents.
BIT , 55(2):459–485, 2015.[21] T. Y. Hou. Numerical solutions to free boundary problems.
Acta Numer. , 4:335–415,1995.[22] C. M. Kuster, P. A. Gremaud, and R. Touzani. Fast numerical methods for Bernoullifree boundary problems.
SIAM J. Sci. Comput. , 29(2):622–634, 2007.[23] A. Laurain and K. Sturm. Distributed shape derivative via averaged adjoint methodand applications.
ESAIM: Math. Model. Numer. Anal. , 50(4):1241–1267, 2016.[24] S. Osher and J. A. Sethian. Fronts propagating with curvature-dependent speed:algorithms based on Hamilton-Jacobi formulations.
J. Comput. Phys. , 79(1):12–49,1988.[25] L. R. Scott and S. Zhang. Finite element interpolation of nonsmooth functions satis-fying boundary conditions.
Math. Comp. , 54(190):483–493, 1990.[26] E. Shargorodsky and J. Toland. Bernoulli free-boundary problems.
Mem. Am. Math.Soc. , 196(914), 2008.[27] J. Sokolowski and J.-P. Zolesio.