# 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

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

[ ]: