diff --git a/FEMII/00_intro_to_fenics.ipynb b/FEMII/00_intro_to_fenics.ipynb
index 84e7ed1..2aabfee 100644
--- a/FEMII/00_intro_to_fenics.ipynb
+++ b/FEMII/00_intro_to_fenics.ipynb
@@ -1,49 +1,656 @@
 {
  "cells": [
   {
    "cell_type": "markdown",
    "metadata": {
-    "collapsed": true,
     "pycharm": {
      "name": "#%% md\n"
     }
    },
    "source": [
-    "# Introduction to FEnICS\n",
+    "# Introduction to FEnICS"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "https://link.springer.com/book/10.1007%2F978-3-642-23099-8\n",
     "\n",
+    "https://link.springer.com/book/10.1007%2F978-3-319-52462-7"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from fenics import *\n",
+    "import matplotlib.pyplot as plt\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Mesh Creation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "FEniCS provides a number of [built-in meshes](https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/built-in_meshes/python/documentation.html) for simple geometries: "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 38,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAD4CAYAAAA0L6C7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAARzElEQVR4nO3d/4td9Z3H8ed7kjEtmSR+m62amM2CplgrG5dBheDKSra4rbQLIjSQBqEhv+yCNQtl+1PoP1D6y/6wQ5XNxpCgSUsWjV0DKjHQambs2NVOXFTMFyzMldqMIZidmPf+8L53nJnMlzP33nM/58vrCUPOzb05981JHjn3nHvnjLk7Sqnq1Jd6AKVUdxNqpSqWUCtVsYRaqYol1EpVrJV5rPTmm2/2TZs25bFqpRQwOjr6ibsPzndfLqg3bdrEyMhIHqtWSgFmdmah+/TyW6mKJdRKVSyhVqpiCbVSFUuolapYmc5+m9lHwGfAF8AVdx/KcyilVPst5y2tv3P3T3KbRCnVlXJ5n3qp3nvvPQ4dOsSWLVtYu3ZtihEK0YkTJwC0HZrb4cEHH8TMEk+TrtZ22Lt3b0fryYragZfNzIF/d/fhuQ8ws93AboCNGzcuurL3338fgLGxseXMWtm0HaLXX3899QiFaGpqiv7+/rb/fFbUW939YzP7C+C4mZ129xMzH9CEPgwwNDS06JUX7rzzTkZGRti1axfr169va/Ayd+HCBfbt28enn34KwJ49e1izZk3iqXrfxMQE+/bt49KlS0Dne6iydubMGQ4cOMDU1BR9fX0dgYaMZ7/d/ePmrxPAr4D7OnrWGtcCfenSJbZs2ZJ6nGS1QPf19XH33XenHidZLdDr1q3jjjvuYGBgoON1LonazFab2ZrWMvAt4J2On7mGzQS9Y8eOWr5Kgdmgn3jiCW666abUIyVpJuidO3d2BTRk21N/DThpZm8DbwIvuvuvu/LsNWou6A0bNqQeKUkCHc0F3c3DryWPqd39Q+Cvu/aMNUygI4GO8gQN+kRZ7gl0JNBR3qBBqHNNoCOBjnoBGoQ6twQ6EuioV6BBqHNJoCOBjnoJGoS66wl0JNBRr0GDUHc1gY4EOkoBGoS6awl0JNBRKtAg1F1JoCOBjlKCBqHuOIGOBDpKDRqEuqMEOhLoqAigQajbTqAjgY6KAhqEuq0EOhLoqEigQaiXnUBHAh0VDTQI9bIS6EigoyKCBqHOnEBHAh0VFTQIdaYEOhLoqMigQaiXTKAjgY6KDhqEetEEOhLoqAygQagXTKAjgY7KAhqEet4EOhLoqEygQaivSaAjgY7KBhqEelYCHQl0VEbQINTTCXQk0FFZQYNQAwLdSqCjMoMGoRboZgIdlR001By1QEcCHVUBNNQYtUBHAh1VBTTUFLVARwIdVQk01BC1QEcCHVUNNNQMtUBHAh1VETTUCLVARwIdVRU0LAO1ma0ws9+Z2Qt5DpRHAh0JdFRl0LC8PfWTwHheg+SVQEcCHVUdNGREbWYbgO8Av8h3nO4m0JFAR3UADdn31D8HfgxcXegBZrbbzEbMbKTRaHRluE4S6Eigo7qAhgyozexRYMLdRxd7nLsPu/uQuw8NDg52bcB2EuhIoKM6gYZse+qtwHfN7CPgEPCwmT2b61QdJNCRQEd1Aw0ZULv7T9x9g7tvAr4PvOLuO3KfrI0EOhLoqI6goULvUwt0JNBRXUEDrFzOg939NeC1XCbpIIGOBDqqM2iowJ5aoCOBjuoOGkqOWqAjgY4EOiotaoGOBDoS6C8rJWqBjgQ6EujZlQ61QEcCHQn0tZUKtUBHAh0J9PyVBrVARwIdCfTClQK1QEcCHQn04hUetUBHAh0J9NIVGrVARwIdCXS2CotaoCOBjgQ6e4VELdCRQEcCvbwKh1qgI4GOBHr5FQq1QEcCHQl0exUGtUBHAh0JdPsVArVARwIdCXRnJUct0JFARwLdecu68km3m5yc5MiRI7UH3Wg0OHLkiEALdFdKuqd+7rnnag8aYP/+/bUHDQh0l0qypz558uT08tq1axkdHWV0dNHLileysbGx6eUbbrhh1napU63tMDU1xfXXX88rr7ySeKI0tbaDu2Nmba8nCepz585NL1++fJkPP/wwxRhJm5ycnHX7woULXLhwIdE06Zq7HSYmJpiYmEg0TbpmbocrV67Q39/f9rqSoN6+fTsHDx5k165drF+/PsUISWudFLt06RIAe/bsqeXLzbNnz/Lss88yNTUFwN69exNPlKbx8XEOHz7M1atXWbVqVUegoQBnv+vWzLPcDzzwQOpxktUCvW7dOu69997U4ySrBfq2227jrrvuYtWqVR2vU6h7mN62imaCrvNJsZmgd+zY0RXQINQ9S6AjgY7yAg1C3ZMEOhLoKE/QINS5J9CRQEd5gwahzjWBjgQ66gVoEOrcEuhIoKNegQahziWBjgQ66iVoyIDazL5iZm+a2dtm9q6Z/TTXiUqeQEcCHfUaNGT7RNll4GF3v2hm/cBJM3vJ3X+b82ylS6AjgY5SgIYMqN3dgYvNm/3NL89zqDIm0JFAR6lAQ8ZjajNbYWZjwARw3N3fmOcxu81sxMxGGo1Gt+csdAIdCXSUEjRkRO3uX7j7FmADcJ+ZfXOexwy7+5C7Dw0ODnZ7zsIm0JFAR6lBwzLPfrv7n4HXgEdymaZkCXQk0FERQEO2s9+DZnZ9c/mrwDbgdN6DFT2BjgQ6KgpoyHb2+1Zgn5mtIP4TeM7dX8h3rGIn0JFAR0UCDdnOfv8eqO83vM5JoCOBjooGGvSJsmUl0JFAR0UEDUKdOYGOBDoqKmgQ6kwJdCTQUZFBg1AvmUBHAh0VHTQI9aIJdCTQURlAg1AvmEBHAh2VBTQI9bwJdCTQUZlAg1Bfk0BHAh2VDTQI9awEOhLoqIygQainE+hIoKOyggahBgS6lUBHZQYNQi3QzQQ6KjtoqDlqgY4EOqoCaKgxaoGOBDqqCmioKWqBjgQ6qhJoqCFqgY4EOqoaaKgZaoGOBDqqImioEWqBjgQ6qipoqAlqgY4EOqoyaKgBaoGOBDqqOmioOGqBjgQ6qgNoqDBqgY4EOqoLaKgoaoGOBDqqE2ioIGqBjgQ6qhtoqBhqgY4EOqojaKgQaoGOBDqqK2ioCGqBjgQ6qjNoqABqgY4EOqo7aCg5aoGOBDoS6Ki0qAU6EuhIoL+slKgFOhLoSKBntyRqM7vdzF41s3Eze9fMnuzFYAsl0JFARwJ9bUv+0HngCvAv7v6Wma0BRs3suLv/IefZrkmgI4GOBHr+ltxTu/sf3f2t5vJnwDiwPu/B5ibQkUBHAr1wyzqmNrNNwL3AG/Pct9vMRsxspNFodGe6ZgIdCXQk0IuXGbWZDQBHgB+5++Tc+9192N2H3H1ocHCwawMKdCTQkUAvXSbUZtZPgD7g7r/Md6QvE+hIoCOBzlaWs98GPA2Mu/vP8h8pEuhIoCOBzl6WPfVW4AfAw2Y21vz6dp5DCXQk0JFAL68l39Jy95OA9WAWQKBbCXQk0MuvUJ8oE+hIoCOBbq/CoBboSKAjgW6/QqAW6EigI4HurOSoBToS6EigOy/LZ79zq9FocPz4cYE+e5ajR48KtEB3paR76qNHj9YeNMDhw4drDxoQ6C6VZE/94osvfjnAypW89NJLKcZI3gcffDC93NfXx9GjRxNOk67Wdrh69SoAzz//fMpxkjVzO/T1tb+/TYL6888/n14eGBjg8uXLKcZI2vnz52fdvu6662q/HdauXQtQ++3g7h2tKwnqxx57jIMHD7Jr1y7Wr+/5d3Em78yZMxw4cICpqSkA9uzZU8uX3a1j6NYe+qmnnko8UZpOnTo1jXr16tWsWLGio/UlP/tdt1qg161bx0MPPZR6nGTNPCl2//33px4nWadOneLYsWNs3ryZe+65p2PQINQ9bSbonTt3MjAwkHqkJOksdzQT9OOPP94V0CDUPWsu6Dq+3AaBbjUX9MqV3TsSFuoeJNCRQEd5ggahzj2BjgQ6yhs0CHWuCXQk0FEvQINQ55ZARwId9Qo0CHUuCXQk0FEvQYNQdz2BjgQ66jVoEOquJtCRQEcpQINQdy2BjgQ6SgUahLorCXQk0FFK0CDUHSfQkUBHqUGDUHeUQEcCHRUBNAh12wl0JNBRUUCDULeVQEcCHRUJNAj1shPoSKCjooEGoV5WAh0JdFRE0CDUmRPoSKCjooIGoc6UQEcCHRUZNAj1kgl0JNBR0UFDth86/4yZTZjZO70YqEgJdCTQURlAQ7Y99X8Aj+Q8R+ES6Eigo7KAhgyo3f0E8KcezFKYBDoS6KhMoKGLx9RmttvMRsxspNFodGu1PU+gI4GOygYauoja3YfdfcjdhwYHB7u12p4m0JFAR2UEDTr7PZ1ARwIdlRU0CDUg0K0EOiozaMj2ltZB4DfA183svJn9MP+xepdARwIdlR00ZPipl+6+vReDpEigI4GOqgAaavzyW6AjgY6qAhpqilqgI4GOqgQaaohaoCOBjqoGGmqGWqAjgY6qCBpqhFqgI4GOqgoaaoJaoCOBjqoMGmqAWqAjgY6qDhoqjlqgI4GO6gAaKoxaoCOBjuoCGiqKWqAjgY7qBBoqiFqgI4GO6gYaKoZaoCOBjuoIGiqEWqAjgY7qChoqglqgI4GO6gwaKoBaoCOBjuoOGkqOWqAjgY4EOiotaoGOBDoS6C8rJWqBjgQ6EujZlQ61QEcCHQn0tZUKtUBHAh0J9PyVBrVARwIdCfTClQK1QEcCHQn04hUetUBHAh0J9NIVGrVARwIdCXS2CotaoCOBjgQ6e4VELdCRQEcCvbwKh1qgI4GOBHr5FQq1QEcCHQl0exUGtUBHAh0JdPsVArVARwIdCXRnJUct0JFARwLdeZlQm9kjZvaemb1vZv/arSc/e/asQAOnT58WaAS6Wy2J2sxWAP8G/APwDWC7mX2jG0/+8ssv1x40wLFjx2oPGhDoLpVly90HvO/uHwKY2SHge8Af2n3SkydPTi9fvHiR/fv3t7uqUtdoNKaXJycnefrppxNOk66Z26HRaDA8PJxwmnTN3A6dlAX1euDcjNvngfvnPsjMdgO7ATZu3LjoCjdu3Mi5c+cYGBhY8rFVrvWXeOONN3LLLbckniZdre2wefPmWu+he4na5vk9v+Y33IeBYYChoaFr7p/Ztm3b2LZtW6YBlVLLK8uJsvPA7TNubwA+zmccpVSnZUF9CrjTzP7KzK4Dvg/8V75jKaXabcmX3+5+xcz+GfhvYAXwjLu/m/tkSqm2ynRWwt2PAcdynkUp1YWSf6JMKdXdhFqpiiXUSlUsoVaqYpn7op8TaW+lZg3gzBIPuxn4pOtP3r00X/sVeTYo9nxZZ/tLdx+c745cUGfJzEbcfSjJk2dI87VfkWeDYs/Xjdn08lupiiXUSlWslKiL/v11mq/9ijwbFHu+jmdLdkytlMonvfxWqmIJtVIVKwnqvC5k2I3M7BkzmzCzd1LPMjczu93MXjWzcTN718yeTD3TzMzsK2b2ppm93Zzvp6lnmpuZrTCz35nZC6lnmZuZfWRm/2NmY2Y20vZ6en1M3byQ4f8Cf09cgOEUsN3d277mWTczs78FLgL/6e7fTD3PzMzsVuBWd3/LzNYAo8A/FmjbGbDa3S+aWT9wEnjS3X+beLTpzGwPMASsdfdHU88zMzP7CBhy944+GJNiTz19IUN3/z+gdSHDQuTuJ4A/pZ5jvtz9j+7+VnP5M2CcuIZcIfLoYvNmf/OrMGdizWwD8B3gF6lnybMUqOe7kGFh/mGWJTPbBNwLvJF2ktk1X96OARPAcXcv0nw/B34MXE09yAI58LKZjTYv5NlWKVBnupChWjgzGwCOAD9y98nU88zM3b9w9y3EtezuM7NCHMKY2aPAhLuPpp5lkba6+98Q19j/p+ah4LJLgVoXMuyg5rHqEeCAu/8y9TwL5e5/Bl4DHkk8SqutwHebx62HgIfN7Nm0I83O3T9u/joB/Io4VF12KVDrQoZt1jwR9TQw7u4/Sz3P3Mxs0Myuby5/FdgGnE47VeTuP3H3De6+ifg394q770g81nRmtrp58hMzWw18C2jrHZieo3b3K0DrQobjwHNFupChmR0EfgN83czOm9kPU880o63AD4i9zFjz69uph5rRrcCrZvZ74j/v4+5euLeOCtrXgJNm9jbwJvCiu/+6nRXpY6JKVSx9okypiiXUSlUsoVaqYgm1UhVLqJWqWEKtVMUSaqUq1v8D7m0r2CWQXKgAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "mesh_2d = RectangleMesh(Point(0,0),Point(5,5), nx=3, ny=3)\n",
+    "plot(mesh_2d)\n",
+    "plt.show()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 39,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[0.        , 0.        ],\n",
+       "       [1.66666667, 0.        ],\n",
+       "       [3.33333333, 0.        ],\n",
+       "       [5.        , 0.        ],\n",
+       "       [0.        , 1.66666667],\n",
+       "       [1.66666667, 1.66666667],\n",
+       "       [3.33333333, 1.66666667],\n",
+       "       [5.        , 1.66666667],\n",
+       "       [0.        , 3.33333333],\n",
+       "       [1.66666667, 3.33333333],\n",
+       "       [3.33333333, 3.33333333],\n",
+       "       [5.        , 3.33333333],\n",
+       "       [0.        , 5.        ],\n",
+       "       [1.66666667, 5.        ],\n",
+       "       [3.33333333, 5.        ],\n",
+       "       [5.        , 5.        ]])"
+      ]
+     },
+     "execution_count": 39,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "mesh_2d.coordinates()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 40,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mesh_3d = BoxMesh(Point(0,0,0), Point(1,1,1), 5, 5, 5)\n",
+    "#plot(mesh_3d)\n",
     "\n",
+    "# Save solution in VTK format\n",
+    "file = File(\"mesh_3d.pvd\");\n",
+    "file << mesh_3d;"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "An interface with [CGAL](https://www.cgal.org) also allows mesh generation for more complex geometries. \n",
+    "See this [example](https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/built-in_meshes/python/documentation.html)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Function Space"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A *Function Space* in FEniCS defines an interpolation scheme over a discretized domain, i.e. a mesh."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 41,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "V = FunctionSpace(mesh_2d, \"Lagrange\", 1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 42,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]"
+      ]
+     },
+     "execution_count": 42,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "V.dofmap().dofs()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Change the degree of interpolation and observe change in degreees of freedom."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The following types of discretizations exist:\n",
+    "- `FunctionSpace` for scalar fields\n",
+    "- `VectorFunctionSpace` for vector fields\n",
+    "- `TensorFunctionSpace` for tensor fields"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 44,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[0,\n",
+       " 1,\n",
+       " 2,\n",
+       " 3,\n",
+       " 4,\n",
+       " 5,\n",
+       " 6,\n",
+       " 7,\n",
+       " 8,\n",
+       " 9,\n",
+       " 10,\n",
+       " 11,\n",
+       " 12,\n",
+       " 13,\n",
+       " 14,\n",
+       " 15,\n",
+       " 16,\n",
+       " 17,\n",
+       " 18,\n",
+       " 19,\n",
+       " 20,\n",
+       " 21,\n",
+       " 22,\n",
+       " 23,\n",
+       " 24,\n",
+       " 25,\n",
+       " 26,\n",
+       " 27,\n",
+       " 28,\n",
+       " 29,\n",
+       " 30,\n",
+       " 31]"
+      ]
+     },
+     "execution_count": 44,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "U = VectorFunctionSpace(mesh_2d, \"Lagrange\", 1)\n",
+    "U.dofmap().dofs()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Functions"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "*Functions* live on a function space. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 45,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "v = Function(V)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Their values are defined at the integration points and interpolated at all other points according to the integrationscheme of their function space."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 46,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])"
+      ]
+     },
+     "execution_count": 46,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# values at integration points\n",
+    "v.compute_vertex_values()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can now discretize analytical functions:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 47,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "scalar_field_expression = Expression(\"exp( - pow(x[0]-x0, 2) / pow(sigma_x, 2) - pow(x[1]-y0, 2) / pow(sigma_y, 2) )\", \n",
+    "                                    x0=3, y0=3, sigma_x=1, sigma_y=1, \n",
+    "                                    degree = 1)\n",
+    "scalar_field_discrete = project(scalar_field_expression, V)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 48,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/Users/dabler/anaconda2/envs/fenicsenv/lib/python3.7/site-packages/dolfin/common/plotting.py:152: UserWarning: The following kwargs were not used by contour: 'scalarbar'\n",
+      "  return ax.tricontourf(mesh2triang(mesh), C, levels, **kwargs)\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.tri.tricontour.TriContourSet at 0x11ce11828>"
+      ]
+     },
+     "execution_count": 48,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD8CAYAAABaQGkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAbBklEQVR4nO2dX6hl113Hv79z5v7LJJORGGtsAhGUQCk00SEWBkRLrbEN1UcD7VPhvqikVCnGF1HwtdSHPnixQUtrSjENSJXqgA0hUNNmYtqmnSqlBowJDDGZzIy999xzz/n5cM5O9tl37b1/a++19vqzfx+4zJ179z1n3XPPZ6+1vvu31iZmhqIoeTIJ3QBFUfyhgitKxqjgipIxKriiZIwKrigZo4IrSsackRxERC8DuAFgAeCEmS/4bJSiKG4QCb7m15n5dW8tURTFOTpEV5SMIUklGxH9F4A3ATCAv2LmA8Mx+wD2AWCK6S/fQrc7bqoyFoio88+OpTLzBr/xOjPf2XacVPCfY+ZXiehnAFwC8AfM/Ezd8ecmd/D7z/ymVYMVBQAme7uhmxCc5eFR6zGX5k9clmRhojk4M7+6/vcqET0F4EEAtYIrSh0Sgcmj5CyQJzRtr5HkBFDQKjgRnQUwYeYb688/BODPxc+gZI3LHten2EM+hy+6nJwkPfi7ADy1nhedAfB3zPx162dSkiHEMPlt8XZ3Bn/utzmahXtuT7QKzsw/BvC+AdqieCTmuW0Ucsfw/HX0OPHYXAdXIiNmaaWckjvm3ymB+XsVFTwycpBWSlJyA8O3z8EJRQUfmDEJXMdG0FWRm3e2ArSoPzSbh26CERV8QFTumvl2Re7l7ubbcnJ0Mkjb+uDyxOTyZKGCD4TK3Twkr5O77mupYXOS4p0tZ5Kn/8olgMpdL3e55ytEXu7I35aTWfy9OyA/Sbkerajgnhm73JL5dvnNbyN3l+NjxOdJKv1XJ2Jyltu6IqxF7rKoi91p/waumR4tnD2WL5Y7Z7xJroJ7IkW5vZRxCsK0Qm6XYhf4eEyX+D4BqeAeiEnuoLXXgjDNJPdiZ5htCqaz5SDPExIV3DFDyJ3EggmLMK2QeyixC4Z+vjJDnVxUcId0lTsJYW2wCNOqci+2u2/24JLpcR4bR6jgjpDKnZ3MZSzm20C8cgP+2mJ14tjb7V2uqoI7QOVG5zCtKvdiJx7Jm5jO0ujhVfCBGIXcFmFaef6bmtxAt7aGOCmo4D0JvQVRcHqEaeVhcCFMSpLXEVPvroL3QOV2E6ZV5W6a/6YQflVPUmXhF9s06O+ggndktHJXdz3pEaaVRZDILfl+TMRwMlLBO5Cd3F22KupYvFKVe0Py4nvb9s0BgOlxt5/LGRXckiTk9r23WMfilTq5N4brHeXu+7MucX6i2d3pvC+bCm5BULlj2RDQcZhm6rWXnjZ1mcS56YpXVHAhXuSORVopPcK0pvl2Ibcvsd9uo+fHj/EEooILcCp3alIXOAjTmubbhXyLQC/PNL8t0QGo4K2MXu4BwrTQcrt4btsTxGJnMsiCExW8AZXbf5hWlXsZSVBWZpJwOq+C16By+yleMYVpMcsNNLcrdvlVcAMq93BhWlVu30GYK0IEarS3a30DQhW8gso9fJhWlXu5vaoAmxynU7UWKyp4iVHJ3fB7uAjTbOfbVbmrn8dC00lnuRXfpTIVfE0WcvcssgkRppWH5IXQyx2Z2JOZ9vBtqOCIXO6Byl5DhGnGXlsot+2xfUn1ZDJ6wYPIHbpWvYLrlWA2YZqp115u9Rd3Mk9DyMXu1OvWyaMWXOUetnhlKLldPg6QzsnCxGgFV7njCNMKuctC8rb/Ci86Drdl8pCMUnCne5cnKLfrMM1mvr36mlnuIcQukDxXrCeByd4uIEzrxYIT0RTA8wD+h5kf7ta08DjdATVhuV2HaTbz7dVzmuWmrfB3G+F5P7EX2/FsPmHTgz8K4AqAc57a4h2VO84wLSa5c0MkOBHdDeAjAP4CwKe8tsgTKnc8YZppvl3IPR1wmF5lEemQvA/SHvyzAD4N4La6A4hoH8A+AOzilv4tc4jK3X2Pcl9hWnkOHIPcofF1C+FWwYnoYQBXmfkyEf1a3XHMfADgAADOTe6IpsZwDHKXQ7M6Yg3TqnJvbfu5T3YT8+P2fm65w0kWu0h68IsAPkpEHwawC+AcEX2RmT/mt2n9SVVuibBSUgjTqnLvbJ+OiGfHiSwzi4xWwZn5MQCPAcC6B/8jlbsGwWO5lLcNX2Fa1+KVpvl2k9xNX+/DGE4aWV4HV7nTCdPKQ/JC4r0t+2H64TzLt3JvrF4VZn4awNNeWuIIlTudMM3Ua3eRu8vPDXFCGPo2RSayOu2p3O1yu14Jtvpas9yS+XZZ0L0tt1Uih3Pznks72/Psh+nZCD52uUNsq+QqTCvkdi12Qflx62TPlSwEj1XuoUg5TKvKfW7bzwbl148j32FnzXL3DCZH7i4VJi94zHIP0Xv7DtO6rgSzCdN8y+2S5TYntVdc0oKr3OG3VeobplXlvnXLreQ352n03L5IVnCVO58wzZfcSqKC5yp3OSgTHR9hmGYz3wZOy33bGbt9v+u4cRJPRrLYIUxnYS6XJSd46nLbSmx8jAB7lNuEabbzbddy27K1fSKqR0+RpH4rlbtbmOZzW6U+cpeH5IXc5xxIfj2i3js0yQiucucZppV77ULun9r6v5pX4B3enJ9tPSZleGcLNOtff5+E4Cr3cGGa6+KV1efNcpd7bYncdccNKX2MdzExEb3gTuWWEpHcIbZVWj3vMGGaSe7bp4fVl6GRtxZ7p75225mjqIK2UEQtuHO5Jb13onKnGKaZhuS2civNRCu4yp1GmGZbvFInd1ns89OfQMK1RVxbg7UR4nJZlIKr3Glsq9Rlvg3Uyy0VW5ETneAq97jCtKrc56ftIdu1hZ8wbbq9zG5n1agEH6vc5bk2MFyY5nolGNAut2m+bSP36NjbBQ671wZEI/hY5K7KXCbFbZX6hmlVue+YNEv+v8u8r3+7JgrBVW5/xSvl78UapknljpHlNjCJ5DZFJoILrnKPO0yryn1+Yl5Rdm3ZbdnnrVuzUS8ZDS64hDHJHese5a6KV0zz7arct5dyrrfGe7MTJwQVfPB7dCcid+zbKgFuwrTykNwkd64MudtqMMFTlltKH7lDVKYBdnum+ZL79skO3lqG2/yBtpa9byEcC3n8FhkSaoOAGKiTu0uCbjP/Lq6B28htE7AVf1Or3rvHJTJABU+aqaCTc33DPMnupJJVXRJZi/l3OWArilyKMtVioUnxnMVa8PJCk6LNxZbJxU0Pij3Ry5s9FHLTWvbJfPX6Fa9jseFieSVZ8XcovjY9Lv5dC204WU9n69HT0WL9+KuRUrGj6sZS0aPVE/Ba9qWF9Cp4ZEh6bskyRZudPyXVW5L9xPtutFAO1Kpyl08IUrnLvXeb3KbXoE3ucu/dJrep926VuyQyd+zJVfBIkAzbpo6vt0qGoq5v8dNWZvrWsl7u8s9K5S6PONrkrvbeQLvc5VGUVO6i9149fovcR5vDNJveG1DBg1P+Y9ceYzFnk/TuJOixbW7pI1l33bTy69py59R17jq5TY/TJnd59CGVu+i9gXa5N4brQrmL3htol7tr7w2o4MEo/4FdIAl7ym9aF0gCLNNmDHUUorfJXe29gXa5yyMRqdzl/KJN7vLoSip30XsDfuQGVHCvTGYnGx+2SObjroM2ye6iLoO24sP0PUAeqgHtcpdHJVK5y1lGm9ymkVab3Bu3KXIsN5BIJVsuFH/UpuvjwErs8nVw42PNNyvZ+rI4nmxcDzdxON/2doNAk+S+E/MCH4l5+eTcJzHvi/bgEeE6aJMk6a6DNkmSfm1xduNDgo/EvPy755CYm9AePAKms+XGQpO+SHp3Op5sLDTpy42T3dYbF1xb3GLctaVJ8nKoNrbE3AXagyeCpHd3HbRJknTboO3a4paNjyaqQ3MgrsS8IJbE3ERrD05EuwCeAbCzPv7vmflPnbZipEyPFhubPBiPEczHp7PNuvS+zI/PbGzyYOL68Y6T2/2aJD8//YmzxLzAR2K+MQSPIDE3IenBZwA+wMzvA3A/gIeI6P3OWzIiJIm661p0SZLuqqKtoCzlW4u9jY8mXCbmXcpQV5+v/rUpQw2dmJto7cGZmQHcXP93a/0x3pUQnpke88aa8CYkc+3JMW2sCTfB88nGmvC+XD/Zbb3HWJ3kxeozl4l5+cQVVWI+AKKQjYimAC4D+AUAn2Pm5wzH7APYB4BdpLVfdWpMjzd3dDHh+jLa4fzMxo4uJkxBm+l6eNPticrix5CYF3hJzC1DNdsyVUAYsjHzgpnvB3A3gAeJ6L2GYw6Y+QIzX9ii/G8ZszHccoCrklWb5YuuSlZNQVvT5bI352dPfZiOKT9OqMS8Sxnq6vGHT8xNWF0mY+ZrRPQ0gIcAvOSlRQlRlbzLtsmugjYbJnPa2HDRhCRoa8Mked3Q3SS5TWJe4CMx71KGCsjl9jX/BgQ9OBHdSUTn15/vAfgggB96a1HCTI5ONj4aj3UctLkqWbXZ+N9Ustq28OT6ye6pjyakiXmXMtTV19b/DpyYD4Wky7kLwN+u5+ETAF9h5q/5bVYeuOjhW5/DUdBmQ1vJap3kdYUwbZInn5gXDCw3IEvRvwvggQHakj3FH71NdEmS7jpokyTpTUFbIWF5b/QqJvGbqt+GTMwLvCTmHofgbWipakRISlZtLqNJkJSszo63Nm580IQpcLOVvu3xfCTmXcpQAX+JuStU8AiQBG02TI4372xiPMZR0Faeh9dVttWVszaJX/65PmWowHgScxMqeEAms5PWpaMFrkpWJzPa2ErZJabQramctam3d7lxA+A/MS+ISW5ABY+emNeGl6kL3fpI7zoxL/CRmDeVoYZEBQ/A5OhksKDt7ed0VLJaF7SZrknbSA9sih/TVsdSuYcuQ5Wgglco/5Fc3uVESgprw8vD5KbyVRvpgdPi+0jMu2zcsPo8bGLepUwVUMEbqZ6RfQpvE7RJendXQVtBXZJet9tLnfh1q9HqxI8pMX/7eSK8HFaHCm7BkMKbCL02vFqX3nTpzCS+bW9fEEtiLipDtcRnmSqggvfCNOeylV6SpLsO2iRJuiRoMy1EcSt9t8S8YPCtjiNEBXeMtJd3FbTZ0DVoM+1EWnd9vG71WZ34bRs6DrnVca/EPFJUcM/4Hta7Klm1DdpspAfse/vyzySZmEdwiQxQwQeHZvNWyV2VrEqCtibqVpbVDd3rbppg09sX0odMzAtirjGXooJHRCxrw9uuh5vEb5qvdx3iNyXmBcG3Oo4cFTwCXAVt7xzbv2TVJNOQ0idxc8AEGI/gh0fAXlxbSUmCNqvH67k2vJCnbi5edxeUOvH7DvGjvjlgIoxHcOD0HyeQ8JJ5eEGIteGmvdqaAjgfvX2UNwdMkHEJXiUS4U0MuTa8PNytq2yr26DRpre3ld5nYh7zVsdVupapAmMXvIrphRxYetdBm22SbioYaSpntentbYf4pnYFuzlgoqjgbQzUy9usDZcgDdoKmgI3n9IDzXc4jeLmgAmjgtsScFjvqmTVFLTV7bhaJ37dTQxthvhthTVRJOYe8V2HDqjg/ekpfIi14cBmONVUvmoS31Vv3zSvd5WYF3RKzDNABXeNcB7vqqLNhrrevSx7mTrxfUoPRLbVcRciKVMFVPBhsOzlfa4Nr24V3DScN4lvIz3QfYg/1sTcNSp4ojjbhNHwfnYlPdCvt48uMU+syAXIRfBiSLTrcKeDAPhcG266KWHd5TOT9EC9+D6kt0nMC8aemJvIQ/CC6twnEeFDrA0H7KQH7Hp7F/P6oW4OmDN5CV4lcuFdB22mJL3upoR1Q3ef0q8e37K393xzQCDf3hvIXfAqpnQzhPSOF76Yeve26+Em8W2kB/wO8U3PH2Vi7pk+ZarA2AQ3EWkv76Jk1bZ3tZEeGKC395yYjwEVvEpg4fsGbcUbv67oxbZ3DSF92+O5SMzHggrexkDDetdBm6maq6nSzaZ3dTGvBwT7tnu6OWAMDFGmCqjg3XDYy9usDW+i3GvVnQR8Sg+47e01MXeDCu4Cz8N625JV08YGNtKvjjd/fcghvo/E3DsRlakCKrgfJMN6QZIuDdoKmgI3G+lXx5/+2pBD/DEm5j6IWvDqPIUi2nHFmrL0LT1817Xhpq1/gXrxY5IeOC2+Jub9aX0XEdE9AL4A4GcBLAEcMPNf+m6YiayEN2CzCWNV5sbe2yC+jfSA3yF+o/SxJOYJ1qEDsh78BMAfMvMLRHQbgMtEdImZf+C5ba2YksjcpK9Sl6TbSNzpeI+9fdulMpvEXNmkVXBmfg3Aa+vPbxDRFQDvBhBccBMp9vK2JatV2Wp7V9she2RD/L5bHUvIuUwVsJyDE9G9AB4A8Jzhe/sA9gFgF7c4aJobohbeVdBmK5rH3t7VEN/0eNEl5gkgFpyIbgXwJIBPMvP16veZ+QDAAQCcm9xRv3ogMFEM649mnYI201C07vKZtWiRDfHLP2OTmBfkEK71rUMHhIIT0RZWcn+Jmb/a+1kjI+pevgUb6QHLa+QOhvhdpe+SmCunkaToBODzAK4w82f8Nyk8IYU3JenliqyCpmG7T+lXj28xZO8xr29KzAtSlHuoMlVA1oNfBPBxAN8johfXX/sTZv4nf82KiyGENwVtTdfDTdID9eLXJc02Q/wh5/U2ZahKPZIU/VkAbrcSwWp+MUloKFzG6Ty+JWgz9VBNRTA+e3sX8/rV48uOTy5Ui6xMFQhcyVYNEVIVHujQyzcEbUVAVFf0UvfGtuntYx/iJyd3pERVqmpKDVOV3sWw3pQEN1W62fT2Lob4vqWvPrfKbU9UgpvIpZcvC18ne7nooq7wxaf0gF1v72JeDzRdI68vQx2URMtUgQQEr5KL8G2YKqxspAfshvih5vVAvfh1j6XISU7wKlkIX/4dGtpvIz1g19u7mNcD7ob45Z/xlZjnXqYKZCB4leTn8XXDQcGwvsCV9EC4Ib5eDnNDdoKbiLaXt9kJRnhTQ6C+Z0phXg/otW7ATZkqMBLBq4QWng+PzEGb7QaPFtIDw8/rge6X7vrInUMduitGKXiVEMN68WW0uuKJOvE9Sg+EuXSXE0OWqQKAu5tPZ4arIZIUPjyy++PbVE0dHlld6qHZ3CqAmhydWPWak9lJp8tdLi+RiX8/6evmsIrN5XtPBU+ZCEsjfSGVW3KiCSW35ATuumNRwcdEogUbErmlo4gxyQ2o4EoDMVwnlsotYWxyAyq44hDb9LpNXpW7P5qiR0btJbRM6ROcqdztaA+eOrZBW6Lz8CoqtwwVXGkkhnl4FZVbjgquOMV3FZnKbYcKHiFDVzulQsxySwuVhi6gUsFzYAQFL7HLLWFouQEVfJxYvtFCz8NV7u6o4IpzXM7DVe5+qOCRovNwldsFKnguZDYPV7ndoJVsY0VwZ9M+DLHpgsrdjvbgiojQQVsVlVuGCq4kh8otRwWPGOugLbN5uAmV2w4VfMxE+IZsQuW2R0M2RYzpFsdDPrcIiWgO68qBeOUGtAdXEkDl7o724A1420bH4vLU2DaAqKJy90N78NzIKGhTufujPfjYsSx4ie16uMrdTGsPTkSPE9FVInppiAaNgZTeIFGjcrciGaL/DYCHPLdDaUAXnhhQuUW0Cs7MzwB4Y4C2KK7IaB5uROUW4yxkI6J9InqeiJ6fc5ovxmhJ6c2rclvhTHBmPmDmC8x8YYvGe1lH8YjKbY1eJgtE6m+cQZHeHVXlPoUKngijXXgSYFtjIA+5AdllsicAfBPAfUT0ChF9wn+zlMGJ8Q2tcvemtdCFmR8ZoiGKsoHK7QStZAvI8vDIqi7dO6m9uVXuVnQOnhBa8FJC5RahgudMLkFbFZVbjAqupIXKbYUKHpjc32BOUbmtUcETY7TzcJW7Eyp47uQwD1e5O6OCK3GjcvdCBVfiReXujQoeAbZvvFHMw1VuJ2gl2xjIYR5uQOVuR3twJUlUbhkquJIcKrccFTwS9M0oQ+W2Q+fgiTLGO55I5FaxN9EeXEkClbsbKrgSPSp3d1TwiNA36WlU7n6o4AmTe8GLyt0fDdkSJ3fJm1C529EeXEkSlVuGCq4kh8otRwWPDH3zNqOvjx0quJIMKrc9KriSBCp3N1RwJXpU7u6o4BGib+h30NeiHyq4Ei0qd39UcCVKVG43aCVbpOgbXHGB9uCKkjEquKJkjAquKBmjgitKxqjgipIxIsGJ6CEi+g8i+hER/bHvRimK4oZWwYloCuBzAH4LwHsAPEJE7/HdMEVR+iPpwR8E8CNm/jEzHwP4MoDf9tssRVFcICl0eTeA/y79/xUAv1I9iIj2Aeyv/zu7NH/ipf7NG4SfBvB66EZYkFJ7U2orkFZ775McJBGcDF/jU19gPgBwAABE9DwzX5A0IDQptRVIq70ptRVIq71E9LzkOMkQ/RUA95T+fzeAV7s0SlGUYZEI/m0Av0hEP09E2wB+F8A/+G2WoiguaB2iM/MJEf0+gH8GMAXwODN/v+XHDlw0biBSaiuQVntTaiuQVntFbSXmU9NpRVEyQSvZFCVjVHBFyRingqdU0kpEjxPRVSKK/no9Ed1DRN8goitE9H0iejR0m5ogol0i+hYRfWfd3j8L3aY2iGhKRP9ORF8L3ZY2iOhlIvoeEb3YdrnM2Rx8XdL6nwB+A6tLa98G8Agz/8DJEziGiH4VwE0AX2Dm94ZuTxNEdBeAu5j5BSK6DcBlAL8T8WtLAM4y800i2gLwLIBHmfnfAjetFiL6FIALAM4x88Oh29MEEb0M4AIztxbluOzBkyppZeZnALwRuh0SmPk1Zn5h/fkNAFewqjCMEl5xc/3frfVHtGkuEd0N4CMA/jp0W1zjUnBTSWu0b8JUIaJ7ATwA4LmwLWlmPeR9EcBVAJeYOeb2fhbApwEsQzdECAP4FyK6vC4Rr8Wl4KKSVqU7RHQrgCcBfJKZr4duTxPMvGDm+7GqfHyQiKKcBhHRwwCuMvPl0G2x4CIz/xJWKzx/bz3dNOJScC1p9ch6LvskgC8x81dDt0cKM18D8DSAhwI3pY6LAD66ntd+GcAHiOiLYZvUDDO/uv73KoCnsJoeG3EpuJa0emIdWn0ewBVm/kzo9rRBRHcS0fn153sAPgjgh2FbZYaZH2Pmu5n5Xqzes//KzB8L3KxaiOjsOmgFEZ0F8CEAtVeCnAnOzCcAipLWKwC+IihpDQYRPQHgmwDuI6JXiOgTodvUwEUAH8eqd3lx/fHh0I1q4C4A3yCi72J14r/EzNFffkqEdwF4loi+A+BbAP6Rmb9ed7CWqipKxmglm6JkjAquKBmjgitKxqjgipIxKriiZIwKrigZo4IrSsb8P+nGjzaM/WqbAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plot(scalar_field_discrete, wireframe=True, scalarbar=False)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 49,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([1.52299797e-08, 2.08579001e-05, 1.10431945e-04, 2.26032941e-06,\n",
+       "       2.08579001e-05, 2.85655008e-02, 1.51239760e-01, 3.09558685e-03,\n",
+       "       1.10431945e-04, 1.51239760e-01, 8.00737403e-01, 1.63895538e-02,\n",
+       "       2.26032941e-06, 3.09558685e-03, 1.63895538e-02, 3.35462628e-04])"
+      ]
+     },
+     "execution_count": 49,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# values at integration points\n",
+    "scalar_field_discrete.compute_vertex_values()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 50,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.09835711945028527"
+      ]
+     },
+     "execution_count": 50,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# and at arbitrary point in domain\n",
+    "scalar_field_discrete(1.2343, 3.4872)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Similar for Vecvtor fields"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 51,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vector_field_expression = Expression( (\"cos(x[0])\", \"sin(x[0])\"), degree = 1)\n",
+    "vector_field_discrete = project(vector_field_expression, U)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 52,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.quiver.Quiver at 0x11cefa828>"
+      ]
+     },
+     "execution_count": 52,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAD4CAYAAAA0L6C7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAUxElEQVR4nO3de3SU9Z3H8fd3JpMEknAJCYhcDAKCYitqRC2gRV3XWttq1a7W2t1q19Nj3Vp7O+u2+0e7l7Pn7Km13d7E1t2qra5W6cVbxRZXRUW5qtwEEeVqglyScMll5rt/JIaIYCbDTH7jj8/rnJzJM/Pk93wg88lzzRNzd0QkHonQAUQkv1Rqkcio1CKRUalFIqNSi0SmpBCD1tTUeF1dXSGGjsbujlbe3LONCZUjSCUK8m2QiC1atGibu9ce7LWCvJvq6upYuHBhIYaOxp3rnuLHrz7G58edxQ2TLggdRz5gzOyNQ72mze9AVjVtAuCPmxbRlukInEZiolIHsrppMwA72nYzb+vywGkkJip1AC3t+9iw5+3u6Qc3LAiYRmKjUgewqmkT1aWVlCVSjBwwlIQZ61reCh1LIqFSBzC0tJL7Zt7E0QOHUlVSzs+m/T3HVBz0QKZIn2V19NvM1gPNQBrocPf6fAfxTDNYJWaW76GLzviqEQDUlFWxrqUBgKTp56vkR1/eSbPcfWohCg1A+0v4tnPJNP0L3voM7m0FWUwxGVZaxY7WFtKeCR2lX3nHIc/GSB4Eu+rB983F257v+Qykt8Geu/A9d4FV4KUzsPJZUPZRLFEdKmrB1JRVkcHZ0babmrKq0HEKzjtew5t/CGbYkB/imSZoXwEdK2DAxVF+j0PIttQOPG5mDtzm7rMPnMHMrgOuAxg7dmzvA7Ythj13vc8Mu6H1T3jHKmhfAxXXYsmaLON+MFw85jQ+OuIEBqUGhI5SUJ7ehLf8GPbOATKQrCPTeC6kN3TPYyXHQtlHg2UsBg89v4LbH1nAr751BUMqc39PZFvq6e6+2cyGA3PNbJW7P9Vzhq6izwaor6/v9c4LVnUjVF6//4nMDnzbRUAbpE7Fys+Bslmd3+xIja2oYWxFXD+oevL0Nnz3z2HPPUD7/hfSG6BkMgy4HEtNgZIpkJoULGexSCYSbGjcyXMr3uBj0ybnPE5WpXb3zV2PDWY2B5gGPPX+X/X+zMrByvcvo305NvhfoewsLDHkcIaWItC5dv4JdKyH5EhINwD7ul5NYoP/tbPQ0u0jU+pImPH0y+sKW2ozqwAS7t7c9fn5wPdyXuKhllN2Rr6HlIAsOQob/O/d0+4O3gKZtyDdCL4nYLriNLiinJPGH8385etpa+/g5fVbOXXi6D6Pk82aegQwp+tUUwnwG3d/rM9LkiOamYFVQaIKSiaEjlN0MhlnzaZGTp88liVrN3Ht9+/jhGNGFKbU7r4OOCmXoCKSnUTCePCZV7j/qWUALH/jLU445qjcxspnMBHJ3VcumcFR1ftPbeZ6HZZKLVIkKspL+c5V53VP53p1pUotUkQ+ckIdnzjzBAByvWBapRYpMl+/9GxqBg3U5rdILAZVlHPzleeS67pad7wTKUKzpk6gumpgTl+rNbVIkTpp/NE5fZ1KLRIZlVokMiq1SGRUapHIqNQikVGpRSKjUotERqUWiYxKLRIZlVokMiq1SGRUapHIqNQikVGpRSKjUotERqUWiYxKLRIZlVokMiq1SGRUapHIqNQikVGpRSKjUotERqUWiUzWpTazpJktMbOHChlIRA5PX9bUNwIrCxVERPIjq1Kb2Wjg48AvChtHRA5XtmvqW4FvAZlDzWBm15nZQjNb2NjYmJdwItJ3vZbazC4CGtx90fvN5+6z3b3e3etra2vzFlBE+iabNfV04JNmth64FzjHzO4uaCoRyVmvpXb3m919tLvXAVcAf3H3zxU8mYjkROepRSJT0peZ3f1J4MmCJBGRvNCaWiQyKrVIZFRqkcio1CKRUalFIqNSi0RGpRaJjEotEhmVOqDF218PHUEipFIH8nZrMz9fMzd0DImQSh3I6qbNLN2xnteat4aOIpFRqQNZ2bQJgAc3vBA4icRGpQ5kddNmAB7ZvIQ9Ha2B00hMVOpAVu3qLPXujlYe3/JS4DQSE5U6gJ1tu9m6b2f39IMbFuDuARNJTFTqAFY1bWZcRS0Dk6UcU1FDXcVwXmt5K3QsiYRKHcDYgTXcPf0fOHpANalECd876TNMqDoqdCyJRNGU2tNbcG8LHaNfHD1wKKlECcPKqni7tTl0HIlM0ZSadAPecDqZHV/G9zyAp7eFTlRww8oq2dG2m45MOnSUfuPueNuS0DGi1qd7lOWT730Q3/fnA57sgNa5eOtcwPDUh7Gyc6BsFpRMwsyCZC2UmrIqALa3tTC8fHDgNIXnrc/jLbdAciykPgQdr0PHCrx9OVZxLZYcETpicB2ZNCWJ5GGNEa7UHW9A27MHPNtz89uh/RXcyjErg0QtJIf1Z8SCu7JuOpeNPYNhXeWOlbe/hDf/ANrmdz6R3oi/dQqwb/9MpWfAEV7q+994jp+++jj3z7yJmvJBOY8TrNSJqpug6qbuaU9vxRvPBRsIZWdhZbM6HxO5/+OKXfRl7liLN98KrY+/+4VMM5TNwFInQMkUSJ0AiSO70NC55bY73cr8ba/yqdH1OY8TrNTvkd6KVf83pE7BrHhiSW48vQnfOwesFFLTINMImQbw3UA7NvBqrOwjoWMWlWk1EymxJPMbV8VRaiudGjqC5JElR2FV33zP855p6Sz4EXKmoy8qSso4tXocC7atpaV9H4u2r+PsESf0eZziOfotRwRLVGIl47DUpNBRikrGMzzX+CrHDTqavek2rn72xzy/bU1OY6nUIkUgYQmW7XyDu15/CoBNe7eT68kelVqkSFwzfhbjKod3Txu5tVqlFikSpYkS/vnES0nkWOZ3qNQiReTEIWO4om46oDW1SDS+NPE8Rg+sLtw+tZmVm9kLZrbMzJab2XdzW5SIZKM8Wcp3Tvw0luM6N5vz1K3AOe7eYmYp4Bkze9Tdn89piSLSq1Oqj2VwamBOX9trqb3zlhwtXZOprg/dpkOkwMbn+Dv2Wa3fzSxpZkuBBmCuuy84yDzXmdlCM1vY2NiYUxgROXxZldrd0+4+FRgNTDOzEw8yz2x3r3f3+tra2nznFJEs9WlP3N13Ak8CFxQkjYgctmyOftea2ZCuzwcA5wGrCh1MRHKTzdHvkcCvzCxJ5w+B+9z9ocLGEpFcZXP0+yXg5H7IIiJ5oCvKRCKjUotERqUWiYxKLRIZlVokMiq1SGRUapHIqNQikVGpRSKjUotERqUWiYxKLRIZlVokMiq1SGRUapHIqNQikVGpRSKjUotERqUWiYxKLRIZlVokMiq1SGRUapHIqNQikVGpRSKjUotERqUWiYxKLRIZlVokMiq1SGRUapHI9FpqMxtjZvPMbKWZLTezG/sjmIjkptc/Og90AF9398VmVgUsMrO57r6iwNlEJAe9rqndfYu7L+76vBlYCYwqdDARyU2f9qnNrA44GVhwkNeuM7OFZrawsbExP+lEpM+yLrWZVQIPAF9196YDX3f32e5e7+71tbW1+cwYJXfH9/05dAyJUFalNrMUnYX+tbs/WNhIR4jMW3jLf+HuoZNIZLI5+m3AL4GV7n5L4SMdIdqXQ8cKaH8pdBKJTDZr6unA1cA5Zra06+PCAueKnre/0vm4957ASSQ2vZ7ScvdnAOuHLEeWjuWdj3sfxqtuxhKDw+aRaOiKslDa3znN3wp75wSNInFRqQPwdANkGvZP77lHB8wkb1TqENqXQ+npYIOgZBKUnwcda0Knkkio1CGkPowNvROSo4AMiapvYqnjQqeSSBRNqV9vaWBba3PoGP3CksMwM0jUQFpX30l+ZfMLHf0i485F8/6DyYNHMbN2MjOHT2Zi1cjON3+sksOhbSfubZiVhk7TL9zboPVprPzc0FGiFazU/7v+WR7dsvRdzyUtwYpdG1mxayO3rX2C4eWDmVE7iZnDj6e++ljKkqlAaQskUdP5mNkOyaPCZikw9zTsewhv+RGUngGlJ3eeAWhfjneswCq/iZWMDh0zuJ07djNkaMVhjRGs1G3ewZ6O1nc9d+Dx333pNvam29jT0Up7Jh1dqa3yeqj8Cp1X4cbJ3aH1Cbzl1v0HA/c24Hvv7zFXCgZcBkd4qe+9cz533DaPO397A0eNHJLzOMFKffW4s7h63Fnd02/s3sbfPP0DxlXUMmP48cyoncSHhoylJJEMFbHgzAaEjlBQ3vos3nLLey+FtSQM+CxWMgVSU6BkwhGz+/F+jp04AndY8OxaPnVpfc7jFM0+dVumnd+e9TVGDxwWOorkgac3QcerUDoNknWd5+XTDZBpBG/GUtOwAbrauKepp9RRXp5iwbNr4ij1xKqRoSNIHllyFFT83UGvL/bMHvAj40xHX5SWlTC1vo5FL6yjsaGJJQtf5/wLT+rzOEVzSkuOHJYYiCVHhI5RVDIZ5/FHllE7fBDtbWm+9Le3s2r55pzGUqlFikAiYTQ37eOPDy4CoGnX3px/jUqlFikSF19+GsdP2X/7v1yv0VCpRYpEMpng69/+BKlU5xmfXK+7UqlFisgxdTV87pqZAFiO298qtUiR+cxVZzJ+4gjtU4vEoqQkyTf+af9meJ+/Ps95RCQPJkw6ispB5Tl9rdbUIkUq1+u/VWqRyKjUIpFRqUUio1KLREalFomMSi0SGZVaJDIqtUhkVGqRyKjUIpHJ5o/O32FmDWb2Sn8EEpHDk82a+n+ACwqcQ0TypNdSu/tTwPZ+yCIieZC3fWozu87MFprZwsZG/dE3kVDyVmp3n+3u9e5eX1tbm69hRaSPdPRbJDIqtUhksjmldQ/wHDDJzDaa2bWFjyUiuer1HmXufmV/BBGR/NDmt0hkVGqRyKjUIpFRqUUio1KLREalFomMSi0SGZVaJDIqtUhkVGqRyKjUIpFRqUUio1KLREalFomMSi0SGZVaJDIqtUhkVGqRyKjUIpFRqUUio1KLREalDsjdQ0eQCKnUgbg7Tzz2cugYEiGVOpCtW3Zyz6/ma20teadSB7Jm1VY2vPk2SxetDx1FIqNSB7Jm9RYAHvrd4sBJJDYqdSBrVm8FYP7/rebtbc2B00hMVOoA3L17TZ1OZ3jsoWWBE0lMVOoAGrbuomnX3u7pR36/mHQ6EzCRxESlDmDNq1s59/wTGTxkIBMnj+Tyq85k45tvh44lkej1T9lK/p067VhmnD2Z67/wC/btbePiy04LHUkiktWa2swuMLPVZrbWzP6xEEFW7NrIyl2bjojztgMGlAJQPaySHdt3B04jsel1TW1mSeAnwF8BG4EXzewP7r4in0GGpCq49OnvU11ayYzaScwYfjynDRtPeTKVz8UUlaHDKmlp3kdraztlZfH+O3tq3tHC/DkvcME154SOEq1sNr+nAWvdfR2Amd0LfAo4rFL/cu1feGDDgvc839jaxJyNLzJn44uUJVKcNmw8M4dPZtaIKQwprTicRRad2toqqodV0ty0j7LauEu9d/c+fvejR7nvP3/P2ZefyYSTx7Fm8TrWLnmdtUte5xt3XM+YSaNCxwzK3VnTvJWJVUdhZjmPk02pRwEbekxvBE4/cCYzuw64DmDs2LG9DlpdVsmxlSPe9dyyHetJ+/6jwKMGDmV85QjGV46gKjUgi6gfLJ//4tl8/otnh45RUO1t7Tw8+wl+828PsOOtXQA8fPsTPHz7E93z1IyqZvuWnUd8qW9b+wR3vDaP+2fcxDGVtTmPk02pD/Yj4z07vu4+G5gNUF9f3+uO8SVjpnHJmGnd06ubNnPNcz/j9GETmFE7mRnDJzNqYHUW8aQYpdNp/nz309z13fvYur7xXa+VV5Rx5c2fZuKpxzLh5HEMHT44UMriUl99LHe8No9nGlcVvNQbgTE9pkcDm3Ne4iFUlJTxp3O/TWVJeb6HlgC2b9lJeUUZn/7qRby9eTvbt+7sfNzS+Vg5pILT/npq6JhFZerQOipKynimcRVXjZuJu+e0GZ5NqV8EJprZOGATcAXw2T4vqRejBw7L95ASUO3oYdReduYhX29va+/HNB8MJYkkZ9Ycx7y3lrO2eSsL336NK+qm93mcXk9puXsHcAPwJ2AlcJ+7L+/zkkR6SJXGfWCwrzKe4devP00qUULaM3zhuZ/y5p5tOY2V1cUn7v4I8EhOSxCRXiUswZDSCh7dvASA1kwHdtDDWVmMlc9gIpK7C48+mTNrjjvscVRqkSJhZtw85WIqkmWd01pTi3zwHTVgCDdMugCAXK8/UalFiswlY07j1OpxOa+prRC/QGFmjcAbvcxWA+R2eK9/KF/uijkbFHe+bLMd4+4HvUKlIKXOhpktdPf6IAvPgvLlrpizQXHny0c2bX6LREalFolMyFLPDrjsbChf7oo5GxR3vsPOFmyfWkQKQ5vfIpFRqUUiE6TU/XEjw1yZ2R1m1mBmr4TOciAzG2Nm88xspZktN7MbQ2fqyczKzewFM1vWle+7oTMdyMySZrbEzB4KneVAZrbezF42s6VmtjDncfp7n7rrRoav0uNGhsCV+b6RYa7M7CygBbjT3U8MnacnMxsJjHT3xWZWBSwCLi6i/zsDKty9xcxSwDPAje7+fOBo3czsa0A9MMjdLwqdpyczWw/Uu/thXRgTYk3dfSNDd28D3rmRYVFw96eA7aFzHIy7b3H3xV2fN9P5++1Fc2Mv79TSNZnq+iiaI7FmNhr4OPCL0FkKKUSpD3Yjw6J5Y35QmFkdcDLw3luyBtS1ebsUaADmunsx5bsV+BZQrH/jyIHHzWxR1408cxKi1FndyFAOzcwqgQeAr7p7U+g8Pbl72t2n0nkvu2lmVhS7MGZ2EdDg7otCZ3kf0939FOBjwJe7dgX7LESp++VGhrHq2ld9APi1uz8YOs+huPtO4EnggsBR3jEd+GTXfuu9wDlmdnfYSO/m7pu7HhuAOXTuqvZZiFJ338jQzErpvJHhHwLk+MDpOhD1S2Clu98SOs+BzKzWzIZ0fT4AOA9YFTZVJ3e/2d1Hu3sdne+5v7j75wLH6mZmFV0HPzGzCuB8IKczMP1e6mK/kaGZ3QM8B0wys41mdm3oTD1MB66mcy2ztOvjwtChehgJzDOzl+j84T3X3Yvu1FGRGgE8Y2bLgBeAh939sVwG0mWiIpHRFWUikVGpRSKjUotERqUWiYxKLRIZlVokMiq1SGT+H0o+CyZw5sL2AAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plot(vector_field_discrete, wireframe=True)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Solving PDEs "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Consider a diffusion problem, such as the transport of heat in an isotropic homogeneous medium:\n",
+    "\n",
+    "$$\\frac{\\partial \\Theta}{\\partial t} = - \\alpha\\, \\nabla^2 \\Theta$$\n",
+    "with \n",
+    "- initial condition $\\Theta(x, t=0) = \\Theta_0$\n",
+    "- isolated boundary condition $\\frac{\\partial \\Theta}{\\partial n} = 0$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To solve this equation, we start by writing the weak form:\n",
+    "    \n",
+    "$$\\int_\\Omega \\left( \\frac{\\partial \\Theta}{\\partial t}\\right ) dx= - \\alpha\\, \\int_\\Omega \\left( \\nabla^2 \\Theta \\right) dx$$ \n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Thermal equilibrium"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We introduce weighting functions $\\nu$:\n",
+    "\n",
+    "$$\\int_\\Omega \\left( \\frac{\\partial \\Theta}{\\partial t}\\right )\\, \\nu\\, dx= - \\alpha\\, \\int_\\Omega \\left( \\nabla^2 \\Theta \\right)\\, \\nu\\, dx $$ "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Using chain rule and divergence theorem, the second term of eq. (8) can be written as:\n",
+    "\n",
+    "\\begin{align}\n",
+    "\\int_\\Omega \\left( \\nabla^2 \\Theta \\right)\\, \\nu\\, dx \n",
+    "&= \\int_\\Omega \\left( \\nabla^2 \\Theta \\right)\\, \\nu\\, dx \\\\\n",
+    "&= \\int_\\Omega \\left( \\nabla \\cdot \\nabla\\, \\Theta \\right)\\, \\nu\\, dx \\\\\n",
+    "&= \\int_\\Omega \\left( \\nabla \\Theta \\right) \\cdot \\, \\left( \\nabla \\nu \\right) \\nu\\, dx - \n",
+    " \\int_\\Omega \\left( \\nabla \\Theta \\right) \\cdot \\nu\\, dx \n",
+    "\\end{align}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "\\begin{align}\n",
+    "\\int_\\Omega \\left( \\nabla^2 \\Theta \\right)\\, \\nu\\, dx \n",
+    "&= \\int_\\Omega \\left( \\nabla^2 \\Theta \\right)\\, \\nu\\, dx \\\\\n",
+    "&= \\int_\\Omega \\left( \\nabla \\cdot \\nabla\\, \\Theta \\right)\\, \\nu\\, dx \\\\\n",
+    "&= \\int_\\Omega \\left( \\nabla \\Theta \\right) \\cdot \\, \\left( \\nabla \\nu \\right) \\nu\\, dx - \n",
+    " \\int_\\Omega \\left( \\nabla \\Theta \\right) \\cdot \\nu\\, dx \n",
+    "\\end{align}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\\int_\\Omega \\left( \\frac{\\partial \\Theta}{\\partial t}\\right )\\, \\nu\\, dV + \\alpha \\int_\\Omega \\left( \\nabla\\nu \\right)\\,\\cdot\\, \\left( \\nabla \\Theta\\right) \\, dV\n",
+    "=  \\int_\\Omega f\\, \\nu\\,  dV\n",
+    "+\\int_{\\partial\\Omega} \\nu\\, \\Theta dS$$ \n",
     "\n"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 68,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n"
+     ]
+    },
+    {
+     "ename": "RuntimeError",
+     "evalue": "\n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n***     fenics-support@googlegroups.com\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a *minimal* running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error:   Unable to define linear variational problem a(u, v) == L(v) for all v.\n*** Reason:  Expecting the left-hand side to be a bilinear form (not rank 1).\n*** Where:   This error was encountered inside LinearVariationalProblem.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2019.1.0\n*** Git changeset:  \n*** -------------------------------------------------------------------------\n",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
+      "\u001b[0;32m<ipython-input-68-623aa4b253e8>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     19\u001b[0m \u001b[0mL\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mT_0\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mdx\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     20\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 21\u001b[0;31m \u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mL\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\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~/anaconda2/envs/fenicsenv/lib/python3.7/site-packages/dolfin/fem/solving.py\u001b[0m in \u001b[0;36msolve\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    218\u001b[0m     \u001b[0;31m# tolerance)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    219\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mufl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclasses\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mEquation\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--> 220\u001b[0;31m         \u001b[0m_solve_varproblem\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\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    221\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    222\u001b[0m     \u001b[0;31m# Default case, just call the wrapped C++ solve function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m~/anaconda2/envs/fenicsenv/lib/python3.7/site-packages/dolfin/fem/solving.py\u001b[0m in \u001b[0;36m_solve_varproblem\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    240\u001b[0m         \u001b[0;31m# Create problem\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    241\u001b[0m         problem = LinearVariationalProblem(eq.lhs, eq.rhs, u, bcs,\n\u001b[0;32m--> 242\u001b[0;31m                                            form_compiler_parameters=form_compiler_parameters)\n\u001b[0m\u001b[1;32m    243\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    244\u001b[0m         \u001b[0;31m# Create solver and call solve\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m~/anaconda2/envs/fenicsenv/lib/python3.7/site-packages/dolfin/fem/problem.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, a, L, u, bcs, form_compiler_parameters)\u001b[0m\n\u001b[1;32m     57\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     58\u001b[0m         \u001b[0;31m# Initialize C++ base class\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 59\u001b[0;31m         \u001b[0mcpp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mLinearVariationalProblem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__init__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mL\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cpp_object\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbcs\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     60\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     61\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;31mRuntimeError\u001b[0m: \n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n***     fenics-support@googlegroups.com\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a *minimal* running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error:   Unable to define linear variational problem a(u, v) == L(v) for all v.\n*** Reason:  Expecting the left-hand side to be a bilinear form (not rank 1).\n*** Where:   This error was encountered inside LinearVariationalProblem.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2019.1.0\n*** Git changeset:  \n*** -------------------------------------------------------------------------\n"
+     ]
+    }
+   ],
+   "source": [
+    "T_0_expression = Expression(\"exp( - pow(x[0]-x0, 2) / pow(sigma_x, 2) - pow(x[1]-y0, 2) / pow(sigma_y, 2) )\", \n",
+    "                                    x0=3, y0=3, sigma_x=1, sigma_y=1, \n",
+    "                                    degree = 1)\n",
+    "T_0 = project(T_0_expression, V)\n",
+    "\n",
+    "u_D = Expression(\"1 + x[0]*x[0] + 2*x[1]*x[1]\", degree=2)\n",
+    "\n",
+    "def boundary(x, on_boundary):\n",
+    "    return on_boundary\n",
+    "\n",
+    "\n",
+    "bc = DirichletBC(V, u_D, boundary)\n",
+    "\n",
+    "T = Function(V)\n",
+    "v = TestFunction(V)\n",
+    "alpha = Constant(1)\n",
+    "\n",
+    "a = inner(grad(v), grad(T)) * dx\n",
+    "L = v*T_0*dx\n",
+    "\n",
+    "solve(a == L, T, [])\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 67,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.0"
+      ]
+     },
+     "execution_count": 67,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD8CAYAAABaQGkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAJSklEQVR4nO3b32tehR3H8c9nTUekOnNRJ52RdRdDJsLsFupFYbDipFPRXVrQKyE326hsIPPSf0C82cWClll0iqCF0TG3ghYp+CuprbPGDXGFlQqhuqC9iKP1s4s8BdEkz0mfc/Lk+fJ+QeiT9lg/SN89zznP0UkEoKZvDHsAgO4QOFAYgQOFEThQGIEDhRE4UNhYk4Nsn5H0maRLki4mmepyFIB2NAq856dJzne2BEDreIsOFOYmT7LZ/rek/0qKpD8kmVnhmGlJ05I0dtXYj6/97rUtTwVw2cfvf3w+yXX9jmsa+HeSnLP9bUlHJf06yaurHb/9B9tz11P3rmswgOYO3XZwrsm9sEZv0ZOc6/24IOmwpN2DzQOwEfoGbnub7Wsuv5Z0h6R3ux4GYHBN7qJfL+mw7cvH/ynJS52uAtCKvoEn+VDSDzdgC4CW8TEZUBiBA4UROFAYgQOFEThQGIEDhRE4UBiBA4UROFAYgQOFEThQGIEDhRE4UBiBA4UROFAYgQOFEThQGIEDhRE4UBiBA4UROFAYgQOFEThQGIEDhRE4UBiBA4UROFAYgQOFEThQGIEDhRE4UBiBA4UROFAYgQOFEThQWOPAbW+x/bbtI10OAtCe9ZzBD0ia72oIgPY1Ctz2pKS7JD3R7RwAbWp6Bn9c0sOSvljtANvTtmdtzy4tLrUyDsBg+gZu+25JC0nm1jouyUySqSRT4xPjrQ0EcOWanMH3SLrH9hlJz0naa/vpTlcBaEXfwJM8kmQyyU5J90l6Ocn9nS8DMDA+BwcKG1vPwUmOSTrWyRIAreMMDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYX0Dtz1u+03bp2yftv3oRgwDMLixBsd8Lmlvkgu2t0o6bvuvSV7veBuAAfUNPEkkXeh9u7X3lS5HAWhHo2tw21tsn5S0IOlokjdWOGba9qzt2aXFpbZ3ArgCjQJPcinJrZImJe22fcsKx8wkmUoyNT4x3vZOAFdgXXfRkyxKOiZpXydrALSqyV3062xP9F5fJel2Se93PQzA4JrcRd8h6SnbW7T8F8LzSY50OwtAG5rcRX9H0q4N2AKgZTzJBhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFNY3cNs32n7F9rzt07YPbMQwAIMba3DMRUm/TXLC9jWS5mwfTfJex9sADKjvGTzJR0lO9F5/Jmle0g1dDwMwuHVdg9veKWmXpDdW+LVp27O2Z5cWl9pZB2AgjQO3fbWkFyQ9lOTTr/56kpkkU0mmxifG29wI4Ao1Ctz2Vi3H/UySF7udBKAtTe6iW9KTkuaTPNb9JABtaXIG3yPpAUl7bZ/sfd3Z8S4ALej7MVmS45K8AVsAtIwn2YDCCBwojMCBwggcKIzAgcIIHCiMwIHCCBwojMCBwggcKIzAgcIIHCiMwIHCCBwojMCBwggcKIzAgcIIHCiMwIHCCBwojMCBwggcKIzAgcIIHCiMwIHCCBwojMCBwggcKIzAgcIIHCiMwIHCCBwojMCBwggcKKxv4LYP2l6w/e5GDALQniZn8D9K2tfxDgAd6Bt4klclfbIBWwC0rLVrcNvTtmdtzy4tLrX12wIYQGuBJ5lJMpVkanxivK3fFsAAuIsOFEbgQGFNPiZ7VtJrkm6yfdb2g93PAtCGsX4HJNm/EUMAtI+36EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGGNAre9z/Y/bX9g+3ddjwLQjr6B294i6feSfi7pZkn7bd/c9TAAg2tyBt8t6YMkHyb5n6TnJN3b7SwAbRhrcMwNkv7zpe/PSrrtqwfZnpY03fv280O3HXx38HkbYruk88MesQ6jtHeUtkqjtfemJgc1Cdwr/Fy+9hPJjKQZSbI9m2SqyYBhG6Wt0mjtHaWt0mjttT3b5Lgmb9HPSrrxS99PSjp3JaMAbKwmgb8l6fu2v2f7m5Luk/TnbmcBaEPft+hJLtr+laS/Sdoi6WCS033+sZk2xm2QUdoqjdbeUdoqjdbeRludfO1yGkARPMkGFEbgQGGtBj5Kj7TaPmh7wfam/7ze9o22X7E9b/u07QPD3rQW2+O237R9qrf30WFv6sf2Fttv2z4y7C392D5j+x+2T/b7uKy1a/DeI63/kvQzLX+09pak/Unea+Vf0DLbP5F0QdKhJLcMe89abO+QtCPJCdvXSJqT9ItN/N/WkrYluWB7q6Tjkg4keX3I01Zl+zeSpiR9K8ndw96zFttnJE0l6ftQTptn8JF6pDXJq5I+GfaOJpJ8lORE7/Vnkua1/IThppRlF3rfbu19bdq7ubYnJd0l6Ylhb2lbm4Gv9Ejrpv1DOKps75S0S9Ibw12ytt5b3pOSFiQdTbKZ9z4u6WFJXwx7SEOR9Hfbc71HxFfVZuCNHmnFlbN9taQXJD2U5NNh71lLkktJbtXyk4+7bW/KyyDbd0taSDI37C3rsCfJj7T8f3j+sne5uaI2A+eR1g71rmVfkPRMkheHvaepJIuSjknaN+Qpq9kj6Z7ede1zkvbafnq4k9aW5FzvxwVJh7V8ebyiNgPnkdaO9G5aPSlpPsljw97Tj+3rbE/0Xl8l6XZJ7w931cqSPJJkMslOLf+ZfTnJ/UOetSrb23o3WmV7m6Q7JK36SVBrgSe5KOnyI63zkp5v8Ejr0Nh+VtJrkm6yfdb2g8PetIY9kh7Q8tnlZO/rzmGPWsMOSa/YfkfLf/EfTbLpP34aEddLOm77lKQ3Jf0lyUurHcyjqkBhPMkGFEbgQGEEDhRG4EBhBA4URuBAYQQOFPZ/+EQxa8CAZuwAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plot(T)\n",
+    "T(2,2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "LinearVariationalProblem()"
+   ]
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 2
+    "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython2",
-   "version": "2.7.6"
-  },
-  "kernelspec": {
-   "name": "python3",
-   "language": "python",
-   "display_name": "Python 3"
+   "pygments_lexer": "ipython3",
+   "version": "3.7.3"
   },
   "pycharm": {
    "stem_cell": {
     "cell_type": "raw",
-    "source": [],
     "metadata": {
      "collapsed": false
-    }
+    },
+    "source": []
    }
   }
  },
  "nbformat": 4,
- "nbformat_minor": 0
-}
\ No newline at end of file
+ "nbformat_minor": 4
+}
diff --git a/FEMII/project_growth_model.ipynb b/FEMII/project_growth_model.ipynb
index 72f05c2..570e4df 100644
--- a/FEMII/project_growth_model.ipynb
+++ b/FEMII/project_growth_model.ipynb
@@ -1,52 +1,51 @@
 {
  "cells": [
   {
    "cell_type": "markdown",
    "metadata": {
-    "collapsed": true,
     "pycharm": {
      "name": "#%% md\n"
     }
    },
    "source": [
     "# Finite Growth \n",
     "using multiplicative Decomposition of Deformation Gradient\n",
     "\n",
     "# Preparation:​\n",
     "- theory overview​\n",
     "- simple example formulation for F_g​\n",
     "\n"
    ]
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 2
+    "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython2",
-   "version": "2.7.6"
-  },
-  "kernelspec": {
-   "name": "python3",
-   "language": "python",
-   "display_name": "Python 3"
+   "pygments_lexer": "ipython3",
+   "version": "3.7.3"
   },
   "pycharm": {
    "stem_cell": {
     "cell_type": "raw",
-    "source": [],
     "metadata": {
      "collapsed": false
-    }
+    },
+    "source": []
    }
   }
  },
  "nbformat": 4,
- "nbformat_minor": 0
-}
\ No newline at end of file
+ "nbformat_minor": 4
+}
diff --git a/FEMII/project_hyperelastic_model.ipynb b/FEMII/project_hyperelastic_model_abaqus.ipynb
similarity index 78%
rename from FEMII/project_hyperelastic_model.ipynb
rename to FEMII/project_hyperelastic_model_abaqus.ipynb
index af7f5db..b5b7e3d 100644
--- a/FEMII/project_hyperelastic_model.ipynb
+++ b/FEMII/project_hyperelastic_model_abaqus.ipynb
@@ -1,56 +1,57 @@
 {
  "cells": [
   {
    "cell_type": "markdown",
    "metadata": {
-    "collapsed": true,
     "pycharm": {
      "name": "#%% md\n"
     }
    },
    "source": [
     "\n",
     "# Hyperelastic Material Model for Brain Tissue\n",
     "\n",
     "## Tasks:​\n",
-    "- implement linear elastic, Neo-Hookean & Ogden model eq 8 in [1]\n",
+    "- implement Abaqus Ogden model \n",
     "- reproduce stress-shear/stretchrelations from Fig 24 in in simple geometry, using parameter values from Table 1 in [1]​\n",
+    "- optimize using FEM model -- better fit possible?\n",
+    "- linearize and compare linearized and non linear models\n",
     "- calibrate linear elastic model to best agree with experimental stress/strain relation in compression​\n",
     "- enable thermal expansion in linear elastic and hyperelastic material; compare temporal stress evolution in test geometry​\n",
     "- apply both models to few seeding scenarios in brain anatomy and compare resulting tumor shapes and stress distribution.​\n",
     "\n",
     "[1] https://link.springer.com/content/pdf/10.1007%2Fs11831-019-09352-w.pdf ​\n"
    ]
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 2
+    "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython2",
-   "version": "2.7.6"
-  },
-  "kernelspec": {
-   "name": "python3",
-   "language": "python",
-   "display_name": "Python 3"
+   "pygments_lexer": "ipython3",
+   "version": "3.7.3"
   },
   "pycharm": {
    "stem_cell": {
     "cell_type": "raw",
-    "source": [],
     "metadata": {
      "collapsed": false
-    }
+    },
+    "source": []
    }
   }
  },
  "nbformat": 4,
- "nbformat_minor": 0
-}
\ No newline at end of file
+ "nbformat_minor": 4
+}
diff --git a/tests/linear_elastics.ipynb b/tests/linear_elastics.ipynb
index acb9dac..dd6361f 100644
--- a/tests/linear_elastics.ipynb
+++ b/tests/linear_elastics.ipynb
@@ -1,407 +1,428 @@
 {
  "cells": [
   {
    "cell_type": "markdown",
    "metadata": {
-    "collapsed": true,
     "pycharm": {
      "is_executing": false,
      "name": "#%% md\n"
     }
    },
    "source": [
     "From https://comet-fenics.readthedocs.io/en/latest/demo/elasticity/2D_elasticity.py.html"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 1,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    },
+    "pycharm": {
+     "is_executing": false,
+     "name": "#%%\n"
+    }
+   },
    "outputs": [],
    "source": [
     "from __future__ import print_function\n",
     "from dolfin import *\n",
     "\n",
     "L = 25.\n",
     "H = 1.\n",
     "Nx = 250\n",
     "Ny = 10\n",
     "mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, \"crossed\")"
-   ],
-   "metadata": {
-    "collapsed": false,
-    "pycharm": {
-     "name": "#%%\n",
-     "is_executing": false
-    }
-   }
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {
+    "pycharm": {
+     "is_executing": false,
+     "name": "#%% md\n"
+    }
+   },
    "source": [
     "Constitutive relation\n",
     "---------------------\n",
     "\n",
     "\n",
     "We now define the material parameters which are here given in terms of a Young's\n",
     "modulus :math:`E` and a Poisson coefficient :math:`\\nu`. In the following, we will\n",
     "need to define the constitutive relation between the stress tensor :math:`\\boldsymbol{\\sigma}`\n",
     "and the strain tensor :math:`\\boldsymbol{\\varepsilon}`. Let us recall\n",
     "that the general expression of the linear elastic isotropic constitutive relation\n",
     "for a 3D medium is given by:\n",
     "\n",
     "$\\boldsymbol{\\sigma} = \\lambda \\text{tr}(\\boldsymbol{\\varepsilon})\\mathbf{1} + 2\\mu\\boldsymbol{\\varepsilon}$\n",
     "\n",
     "for a natural (no prestress) initial state where the Lamé coefficients are given by:\n",
     "\n",
     "$\\lambda = \\dfrac{E\\nu}{(1+\\nu)(1-2\\nu)}, \\quad \\mu = \\dfrac{E}{2(1+\\nu)}$\n",
     "\n",
     "In this demo, we consider a 2D model either in plane strain or in plane stress conditions.\n",
     "Irrespective of this choice, we will work only with a 2D displacement vector $\\boldsymbol{u}=(u_x,u_y)$\n",
     "and will subsequently define the strain operator $\\varepsilon$ as follows:"
-   ],
-   "metadata": {
-    "collapsed": false,
-    "pycharm": {
-     "name": "#%% md\n",
-     "is_executing": false
-    }
-   }
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": 2,
-   "outputs": [],
-   "source": [
-    "def eps(v):\n",
-    "    return sym(grad(v))"
-   ],
    "metadata": {
     "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    },
     "pycharm": {
-     "name": "#%%\n",
-     "is_executing": false
+     "is_executing": false,
+     "name": "#%%\n"
     }
-   }
+   },
+   "outputs": [],
+   "source": [
+    "def eps(v):\n",
+    "    return sym(grad(v))"
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   },
    "source": [
     "which computes the 2x2 plane components of the symmetrized gradient tensor of\n",
     "any 2D vectorial field. In the plane strain case, the full 3D strain tensor is defined as follows:\n",
     "\n",
     "$\\boldsymbol{\\varepsilon} = \\begin{bmatrix} \\varepsilon_{xx} & \\varepsilon_{xy} & 0\\\\\n",
     "\\varepsilon_{xy} & \\varepsilon_{yy} & 0 \\\\ 0 & 0 & 0\\end{bmatrix}$\n",
     "\n",
     "so that the 2x2 plane part of the stress tensor is defined in the same way as for the 3D case\n",
     "(the out-of-plane stress component being given by $\\sigma_{zz}=\\lambda(\\varepsilon_{xx}+\\varepsilon_{yy})$.\n",
     "\n",
     "In the plane stress case, an out-of-plane strain component $\\varepsilon_{zz}$\n",
     "must be considered so that $\\sigma_{zz}=0$. \n",
     "Using this condition in the 3D constitutive relation, one has $\\varepsilon_{zz}=-\\dfrac{\\lambda}{\\lambda+2\\mu}(\\varepsilon_{xx}+\\varepsilon_{yy})$.\n",
     "Injecting into , we have for the 2D plane stress relation:\n",
     "\n",
     "$\\boldsymbol{\\sigma} = \\lambda^* \\text{tr}(\\boldsymbol{\\varepsilon})\\mathbf{1} + 2\\mu\\boldsymbol{\\varepsilon}$\n",
     "where $\\boldsymbol{\\sigma}$, $\\boldsymbol{\\varepsilon}$, $\\mathbf{1}$ are 2D tensors and with\n",
     "$\\lambda^* = \\dfrac{2\\lambda\\mu}{\\lambda+2\\mu}$. \n",
     "Hence, the 2D constitutive relation\n",
     "is identical to the plane strain case by changing only the value of the Lamé coefficient $\\lambda$.\n",
     "We can then have::"
-   ],
-   "metadata": {
-    "collapsed": false,
-    "pycharm": {
-     "name": "#%% md\n"
-    }
-   }
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": 12,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    },
+    "pycharm": {
+     "is_executing": false,
+     "name": "#%%\n"
+    }
+   },
    "outputs": [],
    "source": [
     "E = Constant(1e5)\n",
     "nu = Constant(0.3)\n",
     "model = \"plane_stress2\"\n",
     "\n",
     "\n",
     "mu = E/2/(1+nu)\n",
     "lmbda = E*nu/(1+nu)/(1-2*nu)\n",
     "if model == \"plane_stress\":\n",
     "    lmbda = 2*mu*lmbda/(lmbda+2*mu)\n",
     "\n",
     "def sigma(v):\n",
     "    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)"
-   ],
-   "metadata": {
-    "collapsed": false,
-    "pycharm": {
-     "name": "#%%\n",
-     "is_executing": false
-    }
-   }
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   },
    "source": [
     "Note that we used the variable name ``lmbda`` to avoid any confusion with the\n",
     "lambda functions of Python\n",
     "\n",
     "We also used an intrinsic formulation of the constitutive relation. Example of\n",
     "constitutive relation implemented with a matrix/vector engineering notation\n",
     "will be provided in the :ref:`OrthotropicElasticity` example.\n",
     "\n",
     "Variational formulation\n",
     "-----------------------\n",
     "\n",
     "For this example, we consider a continuous polynomial interpolation of degree 2\n",
     "and a uniformly distributed loading :math:`\\boldsymbol{f}=(0,-f)` corresponding\n",
     "to the beam self-weight. The continuum mechanics variational formulation (obtained\n",
     "from the virtual work principle) is given by:\n",
     "\n",
     "\\text{Find } \\boldsymbol{u}\\in V \\text{ s.t. } \\int_{\\Omega}\n",
     "\\boldsymbol{\\sigma}(\\boldsymbol{u}):\\boldsymbol{\\varepsilon}(\\boldsymbol{v}) d\\Omega \n",
     "= \\int_{\\Omega} \\boldsymbol{f}\\cdot\\boldsymbol{v}  d\\Omega \\quad \\forall\\boldsymbol{v} \\in V\n",
     "\n",
     "which translates into the following FEniCS code::"
-   ],
-   "metadata": {
-    "collapsed": false,
-    "pycharm": {
-     "name": "#%% md\n"
-    }
-   }
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": 13,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    },
+    "pycharm": {
+     "is_executing": false,
+     "name": "#%%\n"
+    }
+   },
    "outputs": [],
    "source": [
     "rho_g = 1e-3\n",
     "f = Constant((0,-rho_g))\n",
     "\n",
     "V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)\n",
     "du = TrialFunction(V)\n",
     "u_ = TestFunction(V)\n",
     "a = inner(sigma(du), eps(u_))*dx\n",
     "l = inner(f, u_)*dx"
-   ],
-   "metadata": {
-    "collapsed": false,
-    "pycharm": {
-     "name": "#%%\n",
-     "is_executing": false
-    }
-   }
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   },
    "source": [
     "Resolution\n",
     "----------\n",
     "\n",
     "Fixed displacements are imposed on the left part of the beam, the ``solve``\n",
     "function is then called and solution is plotted by deforming the mesh::"
-   ],
-   "metadata": {
-    "collapsed": false,
-    "pycharm": {
-     "name": "#%% md\n"
-    }
-   }
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": 11,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    },
+    "pycharm": {
+     "is_executing": false,
+     "name": "#%%\n"
+    }
+   },
    "outputs": [
     {
      "data": {
-      "text/plain": "<matplotlib.collections.PolyCollection at 0x7f410cde6a90>"
+      "text/plain": [
+       "<matplotlib.collections.PolyCollection at 0x7f410cde6a90>"
+      ]
      },
+     "execution_count": 11,
      "metadata": {},
-     "output_type": "execute_result",
-     "execution_count": 11
+     "output_type": "execute_result"
     },
     {
      "data": {
-      "text/plain": "<Figure size 432x288 with 1 Axes>",
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAB6CAYAAACiANjQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAWDklEQVR4nO2debAc1XXGv9Ojp5TLiRMS4QiBnp8EyIYAwfCEWWRjFm2kMMaBACaE4NgKFDgQI8SisJsIJCQFFyREtqGAOBgv2EAlgIXBYGyDpaewi0UICW1sdiWOKy5TvDn5o7fbt9dZemZ65vtVqV7f2/fevt1T+vrM16fviKqCEEJIdXG6PQFCCCGtQSEnhJCKQyEnhJCKQyEnhJCKQyEnhJCKM6EbB500aZKOjIx049CEEFJZxsbG3lHVne36rgj5yMgI1q5d241DE0JIZRGRzUn1tFYIIaTiUMgJIaTidMVaaYXZzonFG0vyfWr1+F1tmg0hhHSfygl5O5gzdHK8MkH0xZHYvgd/c0dZ0yKEkKbobyHXeqJAa11Dkc5pazP3facBdl9xy+KE/R/41a2Nz5cQQpqgv4W8XdgiX9dQzEXi7UUw7/c/52774m628/o+8M6qEiZLCBk0+l/IC0baDWFH5E0g4mD+zmdGxxMJBP/+HTe1fAxCyGDQ/0KeQsP2irnPjMgT22pypB5pUodk3GDm73J2NJo3xrt/yw2ZYxNCBouBFfKGKRDVa70e8cmbwr8JZIwzf/g890Zi3izEwf0br2/t2ISQSkIhN2nGgjGi70wRDwS6oC1TryeLuSne/pjevOfvcYEl7p5N8/J1xY5JCKkkFHKTFv30QhF5ni2TZskk2TVJbQNxD333+XtdHGmrInjg+Wuy50kIqQxtEXIRmQfgBgA1AF9T1WvbMW6ZxPxxtzK9gynyTjzqdbtn9Pfb5UXkBayVWNs86uq+wysC9drP23dxTPAfeOqq/LEIIT1Hy0IuIjUANwGYDWArgDUicq+qvtDq2G0hRZwTH3Z2kKwHnQBCa8UQ2oCi8867efi/1+qL+/6XuULvhHUPrr2i2LEIIV2jHRH5QQA2qOpGABCRbwI4DkBvCHkKDYt4I5ZLVpScZq3kRdZm9O2PkSfuaRG7fypWNowmWDdzZ14ZE/cf/OzS7LkSQjpKO4R8VwBbjPJWAB9rw7idxRLqtkTrzWSwBFG0k1xftL+5bd886gjF3HxYqxoX8wRmH/plwG/mtV/9+OJi8yOEtJ12CHnS/3yNNRJZAGABAAwPDzd9sNX1bzfcZ3btpOhc2hWNF/XHTexj++Wifrc9RlJuuz1OLHpPaIOEiNyeU8oc1QGOPvwfvTHcsX/4yMXFzoUQ0jLtEPKtAKYa5d0AbLcbqeoqAKsAYHR0NCb0ZdLqaodz33da8cam0Nk+dzPZKkXbNGKtmPtiop8/jSIcdeSS8Mbgjfnw6ovaMzghJEI7hHwNgD1FZBqAbQBOBvDZNozbMxRd8XDeB86IVhSN0mNi2iY1tfLWA2E1rZWc/ppk1aSgSbud8NhHzL0ubCOCH92/qMAkCCF5tCzkqvqeiJwD4EG46Ye3qOrzLc+sghRZ8XDepAX5A2VlrDRK0sPVFGsliYhvbkf4IlD7htDAN4vDj1nqGnPe/B6974JCcyKERBHVjrocAFxrhb/ZGWf+Lme7GylrrEQi7EDknXheu3kDMNpHrA6zjbdtZ6ek1auIK8CmkIsE/rj/1MTvFxxXQg/d/ev1dSSI1NUx9gN4/O6FBa4cIYOBiIyp6qhdzzc7e4i0FQ/nTz03u2NSOmKS6GcQsVCMyNt+AKqGUDdEkT52GwFm/fn14Q0CwE++fX4TByekv6GQV4Cs1Q7nT1+Y+qZpqq2S1Nbe1whFMmAaHVK91CdDxFWAQ09aHrFzfnYnhZ0QWit9yPwZF0YzZWxbBYhZMEm2SpIVE4nIfWvFtFUAz0KRyINO01Zx6xD44+ED0Lj1ktTHPSaCOnXCOT15+5eav3CE9Di0VgaIpNUO5/1Jygs7zUbOLUTcidktmR3UFesC4xx0+orgBrDmVoo6GQwYkQ84kfVVgOTX9q0Ftxp+0GlEz+aDzCIReaSfH+E7yI7WzXLwEFUw9tW/b+VSEdJ10iJyCjmJMXfmlQCyRRxAqpCblkzEVom09cqOxEQ3yVqJCHuWkKf+lUgZANbdTGEn1YJCTppmziFXA8jxx726hoTcSRLh5Lo0fzwu+uaxMoQ84rG720/dSGEnvQ2FnLSV2bOuyX3QmZQ/DjT3oBNwBTdXyGFF7OacLCFPbO/VPX0DRZ30HhRyUipHHbEkORoH4iLt1Zm2CpAuzIX98ZS/bj/r2AkRvHmT8KN//3/Hsyso7KT7UMhJxzliznUN2yoxYff7+dE4EtobdfE23lj+8RL2x4TfrENYZ87puWUUdtJ5KOSk63xy/tLGhdz2x5Hcvog/HozppLVD3HYx6vxyELl7dc8voaiTzkAhJz3H4ccuSxZyhHVB/rghoOnZLWmCLqkRNxBPWYztT7tJmDcEY3v91RR2Ug4UctLzzPrM9e5GWoSeJOSmCNttAiEuKOSxB6lRayXLgkm8iTjAi5dT1En7oJCTynHYicszbRXAeyBpCXJq2qG3L88fBywhT3rJKOGGEDwgTdsP4KVLKeykeSjkpPIccsryTD89MaK2hbxFf9y0Vew2qTeFlLavXExRJ41BISd9x8dOW5EooInZKglCHemTFM37B3IQu2nYfXKtl6Q5mGMA2LCIwk6yKUXIRWQZgGMBvAvgVQBnqOp/5/WjkJMymHnGigTxlFShBuJWiNsHTfvjZv9CQh7UaaTtq+dzwS8SpywhnwPgYe/n3q4DAFW9MK8fhZx0ggO/sNLdSLVdEBH0WLQOFPPHU78NhGNFxjBvKsG2RvdLtO/GcynspAPWiogcD+AEVT01ry2FnHSDA85aGctuAVIEGfG6oB+sMYwbgJ1jnpi6mHRTcTT20DQ8pgb1APDaOfwxjUGlE0J+H4C7VPXfUvYvALAAAIaHhw/cvHlzW45LSDPsf44brcdeyy8o5G3xx5NsFWsM9xgavyF4wr/pTP6m6SDRtJCLyEMAJifsWqyq93htFgMYBfAZLXBnYEROeo0/Pc+I1iMiGmbEiIYZK2lt/HLbhNzfB3OfGvvhirrHpi9c0M7LQnqM0iJyETkdwJkAjlLV/yvSh0JOep19z18Zt0yAwmmHMRH3+sbH08h4pnXjrgVfQMhFI167eNubzljU0jUgvUdZDzvnAVgB4HBVfbtoPwo5qRL7XGCIupWHbopvXopjLFPG32/44/YNw85mUT8CN+fgmGLvCbmYbRSb/uqiVi8D6QHKEvINAH4HwC+8qidU9cy8fhRyUmX2vmRlVHhTfPRcW8VrE3vQ6d0wIFa0bggzgoeqGhsnWMI9GNevD4/z2qkXt+dikI7CF4IIKYG9Ll0ZE3LbIklKK4xkq3h9mvLHvbK5T+x2lpCLRAX/1ZMvKeHKkDKgkBPSAT5y5cpsf7yIkGfZKkWFPJKyqBCjLIaPDlE4jrriDuDlEy5t38UgbYdCTkiH+fDVK9NtFTPiTonGgQRbxdFMEQfginbkG0Io1ElCLp6Y++0do/364y9v/4UhTUMhJ6TL7LlkZb6QA/nZKkab0JcvIOTGfjFuFkFE7os8vL/itvOF/ZljryrlupDiUMgJ6SH2WLoyYqvYOeiN5I/HbBUg+qATobhH/HFDrE0hdzyB9/fVnHpQrnnHXHfMl8u5MCQTCjkhPcruy1dEH4YWyVbxy5bP3qg/Htgq3n7TVvGFW4w6v+wgLD8xd0m7LwlJgUJOSEWY/pXljWereG3FzjH3hdy3Vcw8c6/sWPsdI1IH3IjctlkcRMs1px5E8o8dtay8izPgUMgJqSDTblweinXEfmkx7TDFVknzx7OE3Bd6JxKpu3bMD49Y0dkL1udQyAmpOCM3h79par/NmSbkjdgqef64X/bFOigbfSaIIerevwnijvEfn/hKWZdmYKCQE9JnfOhryzJEHG3LVsmzVdKEHAAmeELvwBB2p47vHvrPnbxUfUOakE/oxmQIIa2z+fPhSocjty51fRbfL4e6CSuG921il1vFFG8fu2xy4k/Piuy/65Cb2zqfQYNCTkgfYK90OHL7tWFkXgC7aZ4om9ktQR3y+wBAXSXYdsSN2E978vNB3W0Hfb34xAkAWiuE9D3TvrGkoWwVtxxmqwDJaYdmFO5752JbKpatEq+LWi9DznhwDEfquPnAOzp1mSoBrRVCBhRzpcPd77omkloIhP54WI5H03UV1Iw6sxyLzAvYNmYbX8QBYNyL1h3vAekX150abN/w0Ttzxx1UKOSEDBCvnrQ42J7xnauD7SwP3cxeSSoDxWwVMxp36+pRQZd65GZh7hty6lj49Elw4Pa5dr/v5J/sAEEhJ2RAsVc63Pv7VyRG42kRtl1fNDKvQ2LCn4YfjY973wD8MS999njUJEyFvGyfewuN16+0RchFZCGAZQB2VtV32jEmIaSzvPDpK4Lt/e67LGbBmGmHPna2imOXDdsk+IuwbAp6zS4b7dx94bcA19rx+9Wx5IVjgvaL9r6/qfOvMi0LuYhMBTAbwOutT4cQ0guYKx0e8J//kNm20WwVswxB8MKQ7Y+buPvc7RrqGd8S6rjhxaNR80T/nI88nDn3fqEdEflKAIsA3NOGsQghPYa50uHBD14cicLNl4HMctLfJH/cLKf54z51OHAwnrivZtxM6uq4toso/vWlwwML5nMzHm/yCvQ+LQm5iHwKwDZVfVokO2lVRBYAWAAAw8PDrRyWENIl7JUOZz20KCbaPnk+eCySj0X2YWReE8UEGU9s77ez/9rc8crBwc3is3s8mTm3qpGbRy4iDwGYnLBrMYBLAMxR1f8RkU0ARot45MwjJ6T/OOqRL8XWVzFfy/eFd4IzHonGJzjjEX98gjOOGtQQZsWQuH18y8SvcyPvetA+/FsP2vsReU3qGPJuBjXUMSTvBcJ+3PSnOnehWqDpPHJVPTplwH0BTAPgR+O7AVgnIgep6hstzpcQUjHMlQ7/7LG/AxB9+FmHBALvE/fT4/trsTauKPvUrEi+hnhqpCvm9WC/eRwHddy3cT9MlPfgQDF72vr8k+0xmrZWVPVZAB/0y41E5ISQ/sZc6fD4n5wNIJqtYmILtS3MSWTZMEnH8Bk3sl2icwj7P7ppRiD6Hx/ZkDuXXoB55ISQUvneYTdFyqc8sQATjJeBxlUwJFaUbOx3DJvEjMZrVh/TjnHL9aBtEI37Swog+nC1htB+Ceqkjp9uno4hz6o5YLh3E/PaJuSqOtKusQgh/cudB68Ktk//+d/E9o9Dgt/QAMJslWC/OnACrzvfVgnHqAf7x+GghvFA4MeNY5jHN8d/+vWpGPKEfu+p24qebkdgRE4I6RrmSodnjp0WPMQEwujZz1bxI/ShoJycrRIs5GVG7kYE7j8M9anZZUO8g28NnvD7+1/aMgWOAENe25HddrRyGVqGQk4I6QnMlQ7P/a9TgtfyATfSNpe/rUndyxcvFpnbxB6g2v0tz30irNRHMdsCb2yb4m0Ldp7S+WidQk4I6TnMlQ4veuYEAMCQY4tptlgH7QpG4zU/fdFaCmDI+EYwDgluGpGXkBD+pOovtu+GIa/0gSmd8dUp5ISQnsZe6fCq5z4FIPnhpr3Oix9Zj6tgYuyFopyoHRrLfplo9bEzYBxjQeCaCH6zYxoAYEhqmDC5vAwYCjkhpFKYKx1ev35uZrZKUGdYNJH62MPR/CjfscS7Zv2+kl8eV4UjobjX35jh9p/8cu4xGoVCTgipLAv3ejDYvvHFI4PtmOdtZKdMlGi2yhDGA+Gvw3Hf+PSi8RrqxnZ0ZUYzs8aBG4FHjmmUnUjr9kMhJ4T0BeZKh7e8PCs3W6UIdraKSZat4pbLFW8TCjkhpO+wVzq8a8NMANFslZrUMWRko9SgGJL3AIRphxMtEfdtlTDPXLyxJCg7kKDsboeCXoatAlDICSEDwEl7rAm279m4P4DsFEUz1RFALFulDqDm7cuLuzsRmVPICSEDhbnS4erX9oqlHZr+OOBG53a2ylDKqt22vWJG42VCISeEDCzmSoc/3rRHZtusbBUHYpWduGdekq0CUMgJIQRAdKXDda8PtyVbpVNQyAkhxMJe6fClLVOC7UayVTqVuUIhJ4SQHD48dXuwvWXbLj2TreJDISeEkAaYumu40uEb26bASbFSOplH3vKRROSLIvKSiDwvIkvbMSlCCKkCk3fdjj+ashV/MGVLpL5T2So+LUXkInIEgOMA7KeqvxWRD+b1IYSQfsRc6fC9N8IMmLJtFaB1a+UsANeq6m8BQFXfan1KhBBSbcpc6TCJVuP/GQA+LiJPisijIjIzraGILBCRtSKy9u23327xsIQQQnxyI3IReQjA5IRdi73+OwE4GMBMAN8SkemqGvsJa1VdBWAVAIyOjub/TDYhhJBC5Aq5qh6dtk9EzgJwtyfcPxeROoBJABhyE0JIh2jVI/8+gCMB/EhEZgCYCOCdvE5jY2PviMjmJo85qcgx+oxBO+dBO19g8M550M4XaM85fyipUhJckMKIyEQAtwDYH8C7ABaq6sPZvVpDRNaq6miZx+g1Bu2cB+18gcE750E7X6Dcc24pIlfVdwH8ZZvmQgghpAk6m7VOCCGk7VRRyFd1ewJdYNDOedDOFxi8cx608wVKPOeWPHJCCCHdp4oROSGEEAMKOSGEVJxKCbmIzPNWWtwgIhd1ez5lIyKbRORZEXlKRNZ2ez5lICK3iMhbIvKcUfeHIrJaRF7x/u7UzTm2k5TzvUJEtnmf81Mickw359huRGSqiDwiIuu9VVLP9er78nPOON/SPufKeOQiUgPwMoDZALYCWAPgFFV9oasTKxER2QRgVFX79sUJEfkEgF8DuF1V9/HqlgL4pape692wd1LVC7s5z3aRcr5XAPi1ql7fzbmVhYjsAmAXVV0nIr8HYAzApwH8Nfrwc844379ASZ9zlSLygwBsUNWNXv76N+EuoUsqjKo+BuCXVvVxAG7ztm+D+5+gL0g5375GVXeo6jpv+38BrAewK/r0c84439KokpDvCsBcvX0rSr44PYAC+IGIjInIgm5PpoP8saruANz/FAAGYZ37c0TkGc966QuLIQkRGQHwUQBPYgA+Z+t8gZI+5yoJedLvKVXDF2qew1T1AADzAZztfS0n/ce/ANgd7lIXOwAs7+50ykFEfhfAdwGcp6q/6vZ8yibhfEv7nKsk5FsBTDXKuwHYntK2L1DV7d7ftwB8D669NAi86fmMvt/Y1z9Yoqpvquq4qtYBfBV9+DmLyBBcUfuGqt7tVfft55x0vmV+zlUS8jUA9hSRad5iXScDuLfLcyoNEXm/96AEIvJ+AHMAPJfdq2+4F8Dp3vbpAO7p4lxKxxczj+PRZ5+ziAiArwNYr6orjF19+TmnnW+Zn3NlslYAwEvX+ScANQC3qOo1XZ5SaYjIdLhROOAubvbv/Xi+InIngE/CXeLzTQCXw10e+VsAhgG8DuBEVe2LB4Qp5/tJuF+3FcAmAH/re8f9gIjMAvBjAM8CqHvVl8D1jfvuc84431NQ0udcKSEnhBASp0rWCiGEkAQo5IQQUnEo5IQQUnEo5IQQUnEo5IQQUnEo5IQQUnEo5IQQUnH+H1E382K4oDzvAAAAAElFTkSuQmCC\n"
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAB6CAYAAACiANjQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAWDklEQVR4nO2debAc1XXGv9Ojp5TLiRMS4QiBnp8EyIYAwfCEWWRjFm2kMMaBACaE4NgKFDgQI8SisJsIJCQFFyREtqGAOBgv2EAlgIXBYGyDpaewi0UICW1sdiWOKy5TvDn5o7fbt9dZemZ65vtVqV7f2/fevt1T+vrM16fviKqCEEJIdXG6PQFCCCGtQSEnhJCKQyEnhJCKQyEnhJCKQyEnhJCKM6EbB500aZKOjIx049CEEFJZxsbG3lHVne36rgj5yMgI1q5d241DE0JIZRGRzUn1tFYIIaTiUMgJIaTidMVaaYXZzonFG0vyfWr1+F1tmg0hhHSfygl5O5gzdHK8MkH0xZHYvgd/c0dZ0yKEkKbobyHXeqJAa11Dkc5pazP3facBdl9xy+KE/R/41a2Nz5cQQpqgv4W8XdgiX9dQzEXi7UUw7/c/52774m628/o+8M6qEiZLCBk0+l/IC0baDWFH5E0g4mD+zmdGxxMJBP/+HTe1fAxCyGDQ/0KeQsP2irnPjMgT22pypB5pUodk3GDm73J2NJo3xrt/yw2ZYxNCBouBFfKGKRDVa70e8cmbwr8JZIwzf/g890Zi3izEwf0br2/t2ISQSkIhN2nGgjGi70wRDwS6oC1TryeLuSne/pjevOfvcYEl7p5N8/J1xY5JCKkkFHKTFv30QhF5ni2TZskk2TVJbQNxD333+XtdHGmrInjg+Wuy50kIqQxtEXIRmQfgBgA1AF9T1WvbMW6ZxPxxtzK9gynyTjzqdbtn9Pfb5UXkBayVWNs86uq+wysC9drP23dxTPAfeOqq/LEIIT1Hy0IuIjUANwGYDWArgDUicq+qvtDq2G0hRZwTH3Z2kKwHnQBCa8UQ2oCi8867efi/1+qL+/6XuULvhHUPrr2i2LEIIV2jHRH5QQA2qOpGABCRbwI4DkBvCHkKDYt4I5ZLVpScZq3kRdZm9O2PkSfuaRG7fypWNowmWDdzZ14ZE/cf/OzS7LkSQjpKO4R8VwBbjPJWAB9rw7idxRLqtkTrzWSwBFG0k1xftL+5bd886gjF3HxYqxoX8wRmH/plwG/mtV/9+OJi8yOEtJ12CHnS/3yNNRJZAGABAAwPDzd9sNX1bzfcZ3btpOhc2hWNF/XHTexj++Wifrc9RlJuuz1OLHpPaIOEiNyeU8oc1QGOPvwfvTHcsX/4yMXFzoUQ0jLtEPKtAKYa5d0AbLcbqeoqAKsAYHR0NCb0ZdLqaodz33da8cam0Nk+dzPZKkXbNGKtmPtiop8/jSIcdeSS8Mbgjfnw6ovaMzghJEI7hHwNgD1FZBqAbQBOBvDZNozbMxRd8XDeB86IVhSN0mNi2iY1tfLWA2E1rZWc/ppk1aSgSbud8NhHzL0ubCOCH92/qMAkCCF5tCzkqvqeiJwD4EG46Ye3qOrzLc+sghRZ8XDepAX5A2VlrDRK0sPVFGsliYhvbkf4IlD7htDAN4vDj1nqGnPe/B6974JCcyKERBHVjrocAFxrhb/ZGWf+Lme7GylrrEQi7EDknXheu3kDMNpHrA6zjbdtZ6ek1auIK8CmkIsE/rj/1MTvFxxXQg/d/ev1dSSI1NUx9gN4/O6FBa4cIYOBiIyp6qhdzzc7e4i0FQ/nTz03u2NSOmKS6GcQsVCMyNt+AKqGUDdEkT52GwFm/fn14Q0CwE++fX4TByekv6GQV4Cs1Q7nT1+Y+qZpqq2S1Nbe1whFMmAaHVK91CdDxFWAQ09aHrFzfnYnhZ0QWit9yPwZF0YzZWxbBYhZMEm2SpIVE4nIfWvFtFUAz0KRyINO01Zx6xD44+ED0Lj1ktTHPSaCOnXCOT15+5eav3CE9Di0VgaIpNUO5/1Jygs7zUbOLUTcidktmR3UFesC4xx0+orgBrDmVoo6GQwYkQ84kfVVgOTX9q0Ftxp+0GlEz+aDzCIReaSfH+E7yI7WzXLwEFUw9tW/b+VSEdJ10iJyCjmJMXfmlQCyRRxAqpCblkzEVom09cqOxEQ3yVqJCHuWkKf+lUgZANbdTGEn1YJCTppmziFXA8jxx726hoTcSRLh5Lo0fzwu+uaxMoQ84rG720/dSGEnvQ2FnLSV2bOuyX3QmZQ/DjT3oBNwBTdXyGFF7OacLCFPbO/VPX0DRZ30HhRyUipHHbEkORoH4iLt1Zm2CpAuzIX98ZS/bj/r2AkRvHmT8KN//3/Hsyso7KT7UMhJxzliznUN2yoxYff7+dE4EtobdfE23lj+8RL2x4TfrENYZ87puWUUdtJ5KOSk63xy/tLGhdz2x5Hcvog/HozppLVD3HYx6vxyELl7dc8voaiTzkAhJz3H4ccuSxZyhHVB/rghoOnZLWmCLqkRNxBPWYztT7tJmDcEY3v91RR2Ug4UctLzzPrM9e5GWoSeJOSmCNttAiEuKOSxB6lRayXLgkm8iTjAi5dT1En7oJCTynHYicszbRXAeyBpCXJq2qG3L88fBywhT3rJKOGGEDwgTdsP4KVLKeykeSjkpPIccsryTD89MaK2hbxFf9y0Vew2qTeFlLavXExRJ41BISd9x8dOW5EooInZKglCHemTFM37B3IQu2nYfXKtl6Q5mGMA2LCIwk6yKUXIRWQZgGMBvAvgVQBnqOp/5/WjkJMymHnGigTxlFShBuJWiNsHTfvjZv9CQh7UaaTtq+dzwS8SpywhnwPgYe/n3q4DAFW9MK8fhZx0ggO/sNLdSLVdEBH0WLQOFPPHU78NhGNFxjBvKsG2RvdLtO/GcynspAPWiogcD+AEVT01ry2FnHSDA85aGctuAVIEGfG6oB+sMYwbgJ1jnpi6mHRTcTT20DQ8pgb1APDaOfwxjUGlE0J+H4C7VPXfUvYvALAAAIaHhw/cvHlzW45LSDPsf44brcdeyy8o5G3xx5NsFWsM9xgavyF4wr/pTP6m6SDRtJCLyEMAJifsWqyq93htFgMYBfAZLXBnYEROeo0/Pc+I1iMiGmbEiIYZK2lt/HLbhNzfB3OfGvvhirrHpi9c0M7LQnqM0iJyETkdwJkAjlLV/yvSh0JOep19z18Zt0yAwmmHMRH3+sbH08h4pnXjrgVfQMhFI167eNubzljU0jUgvUdZDzvnAVgB4HBVfbtoPwo5qRL7XGCIupWHbopvXopjLFPG32/44/YNw85mUT8CN+fgmGLvCbmYbRSb/uqiVi8D6QHKEvINAH4HwC+8qidU9cy8fhRyUmX2vmRlVHhTfPRcW8VrE3vQ6d0wIFa0bggzgoeqGhsnWMI9GNevD4/z2qkXt+dikI7CF4IIKYG9Ll0ZE3LbIklKK4xkq3h9mvLHvbK5T+x2lpCLRAX/1ZMvKeHKkDKgkBPSAT5y5cpsf7yIkGfZKkWFPJKyqBCjLIaPDlE4jrriDuDlEy5t38UgbYdCTkiH+fDVK9NtFTPiTonGgQRbxdFMEQfginbkG0Io1ElCLp6Y++0do/364y9v/4UhTUMhJ6TL7LlkZb6QA/nZKkab0JcvIOTGfjFuFkFE7os8vL/itvOF/ZljryrlupDiUMgJ6SH2WLoyYqvYOeiN5I/HbBUg+qATobhH/HFDrE0hdzyB9/fVnHpQrnnHXHfMl8u5MCQTCjkhPcruy1dEH4YWyVbxy5bP3qg/Htgq3n7TVvGFW4w6v+wgLD8xd0m7LwlJgUJOSEWY/pXljWereG3FzjH3hdy3Vcw8c6/sWPsdI1IH3IjctlkcRMs1px5E8o8dtay8izPgUMgJqSDTblweinXEfmkx7TDFVknzx7OE3Bd6JxKpu3bMD49Y0dkL1udQyAmpOCM3h79par/NmSbkjdgqef64X/bFOigbfSaIIerevwnijvEfn/hKWZdmYKCQE9JnfOhryzJEHG3LVsmzVdKEHAAmeELvwBB2p47vHvrPnbxUfUOakE/oxmQIIa2z+fPhSocjty51fRbfL4e6CSuG921il1vFFG8fu2xy4k/Piuy/65Cb2zqfQYNCTkgfYK90OHL7tWFkXgC7aZ4om9ktQR3y+wBAXSXYdsSN2E978vNB3W0Hfb34xAkAWiuE9D3TvrGkoWwVtxxmqwDJaYdmFO5752JbKpatEq+LWi9DznhwDEfquPnAOzp1mSoBrRVCBhRzpcPd77omkloIhP54WI5H03UV1Iw6sxyLzAvYNmYbX8QBYNyL1h3vAekX150abN/w0Ttzxx1UKOSEDBCvnrQ42J7xnauD7SwP3cxeSSoDxWwVMxp36+pRQZd65GZh7hty6lj49Elw4Pa5dr/v5J/sAEEhJ2RAsVc63Pv7VyRG42kRtl1fNDKvQ2LCn4YfjY973wD8MS999njUJEyFvGyfewuN16+0RchFZCGAZQB2VtV32jEmIaSzvPDpK4Lt/e67LGbBmGmHPna2imOXDdsk+IuwbAp6zS4b7dx94bcA19rx+9Wx5IVjgvaL9r6/qfOvMi0LuYhMBTAbwOutT4cQ0guYKx0e8J//kNm20WwVswxB8MKQ7Y+buPvc7RrqGd8S6rjhxaNR80T/nI88nDn3fqEdEflKAIsA3NOGsQghPYa50uHBD14cicLNl4HMctLfJH/cLKf54z51OHAwnrivZtxM6uq4toso/vWlwwML5nMzHm/yCvQ+LQm5iHwKwDZVfVokO2lVRBYAWAAAw8PDrRyWENIl7JUOZz20KCbaPnk+eCySj0X2YWReE8UEGU9s77ez/9rc8crBwc3is3s8mTm3qpGbRy4iDwGYnLBrMYBLAMxR1f8RkU0ARot45MwjJ6T/OOqRL8XWVzFfy/eFd4IzHonGJzjjEX98gjOOGtQQZsWQuH18y8SvcyPvetA+/FsP2vsReU3qGPJuBjXUMSTvBcJ+3PSnOnehWqDpPHJVPTplwH0BTAPgR+O7AVgnIgep6hstzpcQUjHMlQ7/7LG/AxB9+FmHBALvE/fT4/trsTauKPvUrEi+hnhqpCvm9WC/eRwHddy3cT9MlPfgQDF72vr8k+0xmrZWVPVZAB/0y41E5ISQ/sZc6fD4n5wNIJqtYmILtS3MSWTZMEnH8Bk3sl2icwj7P7ppRiD6Hx/ZkDuXXoB55ISQUvneYTdFyqc8sQATjJeBxlUwJFaUbOx3DJvEjMZrVh/TjnHL9aBtEI37Swog+nC1htB+Ceqkjp9uno4hz6o5YLh3E/PaJuSqOtKusQgh/cudB68Ktk//+d/E9o9Dgt/QAMJslWC/OnACrzvfVgnHqAf7x+GghvFA4MeNY5jHN8d/+vWpGPKEfu+p24qebkdgRE4I6RrmSodnjp0WPMQEwujZz1bxI/ShoJycrRIs5GVG7kYE7j8M9anZZUO8g28NnvD7+1/aMgWOAENe25HddrRyGVqGQk4I6QnMlQ7P/a9TgtfyATfSNpe/rUndyxcvFpnbxB6g2v0tz30irNRHMdsCb2yb4m0Ldp7S+WidQk4I6TnMlQ4veuYEAMCQY4tptlgH7QpG4zU/fdFaCmDI+EYwDgluGpGXkBD+pOovtu+GIa/0gSmd8dUp5ISQnsZe6fCq5z4FIPnhpr3Oix9Zj6tgYuyFopyoHRrLfplo9bEzYBxjQeCaCH6zYxoAYEhqmDC5vAwYCjkhpFKYKx1ev35uZrZKUGdYNJH62MPR/CjfscS7Zv2+kl8eV4UjobjX35jh9p/8cu4xGoVCTgipLAv3ejDYvvHFI4PtmOdtZKdMlGi2yhDGA+Gvw3Hf+PSi8RrqxnZ0ZUYzs8aBG4FHjmmUnUjr9kMhJ4T0BeZKh7e8PCs3W6UIdraKSZat4pbLFW8TCjkhpO+wVzq8a8NMANFslZrUMWRko9SgGJL3AIRphxMtEfdtlTDPXLyxJCg7kKDsboeCXoatAlDICSEDwEl7rAm279m4P4DsFEUz1RFALFulDqDm7cuLuzsRmVPICSEDhbnS4erX9oqlHZr+OOBG53a2ylDKqt22vWJG42VCISeEDCzmSoc/3rRHZtusbBUHYpWduGdekq0CUMgJIQRAdKXDda8PtyVbpVNQyAkhxMJe6fClLVOC7UayVTqVuUIhJ4SQHD48dXuwvWXbLj2TreJDISeEkAaYumu40uEb26bASbFSOplH3vKRROSLIvKSiDwvIkvbMSlCCKkCk3fdjj+ashV/MGVLpL5T2So+LUXkInIEgOMA7KeqvxWRD+b1IYSQfsRc6fC9N8IMmLJtFaB1a+UsANeq6m8BQFXfan1KhBBSbcpc6TCJVuP/GQA+LiJPisijIjIzraGILBCRtSKy9u23327xsIQQQnxyI3IReQjA5IRdi73+OwE4GMBMAN8SkemqGvsJa1VdBWAVAIyOjub/TDYhhJBC5Aq5qh6dtk9EzgJwtyfcPxeROoBJABhyE0JIh2jVI/8+gCMB/EhEZgCYCOCdvE5jY2PviMjmJo85qcgx+oxBO+dBO19g8M550M4XaM85fyipUhJckMKIyEQAtwDYH8C7ABaq6sPZvVpDRNaq6miZx+g1Bu2cB+18gcE750E7X6Dcc24pIlfVdwH8ZZvmQgghpAk6m7VOCCGk7VRRyFd1ewJdYNDOedDOFxi8cx608wVKPOeWPHJCCCHdp4oROSGEEAMKOSGEVJxKCbmIzPNWWtwgIhd1ez5lIyKbRORZEXlKRNZ2ez5lICK3iMhbIvKcUfeHIrJaRF7x/u7UzTm2k5TzvUJEtnmf81Mickw359huRGSqiDwiIuu9VVLP9er78nPOON/SPufKeOQiUgPwMoDZALYCWAPgFFV9oasTKxER2QRgVFX79sUJEfkEgF8DuF1V9/HqlgL4pape692wd1LVC7s5z3aRcr5XAPi1ql7fzbmVhYjsAmAXVV0nIr8HYAzApwH8Nfrwc844379ASZ9zlSLygwBsUNWNXv76N+EuoUsqjKo+BuCXVvVxAG7ztm+D+5+gL0g5375GVXeo6jpv+38BrAewK/r0c84439KokpDvCsBcvX0rSr44PYAC+IGIjInIgm5PpoP8saruANz/FAAGYZ37c0TkGc966QuLIQkRGQHwUQBPYgA+Z+t8gZI+5yoJedLvKVXDF2qew1T1AADzAZztfS0n/ce/ANgd7lIXOwAs7+50ykFEfhfAdwGcp6q/6vZ8yibhfEv7nKsk5FsBTDXKuwHYntK2L1DV7d7ftwB8D669NAi86fmMvt/Y1z9Yoqpvquq4qtYBfBV9+DmLyBBcUfuGqt7tVfft55x0vmV+zlUS8jUA9hSRad5iXScDuLfLcyoNEXm/96AEIvJ+AHMAPJfdq2+4F8Dp3vbpAO7p4lxKxxczj+PRZ5+ziAiArwNYr6orjF19+TmnnW+Zn3NlslYAwEvX+ScANQC3qOo1XZ5SaYjIdLhROOAubvbv/Xi+InIngE/CXeLzTQCXw10e+VsAhgG8DuBEVe2LB4Qp5/tJuF+3FcAmAH/re8f9gIjMAvBjAM8CqHvVl8D1jfvuc84431NQ0udcKSEnhBASp0rWCiGEkAQo5IQQUnEo5IQQUnEo5IQQUnEo5IQQUnEo5IQQUnEo5IQQUnH+H1E382K4oDzvAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
      },
      "metadata": {
       "needs_background": "light"
      },
      "output_type": "display_data"
     }
    ],
    "source": [
     "def left(x, on_boundary):\n",
     "    return near(x[0],0.)\n",
     "\n",
     "bc = DirichletBC(V, Constant((0.,0.)), left)\n",
     "\n",
     "u = Function(V, name=\"Displacement\")\n",
     "solve(a == l, u, bc)\n",
     "\n",
     "plot(1e3*u, mode=\"displacement\")"
-   ],
-   "metadata": {
-    "collapsed": false,
-    "pycharm": {
-     "name": "#%%\n",
-     "is_executing": false
-    }
-   }
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   },
    "source": [
     "Validation and post-processing\n",
     "------------------------------\n",
     "\n",
     "The maximal deflection is compared against the analytical solution from\n",
     "Euler-Bernoulli beam theory which is here :math:`w_{beam} = \\dfrac{qL^4}{8EI}`::"
-   ],
-   "metadata": {
-    "collapsed": false,
-    "pycharm": {
-     "name": "#%% md\n"
-    }
-   }
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": 6,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    },
+    "pycharm": {
+     "is_executing": false,
+     "name": "#%%\n"
+    }
+   },
    "outputs": [
     {
      "name": "stdout",
+     "output_type": "stream",
      "text": [
-      "Maximal deflection: 0.00586375434356426\nBeam theory deflection: 0.005859375\n"
-     ],
-     "output_type": "stream"
+      "Maximal deflection: 0.00586375434356426\n",
+      "Beam theory deflection: 0.005859375\n"
+     ]
     }
    ],
    "source": [
     "print(\"Maximal deflection:\", -u(L,H/2.)[1])\n",
     "print(\"Beam theory deflection:\", float(3*rho_g*L**4/2/E/H**3))"
-   ],
-   "metadata": {
-    "collapsed": false,
-    "pycharm": {
-     "name": "#%%\n",
-     "is_executing": false
-    }
-   }
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   },
    "source": [
     "One finds :math:`w_{FE} = 5.8638\\text{e-3}` against :math:`w_{beam} = 5.8594\\text{e-3}`\n",
     "that is a 0.07% difference.\n",
     "\n",
     "\n",
     "The stress tensor must be projected on an appropriate function space in order to\n",
     "evaluate pointwise values or export it for Paraview vizualisation. Here we choose\n",
     "to describe it as a (2D) tensor and project it onto a piecewise constant function\n",
     "space::"
-   ],
-   "metadata": {
-    "collapsed": false,
-    "pycharm": {
-     "name": "#%% md\n"
-    }
-   }
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": 14,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    },
+    "pycharm": {
+     "is_executing": false,
+     "name": "#%%\n"
+    }
+   },
    "outputs": [
     {
      "name": "stdout",
+     "output_type": "stream",
      "text": [
       "Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
       "Stress at (0,H): [ 2.09596028 -0.22361953 -0.22361953  0.77051432]\n"
-     ],
-     "output_type": "stream"
+     ]
     }
    ],
    "source": [
     "Vsig = TensorFunctionSpace(mesh, \"DG\", degree=0)\n",
     "sig = Function(Vsig, name=\"Stress\")\n",
     "sig.assign(project(sigma(u), Vsig))\n",
     "print(\"Stress at (0,H):\", sig(0, H))"
-   ],
-   "metadata": {
-    "collapsed": false,
-    "pycharm": {
-     "name": "#%%\n",
-     "is_executing": false
-    }
-   }
+   ]
   },
   {
    "cell_type": "markdown",
-   "source": [
-    "Fields can be exported in a suitable format for vizualisation using Paraview.\n",
-    "VTK-based extensions (.pvd,.vtu) are not suited for multiple fields and parallel\n",
-    "writing/reading. Prefered output format is now .xdmf::"
-   ],
    "metadata": {
-    "collapsed": false,
     "pycharm": {
      "name": "#%% md\n"
     }
-   }
+   },
+   "source": [
+    "Fields can be exported in a suitable format for vizualisation using Paraview.\n",
+    "VTK-based extensions (.pvd,.vtu) are not suited for multiple fields and parallel\n",
+    "writing/reading. Prefered output format is now .xdmf::"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": 8,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    },
+    "pycharm": {
+     "is_executing": false,
+     "name": "#%%\n"
+    }
+   },
    "outputs": [],
    "source": [
     "file_results = XDMFFile(\"elasticity_results.xdmf\")\n",
     "file_results.parameters[\"flush_output\"] = True\n",
     "file_results.parameters[\"functions_share_mesh\"] = True\n",
     "file_results.write(u, 0.)\n",
     "file_results.write(sig, 0.)\n",
     "\n",
     "\n",
     "\n"
-   ],
-   "metadata": {
-    "collapsed": false,
-    "pycharm": {
-     "name": "#%%\n",
-     "is_executing": false
-    }
-   }
+   ]
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 2
+    "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython2",
-   "version": "2.7.6"
-  },
-  "kernelspec": {
-   "name": "python3",
-   "language": "python",
-   "display_name": "Python 3"
+   "pygments_lexer": "ipython3",
+   "version": "3.7.3"
   },
   "pycharm": {
    "stem_cell": {
     "cell_type": "raw",
-    "source": [],
     "metadata": {
      "collapsed": false
-    }
+    },
+    "source": []
    }
   }
  },
  "nbformat": 4,
- "nbformat_minor": 0
-}
\ No newline at end of file
+ "nbformat_minor": 4
+}