{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Manufactured solutions for Elasticity\n",
    "\n",
    "This is a notebook to play with manufactured solutions for Elasticity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division\n",
    "from sympy import *\n",
    "x, y, z, t = symbols('x y z t')\n",
    "f, g, h = symbols('f g h', cls=Function)\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = symbols('L')\n",
    "a1, a2, a3, b1, b2, b3, c1, c2, c3 = symbols('a1 a2 a3 b1 b2 b3 c1 c2 c3')\n",
    "u0, ux, uy, uz, v0, vx, vy, vz, w0, wx, wy, wz = symbols('u0 u_x u_y u_z\\\n",
    "                                                         v0 v_x v_y v_z\\\n",
    "                                                         w0 w_x w_y w_z')\n",
    "lamda, mu, rho = symbols('lamda mu rho')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = u0 + ux*sin(a1*pi*x/L) + uy*sin(a2*pi*y/L) + uz*sin(a3*pi*z/L)\n",
    "v = v0 + vx*sin(b1*pi*x/L) + vy*sin(b2*pi*y/L) + vz*sin(b3*pi*z/L)\n",
    "w = w0 + wx*sin(c1*pi*x/L) + wy*sin(c2*pi*y/L) + wz*sin(c3*pi*z/L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def laplacian(U, X=[x, y, z]):\n",
    "    lap_U = Matrix([sum([diff(U[k], X[j], 2) for j in range(3)]) for k in range(3)])\n",
    "    return lap_U  \n",
    "\n",
    "\n",
    "def div(U, X=[x, y, z]):\n",
    "    return sum([diff(U[k], X[k]) for k in range(3)])\n",
    "\n",
    "\n",
    "def grad(f, X=[x, y, z]):\n",
    "    return Matrix([diff(f, X[k]) for k in range(3)])\n",
    "\n",
    "\n",
    "def rot(U, X=[x, y, z]):\n",
    "    return Matrix([diff(U[2], X[1]) - diff(U[1], X[2]),\n",
    "                   diff(U[0], X[2]) - diff(U[2], X[0]),\n",
    "                   diff(U[1], X[0]) - diff(U[0], X[1])])\n",
    "\n",
    "def navier(U, X=[x,y,z], lamda=lamda, mu=mu):\n",
    "    return mu*laplacian(U, X) + (lamda + mu)*grad(div(U, X), X)\n",
    "\n",
    "\n",
    "def dt(U, t=t):\n",
    "    return Matrix([diff(U[k], t) for k in range(3)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}\\frac{\\pi^{2}}{L^{2}} \\left(a_{1}^{2} u_{x} \\left(\\lambda + \\mu\\right) \\sin{\\left (\\frac{\\pi a_{1}}{L} x \\right )} + \\mu \\left(a_{1}^{2} u_{x} \\sin{\\left (\\frac{\\pi a_{1}}{L} x \\right )} + a_{2}^{2} u_{y} \\sin{\\left (\\frac{\\pi a_{2}}{L} y \\right )} + a_{3}^{2} u_{z} \\sin{\\left (\\frac{\\pi a_{3}}{L} z \\right )}\\right)\\right)\\\\\\frac{\\pi^{2}}{L^{2}} \\left(b_{2}^{2} v_{y} \\left(\\lambda + \\mu\\right) \\sin{\\left (\\frac{\\pi b_{2}}{L} y \\right )} + \\mu \\left(b_{1}^{2} v_{x} \\sin{\\left (\\frac{\\pi b_{1}}{L} x \\right )} + b_{2}^{2} v_{y} \\sin{\\left (\\frac{\\pi b_{2}}{L} y \\right )} + b_{3}^{2} v_{z} \\sin{\\left (\\frac{\\pi b_{3}}{L} z \\right )}\\right)\\right)\\\\\\frac{\\pi^{2}}{L^{2}} \\left(c_{3}^{2} w_{z} \\left(\\lambda + \\mu\\right) \\sin{\\left (\\frac{\\pi c_{3}}{L} z \\right )} + \\mu \\left(c_{1}^{2} w_{x} \\sin{\\left (\\frac{\\pi c_{1}}{L} x \\right )} + c_{2}^{2} w_{y} \\sin{\\left (\\frac{\\pi c_{2}}{L} y \\right )} + c_{3}^{2} w_{z} \\sin{\\left (\\frac{\\pi c_{3}}{L} z \\right )}\\right)\\right)\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡ 2 ⎛  2               ⎛π⋅a₁⋅x⎞     ⎛  2       ⎛π⋅a₁⋅x⎞     2        ⎛π⋅a₂⋅y⎞ \n",
       "⎢π ⋅⎜a₁ ⋅uₓ⋅(λ + μ)⋅sin⎜──────⎟ + μ⋅⎜a₁ ⋅uₓ⋅sin⎜──────⎟ + a₂ ⋅u_y⋅sin⎜──────⎟ \n",
       "⎢   ⎝                  ⎝  L   ⎠     ⎝          ⎝  L   ⎠              ⎝  L   ⎠ \n",
       "⎢─────────────────────────────────────────────────────────────────────────────\n",
       "⎢                                                  2                          \n",
       "⎢                                                 L                           \n",
       "⎢                                                                             \n",
       "⎢ 2 ⎛  2                ⎛π⋅b₂⋅y⎞     ⎛  2       ⎛π⋅b₁⋅x⎞     2        ⎛π⋅b₂⋅y⎞\n",
       "⎢π ⋅⎜b₂ ⋅v_y⋅(λ + μ)⋅sin⎜──────⎟ + μ⋅⎜b₁ ⋅vₓ⋅sin⎜──────⎟ + b₂ ⋅v_y⋅sin⎜──────⎟\n",
       "⎢   ⎝                   ⎝  L   ⎠     ⎝          ⎝  L   ⎠              ⎝  L   ⎠\n",
       "⎢─────────────────────────────────────────────────────────────────────────────\n",
       "⎢                                                   2                         \n",
       "⎢                                                  L                          \n",
       "⎢                                                                             \n",
       "⎢ 2 ⎛  2                ⎛π⋅c₃⋅z⎞     ⎛  2       ⎛π⋅c₁⋅x⎞     2        ⎛π⋅c₂⋅y⎞\n",
       "⎢π ⋅⎜c₃ ⋅w_z⋅(λ + μ)⋅sin⎜──────⎟ + μ⋅⎜c₁ ⋅wₓ⋅sin⎜──────⎟ + c₂ ⋅w_y⋅sin⎜──────⎟\n",
       "⎢   ⎝                   ⎝  L   ⎠     ⎝          ⎝  L   ⎠              ⎝  L   ⎠\n",
       "⎢─────────────────────────────────────────────────────────────────────────────\n",
       "⎢                                                   2                         \n",
       "⎣                                                  L                          \n",
       "\n",
       "    2        ⎛π⋅a₃⋅z⎞⎞⎞ ⎤\n",
       "+ a₃ ⋅u_z⋅sin⎜──────⎟⎟⎟ ⎥\n",
       "             ⎝  L   ⎠⎠⎠ ⎥\n",
       "─────────────────────── ⎥\n",
       "                        ⎥\n",
       "                        ⎥\n",
       "                        ⎥\n",
       "     2        ⎛π⋅b₃⋅z⎞⎞⎞⎥\n",
       " + b₃ ⋅v_z⋅sin⎜──────⎟⎟⎟⎥\n",
       "              ⎝  L   ⎠⎠⎠⎥\n",
       "────────────────────────⎥\n",
       "                        ⎥\n",
       "                        ⎥\n",
       "                        ⎥\n",
       "     2        ⎛π⋅c₃⋅z⎞⎞⎞⎥\n",
       " + c₃ ⋅w_z⋅sin⎜──────⎟⎟⎟⎥\n",
       "              ⎝  L   ⎠⎠⎠⎥\n",
       "────────────────────────⎥\n",
       "                        ⎥\n",
       "                        ⎦"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = [u,v,w]\n",
    "F = rho*dt(U) - navier(U)\n",
    "\n",
    "simplify(F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Cosserat solid\n",
    "\n",
    "\n",
    "Let's start with the 2D case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = symbols('L')\n",
    "h0, hx, hy = symbols('h0 hx hy')\n",
    "lamda, mu, rho = symbols('lamda mu rho')\n",
    "alpha, mu_c, xi, gamma, J = symbols('alpha mu_c xi gamma J')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = u0 + ux*sin(a1*pi*x/L) + uy*sin(a2*pi*y/L)\n",
    "v = v0 + vx*sin(b1*pi*x/L) + vy*sin(b2*pi*y/L)\n",
    "w = 0\n",
    "h1 = 0\n",
    "h2 = 0\n",
    "h3 = h0 + hx*sin(c1*pi*x/L) + hy*sin(c2*pi*y/L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def coss(U, H, X=[x,y,z], lamda=lamda, mu=mu, mu_c=mu_c, xi=xi, gamma=gamma):\n",
    "    f_u = (lamda + 2*mu)*grad(div(U)) - (mu + mu_c)*rot(rot(U)) + 2*mu_c*rot(H)\n",
    "    f_h = (alpha + 2*gamma + 2*xi)*grad(div(H)) - 2*xi*rot(rot(H)) + 2*mu_c*rot(U) - 4*mu_c*Matrix(H)\n",
    "    return f_u, f_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = [u, v, w]\n",
    "H = [h1, h2, h3]\n",
    "f_u, f_h = coss(U, H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}\\frac{\\pi}{L^{2}} \\left(2 L c_{2} hy \\mu_{c} \\cos{\\left (\\frac{\\pi c_{2}}{L} y \\right )} - \\pi \\left(a_{1}^{2} u_{x} \\left(\\lambda + 2 \\mu\\right) \\sin{\\left (\\frac{\\pi a_{1}}{L} x \\right )} + a_{2}^{2} u_{y} \\left(\\mu + \\mu_{c}\\right) \\sin{\\left (\\frac{\\pi a_{2}}{L} y \\right )}\\right)\\right)\\\\- \\frac{\\pi}{L^{2}} \\left(2 L c_{1} hx \\mu_{c} \\cos{\\left (\\frac{\\pi c_{1}}{L} x \\right )} + \\pi \\left(b_{1}^{2} v_{x} \\left(\\mu + \\mu_{c}\\right) \\sin{\\left (\\frac{\\pi b_{1}}{L} x \\right )} + b_{2}^{2} v_{y} \\left(\\lambda + 2 \\mu\\right) \\sin{\\left (\\frac{\\pi b_{2}}{L} y \\right )}\\right)\\right)\\\\0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡   ⎛                 ⎛π⋅c₂⋅y⎞     ⎛  2                 ⎛π⋅a₁⋅x⎞     2        \n",
       "⎢ π⋅⎜2⋅L⋅c₂⋅hy⋅μ_c⋅cos⎜──────⎟ - π⋅⎜a₁ ⋅uₓ⋅(λ + 2⋅μ)⋅sin⎜──────⎟ + a₂ ⋅u_y⋅(μ \n",
       "⎢   ⎝                 ⎝  L   ⎠     ⎝                    ⎝  L   ⎠              \n",
       "⎢ ────────────────────────────────────────────────────────────────────────────\n",
       "⎢                                                 2                           \n",
       "⎢                                                L                            \n",
       "⎢                                                                             \n",
       "⎢   ⎛                 ⎛π⋅c₁⋅x⎞     ⎛  2                 ⎛π⋅b₁⋅x⎞     2        \n",
       "⎢-π⋅⎜2⋅L⋅c₁⋅hx⋅μ_c⋅cos⎜──────⎟ + π⋅⎜b₁ ⋅vₓ⋅(μ + μ_c)⋅sin⎜──────⎟ + b₂ ⋅v_y⋅(λ \n",
       "⎢   ⎝                 ⎝  L   ⎠     ⎝                    ⎝  L   ⎠              \n",
       "⎢─────────────────────────────────────────────────────────────────────────────\n",
       "⎢                                                 2                           \n",
       "⎢                                                L                            \n",
       "⎢                                                                             \n",
       "⎣                                                0                            \n",
       "\n",
       "          ⎛π⋅a₂⋅y⎞⎞⎞ ⎤\n",
       "+ μ_c)⋅sin⎜──────⎟⎟⎟ ⎥\n",
       "          ⎝  L   ⎠⎠⎠ ⎥\n",
       "──────────────────── ⎥\n",
       "                     ⎥\n",
       "                     ⎥\n",
       "                     ⎥\n",
       "          ⎛π⋅b₂⋅y⎞⎞⎞ ⎥\n",
       "+ 2⋅μ)⋅sin⎜──────⎟⎟⎟ ⎥\n",
       "          ⎝  L   ⎠⎠⎠ ⎥\n",
       "─────────────────────⎥\n",
       "                     ⎥\n",
       "                     ⎥\n",
       "                     ⎥\n",
       "                     ⎦"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(f_u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}0\\\\0\\\\- \\frac{1}{L^{2}} \\left(4 L^{2} \\mu_{c} \\left(h_{0} + hx \\sin{\\left (\\frac{\\pi c_{1}}{L} x \\right )} + hy \\sin{\\left (\\frac{\\pi c_{2}}{L} y \\right )}\\right) + 2 \\pi L \\mu_{c} \\left(a_{2} u_{y} \\cos{\\left (\\frac{\\pi a_{2}}{L} y \\right )} - b_{1} v_{x} \\cos{\\left (\\frac{\\pi b_{1}}{L} x \\right )}\\right) + 2 \\pi^{2} \\xi \\left(c_{1}^{2} hx \\sin{\\left (\\frac{\\pi c_{1}}{L} x \\right )} + c_{2}^{2} hy \\sin{\\left (\\frac{\\pi c_{2}}{L} y \\right )}\\right)\\right)\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡                                                                             \n",
       "⎢                                                                             \n",
       "⎢                                                                             \n",
       "⎢                                                                             \n",
       "⎢ ⎛   2     ⎛           ⎛π⋅c₁⋅x⎞         ⎛π⋅c₂⋅y⎞⎞             ⎛          ⎛π⋅a\n",
       "⎢-⎜4⋅L ⋅μ_c⋅⎜h₀ + hx⋅sin⎜──────⎟ + hy⋅sin⎜──────⎟⎟ + 2⋅π⋅L⋅μ_c⋅⎜a₂⋅u_y⋅cos⎜───\n",
       "⎢ ⎝         ⎝           ⎝  L   ⎠         ⎝  L   ⎠⎠             ⎝          ⎝  L\n",
       "⎢─────────────────────────────────────────────────────────────────────────────\n",
       "⎢                                                                             \n",
       "⎣                                                                             \n",
       "\n",
       "0                                                                             \n",
       "                                                                              \n",
       "0                                                                             \n",
       "                                                                              \n",
       "₂⋅y⎞            ⎛π⋅b₁⋅x⎞⎞      2   ⎛  2       ⎛π⋅c₁⋅x⎞     2       ⎛π⋅c₂⋅y⎞⎞⎞ \n",
       "───⎟ - b₁⋅vₓ⋅cos⎜──────⎟⎟ + 2⋅π ⋅ξ⋅⎜c₁ ⋅hx⋅sin⎜──────⎟ + c₂ ⋅hy⋅sin⎜──────⎟⎟⎟ \n",
       "   ⎠            ⎝  L   ⎠⎠          ⎝          ⎝  L   ⎠             ⎝  L   ⎠⎠⎠ \n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       " 2                                                                            \n",
       "L                                                                             \n",
       "\n",
       "⎤\n",
       "⎥\n",
       "⎥\n",
       "⎥\n",
       "⎥\n",
       "⎥\n",
       "⎥\n",
       "⎥\n",
       "⎥\n",
       "⎦"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(f_h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "- Malaya, Nicholas, et al. \"MASA: a library for verification using manufactured and analytical solutions.\" Engineering with Computers 29.4 (2013): 487-496."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<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",
       "/* Based on Lorena Barba template available at: https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css*/\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:800px;\n",
       "margin-left:16% !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",
       "font-family: 'Fenix', serif;\n",
       "margin-top:12px;\n",
       "margin-bottom: 3px;\n",
       "}\n",
       "h4{\n",
       "font-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: 135%;\n",
       "font-size: 120%;\n",
       "width:600px;\n",
       "margin-left:auto;\n",
       "margin-right:auto;\n",
       "}\n",
       ".CodeMirror{\n",
       "font-family: \"Source Code Pro\";\n",
       "font-size: 90%;\n",
       "}\n",
       "/* .prompt{\n",
       "display: None;\n",
       "}*/\n",
       ".text_cell_render h1 {\n",
       "font-weight: 200;\n",
       "font-size: 50pt;\n",
       "line-height: 100%;\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",
       ".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",
       "\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open('./styles/custom_barba.css', 'r').read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}