{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7.6 Topology Optimization based on the Topological Derivative"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are again interested in minimizing a shape function of the type\n",
    "$$\n",
    "    \\mathcal J(\\Omega) = \\int_\\Omega f(x) \\, \\mbox dx,\\qquad f:\\mathbb{R}^d\\to \\mathbb{R} \n",
    "$$\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# module to generate the 2D geometry\n",
    "from netgen.geom2d import SplineGeometry\n",
    "\n",
    "\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "import numpy as np\n",
    "\n",
    "from interpolations import InterpolateLevelSetToElems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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](interpolations.py) . \n",
    "\n",
    "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:\n",
    "\n",
    "$$\n",
    "\\begin{array}{rl}\n",
    "    \\psi(x) < 0 & \\quad \\Longleftrightarrow  \\quad x \\in \\Omega \\\\\n",
    "    \\psi(x) = 0 & \\quad \\Longleftrightarrow  \\quad x \\in \\partial \\Omega \\\\\n",
    "    \\psi(x) > 0 & \\quad \\Longleftrightarrow  \\quad x \\in \\mathsf{D} \\setminus \\overline \\Omega.\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define the domain $\\mathsf{D} = [-2, 2]^2$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myMAXH = 0.05\n",
    "geo = SplineGeometry()\n",
    "R = 2\n",
    "\n",
    "## add a rectangle\n",
    "geo.AddRectangle(p1=(-R,-R),\n",
    "                 p2=(R,R),\n",
    "                 bc=\"rectangle\",\n",
    "                 leftdomain=1,\n",
    "                 rightdomain=0)\n",
    "\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=myMAXH))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# H1-conforming finite element space\n",
    "fes = H1(mesh, order=1) # C0 conforming finite element space\n",
    "pwc = L2(mesh)   # piecewise constant (on each element) space\n",
    "\n",
    "psi = GridFunction(fes)       # level set function\n",
    "psi.vec[:]=1\n",
    "psi_new = GridFunction(fes)   # grid function for line search\n",
    "\n",
    "a = 4.0/5.0\n",
    "b = 2\n",
    "e = 0.001\n",
    "f = CoefficientFunction( 0.1*( (sqrt((x - a)**2 + b * y**2) - 1) \\\n",
    "                * (sqrt((x + a)**2 + b * y**2) - 1) \\\n",
    "                * (sqrt(b * x**2 + (y - a)**2) - 1) \\\n",
    "                * (sqrt(b * x**2 + (y + a)**2) - 1) - e) )\n",
    "\n",
    "Draw(IfPos(f,0,1), mesh)\n",
    "EPS = myMAXH * 1e-6      #variable defining when a point is on the interface and when not"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Cost(psi):\n",
    "    return IfPos(psi,0,f)*dx\n",
    "\n",
    "print(\"initial cost = \", Integrate(Cost(psi), mesh))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The evolution of the level set function is guided by the $\\textbf{topological derivative}$.\n",
    "We define and draw the topological derivative in $\\Omega$, $\\mathsf{D}\\setminus \\Omega$ and the generalized topological derivative $g$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TDPosNeg_pwc = GridFunction(pwc)\n",
    "TDNegPos_pwc = GridFunction(pwc)\n",
    "TD_pwc = GridFunction(pwc)\n",
    "CutRatio = GridFunction(pwc)\n",
    "\n",
    "TD_node = GridFunction(fes)\n",
    "\n",
    "scene1 = Draw(TD_node, mesh, \"TD_node\")\n",
    "scene2 = Draw(TD_pwc, mesh, \"TD_pwc\")\n",
    "\n",
    "kappaMax = 1\n",
    "kappa = 0.1*kappaMax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define $\\mathsf{NegPos}$ and $\\mathsf{PosNeg}$ as the topological derivative in $\\Omega$ and $\\mathsf{D}\\setminus\\overline{\\Omega}$, respectively:\n",
    "$$\n",
    "\\begin{array}{rl}\n",
    "    \\mathsf{NegPos}(x) :=  -DJ(\\Omega,\\omega)(x) &= f(x) \\quad \\text{ for }x\\in \\Omega \\\\\n",
    "    \\mathsf{PosNeg}(x) :=   DJ(\\Omega,\\omega)(x) &= f(x) \\quad \\text{ for }x\\in \\overline{\\mathsf{D}}\\setminus\\Omega\n",
    "\\end{array}\n",
    "$$\n",
    "and the level set update function (also called generalized topological derivative)\n",
    "$$\n",
    "g_\\Omega(x) := \\chi_\\Omega (x)\\mathsf{NegPos}(x) + (1-\\chi_\\Omega(x))\\mathsf{PosNeg}(x) = f(x). \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TDNegPos_pwc.Set(f)\n",
    "TDPosNeg_pwc.Set(f)\n",
    "\n",
    "InterpolateLevelSetToElems(psi, 1, 0, CutRatio, mesh, EPS)\n",
    "\n",
    "# compute the combined topological derivative using the cut ratio information\n",
    "for j in range(len(TD_pwc.vec)):\n",
    "    TD_pwc.vec[j] = CutRatio.vec[j] * TDNegPos_pwc.vec[j] + (1-CutRatio.vec[j])*TDPosNeg_pwc.vec[j]\n",
    "        \n",
    "TD_node.Set(TD_pwc)\n",
    "\n",
    "NormTD = sqrt(Integrate(TD_node**2*dx, mesh)) # L2 norm of TD_node\n",
    "\n",
    "TD_node.vec.data = 1/NormTD * TD_node.vec        \n",
    "\n",
    "scene1.Redraw()\n",
    "scene2.Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea now is to make a fixed point iteration for $\\ell\\ge 0$\n",
    "$$\n",
    "    \\psi_{\\ell+1} = (1-\\kappa) \\psi_\\ell + \\kappa \\frac{g_{\\psi_\\ell}}{|g_{\\psi_\\ell}|_{L_2(\\mathsf{D})}}, \\qquad \\kappa\\in [0,1]\n",
    "$$\n",
    "Let us do one iteration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scene3 = Draw(CutRatio, mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kappaTest = 0.8\n",
    "for k in range(3):\n",
    "    psi.vec.data =  (1-kappaTest)*psi.vec + kappaTest * TD_node.vec\n",
    "InterpolateLevelSetToElems(psi, 1, 0, CutRatio, mesh, EPS)\n",
    "scene3.Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let us repeat the procedure in the following iterative algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psi.Set(1)\n",
    "InterpolateLevelSetToElems(psi, 1, 0, CutRatio, mesh, EPS)\n",
    "scene_cr = Draw(CutRatio, mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iter_max = 20\n",
    "\n",
    "scene_cr.Redraw()\n",
    "psi_new.vec.data = psi.vec\n",
    "kappa = 0.1*kappaMax\n",
    "for k in range(iter_max):\n",
    "\n",
    "        print(\"================ iteration \", k, \"===================\")\n",
    "\n",
    "        # copy new levelset data from psi_new into psi\n",
    "        psi.vec.data = psi_new.vec\n",
    "        scene_cr.Redraw()\n",
    "        \n",
    "        # current cost function value\n",
    "        J_current = Integrate(Cost(psi),mesh)\n",
    "\n",
    "        print(\"cost in beginning of iteration\", k, \": Cost = \", J_current)\n",
    "        \n",
    "        # level set update function\n",
    "        TDPosNeg_pwc.Set(f)\n",
    "        TDNegPos_pwc.Set(f)\n",
    "\n",
    "        # get fraction of intersecting elements (for psi<0)\n",
    "        InterpolateLevelSetToElems(psi, 1, 0, CutRatio, mesh, EPS)\n",
    "        \n",
    "        # weight the topological derivative with CutRatio\n",
    "        for j in range(len(TD_pwc.vec)):\n",
    "            TD_pwc.vec[j] = CutRatio.vec[j] * TDNegPos_pwc.vec[j] + (1-CutRatio.vec[j])*TDPosNeg_pwc.vec[j]\n",
    "        \n",
    "        TD_node.Set(TD_pwc)\n",
    "\n",
    "        NormTD = sqrt(Integrate(TD_node**2*dx, mesh)) # L2 norm of TD_node\n",
    "        TD_node.vec.data = 1/NormTD * TD_node.vec\n",
    "\n",
    "        psi.vec.data  =  (1-kappa)*psi.vec + kappa * TD_node.vec\n",
    "\n",
    "        for j in range(10): \n",
    "\n",
    "            # update the level set function\n",
    "            psi_new.vec.data  =  (1-kappa)*psi.vec + kappa * TD_node.vec\n",
    "\n",
    "            Jnew = Integrate(Cost(psi_new), mesh)\n",
    "            \n",
    "            Redraw()\n",
    "\n",
    "            print(\"----------- Jnew in  iteration \", k, \" = \", Jnew, \"(kappa = \", kappa, \")\")\n",
    "            print('')\n",
    "\n",
    "            if Jnew <= J_current+1e-10:\n",
    "                print(\"----------------------------\")\n",
    "                print(\"======= step accepted ======\")\n",
    "                print(\"----------------------------\")\n",
    "                kappa = min(1.1*kappa, kappaMax)\n",
    "                break\n",
    "            else:\n",
    "                kappa = kappa*0.8\n",
    "        \n",
    "        print(\"iter\" + str(k) + \", Jnew = \" + str(Jnew) + \"(kappa = \", kappa, \")\")\n",
    "\n",
    "        print(\"end of iter \" + str(k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}