This page was generated from unit-7-optimization/04_Topological_Derivative_Levelset.ipynb.
7.6 Topology Optimization based on the Topological DerivativeΒΆ
- We are again interested in minimizing a shape function of the type :nbsphinx-math:`begin{align*}
mathcal J(Omega) = int_Omega f(x) , mbox dx,qquad f:mathbb{R}^dto mathbb{R}
end{align*}` over all (admissible) shapes \(\Omega \subset \mathsf{D}\) with \(\mathsf{D}\) a given hold-all domain and \(f \in C(\overline{\mathsf{D}})\) a given function.
[1]:
# module to generate the 2D geometry
from netgen.geom2d import SplineGeometry
from ngsolve import *
from ngsolve.webgui import Draw
import numpy as np
from interpolations import InterpolateLevelSetToElems
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:
- :nbsphinx-math:`begin{align*}
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{align*}`
We define the domain \(\mathsf{D} = [-2, 2]^2\):
[2]:
myMAXH = 0.05
geo = SplineGeometry()
R = 2
## add a rectangle
geo.AddRectangle(p1=(-R,-R),
p2=(R,R),
bc="rectangle",
leftdomain=1,
rightdomain=0)
mesh = Mesh(geo.GenerateMesh(maxh=myMAXH))
Draw(mesh)
[2]:
BaseWebGuiScene
We represent the level set function as a GridFunction on the mesh of \(\mathsf{D}\).
[3]:
# H1-conforming finite element space
fes = H1(mesh, order=1) # C0 conforming finite element space
pwc = L2(mesh) # piecewise constant (on each element) space
psi = GridFunction(fes) # level set function
psi.vec[:]=1
psi_new = GridFunction(fes) # grid function for line search
a = 4.0/5.0
b = 2
e = 0.001
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) - e) )
Draw(f, mesh)
EPS = myMAXH * 1e-6 #variable defining when a point is on the interface and when not
[4]:
def Cost(psi):
return IfPos(psi,0,f)*dx
print("initial cost = ", Integrate(Cost(psi), mesh))
initial cost = 0.0
The evolution of the level set function is guided by the \(\textbf{topological derivative}\). We define and draw the topological derivative in \(\Omega\), \(\mathsf{D}\setminus \Omega\) and the generalized topological derivative \(g\).
[5]:
TDPosNeg_pwc = GridFunction(pwc)
TDNegPos_pwc = GridFunction(pwc)
TD_pwc = GridFunction(pwc)
CutRatio = GridFunction(pwc)
TD_node = GridFunction(fes)
scene1 = Draw(TD_node, mesh, "TD_node")
scene2 = Draw(TD_pwc, mesh, "TD_pwc")
scene3 = Draw(CutRatio, mesh, "CutRatio")
kappaMax = 1
kappa = 0.1*kappaMax
Redraw()
- In define \(NegPos\) and \(PosNeg\) as the topological derivative in \(\Omega\) respectively \(\mathsf{D}\setminus\overline{\Omega}\): :nbsphinx-math:`begin{align*}
mathsf{NegPos}(x) := -DJ(Omega,omega)(x) &= f(x) quad text{ for }xin Omega \ mathsf{PosNeg}(x) := DJ(Omega,omega)(x) &= f(x) quad text{ for }xin overline{mathsf{D}}setminusOmega^c
end{align*}` and the level set update function (also called generalised topological derivative) \begin{equation} g_\Omega(x) := \chi_\Omega (x)\mathsf{NegPos}(x) + (1-\chi_\Omega(x))\mathsf{PosNeg}(x) = f(x). \end{equation}
[6]:
TDNegPos_pwc.Set(f)
TDPosNeg_pwc.Set(f)
InterpolateLevelSetToElems(psi, 1, 0, CutRatio, mesh, EPS)
# 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
scene1.Redraw()
scene2.Redraw()
scene3.Redraw()
- The idea now is to make a fixed point iteration for \(\ell\ge 0\) :nbsphinx-math:`begin{align*}
psi_{ell+1} = (1-kappa) psi_ell + kappa frac{g_{psi_ell}}{|g_{psi_ell}|_{L_2(mathsf{D})}}, qquad kappain [0,1]
end{align*}` Let us do one iteration:
[7]:
scene = Draw(psi, mesh)
kappaTest = 0.8
psi.vec.data = (1-kappaTest)*psi.vec + kappaTest * TD_node.vec
scene.Redraw()
Now, let us repeat the procedure in the following iterative algorithm
[8]:
iter_max = 20
scene_psi = Draw(psi)
psi.Set(1)
scene_psi.Redraw()
psi_new.vec.data = psi.vec
kappa = 0.1*kappaMax
# input("Press enter to start optimization")
for k in range(iter_max):
print("================ iteration ", k, "===================")
# copy new levelset data from psi_new into psi
psi.vec.data = psi_new.vec
scene_psi.Redraw()
# current cost function value
J_current = Integrate(Cost(psi),mesh)
print("cost in beginning of iteration", k, ": Cost = ", J_current)
# level set update function
TDPosNeg_pwc.Set(f)
TDNegPos_pwc.Set(f)
# get fraction of intersecting elements (for psi<0)
InterpolateLevelSetToElems(psi, 1, 0, CutRatio, mesh, EPS)
# weight the topological derivative with CutRatio
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
psi.vec.data = (1-kappa)*psi.vec + kappa * TD_node.vec
for j in range(10):
# update the level set function
psi_new.vec.data = (1-kappa)*psi.vec + kappa * TD_node.vec
Jnew = Integrate(Cost(psi_new), mesh)
Redraw()
print("----------- Jnew in iteration ", k, " = ", Jnew, "(kappa = ", kappa, ")")
print('')
if Jnew <= J_current+1e-10:
print("----------------------------")
print("======= step accepted ======")
print("----------------------------")
kappa = min(1.1*kappa, kappaMax)
break
else:
kappa = kappa*0.8
Redraw()
print("iter" + str(k) + ", Jnew = " + str(Jnew) + "(kappa = ", kappa, ")")
print("end of iter " + str(k))
================ iteration 0 ===================
cost in beginning of iteration 0 : Cost = 0.0
----------- Jnew in iteration 0 = 0.0 (kappa = 0.1 )
----------------------------
======= step accepted ======
----------------------------
iter0, Jnew = 0.0(kappa = 0.11000000000000001 )
end of iter 0
================ iteration 1 ===================
cost in beginning of iteration 1 : Cost = 0.0
----------- Jnew in iteration 1 = 0.0 (kappa = 0.11000000000000001 )
----------------------------
======= step accepted ======
----------------------------
iter1, Jnew = 0.0(kappa = 0.12100000000000002 )
end of iter 1
================ iteration 2 ===================
cost in beginning of iteration 2 : Cost = 0.0
----------- Jnew in iteration 2 = 0.0 (kappa = 0.12100000000000002 )
----------------------------
======= step accepted ======
----------------------------
iter2, Jnew = 0.0(kappa = 0.13310000000000002 )
end of iter 2
================ iteration 3 ===================
cost in beginning of iteration 3 : Cost = 0.0
----------- Jnew in iteration 3 = 0.0 (kappa = 0.13310000000000002 )
----------------------------
======= step accepted ======
----------------------------
iter3, Jnew = 0.0(kappa = 0.14641000000000004 )
end of iter 3
================ iteration 4 ===================
cost in beginning of iteration 4 : Cost = 0.0
----------- Jnew in iteration 4 = 0.0 (kappa = 0.14641000000000004 )
----------------------------
======= step accepted ======
----------------------------
iter4, Jnew = 0.0(kappa = 0.16105100000000006 )
end of iter 4
================ iteration 5 ===================
cost in beginning of iteration 5 : Cost = 0.0
----------- Jnew in iteration 5 = 0.0 (kappa = 0.16105100000000006 )
----------------------------
======= step accepted ======
----------------------------
iter5, Jnew = 0.0(kappa = 0.17715610000000007 )
end of iter 5
================ iteration 6 ===================
cost in beginning of iteration 6 : Cost = 0.0
----------- Jnew in iteration 6 = 0.0 (kappa = 0.17715610000000007 )
----------------------------
======= step accepted ======
----------------------------
iter6, Jnew = 0.0(kappa = 0.1948717100000001 )
end of iter 6
================ iteration 7 ===================
cost in beginning of iteration 7 : Cost = 0.0
----------- Jnew in iteration 7 = 0.0 (kappa = 0.1948717100000001 )
----------------------------
======= step accepted ======
----------------------------
iter7, Jnew = 0.0(kappa = 0.2143588810000001 )
end of iter 7
================ iteration 8 ===================
cost in beginning of iteration 8 : Cost = 0.0
----------- Jnew in iteration 8 = 0.0 (kappa = 0.2143588810000001 )
----------------------------
======= step accepted ======
----------------------------
iter8, Jnew = 0.0(kappa = 0.23579476910000013 )
end of iter 8
================ iteration 9 ===================
cost in beginning of iteration 9 : Cost = 0.0
----------- Jnew in iteration 9 = 0.0 (kappa = 0.23579476910000013 )
----------------------------
======= step accepted ======
----------------------------
iter9, Jnew = 0.0(kappa = 0.25937424601000014 )
end of iter 9
================ iteration 10 ===================
cost in beginning of iteration 10 : Cost = 0.0
----------- Jnew in iteration 10 = -0.05669363163218001 (kappa = 0.25937424601000014 )
----------------------------
======= step accepted ======
----------------------------
iter10, Jnew = -0.05669363163218001(kappa = 0.28531167061100016 )
end of iter 10
================ iteration 11 ===================
cost in beginning of iteration 11 : Cost = -0.05669363163218001
----------- Jnew in iteration 11 = -0.09395754397613308 (kappa = 0.28531167061100016 )
----------------------------
======= step accepted ======
----------------------------
iter11, Jnew = -0.09395754397613308(kappa = 0.3138428376721002 )
end of iter 11
================ iteration 12 ===================
cost in beginning of iteration 12 : Cost = -0.09395754397613308
----------- Jnew in iteration 12 = -0.10873311075998884 (kappa = 0.3138428376721002 )
----------------------------
======= step accepted ======
----------------------------
iter12, Jnew = -0.10873311075998884(kappa = 0.3452271214393102 )
end of iter 12
================ iteration 13 ===================
cost in beginning of iteration 13 : Cost = -0.10873311075998884
----------- Jnew in iteration 13 = -0.1139924216640512 (kappa = 0.3452271214393102 )
----------------------------
======= step accepted ======
----------------------------
iter13, Jnew = -0.1139924216640512(kappa = 0.37974983358324127 )
end of iter 13
================ iteration 14 ===================
cost in beginning of iteration 14 : Cost = -0.1139924216640512
----------- Jnew in iteration 14 = -0.11559682782105983 (kappa = 0.37974983358324127 )
----------------------------
======= step accepted ======
----------------------------
iter14, Jnew = -0.11559682782105983(kappa = 0.41772481694156544 )
end of iter 14
================ iteration 15 ===================
cost in beginning of iteration 15 : Cost = -0.11559682782105983
----------- Jnew in iteration 15 = -0.11600240320506246 (kappa = 0.41772481694156544 )
----------------------------
======= step accepted ======
----------------------------
iter15, Jnew = -0.11600240320506246(kappa = 0.45949729863572203 )
end of iter 15
================ iteration 16 ===================
cost in beginning of iteration 16 : Cost = -0.11600240320506246
----------- Jnew in iteration 16 = -0.11608391408168875 (kappa = 0.45949729863572203 )
----------------------------
======= step accepted ======
----------------------------
iter16, Jnew = -0.11608391408168875(kappa = 0.5054470284992942 )
end of iter 16
================ iteration 17 ===================
cost in beginning of iteration 17 : Cost = -0.11608391408168875
----------- Jnew in iteration 17 = -0.11613372171503614 (kappa = 0.5054470284992942 )
----------------------------
======= step accepted ======
----------------------------
iter17, Jnew = -0.11613372171503614(kappa = 0.5559917313492238 )
end of iter 17
================ iteration 18 ===================
cost in beginning of iteration 18 : Cost = -0.11613372171503614
----------- Jnew in iteration 18 = -0.1161363391171118 (kappa = 0.5559917313492238 )
----------------------------
======= step accepted ======
----------------------------
iter18, Jnew = -0.1161363391171118(kappa = 0.6115909044841462 )
end of iter 18
================ iteration 19 ===================
cost in beginning of iteration 19 : Cost = -0.1161363391171118
----------- Jnew in iteration 19 = -0.1161364361153799 (kappa = 0.6115909044841462 )
----------------------------
======= step accepted ======
----------------------------
iter19, Jnew = -0.1161364361153799(kappa = 0.6727499949325609 )
end of iter 19
[ ]: