# 7.4 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|^q ; dx, quad qge 2

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)$$.

[1]:

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

[2]:

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

[3]:

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

[4]:

fes = H1(mesh, order=2, dirichlet=".*")
gfu = GridFunction(fes)
scene_u = Draw (gfu, mesh, "state")

gfp = GridFunction(fes)
scene_p = Draw (gfp, mesh, "adjoint")


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

## Automatic Shape Differentiation¶

The formula for the shape derivative was derived as the partial derivative of the perturbed Lagrangian (brought back to the original domain): :nbsphinx-math:begin{align*}
dmathcal J(Omega; X) = frac{partial}{partial t} left( int_Omega |u - u_d^t|^q ,mbox{det}(F_t) mbox{d} x
• int_{Omega} (F_t^{-top}nabla u) cdot (F_t^{-top} nabla p), mbox{det}(F_t) , dx - int_{Omega}f^t p ,mbox{det}(F_t) ,dx right)|_{t=0}

end{align*} where

• $$T_t(x)=x+tX(x)=y$$

• $$F_t = DT_t = I+t DX$$

• $$u_d^t = u_d \circ T_t$$

• $$f^t = f \circ T_t$$

The integrand depends on the parameter $$t$$ only via $$F_t$$ and $$T_t$$. We define the Lagrangian in this form, involving a parameter t that has the value 0.

[5]:

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

#define symbolic parameter t=0 in order to allow differentiation w.r.t. t
t = Parameter(0)
F = Id(2) + t * grad(X).trans

def Equation(u,v):

q=6
def CostAuto(u):
return (u-ud)**q*Det(F)

def CostAuto2(u):
return (u-ud)**q

Lagrangian = CostAuto(gfu) + Equation(gfu,gfp)


## State equation¶

Equation can also be used to define the bilinear form. The following defines left and right hand side of the PDE in a “BilinearForm”:

[6]:

u, v = fes.TnT()

aAuto = BilinearForm(fes, symmetric=True)
aAuto += Equation(u,v)*dx


Now the PDE can be conveniently solved by calling Newton’s method (which terminates after one iteration since the PDE is linear)

[7]:

gfu.vec[:]=0
Newton(aAuto,gfu,freedofs=fes.FreeDofs())
scene_u.Redraw()

Newton iteration  0
err =  0.13016095740357994
Newton iteration  1
err =  1.3386204426309707e-16


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:

[8]:

p, w = fes.TnT()


[9]:

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

[10]:

fadjoint.Assemble()
scene_p.Redraw()


## Automatic Shape Differentiation¶

Denoting the integrand by $$G^{u,p}$$, the shape derivative is given by :nbsphinx-math:begin{align*}

dmathcal J(Omega; X) =& left( left. frac{partial G^{u,p}}{partial t} + frac{d G^{u,p}}{dy} cdot frac{d T_t}{dt}right)rightrvert_{t=0} \ =& left. frac{partial G^{u,p}}{partial t}rightrvert_{t=0} + frac{d G^{u,p}}{dy} cdot X

end{align*}

[11]:

dJOmegaAuto = LinearForm(VEC)
dJOmegaAuto += Lagrangian.Diff(t) * dx
dJOmegaAuto += (Lagrangian.Diff(x, X[0]) + Lagrangian.Diff(y, X[1])) * dx

[12]:

b = BilinearForm(VEC)

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

[13]:

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

[14]:

b.Assemble()
dJOmegaAuto.Assemble()
SolveDeformationEquationAuto()
Draw(-gfX, mesh, "-gfX")

[14]:

BaseWebGuiScene

[15]:

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

Cost at initial design 6.457232118297974e-13

[16]:

mesh.SetDeformation(gfset)
gfu.vec[:]=0
Newton(aAuto, gfu, fes.FreeDofs())
print('Cost at new design', Integrate (CostAuto2(gfu), mesh))
scene.Redraw()

Newton iteration  0
err =  0.15710464112302447
Newton iteration  1
err =  1.114535141538621e-16
Cost at new design 1.5648982670652649e-13


Thus, the user has to enter the PDE (in its transformed form) only once.

Finally, let us again run the full algorithm:

[17]:

#reset to and solve for initial configuration
scene_u = Draw(gfu)
gfset.Set((0,0))
mesh.SetDeformation(gfset)
gfu.vec[:]=0
Newton(aAuto, gfu, fes.FreeDofs())

LineSearch = True

iter_max = 600
Jold = Integrate(CostAuto2(gfu), mesh)
converged = False
for k in range(iter_max):
print('cost at iteration', k, ': ', Jold)
mesh.SetDeformation(gfset)
scene_u.Redraw()

gfu.vec[:]=0
Newton(aAuto, gfu, fes.FreeDofs(), printing = False)

b.Assemble()
dJOmegaAuto.Assemble()
SolveDeformationEquationAuto()

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

Jnew = Integrate(CostAuto2(gfu), mesh)

if LineSearch:
while Jnew > Jold and scale > 1e-12:
#input('a')
scale = scale / 2

if scale <= 1e-12:
converged = True
break

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

gfu.vec[:]=0
Newton(aAuto, gfu, fes.FreeDofs(), printing = False)
Jnew = Integrate(CostAuto2(gfu), mesh)

if converged==True:
break
Jold = Jnew

Redraw(blocking=True)

Newton iteration  0
err =  0.13016095740357994
Newton iteration  1
err =  1.3386204426309707e-16
cost at iteration 0 :  6.457232118297974e-13
cost at iteration 1 :  5.983807124063546e-13
cost at iteration 2 :  5.335448115442614e-13
cost at iteration 3 :  4.751216080999641e-13
cost at iteration 4 :  4.225554869119028e-13
cost at iteration 5 :  3.7532939439192975e-13
cost at iteration 6 :  3.329631544773859e-13
cost at iteration 7 :  2.950117791617657e-13
cost at iteration 8 :  2.6106378013677934e-13
cost at iteration 9 :  2.307394873259726e-13
cost at iteration 10 :  2.0368937964401986e-13
cost at iteration 11 :  1.7959243304608907e-13
cost at iteration 12 :  1.58154490764688e-13
cost at iteration 13 :  1.3910666052461724e-13
cost at iteration 14 :  1.2220374345252617e-13
cost at iteration 15 :  1.0722269933237963e-13
cost at iteration 16 :  9.396115277268854e-14
cost at iteration 17 :  8.223594470161334e-14
cost at iteration 18 :  7.18817333213215e-14
cost at iteration 19 :  6.274964811991998e-14
cost at iteration 20 :  5.470599957837416e-14
cost at iteration 21 :  4.763104554142932e-14
cost at iteration 22 :  4.141781241702625e-14
cost at iteration 23 :  3.597096478525697e-14
cost at iteration 24 :  3.120570970421562e-14
cost at iteration 25 :  2.7046710732450874e-14
cost at iteration 26 :  2.3426969975031613e-14
cost at iteration 27 :  2.0286613350628315e-14
cost at iteration 28 :  1.7571486957012336e-14
cost at iteration 29 :  1.523145440506158e-14
cost at iteration 30 :  1.3218324534175212e-14
cost at iteration 31 :  1.1483576290842346e-14
cost at iteration 32 :  9.97680187939065e-15
cost at iteration 33 :  8.647382476478587e-15
cost at iteration 34 :  7.452880491841562e-15
cost at iteration 35 :  6.371404569777735e-15
cost at iteration 36 :  5.402186692077566e-15
cost at iteration 37 :  4.549452824120365e-15
cost at iteration 38 :  3.81066808174099e-15
cost at iteration 39 :  3.177057233076808e-15
cost at iteration 40 :  2.637317999005241e-15
cost at iteration 41 :  2.179918202719651e-15
cost at iteration 42 :  1.79405274927764e-15
cost at iteration 43 :  1.4699490993321422e-15
cost at iteration 44 :  1.19890846509175e-15
cost at iteration 45 :  9.732542202716918e-16
cost at iteration 46 :  7.862492466904384e-16
cost at iteration 47 :  6.320100337678094e-16
cost at iteration 48 :  5.054198076362634e-16
cost at iteration 49 :  4.020509582608089e-16
cost at iteration 50 :  3.1808625809796445e-16
cost at iteration 51 :  2.5025610429080207e-16
cost at iteration 52 :  1.9576444454089132e-16
cost at iteration 53 :  1.522494628760481e-16
cost at iteration 54 :  1.176992620352825e-16
cost at iteration 55 :  9.048127696733899e-17
cost at iteration 56 :  6.92514565359858e-17
cost at iteration 57 :  5.3618545468017456e-17
cost at iteration 58 :  4.420112237581565e-17
cost at iteration 59 :  3.8750581046814984e-17
cost at iteration 60 :  3.359801182111596e-17
cost at iteration 61 :  2.936433673760572e-17
cost at iteration 62 :  2.5534777457844717e-17
cost at iteration 63 :  2.2329730829672094e-17
cost at iteration 64 :  1.9455850725269072e-17
cost at iteration 65 :  1.7035596244179587e-17
cost at iteration 66 :  1.4868852810918253e-17
cost at iteration 67 :  1.3038194261638847e-17
cost at iteration 68 :  1.1398948817819407e-17
cost at iteration 69 :  1.0010433195187479e-17
cost at iteration 70 :  8.766375822114432e-18
cost at iteration 71 :  7.709984673730523e-18
cost at iteration 72 :  6.762960604775631e-18
cost at iteration 73 :  5.9567118962218816e-18
cost at iteration 74 :  5.2336035964927505e-18
cost at iteration 75 :  4.616303638134312e-18
cost at iteration 76 :  4.062470646940978e-18
cost at iteration 77 :  3.5883279509185565e-18
cost at iteration 78 :  3.1628228525339155e-18
cost at iteration 79 :  2.7974790061598753e-18
cost at iteration 80 :  2.4695359959173417e-18
cost at iteration 81 :  2.187141765752333e-18
cost at iteration 82 :  1.9335879171612602e-18
cost at iteration 83 :  1.7146446341728737e-18
cost at iteration 84 :  1.5179843170585754e-18
cost at iteration 85 :  1.3477453217895918e-18
cost at iteration 86 :  1.1947378152030797e-18
cost at iteration 87 :  1.0620203568756956e-18
cost at iteration 88 :  9.42619920508546e-19
cost at iteration 89 :  8.389181422472967e-19
cost at iteration 90 :  7.454808121900994e-19
cost at iteration 91 :  6.642993784462727e-19
cost at iteration 92 :  5.90985958686429e-19
cost at iteration 93 :  5.27335517909069e-19
cost at iteration 94 :  4.696598424943588e-19
cost at iteration 95 :  4.196744292065327e-19
cost at iteration 96 :  3.7416417298634275e-19
cost at iteration 97 :  3.3481956568899806e-19
cost at iteration 98 :  2.987729515154402e-19
cost at iteration 99 :  2.676961741158917e-19
cost at iteration 100 :  2.3901932631701643e-19
cost at iteration 101 :  2.143674308478638e-19
cost at iteration 102 :  1.9145906146325793e-19
cost at iteration 103 :  1.7182498907168402e-19
cost at iteration 104 :  1.5347161562700187e-19
cost at iteration 105 :  1.3778886094704812e-19
cost at iteration 106 :  1.2306402738370842e-19
cost at iteration 107 :  1.1051662959335243e-19
cost at iteration 108 :  9.869907499437263e-20
cost at iteration 109 :  8.865194275790789e-20
cost at iteration 110 :  7.916948602902507e-20
cost at iteration 111 :  7.112113509773928e-20
cost at iteration 112 :  6.351466391228181e-20
cost at iteration 113 :  5.706591013342218e-20
cost at iteration 114 :  5.096541175136101e-20
cost at iteration 115 :  4.579752219795809e-20
cost at iteration 116 :  4.090476189550756e-20
cost at iteration 117 :  3.676309228910413e-20
cost at iteration 118 :  3.283816001892287e-20
cost at iteration 119 :  2.9519258267184776e-20
cost at iteration 120 :  2.6369651755044577e-20
cost at iteration 121 :  2.3710859558916776e-20
cost at iteration 122 :  2.118245651360465e-20
cost at iteration 123 :  1.9053531223710627e-20
cost at iteration 124 :  1.702311719373502e-20
cost at iteration 125 :  1.5319576208324058e-20
cost at iteration 126 :  1.368871390655937e-20
cost at iteration 127 :  1.2326625936351742e-20
cost at iteration 128 :  1.1016671614051386e-20
cost at iteration 129 :  9.928597613373438e-21
cost at iteration 130 :  8.876668543488613e-21
cost at iteration 131 :  8.008463698864324e-21
cost at iteration 132 :  7.164182349171988e-21
cost at iteration 133 :  6.47242592758419e-21
cost at iteration 134 :  5.795307483456472e-21
cost at iteration 135 :  5.245202240165372e-21
cost at iteration 136 :  4.702584956631982e-21
cost at iteration 137 :  4.266228090974231e-21
cost at iteration 138 :  3.831655821387504e-21
cost at iteration 139 :  3.4866186367089994e-21
cost at iteration 140 :  3.138584503269148e-21
cost at iteration 141 :  2.8667588806397438e-21
cost at iteration 142 :  2.5877293775522585e-21
cost at iteration 143 :  2.374414102854748e-21
cost at iteration 144 :  2.1500764566304243e-21
cost at iteration 145 :  1.983280157333183e-21
cost at iteration 146 :  1.80197549978233e-21
cost at iteration 147 :  1.6719021565981528e-21
cost at iteration 148 :  1.5242098902184948e-21
cost at iteration 149 :  1.4228792296677806e-21
cost at iteration 150 :  1.3013093903829178e-21
cost at iteration 151 :  1.2222632373813545e-21
cost at iteration 152 :  1.1210093861548509e-21
cost at iteration 153 :  1.0590655932626523e-21
cost at iteration 154 :  9.737821869538259e-22
cost at iteration 155 :  9.248039556722975e-22
cost at iteration 156 :  8.524048119782883e-22
cost at iteration 157 :  8.130498591274212e-22
cost at iteration 158 :  7.51585573047637e-22
cost at iteration 159 :  7.189967316259523e-22
cost at iteration 160 :  6.677866377916812e-22
cost at iteration 161 :  6.391931349253566e-22
cost at iteration 162 :  5.995190482421353e-22
cost at iteration 163 :  5.716750926230614e-22
cost at iteration 164 :  5.45666418195522e-22
cost at iteration 165 :  5.144532349502406e-22
cost at iteration 166 :  4.983804310665939e-22
cost at iteration 167 :  4.652395122518961e-22
cost at iteration 168 :  4.551049856044431e-22
cost at iteration 169 :  4.2240040563249166e-22
cost at iteration 170 :  4.1673280139577306e-22
cost at iteration 171 :  3.8451078411541125e-22
cost at iteration 172 :  3.8319620020642633e-22
cost at iteration 173 :  3.50927637134321e-22

[ ]: