Mircrostructure Optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Microstructure Optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "using Plots\n",
    "import JSON\n",
    "using SparseArrays\n",
    "using Images\n",
    "using StaticArrays, BenchmarkTools\n",
    "using VectorizedRoutines\n",
    "using NLopt\n",
    "using Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "getIndex3d (generic function with 2 methods)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "include(\"./julia/FEM.jl\")\n",
    "include(\"./julia/data_FEM.jl\")\n",
    "include(\"./julia/MBB.jl\")\n",
    "include(\"./julia/topologyOptimization.jl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "densityOptimization (generic function with 3 methods)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function densityOptimization(data_FEM,nelx,nely,totalVolFactor,maxEval,SIMP,sensitivity,p=3,rmin=1.5)\n",
    "\n",
    "    problem=data_FEM(nelx,nely)\n",
    "    E,f,g,idb,ien,iplane,ndf,nel,nen,nnp,nsd,snu,t,xn,Ke=problem;\n",
    "    ρ=ones(nel);\n",
    "    K,F,d,de,Fe=FEM(problem,ρ);\n",
    "    debug=false\n",
    "    if debug\n",
    "        display(\"Degrees of freedom\")\n",
    "        display(heatmap(reshape(idb,2,nelx+1,nely+1)[1,:,:]',c=:YlGnBu,aspect_ratio=:equal))\n",
    "        display(heatmap(reshape(idb,2,nelx+1,nely+1)[2,:,:]',c=:YlGnBu,aspect_ratio=:equal))\n",
    "        display(\"Forces\")\n",
    "        display(heatmap(reshape(f,2,nelx+1,nely+1)[1,:,:]',c=:YlGnBu,aspect_ratio=:equal))\n",
    "        display(heatmap(reshape(f,2,nelx+1,nely+1)[2,:,:]',c=:YlGnBu,aspect_ratio=:equal))\n",
    "        display(\"Initial forces\")\n",
    "        dd=zeros(length(idb[:]))\n",
    "        dd[Bool.(abs.(idb[:] .-1))]=d\n",
    "        display(heatmap(reshape(dd,2,nelx+1,nely+1)[1,:,:]',c=:YlGnBu))\n",
    "        display(heatmap(reshape(dd,2,nelx+1,nely+1)[2,:,:]',c=:YlGnBu))\n",
    "    end\n",
    "    \n",
    "    if !SIMP\n",
    "        X=optimizeDensity(problem,totalVolFactor,maxEval);\n",
    "    else\n",
    "        if sensitivity\n",
    "            X=optimizeDensitySIMPSensitivity(problem,nelx,nely,totalVolFactor,maxEval,p,rmin)\n",
    "        else\n",
    "            X=optimizeDensitySIMP(problem,totalVolFactor,maxEval,p);\n",
    "        end\n",
    "    end\n",
    "    display(heatmap(reshape(X,nelx,nely)',c=:YlGnBu,aspect_ratio=:equal))\n",
    "    return X\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"Minimum compliance problem with OC\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"ndes: 90 x 30\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"volfrac: 0.5 rmin: 1.5 penal: 3\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\" It:1 Obj:1018.7549500315454 Vol:0.5000274820973205 ch:0.2 \""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip6600\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip6600)\" d=\"\n",
       "M0 1600 L2400 1600 L2400 0 L0 0  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip6601\">\n",
       "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip6600)\" d=\"\n",
       "M86.9921 1521.01 L2352.76 1521.01 L2352.76 47.2441 L86.9921 47.2441  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip6602\">\n",
       "    <rect x=\"86\" y=\"47\" width=\"2267\" height=\"1475\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip6600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  86.9921,1521.01 2352.76,1521.01 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  86.9921,1521.01 86.9921,47.2441 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip6602)\">\n",
       "<image width=\"2266\" height=\"755\" xlink:href=\"data:image/png;base64,\n",
       "iVBORw0KGgoAAAANSUhEUgAACNoAAALzCAYAAADN4vhJAAAgAElEQVR4nOzcTW4kR5aFUfMfJgj9\n",
       "zARI69L+FyBAE0HIFEEyIrwGNWlUA40bXWkvnI/nbOCGmQeTTOaXsfz+++/HAAAAAACAD2BZlukb\n",
       "27ZN33h6epq+se/79I0x+pylYqPLXVV8jYzR5yw2cuu6Tt8Yo+Z7ScVGheOY/0/JFRtjjHG9Xqdv\n",
       "vL+/T994e3ubvvHy8jJ949u3by02xhjjcrlM36h4/1aco4ua71YAAAAAAAAAAPDBCW0AAAAAAAAA\n",
       "ACAgtAEAAAAAAAAAgIDQBgAAAAAAAAAAAkIbAAAAAAAAAAAICG0AAAAAAAAAACAgtAEAAAAAAAAA\n",
       "gIDQBgAAAAAAAAAAAkIbAAAAAAAAAAAICG0AAAAAAAAAACAgtAEAAAAAAAAAgIDQBgAAAAAAAAAA\n",
       "AkIbAAAAAAAAAAAICG0AAAAAAAAAACAgtAEAAAAAAAAAgIDQBgAAAAAAAAAAAkIbAAAAAAAAAAAI\n",
       "CG0AAAAAAAAAACAgtAEAAAAAAAAAgIDQBgAAAAAAAAAAAkIbAAAAAAAAAAAICG0AAAAAAAAAACAg\n",
       "tAEAAAAAAAAAgIDQBgAAAAAAAAAAAkIbAAAAAAAAAAAI7I9+AQAAAAAAzLUsS8nOus7/v537Pv/X\n",
       "2l02np6epm+M0ee+bJxro2pn2zYboYo/4yvOMUbNWSq+91Z9f+/iOI7pG9frtcXGzz//PH3j/f19\n",
       "+sbr6+v0jTHGeHt7m77x7du36Rt//vnn9I2K517BJ9oAAAAAAAAAAEBAaAMAAAAAAAAAAAGhDQAA\n",
       "AAAAAAAABIQ2AAAAAAAAAAAQENoAAAAAAAAAAEBAaAMAAAAAAAAAAAGhDQAAAAAAAAAABIQ2AAAA\n",
       "AAAAAAAQENoAAAAAAAAAAEBAaAMAAAAAAAAAAAGhDQAAAAAAAAAABIQ2AAAAAAAAAAAQENoAAAAA\n",
       "AAAAAEBAaAMAAAAAAAAAAAGhDQAAAAAAAAAABIQ2AAAAAAAAAAAQENoAAAAAAAAAAEBAaAMAAAAA\n",
       "AAAAAAGhDQAAAAAAAAAABIQ2AAAAAAAAAAAQENoAAAAAAAAAAEBAaAMAAAAAAAAAAAGhDQAAAAAA\n",
       "AAAABIQ2AAAAAAAAAAAQ2B/9AgAAAAAA/j+WZWmxsW3b9I19r/lVcMVZnp6epm9U3FeXjaodG+fa\n",
       "qPhaH6PPWbr8OV9xjnWt+YyALmep2Ojy81bVTtVZZutyjiqXy2X6xh9//DF94++//56+8f7+Pn2j\n",
       "gk+0AQAAAAAAAACAgNAGAAAAAAAAAAACQhsAAAAAAAAAAAgIbQAAAAAAAAAAICC0AQAAAAAAAACA\n",
       "gNAGAAAAAAAAAAACQhsAAAAAAAAAAAgIbQAAAAAAAAAAICC0AQAAAAAAAACAgNAGAAAAAAAAAAAC\n",
       "QhsAAAAAAAAAAAgIbQAAAAAAAAAAICC0AQAAAAAAAACAgNAGAAAAAAAAAAACQhsAAAAAAAAAAAgI\n",
       "bQAAAAAAAAAAICC0AQAAAAAAAACAgNAGAAAAAAAAAAACQhsAAAAAAAAAAAgIbQAAAAAAAAAAICC0\n",
       "AQAAAAAAAACAgNAGAAAAAAAAAAACQhsAAAAAAAAAAAgIbQAAAAAAAAAAILA/+gUAAAAAAP+2LEub\n",
       "nW3bbIT2ff6vaSs2qnZs5Crev2P0ua8uz6TLXY3R5758v8pV/blVsbOu8z/voGKj4ufGqp+Bu3Bf\n",