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

7.2 PDE-Constrained Shape Optimization in NGSolve

We want to solve the PDE-constrained shape optimization problem :nbsphinx-math:`begin{equation}

underset{Omegasubset mathsf{D}}{mbox{min}} ; J(u) := int_Omega |u-u_d|^2 ; dx

end{equation}` subject to that \((\Omega,u)\) satisfy :nbsphinx-math:`begin{equation}

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

end{equation}` 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 :nbsphinx-math:`begin{align*}

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)

end{align*}` where \(u\) is the solution to the state equation. For \(J(u) = \int_\Omega |u-u_d|^2 \mbox dx\), we get :nbsphinx-math:`begin{align*}

partial_u J(u)(w) = -2 int_Omega (u-u_d)w ,mbox dx.

end{align*}` 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 :nbsphinx-math:`begin{align*}

dmathcal 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 fcdot X + f mbox{div}(X)) p, dx.

end{align*}` 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 :nbsphinx-math:`begin{align*}

  mbox{Find } X in H: qquad b(X, Phi) = dmathcal J(Omega; Phi) ; text{ for all } Phi in H

  end{align*}` 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.0004865769485850591

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

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

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.0004865769485850591
scale =  5.016909006111095
cost at iteration 1 :  0.00048340549129873565
scale =  4.76845478507901
cost at iteration 2 :  0.0004801131572550958
scale =  4.534992340164933
cost at iteration 3 :  0.00047669426791743775
scale =  4.3171503211918045
cost at iteration 4 :  0.0004731439544112973
scale =  4.114867929335298
cost at iteration 5 :  0.0004694581277620631
scale =  3.9276354539113245
cost at iteration 6 :  0.00046563343240193296
scale =  3.754674311081731
cost at iteration 7 :  0.0004616671918286652
scale =  3.5950646910601494
cost at iteration 8 :  0.00045755735242073955
scale =  3.4478319860564635
cost at iteration 9 :  0.00045330242912636285
scale =  3.3120027141480977
cost at iteration 10 :  0.00044890145509412084
scale =  3.1866388283415388
cost at iteration 11 :  0.00044435393620606675
scale =  3.0708572236670593
cost at iteration 12 :  0.0004396598107770552
scale =  2.9638394098922025
cost at iteration 13 :  0.00043481941427441536
scale =  2.8648348493878424
cost at iteration 14 :  0.00042983344869279035
scale =  2.7731603623829457
cost at iteration 15 :  0.00042470295612007266
scale =  2.6881972134102
cost at iteration 16 :  0.00041942929600397456
cost at iteration 17 :  0.0004193665890304989
cost at iteration 18 :  0.00041571542101146927
cost at iteration 19 :  0.0004120020086060284
cost at iteration 20 :  0.00040822701180077015
cost at iteration 21 :  0.0004043911401214325
cost at iteration 22 :  0.00040049515134219215
cost at iteration 23 :  0.0003965398503109536
cost at iteration 24 :  0.000392526087880043
cost at iteration 25 :  0.0003884547599325343
cost at iteration 26 :  0.0003843268064952386
cost at iteration 27 :  0.0003801432109301419
cost at iteration 28 :  0.00037590499919678985
cost at iteration 29 :  0.00037161323917878634
cost at iteration 30 :  0.0003672690400681866
cost at iteration 31 :  0.0003628735518021317
cost at iteration 32 :  0.00035842796454659567
cost at iteration 33 :  0.0003539335082225842
cost at iteration 34 :  0.00034939145207056306
cost at iteration 35 :  0.0003448031042492831
cost at iteration 36 :  0.0003401698114655208
cost at iteration 37 :  0.00033549295863158455
cost at iteration 38 :  0.0003307739685477212
cost at iteration 39 :  0.0003260143016068122
cost at iteration 40 :  0.0003212154555190103
cost at iteration 41 :  0.0003163789650541484
cost at iteration 42 :  0.00031150640179997244
cost at iteration 43 :  0.00030659937393440174
cost at iteration 44 :  0.00030165952601019647
cost at iteration 45 :  0.00029668853875052653
cost at iteration 46 :  0.00029168812885408985
cost at iteration 47 :  0.0002866600488085142
cost at iteration 48 :  0.00028160608671089535
cost at iteration 49 :  0.00027652806609441223
cost at iteration 50 :  0.0002714278457600286
cost at iteration 51 :  0.00026630731961238194
cost at iteration 52 :  0.00026116841649900585
cost at iteration 53 :  0.00025601310005210395
cost at iteration 54 :  0.0002508433685321351
cost at iteration 55 :  0.0002456612546725164
cost at iteration 56 :  0.00024046882552478856
cost at iteration 57 :  0.0002352681823036236
cost at iteration 58 :  0.00023006146023108094
cost at iteration 59 :  0.00022485082837953827
cost at iteration 60 :  0.00021963848951275246
cost at iteration 61 :  0.00021442667992450508
cost at iteration 62 :  0.00020921766927431128
cost at iteration 63 :  0.000204013760419667
cost at iteration 64 :  0.0001988172892443127
cost at iteration 65 :  0.000193630624481986
cost at iteration 66 :  0.0001884561675351279
cost at iteration 67 :  0.0001832963522879963
cost at iteration 68 :  0.0001781536449136035
cost at iteration 69 :  0.0001730305436738912
cost at iteration 70 :  0.00016792957871249737
cost at iteration 71 :  0.0001628533118394438
cost at iteration 72 :  0.0001578043363070056
cost at iteration 73 :  0.00015278527657597594
cost at iteration 74 :  0.00014779878807144148
cost at iteration 75 :  0.0001428475569271144
cost at iteration 76 :  0.00013793429971713816
cost at iteration 77 :  0.00013306176317416345
cost at iteration 78 :  0.00012823272389233774
cost at iteration 79 :  0.00012344998801365307
cost at iteration 80 :  0.00011871639089589221
cost at iteration 81 :  0.00011403479676013763
cost at iteration 82 :  0.00010940809831549105
cost at iteration 83 :  0.00010483921635827813
cost at iteration 84 :  0.00010033109934254142
cost at iteration 85 :  9.588672291807178e-05
cost at iteration 86 :  9.150908943154551e-05
cost at iteration 87 :  8.72012273855055e-05
cost at iteration 88 :  8.296619084889476e-05
cost at iteration 89 :  7.880705881159023e-05
cost at iteration 90 :  7.472693447380333e-05
cost at iteration 91 :  7.072894445924724e-05
cost at iteration 92 :  6.681623793848252e-05
cost at iteration 93 :  6.299198564569548e-05
cost at iteration 94 :  5.9259378768134326e-05
cost at iteration 95 :  5.5621627682217944e-05
cost at iteration 96 :  5.2081960503566696e-05
cost at iteration 97 :  4.8643621409317566e-05
cost at iteration 98 :  4.5309868679294624e-05
cost at iteration 99 :  4.20839723868138e-05
cost at iteration 100 :  3.8969211648465605e-05
cost at iteration 101 :  3.5968871312813714e-05
cost at iteration 102 :  3.3086237926996486e-05
cost at iteration 103 :  3.0324594762422546e-05
cost at iteration 104 :  2.7687215597803575e-05
cost at iteration 105 :  2.5177356836717235e-05
cost at iteration 106 :  2.279824735680523e-05
cost at iteration 107 :  2.055307521423432e-05
cost at iteration 108 :  1.8444969902139046e-05
cost at iteration 109 :  1.647697818480397e-05
cost at iteration 110 :  1.4652030420463887e-05
cost at iteration 111 :  1.2972892412644388e-05
cost at iteration 112 :  1.1442094557637516e-05
cost at iteration 113 :  1.0061824122714005e-05
cost at iteration 114 :  8.833755295603537e-06
cost at iteration 115 :  7.758769642641403e-06
cost at iteration 116 :  6.836474728729002e-06
cost at iteration 117 :  6.0643353642044065e-06
cost at iteration 118 :  5.436047261333094e-06
cost at iteration 119 :  4.938530895039252e-06
cost at iteration 120 :  4.547470864356548e-06
cost at iteration 121 :  4.227139745001233e-06
cost at iteration 122 :  3.9472335639519955e-06
cost at iteration 123 :  3.6925736491469107e-06
cost at iteration 124 :  3.453047267773316e-06
cost at iteration 125 :  3.2240897549672965e-06
cost at iteration 126 :  3.004437836553187e-06
cost at iteration 127 :  2.793742499244892e-06
cost at iteration 128 :  2.591869178683128e-06
cost at iteration 129 :  2.3987316223880632e-06
cost at iteration 130 :  2.2142561206045403e-06
cost at iteration 131 :  2.0383749548855097e-06
cost at iteration 132 :  1.8710252752736256e-06
cost at iteration 133 :  1.7121486793631871e-06
cost at iteration 134 :  1.5616907492087028e-06
cost at iteration 135 :  1.4196004394998868e-06
cost at iteration 136 :  1.2858292478990998e-06
cost at iteration 137 :  1.1603300894808186e-06
cost at iteration 138 :  1.0430557254952852e-06
cost at iteration 139 :  9.339565823154416e-07
cost at iteration 140 :  8.329775913147724e-07
cost at iteration 141 :  7.400537800407698e-07
cost at iteration 142 :  6.551035379830806e-07
cost at iteration 143 :  5.780200281288589e-07
cost at iteration 144 :  5.086559860558595e-07
cost at iteration 145 :  4.4682641048577857e-07
cost at iteration 146 :  3.9253156207588683e-07
cost at iteration 147 :  3.5309272300005164e-07
cost at iteration 148 :  3.4404388882759574e-07
cost at iteration 149 :  3.3158977508157525e-07
cost at iteration 150 :  3.2165404798417593e-07
cost at iteration 151 :  3.106244637076834e-07
cost at iteration 152 :  3.0189490694761225e-07
cost at iteration 153 :  2.9205932192021226e-07
cost at iteration 154 :  2.842737469187099e-07
cost at iteration 155 :  2.754490342916045e-07
cost at iteration 156 :  2.6844078039448904e-07
cost at iteration 157 :  2.6047795170301606e-07
cost at iteration 158 :  2.5412739605791354e-07
cost at iteration 159 :  2.469036509365498e-07
cost at iteration 160 :  2.41118768994629e-07
cost at iteration 161 :  2.345321103381608e-07
cost at iteration 162 :  2.292392155540837e-07
cost at iteration 163 :  2.2320448889999354e-07
cost at iteration 164 :  2.1834315346304976e-07
cost at iteration 165 :  2.1278897447152071e-07
cost at iteration 166 :  2.0830889771670347e-07
cost at iteration 167 :  2.0317520033573583e-07
cost at iteration 168 :  1.9903408964255352e-07
cost at iteration 169 :  1.9427014157173037e-07
cost at iteration 170 :  1.9043217567198546e-07
cost at iteration 171 :  1.859949605060047e-07
cost at iteration 172 :  1.8242962002962153e-07
cost at iteration 173 :  1.7828251332910577e-07
cost at iteration 174 :  1.749636617564355e-07
cost at iteration 175 :  1.7107534526722799e-07
cost at iteration 176 :  1.6798049148532162e-07
cost at iteration 177 :  1.643240613635082e-07
cost at iteration 178 :  1.6143376002282808e-07
cost at iteration 179 :  1.579859935562003e-07
cost at iteration 180 :  1.5528335332060053e-07
cost at iteration 181 :  1.5202410534372859e-07
cost at iteration 182 :  1.4949438339956994e-07
cost at iteration 183 :  1.4640608893901527e-07
cost at iteration 184 :  1.4403635540600838e-07
cost at iteration 185 :  1.4110361939799162e-07
cost at iteration 186 :  1.3888247887497958e-07
cost at iteration 187 :  1.360917372933406e-07
cost at iteration 188 :  1.3400909735588055e-07
cost at iteration 189 :  1.3134833693600504e-07
cost at iteration 190 :  1.2939521534644665e-07
cost at iteration 191 :  1.268537414155168e-07
cost at iteration 192 :  1.2502210531787005e-07
cost at iteration 193 :  1.2259034914166137e-07
cost at iteration 194 :  1.2087298071665722e-07
cost at iteration 195 :  1.18542339326734e-07
cost at iteration 196 :  1.1693272335505479e-07
cost at iteration 197 :  1.1469542609144674e-07
cost at iteration 198 :  1.1318765566728882e-07
cost at iteration 199 :  1.1103665271197933e-07
cost at iteration 200 :  1.0962535000166442e-07
cost at iteration 201 :  1.0755421902984925e-07
cost at iteration 202 :  1.0623446850834506e-07
cost at iteration 203 :  1.042373362815782e-07
cost at iteration 204 :  1.0300462832447782e-07
cost at iteration 205 :  1.0107610462068537e-07
cost at iteration 206 :  9.992628769701268e-08
cost at iteration 207 :  9.806140943978624e-08
cost at iteration 208 :  9.699064945526976e-08
cost at iteration 209 :  9.518483328770217e-08
cost at iteration 210 :  9.418957887984535e-08
cost at iteration 211 :  9.243858074170695e-08
cost at iteration 212 :  9.151553353608065e-08
cost at iteration 213 :  8.981541406000599e-08
cost at iteration 214 :  8.896150306915097e-08
cost at iteration 215 :  8.730859782191078e-08
cost at iteration 216 :  8.652095731025066e-08
cost at iteration 217 :  8.491185107764211e-08
cost at iteration 218 :  8.418780133304511e-08
cost at iteration 219 :  8.261930578824681e-08
cost at iteration 220 :  8.195633633756642e-08
cost at iteration 221 :  8.042547054866628e-08
cost at iteration 222 :  7.982122543433992e-08
cost at iteration 223 :  7.832519876176786e-08
cost at iteration 224 :  7.777746356228528e-08
cost at iteration 225 :  7.631366057483746e-08
cost at iteration 226 :  7.582035090604343e-08
cost at iteration 227 :  7.438631800814476e-08
cost at iteration 228 :  7.394546928698072e-08
cost at iteration 229 :  7.253890280238661e-08
cost at iteration 230 :  7.214866109143434e-08
cost at iteration 231 :  7.076739659173193e-08
cost at iteration 232 :  7.042601037322538e-08
cost at iteration 233 :  6.906801307498906e-08
cost at iteration 234 :  6.877382582790985e-08
cost at iteration 235 :  6.743718191156691e-08
cost at iteration 236 :  6.718862538604912e-08
cost at iteration 237 :  6.587153411356876e-08
cost at iteration 238 :  6.566712221381323e-08
cost at iteration 239 :  6.436788874216229e-08
cost at iteration 240 :  6.420621194307164e-08
cost at iteration 241 :  6.29232407467407e-08
cost at iteration 242 :  6.280296098109233e-08
cost at iteration 243 :  6.153474981053825e-08
cost at iteration 244 :  6.14545957730902e-08
cost at iteration 245 :  6.019973008711453e-08
cost at iteration 246 :  6.015849290998986e-08
cost at iteration 247 :  5.8915640729374976e-08
cost at iteration 248 :  5.891216998966643e-08
cost at iteration 249 :  5.768007712708817e-08
scale =  53.691880265003306
scale =  26.845940132501653
scale =  13.422970066250826
scale =  6.711485033125413
scale =  3.3557425165627066
scale =  1.6778712582813533
scale =  0.8389356291406767
scale =  0.4194678145703383
scale =  0.20973390728516916
scale =  0.10486695364258458
scale =  0.05243347682129229
scale =  0.026216738410646145
scale =  0.013108369205323073
scale =  0.006554184602661536
scale =  0.003277092301330768
scale =  0.001638546150665384
scale =  0.000819273075332692
scale =  0.000409636537666346
scale =  0.000204818268833173
scale =  0.0001024091344165865
scale =  5.120456720829325e-05
scale =  2.5602283604146626e-05
scale =  1.2801141802073313e-05
scale =  6.400570901036657e-06
scale =  3.2002854505183283e-06
scale =  1.6001427252591641e-06
scale =  8.000713626295821e-07
scale =  4.0003568131479104e-07
scale =  2.0001784065739552e-07
scale =  1.0000892032869776e-07
scale =  5.000446016434888e-08
scale =  2.500223008217444e-08
scale =  1.250111504108722e-08
scale =  6.25055752054361e-09
scale =  3.125278760271805e-09
scale =  1.5626393801359025e-09
scale =  7.813196900679512e-10
scale =  3.906598450339756e-10
scale =  1.953299225169878e-10
scale =  9.76649612584939e-11
scale =  4.883248062924695e-11
scale =  2.4416240314623476e-11
scale =  1.2208120157311738e-11
scale =  6.104060078655869e-12
scale =  3.0520300393279345e-12
scale =  1.5260150196639673e-12
scale =  7.630075098319836e-13
No more descent can be found