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
\[\mathcal J(\Omega) = \int_\Omega f(x) \, \mbox dx,\qquad f:\mathbb{R}^d\to \mathbb{R}\]
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
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 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}\]
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))
We represent the level set function \(\psi\) as a GridFunction on the mesh of \(\mathsf{D}\) and visualize the zero level set of the function \(f\).
[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(IfPos(f,0,1), 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")
kappaMax = 1
kappa = 0.1*kappaMax
We define \(\mathsf{NegPos}\) and \(\mathsf{PosNeg}\) as the topological derivative in \(\Omega\) and \(\mathsf{D}\setminus\overline{\Omega}\), respectively:
\[\begin{split}\begin{array}{rl}
\mathsf{NegPos}(x) := -DJ(\Omega,\omega)(x) &= f(x) \quad \text{ for }x\in \Omega \\
\mathsf{PosNeg}(x) := DJ(\Omega,\omega)(x) &= f(x) \quad \text{ for }x\in \overline{\mathsf{D}}\setminus\Omega
\end{array}\end{split}\]
and the level set update function (also called generalized topological derivative)
\[g_\Omega(x) := \chi_\Omega (x)\mathsf{NegPos}(x) + (1-\chi_\Omega(x))\mathsf{PosNeg}(x) = f(x).\]
[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()
The idea now is to make a fixed point iteration for \(\ell\ge 0\)
\[\psi_{\ell+1} = (1-\kappa) \psi_\ell + \kappa \frac{g_{\psi_\ell}}{|g_{\psi_\ell}|_{L_2(\mathsf{D})}}, \qquad \kappa\in [0,1]\]
Let us do one iteration:
[7]:
scene3 = Draw(CutRatio, mesh)
[8]:
kappaTest = 0.8
for k in range(3):
psi.vec.data = (1-kappaTest)*psi.vec + kappaTest * TD_node.vec
InterpolateLevelSetToElems(psi, 1, 0, CutRatio, mesh, EPS)
scene3.Redraw()
Now, let us repeat the procedure in the following iterative algorithm:
[9]:
psi.Set(1)
InterpolateLevelSetToElems(psi, 1, 0, CutRatio, mesh, EPS)
scene_cr = Draw(CutRatio, mesh)
[10]:
iter_max = 20
scene_cr.Redraw()
psi_new.vec.data = psi.vec
kappa = 0.1*kappaMax
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_cr.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
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.05669364161907429 (kappa = 0.25937424601000014 )
----------------------------
======= step accepted ======
----------------------------
iter10, Jnew = -0.05669364161907429(kappa = 0.28531167061100016 )
end of iter 10
================ iteration 11 ===================
cost in beginning of iteration 11 : Cost = -0.05669364161907429
----------- Jnew in iteration 11 = -0.09395828820090642 (kappa = 0.28531167061100016 )
----------------------------
======= step accepted ======
----------------------------
iter11, Jnew = -0.09395828820090642(kappa = 0.3138428376721002 )
end of iter 11
================ iteration 12 ===================
cost in beginning of iteration 12 : Cost = -0.09395828820090642
----------- Jnew in iteration 12 = -0.10873188545899315 (kappa = 0.3138428376721002 )
----------------------------
======= step accepted ======
----------------------------
iter12, Jnew = -0.10873188545899315(kappa = 0.3452271214393102 )
end of iter 12
================ iteration 13 ===================
cost in beginning of iteration 13 : Cost = -0.10873188545899315
----------- Jnew in iteration 13 = -0.11398841292489859 (kappa = 0.3452271214393102 )
----------------------------
======= step accepted ======
----------------------------
iter13, Jnew = -0.11398841292489859(kappa = 0.37974983358324127 )
end of iter 13
================ iteration 14 ===================
cost in beginning of iteration 14 : Cost = -0.11398841292489859
----------- Jnew in iteration 14 = -0.11559719262617797 (kappa = 0.37974983358324127 )
----------------------------
======= step accepted ======
----------------------------
iter14, Jnew = -0.11559719262617797(kappa = 0.41772481694156544 )
end of iter 14
================ iteration 15 ===================
cost in beginning of iteration 15 : Cost = -0.11559719262617797
----------- Jnew in iteration 15 = -0.11600233851477758 (kappa = 0.41772481694156544 )
----------------------------
======= step accepted ======
----------------------------
iter15, Jnew = -0.11600233851477758(kappa = 0.45949729863572203 )
end of iter 15
================ iteration 16 ===================
cost in beginning of iteration 16 : Cost = -0.11600233851477758
----------- Jnew in iteration 16 = -0.11608384153972257 (kappa = 0.45949729863572203 )
----------------------------
======= step accepted ======
----------------------------
iter16, Jnew = -0.11608384153972257(kappa = 0.5054470284992942 )
end of iter 16
================ iteration 17 ===================
cost in beginning of iteration 17 : Cost = -0.11608384153972257
----------- Jnew in iteration 17 = -0.11613377612863536 (kappa = 0.5054470284992942 )
----------------------------
======= step accepted ======
----------------------------
iter17, Jnew = -0.11613377612863536(kappa = 0.5559917313492238 )
end of iter 17
================ iteration 18 ===================
cost in beginning of iteration 18 : Cost = -0.11613377612863536
----------- Jnew in iteration 18 = -0.11613635207766497 (kappa = 0.5559917313492238 )
----------------------------
======= step accepted ======
----------------------------
iter18, Jnew = -0.11613635207766497(kappa = 0.6115909044841462 )
end of iter 18
================ iteration 19 ===================
cost in beginning of iteration 19 : Cost = -0.11613635207766497
----------- Jnew in iteration 19 = -0.11613642191697138 (kappa = 0.6115909044841462 )
----------------------------
======= step accepted ======
----------------------------
iter19, Jnew = -0.11613642191697138(kappa = 0.6727499949325609 )
end of iter 19
[ ]: