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

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

#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.017294037213205
cost at iteration 1 :  0.0004834070682129935
scale =  4.768125361662991
cost at iteration 2 :  0.00048011519350089313
scale =  4.534231415136696
cost at iteration 3 :  0.00047669716218491144
scale =  4.316155256294004
cost at iteration 4 :  0.000473148028814085
scale =  4.113771843070691
cost at iteration 5 :  0.00046946364120922617
scale =  3.9265249630949897
cost at iteration 6 :  0.0004656405916026433
scale =  3.75360331502335
cost at iteration 7 :  0.0004616761602526168
scale =  3.594064615563142
cost at iteration 8 :  0.000457568257567207
scale =  3.4469191638143375
cost at iteration 9 :  0.0004533153684271193
scale =  3.3111835819313975
cost at iteration 10 :  0.00044891650072658145
scale =  3.1859135216674606
cost at iteration 11 :  0.00044437113904143
scale =  3.0702220249912244
cost at iteration 12 :  0.0004396792036418338
scale =  2.9632883913475627
cost at iteration 13 :  0.0004348410146650856
scale =  2.864360957277326
cost at iteration 14 :  0.0004298572610523669
scale =  2.772756118957969
cost at iteration 15 :  0.000424728973761221
scale =  2.6878551588136355
cost at iteration 16 :  0.0004194575027445427
cost at iteration 17 :  0.00041939466281839166
cost at iteration 18 :  0.0004157449599399011
cost at iteration 19 :  0.0004120329975392389
cost at iteration 20 :  0.0004082594338459176
cost at iteration 21 :  0.0004044249768264307
cost at iteration 22 :  0.00040053038287394535
cost at iteration 23 :  0.00039657645561609747
cost at iteration 24 :  0.00039256404482998665
cost at iteration 25 :  0.0003884940454543323
cost at iteration 26 :  0.0003843673966895963
cost at iteration 27 :  0.0003801850811776644
cost at iteration 28 :  0.00037594812425343233
cost at iteration 29 :  0.0003716575932613135
cost at iteration 30 :  0.0003673145969303526
cost at iteration 31 :  0.000362920284802181
cost at iteration 32 :  0.00035847584670661354
cost at iteration 33 :  0.0003539825122801605
cost at iteration 34 :  0.00034944155052317395
cost at iteration 35 :  0.0003448542693917542
cost at iteration 36 :  0.00034022201542090115
cost at iteration 37 :  0.0003355461733757329
cost at iteration 38 :  0.00033082816592787796
cost at iteration 39 :  0.00032606945335443034
cost at iteration 40 :  0.00032127153325708127
cost at iteration 41 :  0.0003164359402992843
cost at iteration 42 :  0.00031156424595946597
cost at iteration 43 :  0.0003066580582985117
cost at iteration 44 :  0.00030171902173988636
cost at iteration 45 :  0.00029674881686089674
cost at iteration 46 :  0.0002917491601937368
cost at iteration 47 :  0.00028672180403506014
cost at iteration 48 :  0.0002816685362629307
cost at iteration 49 :  0.00027659118016009314
cost at iteration 50 :  0.00027149159424258337
cost at iteration 51 :  0.0002663716720927739
cost at iteration 52 :  0.0002612333421960104
cost at iteration 53 :  0.00025607856778005206
cost at iteration 54 :  0.0002509093466565864
cost at iteration 55 :  0.0002457277110641181
cost at iteration 56 :  0.00024053572751158258
cost at iteration 57 :  0.00023533549662206588
cost at iteration 58 :  0.0002301291529760314
cost at iteration 59 :  0.00022491886495349067
cost at iteration 60 :  0.00021970683457456404
cost at iteration 61 :  0.0002144952973378939
cost at iteration 62 :  0.00020928652205638775
cost at iteration 63 :  0.00020408281068976707
cost at iteration 64 :  0.00019888649817340466
cost at iteration 65 :  0.000193699952242926
cost at iteration 66 :  0.0001885255732540454
cost at iteration 67 :  0.00018336579399709017
cost at iteration 68 :  0.00017822307950564944
cost at iteration 69 :  0.00017309992685875998
cost at iteration 70 :  0.00016799886497600075
cost at iteration 71 :  0.00016292245440484255
cost at iteration 72 :  0.00015787328709952612
cost at iteration 73 :  0.00015285398619070933
cost at iteration 74 :  0.00014786720574502425
cost at iteration 75 :  0.00014291563051361736
cost at iteration 76 :  0.00013800197566863134
cost at iteration 77 :  0.00013312898652646714
cost at iteration 78 :  0.00012829943825651012
cost at iteration 79 :  0.00012351613557383916
cost at iteration 80 :  0.00011878191241421155
cost at iteration 81 :  0.0001140996315893716
cost at iteration 82 :  0.00010947218442042489
cost at iteration 83 :  0.00010490249034664924
cost at iteration 84 :  0.00010039349650667475
cost at iteration 85 :  9.594817728841891e-05
cost at iteration 86 :  9.156953384351473e-05
cost at iteration 87 :  8.726059356115525e-05
cost at iteration 88 :  8.302440949529276e-05
cost at iteration 89 :  7.886405973789807e-05
cost at iteration 90 :  7.478264672946268e-05
cost at iteration 91 :  7.078329649600063e-05
cost at iteration 92 :  6.686915779939638e-05
cost at iteration 93 :  6.304340118486987e-05
cost at iteration 94 :  5.930921790538861e-05
cost at iteration 95 :  5.566981869778346e-05
cost at iteration 96 :  5.212843237870052e-05
cost at iteration 97 :  4.8688304219822557e-05
cost at iteration 98 :  4.5352694050244104e-05
cost at iteration 99 :  4.212487401837078e-05
cost at iteration 100 :  3.900812592466782e-05
cost at iteration 101 :  3.600573800762931e-05
cost at iteration 102 :  3.312100102499742e-05
cost at iteration 103 :  3.035720341522574e-05
cost at iteration 104 :  2.7717625242211967e-05
cost at iteration 105 :  2.5205530506518516e-05
cost at iteration 106 :  2.2824157227798994e-05
cost at iteration 107 :  2.0576704431640157e-05
cost at iteration 108 :  1.8466314751504515e-05
cost at iteration 109 :  1.6496050682175127e-05
cost at iteration 110 :  1.4668861414742696e-05
cost at iteration 111 :  1.298753531131632e-05
cost at iteration 112 :  1.1454629801791128e-05
cost at iteration 113 :  1.007236453642305e-05
cost at iteration 114 :  8.84245238994951e-06
cost at iteration 115 :  7.76582080710459e-06
cost at iteration 116 :  6.842130906123434e-06
cost at iteration 117 :  6.0689083805032874e-06
cost at iteration 118 :  5.439917433407755e-06
cost at iteration 119 :  4.9421723157369936e-06
cost at iteration 120 :  4.551597065260479e-06
cost at iteration 121 :  4.233116934342925e-06
cost at iteration 122 :  3.956671493976123e-06
cost at iteration 123 :  3.705866325986535e-06
cost at iteration 124 :  3.4699473518204085e-06
cost at iteration 125 :  3.244291371189636e-06
cost at iteration 126 :  3.0276161030242066e-06
cost at iteration 127 :  2.8195621995451558e-06
cost at iteration 128 :  2.619995179290807e-06
cost at iteration 129 :  2.428834172949116e-06
cost at iteration 130 :  2.2460123446618305e-06
cost at iteration 131 :  2.071468873280284e-06
cost at iteration 132 :  1.9051474673414027e-06
cost at iteration 133 :  1.7469959100276685e-06
cost at iteration 134 :  1.596965635968899e-06
cost at iteration 135 :  1.4550111938767357e-06
cost at iteration 136 :  1.3210895279494578e-06
cost at iteration 137 :  1.1951590291224589e-06
cost at iteration 138 :  1.077178227108375e-06
cost at iteration 139 :  9.671040317734376e-07
cost at iteration 140 :  8.64889191043892e-07
cost at iteration 141 :  7.704788943285349e-07
cost at iteration 142 :  6.83805398339974e-07
cost at iteration 143 :  6.047818021704446e-07
cost at iteration 144 :  5.332885422350769e-07
cost at iteration 145 :  4.6917227869682475e-07
cost at iteration 146 :  4.1221064534116917e-07
cost at iteration 147 :  3.6318560751780237e-07
cost at iteration 148 :  3.3832269759954074e-07
cost at iteration 149 :  3.275058052712108e-07
cost at iteration 150 :  3.18056677854775e-07
cost at iteration 151 :  3.069752795831818e-07
cost at iteration 152 :  2.987317167095616e-07
cost at iteration 153 :  2.88762791106612e-07
cost at iteration 154 :  2.8148128290096083e-07
cost at iteration 155 :  2.7245365256736495e-07
cost at iteration 156 :  2.65971192356009e-07
cost at iteration 157 :  2.5774588357950853e-07
cost at iteration 158 :  2.5194155114741505e-07
cost at iteration 159 :  2.444045295169949e-07
cost at iteration 160 :  2.391833910698957e-07
cost at iteration 161 :  2.3224044943596787e-07
cost at iteration 162 :  2.275255100760796e-07
cost at iteration 163 :  2.2109842023451852e-07
cost at iteration 164 :  2.1682601248735184e-07
cost at iteration 165 :  2.108495213807936e-07
cost at iteration 166 :  2.0696635118587437e-07
cost at iteration 167 :  2.0138579080237439e-07
cost at iteration 168 :  1.978468756940186e-07
cost at iteration 169 :  1.9261624199037055e-07
cost at iteration 170 :  1.8938338020170064e-07
cost at iteration 171 :  1.8446378041135094e-07
cost at iteration 172 :  1.8150436523438398e-07
cost at iteration 173 :  1.7686275872046816e-07
cost at iteration 174 :  1.7414883448150097e-07
cost at iteration 175 :  1.697570093206202e-07
cost at iteration 176 :  1.6726450624787044e-07
cost at iteration 177 :  1.6309824594485587e-07
cost at iteration 178 :  1.6080635297724338e-07
cost at iteration 179 :  1.568447569365115e-07
cost at iteration 180 :  1.547354039001502e-07
cost at iteration 181 :  1.509603324548927e-07
cost at iteration 182 :  1.490177605659399e-07
cost at iteration 183 :  1.4541338104036126e-07
cost at iteration 184 :  1.4362378559478752e-07
cost at iteration 185 :  1.401762004030865e-07
cost at iteration 186 :  1.385274329083882e-07
cost at iteration 187 :  1.3522437432807802e-07
cost at iteration 188 :  1.3370569381075814e-07
cost at iteration 189 :  1.3053627299333508e-07
cost at iteration 190 :  1.29138138104063e-07
cost at iteration 191 :  1.2609263824606106e-07
cost at iteration 192 :  1.2480653326816337e-07
cost at iteration 193 :  1.2187623877206758e-07
cost at iteration 194 :  1.206945278310695e-07
cost at iteration 195 :  1.1787158282715276e-07
cost at iteration 196 :  1.1678738757095116e-07
cost at iteration 197 :  1.140646784185552e-07
cost at iteration 198 :  1.1307177523783776e-07
cost at iteration 199 :  1.1044283263470554e-07
cost at iteration 200 :  1.0953556615558988e-07
cost at iteration 201 :  1.0699448330215464e-07
cost at iteration 202 :  1.061676934332881e-07
cost at iteration 203 :  1.0370905736211274e-07
cost at iteration 204 :  1.0295801763640936e-07
cost at iteration 205 :  1.0057685135480726e-07
cost at iteration 206 :  9.989721668726827e-08
cost at iteration 207 :  9.758893021755039e-08
cost at iteration 208 :  9.697669251815788e-08
cost at iteration 209 :  9.473704127411287e-08
cost at iteration 210 :  9.418849161904558e-08
cost at iteration 211 :  9.201354084489795e-08
cost at iteration 212 :  9.152523712911937e-08
cost at iteration 213 :  8.941133136091306e-08
cost at iteration 214 :  8.89800705379105e-08
cost at iteration 215 :  8.692380723721378e-08
cost at iteration 216 :  8.654660140337719e-08
cost at iteration 217 :  8.45448080676708e-08
cost at iteration 218 :  8.421886377473977e-08
cost at iteration 219 :  8.226857795460033e-08
cost at iteration 220 :  8.199127823801695e-08
cost at iteration 221 :  8.008972999351196e-08
cost at iteration 222 :  7.98586186910983e-08
cost at iteration 223 :  7.800321510320754e-08
cost at iteration 224 :  7.781598311026641e-08
cost at iteration 225 :  7.600429453110333e-08
cost at iteration 226 :  7.585876769743935e-08
cost at iteration 227 :  7.408851547849459e-08
cost at iteration 228 :  7.398264390199792e-08
cost at iteration 229 :  7.225168938487363e-08
cost at iteration 230 :  7.218353789703687e-08
cost at iteration 231 :  7.048987248811821e-08
cost at iteration 232 :  7.045761216052638e-08
cost at iteration 233 :  6.879934834129154e-08
scale =  53.59239427195622
scale =  26.79619713597811
scale =  13.398098567989056
scale =  6.699049283994528
scale =  3.349524641997264
scale =  1.674762320998632
scale =  0.837381160499316
scale =  0.418690580249658
scale =  0.209345290124829
scale =  0.1046726450624145
scale =  0.05233632253120725
scale =  0.026168161265603625
scale =  0.013084080632801812
scale =  0.006542040316400906
scale =  0.003271020158200453
scale =  0.0016355100791002265
scale =  0.0008177550395501133
scale =  0.00040887751977505664
scale =  0.00020443875988752832
scale =  0.00010221937994376416
scale =  5.110968997188208e-05
scale =  2.555484498594104e-05
scale =  1.277742249297052e-05
scale =  6.38871124648526e-06
scale =  3.19435562324263e-06
scale =  1.597177811621315e-06
scale =  7.985889058106575e-07
scale =  3.9929445290532875e-07
scale =  1.9964722645266437e-07
scale =  9.982361322633219e-08
scale =  4.991180661316609e-08
scale =  2.4955903306583047e-08
scale =  1.2477951653291523e-08
scale =  6.238975826645762e-09
scale =  3.119487913322881e-09
scale =  1.5597439566614404e-09
scale =  7.798719783307202e-10
scale =  3.899359891653601e-10
scale =  1.9496799458268005e-10
scale =  9.748399729134003e-11
scale =  4.874199864567001e-11
scale =  2.4370999322835006e-11
scale =  1.2185499661417503e-11
scale =  6.092749830708752e-12
scale =  3.046374915354376e-12
scale =  1.523187457677188e-12
scale =  7.61593728838594e-13
No more descent can be found