PS2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fast PDE/IE course, Skoltech, Spring 2015"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Set 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NOTE:** You can use the code that was presented in lectures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 1 (Image processing) (40pt)\n",
    "\n",
    "Let us consider a problem of de-noising a color image, by solving the following PDE:\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{\\partial}{\\partial t} u _i &= {\\rm div} (\\nabla u _i) + \\lambda (u^{(0)}-u _i) \\\\\n",
    "u _i &= u^{(0)} _i \\\\\n",
    "u _i| _{\\partial\\Omega} &= u^{(0)} _i| _{\\partial\\Omega} \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "where $u_i$, $i=1,2,3$, are the intensity of the red, green, and blue components of the image;\n",
    "$u^{(0)}_i$ represent the original image.\n",
    "The parameter $\\lambda$ should be inversely proportional to the level of the noise, and can be taken, for a start, to be 1/100 or 1/200.\n",
    "The time of integration should be taken so that we can approximate a steady-state solution.\n",
    "\n",
    "* Take a noisy image from <a href=\"https://github.com/oseledets/fastpde/blob/master/fig/color_photo_noisy.png\">here</a> or generate it yourself. (If you generate yourself then make the noise level to be 64 out of 256.)\n",
    "\n",
    "* **(10pt)** Write a code that solves the system\n",
    "* **(10pt)** Set $\\lambda=1/200$ and try getting some results with $\\tau=0.025$, $\\tau=0.25$, and  $\\tau=2.5$. Comment on, and explain the results.\n",
    "\n",
    "Let us consider a more advanced approach:\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{\\partial}{\\partial t} u _i &= {\\rm div} (k \\nabla u _i) + \\lambda (u^{(0)}-u _i) \\\\\n",
    "u _i &= u^{(0)} _i \\\\\n",
    "u _i| _{\\partial\\Omega} &= u^{(0)} _i| _{\\partial\\Omega} \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "where now $k=\\big(1 + \\max_i |\\nabla u_i|^2\\big)^{-1/2}$.\n",
    "\n",
    "* **(10pt)** Write a code that solves the system\n",
    "* **(10pt)** Give results for three values of $\\lambda$, for which the final image is (1) Too noisy, (2) Too blurry, (3) About right"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 2 (Elasticity) (60pt)\n",
    "\n",
    "Consider a domain $\\Omega$ with the boundary consisting of lines between the following 12 points (given in ${\\rm mm}$ in this very order):\n",
    "```\n",
    "(250, -2 S - H), (100 - H, -2 S - H), (100 - H, -S), (150 - H, -S), (150 - H, -H), (0, -H), (0, 0), (150, 0), (150, -S - H), (100, -S - H), (100, -2 S), (250, -2 S)\n",
    "```\n",
    "for some values of $S$ and $H$. The domain is as show here:\n",
    "<img src=\"fig/2Delastic_domain.png\">\n",
    "\n",
    "Next, we assume a piece of elastic material (similar to a beam) of this shape. It is fixed at the leftmost end and the force $P$ is applied to its rightmost end (the illustration below is for different values of $H$ and $S$):\n",
    "<img src=\"fig/2Delastic_domain_ext.png\">\n",
    "We also assume that the material that thickness (width) $W$ in the third dimensions (perpendicular to the plane of the image)\n",
    "\n",
    "To solve this problem, we consider the linear elasticity equations. The elasticity equations describe the displacement $(u,v)$ of the material in the $x$ and $y$ coordiate, respectively. We split the boudary of $\\Omega$ into three parts: $\\Gamma_1 = \\{(x,y)\\in\\partial\\Omega : x = 0\\}$ (leftmost end), $\\Gamma_2 = \\{(x,y)\\in\\partial\\Omega : 0 < x < 250\\}$ (middle part), and $\\Gamma_3 = \\{(x,y)\\in\\partial\\Omega : x = 250\\}$ (rightmost part). It is then easy to write the boundary conditions on $\\Gamma_1$:\n",
    "$$\n",
    "u=v=0\\qquad\\text{on }\\Gamma_1\n",
    "$$\n",
    "\n",
    "To write down the equations and the rest of the boundary conditions, we need the elasticity tensors which are denoted as\n",
    "$$\n",
    "\\begin{align*}\n",
    "C_{11} &= \\begin{pmatrix}\n",
    "\\frac{E (1-\\nu)}{(1-2\\nu)(1+\\nu)} & 0 \\\\ 0 & \\frac{E}{2(1+\\nu)} \\end{pmatrix},\n",
    "&\\qquad\n",
    "C_{12} &= \\begin{pmatrix} 0 & \\frac{E \\nu}{(1-2\\nu)(1+\\nu)}  \\\\ \\frac{E}{2(1+\\nu)} & 0 \\end{pmatrix},\n",
    "\\\\\n",
    "C_{21} &= C_{12}^{\\rm T},\n",
    "&\\qquad\n",
    "C_{22} &= \\begin{pmatrix} \\frac{E}{2(1+\\nu)} & 0 \\\\ 0 & \\frac{E (1-\\nu)}{(1-2\\nu)(1+\\nu)} \\end{pmatrix}\n",
    "\\end{align*}\n",
    "$$\n",
    "Here $E$ is called the modulus of elasticity and $\\nu$ is called the Poisson ratio.\n",
    "\n",
    "The equations are then written as\n",
    "$$\n",
    "\\begin{align*}\n",
    "-{\\rm div} (C_{11} \\nabla u + C_{12} \\nabla v) &= 0\n",
    "\\qquad\\text{in }\\Omega\n",
    "\\\\\n",
    "-{\\rm div} (C_{21} \\nabla u + C_{22} \\nabla v) &= 0\n",
    "\\qquad\\text{in }\\Omega\n",
    "\\\\\n",
    "u = v &= 0 \\qquad\\text{on }\\Gamma_1\n",
    "\\\\\n",
    "C_{11} u_n + C_{12} v_n = C_{21} u_n + C_{22} v_n &= 0 \\qquad\\text{on }\\Gamma_2\n",
    "\\\\\n",
    "C_{11} u_n + C_{12} v_n &= 0 \\qquad\\text{on }\\Gamma_3\n",
    "\\\\\n",
    "C_{21} u_n + C_{22} v_n &= \\frac{P}{H W} \\qquad\\text{on }\\Gamma_3,\n",
    "\\end{align*}\n",
    "$$\n",
    "where $u_n$ and $v_n$ are normal derivatives of the solution.\n",
    "Note that in the last equation $C_{21} \\nabla u + C_{22} \\nabla v$ is the stress on the boundary and we set it equal to the stress $\\frac{P}{H W}$ of the external force $P$. (The stress has the same dimension as the pressure.)\n",
    "\n",
    "Let us now convert this to the variational formulation.\n",
    "We denote the test functions by $p=p(x,y)$ and $q = q(x,y)$.\n",
    "We require that the test functions satisfy the Dirichlet boundary conditions (i.e., $p(0,y)=q(0,y)=0$) and then the weak formulation becomes\n",
    "$$\n",
    "\\int_{\\Omega} \\big[\n",
    "    (C_{11} \\nabla u + C_{12} \\nabla v)\\cdot \\nabla p\n",
    "    +\n",
    "    (C_{21} \\nabla u + C_{22} \\nabla v)\\cdot \\nabla q\n",
    "\\big]\n",
    "= \\frac{P}{H W}  \\int_{\\Gamma_3} q\n",
    "\\qquad\\forall p, q\n",
    "$$\n",
    "\n",
    "* **(15pt)** Write a code that generates a mesh for this domain (you may, but do not need to, use mesh generation packages)\n",
    "\n",
    "* **(35pt)** Write a code that solves this problem by the finite element method and outputs the $y$-component of the deflection of this piece of material, to be precise, $v(250,2 S)$\n",
    "\n",
    "* **(10pt)** For three sets of parameters, $S=20$, $W=3$, $H\\in\\{7,8,9\\}$ (all in ${\\rm mm}$) compute the displacement $v(250,2 S)$, $E = 3.1\\cdot 10^9 {\\rm Pa}$, $\\nu=0.35$, $P = m\\cdot g$, $m=200 {\\rm grams}$, $g=-9810 {\\rm mm}/{\\rm s}^2$.\n",
    "\n",
    "* **(Optional)** Plot the three solutions found. For that you need to change the coordinates of each of the mesh nodes from $(x,y)$ by $(x,y) + (u(x,y),v(x,y))$ and plot the displaced mesh\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "breakpoint": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 100%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 100%;\n",
       "\t\t\tline-height: 110%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 110%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/alex.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "breakpoint": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}