{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Topology optimization using the SIMP method\n", "\n", "This numerical tour investigates the optimal design of an elastic structure through a topology optimization approach. It consists in finding the optimal distribution of a material in a computational domain which minimizes the compliance (or, equivalently, maximizes the stiffness) of the resulting structure under a fixed volume fraction constraint.\n", "\n", "For instance, in the bottom animation, the optimal material density (a continuous field $\\theta(x)\\in[0;1]$) of a cantilever beam evolves with the topoology optimization algorithm iterations. The first 20 iterations correspond to a \"non-penalized\" distribution where intermediate densities (i.e. $0<\\theta(x)<1$) are allowed. However, such a *gray-level* distribution is hard to realize in practice and one would like to obtain a *black-and-white* design of material/void distribution. This is what is achieved in the last iterations where some penalty is applied to intermediate densities, exhibiting finally a truss-like design. Note that the final design is not a global optimum but only a local one. Hence, different designs can well be obtained if changing the initial conditions or the mesh density for instance." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](cantilever_beam.gif \"Title\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The present tour relies on one of the most widely implemented approach in topology optimization, namely the Solid Isotropic Material with Penalisation (SIMP) method developed by Bendsoe and Sigmund [[BEN03]](#References). See also the monograph [[ALL07]](#References) and [Gregoire Allaire's class on *Optimal Design of Structures*](http://www.cmap.polytechnique.fr/~MAP562/index.php) at Ecole Polytechnique. In particular, this tour refers to [lesson 8](http://www.cmap.polytechnique.fr/~MAP562/archives/slides/lesson8.pdf) and is the direct adaptation of the [SIMP Topology Optimization example](http://www.cmap.polytechnique.fr/%7Eallaire/map562/console.simp.edp) of the [CMAP group toolbox written in Freefem++](http://www.cmap.polytechnique.fr/~allaire/freefem.html).\n", "\n", "This tour should also work in parallel.\n", "\n", "## Theory\n", "\n", "### Compliance optimization on binary shapes\n", "\n", "Let us consider a computational domain $D$ in which we seek an elastic material domain $\\Omega$ ($D\\setminus\\Omega$ representing void) so that the compliance under prescribed loading is minimal when subjected to a fixed volume constraint $\\text{Vol}(\\Omega) = \\eta \\text{Vol}(D)$. Replacing the shape $\\Omega$ by its characteristic function $\\chi(x)\\in \\{0,1\\}$, the previous problem can be formulated as the following optimization problem:\n", "\n", "$$\\begin{array}{rl}\\displaystyle{\\min_{\\chi}} & \\displaystyle{\\int_D \\boldsymbol{f}\\cdot\\boldsymbol{u}(\\chi) \\text{ dx}+\\int_{\\Gamma_N} \\boldsymbol{T}\\cdot\\boldsymbol{u}(\\chi)\\text{ dS}}\\\\\n", "\\text{s.t.} & \\int_D \\chi(x) \\text{ dx} = \\eta \\text{Vol}(D) \\\\\n", "& \\chi(x)\\in \\{0,1\\}\n", "\\end{array}$$\n", "\n", "where $\\boldsymbol{u}(\\chi)$ is the solution of the following $\\chi(x)-$dependent elasticity problem:\n", "\n", "$$\\begin{align}\n", "\\text{div}\\, \\boldsymbol{\\sigma} + \\boldsymbol{f}=0 &\\quad\\text{in } D\\\\\n", "\\boldsymbol{\\sigma} = \\chi(x)\\mathbb{C}:\\nabla^s \\boldsymbol{u} &\\quad\\text{in } D\\\\\n", "\\boldsymbol{\\sigma}\\cdot\\boldsymbol{n} = \\boldsymbol{T} &\\quad\\text{on }\\Gamma_N \\\\\n", "\\boldsymbol{u} = 0 &\\quad\\text{on }\\Gamma_D \n", "\\end{align}$$\n", "\n", "where $\\boldsymbol{f}$ is a body force, $\\boldsymbol{T}$ corresponds to surface tractions on $\\Gamma_N$ and $\\mathbb{C}$ is the elastic material stiffness tensor. Owing to the elastic nature of the problem, the compliance can also be expressed using the elastic stress energy density:\n", "\n", "$$\\int_D \\boldsymbol{f}\\cdot\\boldsymbol{u}(\\chi) \\text{ dx}+\\int_{\\Gamma_N} \\boldsymbol{T}\\cdot\\boldsymbol{u}(\\chi)\\text{ dS}=\\int_D \\boldsymbol{\\sigma}(\\chi):(\\chi(x)\\mathbb{C})^{-1}:\\boldsymbol{\\sigma}(\\chi)\\text{ dx}=\\int_D \\nabla^s \\boldsymbol{u}(\\chi):(\\chi(x)\\mathbb{C}):\\nabla^s \\boldsymbol{u}(\\chi)\\text{ dx}\n", "$$\n", "\n", "so that the above optimization problem can be reformulated as:\n", "\n", "$$\\begin{array}{rl}\\displaystyle{\\min_{\\chi,\\boldsymbol{u}}} & \\displaystyle{\\int_D \\boldsymbol{\\sigma}(\\chi):(\\chi(x)\\mathbb{C})^{-1}:\\boldsymbol{\\sigma}(\\chi)\\text{ dx}}\\\\\n", "\\text{s.t.} & \\text{div}\\, \\boldsymbol{\\sigma} + \\boldsymbol{f}=0 \\quad\\text{in } D\\\\\n", "& \\boldsymbol{\\sigma}\\cdot\\boldsymbol{n} = \\boldsymbol{T} \\quad\\text{on }\\Gamma_N \\\\\n", "& \\int_D \\chi(x) \\text{ dx} = \\eta \\text{Vol}(D) \\\\\n", "& \\chi(x)\\in \\{0,1\\}\n", "\\end{array}$$\n", "\n", "### Continuous relaxation and SIMP penalization\n", "\n", "The binary constraint $\\chi(x)\\in \\{0,1\\}$ makes the above problem extremely difficult to solve so that a classical remedy consists in relaxing the constraint by a continuous constraint $\\theta(x)\\in [0,1]$ where $\\theta(x)$ will approximate the characteristic function $\\chi(x)$ by a gray-level continuous function taking values between 0 and 1.\n", "\n", "However, in order to still obtain a final binary density distribution, the binary modulus $\\chi(x)\\mathbb{C}$ will be replaced by $\\theta(x)^p \\mathbb{C}$ with $p>1$ in order to penalize intermediate densities, yielding the so-called *SIMP formulation*:\n", "\n", "$$\\begin{array}{rl}\\displaystyle{\\min_{\\theta,\\boldsymbol{\\sigma}} J_p(\\theta,\\boldsymbol{\\sigma})=\\min_{\\theta,\\boldsymbol{\\sigma}}} & \\displaystyle{ \\int_D \\boldsymbol{\\sigma}:(\\theta(x)^{-p}\\mathbb{C}^{-1}):\\boldsymbol{\\sigma}\\text{ dx}}\\\\\n", "\\text{s.t.} & \\text{div}\\, \\boldsymbol{\\sigma} + \\boldsymbol{f}=0 \\quad\\text{in } D\\\\\n", "& \\boldsymbol{\\sigma}\\cdot\\boldsymbol{n} = \\boldsymbol{T} \\quad\\text{on }\\Gamma_N \\\\\n", "& \\int_D \\theta(x) \\text{ dx} = \\eta \\text{Vol}(D) \\\\\n", "& 0 \\leq \\theta(x) \\leq 1\n", "\\end{array}$$\n", "\n", "### Optimization of the SIMP formulation\n", "\n", "Unfortunately jointly minimizing $J_p$ over $(\\theta,\\boldsymbol{\\sigma})$ is hard to do in practice since this functional is non-convex (except if $p=1$). However, it is convex over each variable $\\theta$ and $\\boldsymbol{\\sigma}$ when fixing the other variable. This makes alternate minimization an attractive method for finding a local optimum. Besides, minimizing directly for a fixed value of $p>1$ does not work well in practice. A better solution consists in performing alternate minimization steps and progressively increase $p$ from 1 to a maximum value (typically 3 or 4) using some heuristic (see later).\n", "\n", "Iteration $n+1$ of the algorithm, knowing a previous pair $(\\theta_n,\\boldsymbol{\\sigma}_n)$ and a given value of the penalty exponent $p_{n}$, therefore consists in:\n", "\n", "* minimizing $J_{p_n}(\\theta_n,\\boldsymbol{\\sigma})$ yielding $\\boldsymbol{\\sigma}_{n+1}=(\\theta_n)^{p_n}\\mathbb{C}:\\nabla^s\\boldsymbol{u}_{n+1}$\n", "* minimizing $J_{p_n}(\\theta, \\boldsymbol{\\sigma}_{n+1})$ yielding $\\theta_{n+1}$\n", "* potentially update the value of the exponent using some heuristic $p_n \\to p_{n+1}$\n", "\n", "Both minimization problems are quite easy to solve since the first one is nothing else than a heterogeneous elastic problem with some known modulus $(\\theta_n(x))^{p_n}\\mathbb{C}$.\n", "\n", "The second problem can be rewritten using the Lagrange multiplier $\\lambda$ corresponding to the toal volume constraint, as follows:\n", "\n", "$$\\begin{equation*}\n", "\\min_{\\theta\\in[0;1]} \\: \\int_D \\theta(x)^{-p_n}e(\\boldsymbol{\\sigma}_{n+1})\\text{ dx} + \\lambda \\left( \\int_D \\theta(x) \\text{ dx} -\\eta \\text{Vol}(D)\\right)\n", "\\end{equation*}$$\n", "\n", "where $e(\\sigma_{n+1}) = \\boldsymbol{\\sigma}_{n+1}:\\mathbb{C}^{-1}:\\boldsymbol{\\sigma}_{n+1}= (\\theta_n)^{p_n} \\boldsymbol{\\sigma}_{n+1}:\\nabla^s \\boldsymbol{u}_{n+1}$. The optimality conditions for this problem yields the following explicit and local condition:\n", "\n", "$$\\begin{equation*}\n", "-p_n \\theta_{n+1}^{-p_n-1} e(\\sigma_{n+1}) + \\lambda = 0\n", "\\end{equation*}$$\n", "\n", "which along with the $[\\theta_{min};1]$ constraint gives:\n", "\n", "$$\\begin{equation*}\n", "\\theta_{n+1} = \\min \\left\\lbrace 1; \\max \\left\\lbrace \\theta_{min}; \\left(\\frac{p_n e(\\sigma_{n+1})}{\\lambda}\\right)^{1/(p_n+1)} \\right\\rbrace\\right\\rbrace\n", "\\end{equation*}$$\n", "\n", "where we replaced the $0$ constraint by a minimum density value $\\theta_{min}>0$ to avoid degeneracy issue with void material.\n", "\n", "Note that the above expression assumes that $\\lambda$ is known. Its value is found by satsifying the volume constraint $\\int_D \\theta(x) \\text{ dx} =\\eta \\text{Vol}(D)$.\n", "\n", "## FEniCS implementation\n", "\n", "We first define some parameters of the algorithm. The most important ones concern the number of total alternate minimizations (AM) `niter` and the parameters controlling the exponent update strategy:\n", "\n", "* `niternp` corresponds to the number of first (AM) for which $p=1$. These are non-penalized iterations yielding a diffuse gray-level field $\\theta(x)$ at convergence (note that we solve this convex problem with AM iterations although one could use a dedicated convex optimization solver).\n", "* `pmax` is the maximum value taken by the penalty exponent $p$\n", "* `exponent_update_frequency` corresponds to the minimum number of AM iterations between two consecutive updates of the exponent" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "%matplotlib notebook\n", "from dolfin import *\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# Algorithmic parameters\n", "niternp = 20 # number of non-penalized iterations\n", "niter = 80 # total number of iterations\n", "pmax = 4 # maximum SIMP exponent\n", "exponent_update_frequency = 4 # minimum number of steps between exponent update\n", "tol_mass = 1e-4 # tolerance on mass when finding Lagrange multiplier\n", "thetamin = 0.001 # minimum density modeling void" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now define the problem parameters, in particular the target material density $(\\eta=40\\%$ here). The mesh consists of a rectangle of dimension $4\\times 1$, clamped on its left side and loaded by a uniformly distributed vertical force on a line of length 0.1 centered around the center of the right side.\n", "\n", "Finally, we initialized the SIMP penalty exponent to $p=1$ and initialized also the density field and the Lagrange multiplier $\\lambda$." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Problem parameters\n", "thetamoy = 0.4 # target average material density\n", "E = Constant(1)\n", "nu = Constant(0.3)\n", "lamda = E*nu/(1+nu)/(1-2*nu)\n", "mu = E/(2*(1+nu))\n", "f = Constant((0, -1)) # vertical downwards force\n", "\n", "# Mesh\n", "mesh = RectangleMesh(Point(-2, 0), Point(2, 1), 50, 30, \"crossed\")\n", "# Boundaries\n", "def left(x, on_boundary):\n", " return near(x[0], -2) and on_boundary\n", "def load(x, on_boundary):\n", " return near(x[0], 2) and near(x[1], 0.5, 0.05)\n", "facets = MeshFunction(\"size_t\", mesh, 1)\n", "AutoSubDomain(load).mark(facets, 1)\n", "ds = Measure(\"ds\", subdomain_data=facets)\n", "\n", "# Function space for density field\n", "V0 = FunctionSpace(mesh, \"DG\", 0)\n", "# Function space for displacement\n", "V2 = VectorFunctionSpace(mesh, \"CG\", 2)\n", "# Fixed boundary condtions\n", "bc = DirichletBC(V2, Constant((0, 0)), left)\n", "\n", "p = Constant(1) # SIMP penalty exponent\n", "exponent_counter = 0 # exponent update counter\n", "lagrange = Constant(1) # Lagrange multiplier for volume constraint\n", "\n", "thetaold = Function(V0, name=\"Density\")\n", "thetaold.interpolate(Constant(thetamoy))\n", "coeff = thetaold**p\n", "theta = Function(V0)\n", "\n", "volume = assemble(Constant(1.)*dx(domain=mesh))\n", "avg_density_0 = assemble(thetaold*dx)/volume # initial average density\n", "avg_density = 0." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now define some useful functions for formulating the linear elastic variational problem." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def eps(v):\n", " return sym(grad(v))\n", "def sigma(v):\n", " return coeff*(lamda*div(v)*Identity(2)+2*mu*eps(v))\n", "def energy_density(u, v):\n", " return inner(sigma(u), eps(v))\n", "\n", "# Inhomogeneous elastic variational problem\n", "u_ = TestFunction(V2)\n", "du = TrialFunction(V2)\n", "a = inner(sigma(u_), eps(du))*dx\n", "L = dot(f, u_)*ds(1)" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "We now define the function for updating the value of :math:`\\theta`. Note that this function will be called many times in each step since we are finding :math:`\\lambda` through a dichotomy procedure. For this reason, we make use of the `local_project` function (see :ref:`vonMisesPlasticity` and :ref:`TipsTricksProjection`) for fast projection on local spaces such as DG0 for the present case." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def local_project(v, V):\n", " dv = TrialFunction(V)\n", " v_ = TestFunction(V)\n", " a_proj = inner(dv, v_)*dx\n", " b_proj = inner(v, v_)*dx\n", " solver = LocalSolver(a_proj, b_proj)\n", " solver.factorize()\n", " u = Function(V)\n", " solver.solve_local_rhs(u)\n", " return u\n", "def update_theta():\n", " theta.assign(local_project((p*coeff*energy_density(u, u)/lagrange)**(1/(p+1)), V0))\n", " thetav = theta.vector().get_local()\n", " theta.vector().set_local(np.maximum(np.minimum(1, thetav), thetamin))\n", " theta.vector().apply(\"insert\")\n", " avg_density = assemble(theta*dx)/volume\n", " return avg_density" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now define a function for finding the correct value of the Lagrange multiplier $\\lambda$. First, a rough bracketing of $\\lambda$ is obtained, then a dichotomy is performed in the interval `[lagmin,lagmax]` until the correct average density is obtained to a certain tolerance." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "def update_lagrange_multiplier(avg_density):\n", " avg_density1 = avg_density\n", " # Initial bracketing of Lagrange multiplier\n", " if (avg_density1 < avg_density_0):\n", " lagmin = float(lagrange)\n", " while (avg_density < avg_density_0):\n", " lagrange.assign(Constant(lagrange/2))\n", " avg_density = update_theta()\n", " lagmax = float(lagrange)\n", " elif (avg_density1 > avg_density_0):\n", " lagmax = float(lagrange)\n", " while (avg_density > avg_density_0):\n", " lagrange.assign(Constant(lagrange*2))\n", " avg_density = update_theta()\n", " lagmin = float(lagrange)\n", " else:\n", " lagmin = float(lagrange)\n", " lagmax = float(lagrange)\n", "\n", " # Dichotomy on Lagrange multiplier\n", " inddico=0\n", " while ((abs(1.-avg_density/avg_density_0)) > tol_mass):\n", " lagrange.assign(Constant((lagmax+lagmin)/2))\n", " avg_density = update_theta()\n", " inddico += 1;\n", " if (avg_density < avg_density_0):\n", " lagmin = float(lagrange)\n", " else:\n", " lagmax = float(lagrange)\n", " print(\" Dichotomy iterations:\", inddico)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, the exponent update strategy is implemented:\n", "\n", "* first, $p=1$ for the first `niternp` iterations\n", "* then, $p$ is increased by some amount which depends on the average gray level of the density field computed as $g = \\frac{1}{\\text{Vol(D)}}\\int_D 4(\\theta-\\theta_{min})(1-\\theta)\\text{ dx}$, that is $g=0$ is $\\theta(x)=\\theta_{min}$ or $1$ everywhere and $g=1$ is $\\theta=(\\theta_{min}+1)/2$ everywhere.\n", "* Note that $p$ can increase only if at least `exponent_update_frequency` AM iterations have been performed since the last update and only if the compliance evolution falls below a certain threshold." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "def update_exponent(exponent_counter):\n", " exponent_counter += 1\n", " if (i < niternp):\n", " p.assign(Constant(1))\n", " elif (i >= niternp):\n", " if i == niternp:\n", " print(\"\\n Starting penalized iterations\\n\")\n", " if ((abs(compliance-old_compliance) < 0.01*compliance_history[0]) and \n", " (exponent_counter > exponent_update_frequency) ):\n", " # average gray level\n", " gray_level = assemble((theta-thetamin)*(1.-theta)*dx)*4/volume\n", " p.assign(Constant(min(float(p)*(1+0.3**(1.+gray_level/2)), pmax)))\n", " exponent_counter = 0\n", " print(\" Updated SIMP exponent p = \", float(p))\n", " return exponent_counter\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, the global loop for the algorithm is implemented consisting, at each iteration, of the elasticity problem resolution, the corresponding compliance computation, the update for $\\theta$ and its associated Lagrange multiplier $\\lambda$ and the exponent update procedure." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "u = Function(V2, name=\"Displacement\")\n", "old_compliance = 1e30\n", "ffile = XDMFFile(\"topology_optimization.xdmf\")\n", "ffile.parameters[\"flush_output\"]=True\n", "ffile.parameters[\"functions_share_mesh\"]=True\n", "compliance_history = []\n", "for i in range(niter):\n", " solve(a == L, u, bc, solver_parameters={\"linear_solver\": \"cg\", \"preconditioner\": \"hypre_amg\"})\n", " ffile.write(thetaold, i)\n", " ffile.write(u, i)\n", "\n", " compliance = assemble(action(L, u))\n", " compliance_history.append(compliance)\n", " print(\"Iteration {}: compliance =\".format(i), compliance)\n", "\n", " avg_density = update_theta()\n", "\n", " update_lagrange_multiplier(avg_density)\n", "\n", " exponent_counter = update_exponent(exponent_counter)\n", "\n", " # Update theta field and compliance\n", " thetaold.assign(theta)\n", " old_compliance = compliance" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The final density is represented as well as the convergence history of the compliance. One can note that the final compliance obtained after the first non-penalized iterations is smaller than the final one. This initially obtained topology is therefore more optimal than the final one, although it is made of large diffuse gray regions (see XMDF outputs or the animation at the beginning of the tour) contrary to the final one which is close to being binary." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "mpl.figure.prototype._root_extra_style = function(el){\n", " var fig = this\n", " el.on(\"remove\", function(){\n", "\tfig.close_ws(fig, {});\n", " });\n", "}\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " event.shiftKey = false;\n", " // Send a \"J\" for go to next cell\n", " event.which = 74;\n", " event.keyCode = 74;\n", " manager.command_mode();\n", " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot(theta, cmap=\"bone_r\")\n", "plt.title(\"Final density\")\n", "plt.show();\n", "\n", "plt.figure()\n", "plt.plot(np.arange(1, niter+1), compliance_history)\n", "ax = plt.gca()\n", "ymax = ax.get_ylim()[1]\n", "plt.plot([niternp, niternp], [0, ymax], \"--k\")\n", "plt.annotate(r\"$\\leftarrow$ Penalized iterations $\\rightarrow$\", xy=[niternp+1, ymax*0.02], fontsize=14)\n", "plt.xlabel(\"Number of iterations\")\n", "plt.ylabel(\"Compliance\")\n", "plt.title(\"Convergence history\", fontsize=16)\n", "plt.show();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "\n", "[ALL07]: Allaire, Grégoire, and Marc Schoenauer. *Conception optimale de structures.* Vol. 58. Berlin: Springer, 2007.\n", "\n", "[BEN03]: Bendsoe, Martin Philip and Sigmund, Ole. *Topology Optimization: Theory, Methods, and Applications*. Springer Science & Business Media, 2003." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Format de la Cellule Texte Brut", "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.8.10" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": false, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }