This page was generated from unit-7-optimization/02_Shape_Derivative_Laplace.ipynb.

7.2 PDE-Constrained Shape Optimization

We want to solve the PDE-constrained shape optimization problem

\[\underset{\Omega\subset \mathsf{D}}{\mbox{min}} \; J(u) := \int_\Omega |u-u_d|^2 \; dx\]

subject to that \((\Omega,u)\) satisfy

\[\int_\Omega \nabla u \cdot \nabla v \; dx = \int_\Omega f v \; dx \; \quad \text{ for all } v \in H_0^1(\Omega),\]

where \(\Omega \subset \mathbb R^2\) for given \(u_d, f \in C^1(\mathbb R^2)\).

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import SplineGeometry

geo = SplineGeometry()
geo.AddRectangle( (0.2, 0.2), (0.8, 0.8), bcs = ('b','r','t','l'))
mesh = Mesh(geo.GenerateMesh(maxh = 0.1))

#given data of our problem (chosen such that \Omega^* = [0,1]^2 is the optimal shape)
f = CoefficientFunction(2*y*(1-y)+2*x*(1-x))
ud = x*(1-x)*y*(1-y)

grad_f = CoefficientFunction( (f.Diff(x), f.Diff(y) ) )
grad_ud = CoefficientFunction( (ud.Diff(x), ud.Diff(y) ) )

State equation

fes = H1(mesh, order=2, dirichlet=".*")
u, v = fes.TnT()
gfu = GridFunction(fes)
scene = Draw (gfu, mesh, "state")

a = BilinearForm(fes, symmetric=True)
a += grad(u)*grad(v)*dx

fstate = LinearForm(fes)
fstate += f*v*dx

def SolveStateEquation():
    rhs = gfu.vec.CreateVector()
    rhs.data = fstate.vec - a.mat * gfu.vec
    update = gfu.vec.CreateVector()
    update.data = a.mat.Inverse(fes.FreeDofs()) * rhs
    gfu.vec.data += update

a.Assemble()
fstate.Assemble()
SolveStateEquation()
scene.Redraw()

Adjoint equation

We set up the adjoint equation

\[\mbox{Find } p \in H_0^1(\Omega): \int_\Omega \nabla w \cdot \nabla p \, \mbox dx = - \partial_u J(u)(w) \quad \text{ for all } w \in H_0^1(\Omega)\]

where \(u\) is the solution to the state equation. For \(J(u) = \int_\Omega |u-u_d|^2 \mbox dx\), we get

\[\partial_u J(u)(w) = -2 \int_\Omega (u-u_d)w \,\mbox dx.\]

However, we can also use the Diff(…) command:

def Cost(u):
    return (u-ud)**2*dx

p, w = fes.TnT()
gfp = GridFunction(fes)
scene = Draw (gfp, mesh, "adjoint")

fadjoint = LinearForm(fes)
fadjoint += -1*Cost(gfu).Diff(gfu,w)

def SolveAdjointEquation():
    rhs = gfp.vec.CreateVector()
    rhs.data = fadjoint.vec - a.mat.T * gfp.vec
    update = gfp.vec.CreateVector()
    update.data = a.mat.Inverse(fes.FreeDofs()).T * rhs
    gfp.vec.data += update

fadjoint.Assemble()
SolveAdjointEquation()
scene.Redraw()

Note that (for linear problems) the operator on the left hand side of the adjoint equation is just the transpose of the state operator.

Shape derivative

The shape derivative of the shape function

\[\mathcal{J}(\Omega):= \int_\Omega |u-u_d|^2\;dx\]

at \(\Omega\) in direction \(X\) is given by the formula

\[\begin{split}\begin{array}{rl} d\mathcal J(\Omega; X) =&\int_{\Omega} \mbox{div}(X) |u-u_d|^2 - 2(u-u_d)\nabla u_d \cdot X \, dx \\ & + \int_{\Omega} (\mbox{div}(X) I - D X - D X^\top )\nabla u \cdot \nabla p \, dx \\ &- \int_\Omega (\nabla f\cdot X + f \mbox{div}(X)) p\, dx. \end{array}\end{split}\]

We now assemble this shape derivaitve in NGSolve as follows:

VEC = H1(mesh, order=2, dim=2)
PHI, X = VEC.TnT()

dJOmega = LinearForm(VEC)
dJOmega += SymbolicLFI(div(X)*( (gfu-ud)**2 - f*gfp + InnerProduct(grad(gfu),grad(gfp))))
dJOmega += SymbolicLFI(-2*(gfu-ud)*InnerProduct(grad_ud,X) - InnerProduct(grad_f,X)*gfp)
dJOmega += SymbolicLFI(-InnerProduct(X.Deriv()*grad(gfu),grad(gfp))-InnerProduct(grad(gfu),X.Deriv()*grad(gfp)))

Descent Direction

• Next, we want to find a vector field \(X\) which yields a decrease of the objective function, i.e. we want to find a vector field \(X\) such that

  \[d \mathcal J(\Omega;X)<0.\]

• This can be achieved by solving an auxiliary boundary value problem of the type

  \[\mbox{Find } X \in H: \qquad b(X, \Phi) = d\mathcal J(\Omega; \Phi) \; \text{ for all } \Phi \in H\]

  for a suitable Hilbert space \(H\) (e.g. \(H=H^1(\Omega)^2\)). Here, \(b(\cdot, \cdot): H \times H \rightarrow \mathbb R\) is a positive definite bilinear form which we are free to choose. Then, \(-X\) is a descent direction since

  \[d\mathcal J(\Omega; -X) = -d\mathcal J(\Omega; X) = - b(X, X) < 0.\]

b = BilinearForm(VEC)
b += InnerProduct(grad(X),grad(PHI))*dx + InnerProduct(X,PHI)*dx

gfX = GridFunction(VEC)

# gfset denotes the deformation of the original domain and will be updated during the shape optimization
gfset = GridFunction(VEC)
gfset.Set((0,0))
scene = Draw(gfset,mesh,"gfset")
SetVisualization (deformation=True)

def SolveDeformationEquation():
    rhs = gfX.vec.CreateVector()
    rhs.data = dJOmega.vec - b.mat * gfX.vec
    update = gfX.vec.CreateVector()
    update.data = b.mat.Inverse(VEC.FreeDofs()) * rhs
    gfX.vec.data += update

b.Assemble()
dJOmega.Assemble()
SolveDeformationEquation()
Draw(-gfX, mesh, "gfX")

WebGLScene

Now let us update the geometry in the direction of \(-X\). Here, we need to choose the step size accordingly (a descent is only guaranteed for this step size being sufficiently small)

print('Cost at initial design', Integrate (Cost(gfu), mesh))
gfset.Set((0,0))
scale = 0.5 / Norm(gfX.vec)
gfset.vec.data -= scale * gfX.vec

Cost at initial design 0.00048657855734481254

Let us now update the design and compute the new value of the cost function:

scene = Draw(gfset)
mesh.SetDeformation(gfset)
scene.Redraw()

a.Assemble()
fstate.Assemble()
SolveStateEquation()
print('Cost at new design', Integrate (Cost(gfu), mesh))

Cost at new design 0.00032579696176024647

#reset the design
gfset.Set((0,0))
mesh.SetDeformation(gfset)

#check if initial value of cost function 0.000486578350214552 is recovered
a.Assemble()
fstate.Assemble()
SolveStateEquation()
print('Cost at new design', Integrate (Cost(gfu), mesh))

Cost at new design 0.0004865785573448126

Finally, this can be applied in an iterative algorithm with a line search for choosing the step size:

scene = Draw(gfset)

gfset.Set((0,0))
mesh.SetDeformation(gfset)
scene.Redraw()
a.Assemble()
fstate.Assemble()
SolveStateEquation()

LineSearch = True

iter_max = 600
Jold = Integrate(Cost(gfu), mesh)
converged = False

# input("Press enter to start optimization")
for k in range(iter_max):
    mesh.SetDeformation(gfset)
    scene.Redraw()

    print('cost at iteration', k, ': ', Jold)

    a.Assemble()
    fstate.Assemble()
    SolveStateEquation()

    fadjoint.Assemble()
    SolveAdjointEquation()

    b.Assemble()
    dJOmega.Assemble()
    SolveDeformationEquation()

    scale = 0.01 / Norm(gfX.vec)
    gfsetOld = gfset
    gfset.vec.data -= scale * gfX.vec

    Jnew = Integrate(Cost(gfu), mesh)

    if LineSearch:
        while Jnew > Jold and scale > 1e-12:
            #input('a')
            scale = scale / 2
            print("scale = ", scale)
            if scale <= 1e-12:
                converged = True
                break

            gfset.vec.data = gfsetOld.vec - scale * gfX.vec
            mesh.SetDeformation(gfset)

            a.Assemble()
            fstate.Assemble()
            SolveStateEquation()
            Jnew = Integrate(Cost(gfu), mesh)

    if converged==True:
        print("No more descent can be found")
        break
    Jold = Jnew

    Redraw(blocking=True)

cost at iteration 0 :  0.0004865785573448126
scale =  5.0172940372132
cost at iteration 1 :  0.0004834070682129935
scale =  4.76812536166299
cost at iteration 2 :  0.00048011519350089313
scale =  4.534231415136696
cost at iteration 3 :  0.0004766971621849114
scale =  4.316155256294001
cost at iteration 4 :  0.0004731480288140851
scale =  4.113771843070691
cost at iteration 5 :  0.00046946364120922617
scale =  3.9265249630949883
cost at iteration 6 :  0.00046564059160264314
scale =  3.7536033150233488
cost at iteration 7 :  0.000461676160252617
scale =  3.5940646155631444
cost at iteration 8 :  0.0004575682575672066
scale =  3.446919163814335
cost at iteration 9 :  0.00045331536842711914
scale =  3.3111835819313984
cost at iteration 10 :  0.0004489165007265815
scale =  3.1859135216674597
cost at iteration 11 :  0.00044437113904143
scale =  3.070222024991224
cost at iteration 12 :  0.00043967920364183376
scale =  2.9632883913475627
cost at iteration 13 :  0.0004348410146650856
scale =  2.8643609572773254
cost at iteration 14 :  0.00042985726105236724
scale =  2.7727561189579695
cost at iteration 15 :  0.00042472897376122105
scale =  2.6878551588136355
cost at iteration 16 :  0.0004194575027445427
cost at iteration 17 :  0.0004193946628183917
cost at iteration 18 :  0.000415744959939901
cost at iteration 19 :  0.00041203299753923877
cost at iteration 20 :  0.0004082594338459174
cost at iteration 21 :  0.0004044249768264305
cost at iteration 22 :  0.0004005303828739452
cost at iteration 23 :  0.00039657645561609736
cost at iteration 24 :  0.00039256404482998654
cost at iteration 25 :  0.00038849404545433243
cost at iteration 26 :  0.0003843673966895961
cost at iteration 27 :  0.0003801850811776647
cost at iteration 28 :  0.00037594812425343206
cost at iteration 29 :  0.0003716575932613135
cost at iteration 30 :  0.00036731459693035263
cost at iteration 31 :  0.0003629202848021809
cost at iteration 32 :  0.0003584758467066134
cost at iteration 33 :  0.0003539825122801603
cost at iteration 34 :  0.00034944155052317395
cost at iteration 35 :  0.00034485426939175415
cost at iteration 36 :  0.00034022201542090115
cost at iteration 37 :  0.00033554617337573314
cost at iteration 38 :  0.0003308281659278783
cost at iteration 39 :  0.00032606945335443
cost at iteration 40 :  0.00032127153325708154
cost at iteration 41 :  0.00031643594029928424
cost at iteration 42 :  0.000311564245959466
cost at iteration 43 :  0.00030665805829851187
cost at iteration 44 :  0.00030171902173988636
cost at iteration 45 :  0.0002967488168608968
cost at iteration 46 :  0.0002917491601937371
cost at iteration 47 :  0.0002867218040350602
cost at iteration 48 :  0.0002816685362629307
cost at iteration 49 :  0.00027659118016009325
cost at iteration 50 :  0.00027149159424258315
cost at iteration 51 :  0.0002663716720927741
cost at iteration 52 :  0.0002612333421960104
cost at iteration 53 :  0.00025607856778005216
cost at iteration 54 :  0.00025090934665658636
cost at iteration 55 :  0.0002457277110641179
cost at iteration 56 :  0.00024053572751158274
cost at iteration 57 :  0.00023533549662206604
cost at iteration 58 :  0.0002301291529760313
cost at iteration 59 :  0.0002249188649534908
cost at iteration 60 :  0.00021970683457456404
cost at iteration 61 :  0.00021449529733789382
cost at iteration 62 :  0.00020928652205638766
cost at iteration 63 :  0.00020408281068976682
cost at iteration 64 :  0.00019888649817340436
cost at iteration 65 :  0.0001936999522429261
cost at iteration 66 :  0.00018852557325404546
cost at iteration 67 :  0.00018336579399709006
cost at iteration 68 :  0.0001782230795056496
cost at iteration 69 :  0.00017309992685875982
cost at iteration 70 :  0.00016799886497600097
cost at iteration 71 :  0.00016292245440484253
cost at iteration 72 :  0.0001578732870995265
cost at iteration 73 :  0.00015285398619070944
cost at iteration 74 :  0.00014786720574502403
cost at iteration 75 :  0.0001429156305136172
cost at iteration 76 :  0.0001380019756686315
cost at iteration 77 :  0.00013312898652646733
cost at iteration 78 :  0.0001282994382565103
cost at iteration 79 :  0.00012351613557383932
cost at iteration 80 :  0.00011878191241421173
cost at iteration 81 :  0.0001140996315893717
cost at iteration 82 :  0.00010947218442042487
cost at iteration 83 :  0.00010490249034664921
cost at iteration 84 :  0.00010039349650667474
cost at iteration 85 :  9.594817728841897e-05
cost at iteration 86 :  9.156953384351481e-05
cost at iteration 87 :  8.726059356115537e-05
cost at iteration 88 :  8.302440949529279e-05
cost at iteration 89 :  7.886405973789805e-05
cost at iteration 90 :  7.478264672946267e-05
cost at iteration 91 :  7.07832964960006e-05
cost at iteration 92 :  6.686915779939642e-05
cost at iteration 93 :  6.304340118486994e-05
cost at iteration 94 :  5.930921790538869e-05
cost at iteration 95 :  5.5669818697783556e-05
cost at iteration 96 :  5.2128432378700646e-05
cost at iteration 97 :  4.868830421982279e-05
cost at iteration 98 :  4.5352694050244125e-05
cost at iteration 99 :  4.2124874018370704e-05
cost at iteration 100 :  3.900812592466795e-05
cost at iteration 101 :  3.6005738007629406e-05
cost at iteration 102 :  3.312100102499749e-05
cost at iteration 103 :  3.0357203415225813e-05
cost at iteration 104 :  2.7717625242212102e-05
cost at iteration 105 :  2.5205530506518506e-05
cost at iteration 106 :  2.2824157227799025e-05
cost at iteration 107 :  2.0576704431640238e-05
cost at iteration 108 :  1.8466314751504556e-05
cost at iteration 109 :  1.649605068217509e-05
cost at iteration 110 :  1.4668861414742735e-05
cost at iteration 111 :  1.2987535311316407e-05
cost at iteration 112 :  1.1454629801791167e-05
cost at iteration 113 :  1.0072364536423074e-05
cost at iteration 114 :  8.842452389949546e-06
cost at iteration 115 :  7.765820807104616e-06
cost at iteration 116 :  6.84213090612346e-06
cost at iteration 117 :  6.0689083805033095e-06
cost at iteration 118 :  5.4399174334077805e-06
cost at iteration 119 :  4.942172315737003e-06
cost at iteration 120 :  4.5515970652604826e-06
cost at iteration 121 :  4.2331169343429204e-06
cost at iteration 122 :  3.956671493976132e-06
cost at iteration 123 :  3.70586632598654e-06
cost at iteration 124 :  3.4699473518204077e-06
cost at iteration 125 :  3.244291371189643e-06
cost at iteration 126 :  3.027616103024207e-06
cost at iteration 127 :  2.819562199545164e-06
cost at iteration 128 :  2.6199951792908104e-06
cost at iteration 129 :  2.4288341729491163e-06
cost at iteration 130 :  2.2460123446618343e-06
cost at iteration 131 :  2.071468873280283e-06
cost at iteration 132 :  1.9051474673414052e-06
cost at iteration 133 :  1.7469959100276731e-06
cost at iteration 134 :  1.5969656359689012e-06
cost at iteration 135 :  1.4550111938767366e-06
cost at iteration 136 :  1.3210895279494625e-06
cost at iteration 137 :  1.195159029122459e-06
cost at iteration 138 :  1.077178227108377e-06
cost at iteration 139 :  9.67104031773439e-07
cost at iteration 140 :  8.648891910438971e-07
cost at iteration 141 :  7.704788943285323e-07
cost at iteration 142 :  6.838053983399858e-07
cost at iteration 143 :  6.047818021704172e-07
cost at iteration 144 :  5.332885422352927e-07
cost at iteration 145 :  4.691722786979197e-07
cost at iteration 146 :  4.122106453794152e-07
cost at iteration 147 :  3.6318560846579926e-07
cost at iteration 148 :  3.383227077800786e-07
cost at iteration 149 :  3.275058122821268e-07
cost at iteration 150 :  3.180566846079047e-07
cost at iteration 151 :  3.069752848491972e-07
cost at iteration 152 :  2.9873172207063575e-07
cost at iteration 153 :  2.8876279537055566e-07
cost at iteration 154 :  2.8148128732694684e-07
cost at iteration 155 :  2.7245365618681053e-07
cost at iteration 156 :  2.659711961069519e-07
cost at iteration 157 :  2.577458867487123e-07
cost at iteration 158 :  2.519415543854804e-07
cost at iteration 159 :  2.4440453235150207e-07
cost at iteration 160 :  2.391833939037126e-07
cost at iteration 161 :  2.322404520096468e-07
cost at iteration 162 :  2.275255125825601e-07
cost at iteration 163 :  2.2109842259744485e-07
cost at iteration 164 :  2.1682601472331654e-07
cost at iteration 165 :  2.108495235685734e-07
cost at iteration 166 :  2.069663531946836e-07
cost at iteration 167 :  2.0138579284132023e-07
cost at iteration 168 :  1.9784687750959574e-07
cost at iteration 169 :  1.92616243900553e-07
cost at iteration 170 :  1.893833818511074e-07
cost at iteration 171 :  1.8446378220846895e-07
cost at iteration 172 :  1.815043667395451e-07
cost at iteration 173 :  1.768627604170693e-07
cost at iteration 174 :  1.7414883586040645e-07
cost at iteration 175 :  1.6975701092691235e-07
cost at iteration 176 :  1.6726450751544135e-07
cost at iteration 177 :  1.6309824746926869e-07
cost at iteration 178 :  1.6080635414597253e-07
cost at iteration 179 :  1.5684475838610144e-07
cost at iteration 180 :  1.5473540498058487e-07
cost at iteration 181 :  1.5096033383564127e-07
cost at iteration 182 :  1.4901776156705669e-07
cost at iteration 183 :  1.4541338235738942e-07
cost at iteration 184 :  1.4362378652428317e-07
cost at iteration 185 :  1.401762016608254e-07
cost at iteration 186 :  1.385274337729063e-07
cost at iteration 187 :  1.3522437553039838e-07
cost at iteration 188 :  1.337056946160749e-07
cost at iteration 189 :  1.305362741436459e-07
cost at iteration 190 :  1.2913813885523297e-07
cost at iteration 191 :  1.260926393473889e-07
cost at iteration 192 :  1.2480653396964013e-07
cost at iteration 193 :  1.2187623982712219e-07
cost at iteration 194 :  1.2069452848679957e-07
cost at iteration 195 :  1.1787158383837391e-07
cost at iteration 196 :  1.1678738818445497e-07
cost at iteration 197 :  1.1406467938815515e-07
cost at iteration 198 :  1.1307177581227295e-07
cost at iteration 199 :  1.1044283356470267e-07
cost at iteration 200 :  1.0953556669380213e-07
cost at iteration 201 :  1.0699448419439799e-07
cost at iteration 202 :  1.0616769393786276e-07
cost at iteration 203 :  1.0370905821830406e-07
cost at iteration 204 :  1.0295801810969967e-07
cost at iteration 205 :  1.0057685217652089e-07
cost at iteration 206 :  9.989721713142769e-08
cost at iteration 207 :  9.758893100624682e-08
cost at iteration 208 :  9.697669293517102e-08
cost at iteration 209 :  9.473704203115299e-08
cost at iteration 210 :  9.418849201073959e-08
cost at iteration 211 :  9.201354157155047e-08
cost at iteration 212 :  9.152523749719007e-08
cost at iteration 213 :  8.94113320583662e-08
cost at iteration 214 :  8.898007088393454e-08
cost at iteration 215 :  8.692380790658121e-08
cost at iteration 216 :  8.654660172882282e-08
cost at iteration 217 :  8.454480870999915e-08
cost at iteration 218 :  8.421886408098096e-08
cost at iteration 219 :  8.226857857087433e-08
cost at iteration 220 :  8.199127852634141e-08
cost at iteration 221 :  8.008973058465914e-08
cost at iteration 222 :  7.985861896271571e-08
cost at iteration 223 :  7.80032156701013e-08
cost at iteration 224 :  7.781598336631388e-08
cost at iteration 225 :  7.600429507457169e-08
cost at iteration 226 :  7.585876793898976e-08
cost at iteration 227 :  7.408851599931699e-08
cost at iteration 228 :  7.39826441300638e-08
cost at iteration 229 :  7.225168988378816e-08
cost at iteration 230 :  7.218353811257508e-08
cost at iteration 231 :  7.048987296582322e-08
cost at iteration 232 :  7.045761236444214e-08
cost at iteration 233 :  6.879934879844687e-08
scale =  53.59239428273243
scale =  26.796197141366214
scale =  13.398098570683107
scale =  6.6990492853415535
scale =  3.3495246426707768
scale =  1.6747623213353884
scale =  0.8373811606676942
scale =  0.4186905803338471
scale =  0.20934529016692355
scale =  0.10467264508346177
scale =  0.05233632254173089
scale =  0.026168161270865443
scale =  0.013084080635432722
scale =  0.006542040317716361
scale =  0.0032710201588581804
scale =  0.0016355100794290902
scale =  0.0008177550397145451
scale =  0.00040887751985727255
scale =  0.00020443875992863628
scale =  0.00010221937996431814
scale =  5.110968998215907e-05
scale =  2.5554844991079535e-05
scale =  1.2777422495539767e-05
scale =  6.388711247769884e-06
scale =  3.194355623884942e-06
scale =  1.597177811942471e-06
scale =  7.985889059712355e-07
scale =  3.9929445298561773e-07
scale =  1.9964722649280886e-07
scale =  9.982361324640443e-08
scale =  4.9911806623202216e-08
scale =  2.4955903311601108e-08
scale =  1.2477951655800554e-08
scale =  6.238975827900277e-09
scale =  3.1194879139501385e-09
scale =  1.5597439569750693e-09
scale =  7.798719784875346e-10
scale =  3.899359892437673e-10
scale =  1.9496799462188366e-10
scale =  9.748399731094183e-11
scale =  4.8741998655470914e-11
scale =  2.4370999327735457e-11
scale =  1.2185499663867729e-11
scale =  6.092749831933864e-12
scale =  3.046374915966932e-12
scale =  1.523187457983466e-12
scale =  7.61593728991733e-13
No more descent can be found