{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 7.5 PDE-Constrained Shape Optimization (fully automated)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to solve the PDE-constrained shape optimization problem\n", "$$\n", " \\underset{\\Omega\\subset \\mathsf{D}}{\\mbox{min}} \\; J(u) := \\int_\\Omega |u-u_d|^q \\; dx, \\quad q\\ge 2\n", "$$\n", "subject to that $(\\Omega,u)$ satisfy\n", "$$\n", " \\int_\\Omega \\nabla u \\cdot \\nabla v \\; dx = \\int_\\Omega f v \\; dx \\; \\quad \\text{ for all } v \\in H_0^1(\\Omega),\n", "$$\n", "where $\\Omega \\subset \\mathbb R^2$ for given $u_d, f \\in C^1(\\mathbb R^2)$.\n", "\n", "Again, we want to compute the shape derivative by differentiation of a suitably defined perturbed Lagrangian using automated differentiation. Here this is accounted for automatically by NGSolve. For details we refer to \n", "\n", "P. Gangl, K. Sturm, M. Neunteufel, J. Schöberl.\n", "Fully and Semi-Automated Shape Differentiation in NGSolve,\n", "Struct. Multidisc. Optim., 63, pp.1579-1607, 2021." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from ngsolve import *\n", "from ngsolve.webgui import Draw\n", "from netgen.geom2d import SplineGeometry\n", "from ngsolve.solvers import *" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "geo = SplineGeometry()\n", "geo.AddCircle(c=(0.5,0.5), r=0.5, bc = 'circle')\n", "mesh = Mesh(geo.GenerateMesh(maxh = 0.08))\n", "Draw(mesh)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#given data of our problem (chosen such that \\Omega^* = [0,1]^2 is the optimal shape)\n", "f = CoefficientFunction(2*y*(1-y)+2*x*(1-x))\n", "ud = x*(1-x)*y*(1-y)\n", "\n", "grad_f = CoefficientFunction( (f.Diff(x), f.Diff(y) ) )\n", "grad_ud = CoefficientFunction( (ud.Diff(x), ud.Diff(y) ) )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fes = H1(mesh, order=2, dirichlet=\".*\")\n", "gfu = GridFunction(fes)\n", "scene_u = Draw (gfu, mesh, \"state\")\n", "\n", "gfp = GridFunction(fes)\n", "scene_p = Draw (gfp, mesh, \"adjoint\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that (for linear problems) the operator on the left hand side of the adjoint equation is just the transpose of the state operator." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Automatic Shape Differentiation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The formula for the shape derivative was derived as the partial derivative of the perturbed Lagrangian (brought back to the original domain):\n", "$$\n", " d\\mathcal J(\\Omega; X) = \\frac{\\partial}{\\partial t} \\left. \\left( \\int_\\Omega |u - u_d^t|^q \\,\\mbox{det}(F_t) \\mbox{d} x \n", " + \\int_{\\Omega} (F_t^{-\\top}\\nabla u) \\cdot (F_t^{-\\top} \\nabla p)\\, \\mbox{det}(F_t) \\, dx - \\int_{\\Omega}f^t p \\,\\mbox{det}(F_t) \\,dx \\right) \\right \\rvert_{t=0} \n", "$$\n", "where \n", "