joule-2d/joule.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Joule heating model\n",
    "The thoretical information comes from [here](https://reference.wolfram.com/language/PDEModels/tutorial/Multiphysics/ModelCollection/JouleHeating.html). \n",
    "\n",
    "## Electroestatic model\n",
    "In this model the electric potential field $V(x,y,z)$ is assumed to be independent of the temperature $T$. That is, the electrical conductivity $G$ is kept at a constant value at all time. By Guass's law the stationary potential field satisfies Poisson's equation:\n",
    "\n",
    "$$\\nabla \\cdot (-G \\cdot \\nabla V)=0$$\n",
    "\n",
    "## Heat transfer model\n",
    "The heat equation is used to solve for the temperature field in a heat transfer model:\n",
    "\n",
    "$$\\rho C_p \\frac{dT}{dt}+\\rho C_p v\\cdot \\nabla T + \\nabla \\cdot(-k \\nabla T) = Q$$\n",
    " \n",
    "For a steady-state heat transfer model the transient term in the equation are set to zero. Since a solid is modeled, the internal velocity also vanishes and the heat equation simplifies to:\n",
    "\n",
    "$$ \\nabla \\cdot(-k \\nabla T) = Q $$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving the PDE Model\n",
    "In order to solve the PDE model, first the electroestatic model will be solved first to simulate/obtain the potential field $V$ of the wire. The heat transfer model is then constructed to show the heating effect of the electric current.\n",
    "\n",
    "## The electroestatic model\n",
    "There are two types of the boundary conditions involved in the electrostatics model. At both ends of the wire an electric potential difference $V_0=0.2$ is applied.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fenics import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f89569e9160>,\n",
       " <matplotlib.lines.Line2D at 0x7f89569e92e8>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAABJCAYAAAA9vFENAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHCpJREFUeJztnXtsXNd95z+Hj6E4FIfkDMnh25T1sh62Y0u2bLcpHHsLO01qB029dopFtoESx9lku2kehpwCeaFBWzQu0sJFF4Gd3U22eW0cLAw3hYu6zUOWZEl0HMuSLVu2LL7fFN/DGY1O/7j3jO8dDjmXc88lZ6TzAQjM48w9v3N57++c8zvf87tCSonBYDAYDGUbbYDBYDAYigPTIRgMBoMBMB2CwWAwGGxMh2AwGAwGwHQIBoPBYLAxHYLBYDAYANMhGAwGg8HGdAgGg8FgAEyHYDAYDAabio02YC00NjbK7u7ujTbDYDAYSoaenp5xKWWTl7Il1SF0d3dz8uTJNf/ua1/7mut9NBqloiL4po+Ojrreh8NhNm/evO71AjQ2NlJWFvyEMLvu2tpaqqur171eIQSNjY0IIda97vr6ekKh0LrXW1lZSX19/Ya0eaPuqU2bNhGJRAKvN5FIMDMz4/osFotRXl4eeN2qzV/5ylcK+r0Q4oLXsiXVIRTCpUuXln02OTmZeV1eXk5nZ2cgTiv74l1YWGBhYSHzPhaL0dTUpP0Gnpubc9UDMD4+nnkdCoXo6OigqqpKa72wvM2zs7PMzs5m3sfjcRoaGrS3ObteKSVjY2OZ9+FwmLa2NiorK7XWm6vuixcvut63t7cH4rSy602lUq42RyIRWlpatDstKeWyup33lBCCzs5OampqtNYLy9ucSCRIJBKZ99FolKamJu2Dn+np6WUdwsTEROZ1ZWUlHR0dbNq0SWu9kHuAFxRXfIdQUVFBZWUl+/fv54477qCvr4/e3l76+voYGhoinU7zzjvv0NjYSGdnJ11dXXR2dhKNRn07reeff54jR47wpS99iaGhoUy9fX19zM/PMzExwfz8PJ2dnZm/9vZ2305rfn6eb37zm7z//e9nx44drjaPjo6STCY5f/488Xg8U29XVxd1dXW+6gV4+umnGRwc5JFHHmFgYMDV5qWlJUZGRpibm8uc587OTlpbW307rYGBAZ588kk+8pGPEIvFXPWOj4+zsLDA22+/TVtbm6vNOpzWU089RSgU4oEHHqC/vz9T98DAAKlUioGBAebm5lz1Njc3+3Zar776Kk8//TSf+tSnkFK62nzx4kVmZmZYWFigvb09U29HR4eWwc/jjz/Ojh07uOuuuzJ19vX1MTg4SDqdpre3l2g0mvk/d3V1EYvFfN9Thw8f5vnnn+fQoUOMj49n6u3t7WVubo7JyUnm5+fp6OjI1Nve3u57xnbp0iW+8Y1vcPfdd3PDDTe47qnh4WFSqRTnz5+nubnZ5Ud0zNieffZZXn/9dV/H8MoV3yE42bx5M7t27WLXrl2ANaIaHBzM/GNfe+01fv3rXwNQU1PjctStra0FT4nLy8vp6Oigo6MDsEZYU1NT9Pb2Zup+8803ASgrK6O1tdXlPAoNMwkhqK+vp76+nuuvvx6ApaUll9N6+eWXOXHiBGCNKJ31xuPxgp1WZWUl3d3dqDUfNWJ3Oq3XXnsNsDpt5bTUX6FOSwhBLBYjFotx0003AdbMzOm0jh8/ztGjRwFrROlss58w06ZNm9i2bRvbtm0DIJ1OMzIykmnzhQsXePXVV4F3Z2lOp1XojE0IQXNzM/F4nFtuuQWwZmbOc/3CCy9w+PBhAJqamjL1dnV1+XJaNTU1XHfddVx33XWA5TgHBwcz9Z49e5aXX34ZgOrqate5bmtrK/ieKisro729nfb2dm677TaklFy8eNHlqH/+859nzk9LS4vLUfuZsUUiEfbs2cOePXsASCaT9Pf3Z9p86tQpenp6AMvnOAc/QczYdHJVdQjZVFZWcs0113DNNdcAltMaHx/PXFC9vb2Znrm8vHyZ0wqHwwXVK4QgGo0SjUZ5z3veA1hOy+moT548ybFjxwBoaGhw3Uh+wkxVVVVs3bqVrVu3AnD58mWX0+rr6+P06dOA22l1dnb6CjMpp9Xc3Mz+/fsBy2k5HfWRI0e4fPky8K7TUm32E2YKh8Ps3LmTnTt3ApbTcs7Y3nzzTX7zm98AllPPdlqFztjKy8tpa2ujra0t47Smp6dd5/oXv/hF5vyoGZty1H6cVm1t7TKnNTAw4Pofv/TSS4DltJxt9uO0KioqMvaDdU9NTEy4HPUbb7yROT9q8KOcZqEzNiEEDQ0NNDQ0cMMNNwBWOMl5fb300kscP34cgLq6Opej9jNjC4VCXHvttVx77bWAdU+Njo66Zi9nzpwBLJ+TPWMLIsxUKFd1h5CNEIKmpiaamprYt28fYMXjnRfz0aNHeeGFFwC0hpnC4TA7duxgx44dgDW6dDqtt956i1deeQVwOy2/YSY1I2ltbeXAgQMAGael6v7lL3+JlNLltHSEmWpra9m9eze7d+8GyIRXVL1nzpzJOC3njK2rq8tXmKmioiJzLLCc1uTkpMtR55qxqf+zjhmb02k5R5fOGVtdXZ22MFMoFGLLli1s2bIFsJzW2NiYy2nlmrH5DTOpxf3GxsbMjG1+fn7VGZvTUfudsW3fvp3t27cD1j01PDycqff8+fOcOnUKsAZKusJMZWVltLS00NLSkpmxzczMuPzI4cOHUc+iCSLMVCimQ8jDWsJM4XDYdTGn0+mC6/UTZorFYr7aXFdXx/XXX19QmEmN8AvBS5hJzdiyw0x+4vFewkwnTpzIzNiyw0x+HjKVHWa6fPkyw8PDnsJMyWSy4HrLysqIx+PE43HXjM1rmMlPm/2EmXKJRLyiZvmFhJlaW1sLrhdyh5mcg5/Vwkx+/MhaEaX0xLT9+/dLHbJTP7HxbFSYycuFWldXV3CYKReLi4vLFC25KC8v16pmUmEmL0SjUa1qptnZWebm5vKWq6qq0qpmunTpkkvBsxqNjY3a1EwqzLS4uJi3bE1NjVY1UzKZdClpVqO5uVlbbFyFmVKpVN6ykUhEq5opkUgwNTWVt5wKgeryIyrMtJo//vKXv1zQ9SyE6JFS7vdU9krvEJQ6YCXUSFyn05qbm2NwcDBvuVgspkXNpFBhpmzJaTahUEiLmsnJ9PS0p04iHo9rUTMpVJgp32g5HA5rUTMpVJjJi8Nsb2/X7rT6+/vzzsYikYhWp6XkptPT06uWE0JoUzMp1BpbPqLRqBY1k0KFmebn51ctp9YGdO4/mZmZYXh4OPPexz4Ezx2Cp5CREOJe4G+BcuBJKeVfZn1fBXwX2AdMAA9KKd+xv3sMOAikgT+RUj5nf/4d4IPAqJRyrxc7CsEpOz1w4IBrejgyMpKRyGXH8erq6rTJTg8dOuQKM/X19ZFIJJiYmCCRSGhTMymcstOtW7e62jw+Pk4ymeTChQva1ExOlOz04YcfdoWZ+vv7SaVSjIyMsLi46JoS65ixOWWn9fX1rnM9NTXFwsICFy5c0KZmcqJkpx/+8Idd51pJMAcGBrSqmRRKdvrII49krmPV5tnZWWZmZkgkEq4wU0dHhxanpWSnd955p6veoaEhpJT09fVpVTMplOz00UcfXbYGsri4yOTkJIuLi9rUTAo1sLzrrrvYu3fvMil3KpXiwoULWtVMiqKSnQohyoG/B34X6AdOCCGekVKecRQ7CExJKbcJIR4C/gp4UAixG3gI2AO0Af8qhNghpUwD/xt4AqsjWRfq6uqoq6tj716r/1laWnLF8V555ZXMTuja2tplcjE/EsyNUjOp2HixqZl6e3tzxsaNmmntqPUAo2a6stVM64GXbvNW4JyU8m0AIcQPgfsBZ4dwP/BV+/VPgCeEdVfdD/xQSrkEnBdCnLOPd1RK+UshRLeORhRKVVVVTrmY02kpCabaieh01H6cViFqJnVBX2lqJnUjGTXTladmUmEmo2bSr2YKAi8dQjvQ53jfDxxYqYyU8pIQYhqI2Z8fy/pte8HWBox
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Creating mesh and define function space\n",
    "Nx = 0.013\n",
    "Ny = 0.0015\n",
    "eleX = 10\n",
    "eleY = 5\n",
    "mesh = RectangleMesh(Point(0,0),Point(Nx,Ny),eleX, eleY,'left')\n",
    "V = FunctionSpace(mesh, 'P', 1)\n",
    "plot(mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Boundary conditions\n",
    "tol = 1E-14 # tolerance for coordinate comparisons\n",
    "#at y=1\n",
    "def Dirichlet_boundary1(x, on_boundary):\n",
    "    return on_boundary and abs(x[1] - Ny) < tol\n",
    "#at y=0\n",
    "def Dirichlet_boundary0(x, on_boundary):\n",
    "    return on_boundary and abs(x[1] - 0) < tol\n",
    "#at x=0\n",
    "def Dirichlet_boundarx0(x, on_boundary):\n",
    "    return on_boundary and abs(x[0] - 0) < tol\n",
    "#at x=20\n",
    "def Dirichlet_boundarx1(x, on_boundary):\n",
    "    return on_boundary and abs(x[0] - Nx) < tol\n",
    "\n",
    "#bc0 = DirichletBC(V, Constant(0), Dirichlet_boundary0)\n",
    "#bc1 = DirichletBC(V, Constant(0), Dirichlet_boundary1) \n",
    "bc2 = DirichletBC(V, Constant(0), Dirichlet_boundarx0)\n",
    "bc3 = DirichletBC(V, Constant(0.003), Dirichlet_boundarx1)\n",
    "bcs = [bc2,bc3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f8956953eb8>,\n",
       " <matplotlib.lines.Line2D at 0x7f8955f23278>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAABWCAYAAADPPOFDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsvXmQI9l93/nJO3HfBaAKVQCqu7q7+j5nOAcpahg8JYqiZZlcybosWbLWWm/YG7thaSN2ZXnDthwbYVsryRYlSpYoLQ8dvCSKx5Aczgw5w+lz+p6+gDq67gMoAIUbuX+8bFQfVQVMN4bkzOIbgYgq4GXi5UNm/vL93vf7/UmWZdFHH3300Ucf8ve7A3300UcfffxgoB8Q+uijjz76APoBoY8++uijDxv9gNBHH3300QfQDwh99NFHH33Y6AeEPvroo48+gH5A6KOPPvrow0Y/IPTRxxaQJCkrSVJZkqTiXa8nJUmy7nuvKEnSR+xt/rv9+Yfu29d/st//+e/LwfTRRxdQv98d6KOPH3B80LKsZ+/8I0lSyv7Tb1lWY4ttrgE/C3ze3kYF/hFw843rZh99PDr6M4Q++ug9vgg8LUlSwP7/fcB5YO7716U++uiMfkDoo4/eo4KYHXzU/v9ngT/7/nWnjz66Qz8g9NHH9vicJEk5+/W5u95fuuv9nCRJ4/dt92fAz0qS5Ad+CPgcffTxA47+GkIffWyPH99iDSG8zRoClmW9KElSBPjfgb+1LKssSdIb2tE++nhU9ANCH328cfhz4P8Afvj73ZE++ugG/YDQRx9vHH4HeAF4/vvdkT766Ab9NYQ++ng45O7TIfyr+xtYlrViWdbXrX7RkT7eJJD652offfTRRx/QnyH00UcfffRhox8Q+uijjz76APoBoY8++uijDxv9gNBHH3300QfQDwh99NFHH33YeFPpECRJ2pISJfHGq0AttmdkPdiHHvXprt1YVqtj460UsT3lk3Xqh/Q9eNbo2AcJkHr2M2yJVod+yN+DsejYB3ssHhZbX3obTQCr2eEaUb4H12mnPsgS3YvGH/6qaTW3/1xRHnrX2+LuQ2uIPixZlhXpZts3VUAwcLBDOUCutUTOWmSdAgASMl4piF8KEzCH8KtRDNmx6T4k0+j8RY7Ntz09/1lKtRX2ht/FamOe1fI0+fXbNK266J/qxu9M4HcmCDiGcQWGkDvcGFtG55+gaW60yecnOXP6v7Jj5wdQVZN8boJ8Pku5vCyOT1bw+BJ4A0l8gRTeQBJNd9n76Xxjahqdr5SGIfHaV/4b9XKB4eMfpLiQobiQpbQ0SatRA0B3+XDF0vYrhSMQR7rvxtjUO34VrW3aFCevk/3L/0r0Pf8ASZYpT2dZn85Qz9ljoSiY8RHMZArHSBrHcArF6dr8e7Qu+mJsfnO4/Tu/A60Wgfe+l/JkhsqtDLWJSay6cLZQA36M0RTGaApzNI2ajDwwFg9A7xT4QdY37jjrZ64y/39/gvAvfxir3qRybYLqtUkai6sASKqKOTaEY88Ijj0jmLsTKC5xnpt6veN3ObQtXTracOlVXvmlT6C4DJIfOUbu4gz5SzMUrs7Tqou+OmMeggfihA7ECB6M40kGkOR7zzm3Wt32e/xaedvPr391gq/8+rf58X9/jEqhztTZFabOLbM2K7ZTDZmhAwHGjnlIH/GTPOTD4dn8BPAp6x2PO6QWN33/V5+5wsCwzgf+cZjsmRwXTld47UKFurhEGBzROHDcZP8xkwPHTY6MNbZ8mLsDv9w5QAXkjXvcJz5d4Ff+5fJEx41svKkCgozCkLyDIXkHADWrSs5aJG+JADHZusbE+lUAnLIXvxoloETxq1Fcsq/jYHcDVdYJO1KEHLsBaFktipUFVstT5Nanya1PM7d2BQBlQsfnSuB3DxNwD+NzJVCVLgJSF3C5BgiFdjM4eAKAslRibXWCtdUs+dUJbmdeZPrWtwBwuCJ4A0k8A2k8oRQOd6Q3Y6Gb+BJ78CX2ANBqNSmvzFBcyFBYylKcucHqjTMAyLqJayCJK5bCHU3jjCZB781YGJEYruROAkefBKBeXKM8nWkHiNXvPMfqi98AQI9EcYykMe0AoQXDPRkL2eHAOT6OeUiMhdVoUJu+TeVWluqtLOWr1ymdFGMhOUyM0STGzhTmjiR6egTZ6CI6dgF9JIYxmsD77sfFd5VWKF+donx1kvLVSVY++6KYTUgS+nAEx54RvPsGce4dRo/6ezIWmtck/OQOwk+K67RVa1C4vkD+4gyly9PMvzzJ1JdfE209BsEDMREgDsQJjA/07K4U3eMjssPL8Y+OApCfW2faDg7TZ1d49g8yWGIoiO1ykz7iF6+jfgKDZk/GwhdUefzdPt7/XjEdqFZbXLtQ5cKpChdOV/j2syW+9JfiwdYfkDl8TOPIcZ0jx3X2H9QxzO+t/9WbKiDcD10yGJASDJAAoGU1KZjr5Brz5BrzLNanmKldB0CTDPzKAIHGEH49hk8bQJEe/fBlScbriOF1xEgG7ZtzPU9ufZrV2gy54hS3Zp9HTD0lPI4ofvew/RrBMEKP3AcA3fAQju0nHNsPQLNZp5ifJr+aZW11guX5y8xPnwJA1V14Qyk8oRTeUAq3P4Gs9GAsZAVXeBhXeJiw8UNYlkWtuEpp9hal+SzFuQxzp74KWCDJOMKDOOMpXPE0zsE0utv/yH0A0NxetD2H8O45BEBdqlGZmaQ8maUymaFw6VXyp18GQHG5cYykMVJpHMkUxmACSX30sZBUFSOVxEgl4RkxFo2lZaq3slQmMlRvZsl//ivkAWQZfXgQY2dKvHakUAfcj9wHADXoxfPkPjxP7gOgVa5SuXHbDhBTFF68QP6r9nkRcOMaT+AcH8a1N4FjNIakPnpeQ9ZVfPsG8e0bxKsfwLIsStN5li/MsXJhlpULc8x/RzzESqpMaHeIyMFo++UIbT5jf73wxZz43u9k3/vF/cJZzTP5ap7M2RyZszlOf3GW73xqWrSNGqSP+Bk/5mTnUS/De1wo6qPfnA1D5sBxBweOi2OyLIvJW3UunKpw/UyRs6dqPPesmCFpOuw7IALE0RM6h4/p+CNvbPrxTR0Q7ocsKQTUKAE1CojBXm/lWW0skGvaQaI4BdhpJi1CQIsT0GP4tRiG4uxJPxyaD4fPR8xxGIBGs0q+NM1qcZJccYqZ5XNMLZ4EwNC9+D0j+LxJ/J4R3K4osvToF6GiaPiCaXzBNCDWHor1ZdaWsxSWs6ytZFmZvQSAJKu4Awm8oRSuaBpPJIVmbJ5aeT2QJAnDE8TwBAnuOg5Ao1pmfT5LcT5LcSHDyuVXWD7/IgCaJ4ArnsIZH8UVT2GG4vSC9yDrOs7UTpypnQBYrRa1xXnKkxnKUyJIFK9cEH1WVYyEnWZKpjFHtk4zvR5IkoQWCaNFwrjecRSAZmmd2q0JKjcmqN7MUnz+ZQpfF2OhRoIYO5MYY0mMsRTa0EDnNFMXkB0GzgOjOA+Ip2ar2YK5GdavTFO6MkXp8hT574hZtqSrOHcN4hofJnBwEPf4EKrbfOQ+SJKEe9iPe9hP8gNiRlXLV9oBIndxhtf++gpXPnkRAHfCw8BdAcI3ZjyQZnoYmC6VXU+G2PWkeChrNS1mrxXbASJ7Jse5L88DYDhl0gc9jB31sPOol9HDHpyeHjw4SBLJHTrJHTqhnxJpq5XlJudO1zhzssbZUzX+/E9K/MkflAAYHVU4cULnxHGdEyd0duxQejKTuYO3VEC4H5Ik4VL8uBQ/CXYBUNdb5GpzrNbnWK3NMrF+nuz6OQCcio+AHsffGCZgDOHSgr1JrSgGIe8OQl57Cm21KJbnyRUnWV2fJleYZH5ZnPyKrOPzJPB5RvB7k/g8w0g9+JkkScbpjeL0RomlRTqhVilSWMmytpyhsJxl5voLWNeeA8D0RPBE0njCKTyRFKanR2kmw4F3ZBzvyDhNHaxmk/LSDOuzGUqzGYq3b5K7dhYQaSZnPIlzKIVzMI0jnkTpQZpJkmWMaBwjGsd/QqSZapU1KhMZyhPiCT73wnPkviXSTFokKmY
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#4 Defining variational problem and its solution\n",
    "rho = 1.43E-7 #ohm*m\n",
    "G = 1/rho\n",
    "voltage = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "f = Constant(0)\n",
    "a = dot(G*grad(voltage), grad(v))*dx\n",
    "L = f*v*dx\n",
    "\n",
    "# Compute solution\n",
    "voltage = Function(V)\n",
    "solve(a == L, voltage, bcs)\n",
    "\n",
    "# Plot solution and mesh\n",
    "plot(voltage,)\n",
    "plot(mesh, color='k', title=\"FEM\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#saving the solution to a VTK file\n",
    "vtkfile = File('solution/sol.pvd')\n",
    "vtkfile << voltage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.003 ,  0.003 ,  0.0027,  0.003 ,  0.0027,  0.0024,  0.003 ,\n",
       "        0.0027,  0.0024,  0.0021,  0.003 ,  0.0027,  0.0024,  0.0021,\n",
       "        0.0018,  0.003 ,  0.0027,  0.0024,  0.0021,  0.0018,  0.0015,\n",
       "        0.0027,  0.0024,  0.0021,  0.0018,  0.0015,  0.0012,  0.0024,\n",
       "        0.0021,  0.0018,  0.0015,  0.0012,  0.0009,  0.0021,  0.0018,\n",
       "        0.0015,  0.0012,  0.0009,  0.0006,  0.0018,  0.0015,  0.0012,\n",
       "        0.0009,  0.0006,  0.0003,  0.0015,  0.0012,  0.0009,  0.0006,\n",
       "        0.0003,  0.    ,  0.0012,  0.0009,  0.0006,  0.0003,  0.    ,\n",
       "        0.0009,  0.0006,  0.0003,  0.    ,  0.0006,  0.0003,  0.    ,\n",
       "        0.0003,  0.    ,  0.    ])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "volt_nodal = voltage.vector()\n",
    "volt_array = volt_nodal.get_local()\n",
    "volt_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "volt_array(       0,       0) = 0.003\n",
      "volt_array(  0.0013,       0) = 0.003\n",
      "volt_array(  0.0026,       0) = 0.0027\n",
      "volt_array(  0.0039,       0) = 0.003\n",
      "volt_array(  0.0052,       0) = 0.0027\n",
      "volt_array(  0.0065,       0) = 0.0024\n",
      "volt_array(  0.0078,       0) = 0.003\n",
      "volt_array(  0.0091,       0) = 0.0027\n",
      "volt_array(  0.0104,       0) = 0.0024\n",
      "volt_array(  0.0117,       0) = 0.0021\n",
      "volt_array(   0.013,       0) = 0.003\n",
      "volt_array(       0,  0.0003) = 0.0027\n",
      "volt_array(  0.0013,  0.0003) = 0.0024\n",
      "volt_array(  0.0026,  0.0003) = 0.0021\n",
      "volt_array(  0.0039,  0.0003) = 0.0018\n",
      "volt_array(  0.0052,  0.0003) = 0.003\n",
      "volt_array(  0.0065,  0.0003) = 0.0027\n",
      "volt_array(  0.0078,  0.0003) = 0.0024\n",
      "volt_array(  0.0091,  0.0003) = 0.0021\n",
      "volt_array(  0.0104,  0.0003) = 0.0018\n",
      "volt_array(  0.0117,  0.0003) = 0.0015\n",
      "volt_array(   0.013,  0.0003) = 0.0027\n",
      "volt_array(       0,  0.0006) = 0.0024\n",
      "volt_array(  0.0013,  0.0006) = 0.0021\n",
      "volt_array(  0.0026,  0.0006) = 0.0018\n",
      "volt_array(  0.0039,  0.0006) = 0.0015\n",
      "volt_array(  0.0052,  0.0006) = 0.0012\n",
      "volt_array(  0.0065,  0.0006) = 0.0024\n",
      "volt_array(  0.0078,  0.0006) = 0.0021\n",
      "volt_array(  0.0091,  0.0006) = 0.0018\n",
      "volt_array(  0.0104,  0.0006) = 0.0015\n",
      "volt_array(  0.0117,  0.0006) = 0.0012\n",
      "volt_array(   0.013,  0.0006) = 0.0009\n",
      "volt_array(       0,  0.0009) = 0.0021\n",
      "volt_array(  0.0013,  0.0009) = 0.0018\n",
      "volt_array(  0.0026,  0.0009) = 0.0015\n",
      "volt_array(  0.0039,  0.0009) = 0.0012\n",
      "volt_array(  0.0052,  0.0009) = 0.0009\n",
      "volt_array(  0.0065,  0.0009) = 0.0006\n",
      "volt_array(  0.0078,  0.0009) = 0.0018\n",
      "volt_array(  0.0091,  0.0009) = 0.0015\n",
      "volt_array(  0.0104,  0.0009) = 0.0012\n",
      "volt_array(  0.0117,  0.0009) = 0.0009\n",
      "volt_array(   0.013,  0.0009) = 0.0006\n",
      "volt_array(       0,  0.0012) = 0.0003\n",
      "volt_array(  0.0013,  0.0012) = 0.0015\n",
      "volt_array(  0.0026,  0.0012) = 0.0012\n",
      "volt_array(  0.0039,  0.0012) = 0.0009\n",
      "volt_array(  0.0052,  0.0012) = 0.0006\n",
      "volt_array(  0.0065,  0.0012) = 0.0003\n",
      "volt_array(  0.0078,  0.0012) = 0\n",
      "volt_array(  0.0091,  0.0012) = 0.0012\n",
      "volt_array(  0.0104,  0.0012) = 0.0009\n",
      "volt_array(  0.0117,  0.0012) = 0.0006\n",
      "volt_array(   0.013,  0.0012) = 0.0003\n",
      "volt_array(       0,  0.0015) = 0\n",
      "volt_array(  0.0013,  0.0015) = 0.0009\n",
      "volt_array(  0.0026,  0.0015) = 0.0006\n",
      "volt_array(  0.0039,  0.0015) = 0.0003\n",
      "volt_array(  0.0052,  0.0015) = 0\n",
      "volt_array(  0.0065,  0.0015) = 0.0006\n",
      "volt_array(  0.0078,  0.0015) = 0.0003\n",
      "volt_array(  0.0091,  0.0015) = 0\n",
      "volt_array(  0.0104,  0.0015) = 0.0003\n",
      "volt_array(  0.0117,  0.0015) = 0\n",
      "volt_array(   0.013,  0.0015) = 0\n"
     ]
    }
   ],
   "source": [
    "coor = mesh.coordinates()\n",
    "if mesh.num_vertices() == len(volt_array):\n",
    "    for i in range(mesh.num_vertices()):\n",
    "        print('volt_array(%8g,%8g) = %g' % (coor[i][0], coor[i][1], volt_array[i]))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the projection of the voltage function; the soltion. This is:\n",
    "$$\\nabla u$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.tri.tricontour.TriContourSet at 0x7fc67e808278>"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAABWCAYAAADPPOFDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHCtJREFUeJztnX20JEd5n5/fzOzM3bva1bKrHNv6QBIW5iDZCQiM4wTHGGwj4YD8h4lF4g8ZCZkYbBwSHBQSGzgQ44ScCBsIAYHxFxICbM4Gg+UQQTAQSSiQYK2EYn2srUUCa3e12tVd7czOnTd/VPfcmr7d090z3fNxt55z7rk93dU1Vd019ev3fauqZWYEAoFAINCYdwECgUAgsBgEQQgEAoEAEAQhEAgEAhFBEAKBQCAABEEIBAKBQEQQhEAgEAgAQRACgYVDjt+V9JikO+ZUBpN0kff5xZI+WfDcl0r6aH2lC9RFEIRAoGIkXSXpi1Nk8Xzgx4Bzzex5FRVrWt4OvKNIQjP7b8Alkv5uvUUKVE0QhEBg8TgfOGBma9NmFFkbU/3OJX0/cKaZ3VbitBuBa6f53sDsCYIQKISk75Z0RNKl0eezJT0q6QVjznmVpHskHZd0t3fuMyV9XtJRSfslvcw758OS3ivpM5KekPQlSd8p6frIhfINSc/20h+QdF2U/2ORq2UlUYb7orLvk3S2d8wkvVrSX0VleY8kecdfGZX/MUm3SDo/71xJzwTeB/xgVP6jGdfm7Kg8R6LyvSrafzVwg3f+W1LObUr6T5IOSXpQ0muj8rSi45+X9HZJXwJOAE+T9AvevXhA0i8m8nyDpEckPSzplYmvvBz4n17aC/zv877zGu+czwM/kVb3wAJjZuEv/BX6A14F3A2sArcA7xyT9uXAN4HvBwRchHvy3QbcB/wboA28EDgOPCM678PAIeA5wApwK/Ag8HNAE3gb8Dnvew4AdwHnAXuALwFvi469MMrrUqAD/A7wBe9cAz4F7AaeCjwKXBYduyIq5zOBFvBvgS8XPPcq4Is51/ILwHujOj4rOv+FRc4HXh3dh3OBpwCfjcrTio5/Hvgb4JKo7NtwnfN3R/fih3FCcWmU/jLg28D3AjuAj0T5XRQd/xjwBu/7L/C/z/vOa7zPe6I0u+bdbsNfid/4vAsQ/pbrD9gH/CXwdaAzJt0twOtS9v8Q8C2g4e27EXhztP1h4APesV8G7vE+fx9w1Pt8AHi19/klwP3R9geB/+AdOwM4BVwQfTbg+d7xm4E3RtufAa72jjWiTvT8AufmdejnAevATm/fbwIfLnj+rcAvep9/NEUQ3ppzHz8Z3x/gQ8A7vGPfkxCE/564xkUEYVuU5qnzbrPhr/hfcBkFyvIB3JPk75hZF0DSD0XujSck7Y/SnQfcn3L+2cBDZjbw9v01cI73+dve9pMpn89I5PlQIq/YLXR29BkAM3sCOJz4rm952ye8vM8H3hW5g44CR3BP10XOzeNs4IiZHU+U+5yM9Gnn+3V+KCXNyD5Jl0u6LXJRHcUJ51kZ+f01ozwG7CxYtpg4farLLLCYBEEIFEbSGcD1uCfvN0vaA2Bmf2FmZ0R/l0TJH8K5KJI8DJyXCHQ+FedempTzEnk97H2X7/ffAewt+F0P4Z7Cd3t/283sywXOzVtC+GFgjyS/ky1zDR7BuYtizktJMyyDpA7wCeCdwHeY2W7g0ziBi/NLXkOfr+Oshpg42L3q7fvOxDnPxAXGj2VXI7BoBEEIlOFdwJ1mdg3wp7jgaRY3AP9K0nOiYOtFUVD2dtzT9K9J2hYFpV8K3DRFuV4j6dxIoN4ExGPgbwR+QdKzok7x3wO3m9mBAnm+D7hO0iUAks6U9PKC5fk2cK6kdtpBM3sI+DLwm5JW5IZnXg38YcH8bwZeJ+kcSbuBf52Tvo2LoTwK9CVdDvx4Ir+rJF0saRX4jcT5n8bFHeLyP4oTr5+JAtyvZLP4/zDO7RZYIoIgBAoh6Qpc8PGfR7teD1wq6Z+lpTezj+HGrn8EFzT+JLDHzHo4AbgcF/B9L/BzZvaNKYr3EeDPgQdwbqq3RWX4LPDvcE/Hj+A6rSuLZGhmfwL8FnCTpGO4wPXlBctzK7Af+JakQxlpXoHzxT8M/AnwG1F5i/ABXH2/DnwN12H3cXGJTUSuqV/BdfyPAf8UFwuKj38GZ/ndiguk35o4/6vA45J+wNv9KuANOBfcJTiBS9bvvxasT2BBkFl4QU5geZF0ABfMLNqZbjmiJ/73mdn5uYkn/44fB37JzH6yQNqXAj9rZv+krvIE6iEIQmCpOR0FQdJ24EdwVsJ34Cyg28zsV+dasMDSE1xGgcDyIeAtOPfP14B7gF+fa4kCW4JgIQQCgUAACBZCIBAIBCKCIAQCgUAAcOucLA3bOjuss7pn3sVAZdxsg/wkE+edg22s05ZPwUeDrDwLlbvAtcjLR4Pi18caJeofn1Pmmvkkrl9mPinX2VKSWtb9SEvr72t4++TlFaeRQQMkdx0loxFtNxpGU+4mNeVvD4bbLdZpeftb0U1tCZophVvH6Ee3rB8Vbt0a9KMK9mmybhv7N7Y13DcYuO2BCYv2m2mjPZmG0/A03OeqOiSj7SmtOWU0MaXkkXp+yd98KSZ4hF87cvCQmf2dImmXShA6q3v4ey963byLQatbrtNudMu1kObJfqn041hfKX6LB51ira3f2fzDL3JNilyHInVvnjiVm8ZnfXVbqfRQ7rrB5muXdo1cvpuvcb+Tki7r/LS03vS3QWd036BjDDru3gy2GdZx96DRXqfVcdMW2u0+q+0eADvaPXa2uwDsap/kzG1Psqt1EoA9rTXOarnVNnY319jbdBOW9zRPsLfh8n1KY6OAjw26w+3DA1fvI+tucvPh9R0cXd8BwKH+To703faxvluo9vFT2znWc9vHex3Weq5CJ3ptej13b/rdpqtXr4m6Lv/GKdHoumsX/2+6qtGIihN/9ml2N+9z+ze361ZW2pPVKkFaW5mE2//oXyaXIslkqQRhWR1cg06jtChURfNkv3Dn1ugOColCq2sjHV5RgSxyHdZXWpUKIowKSFFx8MtQ5Polr13yGrl8iolBFnlikF02MegYjVNiQAPrDBj0mvSBVmd92MGutnvDjjcWhcdPbQdgV+vksNM+q3V82Jnvba5FnbwTBV8EYEMIYEMMYnY31zi6vmMoMjHH+iucue1Jt91bGZZlrdceCldMv9uk0V5nAKjbYLDNRuvcFettJwKDjhOF+JrFwpAUgzQRgHQhyBKB5G8i6wHBpyoBmIblEgTcRataiReNOjrFopQRhXmxvrqttJUQE59XxmqI70WeMOSJQvPkoNCPPss6mJa4s/Tpd5tDS+GE1+Ee7zn12dU+WVoUIFsIDkfnAMPzk+xqnRwRhTRO9Nq02/1hHdJEYZTRaxoLQ7PnhDYWhaQYTCMCWceKiMO8WDpBgPmLQr+j2jvEqkShrOsDiotC2TyLMk9BHEdRYRifx6gotLrFrIQ066AosctoHL1ei3a7z4nIQtjR7pUWBRgVAigmBof6mxdS9UXh8VPb2dV2bqvjvQ47ItHyRQHc2h1xKxvQcBZRZCVsXIO8zjhyMyV+32VEIKutx7+pLHFIfsc8LIalFASYvyiUZRK30bQd4zQdV5WiUIe7bFIrYZJ4QpJx7qQirqOqyHMXFRGCLNZ67dKiAC6eAOWFIM4PNuIIsOGyitnZ7o6IQkyv16LVWXei0Gs6t1jCx5znRvKthQ0rLb6GLq+4z0kKQZE27qfJEwf/u2B24rC0ggDLJwqzYhohqJpJxKBq66AKEcgizWoo6zpKWgnNrk3lNioSVyhK7L/3OdZfGQaaD/V3DmMAhxMuoKqEIA4uw4YrK0m73R8RBSDHhQR51kIsDM2u0e+4+xT3Ob6XINnGs9pvso3EpIkDZFsPZcShTIwKllwQYOuLQlkroSoxqMI6qFsM8qyEOoUgSTJ4X7UouM5pfBkGU7iVYvzRRnkc6e9gT2tt2MHvbq5tEoWqhABGxSAOfvuMiIIXV4ANYUi6kcYFnTd
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "V_g = VectorFunctionSpace(mesh, \"Lagrange\", 1)\n",
    "w = TrialFunction(V_g)\n",
    "v = TestFunction(V_g)\n",
    "a = inner(w, v)*dx\n",
    "L = inner(grad(voltage), v)*dx\n",
    "grad_u = Function(V_g)\n",
    "solve(a == L, grad_u)\n",
    "plot(grad_u, title=\"grad(u)\")\n",
    "grad_u_x, grad_u_y = grad_u.split(deepcopy=True) # extract components\n",
    "plot(grad_u_y, title=\"y-component of grad(u)\")\n",
    "plot(grad_u_x, title=\"x-component of grad(u)\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous code can also be substituted by the method project, like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Calling FFC just-in-time (JIT) compiler, this may take some time.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.quiver.Quiver at 0x7fc67efcbcc0>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAABJCAYAAAA9vFENAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXd8HNW5979nd7V91XvfpmpbXbuqljEGQkIg3FwguTfvzQ2ENz33pgNJwCSkXW4SSAcDoYUSejHVNs2AG2Djgptsy92yZXVpV9LO+8cuRitmzWq0Wtsv8/18/LH2nJk5Z3ZnTnme33mOkCQJFRUVFRUVzamugIqKiorK6YHaIaioqKioAGqHoKKioqISQu0QVFRUVFQAtUNQUVFRUQmhdggqKioqKoDaIaioqKiohFA7BBUVFRUVQO0QVFRUVFRC6E51BaZDenq6VFxcPO3zJgISw6N++UwJDh7tIxAIX7GtT9BiNRuwmg2YjHqEgvoCDI/4P3Tt9xkYGqVvYDQsTaMRWEx6rJZg2Vqtsj57fHyC0dEx2TwpIHHoYB9TV6kbjAlYrQYsViNGY4KicgGGhn1EWgB/vH+YoZHw30KrEVhC37XVbECjUfZtj/nH8fnGZfPGxyc4crg/LE0ARlMCFqsRq9WI3qD8dRgc8kXMO3psEJ8v/LfQ6bRYLAasFgMWsx4hlN3z6Ng4Y+MTsnkjY2N09w+FpQkhsBgSsBkNWI0G9DqtonIhANJg5OyJAyBN+S2EHoQNNDYQFlD4Vg37xxgPBGTzBn0+eoZGwtI0QmA16LEZDdiMBnQahe/URICRCO8UksSB7v4Pv1N63aR2ZAbv1IifQISXqm9ghIHh8OdPqxFYTAZsZgPbt246KklSRjTliDMpdEV9fb20du3aaZ+3eechvvTTexWXm2gx4q0qprXGQWutE7NRH/W5X736H7y79YCicoWACncOLfVO5nvcFOWnRX3uy8s28/NrH1VULkBqmgVPsxtvs4sGr4uEhOgbjk9f+VeOHR/66ANl0Go1VJfn01LvYH6jm+yMxKjPffDulSz50zJF5QJk5ybjaSmhqa2E6np71B3T2NgEZ1/0v4rL1et11FUX0dzopK3JTUqyJepzr390GQ+8uUFx2c7MVDoqHJxV4aSqMCfqjkka70Q6ep7ichFWMLQhDAvAsBChsUV96pfveZRXd+xWXPSc3CwWlDg4u9xFaVZ61Oe9+c4uvvOLRxSXm5xoornGQUudg+ZaBwZ99AOQS35wB3sOHldU7up7vrtOkqT6aI5VTUZRkJeVRFFOKkW5qZgMynv56ZJkM1GUl0phbioZadG/MDNFoxXkFaRRUJBGQVE6Ol38HpP0ZAuFuSkU5qaSkmSOW7kJCVryC9MoKEojryBV8SxFCVkZNgrzUinMT8VmM8WtXJM+AXtmKvb0VPJTkxTPUhShLQCtA7T20GwhPiQZDdjTUrCnp5CTaI1buUJAQU4KhXnBZ1s/jQFWPPlYzBBG/WMc7O6XzRsbm+DrNzwYNuUy6HU0zimitdZBc7WDjBTlD87BI334/PJmjIeefovHnl8fluYoTKel3klLvZNyV7Zik9HQkI+j3QOyeX29w3z/2/cQmPjgt7dYDTR4nHhb3DR4nSQmKm+Y9h48zsSE/JT+z/e8wsp1nWFpFa7s4D3XOXAVZShumPr7huntkZ+Z7O7s5ufXPBSWlpxiwdPqxtNSQl2jA5M5+pnfZCRJYs/enoj5N9z4FNt2Hj7xWasRzK3Mp9njornRSUFeqqJyAbr7B+kfkTdXrdy2h18/9XJYWk6yjY5yBx3lDhqd+eh1ysxkkuSHia4ImQGk41+CwJFJiXowNIVmBAsQ2hxF5QIc6O1nZEzedPPw25u4/fV1YWmO9FQWlNhZUOqgOj8XncJ3anjUz5Gj8u/U0Iifr1/3IP6xD951szGBxqpiWuucNNXaSUlUPsDZd7g3omnwtsfe5IU3t4allRZn0lbjoLXGSYUjO+oZwhnlQ1CKUZ+APU/e3PLIi+sZGPaRkWI9YRKqqyzAqI/NTCAnM0k2fXjEz4o3tpGg01JTWUBLvYPmemfE46eLxWLAYjHI5i35y3ICExK5+Sk0tbjxtriZM68AnWJ7cjgFOSmy6Ye6+3nznd0YDToa5xXTUu+gqcZBWkpsRoiJSWYSI8wq7l4SbBgdriw8rW68rSWUVuTFZCYghKC4UP752rh5P9t2HsZmNeKpt9PU4MRT58BmM864XICMRCsZMiNdSZK4+sHnEALmFeSc6ATc2WkxmQkIoQedSzZPGnky2Blo0sHQgTCcBfpmhCY2M77cZHkzon98nKff3YpOo6G+KI8FJQ46ShwUpSXHpFyzUU9xBLPtnY+swj82Tk5GIi11TlrrHVSX56NPiE0Tm58lfw89fcO88tZODAlaGiqLaK1x0FJtJzNVmUXhY9EhnIy1W7r4xufbuWjBPKxm+QZ0NnjpjW14a+386yfrKHNmxa1cv3+czl1H+PYPz+ecc+bOyJE6XZ57eTMXnjWHSz5VR0Gu8lHxdOk5NsjgkI+rb/gX2uaXo42jCWzFS1v4t880cPFF9aSnx8/st2n/EdITLfzlnM/Q4i6Kqwns0PDrCNP3SbNcQoIuNgOcaFi2awc1Fdl8oaaW2qzcuJnAAoEAW9/dx/c/P5/zzp2HyaRspqmE51Zt4dyWMi47txZnXvT+kEh8LExGJ+OFN7fy4z89TYJOS215Pm01DlqqHeRmzO6D7B8b57Jv3sbhowMU5KTQXOegpd5JVVlezEbqkbjv3tdZcutLGI0J1NXbaWp24/E4SU2bXZvq8b5hLv2/tzAyOoazOIOWBifNDU7KXTmz3mD9/jdP8/Tjb2OzGWnwOvG0uGnwOLHNwDQWDbs6j/Dly5cAUFaeR1Ozi6ZmN3a7ctNYtHz5jkdYuWMP6VYz80sddJQ5aHIVYo7R7DcSB4dW8tqhbyHQkmGqJcfcRq6lHWtCwayWOxEIcM7Tt9LZ30OBNZmz8pwszHPjySxEr53dd+qJh9bwxxufQW/QUVNvx9NagrfFTXpm9KIIJQyO+Ljwh7fRNzSKIzeNtioHbVUO5jpz0IbUVEKIqE1GH4sOYWBolA3b5ZU+gYDE9bc8S/8U2aAzP43WGidtNQ4qnNknvtzpsnHrAfqHRmXzlq18j2df3hyWZrMY8FTbaal34q0uJlGhk/HYsUF2bD8kmzc87OeXNzzxITt/WXkuTU3BBsvhzFTcYK1dv4exMXl750NPr2P127vD0lKSzDTVO2hpcFJfVYxZ4QjrwL4e9nbJ2/KPHO7j5hufCUvTaAVz5xXibXXjbSkhv0DZrCUQkFizemfE/NuXvMyOHYfD0rKykvCGvuuq6kL001CcTGb74aMc6JW3a7+77xB/Xv5mWJpep8XrKKSjzE5HmYPsJGWzlvHAMN0jb0XIlVjb/XNGJ46GpdoSism1tJNjbifNOBeNUHbP648doGd0WDZv2f4d3Lv97bA0q05PW66DhXlOOnJdpBmVma56jw+xbYt8O+L3T/Crax/BP0X27CrNCT5fbSW4SpQPfNZt3cuIT95v8virG1nx1o6wtCSrkZa5dtqrnCxqLFU7hMls7jzEf177D8XlpthMNFc7aK2x0zTXPi098VeuUS471WoEc8vyaKlz0tboimibl+OlFZv52eLHFJULkJmZiLfZhbfJTW1t8bRkpxf9558Vy071CVpq5hYGncyNLjKnYWZ54N7XWfLn5YrKBSgoTMMT8qnMrSqcluz0vEW/VlyuyaSnvsGOt8lNU5OLpOToG6zFjy/jgdXKZadlORksKHOwoMxBZV5W1IOAAf9unt37L4rL1WuSyDY3k2tpJ9vcQoImej/SF5c/wMsHOz/6QBkEUJuex8J8N4vy3biSojezrHljB9f8t/J2JC3DRmOzG2+rm7pG57TMtf9yzR3sOaRMdrrudlV2GlPSki1kpljJTLFNSzs8U8xmAxmpNjLSrCTFUY4oBKSl20hPt5GRYYur7NRqMZCeaiU9zYrNGj+fjlarIS0jeL/
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "grad_u = project(grad(voltage), VectorFunctionSpace(mesh, \"Lagrange\", 1))\n",
    "plot(grad_u,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heat Transfer Model\n",
    "The electrical heat generation is coupled to the stationary heat equation as the source term $Q$ on the right hand side, which is known as the electromagnetic heat source.\n",
    "\n",
    "$$\\nabla \\cdot(-k\\nabla T)=0$$\n",
    "\n",
    "Joule's first law states that the heat generated within an object is equivalent to the product of its conductivity $G$ and the square of the potential gradient:\n",
    "\n",
    "$$Q=G|\\nabla V|^2$$\n",
    "\n",
    "\n",
    "**Cambiar el valor de Q por uno numerico obtenido de las mediciones y también puedo tratar de calcular el valor de Q del Vector**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.quiver.Quiver at 0x7fc67e34d1d0>"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAABJCAYAAAA9vFENAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXmUG9Wd7z9Xu1pS7/uurVe3e29J3W27bfaEAMkkgeQN5/EGkskCk8kkk0AIwTBAZrINCWTgAYYQdsJqFrPZ7BiMjTE2Nl7a+97tdu+LpFa9PyScViMZdbUk24/6nNPntO6tqvsrqapu3d/ve39XSJKEgoKCgoKC6kQboKCgoKBwcqB0CAoKCgoKgNIhKCgoKCiEUDoEBQUFBQVA6RAUFBQUFEIoHYKCgoKCAqB0CAoKCgoKIZQOQUFBQUEBUDoEBQUFBYUQmhNtwEzIzs6WysvLZ7zfZEBidNwbuVKCA70DBALhM7Z1WjXmFD3mFD1Ggw4hw16A0THvZ479KUMj4wwMjYeVqVQCk1GH2RRsW62W12f7/ZOMj/si1kkBiYMHBpg+S11v0GI26zGZDRgMWlntAoyMThBtAvzRwVFGRsN/C7VaYAp91+YUPSqVvG/b5/UzMeGPWOf3T3L40GBYmQAMRi0mswGz2YBOL/92GB6ZiFrXe2SYiYnw30KjUWMy6TGb9JhSdAgh75zHfX58/smIdWM+Hz1DI2FlAoFJr8Vi0GM26NFp1LLahQBIw9GrJ/eDNO23EDoQFlBZQJhA5l016vXhDwQi1g1PTNA3MhZWphICs16HxaDHYtCjUcm8pyYD0e8pSWJ/z+Bn7ymt5ti9bJzNPTXmJRDlphoYGmNoNPz6U6sEJqMeS4qerZs/7pUkKSeWdsSplLqipaVFWr169Yz329h9kH/61QOy2001GXDXl9PZaKOzyU6KQRfzvt+/6kHWb94vq10hoMZZQEeLnQVuJ2XFWTHv+/ryjdxw7ZOy2gXIzDLhanfibnfQ6nag1cb+4DjvO7dx5OjI528YAbVaRUN1cfCcXU7yc1Jj3vfR+9/hrj8vl9UuQH5hOq4OJ57OChparDF3TD7fJKdf8HvZ7ep0GpobymhvszPP4yQj3RTzvtc/tZxH3vtIdtv23Ey6qm0sqrZTX1oQc8ck+bcj9Z4tu12EGfTzEPqFoD8NobLEvOt37n+SN7ftlN30nMI8FlbYOL3aQWVedsz7vbt2Bz+56QnZ7aanGvE02uhosdHeZEOvi/0F5Js/v4ddB47KanfVfT9ZI0lSSyzbKi6jGCjKS6OsIJOywkyMevm9/ExJsxgpK8qktCiTnKzYb5jZolILikqyKCnJoqQsG40meZdJdoaJ0qIMSgszyEhLSVq7Wq2a4tLg+RaVZMoepcghL8dCaVEmpcWZWCzGpLVr1Gmx5mRizc6kODNN9ihFFuoSUNtAbQ2NFpJDmkGPNSsDa3YGBanmpLUrBJQUZFBWlEFZYSa6GbxgJZMvxAhh3OvjQM9gxDqfb5If3vho2JBLr9PQNqeMziYb7Q02cjLkXzgHDg8w4Y3sxnjsuQ946sV1YWW20mw6Wux0tNqpduTLdhmNjEzQ2zMUsW6gf5R//9H9BCb//tubzHpaXXbcHU5a3XZSU+U/mPYcOMrkZOQh/f/c/wZvr+4OK6tx5AfPucWOoyxH9oNpcGCU/igjk53dPdzwy8fDytIzTLg6nLg6nDS32TCmxD7ym4okSeza0xe1/sbfPcuW7kPHPqtVgrraYtpdDtrb7JQUZcpqF6BncJjB8cjuqre37OK/nns9rKwg3UJXlY2uahtttmJ0GnluMknywuTuKJUBpKP/BIHDUwp1oPeERgQLEeoCWe0C7O8fZMwX2XXz+NqPufudNWFltuxMFlZYWVhpo6G4EI3Me2p0zMvhI5HvqZExL5df+yhe39/vdaNBi6u+nI4WO55G66xecPYe6o/qGlzy9Lu8/O7msLLK8lzmNdjobLRTY8uPeYRwSsUQ5GLQabEWRXa3PPHKOoZGJ8jJMB9zCTXXlmDQxWckUJCbFrF8dMzLq+9sQatR0zinJDiMbLFH3X6mmEx6TCZ9xLq7bltBYFKisDgDT4cTd4eTOXNL0Mj2J4dTUpARsfxgzwDvrt2BQa+hrb6cjmY7niYbWRnxeUNMTUshNcpNd99dbwBgc+Ti6qzA3VlBZXVhXEYCQgjKSyNfXxs27mNL9yEsZgOuFiueVjuuZhsWi2HW7QLkpJrJifCmK0kSv/jbiwgBc4sL6KoOdgLOvKy4jASE0IHGEbFOGnsm2BmoskHfhdAvAl07QhWfEV9hemQ3otfv57n1m9GoVLSUFbGwwkZXhY2yrPS4tJti1FEexW177xPv4fX5KchJDb7cNNtoqClGp43PI7Y4L/I59A2O8sYH3ei1alpry+hstNFRbyU3U55H4QvRIRyP1Zt2c/m353PBwrmYUyI/QBPBayu34G628o0vN1Nlz0tau16vn+07evjRz7/EmWfVoZuBH3O2vPj6Js5fVMc3v9JMSaH8t+KZ0ndkmOGRCX5xw9eYt6AadRJdYK++volvf62Nfzi/mezs5Ln9Pt5/mKxUE7ed8VU6HGVJdYEdGF2JMPyMLPM30WlijwHNluU7ttFYk8/FjU005RUmzQUWCATYvH4vP/3HLs4+sw6jUd5IUw4vvruJs9qruOisJuxFscdDovGFcBkdj5ff28zVtz2HVqOmqao4OMyqt1GYE5839Wh4fX4uumIJh3qHKCnIoL3ZRkeLnfqqori9qUfjwYdWctddr2MwaGluLqfd48DlspOZmVif6tGBUS785zsYG/dhL8+ho9VOe6udakdBwh9YN//meZ5buhaLxUCr246r3Umry4ZlFq6xWNixo4fLvrMEgOqqQjweBx6PA6tVvmssVi679wne3raLbHMKCyptLKy04bGXkhKn0W809o+8w2sHfoxATa6xkSJTJ0UpnVh0JQltdzIQ4Mxn72T7YB8l5nQWFdk5rdiJK7cUnTqx99TSJ1Zzy+9fQKfT0NRqxdXuwN3hJHsGogg5DI9NcP7PlzAwMo6tMIt59Tbm1duosxegDqmphBAxu4y+EB3C0Mg4H22LrPQJBCSuu+sFBqfJBu3FWcxrsNPZYKPWln/sy50pGzbvZ3BkPGLd8rc/4YXXN4aVWUx6XA1WOlrsuBvKSZUZZDxyZJit2w5FrBsb9XLTr5/5jJ+/qqoAj8dBu8eBzZYr+4G1et0ufL7I/s7HnlvDqrU7w8oy0lLwtNjoaLXTUl9Oisw3rP37jrJn95GIdYcPDvCn378QVqZSC+rmluIOuc2KS+SNWgIBiVXvb49af/fdr7Nt2+Gwsry8VNzu4HddX18qe6S25VAvBwYi+7XX7z3In199N6xMp1HjtpWysNJKV6WN/DR5oxZfYJTDY2uj1EqsOvxrxiZ7w0pTteXBzsHUSbahDpWQd87revfTNzEasW753m08sDXcLrNWx7wCG6cV2ekqcpBlkOe66j86wuZNByLW+Xx+fr34KbzT4oXOyvzQ9VWBoyJf9ovPms17GJuIHDd5+s0NvPrBtrCyNLOBjjor8+vtnNFWqXQIU9m4/SCXXP+g7HYzLEY66m101lvxzLXOSGn0vavly07VKkFdVREdzXbmtTmi+uYj8dprm7j+P56W1S5Abm4qHrcDt9tOU1P5jGSnF/yf/5EtO9Vp1TTWlQb9sG0OcmfgZnnkgZXcddsKWe0ClJRmHXuzq6svnZHs9Kyzfyu7XaNRR0tLOR53cPSQNoPg4+Kly3nkffmy0+r8HLqqbCyqslFbmBfzS8CgdxfP7v6m7HZ1qlQKU9opMnVSaGpHq4o9jnTJikd4fX/0Dvh4CKApp4jTipycUeLEkRa7m+X9d7v5xU8ektUuQFa25dj11dxqm9G8l3+4+h52HZQnO11ztyI7jStZ6SZyMszkZlrQxylIFAspKXpyMi3kZJlJS6IcUQjIzjKTnW0mJ8eSVNmp2aQnO9NMdpYZizl5MR21WkVWjoWc3FSycyxJ9bmnpRn
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "Q = 1/rho*grad_u\n",
    "plot(Q,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'assamble' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-47-080bf2ec9c68>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0menergy\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mQ\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mE\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0massamble\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0menergy\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'assamble' is not defined"
     ]
    }
   ],
   "source": [
    "energy = (Q)\n",
    "E = assamble(energy)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}