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": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fenics import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f2daacc7f98>,\n",
       " <matplotlib.lines.Line2D at 0x7f2daacb1128>]"
      ]
     },
     "execution_count": 27,
     "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": 28,
   "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": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.tri.tricontour.TriContourSet at 0x7f2daaa6a390>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAABJCAYAAAA9vFENAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAACO1JREFUeJzt3WuIXGcdx/Hvb++bxCZx9UXbBLulobB5Y2uo9YJUKySpmviiQgJCCylBTEXtC0kpiAYLFoWGklYJbTVWzcVYJNRqFVtQQZPGVmsujW4u2o2FXo0mNdf9++I8KZPJZs+zyZzZneH3gYUz5/yf5zz/PWfnP2eembOKCMzMzDomewBmZjY1uCCYmRnggmBmZokLgpmZAS4IZmaWuCCYmRnggmBmZokLgpmZAS4IZmaWdE32ACaiR73RpxnjxqijpMZJ5Tsq6wOgI6ef8pjIGk9OP+XdRNaYM/rJ2VdGTEP2lZV3xlgalVNW3hl3B8jaV3k/yjqVy/vJixktjelswL46M/bTlTWWM+UxNGZfWTEZ++rMOC86606e5184+VpEvLu8ZYsVhD7N4MauhePGdEyfNu529fWW76i/vzQk+sv7if7u0pjR3vJDcKYvJ6b8Lz8rprf8jDudEZPTz5me0hBGS2LOZBzOnJjR8kPFaG/5k9loT3lM9Jb/4dNTHtPRU/6E1t17ujSmr+dUaUx/d3k/03tOlMZclhEzo2v8mFnd/yvtY1b3W6Ux7+w6Vhozs7O8n4Guo+X76iiPGejMyCvjxcTsjnNP+GlXHPpHaaPEbxmZmRnggmBmZokLgpmZAS4IZmaWuCCYmRnggmBmZokLgpmZAS4IZmaWZBUESYsk7ZM0LGn1GNt7JW1O27dLuqpm291p/T5JC2vWPyrpFUm7GpGImZldmtKCIKkTeBBYDAwByyUN1YWtAN6MiGuA+4H7UtshYBkwH1gEPJT6A/h+WmdmZlNAzhXCDcBwRByIiJPAJmBpXcxSYENa3grcLElp/aaIOBERB4Hh1B8R8VvgjQbkYGZmDZBTEK4EXqp5PJLWjRkTEaeBI8BAZttxSVopaaeknafi+ESampnZBEz5SeWIWB8RCyJiQbf6Jns4ZmZtK6cgHAbm1jyek9aNGSOpC5gJvJ7Z1szMpoCcgvAsME/SoKQeiknibXUx24Db0vKtwNMREWn9svQppEFgHrCjMUM3M7NGKi0IaU7gTuApYC+wJSJ2S1ojaUkKewQYkDQM3AWsTm13A1uAPcAvgVURcQZA0kbgD8C1kkYkrWhsamZmNhFZ/yAnIp4Enqxb99Wa5ePAZy7Q9l7g3jHWL5/QSM3MrFJTflLZzMyawwXBzMwAFwQzM0tcEMzMDHBBMDOzxAXBzMwAFwQzM0tcEMzMDHBBMDOzxAXBzMwAFwQzM0tcEMzMDHBBMDOzxAXBzMwAFwQzM0tcEMzMDMgsCJIWSdonaVjS6jG290ranLZvl3RVzba70/p9khbm9mlmZs1VWhAkdQIPAouBIWC5pKG6sBXAmxFxDXA/cF9qO0TxP5jnA4uAhyR1ZvZpZmZNlHOFcAMwHBEHIuIksAlYWhezFNiQlrcCN0tSWr8pIk5ExEFgOPWX06eZmTVRTkG4Enip5vFIWjdmTEScBo4AA+O0zenTzMyaqGuyB1BG0kpgZXp44tenNu4at8G/Kx9So7wLeG2yB9FA7ZRPO+UC7ZVPO+UCzcnnPbmBOQXhMDC35vGctG6smBFJXcBM4PWStmV9AhAR64H1AJJ2RsSCjDFPee2UC7RXPu2UC7RXPu2UC0y9fHLeMnoWmCdpUFIPxSTxtrqYbcBtaflW4OmIiLR+WfoU0iAwD9iR2aeZmTVR6RVCRJyWdCfwFNAJPBoRuyWtAXZGxDbgEeAxScPAGxRP8KS4LcAe4DSwKiLOAIzVZ+PTMzOzXCpeyLcGSSvTW0gtr51ygfbKp51ygfbKp51ygamXT0sVBDMzq45vXWFmZsAkFoR2ux1Go/ORNFfSM5L2SNot6YutmkvNtk5Jz0t6ovosztlvFefaLElbJb0oaa+kD7RwLl9O59guSRsl9TUjl7Tvi8pH0kD6+zgqaV1dm/dJ+mtq84AktWIukqZJ+nk6x3ZL+mblSURE038oJpL3A1cDPcBfgKG6mM8D303Ly4DNaXkoxfcCg6mfzpw+Wyyfy4HrU8w7gL81I58qcqlpdxfwY+CJVj7X0rYNwB1puQeY1Yq5UHwh9CDQn+K2ALe3wLGZDnwY+Bywrq7NDuBGQMAvgMWtmAswDfhozTn2u6pzmawrhHa7HUbD84mIlyPiOYCI+C+wl+Z8m7uKY4OkOcAngIebkEOthucjaSbwEYpP1xERJyOiGV+JrOTYUHzasF/Fd4imAf+qOI+zLjqfiDgWEb8HjtcGS7ocuCwi/hjFM+kPgE9XmkWh4blExFsR8UxaPgk8R/GdrcpMVkFot9thVJHP29Kl5XXA9gaO+UKqymUt8BVgtPFDHlcV+QwCrwLfS2+BPSxpejXDH3ucdeMZMyYnl4g4DHwb+CfwMnAkIn5VyejPdyn5jNfnSEmfVagil7dJmgV8CvjNJY90HJ5UnuIkzQB+CnwpIv4z2eO5GJI+CbwSEX+a7LE0SBdwPfCdiLgOOAa05C3cJc2meOU6CFwBTJf02ckdldVKV24bgQci4kCV+5qsgjCR22Gc/YWU3Q4jp8+qVJEPkropisGPIuLxSkZ+vipy+RCwRNIhikvpj0n6YRWDH0MV+YwAIxFx9optK0WBqFoVuXwcOBgRr0bEKeBx4IOVjP58l5LPeH3Wvq3SrOeBKnI5az3w94hY24Bxjq/qyZYLTMB0AQcoXpWcnYCZXxezinMnYLak5fmcOzl2gGJCp7TPFstHFO9/rm31Y1PX9iaaO6lcST4UE3zXpuWvAd9qxVyA9wO7KeYORPEe9xem+rGp2X475ZPKt7RwLt+geFHY0ZRj0oydXOAXeAvFJ2f2A/ekdWuAJWm5D/gJxeTXDuDqmrb3pHb7qJl1H6vPVs2H4lMHAbwA/Dn9VH5iV3VsarbfRBMLQoXn2nuBnen4/AyY3cK5fB14EdgFPAb0tsixOURxq5yjFFdtQ2n9gpTLfmAd6Qu4rZYLxVVGUHyg5OxzwB1V5uBvKpuZGeBJZTMzS1wQzMwMcEEwM7PEBcHMzAAXBDMzS1wQzMwMcEEwM7PEBcHMzAD4P2YY2EixWfAwAAAAAElFTkSuQmCC\n",
      "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)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "#saving the solution to a VTK file\n",
    "vtkfile = File('solution/sol.pvd')\n",
    "vtkfile << u"
   ]
  },
  {
   "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$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ Product(FloatValue(6993006.993006993), Abs(Indexed(Grad(Product(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 156), FiniteElement('Lagrange', triangle, 1)), 196), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 156), FiniteElement('Lagrange', triangle, 1)), 196))), MultiIndex((FixedIndex(0),))))),\n",
       "       Product(FloatValue(6993006.993006993), Abs(Indexed(Grad(Product(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 156), FiniteElement('Lagrange', triangle, 1)), 196), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 156), FiniteElement('Lagrange', triangle, 1)), 196))), MultiIndex((FixedIndex(1),)))))], dtype=object)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "Q = 1/rho*np.abs(grad(dot(voltage,voltage)))\n",
    "Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "ename": "RuntimeError",
     "evalue": "Don't know how to plot given object:\n  [ Product(FloatValue(6993006.993006993), Abs(Indexed(Grad(Product(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 156), FiniteElement('Lagrange', triangle, 1)), 196), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 156), FiniteElement('Lagrange', triangle, 1)), 196))), MultiIndex((FixedIndex(0),)))))\n Product(FloatValue(6993006.993006993), Abs(Indexed(Grad(Product(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 156), FiniteElement('Lagrange', triangle, 1)), 196), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 156), FiniteElement('Lagrange', triangle, 1)), 196))), MultiIndex((FixedIndex(1),)))))]\nand projection failed:\n  'numpy.ndarray' object has no attribute 'ufl_domain'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m/usr/local/lib/python3.6/dist-packages/dolfin/common/plotting.py\u001b[0m in \u001b[0;36mplot\u001b[0;34m(object, *args, **kwargs)\u001b[0m\n\u001b[1;32m    427\u001b[0m                          \"piecewise linears.\")\n\u001b[0;32m--> 428\u001b[0;31m             \u001b[0mobject\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mproject\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobject\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmesh\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmesh\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    429\u001b[0m             \u001b[0mmesh\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mobject\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfunction_space\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmesh\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.6/dist-packages/dolfin/fem/projection.py\u001b[0m in \u001b[0;36mproject\u001b[0;34m(v, V, bcs, mesh, function, solver_type, preconditioner_type, form_compiler_parameters)\u001b[0m\n\u001b[1;32m     94\u001b[0m             \u001b[0;31m# Otherwise try extracting function space from expression\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 95\u001b[0;31m             \u001b[0mV\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_extract_function_space\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmesh\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     96\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.6/dist-packages/dolfin/fem/projection.py\u001b[0m in \u001b[0;36m_extract_function_space\u001b[0;34m(expression, mesh)\u001b[0m\n\u001b[1;32m    150\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mmesh\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 151\u001b[0;31m         \u001b[0mdomain\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mexpression\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mufl_domain\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    152\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mdomain\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'numpy.ndarray' object has no attribute 'ufl_domain'",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-34-86a109355051>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mplot\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[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/usr/local/lib/python3.6/dist-packages/dolfin/common/plotting.py\u001b[0m in \u001b[0;36mplot\u001b[0;34m(object, *args, **kwargs)\u001b[0m\n\u001b[1;32m    432\u001b[0m             \u001b[0mmsg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"Don't know how to plot given object:\\n  %s\\n\"\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    433\u001b[0m                   \u001b[0;34m\"and projection failed:\\n  %s\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobject\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0me\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 434\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    435\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    436\u001b[0m     \u001b[0;31m# Plot\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mRuntimeError\u001b[0m: Don't know how to plot given object:\n  [ Product(FloatValue(6993006.993006993), Abs(Indexed(Grad(Product(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 156), FiniteElement('Lagrange', triangle, 1)), 196), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 156), FiniteElement('Lagrange', triangle, 1)), 196))), MultiIndex((FixedIndex(0),)))))\n Product(FloatValue(6993006.993006993), Abs(Indexed(Grad(Product(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 156), FiniteElement('Lagrange', triangle, 1)), 196), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 156), FiniteElement('Lagrange', triangle, 1)), 196))), MultiIndex((FixedIndex(1),)))))]\nand projection failed:\n  'numpy.ndarray' object has no attribute 'ufl_domain'"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Collecting vtkplotter\n",
      "\u001b[31m  ERROR: Could not find a version that satisfies the requirement vtkplotter (from versions: none)\u001b[0m\n",
      "\u001b[31mERROR: No matching distribution found for vtkplotter\u001b[0m\n",
      "\u001b[33mWARNING: You are using pip version 19.1, however version 21.2.1 is available.\n",
      "You should consider upgrading via the 'pip install --upgrade pip' command.\u001b[0m\n"
     ]
    }
   ],
   "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
}