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")

BaseWebGuiScene

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.0004865785615789076

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.00032579700603296694

#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.00048657856157890744

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.0004865785615789076
scale =  5.0172940372132
cost at iteration 1 :  0.0004834070731679914
scale =  4.76812536166299
cost at iteration 2 :  0.00048011519918462055
scale =  4.534231415136696
cost at iteration 3 :  0.0004766971686029877
scale =  4.316155256294001
cost at iteration 4 :  0.00047314803597050544
scale =  4.113771843070691
cost at iteration 5 :  0.0004694636491068562
scale =  3.9265249630949883
cost at iteration 6 :  0.00046564060024363287
scale =  3.7536033150233488
cost at iteration 7 :  0.00046167616963875113
scale =  3.5940646155631444
cost at iteration 8 :  0.00045756826770020266
scale =  3.446919163814335
cost at iteration 9 :  0.00045331537930887715
scale =  3.3111835819313984
cost at iteration 10 :  0.00044891651235940577
scale =  3.1859135216674597
cost at iteration 11 :  0.0004443711514282213
scale =  3.070222024991224
cost at iteration 12 :  0.0004396792167862605
scale =  2.9632883913475627
cost at iteration 13 :  0.00043484102857173975
scale =  2.8643609572773254
cost at iteration 14 :  0.0004298572757269101
scale =  2.7727561189579695
cost at iteration 15 :  0.0004247289892105183
scale =  2.6878551588136355
cost at iteration 16 :  0.00041945751897679365
cost at iteration 17 :  0.00041939467924363087
cost at iteration 18 :  0.0004157449768855942
cost at iteration 19 :  0.00041203301501048864
cost at iteration 20 :  0.00040825945184834076
cost at iteration 21 :  0.0004044249953661821
cost at iteration 22 :  0.0004005304019577407
cost at iteration 23 :  0.0003965764752512381
cost at iteration 24 :  0.00039256406502438353
cost at iteration 25 :  0.0003884940662165308
cost at iteration 26 :  0.00038436741802880004
cost at iteration 27 :  0.00038018510310376143
cost at iteration 28 :  0.0003759481467770182
cost at iteration 29 :  0.00037165761639372004
cost at iteration 30 :  0.00036731462068367155
cost at iteration 31 :  0.0003629203091892909
cost at iteration 32 :  0.0003584758717412077
cost at iteration 33 :  0.0003539825379767731
cost at iteration 34 :  0.0003494415768972085
cost at iteration 35 :  0.00034485429645951125
cost at iteration 36 :  0.00034022204319960784
cost at iteration 37 :  0.00033554620188357094
cost at iteration 38 :  0.0003308281951840144
cost at iteration 39 :  0.0003260694833790455
cost at iteration 40 :  0.0003212715640714041
cost at iteration 41 :  0.0003164359719256182
cost at iteration 42 :  0.00031156427842122375
cost at iteration 43 :  0.00030665809162024613
cost at iteration 44 :  0.0003017190559473237
cost at iteration 45 :  0.00029674885198096917
cost at iteration 46 :  0.0002917491962546171
cost at iteration 47 :  0.00028672184106619454
cost at iteration 48 :  0.0002816685742950744
cost at iteration 49 :  0.0002765912192253453
cost at iteration 50 :  0.00027149163437442225
cost at iteration 51 :  0.0002663717133260934
cost at iteration 52 :  0.00026123338456715594
cost at iteration 53 :  0.0002560786113268585
cost at iteration 54 :  0.0002509093914184146
cost at iteration 55 :  0.0002457277570818929
cost at iteration 56 :  0.00024053577482783175
cost at iteration 57 :  0.0002353355452809573
cost at iteration 58 :  0.00023012920302341262
cost at iteration 59 :  0.0002249189164369292
cost at iteration 60 :  0.00021970688754338494
cost at iteration 61 :  0.00021449535184322124
cost at iteration 62 :  0.00020928657815118374
cost at iteration 63 :  0.00020408286842887272
cost at iteration 64 :  0.00019888655761358023
cost at iteration 65 :  0.00019370001344289143
cost at iteration 66 :  0.00018852563627452015
cost at iteration 67 :  0.00018336585890083392
cost at iteration 68 :  0.00017822314635750378
cost at iteration 69 :  0.00017309999572568617
cost at iteration 70 :  0.00016799893592712244
cost at iteration 71 :  0.00016292252751148233
cost at iteration 72 :  0.00015787336243524718
cost at iteration 73 :  0.00015285406383135154
cost at iteration 74 :  0.00014786728576874348
cost at iteration 75 :  0.00014291571300092203
cost at iteration 76 :  0.00013800206070241748
cost at iteration 77 :  0.0001331290741920521
cost at iteration 78 :  0.0001282995286416666
cost at iteration 79 :  0.0001235162287688252
cost at iteration 80 :  0.00011878200851179973
cost at iteration 81 :  0.00011409973068487517
cost at iteration 82 :  0.0001094722866117213
cost at iteration 83 :  0.00010490259573420029
cost at iteration 84 :  0.00010039360519354304
cost at iteration 85 :  9.59482893802794e-05
cost at iteration 86 :  9.156964944866149e-05
cost at iteration 87 :  8.726071279050268e-05
cost at iteration 88 :  8.302453246236891e-05
cost at iteration 89 :  7.886418655883118e-05
cost at iteration 90 :  7.478277752295755e-05
cost at iteration 91 :  7.078343138330443e-05
cost at iteration 92 :  6.686929690425188e-05
cost at iteration 93 :  6.304354463345377e-05
cost at iteration 94 :  5.9309365826232125e-05
cost at iteration 95 :  5.566997122167291e-05
cost at iteration 96 :  5.212858963855242e-05
cost at iteration 97 :  4.868846635053017e-05
cost at iteration 98 :  4.5352861188487664e-05
cost at iteration 99 :  4.2125046302386385e-05
cost at iteration 100 :  3.9008303493963315e-05
cost at iteration 101 :  3.600592100263525e-05
cost at iteration 102 :  3.312118958664038e-05
cost at iteration 103 :  3.0357397684400802e-05
cost at iteration 104 :  2.7717825359129246e-05
cost at iteration 105 :  2.520573660988727e-05
cost at iteration 106 :  2.28243694538019e-05
cost at iteration 107 :  2.0576922912627565e-05
cost at iteration 108 :  1.8466539614314545e-05
cost at iteration 109 :  1.6496282045939493e-05
cost at iteration 110 :  1.4669099387971284e-05
cost at iteration 111 :  1.2987779987943602e-05
cost at iteration 112 :  1.1454881255687179e-05
cost at iteration 113 :  1.0072622813555762e-05
cost at iteration 114 :  8.842717496818089e-06
cost at iteration 115 :  7.7660926930044e-06
cost at iteration 116 :  6.842409434733028e-06
cost at iteration 117 :  6.06919328215366e-06
cost at iteration 118 :  5.440208221421105e-06
cost at iteration 119 :  4.9424681388560385e-06
cost at iteration 120 :  4.551896484353188e-06
cost at iteration 121 :  4.23341783311608e-06
cost at iteration 122 :  3.956971777491843e-06
cost at iteration 123 :  3.7061649322349924e-06
cost at iteration 124 :  3.470244162683969e-06
cost at iteration 125 :  3.244586663284191e-06
cost at iteration 126 :  3.0279102205696687e-06
cost at iteration 127 :  2.8198554430474035e-06
cost at iteration 128 :  2.6202877787151974e-06
cost at iteration 129 :  2.4291262848947427e-06
cost at iteration 130 :  2.246304055024347e-06
cost at iteration 131 :  2.0717602007679703e-06
cost at iteration 132 :  1.9054383669070419e-06
cost at iteration 133 :  1.7472862760882244e-06
cost at iteration 134 :  1.59725530543601e-06
cost at iteration 135 :  1.455299949046631e-06
cost at iteration 136 :  1.321377099328659e-06
cost at iteration 137 :  1.1954450982857058e-06
cost at iteration 138 :  1.0774624297453897e-06
cost at iteration 139 :  9.67385961114444e-07
cost at iteration 140 :  8.651684020082311e-07
cost at iteration 141 :  7.707549087236222e-07
cost at iteration 142 :  6.840777119748682e-07
cost at iteration 143 :  6.050498938947853e-07
cost at iteration 144 :  5.335518884441783e-07
cost at iteration 145 :  4.694303650477785e-07
cost at iteration 146 :  4.124630156896632e-07
cost at iteration 147 :  3.634318283779554e-07
cost at iteration 148 :  3.3856406608880464e-07
cost at iteration 149 :  3.277442264944354e-07
cost at iteration 150 :  3.182937313931422e-07
cost at iteration 151 :  3.0721010343773706e-07
cost at iteration 152 :  2.989652581765771e-07
cost at iteration 153 :  2.889941666590611e-07
cost at iteration 154 :  2.817114464355738e-07
cost at iteration 155 :  2.726817143271275e-07
cost at iteration 156 :  2.661981044975951e-07
cost at iteration 157 :  2.5797075848524127e-07
cost at iteration 158 :  2.5216533373766275e-07
cost at iteration 159 :  2.446263392803478e-07
cost at iteration 160 :  2.394041620933194e-07
cost at iteration 161 :  2.324593113446515e-07
cost at iteration 162 :  2.2774338378630533e-07
cost at iteration 163 :  2.2131444737105525e-07
cost at iteration 164 :  2.1704109930712592e-07
cost at iteration 165 :  2.110628226381814e-07
cost at iteration 166 :  2.0717875757032914e-07
cost at iteration 167 :  2.0159647085373775e-07
cost at iteration 168 :  1.9805670402883647e-07
cost at iteration 169 :  1.92824401275835e-07
cost at iteration 170 :  1.8959072875850023e-07
cost at iteration 171 :  1.8466951516401852e-07
cost at iteration 172 :  1.817093281847271e-07
cost at iteration 173 :  1.7706616103501243e-07
cost at iteration 174 :  1.7435150196244062e-07
cost at iteration 175 :  1.6995816726006497e-07
cost at iteration 176 :  1.6746496446916387e-07
cost at iteration 177 :  1.632972436781887e-07
cost at iteration 178 :  1.610046843600991e-07
cost at iteration 179 :  1.570416748992157e-07
cost at iteration 180 :  1.5493168723828189e-07
cost at iteration 181 :  1.5115524752500243e-07
cost at iteration 182 :  1.4921207120100335e-07
cost at iteration 183 :  1.4560636672360343e-07
cost at iteration 184 :  1.4381619559965773e-07
cost at iteration 185 :  1.403673270218137e-07
cost at iteration 186 :  1.3871801127306833e-07
cost at iteration 187 :  1.3541370920973796e-07
cost at iteration 188 :  1.3389450662687736e-07
cost at iteration 189 :  1.3072388065520798e-07
cost at iteration 190 :  1.293252487449075e-07
cost at iteration 191 :  1.2627858057384486e-07
cost at iteration 192 :  1.2499200256204807e-07
cost at iteration 193 :  1.2206057519058468e-07
cost at iteration 194 :  1.2087841422651111e-07
cost at iteration 195 :  1.1805437046203054e-07
cost at iteration 196 :  1.1696974729268435e-07
cost at iteration 197 :  1.1424597224821e-07
cost at iteration 198 :  1.1325266243311949e-07
cost at iteration 199 :  1.1062268563230833e-07
cost at iteration 200 :  1.0971503303058768e-07
cost at iteration 201 :  1.0717294656753895e-07
cost at iteration 202 :  1.0634579037961123e-07
cost at iteration 203 :  1.0388618024380602e-07
cost at iteration 204 :  1.0313479334798823e-07
cost at iteration 205 :  1.0075268156251197e-07
cost at iteration 206 :  1.0007271826797578e-07
cost at iteration 207 :  9.77635139255019e-08
cost at iteration 208 :  9.715096558062494e-08
cost at iteration 209 :  9.491042321582604e-08
cost at iteration 210 :  9.436158037516989e-08
cost at iteration 211 :  9.218576439983938e-08
cost at iteration 212 :  9.169718447284387e-08
cost at iteration 213 :  8.958243863367349e-08
cost at iteration 214 :  8.915091812080677e-08
cost at iteration 215 :  8.709383912973698e-08
cost at iteration 216 :  8.671638970360156e-08
cost at iteration 217 :  8.471380434510504e-08
cost at iteration 218 :  8.438763215990687e-08
cost at iteration 219 :  8.243657730530195e-08
cost at iteration 220 :  8.215906502250728e-08
cost at iteration 221 :  8.025677008371996e-08
cost at iteration 222 :  8.002546118829617e-08
cost at iteration 223 :  7.816933262686722e-08
cost at iteration 224 :  7.798191768019317e-08
cost at iteration 225 :  7.616952525530332e-08
cost at iteration 226 :  7.60238297902203e-08
cost at iteration 227 :  7.425289428491623e-08
cost at iteration 228 :  7.414686809755577e-08
cost at iteration 229 :  7.241525030766198e-08
cost at iteration 230 :  7.234695794138852e-08
cost at iteration 231 :  7.065264874851825e-08
cost at iteration 232 :  7.062026099901922e-08
cost at iteration 233 :  6.896137237937607e-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