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

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

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.0004834070731679915
scale =  4.76812536166299
cost at iteration 2 :  0.00048011519918462055
scale =  4.534231415136698
cost at iteration 3 :  0.0004766971686029878
scale =  4.316155256294001
cost at iteration 4 :  0.00047314803597050544
scale =  4.11377184307069
cost at iteration 5 :  0.0004694636491068565
scale =  3.926524963094991
cost at iteration 6 :  0.000465640600243633
scale =  3.75360331502335
cost at iteration 7 :  0.0004616761696387511
scale =  3.594064615563144
cost at iteration 8 :  0.00045756826770020255
scale =  3.4469191638143335
cost at iteration 9 :  0.0004533153793088772
scale =  3.3111835819313975
cost at iteration 10 :  0.0004489165123594056
scale =  3.1859135216674592
cost at iteration 11 :  0.00044437115142822135
scale =  3.070222024991224
cost at iteration 12 :  0.00043967921678626044
scale =  2.9632883913475636
cost at iteration 13 :  0.0004348410285717398
scale =  2.8643609572773268
cost at iteration 14 :  0.00042985727572691
scale =  2.772756118957971
cost at iteration 15 :  0.0004247289892105181
scale =  2.687855158813637
cost at iteration 16 :  0.0004194575189767934

cost at iteration 17 :  0.00041939467924363054
cost at iteration 18 :  0.0004157449768855941
cost at iteration 19 :  0.00041203301501048853
cost at iteration 20 :  0.00040825945184834087
cost at iteration 21 :  0.000404424995366182
cost at iteration 22 :  0.0004005304019577407
cost at iteration 23 :  0.000396576475251238
cost at iteration 24 :  0.00039256406502438375
cost at iteration 25 :  0.00038849406621653087
cost at iteration 26 :  0.00038436741802880004
cost at iteration 27 :  0.00038018510310376143
cost at iteration 28 :  0.0003759481467770181
cost at iteration 29 :  0.00037165761639371993
cost at iteration 30 :  0.00036731462068367117
cost at iteration 31 :  0.0003629203091892907
cost at iteration 32 :  0.00035847587174120746
cost at iteration 33 :  0.00035398253797677324
cost at iteration 34 :  0.0003494415768972088
cost at iteration 35 :  0.00034485429645951136
cost at iteration 36 :  0.0003402220431996079
cost at iteration 37 :  0.00033554620188357083

cost at iteration 38 :  0.0003308281951840142
cost at iteration 39 :  0.00032606948337904595
cost at iteration 40 :  0.00032127156407140415
cost at iteration 41 :  0.0003164359719256182
cost at iteration 42 :  0.00031156427842122364
cost at iteration 43 :  0.00030665809162024635
cost at iteration 44 :  0.00030171905594732364
cost at iteration 45 :  0.00029674885198096917
cost at iteration 46 :  0.00029174919625461694
cost at iteration 47 :  0.00028672184106619454
cost at iteration 48 :  0.0002816685742950742
cost at iteration 49 :  0.0002765912192253451
cost at iteration 50 :  0.0002714916343744222
cost at iteration 51 :  0.00026637171332609324
cost at iteration 52 :  0.0002612333845671558
cost at iteration 53 :  0.0002560786113268586
cost at iteration 54 :  0.00025090939141841456
cost at iteration 55 :  0.000245727757081893
cost at iteration 56 :  0.0002405357748278317
cost at iteration 57 :  0.00023533554528095698
cost at iteration 58 :  0.00023012920302341265
cost at iteration 59 :  0.00022491891643692898
cost at iteration 60 :  0.00021970688754338494
cost at iteration 61 :  0.00021449535184322118
cost at iteration 62 :  0.00020928657815118376

cost at iteration 63 :  0.000204082868428873
cost at iteration 64 :  0.00019888655761358032
cost at iteration 65 :  0.00019370001344289149
cost at iteration 66 :  0.00018852563627451988
cost at iteration 67 :  0.00018336585890083409
cost at iteration 68 :  0.00017822314635750365
cost at iteration 69 :  0.00017309999572568627
cost at iteration 70 :  0.00016799893592712246
cost at iteration 71 :  0.00016292252751148238
cost at iteration 72 :  0.00015787336243524693
cost at iteration 73 :  0.0001528540638313514
cost at iteration 74 :  0.00014786728576874343
cost at iteration 75 :  0.000142915713000922
cost at iteration 76 :  0.00013800206070241727
cost at iteration 77 :  0.00013312907419205217
cost at iteration 78 :  0.00012829952864166661
cost at iteration 79 :  0.0001235162287688251
cost at iteration 80 :  0.00011878200851179958
cost at iteration 81 :  0.00011409973068487521
cost at iteration 82 :  0.00010947228661172123
cost at iteration 83 :  0.00010490259573420047
cost at iteration 84 :  0.00010039360519354301
cost at iteration 85 :  9.594828938027915e-05
cost at iteration 86 :  9.156964944866153e-05
cost at iteration 87 :  8.726071279050268e-05
cost at iteration 88 :  8.302453246236872e-05

cost at iteration 89 :  7.88641865588312e-05
cost at iteration 90 :  7.478277752295763e-05
cost at iteration 91 :  7.078343138330423e-05
cost at iteration 92 :  6.686929690425183e-05
cost at iteration 93 :  6.304354463345361e-05
cost at iteration 94 :  5.9309365826231996e-05
cost at iteration 95 :  5.56699712216729e-05
cost at iteration 96 :  5.212858963855236e-05
cost at iteration 97 :  4.868846635053011e-05
cost at iteration 98 :  4.535286118848774e-05
cost at iteration 99 :  4.212504630238634e-05
cost at iteration 100 :  3.900830349396327e-05
cost at iteration 101 :  3.6005921002635226e-05
cost at iteration 102 :  3.3121189586640336e-05
cost at iteration 103 :  3.0357397684400822e-05
cost at iteration 104 :  2.7717825359129206e-05
cost at iteration 105 :  2.5205736609887297e-05
cost at iteration 106 :  2.2824369453801957e-05
cost at iteration 107 :  2.057692291262757e-05
cost at iteration 108 :  1.8466539614314582e-05
cost at iteration 109 :  1.6496282045939547e-05
cost at iteration 110 :  1.4669099387971241e-05
cost at iteration 111 :  1.298777998794354e-05
cost at iteration 112 :  1.145488125568711e-05
cost at iteration 113 :  1.007262281355575e-05
cost at iteration 114 :  8.842717496818089e-06

cost at iteration 115 :  7.766092693004397e-06
cost at iteration 116 :  6.842409434732994e-06
cost at iteration 117 :  6.069193282153669e-06
cost at iteration 118 :  5.440208221421094e-06
cost at iteration 119 :  4.9424681388560335e-06
cost at iteration 120 :  4.551896484353182e-06
cost at iteration 121 :  4.233417833116078e-06
cost at iteration 122 :  3.956971777491831e-06
cost at iteration 123 :  3.7061649322349903e-06
cost at iteration 124 :  3.470244162683972e-06
cost at iteration 125 :  3.2445866632841836e-06
cost at iteration 126 :  3.027910220569668e-06
cost at iteration 127 :  2.8198554430473946e-06
cost at iteration 128 :  2.6202877787151928e-06
cost at iteration 129 :  2.429126284894742e-06
cost at iteration 130 :  2.2463040550243446e-06
cost at iteration 131 :  2.0717602007679695e-06
cost at iteration 132 :  1.9054383669070406e-06
cost at iteration 133 :  1.7472862760882225e-06
cost at iteration 134 :  1.5972553054360095e-06
cost at iteration 135 :  1.4552999490466279e-06
cost at iteration 136 :  1.32137709932866e-06
cost at iteration 137 :  1.1954450982857e-06
cost at iteration 138 :  1.0774624297453983e-06
cost at iteration 139 :  9.673859611144226e-07
cost at iteration 140 :  8.651684020082804e-07

cost at iteration 141 :  7.70754908723486e-07
cost at iteration 142 :  6.840777119752722e-07
cost at iteration 143 :  6.050498938933369e-07
cost at iteration 144 :  5.33551888448009e-07
cost at iteration 145 :  4.6943036500608073e-07
cost at iteration 146 :  4.1246301529351296e-07
cost at iteration 147 :  3.6343181510421024e-07
cost at iteration 148 :  3.3856392909400196e-07
cost at iteration 149 :  3.2774413393563646e-07
cost at iteration 150 :  3.1829363965946713e-07
cost at iteration 151 :  3.0721003489110865e-07
cost at iteration 152 :  2.989651852576779e-07
cost at iteration 153 :  2.8899411214234657e-07
cost at iteration 154 :  2.8171138616899657e-07
cost at iteration 155 :  2.7268166895265493e-07
cost at iteration 156 :  2.661980533481267e-07
cost at iteration 157 :  2.57970719548692e-07
cost at iteration 158 :  2.5216528948974876e-07
cost at iteration 159 :  2.446263051437108e-07
cost at iteration 160 :  2.3940412326028583e-07
cost at iteration 161 :  2.324592809429861e-07
cost at iteration 162 :  2.2774334931668095e-07
cost at iteration 163 :  2.21314419971787e-07
cost at iteration 164 :  2.1704106842657958e-07
cost at iteration 165 :  2.110627977140854e-07
cost at iteration 166 :  2.0717872969058573e-07

cost at iteration 167 :  2.0159644801141513e-07
cost at iteration 168 :  1.980566786919977e-07
cost at iteration 169 :  1.92824380213092e-07
cost at iteration 170 :  1.8959070560123464e-07
cost at iteration 171 :  1.8466949564320423e-07
cost at iteration 172 :  1.8170930691420843e-07
cost at iteration 173 :  1.7706614286558312e-07
cost at iteration 174 :  1.7435148233958673e-07
cost at iteration 175 :  1.6995815028675377e-07
cost at iteration 176 :  1.674649462965877e-07
cost at iteration 177 :  1.632972277727086e-07
cost at iteration 178 :  1.6100466747324407e-07
cost at iteration 179 :  1.5704165995428115e-07
cost at iteration 180 :  1.5493167149874422e-07
cost at iteration 181 :  1.511552334499089e-07
cost at iteration 182 :  1.4921205649143035e-07
cost at iteration 183 :  1.4560635344090982e-07
cost at iteration 184 :  1.4381618181979e-07
cost at iteration 185 :  1.4036731446480155e-07
cost at iteration 186 :  1.3871799833663546e-07
cost at iteration 187 :  1.3541369732043512e-07
cost at iteration 188 :  1.3389449445913393e-07
cost at iteration 189 :  1.3072386938285008e-07
cost at iteration 190 :  1.293252372806717e-07
cost at iteration 191 :  1.262785698736442e-07
cost at iteration 192 :  1.249919917441193e-07

cost at iteration 193 :  1.2206056502275454e-07
cost at iteration 194 :  1.2087840400438844e-07
cost at iteration 195 :  1.1805436079099762e-07
cost at iteration 196 :  1.1696973762152907e-07
cost at iteration 197 :  1.1424596304197864e-07
cost at iteration 198 :  1.1325265327289772e-07
cost at iteration 199 :  1.1062267686193647e-07
cost at iteration 200 :  1.0971502434537597e-07
cost at iteration 201 :  1.0717293820671682e-07
cost at iteration 202 :  1.0634578213701091e-07
cost at iteration 203 :  1.0388617226850391e-07
cost at iteration 204 :  1.0313478551864953e-07
cost at iteration 205 :  1.0075267395068628e-07
cost at iteration 206 :  1.0007271082519713e-07
cost at iteration 207 :  9.776350665685512e-08
cost at iteration 208 :  9.715095850001704e-08
cost at iteration 209 :  9.491041627160415e-08
cost at iteration 210 :  9.436157363437631e-08
cost at iteration 211 :  9.218575776265983e-08
cost at iteration 212 :  9.169717805129789e-08
cost at iteration 213 :  8.95824322873811e-08
cost at iteration 214 :  8.915091199953692e-08
cost at iteration 215 :  8.709383305927959e-08
cost at iteration 216 :  8.671638386505517e-08
cost at iteration 217 :  8.471379853642762e-08
cost at iteration 218 :  8.438762658780106e-08

cost at iteration 219 :  8.243657174525482e-08
cost at iteration 220 :  8.215905970170576e-08
cost at iteration 221 :  8.025676475998065e-08
cost at iteration 222 :  8.002545610469774e-08
cost at iteration 223 :  7.816932752786989e-08
cost at iteration 224 :  7.798191282063738e-08
cost at iteration 225 :  7.616952037017834e-08
cost at iteration 226 :  7.602382514240063e-08
cost at iteration 227 :  7.4252889603437e-08
cost at iteration 228 :  7.414686364995029e-08
cost at iteration 229 :  7.241524582019547e-08
cost at iteration 230 :  7.234695368319328e-08
cost at iteration 231 :  7.065264444598434e-08
cost at iteration 232 :  7.062025692009145e-08
cost at iteration 233 :  6.896136825320757e-08
scale =  53.59239366683971
scale =  26.796196833419856
scale =  13.398098416709928
scale =  6.699049208354964
scale =  3.349524604177482
scale =  1.674762302088741
scale =  0.8373811510443705
scale =  0.41869057552218525
scale =  0.20934528776109262
scale =  0.10467264388054631
scale =  0.052336321940273156
scale =  0.026168160970136578
scale =  0.013084080485068289
scale =  0.0065420402425341445
scale =  0.0032710201212670723
scale =  0.0016355100606335361
scale =  0.0008177550303167681
scale =  0.00040887751515838403
scale =  0.00020443875757919202
scale =  0.00010221937878959601
scale =  5.1109689394798004e-05
scale =  2.5554844697399002e-05
scale =  1.2777422348699501e-05
scale =  6.3887111743497505e-06
scale =  3.1943555871748752e-06
scale =  1.5971777935874376e-06
scale =  7.985888967937188e-07
scale =  3.992944483968594e-07
scale =  1.996472241984297e-07
scale =  9.982361209921485e-08
scale =  4.9911806049607426e-08
scale =  2.4955903024803713e-08
scale =  1.2477951512401856e-08
scale =  6.238975756200928e-09
scale =  3.119487878100464e-09
scale =  1.559743939050232e-09
scale =  7.79871969525116e-10
scale =  3.89935984762558e-10
scale =  1.94967992381279e-10
scale =  9.74839961906395e-11
scale =  4.874199809531975e-11
scale =  2.4370999047659876e-11
scale =  1.2185499523829938e-11
scale =  6.092749761914969e-12
scale =  3.0463748809574845e-12
scale =  1.5231874404787422e-12
scale =  7.615937202393711e-13
No more descent can be found