# 7.7 PDE-Constrained Topology OptimizationΒΆ

We are again interested in minimizing a shape function of the type

$\mathcal J(\Omega) = \int_{\mathsf{D}} (u-u_d)^2 \; dx$

subject to $$u:\mathsf{D} \to \mathbb{R}$$ solves

$\int_{\mathsf{D}} \beta_\Omega\nabla u\cdot \nabla \varphi\;dx = \int_{\mathsf{D}} f_\Omega\varphi \quad \text{ for all } \varphi \in H^1_0(\mathsf{D}),$

where

• $$\beta_\Omega = \beta_1\chi_\Omega + \beta_2 \chi_{\mathsf{D}\setminus \Omega}$$, $$\qquad \beta_1,\beta_2>0$$,

• $$f_\Omega = f_1\chi_\Omega + f_2 \chi_{\mathsf{D}\setminus \Omega}$$, $$\qquad f_1,f_2\in \mathbb{R}$$,

• $$u_d:\mathsf{D} \to \mathbb{R}$$.

[1]:

from netgen.geom2d import SplineGeometry  # SplieGeometry to define a 2D mesh

from ngsolve import * # import everything from the ngsolve module
from ngsolve.webgui import Draw

import numpy as np

from interpolations import InterpolateLevelSetToElems # function which interpolates a levelset function

myMAXH = 0.05
EPS = myMAXH * 1e-6      #variable defining when a point is on the interface and when not


The file interpolations.py can contains a routine to compute the volume ratio of each element with respect to a given level set function. The file can be downloaded from interpolations.py .

We begin with defining the domain $$\mathsf{D}$$. For this we define a rectangular mesh with side length $$R$$ using the NGSolve function SplineGeometry() with maximum triangle size maxh. The boundary name of $$\partial\mathsf{D}$$ is "rectangle".

[2]:

geo = SplineGeometry()

R = 2

## add a rectangle
p2=(R,R),
bc="rectangle",
leftdomain=1,
rightdomain=0)

geo.SetMaterial (1, "outer") # give the domain the name "outer"

mesh = Mesh(geo.GenerateMesh(maxh=myMAXH)) # generate ngsolve mesh
Draw(mesh)

[2]:

BaseWebGuiScene

• fes_state,fes_adj - $$H^1$$ conforming finite element space of order 1

• pwc - $$L^2$$ subspace of functions constant on each element of mesh

• $$u,v$$ are our trial and test functions

• gfu, gfp are Gridfunctions storing the state and adjoint variable, respectively

• gfud stores $$u_d$$ - the target function

[3]:

# H1-conforming finite element space

fes_state = H1(mesh, order=1, dirichlet="rectangle")
fes_adj = H1(mesh, order=1, dirichlet="rectangle")
fes_level = H1(mesh, order=1)

pwc = L2(mesh)   #piecewise constant space

## test and trial functions
u, v = fes_state.TnT()

p, q = fes_adj.TnT()

gfu = GridFunction(fes_state)

gfud = GridFunction(fes_state)


We represent the design $$\Omega \subset \mathsf{D}$$ by means of a level set function $$\psi: \mathsf{D} \rightarrow \mathbb R$$ in the following way:

$\begin{split}\begin{array}{rl} \psi(x) < 0 &\quad \Longleftrightarrow\quad x \in \Omega \\ \psi(x) = 0 &\quad \Longleftrightarrow\quad x \in \partial \Omega \\ \psi(x) > 0 &\quad \Longleftrightarrow\quad x \in \mathsf{D} \setminus \overline \Omega. \end{array}\end{split}$
• psi - gridfunction in fes_state defining $$\psi$$

• psides - gridfunction describing the optimal shape

• psinew - dummy function for line search

[4]:

psi = GridFunction(fes_level)
psi.Set(1)

psides = GridFunction(fes_level)
psinew = GridFunction(fes_level)


Next we solve the state equation $$\int_{\mathsf{D}} \beta_\Omega \nabla u \cdot \nabla \varphi \;dx = \int_{\mathsf{D}} f_\Omega\varphi \;dx \quad \text{ for all } \varphi \in H^1_0(\mathsf{D}).$$

• beta, f_rhs - piecewise constant gridfunction

• B - bilinear form defining $$\int_{\mathsf{D}} \beta_\Omega\nabla \psi \cdot \nabla \varphi \;dx$$

• B_adj - bilinear form defining $$\int_{\mathsf{D}} \beta_\Omega\nabla \psi \cdot \nabla \varphi \;dx$$

• L - right hand side $$\int_{\mathsf{D}} f_\Omega \varphi \;dx$$

[5]:

# constants for f_rhs and beta
f1 = 10
f2 = 1

beta1 = 10
beta2 = 1

# piecewise constant coefficient functions beta and f_rhs

beta = GridFunction(pwc)
beta.Set(beta1)

f_rhs = GridFunction(pwc)
f_rhs.Set(f1)

# bilinear form for state equation
B = BilinearForm(fes_state)
B += beta*grad(u) * grad(v) * dx

L = LinearForm(fes_state)
L += f_rhs * v *dx

# solving
psides.Set(1)
InterpolateLevelSetToElems(psides, beta1, beta2, beta, mesh, EPS)
InterpolateLevelSetToElems(psides, f1, f2, f_rhs, mesh, EPS)

B.Assemble()
L.Assemble()

inv = B.mat.Inverse(fes_state.FreeDofs(), inverse="sparsecholesky") # inverse of bilinear form
gfu.vec.data = inv*L.vec

scene_u = Draw(gfu)


As an optimal shape we define $$\Omega_{opt} := \{f<0\}$$, where $$f$$ is the function from the levelset example describing a clover shape. We construct $$u_d$$ by solving

$\int_{\mathsf{D}} \beta_{\Omega_{opt}} \nabla u_d \cdot \nabla \varphi \;dx = \int_{\mathsf{D}} f_{\Omega_{opt}}\varphi \;dx \quad \text{ for all } \varphi \in H^1_0(\mathsf{D}).$
[6]:

a = 4.0/5.0
b = 2
f = CoefficientFunction( 0.1*( (sqrt((x - a)**2 + b * y**2) - 1) \
* (sqrt((x + a)**2 + b * y**2) - 1) \
* (sqrt(b * x**2 + (y - a)**2) - 1) \
* (sqrt(b * x**2 + (y + a)**2) - 1) - 0.001) )

psides.Set(f)
InterpolateLevelSetToElems(psides, beta1, beta2, beta, mesh, EPS)
InterpolateLevelSetToElems(psides, f1, f2, f_rhs, mesh, EPS)

B.Assemble()
L.Assemble()
inv = B.mat.Inverse(fes_state.FreeDofs(), inverse="sparsecholesky") # inverse of bilinear form

gfud.vec.data = inv*L.vec

Draw(gfud, mesh, 'gfud')

[6]:

BaseWebGuiScene


The function SolvePDE solves the state and adjoint state each time it is called. Note that since the functions beta and f_rhs are GridFunctions we do not have to re-define the bilinear form and linear form when beta and f_rhs change.

[7]:

def SolvePDE(adjoint=False):
#Newton(a, gfu, printing = False, maxerr = 3e-9)

B.Assemble()
L.Assemble()

inv_state = B.mat.Inverse(fes_state.FreeDofs(), inverse="sparsecholesky")

# solve state equation
gfu.vec.data = inv_state*L.vec

if adjoint == True:
# solve adjoint state equatoin
duCost.Assemble()
gfp.vec.data = -inv_adj * duCost.vec
scene_u.Redraw()

InterpolateLevelSetToElems(psi, beta1, beta2, beta, mesh, EPS)
InterpolateLevelSetToElems(psi, f1, f2, f_rhs, mesh, EPS)

SolvePDE()

# define the cost function
def Cost(u):
return (u - gfud)**2*dx

# derivative of cost function
duCost += 2*(gfu-gfud) * q * dx

print("initial cost = ", Integrate(Cost(gfu), mesh))

initial cost =  39.29567644495182


It remains to implement the topological derivative: In define $$NegPos$$ and $$PosNeg$$ as the topological derivative in $$\Omega$$ respectively $$D \setminus \overline{\Omega}$$:

$\begin{split}\begin{array}{rl} \mathsf{NegPos}(x) := -DJ(\Omega,\omega)(x) &= (-1) * \left( 2 \beta_2 \left(\frac{\beta_2-\beta_1}{\beta_1+\beta_2}\right)\nabla u(x)\cdot \nabla p(x) - (f_2-f_1)p(x) \right) \\ &= 2 \beta_2 \left(\frac{\beta_1-\beta_2}{\beta_1+\beta_2}\right)\nabla u(x)\cdot \nabla p(x) - (f_1-f_2)p(x) \quad \text{ for }x\in \Omega \\ \mathsf{PosNeg}(x) := DJ(\Omega,\omega)(x) &= 2 \beta_1 \left(\frac{\beta_1-\beta_2}{\beta_1+\beta_2}\right)\nabla u(x)\cdot \nabla p(x) - (f_1-f_2)p(x) \quad \text{ for }x\in D \setminus \overline{\Omega} \end{array}\end{split}$

and the update function (also called generalised topological derivative)

$g_\Omega(x):= -\chi_\Omega (x)DJ(\Omega,\omega)(x) + (1-\chi_\Omega)*DJ(\Omega,\omega)(x) .$
[8]:

BetaPosNeg = 2 * beta2 * (beta1-beta2)/(beta1+beta2)    # factors in TD in positive part {x: \psi(x)>0} ={x: \beta(x) = \beta_2}
BetaNegPos = 2 * beta1 * (beta1-beta2)/(beta1+beta2)    # factors in TD in positive part {x: \psi(x)>0} ={x: \beta(x) = \beta_2}
FPosNeg = -(f1-f2)

psinew.vec.data = psi.vec

## normalise psi in L2
normPsi = sqrt(Integrate(psi**2*dx, mesh))
psi.vec.data = 1/normPsi * psi.vec

kappa = 0.05

# set level set function to data
InterpolateLevelSetToElems(psi, beta1, beta2, beta, mesh, EPS)
InterpolateLevelSetToElems(psi, f1, f2, f_rhs, mesh, EPS)

Redraw()

TD_node = GridFunction(fes_level)
#scene1 = Draw(TD_node, mesh, "TD_node")

TD_pwc = GridFunction(pwc)
#scene2 = Draw(TD_pwc, mesh, "TD_pwc")

TDPosNeg_pwc = GridFunction(pwc)
TDNegPos_pwc = GridFunction(pwc)

cutRatio = GridFunction(pwc)

# solve for current configuration

Finally we can define a loop defining the levelset update: :nbsphinx-math:begin{align*}

psi_{k+1} = (1-s_k) frac{psi_k}{|\psi_k|} + s_k frac{g_{psi_k}}{|g_{\psi_k}|}

end{align*}

[9]:

scene = Draw(psi)
scene_cr = Draw(cutRatio, mesh, "cutRatio")

[10]:

iter_max = 20
converged = False

xm=0.
ym=0.
psi.Set( (x-xm)**2+(y-ym)**2-0.25**2)
psinew.vec.data= psi.vec
scene.Redraw()

J = Integrate(Cost(gfu),mesh)

for k in range(iter_max):
print("================ iteration ", k, "===================")

# copy new levelset data from psinew into psi
psi.vec.data = psinew.vec
scene.Redraw()

J_current = Integrate(Cost(gfu),mesh)

print( Integrate( (gfu-gfud)**2*dx, mesh) )

print("cost in beginning of iteration", k, ": Cost = ", J_current)

# compute the piecewise constant topological derivative in each domain
TDPosNeg_pwc.Set( BetaPosNeg * grad(gfu) * grad(gfp)  + FPosNeg*gfp )
TDNegPos_pwc.Set( BetaNegPos * grad(gfu) * grad(gfp)  + FPosNeg*gfp )

# compute the cut ratio of the interface elements
InterpolateLevelSetToElems(psi, 1, 0, cutRatio, mesh, EPS)
scene_cr.Redraw()

# compute the combined topological derivative using the cut ratio information
for j in range(len(TD_pwc.vec)):
TD_pwc.vec[j] = cutRatio.vec[j] * TDNegPos_pwc.vec[j] + (1-cutRatio.vec[j])*TDPosNeg_pwc.vec[j]

TD_node.Set(TD_pwc)

normTD = sqrt(Integrate(TD_node**2*dx, mesh)) # L2 norm of TD_node
TD_node.vec.data = 1/normTD * TD_node.vec # normalised TD_node

normPsi = sqrt(Integrate(psi**2*dx, mesh)) # L2 norm of psi
psi.vec.data = 1/normPsi * psi.vec  # normalised psi

linesearch = True

for j in range(10):

# update the level set function
psinew.vec.data = (1-kappa)*psi.vec + kappa*TD_node.vec

# update beta and f_rhs
InterpolateLevelSetToElems(psinew, beta1, beta2, beta, mesh, EPS)
InterpolateLevelSetToElems(psinew, f1, f2, f_rhs, mesh, EPS)

# solve PDE without adjoint
SolvePDE()

Redraw(blocking=True)

Jnew = Integrate(Cost(gfu), mesh)

if Jnew > J_current:
print("--------------------")
print("-----line search ---")
print("--------------------")
kappa = kappa*0.8
print("kappa", kappa)
else:
break

Redraw(blocking=True)
print("----------- Jnew in  iteration ", k, " = ", Jnew, " (kappa = ", kappa, ")")
print('')
print("iter" + str(k) + ", Jnew = " + str(Jnew) + " (kappa = ", kappa, ")")
kappa = min(1, kappa*1.2)

print("end of iter " + str(k))

================ iteration  0 ===================
39.29567644495183
cost in beginning of iteration 0 : Cost =  39.29567644495183
----------- Jnew in  iteration  0  =  11.859500825448178  (kappa =  0.05 )

iter0, Jnew = 11.859500825448178 (kappa =  0.05 )
end of iter 0
================ iteration  1 ===================
11.859500825447913
cost in beginning of iteration 1 : Cost =  11.859500825447913
----------- Jnew in  iteration  1  =  6.3662170311137425  (kappa =  0.06 )

iter1, Jnew = 6.3662170311137425 (kappa =  0.06 )
end of iter 1
================ iteration  2 ===================
6.366217031113691
cost in beginning of iteration 2 : Cost =  6.366217031113691
----------- Jnew in  iteration  2  =  6.219404002334578  (kappa =  0.072 )

iter2, Jnew = 6.219404002334578 (kappa =  0.072 )
end of iter 2
================ iteration  3 ===================
6.219404002334453
cost in beginning of iteration 3 : Cost =  6.219404002334453
--------------------
-----line search ---
--------------------
kappa 0.06912
--------------------
-----line search ---
--------------------
kappa 0.055296000000000005
--------------------
-----line search ---
--------------------
kappa 0.04423680000000001
--------------------
-----line search ---
--------------------
kappa 0.03538944000000001
--------------------
-----line search ---
--------------------
kappa 0.028311552000000007
--------------------
-----line search ---
--------------------
kappa 0.02264924160000001
--------------------
-----line search ---
--------------------
kappa 0.018119393280000007
--------------------
-----line search ---
--------------------
kappa 0.014495514624000005
----------- Jnew in  iteration  3  =  5.650460303043148  (kappa =  0.014495514624000005 )

iter3, Jnew = 5.650460303043148 (kappa =  0.014495514624000005 )
end of iter 3
================ iteration  4 ===================
5.650460303043158
cost in beginning of iteration 4 : Cost =  5.650460303043158
--------------------
-----line search ---
--------------------
kappa 0.013915694039040007
----------- Jnew in  iteration  4  =  5.538240999806068  (kappa =  0.013915694039040007 )

iter4, Jnew = 5.538240999806068 (kappa =  0.013915694039040007 )
end of iter 4
================ iteration  5 ===================
5.538240999805812
cost in beginning of iteration 5 : Cost =  5.538240999805812
--------------------
-----line search ---
--------------------
kappa 0.013359066277478408
--------------------
-----line search ---
--------------------
kappa 0.010687253021982727
--------------------
-----line search ---
--------------------
kappa 0.008549802417586181
--------------------
-----line search ---
--------------------
kappa 0.006839841934068946
----------- Jnew in  iteration  5  =  5.4770935164670815  (kappa =  0.006839841934068946 )

iter5, Jnew = 5.4770935164670815 (kappa =  0.006839841934068946 )
end of iter 5
================ iteration  6 ===================
5.47709351646705
cost in beginning of iteration 6 : Cost =  5.477093516467048
--------------------
-----line search ---
--------------------
kappa 0.006566248256706188
----------- Jnew in  iteration  6  =  5.085963722546003  (kappa =  0.006566248256706188 )

iter6, Jnew = 5.085963722546003 (kappa =  0.006566248256706188 )
end of iter 6
================ iteration  7 ===================
5.085963722546078
cost in beginning of iteration 7 : Cost =  5.085963722546078
--------------------
-----line search ---
--------------------
kappa 0.00630359832643794
--------------------
-----line search ---
--------------------
kappa 0.005042878661150352
--------------------
-----line search ---
--------------------
kappa 0.004034302928920282
----------- Jnew in  iteration  7  =  5.024676369417312  (kappa =  0.004034302928920282 )

iter7, Jnew = 5.024676369417312 (kappa =  0.004034302928920282 )
end of iter 7
================ iteration  8 ===================
5.0246763694173495
cost in beginning of iteration 8 : Cost =  5.024676369417351
--------------------
-----line search ---
--------------------
kappa 0.003872930811763471
----------- Jnew in  iteration  8  =  4.680164002928751  (kappa =  0.003872930811763471 )

iter8, Jnew = 4.680164002928751 (kappa =  0.003872930811763471 )
end of iter 8
================ iteration  9 ===================
4.6801640029288505
cost in beginning of iteration 9 : Cost =  4.6801640029288505
--------------------
-----line search ---
--------------------
kappa 0.003718013579292932
--------------------
-----line search ---
--------------------
kappa 0.0029744108634343455
--------------------
-----line search ---
--------------------
kappa 0.0023795286907474767
--------------------
-----line search ---
--------------------
kappa 0.0019036229525979814
----------- Jnew in  iteration  9  =  4.548013307504089  (kappa =  0.0019036229525979814 )

iter9, Jnew = 4.548013307504089 (kappa =  0.0019036229525979814 )
end of iter 9
================ iteration  10 ===================
4.548013307504103
cost in beginning of iteration 10 : Cost =  4.548013307504104
----------- Jnew in  iteration  10  =  4.297138500840764  (kappa =  0.0022843475431175778 )

iter10, Jnew = 4.297138500840764 (kappa =  0.0022843475431175778 )
end of iter 10
================ iteration  11 ===================
4.2971385008408545
cost in beginning of iteration 11 : Cost =  4.2971385008408545
--------------------
-----line search ---
--------------------
kappa 0.002192973641392875
--------------------
-----line search ---
--------------------
kappa 0.0017543789131143001
--------------------
-----line search ---
--------------------
kappa 0.00140350313049144
----------- Jnew in  iteration  11  =  4.240506681973666  (kappa =  0.00140350313049144 )

iter11, Jnew = 4.240506681973666 (kappa =  0.00140350313049144 )
end of iter 11
================ iteration  12 ===================
4.2405066819736446
cost in beginning of iteration 12 : Cost =  4.240506681973644
----------- Jnew in  iteration  12  =  4.116140835706864  (kappa =  0.0016842037565897282 )

iter12, Jnew = 4.116140835706864 (kappa =  0.0016842037565897282 )
end of iter 12
================ iteration  13 ===================
4.116140835706807
cost in beginning of iteration 13 : Cost =  4.116140835706807
----------- Jnew in  iteration  13  =  4.105202561884528  (kappa =  0.002021044507907674 )

iter13, Jnew = 4.105202561884528 (kappa =  0.002021044507907674 )
end of iter 13
================ iteration  14 ===================
4.105202561884544
cost in beginning of iteration 14 : Cost =  4.105202561884544
----------- Jnew in  iteration  14  =  4.030283152748079  (kappa =  0.0024252534094892086 )

iter14, Jnew = 4.030283152748079 (kappa =  0.0024252534094892086 )
end of iter 14
================ iteration  15 ===================
4.03028315274824
cost in beginning of iteration 15 : Cost =  4.03028315274824
--------------------
-----line search ---
--------------------
kappa 0.00232824327310964
--------------------
-----line search ---
--------------------
kappa 0.0018625946184877122
----------- Jnew in  iteration  15  =  4.01335836417223  (kappa =  0.0018625946184877122 )

iter15, Jnew = 4.01335836417223 (kappa =  0.0018625946184877122 )
end of iter 15
================ iteration  16 ===================
4.013358364172229
cost in beginning of iteration 16 : Cost =  4.013358364172229
----------- Jnew in  iteration  16  =  3.8845675205276464  (kappa =  0.0022351135421852545 )

iter16, Jnew = 3.8845675205276464 (kappa =  0.0022351135421852545 )
end of iter 16
================ iteration  17 ===================
3.8845675205277166
cost in beginning of iteration 17 : Cost =  3.8845675205277166
--------------------
-----line search ---
--------------------
kappa 0.0021457090004978444
--------------------
-----line search ---
--------------------
kappa 0.0017165672003982757
----------- Jnew in  iteration  17  =  3.836992769594235  (kappa =  0.0017165672003982757 )

iter17, Jnew = 3.836992769594235 (kappa =  0.0017165672003982757 )
end of iter 17
================ iteration  18 ===================
3.8369927695942376
cost in beginning of iteration 18 : Cost =  3.8369927695942376
----------- Jnew in  iteration  18  =  3.736383183809279  (kappa =  0.0020598806404779307 )

iter18, Jnew = 3.736383183809279 (kappa =  0.0020598806404779307 )
end of iter 18
================ iteration  19 ===================
3.7363831838092887
cost in beginning of iteration 19 : Cost =  3.7363831838092887
--------------------
-----line search ---
--------------------
kappa 0.0019774854148588133
----------- Jnew in  iteration  19  =  3.71813921466786  (kappa =  0.0019774854148588133 )

iter19, Jnew = 3.71813921466786 (kappa =  0.0019774854148588133 )
end of iter 19


After 999 iterations: $J =$

[ ]: