{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7.1 Shape optimization via the shape derivative\n",
    "In this tutorial we discuss the implementation of shape optimization algorithms to solve\n",
    "\n",
    "$$ \\min_{\\Omega} J(\\Omega) := \\int_{\\Omega} f(x)\\;dx, \\quad f\\in C^1(\\mathbf{R}^d). $$\n",
    "The analytic solution to this problem is $$\\Omega^*:= \\{x\\in \\mathbf{R}^d:\\; f(x)\\le 0\\}.$$ \n",
    "\n",
    "This problem is solved by fixing an initial guess $\\Omega_0\\subset \\mathbf{R}^d$ and then \n",
    "considering the problem\n",
    "$$ \\min_{X} J((\\mbox{id}+X)(\\Omega_0)), $$\n",
    "where $X:\\mathbf{R}^d \\to \\mathbf{R}^d$ are (at least) Lipschitz vector fields. We approximate $X$ by a finite element function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initial geometry and shape function\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We choose $f$ as \n",
    " $$ \n",
    "       \\begin{array}{rl}\n",
    "       f(x,y) =& \\left(\\sqrt{(x - a)^2 + b y^2} - 1 \\right) \\left(\\sqrt{(x + a)^2 + b y^2} - 1 \\right) \\\\\n",
    "                & \\left(\\sqrt{b  x^2 + (y - a)^2} - 1 \\right)  \\left(\\sqrt{b  x^2 + (y + a)^2} - 1 \\right) - \\varepsilon.\n",
    "       \\end{array}\n",
    "       $$\n",
    "where $\\varepsilon = 0.001, a = 4/5, b = 2$. The corresponding zero level sets of this function are as follows. The green area indicates where $f$ is negative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "    <img src=\"opt_vs_init.png\" width=\"800\" aligh=\"center\"/>\n",
    "</center>\n",
    "\n",
    "Let us set up the geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate the geometry\n",
    "from ngsolve import *\n",
    "\n",
    "# load Netgen/NGSolve and start the gui\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.geom2d import SplineGeometry\n",
    "\n",
    "geo = SplineGeometry()\n",
    "geo.AddCircle((0,0), r = 2.5)\n",
    "\n",
    "ngmesh = geo.GenerateMesh(maxh = 0.08)\n",
    "mesh = Mesh(ngmesh)\n",
    "mesh.Curve(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#gfset.vec[:]=0\n",
    "#gfset.Set((0.2*x,0))\n",
    "#mesh.SetDeformation(gfset)\n",
    "#scene = Draw(gfzero,mesh,\"gfset\", deformation = gfset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the function f and its gradient\n",
    "a =4.0/5.0\n",
    "b = 2\n",
    "f = CoefficientFunction((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) - 0.001)\n",
    "\n",
    "# gradient of f defined as vector valued coefficient function\n",
    "grad_f = CoefficientFunction((f.Diff(x),f.Diff(y)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# vector space for shape gradient\n",
    "VEC = H1(mesh, order=1, dim=2)\n",
    "\n",
    "# grid function for deformation field\n",
    "gfset = GridFunction(VEC)\n",
    "gfX = GridFunction(VEC)\n",
    "\n",
    "Draw(gfset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Shape derivative\n",
    "---\n",
    "\n",
    "$$DJ(\\Omega)(X) = \\int_\\Omega f \\mbox{div}(X) + \\nabla f\\cdot X\\;dx$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test and trial functions\n",
    "PHI, PSI = VEC.TnT()\n",
    "\n",
    "# shape derivative\n",
    "dJOmega = LinearForm(VEC)\n",
    "dJOmega += (div(PSI)*f + InnerProduct(grad_f, PSI) )*dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bilinear form\n",
    "---\n",
    "$$ (\\varphi,\\psi) \\mapsto b_\\Omega(\\varphi,\\psi):= \\int_\\Omega (\\nabla \\varphi+\\nabla \\varphi^\\top): \\nabla\\psi+\\varphi\\cdot \\psi\\; dx.$$\n",
    "to compute the gradient $X:= \\mbox{grad}J(\\Omega)$ by\n",
    "$$ b_\\Omega(X,\\psi) = DJ(\\Omega)(\\psi)\\quad \\text{ for all } \\quad \\psi \\in H^1(\\Omega)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# bilinear form for H1 shape gradient; aX represents b_\\Omega\n",
    "aX = BilinearForm(VEC)\n",
    "aX += InnerProduct(grad(PHI) + grad(PHI).trans, grad(PSI)) * dx\n",
    "aX += InnerProduct(PHI, PSI) * dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The first optimisation step\n",
    "Fix an initial domain $\\Omega_0$ and define\n",
    "$$ \\mathcal J_{\\Omega_0}(X):= J((\\mbox{id} + X)(\\Omega_0)). $$\n",
    "\n",
    "Then we have the following relation between the derivative of $\\mathcal J_{\\Omega_0}$ and the shape derivative of $J$:\n",
    "$$\n",
    "D\\mathcal J_{\\Omega_0}(X_n)(X) =   DJ((\\mbox{id}+X_n)(\\Omega_0))(X\\circ(\\mbox{id}+X_n)^{-1}).\n",
    "$$\n",
    "Here\n",
    "\n",
    "- $(\\mbox{id}+X_n)(\\Omega_0)$ is current shape\n",
    "\n",
    "Now we note that $\\varphi \\mapsto \\varphi\\circ(\\mbox{id}+X_n)^{-1}$ maps the finite element space on $(\\mbox{id}+X_n)(\\Omega_0)$ to the finite element space on $\\Omega_0$. \n",
    "Therefore the following are equivalent:\n",
    "\n",
    "- assembling $\\varphi \\mapsto D\\mathcal J_{\\Omega_0}(X_n)(\\varphi)$ on the fixed domain $\\Omega_0$ \n",
    "- assembling $\\varphi \\mapsto DJ((\\mbox{id}+X_n)(\\Omega_0))(\\varphi)$ on the deformed domain $(\\mbox{id}+X_n)(\\Omega_0)$. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# deform current domain with gfset\n",
    "mesh.SetDeformation(gfset)\n",
    "\n",
    "# assemble on deformed domain\n",
    "aX.Assemble()\n",
    "dJOmega.Assemble()\n",
    "\n",
    "mesh.UnsetDeformation()\n",
    "# unset deformation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's make one optimization step with step size $\\alpha_1>0$:\n",
    "\n",
    "$$\\Omega_1 = (\\mbox{id} - \\alpha_1 X_0)(\\Omega_0).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute X_0\n",
    "gfX.vec.data = aX.mat.Inverse(VEC.FreeDofs(), inverse=\"sparsecholesky\") * dJOmega.vec\n",
    "\n",
    "print(\"current cost \", Integrate(f*dx, mesh))\n",
    "print(\"Gradient norm \", Norm(gfX.vec))\n",
    "\n",
    "alpha = 20.0 / Norm(gfX.vec)\n",
    "\n",
    "gfset.vec[:]=0\n",
    "scene = Draw(gfset)\n",
    "# input(\"A\")\n",
    "gfset.vec.data -= alpha * gfX.vec\n",
    "mesh.SetDeformation(gfset)\n",
    "#draw deformed shape\n",
    "scene.Redraw()\n",
    "# input(\"B\")\n",
    "\n",
    "print(\"cost after gradient step:\", Integrate(f, mesh))\n",
    "mesh.UnsetDeformation()\n",
    "scene.Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimisation loop\n",
    "Now we can set up an optimisation loop. In the second step we compute\n",
    "$$ \\Omega_2 = (\\mbox{id} - \\alpha_0 X_0 - \\alpha_1 X_1)(\\Omega_0)$$\n",
    "by the same procedure as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "\n",
    "iter_max = 50\n",
    "gfset.vec[:] = 0\n",
    "mesh.SetDeformation(gfset)\n",
    "scene = Draw(gfset,mesh,\"gfset\")\n",
    "\n",
    "converged = False\n",
    "\n",
    "alpha =0.11#100.0 / Norm(gfX.vec)\n",
    "# input(\"A\")\n",
    "for k in range(iter_max):\n",
    "    mesh.SetDeformation(gfset)\n",
    "    scene.Redraw()\n",
    "    Jold = Integrate(f, mesh)\n",
    "    print(\"cost at iteration \", k, ': ', Jold)\n",
    "    \n",
    "    # assemble bilinear form\n",
    "    aX.Assemble()\n",
    "    \n",
    "    # assemble shape derivative\n",
    "    dJOmega.Assemble()\n",
    "        \n",
    "    mesh.UnsetDeformation()\n",
    "\n",
    "    gfX.vec.data = aX.mat.Inverse(VEC.FreeDofs(), inverse=\"sparsecholesky\") * dJOmega.vec\n",
    "    \n",
    "    # step size control\n",
    "    \n",
    "    gfset_old = gfset.vec.CreateVector()\n",
    "    gfset_old.data = gfset.vec    \n",
    "    \n",
    "    Jnew = Jold + 1\n",
    "    while Jnew > Jold:\n",
    "\n",
    "        gfset.vec.data = gfset_old\n",
    "        gfset.vec.data -= alpha * gfX.vec\n",
    "\n",
    "        mesh.SetDeformation(gfset)\n",
    "        \n",
    "        Jnew = Integrate(f, mesh)\n",
    "        \n",
    "        mesh.UnsetDeformation()\n",
    "\n",
    "        if Jnew > Jold:\n",
    "            print(\"reducing step size\")\n",
    "            alpha = 0.9*alpha\n",
    "        else:\n",
    "            print(\"linesearch step accepted\")\n",
    "            alpha = alpha*1.5\n",
    "            break\n",
    "            \n",
    "    print(\"step size: \", alpha, '\\n')\n",
    "\n",
    "    time.sleep(0.1)\n",
    "    Jold = Jnew\n",
    "        \n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n",
    "## Using SciPy optimize toolbox\n",
    "\n",
    "We use the package scipy.optimize and wrap the shape functions around it. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first setup the initial geometry as a circle of radius r = 2.5 (as before) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The code in this cell is the same as in the example above.\n",
    "\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.geom2d import SplineGeometry\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "ngsglobals.msg_level = 1\n",
    "\n",
    "geo = SplineGeometry()\n",
    "geo.AddCircle((0,0), r = 2.5)\n",
    "\n",
    "ngmesh = geo.GenerateMesh(maxh = 0.08)\n",
    "mesh = Mesh(ngmesh)\n",
    "mesh.Curve(2)\n",
    "\n",
    "# define the function f\n",
    "\n",
    "a =4.0/5.0\n",
    "b = 2\n",
    "f = CoefficientFunction((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) - 0.001)\n",
    "\n",
    "# Now we define the finite element space VEC in which we compute the shape gradient\n",
    "\n",
    "# element order\n",
    "order = 1\n",
    "\n",
    "VEC = H1(mesh, order=order, dim=2)\n",
    "\n",
    "# define test and trial functions\n",
    "PHI = VEC.TrialFunction()\n",
    "PSI = VEC.TestFunction()\n",
    "\n",
    "\n",
    "# define grid function for deformation of mesh\n",
    "gfset = GridFunction(VEC)\n",
    "gfset.Set((0,0))\n",
    "\n",
    "# only for new gui\n",
    "#scene = Draw(gfset, mesh, \"gfset\")\n",
    "\n",
    "#if scene:\n",
    "#    scene.setDeformation(True)\n",
    "\n",
    "# plot the mesh and visualise deformation\n",
    "#Draw(gfset,mesh,\"gfset\")\n",
    "#SetVisualization (deformation=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Shape derivative"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$DJ(\\Omega)(X) = \\int_\\Omega f \\mbox{div}(X) + \\nabla f\\cdot X\\;dx$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fX = LinearForm(VEC)\n",
    "\n",
    "# analytic shape derivative\n",
    "fX += f*div(PSI)*dx + (f.Diff(x,PSI[0])) *dx + (f.Diff(y,PSI[1])) *dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bilinear form for shape gradient\n",
    "Define the bilinear form\n",
    "$$ (\\varphi,\\psi) \\mapsto b_\\Omega(\\varphi,\\psi):= \\int_\\Omega (\\nabla \\varphi+\\nabla \\varphi^\\top): \\nabla\\psi + \\varphi\\cdot\\psi\\; dx $$\n",
    "to compute the gradient $X =: \\mbox{grad}J$ by\n",
    "$$ b_\\Omega(X,\\psi) = DJ(\\Omega)(\\psi)\\quad \\text{ for all } \\psi\\in H^1(\\Omega)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cauchy-Riemann descent CR\n",
    "CR = False\n",
    "\n",
    "# bilinear form for gradient\n",
    "aX = BilinearForm(VEC)\n",
    "aX += InnerProduct(grad(PHI)+grad(PHI).trans, grad(PSI))*dx + InnerProduct(PHI,PSI)*dx\n",
    "\n",
    "## Cauchy-Riemann regularisation\n",
    "\n",
    "if CR == True:\n",
    "    gamma = 100\n",
    "    aX += gamma * (PHI.Deriv()[0,0] - PHI.Deriv()[1,1])*(PSI.Deriv()[0,0] - PSI.Deriv()[1,1]) *dx\n",
    "    aX += gamma * (PHI.Deriv()[1,0] + PHI.Deriv()[0,1])*(PSI.Deriv()[1,0] + PSI.Deriv()[0,1]) *dx\n",
    "\n",
    "\n",
    "aX.Assemble()\n",
    "invaX = aX.mat.Inverse(VEC.FreeDofs(), inverse=\"sparsecholesky\")\n",
    "\n",
    "gfX = GridFunction(VEC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wrapping the shape function and its gradient\n",
    "\n",
    "Now we define the shape function $J$ and its gradient grad(J)  use the shape derivative $DJ$. The arguments of $J$ and grad(J) are vector in $\\mathbf{R}^d$, where $d$ is the dimension of the finite element space. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def J(x_):\n",
    "        \n",
    "    gfset.Set((0,0))\n",
    "    # x_ NumPy vector \n",
    "    gfset.vec.FV().NumPy()[:] += x_\n",
    "\n",
    "    mesh.SetDeformation(gfset)\n",
    "\n",
    "    cost_c = Integrate (f, mesh)\n",
    "\n",
    "    mesh.UnsetDeformation()\n",
    "\n",
    "    Redraw(blocking=True)\n",
    "\n",
    "    return cost_c\n",
    "\n",
    "def gradJ(x_, euclid = False):\n",
    "\n",
    "    gfset.Set((0,0))\n",
    "    # x_ NumPy vector \n",
    "    gfset.vec.FV().NumPy()[:] += x_\n",
    "\n",
    "    mesh.SetDeformation(gfset)\n",
    "\n",
    "    fX.Assemble()\n",
    "\n",
    "    mesh.UnsetDeformation()\n",
    "\n",
    "    if euclid == True:\n",
    "        gfX.vec.data = fX.vec\n",
    "    else: \n",
    "        gfX.vec.data = invaX * fX.vec\n",
    "\n",
    "    return gfX.vec.FV().NumPy().copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient descent and Armijo rule\n",
    "\n",
    "We use the functions $J$ and grad J to define a steepest descent method.We use scipy.optimize.line_search_armijo to compute the step size in each iteration. The Arjmijo rule reads\n",
    "$$ J((\\mbox{Id}+\\alpha_k X_k)(\\Omega_k)) \\le J(\\Omega_k) - c_1 \\alpha_k \\|\\nabla J(\\Omega_k)\\|^2, $$\n",
    "where $c_1:= 1e-4$ and $\\alpha_k$ is the step size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import scipy linesearch method\n",
    "from scipy.optimize.linesearch import line_search_armijo\n",
    "\n",
    "def gradient_descent(x0, J_, gradJ_):\n",
    "\n",
    "    xk_ = np.copy(x0)\n",
    "    \n",
    "    # maximal iteration\n",
    "    it_max = 50\n",
    "    # count number of function evals\n",
    "    nfval_total = 0\n",
    "\n",
    "    print('\\n')\n",
    "    for i in range(1, it_max):\n",
    "        \n",
    "        # Compute a step size using a line_search to satisfy the Wolf\n",
    "        # compute shape gradient grad J \n",
    "        grad_xk = gradJ_(xk_,euclid = False)\n",
    "        \n",
    "        # compute descent direction\n",
    "        pk = -gradJ_(xk_)\n",
    "        \n",
    "        # eval cost function\n",
    "        fold = J_(xk_)\n",
    "        \n",
    "        # perform armijo stepsize\n",
    "        \n",
    "        if CR == True:\n",
    "            alpha0 = 0.15\n",
    "        else:\n",
    "            alpha0 = 0.11\n",
    "            \n",
    "            \n",
    "        \n",
    "        step, nfval, b = line_search_armijo(J_, xk_, pk = pk, gfk = grad_xk, old_fval = fold, c1=1e-4, alpha0 = alpha0)\n",
    "        nfval_total += nfval\n",
    "        \n",
    "        # update the shape and print cost and gradient norm\n",
    "        xk_ = xk_ - step * grad_xk\n",
    "        print('Iteration ', i, '| Cost ', fold, '| grad norm', np.linalg.norm(grad_xk))\n",
    "        \n",
    "        mesh.SetDeformation(gfset)\n",
    "        scene.Redraw()\n",
    "        mesh.UnsetDeformation()\n",
    "        \n",
    "        if np.linalg.norm(gradJ_(xk_)) < 1e-4:\n",
    "\n",
    "            #print('#######################################')\n",
    "            print('\\n'+'{:<20}  {:<12} '.format(\"##################################\", \"\"))\n",
    "            print('{:<20}  {:<12} '.format(\"### success - accuracy reached ###\", \"\"))\n",
    "            print('{:<20}  {:<12} '.format(\"##################################\", \"\"))\n",
    "            print('{:<20}  {:<12} '.format(\"gradient norm: \", np.linalg.norm(gradJ_(xk_))))\n",
    "            print('{:<20}  {:<12} '.format(\"n evals f: \", nfval_total))\n",
    "            print('{:<20}  {:<12} '.format(\"f val: \", fold) + '\\n')\n",
    "            break\n",
    "        elif i == it_max-1:\n",
    "\n",
    "            #print('######################################')\n",
    "            print('\\n'+'{:<20}  {:<12} '.format(\"#######################\", \"\"))\n",
    "            print('{:<20}  {:<12} '.format(\"### maxiter reached ###\", \"\"))\n",
    "            print('{:<20}  {:<12} '.format(\"#######################\", \"\"))\n",
    "            print('{:<20}  {:<12} '.format(\"gradient norm: \", np.linalg.norm(gradJ_(xk_))))\n",
    "            print('{:<20}  {:<12} '.format(\"n evals f: \", nfval_total))\n",
    "            print('{:<20}  {:<12} '.format(\"f val: \", fold) + '\\n')\n",
    "        \n",
    " \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now we are ready to run our first optimisation algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfset.vec[:]=0\n",
    "x0 = gfset.vec.FV().NumPy() # initial guess = 0 \n",
    "scene = Draw(gfset, mesh, \"gfset\")\n",
    "gradient_descent(x0, J, gradJ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### L-BFGS method\n",
    "Now we use the L-BFGS method provided by SciPy. The BFGS method requires the shape function $J$ and its gradient grad J. We can also specify additional arguments in options, such as maximal iterations and the gradient tolerance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the BFGS method we replace\n",
    "$$ b_\\Omega(X,\\varphi) = \\nabla J(\\Omega)(\\varphi) $$\n",
    "by \n",
    "$$ H_\\Omega(X,\\varphi) = \\nabla J(\\Omega)(\\varphi),$$\n",
    "where $H_\\Omega$ is an approximation of the second shape derivative at $\\Omega$. On the discrete level we solve\n",
    "$$\n",
    "\\begin{array}{rl}\n",
    " \\text{ solve }\\quad  B_nX_k & = \\nabla J(\\Omega_k)  \\\\\n",
    "    & \\\\ \n",
    " \\text{ update } \\quad s_k & := -\\alpha_k p_k \\\\\n",
    "   y_k & := \\nabla J_h(\\Omega_{k+1}) - \\nabla J_h(\\Omega_k) \\\\\n",
    "B_{k+1} &:= B_k + \\frac{y_k\\otimes y_k}{y_k\\cdot s_k} + \\frac{B_ks_k\\otimes B_ks_k}{B_ky_k\\cdot B_ks_k}.   \n",
    "\\end{array} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "x0 = gfset.vec.FV().NumPy()\n",
    "\n",
    "# options for optimiser\n",
    "\n",
    "options = {\"maxiter\":1000,\n",
    "              \"disp\":True, \n",
    "              \"gtol\":1e-10}\n",
    "\n",
    "# we use quasi-Newton method L-BFGS\n",
    "minimize(J, x0, method='L-BFGS-B', jac=gradJ, options=options)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Improving mesh quality via Cauchy-Riemann equations\n",
    "\n",
    "In the previous section we computed the shape gradient grad J:= X of $J$ at $\\Omega$ via \n",
    "$$ \\int_{\\Omega} \\nabla X : \\nabla \\varphi + X\\cdot \\varphi\\; dx = DJ(\\Omega)(\\varphi) \\quad \\forall \\varphi \\in H^1(\\Omega)^d. $$\n",
    "This may lead to overly stretched or even degenerated triangles. One way to improve this is to modify the above equation by \n",
    "$$ \\int_{\\Omega} \\nabla X : \\nabla \\varphi + X\\cdot \\varphi +  \\color{red}\\gamma   \\color{red}B\\color{red}X\\cdot \\color{red}B\\color{red}\\varphi\\;dx = DJ(\\Omega)(\\varphi) \\; \\forall \\varphi \\in H^1(\\Omega)^d, $$\n",
    "where\n",
    "$$\n",
    "B := \\begin{pmatrix}\n",
    "-\\partial_x & \\partial_y \\\\\n",
    "\\partial_y & \\partial_x\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "The two equations $BX = 0$ are precisely the Cauchy-Riemann equations $-\\partial_x X_1 + \\partial_y X_2=0$ and $\\partial_y X_1 + \\partial_x X_2$. So the bigger $\\gamma$ the more angle preserving is the gradient. So by adding the B term we enforce conformaty with strength $\\gamma$. \n",
    "\n",
    "\n",
    "This only means we need to add the $\\alpha$ term to the above bilinear form aX:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 100\n",
    "\n",
    "aX += alpha * (PHI.Deriv()[0,0] - PHI.Deriv()[1,1])*(PSI.Deriv()[0,0] - PSI.Deriv()[1,1]) *dx\n",
    "aX += alpha * (PHI.Deriv()[1,0] + PHI.Deriv()[0,1])*(PSI.Deriv()[1,0] + PSI.Deriv()[0,1]) *dx"
   ]
  },
  {
   "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
}
