fenics_demos.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# FEniCS Demos\n",
    "\n",
    "### Matt Hennessy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Recap of last week...\n",
    "\n",
    "### Steps needed when using FEM\n",
    "\n",
    "* Derive the weak form of the problem\n",
    "* Create a mesh for the domain\n",
    "* Choose finite elements (cell geometry + basis functions)\n",
    "* Assemble the equations locally and then globally\n",
    "* Solve!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# What is (and isn't) FEnICS?\n",
    "\n",
    "An open-source Python/C++ library for solving ODEs and PDEs using finite elements\n",
    "* Very high level\n",
    "* Parallelised\n",
    "* Well maintained, lots of documentation, large user community\n",
    "\n",
    "FEniCS is **not** a GUI-based program with pre-defined models like COMSOL or Abaqus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Example 1 - Poisson's equation\n",
    "\n",
    "Let's solve\n",
    "\\begin{align}\n",
    "\\nabla^2 u + f = 0 \\quad &\\text{in $\\Omega$} \\\\\n",
    "u = 0 \\quad &\\text{on $\\partial \\Omega$}\n",
    "\\end{align}\n",
    "\n",
    "For simplicity, let's take $\\Omega$ to be the unit square"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# The weak form\n",
    "\n",
    "Recall that the weak form is given by\n",
    "\n",
    "\\begin{align}\n",
    "\\underbrace{-\\int_{\\Omega} \\nabla u \\cdot \\nabla v\\, dx + \\int_{\\Omega} f v\\, dx}_F = 0 \\quad \\text{for all $v \\in V_h$}\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# FEniCS Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from fenics import *\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "N = 20\n",
    "mesh = UnitSquareMesh(N, N)\n",
    "\n",
    "P1 = FiniteElement('CG', 'triangle', 1)\n",
    "Vh = FunctionSpace(mesh, P1)\n",
    "\n",
    "u = Function(Vh)\n",
    "v = TestFunction(Vh)\n",
    "\n",
    "def boundaries(x):\n",
    "    return near(x[0], 0) or near(x[0], 1) or near(x[1], 0) or near(x[1], 1)\n",
    "\n",
    "bc = DirichletBC(Vh, 0.0, boundaries)\n",
    "\n",
    "f = 1.0\n",
    "F = -dot(grad(u), grad(v)) * dx + f * v * dx\n",
    "\n",
    "solve(F == 0, u, bc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0735267092333906\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.tri.tricontour.TriContourSet at 0x7fe8b40d3850>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAAD8CAYAAACVSwr3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2dbahm13Xff+vemTuakWXrxU4wkkzVICdxSwy2YpuQtk6Ca8mhqAFDHYeamAShEoV+tCk0KeSLQ74kIbaFMMYEikUgJlWKYmNSEhtc1XKK32TjMJVba6yAsaX4ZWY0L/fufnjuuXffc/fLWvvsfV7unD8IzX2e85yzn3P2/j1rrb322uKcY9WqVat8bU3dgFWrVs1PKxhWrVp1TCsYVq1adUwrGFatWnVMKxhWrVp1TCsYVq1adUxZMIjIR0XkOyLy1cj7IiJ/LCLnReTLIvKG+s1ctWrVmNJYDB8D7k+8/wBw7/5/DwEfHt6sVatWTaksGJxznwFeSBzyIPCnbqOngFtF5NW1Grhq1arxdarCOe4EnvP+vrD/2j/0DxSRh9hYFWyz/cZz8ooKl1+1alVMP3QvfNc59yrr52qAQQKvBfOsnXOPAY8BvHzrDveWU2+vcPllauvmc1M3YdW+9i5emroJzfTpax//fyWfqwGGC8Dd3t93Ac9XOG8zrYNyla+p+8McwVQDDE8Aj4jI48Cbge875465Ea019cOVc2cnvf6q9nKXLjc5r7bvjgmQLBhE5OPAW4FXisgF4HeB0wDOuUeBJ4F3AOeBS8B7WzU2pFIgLGIgn71p6hYsV5dfqn7Kmn2mBDJdXx8DEFkwOOd+NfO+A36rWot6qmUJNAHBggeuO7sz2bXl8tX2F2n1bCoBx++PVkjExkRNYNRwJZrJAoUqA7/hQJ9yIM5Nc7oXZkhZ+ogSIqm+28p9yWm2YNBAwQyDwoE/RUd2N81n8CxF8pLdEhn6bJNgSfU3IzTGBsQswZCDghoIRhAM7SRTD+a9M7N8nGZtXble9Lla998CmFCfUVkh/b6ZAYWcOzsqHGbXk1JQSAJBCQHr4K892E/K4G2pVvdICxztM48BJNbH1NZFBBJjWg+z6qUhKAyFgRYEpQBo0Yl3z25XP+dJ1vblXdVxlmelgUioz6SsjX5fjIKi69cJQLSGw6zAYFICChoYWEAwdPCPNdB3b1o2ULZf0g3wvkrubw4mqWeegka/X2lBEYREAhAhOGzdfK7azMSswRC1FiJQqAUEKwhqDPylD+oaqnUPNIDJPbMUOEL9IwYLLSjc2Z20BdEgLyOlWYMhqAIo5GCgBUEpAGoP+utnQstTTrZOXdFvc6C93ymAxJ51DBj9PqQBRR8SXR+OWg89OLS0GpYFBiMUUkDQwMACgtLBP9Ug371p3OJd2y/tDfq89T5pQJJ6ZjFohPpECBZ+/8pBQg0IJRxqaLZgOOZGKKEw1DrQwMAKgaGDf+xB3EK1v0MONJp7noJH6BlrYFECCXfTTtDFCAJiJLditmA4ogpQiAEhBwItBEoGf4sBf/1M9VOOplNX9Mdq710KILFnFgNGvy+EQKGFRB8QMetBoxZWw2zA4E9VHrEWBkIhBIQaMLCAoBQASx7kJSr9vimgpO59DBqhZxuCRQ4UKUikABFyLca2GmYDBos0UKgJBA0ErIO/1qDfPcGByG1lwDF3L2PgiD2zEDD8PpCzKFKWRAgQGjiMrXmDIWAt5KBgcRlSlkEOBloQlACg5WDfndAK2Ta4CmC7DymIxJ6BBhg5SGzOc/TaOUCUwGFsq2EWYIi6ET2VQMEChBQMNCDQQmDIwJ9yYA/VkLbnoJK6pzFohJ5XHxb9554ChRYQMTjAUddiSsthFmDQyAqFGkDIwSAHAisAag363YUtzNxW9H3NvYnBI/QcNLAIWRQpa8ICiJhrkVIyCaqyFgGGVlAoAUItGFghsLTBbpHlu6UgErunIWBoYNF/1jFrwgKInPXQdyuSLkVDd2IRYLBoCBRiQBgKAwsEhgJgb4GuxpZlmjJyfyzA0FgWIYui6wdDABGzHixwCMmfsqyR/ThfMCiXUfvWQm0opICQgoEWBFYILHHQa6T9XimAWIChAUUKEiWAyFkPOTiMrfmCYV/WzEZfGihYrIShMLCAoBQCS3U5NHGG1D2JQaN/P3KgSEHCAogSOKTkWw1jxBpmD4aUUtZCDgpjAEEzSK0AmGLg7505blZvXak7pZr7Xjlw9O9jKShS1oQFECHrIeZadHAoshq8OEPNDMhZgSFXsk1rLdSCghUINUFQCoDQIG6h2tfJgSZ1P0LQCN3nECz886asiSGA0FgPB8cHpjKn0KzAcCBFfCFmLZRAQWslxKyD3CDOwcAea2g/+C3XqGE95K6XuobGXdhcwz9f+jwxS0ILiJz1cOTciXiDbzWMmdcwTzDsy48vxKyFIVBoBYRaICgBwFgWQ+1r5uASu0boc6H72x/oWkiUAqIPBzhqPfStBku8AdrHGSYHQ8mGMtrCKkOhYAHCUBhYBlfLwb93uvzcW9fKLYfUd0pBo/+52LEpayAFiRQg+oHK3TOShYOvObsUk4Ohk7YkvNaFsEBBAwSrdVADBEUWw4CBPVStoGKxFjSBUg0kNICIWQ8pOPRjDj4cYlbDFFOXswHDgUILp4wVnLVQ0LgNFiCkYFATBLUHvzszrLpSX3LFvsw89p0swIhZC91xObfDH/RWQAyBg0VjxRnmBwalYtZCaiFUDShYLYTUYG8JgtqDvdW1cxAJffcYLHKg8N+PQcJiQQyBg6+Y1aB2JxqkRi8KDLnYgjZ5KaRSKLQAggUEUwKghmLtTwGjf380oNBCIhZTCAFiCBysVsPY7sTswVC6EYzFhchBwWIlxAZ9DRgMgcDWzvRz4wB7V3Wl8kLfNQYLDShS7sTeGReNQ4QA0QIOmlhDXy1nJmYLBu0OUroybN7xFaBQAwgtQDCXwZ9Sro0pcPTvhwYUfUjELIUYOEKACMHBPyYHB4tC7sQYcYbZgsEqTWZjX1YoWNyGUiBoYTAUAqfOjAeR61f0VbVD3ysGCw0ouvsdsyJiLkYIECk49I9JwaGm1dBKiwFDbhVlTDFroQYUagEhBwMrBMYc9Brl2pMDR//7a0DRh0TMioiBIPR6CA6bY8LtDuU6DNGYcYbFgCGkkBuhsRZaQMEKhFowqAWBnZ06He7qVXuXin2HGDA0oOjur9aKSAEiBYfNMYdwSMUcrFbDwTkmSHZaNBg6xaYoc9OTYIeC1kooAUIOBiUQqDXga19PA5D+99WAog+JnBURi0H04w+amQv/fc0y8hpqFYBUOeMicr+IfENEzovI+wPvv0JE/lJEviQiz4jIe6u3dF8pN8K3FjQuxFAo7J1xx1877YJQcGf2glDY2tk9+C+kU2d2D/5LaWfnevC/uaqkvf69iN2P1L0MPYPY89JYhH7/SPWd3UhfPNJHA5ZuKrBeOlunVRbbIrINfBB4G3ABeFpEnnDOfc077LeArznn/o2IvAr4hoj8V+dcEcq6GYnUly/ZK7IUChorwWIhDLUMSgf8uZ1p9yq4dDXfmUPfLWZd+Pepb03krIhYDKLvXmhmLg7fj09l1o43tJbGlXgTcN459yyAiDwOPAj4YHDALSIiwMuAFwB7740st04lNoXcCI0LcfT8R/9OQUHrNliBkIKBBQRTD/6UUm1LQaP//UOg6O5fyN3o7rsPiFgMou9e5FwLzWxFX7FYQ6fU4qqxApAaMNwJPOf9fQF4c++YPwGeAJ4HbgH+nXPu2MgQkYeAhwBuQreXhC+tGxF8P2ItaKckN8fmoRCLIYSgUMMyGAqCmxuB5KLCOvDV/x4aUKQAAWErIhaD8AERiz2kVm7G4GCxGlLZkNkAZOW0aA0YQnejP0LeDnwR+EXgJ4BPi8hnnXM/OPIh5x4DHgN4+dYdTZYBWq0FX6mYQikUUrGDmHJAsMKg1eAvvaYGGv53jEHCv08xSMRcjJx7kYODZrYipZzVcHCcMZ+hVnk3DRguAHd7f9/FxjLw9V7gA845B5wXkW8CPwV8fnALAyqJLxx8NmEt+GoFhRgQasFgCghYFWpjChYWSPQBEXMxYu5FHw4Qjzuk4FBqNXRKuRNjSAOGp4F7ReQe4NvAu4B39475FvBLwGdF5MeBnwSerdnQIapewTkDBYuVUAMIpTC4ZadtNOyHV/WVbvvfIQaKHCRKAJGCA4Sth4P3Ei5GTkNSpVsrCwbn3HUReQT4FLANfNQ594yIPLz//qPA7wEfE5GvsHE93uec+27DdgOHgccuvmBxI/aSMxRxa6EEClYg5GBgBUFrAFivqwGG/x1zkLACogQOvkKzFaFgZInVMKRWQ02pEpycc08CT/Zee9T79/PAv67btI20Zdw00lVxng4KKSBoYTAVBCwKtTEFixwkcoDQWA85OIRmK2qW0M+VgRtbi8l8tKyPOPK5QH9LWQsxWaEwJhBKYPDynXZbqPv6wVXdjmL97xADRXc/LIDQWg9D4GC1Gjpp3IkuANnNTIwxZTlbMAzJ7ApWZxpgLbSCwhAgWGAwFgS019bAwv9+IUikrIgUIGrDoYU6dyIXgGy5/Hq2YIipm5FIxRcOjjVYC2NCoQQIWhC0gsArTm+mwL5/TZdzkpIVFhpIaC2IkPUwBA5Wq6FTKM4wJ3diUjCUlI7PSWst5FKeQ7JCwWIlDAFCCQy6gd76c1qQ+N9BA4k+IHIuRs56sMJBo6DrUOhOjK3ZWQzayk0l0mwPF7MWWkGhFAgWGJRCoIZC187BQgMJKyBqw0FjNSxZswNDLWm3os8lMkEbKKRiCDEoaGEwJQg06rcvBYocJG7ZuaJ2L2rAQaNwZSidO1FLWzefY+/ipeLPn1gwpKQKRCorNdeEwhAglMLgtlPlnSelF6/r3US/7RpI9AFhsR6GwkFjNfjSuhND1KImwyLA0E1VxlKhu8Bj6WxEStqsRpgGClYgtAKB5joaWGiCnC/feSlqPUAYEH04wPGgZAoOMWngYFEXgOzPTIxdA3IRYOirPyOhUSi+oJ2J6KSZgfClgUJLIIwFAa367UmBImdFxKwHCLsXsZkLX6HpzE4xqyEmrTvRSRuAHKvM2yLBkFMqvmCxIFJl2HIuRAsoTAGD205fjL734rWbh53ba6sGEjFAlMIh5Fb4SlkNMZcipJTrMNcCLrMGQ8106OPnHmYt+GoNhRZASA34mufQwqNrfwkgUrEHKxxSLoUmEKl1J8asC1miWYNhqIbGFzqlXIhaUCi1ErRAqAGCEvWvmwOFxop4xenLauuhBhx8WVwKqzthUeu06BMDhtRqyi6+kN2aPpG30CkVV6gJhRQQNDCYCgQ5WUCRsiJS1kMNOPjSBCJzQci5Wwh9LR4M/VRoTf6CJnehk9aFKIGC1UrIAWEoDG7dHhaf+Mddeyar3+YYJHKA0LgWVjhorYacNK5FF2foApBzSI2eDxgihWBjytZ4NLoRVmshVU+hNhRqA2EoAKzn1QKj+x5WQGhdi1gyVEw+HGJWgyUI2VfOipiyNsN8wKBUbHOZmLTl20LyrQWLC9GXBgolVoIWCK1AoJV/fQ0kNIAIwQHyrkUfDrVdik5LT422rQqZUEPqPILegtBsKmuxFjSLoFpB4dbtS5NDoa+uTZp2pb5j7N5oZnD6zySVnh571v6PhnVXcs2aHV/qH0Oj1Z3S4iyGFsqlP2utBU1Rlb61YIVCDgg1QHD79o9Mx7+w+7Ki62gsiZT1ELIc4LhroQ1IdsrlN6SUcieGBiDHzH6cPRgslZssgUftUuu+LJu/5FyIqaFgBYD1PBZg3Lp9Kelm3Hb6YnU4+EplRsZiDdog5BLdilmCofW+fBppdptOWQs5F8IChZpAqAUD67U0kOi+h9V6iMGhL0sw0mo1BKtLG4BQmgHZqorTYmIMKfVzGEoTm0K+Yq1t5jUrJFtB4fbtHx38N5UsbcjFIEL3JHTvrAvMxtyTI5ZbY90kqZVOBBhiUhVmUS6vhnTegq8SF8IKBU0Ab2oYxGQBREylcOg/i1Qg0n/G/rP3fyxClmWqT2l/tEKlCg/OUVgY2aJZuhJapW6er+5haOMKGjfCl9+ZxoJCTENBcMd2eZLU93Zti6o0rkYq9hCKO8SmMlPxBmt+Q4mWFmdYNBhKFQJEzo3QWgtWWaCgsRC0GgIAyzm1sLh9+0dJOEA49qCFQ6li2ZChnIahcYY5aZGuhKUOQ0oWNyKmIdZCLShYXIY7ti82gUKN6+W+R+weaNwKi0uxhP0/W2uRYGipkBtRYi2MCYWcusE5JhCGtGEqOGgUmq62xhksqvUjaNWiwJDKANvtRXkt2WWlsxEpa6GWYgNBYyVMDYOYNO0qgcNQxayGmq5jSLXKA9TUiY8xTHHTa1kLKSikNBQGtxsG3gsFKyo7de2MxSG67xmKPYSCkpp4gzXxKSVtnCGmWrUZWujEg6HT0IzHlrKsjqwNBQsELJ+3AOOO7YvJIGUqMNlXLENyDuoCkbHU6FSSU267utpalCtRUyEf0BJfiLkRJdZCSCFroSYUbt++NBgKNc9f4lpoXYpcrMGXxp3Qxhk6zfHHKKfFg2EumWKl0roQtaDQGghDrlfiAoXulbU+Rav9PmsFIKfQbMHQshBsTrXSoPvq/3LVKr+mi/SPC4TS66e+SwyOGji0KqXfoq/M4cdutmDwNbQWQ0ratfQ13AiNrNaCFgpzUSs4WDXUnUjJWp9hjloEGKaSZYm1VTVciLxfXs9KuGNrjzu26nR4TbuscBjqUrRyJ1qqpVWtAoOI3C8i3xCR8yLy/sgxbxWRL4rIMyLyt3WbeVzadRJaWddH5FQadNQqPXCGAaGDgP+f5j2rcu2snYNhCULmZA1ALk1Z5IjINvBB4G3ABeBpEXnCOfc175hbgQ8B9zvnviUiP1bSGHd2WNLBsSWsvcVTB/9vEBSy/OIMtRZyUNCqlgUQOs/39vTG6O3bl6LTm7GpzNAUpja3Yckaq4qT5um9CTjvnHvWOXcVeBx4sHfMu4FPOOe+BeCc+07dZtpkrannKxRM0sQXfFmthSVDIXV+yzVKLIeSKczUs/DhXisLsvsRKpmyrG0VW6QBw53Ac97fF/Zf8/Va4DYR+RsR+TsReU/oRCLykIh8QUS+cM1N79MtLUhUAwo1YwXa62k1plvRaYg70ckyM9EBop+6P7e0aA0YQnMnffydAt4I/DLwduA/i8hrj33Iucecc/c55+47LfUq2i5JfTeiJGehVGMCYQ7X7TTFrlxL+9HpSwOGC8Dd3t93Ac8Hjvmkc+6ic+67wGeA19dp4jSqNSNRK+g41FqYenBqr2+1Gmq7E6s20oDhaeBeEblHRHaAdwFP9I75b8C/EJFTInIOeDPw9bpNnZdi+Qs1TNPamhoKnebSjlpqOZ09tbJgcM5dBx4BPsVmsP+Zc+4ZEXlYRB7eP+brwCeBLwOfBz7inPtqiwaH1qdrysanpJ1mql3AQ+tGDLEW5jYYNe2pYTVY5MNcE4Cco4bO6PWlypBwzj0JPNl77dHe338A/EG9ps1DVcu4VfZ1lwaFTnds7WWnM1NTmBr1py7Hmra0bHg7Z90wy66nUA1ftjQaXwsKt22FzbEX94YVEtDAoaWG1IUcslNVTqX7S9TWmhLtqdXiqZCGzkYkze2BULht68zBf0OOySnXzth3bOFOzFlTlHdbLYaKmmPgUaMhg7v/easl0dJyyG17V1OWnbCXoNViKFBqRWVMufiCJehYy1oY+oufOqflvKk2W6yGbNsSzyAWgLxRdcOBoUtRTSWghKahrFHpKebKrVBorTGu4auFO9GqyO/cdcOAYYnltUKaU22FmiqxGobKCu+5T1nW1KzBMMYefVOoVflzrcb+Jb9RteS06FmDoYaWainMcT+IEp10CJ3U7MdFg2EOtfFWTScNPKe2zmCZRWEXDYaU5raMNaUaQTNt4PGk/4KvqqMTC4aTqKUGHrUwqpGtOedEpyW5tSsYRtAU9QBOmpYKxaXqxIHBWtYttbKy9WamY2t1I1ZIa3XiwDCVlpoOPZZWKB1qCS7F5GCQc2fzB61Kaq7Lq1ctV5ODYcmK5dS3SIce6mMv5Rd7hdw8tBgwXL+BchZOSnLT0rS6g4daDBhOXZmfX/aDq+FK16UFQDqFNlgZUs0IhhdWGUtTFm/5/rXVre20GDCsimvKwbTKrq0r87d+1x61ahTNxWI5SdvVtdSJA8PWRP1vCWboXAbnqvnrxIHBqlRF31YFP1fZNTTGMgctwYXodMODYQzd6Oar1lKpESvp74C9qkwrGGYgbWdO/WpqB9VJcidCszdz1Na15VgKnU4sGLZP1jKHVYUaq0r0SdOJBcPStZRfwxtdV6+ezB0YJgeDu1SebbatSHqaY8Bn6l+xMd2Jk+S6WCVXJh9exVpuy1edKKViJK1mJKwZqhdvoFmqWYNh+3K9LeP6lkOK5iHz0NophqZFx1QjAAnj/JKPbS20mJH44dVlLD6rrVmDoYVqRIhbdJZQp47FGWr9gr64d6XJ4LWet8RaKInBpKaN/QS12BqYlFpuT7f9Un7FqVyuG22/4cBQW34nSmU/jpXLUJIL0A3kIZCocY7amjqWs2StYJhIQzttLZeiL8sArwGDWtaC1Y1o5eoN1fZMuLqCwVNtc9Da+SzuRE41sghDA7+mZbDUVaE3Qqr8Mp/MiJpzJ8jFGmoOvNpuQq5ttWcihrpy2uBzau3NkrSCYUT1O2fInagdWZ/jr/IQKGjdCIurFgs83qgzEqAEg4jcLyLfEJHzIvL+xHE/KyK7IvLOek3Ma6hfpqW8/6vhdxptALJUKXdC88s6JzjMoS214gsnNesRFGAQkW3gg8ADwOuAXxWR10WO+33gU7Ub6UszdTMnjRHkWgocNG0Y21pYFZamt7wJOO+ce9Y5dxV4HHgwcNxvA38OfKdi+yZTq1+DUneixtqJKeEw5bWnWPa+5HRo0IHhTuA57+8L+68dSETuBH4FeDR1IhF5SES+ICJfuObCpdfH1NIe3lCXAsYfoN/b21Jf02ot1FAN12/IbFZXcWxuq4E1TyyUKthfvfSHwPucc8kcZufcY865+5xz950We3bZGKoxZWntbJYgZC04jAEIyzVKoFDiRqRcu1jg0Y8tWWepukzbLiW/+78GBKcmzGnQPLkLwN3e33cBz/eOuQ94XET+L/BO4EMi8m+rtNCg7mb3Kdx/KNa0aL8zaAKQffU741DTtgYc4BAQ/n/FbRpwrpZQGNONGGOqcvulw9/erSvXm11H40g/DdwrIvcA3wbeBbzbP8A5d0/3bxH5GPDfnXN/UbGdi9P3r501bWDyj7vnuLW329QLuy8r2ta9G2glu1eFBnR/d6ha1kY2D6Oy+9AH9BA34iTPSIDCYnDOXQceYTPb8HXgz5xzz4jIwyLycOsGjqUQ7Ws+fI3VUMulOPxsneh8Laui0wu75wZBoYW1oHEjbiSper5z7kngyd5rwUCjc+7Xhzcrr1NX4PqI+SeXru5wbmfjm1y8usPN+//+4dUz3LKz8V1+cPWm6H6WLfS93Zuz29kNsR5aSDW1WgEKfdWwFkLxhVBMamhQW1OAqLUWEZb3/araCj3EWmsm+p2vttUAenNb8yvdUtrr13IfbvTK3EM1KzDUXlO+VJXAYa6AsFwv9x1quRB9YNdIgz4payQ6zQoMvrQR1zmYXanZiRKrIabcOgrLr21rOFgBVAsKIZVmn8amKbWxp9hUZYnG3tT5ZIdWA9q6IuydcWxdE/ZOH73Ze1e32do56rZcvXqKnZ0NpGJxhqF68drN3Hb6aKwgNEsB+ZkKTdzh8FzhAWONRwyFTAkUYlqCCxHLYZhLLQa4gcDQAWEM9YOQ/anLF6+f47ZT+cE3BhxCCg30DhY1LQ3dzEoYCloXIhd0LCnj1mlo4HGqfVY1mq0rEVLKnOrT1pJqWiMAOXSJrvWXTuNW1MwDqBmb0LTthd2XmaBQQ7WyHa2aWzo0LAwMGrW4yb5PGcuC7MsaawD9LEWn1ODp1A3COWxgo21H6jvF7sdY1kIovhAKPNbalm6q1cQnDgwWhR7ekOiy1WqwwCEHCI2mgITlmjnQtYSC1VpoWRV6DrqhwZBT7OHXshpADwcYbj348gdsTVCUnjcHhCFQqKXcbETIJS2dkShaQHW5XnLdIoOP2y/tsXtTOdNCgUi5soU7Ezfb/NmJlPxMSMgHIiEcjAzNVMAhHEJBSTgcYNY1FrlBfMf2xSaWRg5mKRhqoVBqLWhktTD7C/2GSl5qE6C4oSwGLbWtD9vamYZaDpAPwHUWRK0aknUDmfm25dynFlDoK5a7kHMjhsQXNHk5NXdoi2nRYNCaW6mAZO4h+p0gFoTsq9/ZhkyJ5eCgidDXBESpLKDK1lSoBIW+hiyYKl0fMccZCVg4GDrFKGs114Ysful3qhwctFYD5KcytdN3ta0IyzU10lgJNWMKWmshppxlaYkrzCm5CRYQY9i+vMvuWZ1pv30Fds9sKLybea65hKdQFiTEMyE10sYbgGDMAQjGHSAfe+grNlhL6j+kzqdViZUAcShYXYiUtTDUjRiSCu2r5WLCvmYPhtoKASGUHu3r+pVtTp3JP5R+mnQ/EBlSrKBLLDsyFpTsZAVEX2O7HKVAgHIo5GRNaNJamnN1G0KapSvRItKqfSiah6yNNYAu3hCrDZByLTTuhTYGMba0basBhZBKrQVfQ/JdrC7u2AuoYKZgCKnUjEo9hJyJ5z/8lAmZy2toAQfQAQLmAQlLG3LfywIFyywEpEFvSWrSuA+5H6sp91BZnCtx6orj+hndTc/FGXzl3Im+UnkNmpWXoWpPXcfWxh0O3s/EH3yFBmap22G5hkY5yFmtBA0UhsYWfAuzNL7QDzzWqA69d3HYM10cGDoNSXLKBR79ZCc/CJmKNfQDkZp4Q6wUXCruAHlAgA4Snaa0JDQWT8pq0lgJoIOCtuZCyWyExpWdQ22RTosFQ6eu9uP2FcfuGTmYmdDKuhzbh0PfaiiFAxC0HmJVpnOAgHJItJZlFakVCFAHCn2VWgsh1cp4bFk6Hk4AGCzKuRe+O2IlNNIAABJySURBVBGzGvrKpUprZyosrkUnDSDg+GAcExTW5eSanAStlQBlUEi5ELVKuM19huKGAMPWFdg703+tvHCLxaWAMBwAk2sBeUBAHhKQH6wWcNSqmFQKBGgHhb58KMSsBd+NsOw6NUQtaqXOGgxbV66zd6a8idoApA+JGlaDBg5gcy0gDwiwQyJ4jpHKo2kzFlNTkC2hUHtpdciNiGU8djMS3VRlNys3xjoJmBMYLr8EZ+vsZ5mKM1hnK3ylApE14QDpPSo0gIDjA68UFLVkSV3O5SNYgAC2mAKkXQiLtWDVlPtV+poPGCqoC0CGlHMnNFZDXzXhAGHXAsLWAxwdPJrt8MYGRckaBk2CUgsoDJmFsCrkWsxpRgIWAoZuvcT2S7vs3rR9kMvQTVmmdqUKWQi1rAaoBwdIWw8QBwTorQhfloHbQaRVEZQWQIAyKKRcCKu1MPcgY0yLAENLlVgNreEAx60HsAGikwUUKdUGgnaLuNw6B62VAGVQqGUtWOILc9CJBUMozhByJ7TKwaGvEByAY3kOQDEgIA0JCA/AWrDQqmSfyFIgQDsotIotzM2NgBMIhlCcIedOaKwGsMUbIJzjELMewA4IsEGik3agpgAyZAv5mDSrIK1AAN2UpAUKWgVjCQnXwhp4bFXWDU4gGGorBYeQSwEk3QoIWw+gAwToINFpyO7bLQZ/J8tS6BIYgD7IaIWCNbZgdSNiU5UauUt1rMETA4ZQADLnTmisBkgvsAq5FKGYA6CyHiANCNBBolNtWFhVWtYutwrSCgQYDwoxDZ2NGCuHAWYIBrl8FXdWP2UQWkyldSdSSmVGauINFtcCjlsPkAcE2CDRqWSwdjAZUr8yJe2eHK2AADYopDTEWtCo9ToJmCEYOslLV3E3HT7o2JRliTRWQ//vXLyh61Sp2QpIWw+QBgToIdFJC4ucagPBsjlPrkhrKlmpFRSs05O5acvBiU0V95SAGYOhkyUtOudODK0FqQlGauIOkAcE5CEBaVBAfADWAkZOpXt6aio21wAC1INCXzlrIeRG9OMLU0k14kTkfuCPgG3gI865D/Te/zXgfft//gj4D865L9VsaEoltRliVsPmPX28oetEGtcC9ICAtBXRSWtN9DV0E97aspRutwIBdFYClLsPUG4tJM85YgFYX1kwiMg28EHgbcAF4GkRecI59zXvsG8C/8o596KIPAA8Bry5RYN9pdwJP84QsxpScPClCUZqXAuIr8zUAAL0kOhkgcVYsu7fkFvXkEpn1loJoIPCEGshprmsj/ClsRjeBJx3zj0LICKPAw8CB2Bwzn3OO/4p4K6ajeyrizOElEqPzp5XGW8AHRzAZj1AGhBwfJDkytenBmFraAzZwEVTqTm3l6TWSoDhUMgFHCHvRqQ05owE6MBwJ/Cc9/cF0tbAbwB/FXpDRB4CHgK4iXPsXbzE1s3tyopZrYa+NHAAVHEHCFsPkAYExCEBemsipCEDt7YsJdtLgAD1oNBXylWwuBFTFn/tSwOG0B0J4k5EfoENGH4+9L5z7jE2bgYv37rDHF1JbT7jxxmsVoM23hD6G3RxB8gDAvKQAL014csKjVay7tkAeRBAeuGTFgiQDjSGXusPfKu1EHIjUoHH/lRliyItoAPDBeBu7++7gOf7B4nIzwAfAR5wzn2vTvM2Ss1M5KYtNVbD5ho2OAAq1wL0gIA8JEAPCl/WAakBSckg10gDAsgXUhkCBDgOhZT7sHmf5PtWpQKPLdOhQQeGp4F7ReQe4NvAu4B3+weIyGuATwD/3jn399VbqZR1diIVcMzBIfRayLWAPCCgHBIQH0gaYMTUatD70gLAV4llAPG1Dlq3YQgUNNZCyo0YO74ACjA4566LyCPAp9hMV37UOfeMiDy8//6jwO8AdwAfEhGA6865+4Y2rp/k5CsVgISj7kTMaoB0vCEEB0DlWoAeEGCDxMF3zGyblxt4Q8ChUcnA76Qtq5ZbFm0BAuStBCi3FHIuxOb1eay0VD0559yTwJO91x71/v2bwG/WapQ1LdqaBZkq/ZaDw+aYvGsBeUBAGSQ6lcDC15CBW1uW+oqaGglDgRB7zQoFTfpzyFoIuRFjpEJ3mk/PAFXdx9zu17EgZH/9hDbe0L0POtcC9IAAGyQgDQpIDzALNGprSGHVISDolEpU0lgJMAwKGmsheM0J3AiYCRjcpcvIufQS31AA0ncnNFZDqiZkCA5gdy2Sr3udswQSEB4kOVh0ql31uIUsFZM0NRIsFkLs9ZBrUAsKvrVgXWZ9MCNReZ0EzAQMtZWautTGGzppXQvQAwJskIA4KCA9mLTQGFNDyqUNgQHkshXrQyGl0ryF1jMSsAAwhAKQvjuRC0JCOq8hBAewuxab49KACL0HeUhAeECkYHFwvsJBqAFK7erJviwVk3LrGXK1Eix7TVqhoHUhckHHMeMLMGMwWAOQcNSd6E9dauMNB68pXQvQA8J/L/p+r5OnduCODR4NMHJqOej7KimbplncZLUOQA+E0LGWOgsaa2Gq+ALMGAwh+XGGXBAS7HAAu2vRHQdxQEAZJMAGik7agVYDIDXakZNlhWOJddBJ4zbEjg1BocRamGo1ZV+LAkNMqSBkKukpWOlJ6VpAHBD94zef0UMidgykB4kGGr5qDdxasgAA8hDQHGexEGLHp9wH0AUcU/LdiFR8oVa9R1gIGFKJTiGlZihCwUhA5VqADhD+8f3PbD6Xh0CsM6c24rUMLCtESmQd6MlzVYBAp1Sq8hAgbF4rh4JvLQxxI/YuDt9hbBFg8BVzJ/pByFy8AeKzFZu/N/+3AALSkEjVe/CVHPyJjm/ZvbvmoK2lkj0ZtJ/JrVvQxBCOvJexEqDcUvChELMWWi2e6rQ4MKSUmqEIuRSpqczD1zb/1wAC0pAIdbQQLEqshdTnYrKAxKrSTV1rna8EBLnPaiwECMcTclAoji00yGGAOYLBy370ZyZi7kQ/CGmNN2hci83r7L/eu37SZTj8d6zmQ78TpmpS5gaHdaDXHrxDVNoW7QpGTUUlCxA2r9eHQsxayKlmfAHmCAaF+lmQqRmKEByA6GzFwXFGQIAeEp1ysQlfmtL3pYOrpuXQEjbWZcw5GFhchaPv6YAA9SyFMZKafC0SDDml4g2Hx+TjDpAHxOa9QBsSkOjU77ipfTU1g6J0B+85WA5Dahdo6ytqrtECCJvj8lDQWgut4wswIzAcWS8RcSd85ayGEBwAlWuxOe7o9WKA2Lx3+O8cJA6OC8YWjr/WSbMZ79DCIKVgqXFtrSxFVqENDDrVhEJKUWuhUXwBZgQGjfpxBiscQOdabI47/HcoBtEpZUls3g9/F0tsAXQDonQn705jDe6QrAO+L1NtxZyLkckt0AJhc6wOCqWxhVZaFBigDA6AyrXYHHc86SdmRcDRTpSzJjbHHD8HpDu29pd86ODqKwea2tfTqARe6kVNFYGwOd4OhZzGcCNgBmDwK0Vr3QkrHEDnWmyOO3zIWivi4LMZa2JzzPHPbY4Nvw76wTDEFQhp7IFfw2IxrVdQZB2mUplT6x1KoXCs2OvIQcdOk4OhlTRwgKMPsMSK2Hwu0oZA54jWg8h06BQ4Ds4xj0LQTVW6Iax2H4dcERUrDDaf0VkJOSiMZS3AgsCgsRqAY5YDkHUtOmmsiM1n45ZEJwssDs6brHYdfSspDVDG1tDdnsPntE25aqso5VZBJjMYDW6DGQpe4NHPYaiRDg0LAoNWwUpPEeuhk8WK2Hw27m4cniPcvtR+F9rOnQLI8XOqD52lrAM+JFMpNcVyaCsMDt5TWAkwraXQad5g6NWAzFkNnWJwAIKJUBorolPOmuhkBcbh+dPvQ53B0pcFNjG1aFdOpfs+aqsnaVY/5qYetUCAYVCoZS3A3MEAKjgAxwAR26TGf0gWK6JTzpo4PFe846X2vrB29NJ9OvuaYlCHVHuDV2v5NG35dk0eggUIYIBCxI2oqfmDQakQIEJxB18aK8JXyuXoqxQah9fS1UqY407JY2joPo/W/RtKrAJfudyEEii01OzAEKwYHbAagGBGZMi98B+K1Yo4clyvc6Q3vEl3vGxF6wobnFp25RpTY23eWrJ5iyUjUZt/oHUbDl5XQqGVtQAzBAPoyslDPF06VdjFAolO8Y104x0jX6DW1mktG+p0mtPuybVUc6emkqXOQ2EA6dyEqS2FTrMEQ1CRzWhi1kP/5seClL5yLkdfFusi+PkMPHxNvXXZ9TMyeRssGlI70Vo9SZvCnEtWskChpbUAMwGDn/3YKepSgAkQB+97D0VjTRy0LQILSHegXKFaGNZ5wQaWoRobCmMURS0tn2ZZy6DNXJyLpdBpFmAAAxxABQjQQQLioIB4J0gBA+ydTgOSY9eYSUXhOalmyfXSxUyDYdBpIijAjMAQUzLekNnrUgMJCD/IXPFZS6fJQQTqdmitSmCk0ZT7IeRUa+ViyRoGdU5CAggxF6JmDgMsAAyggIOvCChiD0VrVQQ/q6xcXXsZrQY0Gs15AMN0y4+HLFwyZykqrIJcPKE2FGAhYIDDm5OdrfBvdGbnbAg/SO0OWDVWvlnK4neaw3r9OWrMlYjFacoG90ATYGwBBVgQGDpppzIBMyQ6WR66dRu9Y9eaaFntqqNqth7BGCewzDa0ggIsEAyQvnlqlyMmA0BgmgUuMBxIc9RU99KsgUHBGlONLaEAMwND7Mv2ZytS8m+62rLwVTsSbASNVosZRFNrwsh+p1IQtB78Kc0KDDGFpjI1sj6QIpDkNIOOuaqeWicWTQkDXyowiMj9wB8B28BHnHMf6L0v+++/A7gE/Lpz7n/XbKj2hpUApFPrhz5XWYF4o96nIZrLgNcqCwYR2QY+CLwNuAA8LSJPOOe+5h32AHDv/n9vBj68///R1eoBDAHO3LUO9OUN3NbSWAxvAs47554FEJHHgQcBHwwPAn/qnHPAUyJyq4i82jn3D9VbPJHWjrPqRpIGDHcCz3l/X+C4NRA65k7gCBhE5CHgof0/r3z62se/amrttHol8N2pG6HUktoKy2rvktoK8JMlH9KAIbTet7+iRnMMzrnHgMcAROQLzrn7FNefhZbU3iW1FZbV3iW1FTbtLfmcppLHBeBu7++7gOcLjlm1atVCpAHD08C9InKPiOwA7wKe6B3zBPAe2egtwPdPUnxh1aobTVlXwjl3XUQeAT7FZrryo865Z0Tk4f33HwWeZDNVeZ7NdOV7Fdd+rLjV02hJ7V1SW2FZ7V1SW6GwvbKZSFi1atWqQ82zWuiqVasm1QqGVatWHVNzMIjI/SLyDRE5LyLvD7wvIvLH++9/WUTe0LpNMSna+mv7bfyyiHxORF4/RTu99iTb6x33syKyKyLvHLN9vTZk2yoibxWRL4rIMyLyt2O3sdeWXF94hYj8pYh8ab+9mrhaE4nIR0XkOyISzAsqGmPOuWb/sQlW/h/gnwI7wJeA1/WOeQfwV2xyId4C/K+WbRrY1p8Dbtv/9wNTtVXbXu+4/8EmQPzOubYVuJVNNu1r9v/+sTnfW+A/Ab+//+9XAS8AOxO1918CbwC+GnnfPMZaWwwH6dTOuatAl07t6yCd2jn3FHCriLy6cbtCyrbVOfc559yL+38+xSZfYypp7i3AbwN/DnxnzMb1pGnru4FPOOe+BeCcm3t7HXDL/gLCl7EBwySltZxzn9m/fkzmMdYaDLFUaesxY8jajt9gQ+GplG2viNwJ/Arw6IjtCklzb18L3CYifyMifyci7xmtdcelae+fAD/NJpHvK8B/dM7NdYcf8xhrXY+hWjr1CFK3Q0R+gQ0Yfr5pi9LStPcPgfc553Y3P2yTSdPWU8AbgV8CzgL/U0Secs79fevGBaRp79uBLwK/CPwE8GkR+axz7getG1cg8xhrDYYlpVOr2iEiPwN8BHjAOfe9kdoWkqa99wGP70PhlcA7ROS6c+4vxmnigbT94LvOuYvARRH5DPB6YAowaNr7XuADbuPEnxeRbwI/BXx+nCaaZB9jjYMip4BngXs4DOL8s94xv8zRwMjnJwrgaNr6GjbZnT83RRut7e0d/zGmCz5q7u1PA3+9f+w54KvAP59xez8M/Jf9f/848G3glRP2h39CPPhoHmNNLQbXLp16qrb+DnAH8KH9X+HrbqKVdsr2zkKatjrnvi4inwS+DOyxqRQ2ybJ85b39PeBjIvIVNgPufc65SZZji8jHgbcCrxSRC8DvAqe9tprH2JoSvWrVqmNaMx9XrVp1TCsYVq1adUwrGFatWnVMKxhWrVp1TCsYVq1adUwrGFatWnVMKxhWrVp1TP8fKfvpM7cBWakAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(u(0.5, 0.5))\n",
    "plot(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Shape of a liquid drop\n",
    "\n",
    "The shape of a liquid drop resting on a flat surface is given by\n",
    "\\begin{align}\n",
    "\\nabla \\cdot \\left(\\frac{\\nabla h}{\\sqrt{1 + |\\nabla h|^2}}\\right) - B h + p = 0.\n",
    "\\end{align}\n",
    "\n",
    "Here, $h(x,y)$ is the height of the drop at $(x,y)$, $B$ is the Bond number, and $p$ is the liquid pressure. \n",
    "\n",
    "We impose $h = 0$ when $x \\in \\partial \\Omega$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Weak form\n",
    "The weak form is given by\n",
    "\\begin{align}\n",
    "\\int_{\\Omega} \\vec{n} \\cdot \\nabla v\\, dx - \\int_{\\Omega} B h v\\, d x + \\int_{\\Omega} p v\\, dx = 0 \\quad \\text{for all $v \\in V_h$}\n",
    "\\end{align}\n",
    "where\n",
    "\\begin{align}\n",
    "\\vec{n} = -\\frac{\\nabla h}{\\sqrt{1 + |\\nabla h|^2}}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# FEniCS implementation\n",
    "\n",
    "* The ```solve``` function can handle nonlinear problems (Newton's method)\n",
    "* The Jacobian can be computed **analytically** using the ```derivative``` function\n",
    "* The ```interpolate``` method can be used to create a non-zero initial guess"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# FEniCS code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fenics import *\n",
    "import matplotlib.pyplot as plt\n",
    "from timeit import default_timer as timer\n",
    "\n",
    "# define the Bond number and drop pressure\n",
    "B = Constant(1)\n",
    "p = Constant(1)\n",
    "\n",
    "# create a mesh for the unit square\n",
    "N = 30\n",
    "mesh = RectangleMesh(Point(0, 0), Point(1, 1), N, N)\n",
    "\n",
    "# choose elements\n",
    "P1 = FiniteElement('CG', 'triangle', 1)\n",
    "Vh = FunctionSpace(mesh, P1)\n",
    "\n",
    "# define boundaries and create h = 0 boundary condition\n",
    "def boundaries(x):\n",
    "    return near(x[0], 0) or near(x[0], 1) or near(x[1], 0) or near(x[1], 1)\n",
    "\n",
    "bc = DirichletBC(Vh, 0.0, boundaries)\n",
    "\n",
    "# create solution (h) and test function (v)\n",
    "h = Function(Vh)\n",
    "v = TestFunction(Vh)\n",
    "\n",
    "# weak form and Jacobian\n",
    "n = -grad(h) / sqrt(1 + dot(grad(h), grad(h)))\n",
    "F = dot(n, grad(v)) * dx + p * v * dx - B * h * v * dx\n",
    "J = derivative(F, h)\n",
    "\n",
    "# initial guess\n",
    "h_0 = Expression(\"x[0] * (1-x[0]) * x[1] * (1-x[1])\", degree=1)\n",
    "h.interpolate(h_0)\n",
    "\n",
    "# create function to solve the problem and time the calculation\n",
    "def solve_and_time():\n",
    "    start = timer()\n",
    "    solve(F == 0, h, bc, J = J)\n",
    "    end = timer()\n",
    "    print(\"Problem assembled and solved in\", (end - start) * 1000, 'ms')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Problem assembled and solved in 3.0049079996388173 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.tri.tricontour.TriContourSet at 0x7fe8905d3580>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAAD8CAYAAACVSwr3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2da6xtV3Xff+Oc+/I1Nn5BhGwj3Mh50CpI4ACK0pYEUWyqyoqEFEJUFJTIshpH/Qjqh6RSvhClH5IogHWFKIpaYUUNSl3JAaFUCUjEwaQiBoNAt0bFF0eywDyMj+89954z+2GfdbzO2vMxxlxzrjX3OesvHd2791577bnXmvO3xxhzzDHFOceiRYsW9bU1dwMWLVrUnhYwLFq0aE0LGBYtWrSmBQyLFi1a0wKGRYsWrWkBw6JFi9aUBIOIfFxEnhORrwZeFxH5ExG5KCJPisgbyzdz0aJFU0pjMXwCuDfy+n3A3Qd/DwAfHd+sRYsWzakkGJxznwOejxxyP/BnbqXHgZtE5DWlGrho0aLpdarAOW4Hnuk9vnTw3D8NDxSRB1hZFWyz/abz8soCH79o0aKQXnDPf9c59yrr+0qAQTzPefOsnXMXgAsAN27d6t566p0FPr5dbV1/fu4mLDJq/8WduZtQVJ+9+sn/l/O+EmC4BNzZe3wH8GyB8xbTMkAXaTVHX2kRRiXA8CjwkIg8ArwF+KFzbs2NqKW5Br2cv26Wz100r9zOS8XPmerDc4AjCQYR+STwNuA2EbkE/B5wGsA59zDwGPAu4CKwA7y/RkNLA2D2gX3duXk//yTqpcujTzGm3+RCxdf3a8MiCQbn3K8lXnfAbxdr0UA5QCg66DdgALvzZ+duQnHJzpXyJy19L42gifVLKzS2rj9fFQ4lXInissAgGwKVBvxxHKRzaI7raIaRpg8p4eHrxylYdOOkBiCaA0MKCmYQjADAFJ1z/9yZ6p9xErV1edf8Huv9VoEk1P8UwOj39RgkalgPzYEhJhUUMkAwFgBzDu796zbqFpq09dK17PeOuSdaqMT6TRIaw36aAIWcv65K4DOkjelVSSgogJALgNID/zgP5pIqfZ20oNHc7xQ8fH0tCot+/w1AohsDPkCUthqa6qEhNyIKhUJAGDP4aw30vXPbVc57XLR9ec90vPY+aQAS6i8xYAz7YRAUXZ+OAKI2HJoCg09BKCSAUBIGJQb+VIN879xmrqTfvrxvfk/uNU0BJXW/Y+AY9ikNKHIAUdu1aAYMPmvBC4UIEGIwsFgEFhCUGvCbOqBLqdT31wBGc89i8PD1jxAsNKDo91svJK47VyQHw6JmwDCUBQpjgKCFQA4Aagz2a2dPLkBOXdEMev31iUEkdL9DwBj2oxQoQpaEO39WDYeaVkOTYNBCIQSEsTCwQGDM4J9ikO+d861xm1fbl/P2MrFerxRIYvcuBI1h39CAwgeJGCBagEMTYMjJbrRCoQQMrBAoMfBbHNhjVeI7aeCSuv4xcPjutQ8WGlDEIBECRDD+kHArSgUgmwBDUgNroSQUUkDQwCAHALUG/LWzmwuSU1f0loT2+sUAErpvIWAM+0IKFDFIWABRJT08oebAkMpXsEAhBIQpYZALgE0e4LnK+c4pmISuvxUYPlj0+0kMElZAJBOsBlZDDXeiOTCsSZGnoIVCCAgpENSAQMmBv3cMl2dsK38kU9cxBA4rMPp9wAqJFCBScJjDamgfDD35rIUhFEoBIQUDLQRyADDVQJ/ic7QDfChL22KfEbr+WmD4QKGFhBYQWXCoPIXZNhgS1kIuFHKBkIKBBQJjB+WmWAlj26kBS+wzQu/33SsfLPr33AqJGCCKwKGimgJDdL36wFpIQcEChFwYaECQMzBKD/rWIGKxIjRtj53P934tLIag0EJCAwif9eCDw1xqCgxHFLEWSkEhBISpYJA7YFsb6FZZ258CSeh8ofdpYdG/1xZIxACRsh6GcIgGI3vuROkAZLtg6KlvLaSSl8ZAIQSEFAw0Hd06GEoN/v1GIbJV0GoYC4DhsTFIxNwNHyBKwKEvjTtRIpdhI8AQU99a0EDBAoTV8f7XUp21ptXQ6mC3yPIdUhDJBUDouBgkUlaExnqwwkE1hVlYbYJBWWwlnc2YD4UcIJS0HEoN/r1GC0RtG/p57FqEoKGFRcpS0ECiDwit9eCDg1ZTBCFnB0MqHVrrRgythVwo+ICQC4OaEGh1wGulbX8KIL5rp4XFcKBrIeEDRCgG0QdECg5ZVkOlacvZwZCrmAuRgsJYIIyFgRYEJQb//tm8BUu1tXVFN7UbuwYhaGhhEQOB5rX+86EYxNC9sMJhLjUDhlQqdN9aiK99sEGhBBBSMNCAwAqBVge8Vtr2xwDiu2YaWORCwve81r2wwOFIu2eKNTQDBp80VZj61kItKOQAIQUDLQhKAGD/TFsQ2drVJ4KFvn8IGMPr6gOFFhIaQGjdCwscNFaDL85QcsqyaTD4pC2sYoWCZfCHji0BAysIWhv0GmnanILH8DppQGGFhMVSCFkPFji0pPbAYCj/HrMWjh5nh0IpIKRgYAHBWAi4s/N2QrmiX54e+q4hYGhAoYWEBhDa51JwOPI5AathjozI9sDgkbWCcyzFOQcKViCUgEEOBOYe+Clp2peCx/C65IIiBokYIHKshxgcrFaDN85QYWZiI8CQUt9aiLkQKShorYQcIKRgoAVB0cFfGyQGC6FT6PuFgGEFRcyS0ADCYj3E4NBXHw7WWEMtNQsG7xLrRJbjUFNBIQSEGAyqgaAlq8HSlgREhtchFxT9exKyInyAyLEeYnCIuRQ++dyJmolOzYJBq5C1EEtzTkFhTiCoQVAYANtnyvwy7e1mltMPfZ8AAKyg8FkSISvCBwiL9RB7HFLIaphLzYPBF1+wlnLvWwsloGAFQgwGNUFQarDX+kwVRIbfOxMU/XsQsiJ8gMixHmIwsFoN3toNE+QzNA8GrbTWQl+loFAFCAYQzAGAEgq1OwoMIyh8lkTIivABwmI9hOAQA8WYQGRNbQwY/NWZNEVbw9bCkXMVhEI2EBQwGAuB02fmKQRydVff1XzfMQiL/jXzQKB/vUNWhA8QOdZDFAC912KByMPjPe7ElNOWGwOGTukKz7pVkf3HKSiMtRKmhMFcAz+lVLtS4Ohfg7GQ0AAiZD1Y4FDSpZhaTYLBul19KG8ht5ZCDhRKAmEOEJw/V8Zn3bmct/LL9z1CsBheHy8ouutrBETKegi5FrEVnRpQdPIXd5l+2lIFBhG5F/hjYBv4mHPuQ4PXXwn8N+C1B+f8L865/2pujSHr0afQ9KTFhbBCIeQ2BKEwAgi5ECg16Et+lgYgw++bAkUUELAGCR8gQtZDyrXow6GUS6GZnTicsiyc5JQEg4hsAx8G3gFcAp4QkUedc1/rHfbbwNecc/9ORF4FfENE/rtzTtVLQisr4/UX8rZ/i7kQMSiMshImgsGUABgrX1tTsOhfCx8kku7G2f2gBZGyHrSuxeFrRkvBnNdQeWZCYzG8GbjonHsaQEQeAe4H+mBwwA0iIsArgOeBYs5uF3j0xRd8bkTIWqgJhZJA0MBgDATOn5kWIDu7OvfCAgstJNYAEXAxQtaDJjD58mtlXIrY7MRUAUgNGG4Hnuk9vgS8ZXDMnwKPAs8CNwC/6pxb+2Yi8gDwAMA57BvZxjRmL8jYIqjSUBgDBCsMpgZASKF2aIDR/84pSJQChAUOuS5FJ80MBejjDKWWXmvA4Btxw2/yTuDLwC8DPwl8VkQ+75z70ZE3OXcBuABw49at2WHYHDciNQtx+FpqAVQmFEJAKAWDViBgka/NMVikIBGzIqKAyIBDXzkuhTel2uNOzJUFqQHDJeDO3uM7WFkGfb0f+JBzzgEXReRbwM8AXyzSSoO0QcdOFhciBwo5QNDAIBcEN5yZfudkgBd2dTNNw+8VAkV3jWJWRMiC0FgPKTj4Zis6OMRcCp98VsPcyU4aMDwB3C0idwHfAd4DvHdwzLeBtwOfF5GfAH4aeLpkQ31SbzYbsBaKQqGAlZACghUGc0HAJ19bNLDof2cfJGJWRMjF0FoPVjj0FXIptGsn5lYSDM65ayLyEPAZVtOVH3fOPSUiDx68/jDw+8AnROQrrFyPDzjnvlux3Ufkiy+M2bClFhRqAqEmBG48ffTcP7o64uL2ZIVFdx2sVkQMECXhEHMpfNK6E3NIlcfgnHsMeGzw3MO9/z8L/JuyTTsqzcKpVA1HrbXQlxUKFishBoRaMBgO8hzlnEMLk/53CkFCa0X4AFEDDn2FXIr0TIQuCDmVZs18jO0p0eUwrO9LmZe/oFXfWigBhdJAsMCgBARKydeWFCwskAgBQmM9eF0LAxxiLoVPVneiC0BOmQHZXEq0Nh06ukN1ouhKyFqI1VGoDYUYELQwqAmCV55eTYH98Gq8zL9FFhclBYkQIEZZDwk49BVyKVJWQ3w9hT8A2eUy1Exyag4MVo3JX4ipby3UhMJYIOTCoBvotd9nAUn/u2ggYQVEaTikpjFzVSLOMHZj240HQyftHhEaayFaZWkCKKSAYIVBLgRKyPfZGlhoIJECRC049BVyKVLTl5py81rVKPF2bMDQyVL2XaPYCkkNFKyuQwwKWiDMCQKNhu1LgaL73jFAWOAAR12LUM5DSDGXIiWrOzGXjh0YtNJYC5ppyb40UKgFhDEwuOl0vsnp0w+u2tLd+22PQSIGiJD1oHUthnDIdSnGTF+2pI0Cg2VGwlL2XaWEC5ELhTFAsMKgNACsn6MBhgYSuYAoDQeLtO5EK9ooMHTqZiS6wGOoIAsEXAurtTAxFGJAsMBgKhBoNWxPChQpSKQAUQUOPWmsBp8sQOhmJqZeM9EOGCJFWrRVocdsMqvVXFDQAqEkDG4+/WL09e9fvX7U+ftt1UAiZkHUhMMRGa0GqzvRBSDnzoBsBwwFlRNs1FoLfZWEwhgrIRcGqYE/9v0WcGggEculCFkPWjjElGM1+KSt1xDTVElOTYNBu7O1RdqEJp+sFZrHQqEkEMZCIEe+z9TAovteuYDIgYPJpVBYDdpsyFbjDHXzi2eW1nKoYS3UhMJNp3dUULj59IuHf63I0qbU9wxdI9819V374T0a3sNorYxe30jtEeJbizNmCh0CP5oja6b21bTFUEq++EJJayEHCrlASKkkBG7Zjp/r+b1xMYZ+W2OWRMyCCFkPpSyHvmIuRadUENJaq2EuNQkG67b3MJ7Aa+r9CsQCjlNAoRYQUgN/7Pst4OjanwJEzL3IhUNM1uSn46KN+8apqUpr2bZoQpNH1hLuNaFggcFYCORo+JkaUKQAkbIecuCQbTX0Yg2pIKRmdqKLN7QwM3FsYwzxqcsyFztlLWhXRfqgEPOvtT76LdsvHv61IEt7Ut8xdH1811ITc4jFG2rv7pWydrWVykpqY8CQW4dBm78QCjr23YhonUYFFHwdNASFkCxAaFlaSGgAMVQpOPTVv/dHXMtEEDL2I1TC/c1xuzXaOFciV6nqz1ZZS7lPAYUxMLjl1I+z39vX89deYf/sg3bHXI2bT78YdS+GroXWrYgp5lKM0SYEIE8MGHyKLq+mnLUwFgqlgVAKAtpza2GRAkQs/pALh1S8IVe5cYaYpkyL3hhXopR8pp02d6GTxVpoBQq3nPrx4d/Usn62xr3wSetWWBSKNYTcieOiYwGGzlcrPmWpkDXgWAoK2jjCXDAISdue1PcbA4chrC2xhpSscYZW1bwrEVtAlbvNPeS7ETFrQTs12ZcVCloYlNSrt1848vi5vRtGn7PfxpirEXMvQnGHWL5Dp5RLEVI/r0GT8JRSN0U5/HduNQ8Gn7R1HrsZCXXg0WgSxn5ZNC5EaSjkAmE48McebwVH1+4UIMbAIbY606d+rEEVhPSsn4gtqmo9AHksXIkxSiU1aeewc1wIn6aCwqu3Xzj8K63cc6e+Q+gaaBO9rC5FDWl/pKI1RpRlCMboRIFB6+uF1kX03Yixm8j6rIUcKFhiCDVhUOozS8KhZG0KS8JTyThDrSroKW0UGKpmgI2ILKdiCxoXIhcKGk0Ng7HtSMFujOVgsRo08Lcuxd8UNQUG7WYzY5UKPJbUmOmysVBoBQhDWQARfE05TTtpeTufpTBhXyuppsDQmkKzEf1fEutMhNZaGAOFUkC4ZXsn+jdWmnZa4aBxKVJWQ1/9ez3l+omhpl4v0eyshLV60zCXQbNGwucLljQNNdZCDShYlTvIQ+97fs9WPv7V2y9EZzJuOfXj4IyFb7Yilj7dKTZLYV2a7Zu2HLP3RAtqFgyl1EWBo8GfQplrVmvBWkMhBgXzDECBX3zNubWQ6NofAkQMDhqlchus6yhKqJuybCV3oa/NRdqBpsh21LoRJRSyFkpBoZQbUOvzYt8ldA20LkWuQu6EV5Efme7HKTRlOUfmbkgbDwartMGgsT5kaibC4kKUgMLUQBjz+Tlw0CgVa+hLA/2p4wxD1SiW3OnEgSGmMfGFWlvQl4JCKyoBB/95x1sNte6hdWYix3IoPaO3EWDwFWmJZYbNrZLWQkiq6b6ZrYSQtO0KfUeLSzFUrenL45bPoAKDiNwrIt8QkYsi8sHAMW8TkS+LyFMi8rdlm1lemvqOMG18YajQANBCoZRu3do/8ldKNeAwVMpq0LoTOZmu2j7WopJOiohsAx8G3gFcAp4QkUedc1/rHXMT8BHgXufct0Xk1bUaXFwFbl7MBNX8Qvl+6XJ96TFA0A760HHf27cboLds75inNw/f65mpCC220kozO1GrstNQXVHYOaS5k28GLjrnnnbO7QKPAPcPjnkv8Cnn3LcBnHPPaRsg5/Ur3qbSmKBSKndhbLQ8GpjLgEJJSyD3XKl2j03WGl7zku6Et69UtBRya59apfmU24Fneo8vHTzX108BN4vI34jIP4jI+3wnEpEHRORLIvKlq+5yVoPnWlRSSxZroWR6c2m3YOz5c+Hgu1bWeM3YKk/HUZr5Dt9IHNo3p4A3AW8HrgP+TkQed85988ibnLsAXAC4cevW4jbSMOtRldx0oFTwKBRfGOtGlJLFWqgJBN9naV2MMW7ForLS3LFLwJ29x3cAz3qO+bRz7kXn3HeBzwFvKNPEcpp7QcvQpC1lLWihUNtKKPG5se8yxmKyuBN92FsCkMdpZkIDhieAu0XkLhE5A7wHeHRwzP8E/qWInBKR88BbgK+Xber0spaIb9kknQMIU7WhhDthlbVvdBpmP2r3PZlaSTA4564BDwGfYTXY/9w595SIPCgiDx4c83Xg08CTwBeBjznnvlqv2SdLY62FFqDQSdOWWlZDTLWhnrJWo1srzhBXU+VUOuceAx4bPPfw4PEfAn9Yrmn1NMX8ciqpqcQU5aZBoZMm7mCNN2gWWWlWXS5aaSMyH1MqufjEMlVZK4W2rzG/kC1CodOYtulTwePuxBTB4U1NcjoWYMiWNvtxgozH0tZCy1DQas50bgv0515MVUMnGwwFVcNHbaEs281bZ71/JVQaXrU21hmbGr2JWsAwgTTxBYtqWgtaAJSCRKqtoe9aA5otzypNrQUMiw6VO8hLWhGlNYRwyQIuUWXGFlop1rKAAX9iSu48NUxcmXigXGuhxMDOPcdxiIccNy1gaEA+3zi4NqBCQK7kr30Ny6GmOzEnxFtWc2DYP1d/OetxVc4vb42BPKdbMeXO3j6r8rikRTcHhr5Se/S14o9tqmoOYOu5F3eiLTUNhpalneeeLNhlVKvBwkVtaAHDxKq9uKclLfDZXC1gKKCp5r9bLOy66HhqAYNBUxSAhfoZj5v4Sz4GirPlMmywNgYMU2/quWlqNXi3iRBatEFgOK6acnqtdeXAraZ1NcXq2Va1gOGEafkFj2tZL7HSAoZFixataQHDCdP390+ueazRD6+2t8/JHFrAMLNS5chOknJ2snpu74YKLVkptSPVcdYChmOinEE1hRYLZTPVZm/y6NSVNqfjaqjmryBs5mAdsxHNcC/LpSBsWhsDhhb0wq7ftJzKL112aVo0lRYwLFq0aE0LGDKlDUwNzdYxW7SX1BTuxCa6LItWWsCwqIqsUGg1eHpStdF3Y3v5QToi6+Cq9Ys+p6Uw5fTvzuX1amN7u/HiQpui5sCwdfnlclnbl49HmayULJ25dACy9CCuAYXQdy4xe/ODq0tA16fmwDCHfJT3/RpoVaKz5Xb6HJO81GDOPc/iRrSn5Y4sAlaDegwgWg00zpbDcCVvaLXiHi9gmEBTz0yM+QXuAKEZ6JZjQ0q1taYbMdSyTuJlLWAopBqdKtT5p0p08g38EjDoVNqFmCLwuLN7MrY3ONlgUJp7U3SG0p269KArBQOL5sz07OephDJeO13dPVW7OZPrZIMhQ/1OMsXquzFWQ8tBvTFt07oRKZdtihkJyYw1zK3NbPVAJQM2peg/7HSaOIPVathUOGjaZLUWNNduWTylV3u9ZgJtGsXHBtpagsNYKNRaeWqNEVmns7d2Jfp67Mdt+7IzfVYJqXqMiNwrIt8QkYsi8sHIcT8vInsi8u5yTdwczRHV1v6ytgCHWm3YpGI3W1dWgNje7R7P2JiIkndKRLaBDwP3Aa8Hfk1EXh847g+Az5RuZCmlqD3UmCQnjSzuROyXchPgoP3sUtZCKn8hFl8IBR5TQejjkg4NOovhzcBF59zTzrld4BHgfs9xvwP8BfBcwfaZ1JljHYVfpnIaCKmb2u8U2gBkKs5QUhY4TAkIy+ct9SbakeaO3Q4803t86eC5Q4nI7cCvAA/HTiQiD4jIl0TkS1fdZWtbj6VKWQ1W1QZE6fMHZ2c818qaQLYkNq1Lc+d8P7fDaMgfAR9wzkVXPTnnLjjn7nHO3XNazmnbeERzBGIsmrOT5fzilh7AueerGXC0uBFWeWexCgW3T11Z7+vbl6cpcaiZm7sE3Nl7fAfw7OCYe4BHRATgNuBdInLNOfeXRVpZU1e24OzRi3119xSnz1wrcvofXD3PTadf3nfx+1evX9s78fm969f2V3z+2iu8u1Q9t3dDdPelboBZ93oMDebQ7lClYDKqlmMBa2EoTW6KNva0abNffWla/gRwt4jcJSJngPcAj/YPcM7d5Zx7nXPudcD/AP7DRkBBoVAnmCLRaYxLUcpf7yyA4V8JadpocSF8SsV1YhaeJfA4t2Sn7PRG8g47564BD7Gabfg68OfOuadE5EERebBoayZUDs21nWPY2TRBSOsv3ZRwKK3n986PgkL4vOlrWCvbseSMRAsrLFVpfs65x4DHBs95A43Oud+wNMDtvIScj/vl25f32Tt3dCCfuuK4dtY2/ZjS3u4222faKg4Tcikg7VZAvmtRS1pYRadnK7gQYLf8tFmy1mlyrbZeKuPues9d7cyNSnuTrKnR1k5lsRpiZrN+3cC81oPWSgA7FEJKBR21bsRJ1MaDYWqzK5TPMFTKnSgpCxymBoT1M3OgULu+Rf+eJwOPERd1mPU4VAsuRKeNB4NW0SSnRLxBG4Wey2oAa1bg+SN/pZVz7uf2bihmKYB9inKOfSoPE/EaAkKnY7OQfPsK7PXu7dYV2D+7ovOeIaBcMs7ww6vX8crTLx0+Hk5dgn76EuLxBngZDqm4w9p5PQNYG5MoAZYU1GJQ8IFUk2Fa0o3wBR5LT1VOvUXjsQFDpyEgYpIrW7iz4QseymfY2T3D+TMre/CF3bPccEaPfB8cfMqFA+QD4ujn13c5VDMrBaAwxloIuRFzBx5rq1lXombENUe57oQmEzL0CxdzKzSmdco8n0vadlmhoFW17NQNTmgaaqO+SSlzykvxETfVanr6fsGscAC9390KICztyIGCxloYQmEIcmtS03FaUdlXU2Aonb2VK83NjnUajdVQEg6tA8IKhFpQGCOLG+GLL2hW+Po019qgpsBQW76bkwoS9TtBzJ1IWQ1a8zUHDmCL2ncDtRYk+uc3zZYkvoMFCj5ZrIVaCk1Vms4xwQ5tGxl83L7s2DunJ7B1ZkKrfhByqB9dPcuNp+MWkDYQ2SkUkDx8/WBgpYKTfcUGbip4WQosGqhZoWBJZvLJvDbCZyVEAo+pqUrfysoptZFg6JRKi+6mLL2v7Qr7ZwYXv7fSMjRtuXP5DOfP+WEwnKEYwmE4fQnhKUxgbRoT0nCAPED4VNvtGAME0EPBJ4u1EHIjaqyPaCWn4Vi4EsOLWfrixnxK6y+LNt4AcbdCE5Xv/PWWaiJa2lQKClYXYsxKypLxhTnVPBgs/pQGCDlxhqEssQbfHHkJOIAeEGALUtaQ5fNT32sMFKxSWQvHaJqy00a7ElUUcCdixVuGsYaUSxFSKObgy47sqxtEKRcD/Ob7WJdD8xmq9yUgF4OkFgo1rYW+fPGFsYHHqao3wQkCgy8A6Y0zKBWLNfikiTdAHA7gjzt0sgDiyPsSA3kIjtKWh8bqqQGFlGpVCW+1ZHxfG2MD5dJSexOOuBMB09ASa9BMfYXM3FjwTDM1Z3ExNOrHBUpCQdvOWlCwWAsaN8IaXygRC9u6XGD+06NjazHE1kxsXRH2z+oshdiiqqHVkONSxCwHINt6gKO/xFYropYswLICAfQxhRQUau8p0lcIEF1y09QLqOAYgsGyiArs7sTYQrEhOAAm1wL0gID1ATkVKHIsl5RVZIXC2LiC1VrIjS90gMjOYXip3JYMTYJh6/Iu++f8N+vUlX2und06THLqchk0QEglOh1ZbanIaQC71QDhYGSO9QA2QHQKDdgcYJRwWzQuUszFyoWCTzmrKH3yuRGl4gu1Fxk2CYbS8iU6WdyJoYZWwxRwgHSmZH9wWSDRV+1qSH1pU5lLAAH8ULC4ELnWgk+tJDKFdCLAkFLfndBaDSmXQgsHwOxaQNh66DQcdLmgKC3LNn2pDMbaUBhjLfRVYn3E1DrWYPC5F2PWTVhcCp9CRV1yrYdOmvUWc4EiZ7/OXCBAPShY05+12Y5Wy2GKBVTQOBi2XrrG/nX2JmoDkCF3ImQ1DGV1KSAPDuC3HjpprYi+QgM2FxglNuzVrHGoAYWU1qCQ6UZsQv5Cp6bBYFV8ijK8oAr0sxMpl0ILB0DtWoANEGCDRF81d+T2SVszwQoE0EOhlgsBejdiOCMx9x6tGwGG7ct77J0LmzLcM+AAABE4SURBVHKxVZYWd0ITa4A8OACjrQfQAQLWB1wuKErLUjxFk5NQGwo51kLKjbC4D10OgzXBb//Fcfe7HTC8dBmui++A3e1INZyyHCP17ETEpfDJF3MoZT3A0UGTggSEB2QtYIypnpQCQiy1OTQdWRMKGvnciJZnJtoBQ2H5LIW+O2G2GgbSzFJo4QBx6wHCgAA7JPqquRGORWOsA0jUUyjhPsQ2kglYC1Y3oiUdOzB4XYdEMDJmNYxxKSAMB/C7FrBuPcDRQaGFBNhBMZUsy6FTi58sVgLooBBbVp2Tt9CpBgRq1Eo9dmBIKRWEhEQgshAcIG49gB8QoLMiOvkG4NSwyKmJoFkJOdZKAAUUDApZC1o3opXAIxwjMPgCkClLoe9OqK0GUMEBWAtIAmrrAfSAAB0kOmkGqhYeJfdo0C6LtgIBRkChoLWwSWoODLJzBXdevwpKG4DsQyJkNfThMLQarHAAu/UAeYCA9UFlAYVP1TZlGchSIyEHCFAHCkNprAWflRBzLUKrKrt1ErWWXEODYNCqm5nwybLCMpYJWRMOsG49gA4QEIcE+AfcWFiMVc7GsalFT1YgQBkoWKYnh7JnO3ZTltNkPULDYOhWWHbZj10uQzdl6VPKnQgBI+ZSJBOflHAAigEC1geMZv9MzcDMhUep3aK1eztUAwKYoDCUxVroa+5y8UM1C4baipaWt8QbwAsHQGU9gA4QEIYE5IHCp6m3g7ds8pJKZbYAAcpAYUwFaK9rMbI4i9spE1hWZWmIyL0i8g0RuSgiH/S8/usi8uTB3xdE5A1FWqeUL4q7raD1EaJH3LVhZ/D6mp7nfB3v6u6pYEfduXwmuux3Z/fM4V9KL+yeXfubU772qGojKL5z7LpNCQWNtdBizoJPSYtBRLaBDwPvAC4BT4jIo865r/UO+xbwr51z3xeR+4ALwFtqNLgvX5whtQlNLP5gcSm6ThOLOUDceoB19wLiFsThMUpLoi8LHDQWRw3YaBc4pUqvjQEClINCSH1AtOZGgM6VeDNw0Tn3NICIPALcDxyCwTn3hd7xjwN3lGxkrkLxBc0Mxeq1dLzBG5AELyB8S7Y1gAA9JA6PV8IipKksDMtKR00dRhMQIAsKQ8WsTYu14LV8JywZ35cGDLcDz/QeXyJuDfwm8Fe+F0TkAeABgHOcZ//FHbaut6fk+gKQY9ZNxFKls+AAJusB4oAAPSQOjw8MuLHAyFXufg1jYNCpNBRicQXr0mqttTDljATowOC7Ct5vIyK/xAoMv+h73Tl3gZWbwY1bt5rtp9Qqy059d0JjNax9jgIOgM61gCxAgA4SnbR7XJTaUKWWLNWZY0CIZi8qgAA6KMRciJC1kJ6hCFsJU+QwgA4Ml4A7e4/vAJ4dHiQiPwd8DLjPOfe9Ms1bKVawJZbPkFLMpUjBAYzWAwQBAWlIQBgUMA4WcyinPLumToIVCKCDgs9KyIFCX31rYbQbUbBCNOjA8ARwt4jcBXwHeA/w3v4BIvJa4FPAv3fOfbNU42LVolPuhMZqWPu8gnAA1ICANCTABgrQDb5a8CixL4O2YEpybcMIIMA4KKy9r/CsRI0FVKAAg3Pumog8BHwG2AY+7px7SkQePHj9YeB3gVuBj4gIwDXn3D25jbKmRVsVA4UGDkDStQAFICAbEmAHhU9TbqwSk7VqUi4MIJzaXAIKa8cqrIWjz88TbBxKdTecc48Bjw2ee7j3/98Cfiu3EW7nJeT8dcpiLetxhr47obEa1s5phMPqGL/1AAZAQNSKgPUBYAFFX2M2ySmlMSXTVCseE8VTSgJhdcz6uWIuRGrHqaPP7ff+v7rntfeS6OvEZj4OYaCBA6StBxgJCIhWirKAoq+SdQxry7TsORMGEJ6GLAGF2GutWwuwQWDwBSD7cYYcq8EKB9BbDxBeZ9HvrElIgAkUh+1UAmNuZdU+UJRVKwUEyIOCNpbQQu0FnzYGDKWUAwcYbz3AugUBCkiACRSdtAOuFkByi514payvmKrDOBYIq2M9xyWgMCbL0Ze/UHuqEjYADL6ZiX6cwWo1gB0OYLMeYB0QYIMEGEDRyVCstlPRATxGxuKqoCvIGstaLA0E33NR96JnLfTdCO00Za0ZCdgAMJRSDhwgz3pYHduDQAYkwACKTprBlQGPIsoY+H1pqzKnN33RLZk++h7PsRlQSOUtrJ2vZy1MGXiEhsHgm7LU7EwVshpWr9ngEHouZD10GgMJ0IOiUxIYfY0coLVlLcsOuhJrpYAAOqvAAgWLtTCFGwENgyGlkDuxfpwdDqC3HiAMCW95+gQkQA+KTqnBZAJHZeUM/L4sdRaTG7+MBIL2uVhcochMxEHWY6laDLAhYIhlQHYKxRpWr9ngYHkO/C4GxK2I1fsGAFCC4vB4xZZ6MH4wziVrsVVN0ZRo3Y2RQPA9P4RCzIXw5S7A9G4EtAiGRJJT352ILaoqBQfwWw/D5yEMCEhDYvV+HSgOj08MHC045tKYKsvq3aQTlncsc1E7+EPPpaDQUt7CUO2BoSd7xeijLkUOHEAHiNjzMTcD1jurFhQvn1M34I9LeXNzsVWFG14KCKHnrVAIWQsh1ZyRgMbBEFLMaojFG1avx+EA+YDwvTbsgBpQQBgWq3OmB4oWHnNrTM1E0EFg9TmJ8xQChS+eYIHCUH03YqrAI2wQGDRxBp98y7J9cABGA6L/Wuh1DSgg3OFjwDj6OcfDWgD94O/LWl5t7GsaIKyOi0NBE1tIWQtjd7qGhsBwuJBKKYvV0N2MmFuxOi4PEL7XNK+DvwPHttDTDhItQOZWzqDvy1IxKVlaLQMUoRmHsVAYKmgtFK7D0KkZMBxRIAAZr8+Qdil8MQdgNCD6r2lejx0H8c6e2nfz8PPardGSJWvJNEiDQHNMfKeoclAYao6ZiL7aBENPsQDkMOEpBw6r4/IAAToIxJZ7D5XaQUs7QLQAmVs5A34o885OI2AA44CwOs73XNiF6FsLtYOOnZoHw1BDq0EDByDpWqyOjQNi9Z4wJGC8tRDrlNpt96DMgGtJuZWP1KscM2EA4dyEHCjkqmRyEzQABk2l6KHVYIXD6jm/9QB6QKzeE4YE1LUWLIPDApE5VarUWWmrAdIrIS1AWB2fhkLMWlhTpfgCNACGoBTVnGIKwQEwA2L1HjskoJy1kHpvzrk2SWO+j/a9miXRsczFMUCANBSmciOgMTDEZiY0VgOQtBxWz/tzHUKAWL1HD4lOFlgcnrug1WA57xSqCaucc48Fweoc4SBiKMCYA4U1VbQWoDEwWOWbpfC5FUDQeli95rcgIA2J1fv9AAh1vNQWehpZB/qmWhBj220vjJI+PpXKPAYI4IfClNYCtA6GgTvhXYqtgAOEAbF6ze9iwHon0IBida7w4E911hg4Dj+zcD8pYVHMBZ/cvR+1ZdU0axpi04+h4KIWCmsaWAv9wGOJ5CZoEAypRCcLHIAsQHTKBcXqXOFOl9pKz9LRNRDRqCWLouQmrzk1FceCYPV6JFkpkKMQgsLU1gI0CIY1eYKQITgAZkCAHxKrY8LuRqdQJ4rtjqXprNp9OFvcKXlKjSmmal3dqCm5lpp6HA2FyrGFTu2DAdRwgHB2ZAgQsH4zNdbE6rjw4I91Os2WerkdPndj3zlVs1Jy7tJmy/ZwmjyEWCZjS5ZCpybB4HUnjHCAdesBjt6gUJk4DShWx0XKumdCwyfL3pytliMvqZJ1DHK2mR8LAojHEoJA8FgLpRObOjUJhqACcADMgID1m6cFRafYztvaDhcDSKc5C3r4oNRygZFOOQP+6Ptt2YjatQ1ZQIBJoQANgyEYhAwkPmkA0UkLCgjDAtKdJwaOl88xrgNrwDJGc0Bg7DXRfUZ+GrJlgZNmlsEKBPBDodSMBDQCBk1a9BFFsiJjgOikBQXEO0GqYrW182lAsv4Z7f9611aJtQZD5a5u1BZTScYPIkHGmpZCpybAEFJ06rK7cAlAHDlfaJVm4GamCsNoO08KIJ1qdHCNcoDk01ztt6jUcuacakqqYGJi1iEEhZLWAjQEhpDV0F2IJCAgubZieGNS9SRTN19bUapUZ9QCxqqTNKDVnzeyjJp5RkExDTmFpdCpGTBA3KVIAgLWL64RFEc+T1GENrfz5JSog/mLd2yqatZKHDWlaMhJmBIK0BgYNOpfoGQpON+FV67YtNxwSyVrmLaop1a5sBqqxe9mVfH8gcykJC0MSrsR0CAYui+pCUaaINFJc5OMy71rJ6JYwZOjTRnQcyb9rKlwFqLFKqgBg76aA0Mn60xF6KJaCsweqsQNH1FLYqimBsNx0USpxUONdQlqA6FTs2AAm/UQUs6NyILJUDN1vEXTaEqffyoY9KUCg4jcC/wxsA18zDn3ocHrcvD6u4Ad4Decc/+nVCMtF2YMRDpNHegZoyIQm1mbdL1LaY7BblESDCKyDXwYeAdwCXhCRB51zn2td9h9wN0Hf28BPnrw7+SqdcFLAKeGTuKgmkKtD9za0lgMbwYuOueeBhCRR4D7gT4Y7gf+zDnngMdF5CYReY1z7p+Kt3gmnfSOsuhkSQOG24Fneo8vsW4N+I65HTgCBhF5AHjg4OGVz1795FdNrZ1XtwHfnbsRSm1SW2Gz2rtJbQX46Zw3acDgW+A/XNurOQbn3AXgAoCIfMk5d4/i85vQJrV3k9oKm9XeTWorrNqb8z7N0rxLwJ29x3cAz2Ycs2jRog2RBgxPAHeLyF0icgZ4D/Do4JhHgffJSm8Ffnic4guLFp00JV0J59w1EXkI+Ayr6cqPO+eeEpEHD15/GHiM1VTlRVbTle9XfPaF7FbPo01q7ya1FTarvZvUVshsr6wmEhYtWrToZdUt/7No0aKN1AKGRYsWrak6GETkXhH5hohcFJEPel4XEfmTg9efFJE31m5TSIq2/vpBG58UkS+IyBvmaGevPdH29o77eRHZE5F3T9m+QRuSbRWRt4nIl0XkKRH526nbOGhLqi+8UkT+l4j840F7NXG1KhKRj4vIcyLizQvKGmPOuWp/rIKV/xf4Z8AZ4B+B1w+OeRfwV6xyId4K/H3NNo1s6y8ANx/8/7652qptb++4/80qQPzuVtsK3MQqm/a1B49f3fK1Bf4T8AcH/38V8DxwZqb2/ivgjcBXA6+bx1hti+Ewndo5twt06dR9HaZTO+ceB24SkddUbpdPybY6577gnPv+wcPHWeVrzCXNtQX4HeAvgOembNxAmra+F/iUc+7bAM651tvrgBsOFhC+ghUYZimx5Zz73MHnh2QeY7XBEEqVth4zhazt+E1WFJ5LyfaKyO3ArwAPT9gunzTX9qeAm0Xkb0TkH0TkfZO1bl2a9v4p8LOsEvm+AvxH51yrJbvNY6x2PYZi6dQTSN0OEfklVmD4xaotikvT3j8CPuCc21v9sM0mTVtPAW8C3g5cB/ydiDzunPtm7cZ5pGnvO4EvA78M/CTwWRH5vHPuR7UblyHzGKsNhk1Kp1a1Q0R+DvgYcJ9z7nsTtc0nTXvvAR45gMJtwLtE5Jpz7i+naeKhtP3gu865F4EXReRzwBuAOcCgae/7gQ+5lRN/UUS+BfwM8MVpmmiSfYxVDoqcAp4G7uLlIM4/HxzzbzkaGPniTAEcTVtfyyq78xfmaKO1vYPjP8F8wUfNtf1Z4K8Pjj0PfBX4Fw2396PAfz74/08A3wFum7E/vI5w8NE8xqpaDK5eOvVcbf1d4FbgIwe/wtfcTCvtlO1tQpq2Oue+LiKfBp4E9llVCptlWb7y2v4+8AkR+QqrAfcB59wsy7FF5JPA24DbROQS8HvA6V5bzWNsSYletGjRmpbMx0WLFq1pAcOiRYvWtIBh0aJFa1rAsGjRojUtYFi0aNGaFjAsWrRoTQsYFi1atKb/DzoJfPghRbP0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "solve_and_time()\n",
    "plot(h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# If you want to customise the nonlinear solver...\n",
    "* Use the ```NonlinearVariationalProblem``` class\n",
    "* Can invoke the PETSC SNES solver\n",
    "* Can choose the linear solver (e.g. MUMPS, GMRES, preconditioners, etc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "```python\n",
    "problem = NonlinearVariationalProblem(F, X, bcs = bcs, J = J)\n",
    "solver = NonlinearVariationalSolver(problem)\n",
    "\n",
    "solver.parameters['nonlinear_solver'] = 'snes'\n",
    "solver.parameters['snes_solver']['absolute_tolerance'] = 1e-8\n",
    "solver.parameters['snes_solver']['linear_solver'] = 'mumps'\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Automatic derivation of the weak form\n",
    "\n",
    "FEniCS can derive the weak form if the PDE can be \n",
    "obtained from minimisation of a functional\n",
    "\n",
    "The total energy of the drop is\n",
    "\\begin{align}\n",
    "E[h,\\nabla h] = \\int_{\\Omega} \\left(\\sqrt{1 + |\\nabla h|^2} - p h + \\frac{B}{2}h^2\\right)\\, dx\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The weak form can be obtained using the ```derivative``` function:\n",
    "```python\n",
    "E = (sqrt(1 + dot(grad(h), grad(h))) - p * h + B * h**2 / 2) * dx\n",
    "F = derivative(E, h)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# What happens if $B \\gg 1$?\n",
    "\n",
    "* The Bond number is defined as $B = \\rho g H / \\gamma$\n",
    "* $B \\gg 1$ corresponds to a situation where gravity is very strong"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Problem assembled and solved in 5.312493000019458 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fe8904010d0>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVEAAAD8CAYAAADOg5fGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3df5BdZZ3n8ffn3v4Rwo8kEGEVnAHXOA7OCAMRqBnHdaRUws4QXGEELflRWClmYdzdqq0VdsqZqR3/wJrZ0WFV2Iho4uzwo9AdYolSGMcfOxoVFAOISkQGIwyQAB3SMX27+373j/Pczsnte849nds/bnc+r6pbfe855zn3uZ3bnzznec45jyICMzM7NLWFroCZ2WLmEDUz64FD1MysBw5RM7MeOETNzHrgEDUz60HXEJV0q6RnJT1csF6SbpS0Q9J2SWfMfjXNbCmQdJ6kn6S8uK7D+sI8KSor6WJJj0hqSlrbtr/r0/Y/kfT23PIzJT2U1t0oSWn5sKQ70vLvSDq522eq0hL9DHBeyfp1wJr02ADcVGGfZnaYkVQHPk6WGacCl0o6tW2zjnnSpezDwH8AvtH2fqcClwCvI8uwT6T9kPa7IfderYy7CnghIl4NfAT4cLfP1TVEI+IbwPMlm6wHNkdmG7BS0su77dfMDjtnATsi4vGIaAC3k+VHXlGeFJaNiEcj4icd3m89cHtEjEXEz4EdwFlpf8dExLcju9poM3Bhrsym9Pwu4NxWK7XIwIx+BZ2dCPwi93pnWvZ0+4aSNpClP3XqZy7Xio47LK1zzd24ZlOazdLVe5q7d0XEyw5192/9gyNi9/Pl79Hyg+2NR4D9uUUbI2Jj7nWnrDi7bTdFeVKlbLsTgW0d9jWenrcvP+j9I2JC0ghwHLCr6E1mI0Q7JV7Ha0nTL3QjwDG14+KcgbdP26Z+3CpYeUzhmzWXDR1aLYFYNhsf12z2aP9ET+Vr+xuF6+KpZ7h3z6f/pZf9736+yT9/+RWVtl3+iif2R8Takk2qZEXRNpVzpsd9zfh9ZiNVdgKvzL0+CXjqUHZUP24VzX+zmomVyzqunxyu0Rx2S9QOL7Wxzi3B+lgTWEa9QxDXX9w3x7U6JFWyomiboQplq77fzvS8075aZXZKGgBWUN6dOSshugW4VtLtZM3rkYiYdijfTStAx05YPi0oJ4ZrTA5n/0FMDhfvY3KotOvCbNGoN/KNn/rB68ZaP4OBsSbjxxz8Z7zs2f30qe8BaySdAvySbNDn3W3bdMwTSc9VKNtuC/APkv4WeAXZANJ3I2JS0kuSzgG+A1wG/K9cmcuBbwMXAV+NLndp6hqikm4D3gyslrQT+AtgECAibgbuAc4n67TdB1zZbZ/t8gH6q5cNHhSUrWBsDsFkhyP5ZsWj+7LwNZtPrRAsl33vax2O1uuNbHl9GBq5gK2PwfCeydmp5BxIfYzXAveS/c9wa0Q8IunqtL4wT4rKAkh6B1kIvgz4oqQHI+Ltad93Aj8CJoBrIqL1C/oTsjOPjgC+lB4AnwI+K2kHWQv0km6fSwt1K7xWn2h7gDaOFuNHHdiuOXQgAJuD0BzqXN+i5WaLTa2hwte18exnfSwFaS5kB/cGy5+bZNmz+6m/uK/VJ/pAl37KUmecNhwz6BPt6b0WqwUdaekUoGPHZqGZD8zmUBBDTRg+0DdUGzz4f1wBA0P9+7+wWScTjfq0ZQE0xw8sn/pWj9VQo0atIZqDqaWaQnXZ7rmp3yTBSLNS0/mwtWAhqnodVh7DxMplNI4ZYHIYxo/KHhNHZWHZCs7a4ORUSA4PHdyJvmxwvOt7HTlUPIJpNp9GG239T0ce/HL/+OC0MmON7M90YrBOc7zO5FiNZiMbN6g1xMDowS3XyWUDTI9mmyt9c85P+6BQPkDhQIC2h2ZRQB456OC0/tP6Xo6Od+7MP3KoMT1ocyaAJhDDTRir0aRGsyEmh7KBionhGoNkpwLWlh8Be2b9I1ibvglRyPo/84funQK0PTTzYXmUg9MWifx3dW9boB452CgMWciCFLIwJbVIs4FXMTmcnQrYV3/YS1zf/a6bQ1EaoO0tzPyXccXgr+a1rmazYcXgrxgZP2JGZVp9qTHUpDkkJodV+UwVm119E6KtUfgYKg/QotA8ZnDm58atGhjtveJmBV6YOLL7Rskxg/vZM37gIpN8sB412JjWWoUDh/bNRg1Gs0N6n8o3/xY8RCdzJ9I3c33qRQFaFJzdAnFFvS+v4LAlrNN3bmRyeeH2qwZGZxS8kIJ0qEZzcG6u5JsMGKl26fxha0FDtHUd/ORwdiJ961C+LEBbwZkPzaKAdHBav2l9J8vCdKYawx6LX0gLGqKxbGDqEs9Wf05tsDxAW+GZD8hOYXncwN65/wBmC6Bs0Mnm34IfzsPMr3kvCtBegnNlza1W682Lzeqty/bv6u6J7DK9FfV9U63U/PO89sHVxuhQujAldYv5Jj3zqi9CtIpWK/RQAtQBafOhyvesKGiPG9g7oyDtxKPzC6OvQrQ5FAhKB5NaygLUoWn9qv27mQ/VsiDN6zRSP1cmELubnW9NaZm+CtGqigK0SngeV+vb24TZElUWQvnv7IvN5VPf590TRxW2Qlun+bWubMpOCcxObfHtIOdf34Vo0U1EWofyMw1Qh6YttE7fwU7BurK2b9rhftnh/JFDjalr7fN9oja/+iZEy04S7nQifVGAHmpornBfvPVoJudT5r+n+UBtBWn7oX1L/oR8yLq+RvEZ9gtpwUN0omQksag/9FAC1CFpc63Kd6xT0La+t60w7RSkLccM7p/xJaI2t/o2WmZ6F6aiAF1Rc4Ba/5jpd7HTWSe+R0R/WfCWaF4Mzez6sm4DSQ5P60crasUt0qoj4Ud1udPTbJmM2ozOfz0c9WXMFN0jdEV9X+UT6h2g1s+6fT9ncppep8HYfp0eXNJ5kn4iaYek6zqsl6Qb0/rtks7oVlbSsZLuk/RY+rkqLR+UtEnSQ5IelXR9Wn60pAdzj12SPprWXSHpudy693X7TP35m56h9kP52QjQFTV31lu5XqfNKGqRLlWS6sDHgbeSTU38PUlbIuJHuc3Wkc3KuYZsts+bgLO7lL0O2BoRN6RwvQ74AHAxMBwRvy1pOfAjSbdFxBPA6bl6PQB8PleHOyLi2qqfa9GEaNWbicw0QB2WdqjKvjtVA7Y9SDsd0reOvtpPdTpysPwu+H3oLGBHRDwOkKZFXk82G2fLemBzmqZ4m6SVkl4OnFxSdj3ZjMQAm4CvkYVoAEem+eOPABq03etf0hrgeOCbh/qhFu1Bb6fDnaoBuqI2PPUwmwsz+Y7NVtfTIriX6InAL3Kvd6ZlVbYpK3tCRDwNkH4en5bfBYwCTwNPAn8TEc+3vd+lZC3P/HTB70xdCXdJemW3D7UoWqKzcfNkB6YtlNZ3r6x12unQvv3k+xX1fdPuN3rkUIOR0bk75WmC+rTTrEqslnR/7vXGiNiYe93paoD2uc6LtqlStt1ZZJOlvgJYBXxT0ldardnkEuC9uddfAG6LiDFJV5O1bN9S9iaLIkTL+IokWyxW1IaX+vTDu7rMO78TyLfsTgKeqrjNUEnZZyS9PCKeTof+z6bl7wa+HBHjwLOS/hlYC7S6BE4DBiLigdZOIyI/+fQngQ+XfB6gDw/n26dEng1uhZr1he8BaySdImmIrBW4pW2bLcBlaZT+HGAkHaKXld0CXJ6eXw7cnZ4/Cbwl7etI4Bzgx7n3uhS4Lf/mKYRbLgAe7fahFn1L1GwxqdIancn5ootJRExIuha4F6gDt0bEI+mwmYi4GbgHOB/YAewDriwrm3Z9A3CnpKvIgvPitPzjwKeBh8m6Az4dEdtzVfrj9F5575d0AdnMK88DV3T7XA5RM5s3EXEPWVDml92cex7ANVXLpuW7gXM7LN/LgUDttL9XdVh2PXB98SeYziFqZoUmozar80EtRX3XJ9qu0x2cZnJ6k/tDrd/08p30NN/9py9C1NMamNli1Rchama2WDlEzfpEpy4pzxfW/xyiZmY9WNKj8x5UMuuNR+e7c0vUzKwHlUK0wo1UV0j6gqQfSnpE0pWzX1WzpcNHSUtH1xDN3Qx1HXAqcKmkU9s2uwb4UUScRnZfv/+Zrm81M1vSqrREp26kGhENoHUz1LwAjpYk4Ciya05n/04iZmZ9psrAUqeboZ7dts3HyO6k8hRwNPCuiJg28YGkDcAGgGW1yvcoNLMFMklt2j1M7WBVWqJVbob6duBBspufng58TNIx0wpFbIyItRGxdqjmubPNbPGrEqJVbqR6JfD5yOwAfg68dnaqaGbWv6qEaJUbqT5JuhWVpBOA3yDdPdrMbCnr2ida8UaqfwV8RtJDZIf/H4iIXXNYbzOzvlDpiqUKN1J9Cnjb7FbNzKz/LenLPs2sN5NRY8/40puqZDb5sk8zsx44RM1s3lS4hFySbkzrt0s6o1tZScdKuk/SY+nnqrT8PZIezD2akk5P676W9tVad3xaPizpjvQe35F0crfP5BA1s3lR8RLydcCa9NgA3FSh7HXA1ohYA2xNr4mI/xMRp0fE6cB7gSci4sHce72ntT4iWnPVXwW8EBGvBj7CYpx33syWrCqXkK8HNqdzzrcBK9Nc8GVl1wOb0vNNwIUd3nvaHPMF8vu6Czg3Xc5eyANLZlZoMsTIeOWrC1dLuj/3emNEbMy9rnIJeadtTuxS9oSIeBogIp5uHZq3eRfTA/vTkiaBzwEfStM1T71POr1zBDgOKDxl0yFqZrNlV0SsLVlf5RLyom2qlO38ptLZwL6IeDi3+D0R8UtJR5OF6HuBzYfyPj6cN7P5UuUS8qJtyso+kw75ST+f5WCX0HYoHxG/TD9fAv6BrLvgoPeXNACsILsrXSGHqJnNlyqXkG8BLkuj9OcAI+lQvazsFuDy9Pxy4O7WziTVgIvJ+lBbywYkrU7PB4E/BB7usK+LgK+mw/xCPpw3s3lR8RLye4DzgR3APrKbGxWWTbu+AbhT0lVk9/G4OPe2bwJ2RkT+Xh7DwL0pQOvAV4BPpnWfAj4raQdZC/SSbp/LIWpm86bCJeRBNlNGpbJp+W7SDZA6rPsacE7bslHgzILt93NwCHflEDWzQpNRY++4Z/op4z5RM7MeOETNzHrgEDUz64FD1MysBx5YMrNCzRCjHlgq5ZaomVkPHKJmZj1wiJqZ9cAhambWA4eomVkPPDpvZoWaTTHa8Oh8GbdEzcx64BA1M+uBQ9TMrAcOUTOzHnhgycwKNUPsHx9c6Gr0NbdEzcx64BA1s3kj6TxJP5G0Q9J1HdZL0o1p/XZJZ3QrK+lYSfdJeiz9XJVb93pJ35b0iKSHJC2TtFzSFyX9OC2/Ibf9FZKek/Rgeryv22dyiJrZvJBUBz4OrANOBS6VdGrbZuuANemxAbipQtnrgK0RsQbYml63pjz+e+DqiHgd8GZgPJX5m4h4LfA7wO9JWperwx0RcXp63NLtczlEzWy+nAXsiIjHI6JBNo3x+rZt1gObI7MNWJnmki8rux7YlJ5vAi5Mz98GbI+IH0I2oV1ETEbEvoj4p7SsAXyfbB77Q+IQNbNCEWKsMVDpAayWdH/usaFtdycCv8i93pmWVdmmrOwJaW560s/j0/LXACHpXknfl/Tf2j+fpJXAH5G1YFvemboS7pL0ytJfEB6dN7PZsysi1pasV4dlUXGbKmXbDQBvBN5ANof9VkkPRMRWmDrcvw24MTcv/ReA2yJiTNLVZC3bt5S9SaWWaLfO4LTNm1NH7COSvl5lv2Z2WNkJ5Ft2JwFPVdymrOwz6ZCf9PPZ3L6+HhG7ImIf2Zz1Z+T2sRF4LCI+2lqQDvnH0stPUjA/fV7XEK3SGZyaxJ8ALkgduBd326+ZHXa+B6yRdIqkIeASYEvbNluAy9Io/TnASDpELyu7Bbg8Pb8cuDs9vxd4fRqNHwD+HfAjAEkfAlYA/zn/5q0wTi4AHu32oaoczk916KY3aXXo/ii3zbuBz0fEkwAR8ey0vZjZYS0iJiRdSxZudeDWiHgkHTYTETeTtRbPB3aQHYJfWVY27foG4E5JVwFPkhpxEfGCpL8lC+AA7omIL0o6Cfgz4MfA9yUBfCyNxL9f0gXABPA8cEW3z1UlRDt16J7dts1rgEFJXwOOBv4uIja37yh1NG8AWFY7qsJbm9lSEhH3kAVlftnNuecBXFO1bFq+Gzi3oMzfk53mlF+2k859rETE9cD1pR+iTZUQrdKhO0DWd3AucATwbUnbIuKnbRXcSNYPwYrB47t1CpvZAouAiUZ9oavR16qEaNXO4F0RMQqMSvoGcBrwU8zMlrAqo/NVOoPvBn5f0oCk5WSH+107ZM3MFruuLdEqncER8aikLwPbgSZwS0Q8PJcVNzPrB5VOtu/WGZxe/zXw17NXNTOz/ucrlsysWIjmuAeWyvjaeTOzHjhEzcx64BA1M+uBQ9TMrAcOUTOzHnh03syKNYExt7XK+LdjZtYDh6iZWQ8comZmPXCImpn1wANLZlYshBpua5Xxb8fMrAcOUTObN91mDk4T1N2Y1m+XdEa3spKOlXSfpMfSz1Vp+cmSfpVmIX5Q0s25MmdKeijt60aliZYkDUu6Iy3/jqSTu30mh6iZzYsqMwendWvSYwNwU4Wy1wFbI2INsDW9bvlZRJyeHlfnlt+U9t96r/PS8quAFyLi1cBHgA93+1wOUTObL1MzB0dEA2jNHJy3HtgcmW3AyjSNcVnZ9cCm9HwTcGFZJdL+jomIb6eJ8TbnyuT3dRdwbquVWsQhamazZbWk+3OPDW3rO80cfGLFbcrKnpDmpif9PD633SmSfiDp65J+P/ceOwv2NfU+ETEBjADHlXxmj86bWYmAWqO0IZa3KyLWlqyvMnNw0TZVyrZ7Gvi1iNgt6UzgHyW9rsu+Zvw+boma2XypOnNwp23Kyj6TDtFbh+rPAkTEWJqTnoh4APgZ8Jq0r5MK9jX1PpIGgBXA82UfyiFqZvOlyszBW4DL0ij9OcBIOkQvK7sFuDw9v5xs9mEkvSwNSCHpVWQDSI+n/b0k6ZzU33lZq0zbvi4Cvpr6TQv5cN7M5kWVmYPJJsQ8H9gB7AOuLCubdn0DcKekq4AngYvT8jcB/0PSBDAJXB0RrVblnwCfAY4AvpQeAJ8CPitpB1kL9JJun8shambzptvMwanVd03Vsmn5buDcDss/B3yuYF/3A7/VYfl+DoRwJQ5RMyukmQ0sHZbcJ2pm1gOHqJlZDxyiZmY9cIiamfXAA0tmViygNr7QlehvbomamfXAIWpm1gOHqJlZDxyiZmY9cIiamfXAo/NmVkhNqI8tdC36m1uiZmY9qBSi3Wboy233BkmTki6avSqamfWvriFacYa+1nYfJrvfn5nZYaFKS7TKDH0Af0p2775nZ7F+ZmZ9rcrAUqdZ9s7ObyDpROAdwFuANxTtKM3+twFgWe2omdbVzOZbQK2x0JXob1VaolVmv/so8IGImCzbUURsjIi1EbF2qHZE1TqamfWtKiFaZYa+tcDtkp4gm9zpE5IunJUamtmS0W2QOk1Qd2Nav13SGd3KSjpW0n2SHks/V6XlZ0l6MD1+KOkdafnRueUPStol6aNp3RWSnsute1+3z1QlRLvO0BcRp0TEyRFxMnAX8B8j4h8r7NvMDhMVB6nXkc3KuYas6++mCmWvA7ZGxBpga3oN8DCwNiJOB84D/rekgYh4KSJObz2AfwE+n6vDHbn1t3T7XF1DNCImgNYse48Cd7Zm6GvN0mdmVkGVQer1wObIbANWprnky8quBzal55uACwEiYl/KL4BlTO+GRNIa4Hjgm4f6oSpdsdRthr625VccamXMbFFbLen+3OuNEbEx97rrIHXBNid2KXtCmkueiHha0vGtjSSdDdwK/Drw3lyotlxK1vLMB+w7Jb0J+CnwXyLiF5TwZZ9mVkgB9eqj87siYm3Z7josa28dFm1Tpez0DSK+A7xO0m8CmyR9KU2L3HIJ8N7c6y8At0XEWDrS3kR21lEhX/ZpZvOlyiB10TZlZZ9Jh/ykn9POVY+IR4FRcnPNSzoNGIiIB3Lb7Y6I1t0CPgmc2e1DOUTNbL50HaROry9Lo/TnACPpUL2s7Bbg8vT8cuBugLTtQHr+68BvAE/k3utS4Lb8m7fCOLmAbByolA/nzWxeRMSEpNYgdR24tTVIndbfTDb2cj6wA9gHXFlWNu36BuBOSVcBTwIXp+VvBK6TNA40yc4a2pWr0h+n98p7v6QLgAngeeCKbp/LIWpm86bbIHUa4Lmmatm0fDdwboflnwU+W1KXV3VYdj1wffEnmM4hambFfNlnV+4TNTPrgUPUzKwHDlEzsx44RM3MeuAQNTPrgUfnzaxQdtln16srD2tuiZqZ9cAhambWA4eomVkPHKJmZj3wwJKZFWtCfaz7Zoczt0TNzHrgEDUz64FD1MysBw5RM7MeOETNzHrg0XkzK6SA+pgv+yzjlqiZzRtJ50n6iaQdkq7rsF6Sbkzrt0s6o1tZScdKuk/SY+nnqrT8rZIekPRQ+vmWXJmvpX09mB7Hp+XDku5I7/EdSSd3+0wOUTObF5LqwMeBdcCpwKWSTm3bbB2wJj02ADdVKHsdsDUi1gBb02uAXcAfRcRvk80C2j7f0nsi4vT0aE2zfBXwQkS8GvgI8OFun8shambz5SxgR0Q8HhEN4HZgfds264HNkdkGrEzTGJeVXQ9sSs83ARcCRMQPIqI1N/0jwDJJw13qmN/XXcC5klRWwCFqZrNltaT7c48NbetPBH6Re70zLauyTVnZE9Lc9KSfx3eo2zuBH0RE/vqrT6dD+Q/mgnLqfSJiAhgBjiv8xHhgycxKqBkMjDWrbr4rItaW7a7DsvZRq6JtqpTt/KbS68gOy9+WW/yeiPilpKOBzwHvBTYfyvu4JWpm82Un8Mrc65OApypuU1b2mXTIT/rZ6t9E0knA/wUui4iftZZHxC/Tz5eAfyDrLjjo/SUNACuA58s+lEPUzObL94A1kk6RNARcAmxp22YLcFkapT8HGEmH6GVlt5ANHJF+3g0gaSXwReD6iPjn1htIGpC0Oj0fBP4QeLjDvi4CvhoRpS1RH86b2byIiAlJ1wL3AnXg1oh4RNLVaf3NwD3A+cAOYB9wZVnZtOsbgDslXQU8CVycll8LvBr4oKQPpmVvA0aBe1OA1oGvAJ9M6z8FfFbSDrIW6CXdPpdD1MzmTUTcQxaU+WU3554HcE3Vsmn5buDcDss/BHyooCpnFrzHfg6EcCUOUTMrFlCrPrB0WHKfqJlZDyqFaIVLtd6TLtHaLulbkk6b/aqamfWfriFa8VKtnwP/LiJeD/wVsHG2K2pm1o+qtES7XqoVEd+KiBfSy21k53CZmS15VUK0yqVaeVcBX+q0QtKG1iVhjeavqtfSzKxPVRmdr3wZlKQ/IAvRN3ZaHxEbSYf6KwaP900KzfqcmlD36HypKiFa5VItJL0euAVYl87bMjNb8qoczne9VEvSrwGfB94bET+d/WqamfWnri3Ripdq/TnZ7aI+ke4oNdHlbi5mh7WR5lj3jWxRqHTFUoVLtd4HvG92q9a7keYYK2rd7sFqZnbofNmnmRWLoL5/YqFr0dd82aeZWQ8comZ9YqTDmUQvNpfPf0VsRpZMiHb6ApqZzbW+D9E948umLZvJ/84eBbV+4+/k0tL3IWpmB7wwceRCV8Ha9MXofK1x4PlYYwD8PTHrC2oG8uh8qcOiJerDJ+sXh/t3scK9iSXpxrR+u6QzupWVdKyk+yQ9ln6uSsuPk/RPkvZK+lhu++WSvijpx5IekXRDbt0Vkp5L89E/KKnr+e+HRYiaLSa7m9PHAZaCivcmXgesSY8NwE0Vyl4HbI2INcDW9BpgP/BB4L92qM7fRMRrgd8Bfk/Suty6OyLi9PS4pdvnWvQhulS/cGa96NND8K73Jk6vN0dmG7AyzSVfVnY9sCk93wRcCBARoxHx/8jCdEpE7IuIf0rPG8D36eEeyIsiRKt2pped5jTSHJt6mM23bt+9KqfojUyWn5VS7/+vdpV7ExdtU1b2hDQ3Penn8VUrlOam/yOyFmzLO1NXwl2SXllQdEpfDCwdiheby1lZ2zdt+UgTVnT5ryH/Zfa19TYXZuM/6yqn8o02hnp+n1IR1PY3um+XWS3p/tzrjekewi1V7k1ctE3l+xpXJWkAuA24MSIeT4u/ANwWEWPpJkubgLeU7acvQ3S0McSRgwf/w41MLmdFfXpo9qrTl93BajPRa2DmW6GLvHtqV5e7t1W5N3HRNkMlZZ+R9PKIeDod+j9bsb4bgcci4qOtBW33Qv4k8OFuO+nLEC2ze+IojhvYe/Cy5jKOqx3o9qjSGi3jQ36bL1WvtOt2KL9ITN2bGPgl2b2J3922zRbgWkm3A2cDIykcnyspuwW4HLgh/by7W0UkfQhYQdvd51phnF5eADzabV99FaJqlCffTFqjvQap2VyrGqC7J47quHx0fI4P5WdZxXsT3wOcD+wA9gFXlpVNu74BuFPSVcCTwMWt95T0BHAMMCTpQuBtwB7gz4AfA99P90D+WBqJf7+kC4AJ4Hngim6fq69CNG90fIijBsv7Yor6RVtaX1KHqfWbbgG6VG88UuHexAFcU7VsWr4bOLegzMkFVenUx0pEXA9cX1Cmo0UZL53+Zy7rS/LNSayfFH0fF3l/6GGrb1uiAHvT4coxg/sLt8m3Rtv7RvPKgtQtVZttM/2POx+g+VZopwbDyPgRU38bABON+tTzemOWJ9FtNmG/xwjK9HWItuwZX8aqgVHgQL9ofoCpPUjzikI1zy1VWwidWp5lAfrCxJEd72oGUGt0PDq1edB3ITrRqLN/aJAjhw7uD22dcL+ivq9rkOblv6hVAtVsrnQ7XC8K0JHJ5UtldH5J6psQbb/aotO5op20B2lLt0A16xftg0jtAdpu7/gQo+NDjDaG2D8+OLW8VvmceJtNCx6iA2NNWuNbtYYOugQhfwpHq1+01QrNn+7U6dzRboFqtlDKRt7LAnTP+DJGxo+Ys3rZoVnwEK2q1Re0amD0oCAFpg7tgWlhCt1PF3HI2mw7lFOUigJ0ZHJ56f0jmuEqlRUAAAjPSURBVON16oVre9QMDyx10Rchmo0oZh3jzfFqX4d8S7S9VdrSKVA7Warn5Fn/6zj63hag1t/6IkTbjTW6V6vVIgUOapW25FunZv2qU0jml7VaoK1D+fb+0NbfSmt0vj4GtbEmtf0NYt+v5uETWP+F6FiNicHeD07yAWvWb4pamO2H7a1urKIAnWjUYaxGbRzqaWCpPuZz9uZT34RorZH9bzpJdkjfuqXsssFxYPp1wiPjR7BiMPufds/4ssIT8ltfytZ5pmYLqdu9cfPngeYHkQoDlAP3nKg1oD42yyfbW1d9E6IAtfHsCxHDB/4nzZ/CUSYfqp0UnaRcdjWU2aEq+r6VaR95z1+VVBSgzfE69YamThEcmO1WaLPpboEu+iJE62PZoUhzLLVGx2o0Yao1OtYYYHhobqY78Ckj1i/2th1t5Y++igKUsVwrdLYv+bRKFjREtX+C2liT+lhtqmO81hDNRo2Ag4K0k9HG0EFXNrW+dFVO0jfrJ51ua5e/a31RgKpR8yWfC2xBQ7S2v0F9LDvsqTeCWkPUxlOQciBIoThM948PTvWbtuS/fO2Xj5r1k7LpPfJdWUUBClk3mC2cvjicHxhr0sidLjwwKpoN0RxKrdKhJs2xGo3hOo3RIWqDkweVHxsq/hgjo8WH63PVRWDWrsppey35uzK15A/dW63P/GlNg3sPnN408OJ+eHEPzV+5v38+9EWIttQbwEvZ8+YQgJgcFs3BGs2hgL2tL9dg2ibrA+r0H3EMde9gdxvV+kHRjA75w/R6el0bP3Cfida18oN7Yeil4Ijnxhl+Zh+1f93F5O4X5rjW1tJXITq4N2gdwEwOiebQgS9Kc6hTv4+YLJxTzjcJtcWnaNrjbODowPOp7RtBfYw5C9CIcIu2i74J0eXPHThEnxiuAcHkcOeQnOwYqDPTXFzT09giMtt3U2qNutfHOp8HOrRnwi3QBVQpRCWdB/wd2VHFLRFxQ9t6pfXnk00udUVEfL9qJQb3HNw3OQhMDtdoDpe3Jie6rDdbTMrO8ayldZ2uRhp4cf+iCdBesqSorKRjgTuAk4EngD+OiBfSuuuBq4BJ4P0RcW9afibwGeAIsnmb/lNEhKRhYDNwJrAbeFdEPFH2mbqGqKQ68HHgrWRzQn9P0paI+FFus3XAmvQ4G7gp/exq6F/3dFzeXDZELOubhrLZgtL+CWr7C5q4L+5ZLAF6yFnSpex1wNaIuEHSden1BySdSja18uuAVwBfkfSaiJhM+90AbCML0fOAL5EF7gsR8WpJl5DNO/+uss9VJaXOAnZExOPpF3E7sB7If/D1wOY0U982SSvb5m+ertmEFzsHKLhH06yqxRCgySFnCVkrs6jseuDNqfwm4GvAB9Ly2yNiDPi5pB3AWa1plCPi22lfm4ELyUJ0PfCXaV93AR+TpFSfjqqE6InAL3KvdzK9ldlpmxOBg0JU0gay9AcY+/K/fuLhCu/fL1YDuxa6EhUtprrC4qrvYqorwG/0UvileP7e+8ZvW11x82WS7s+93hgRG3Ove8mSsrIntBpsEfG0pONz+9rWYV/j6Xn78oPeP811PwIcR8m/eZUQ7TSK057KVbYh/UI3Aki6PyLWVnj/vrCY6ruY6gqLq76Lqa6Q1beX8hFx3mzVhd6ypFLGzMK+Zvw+VY6adwKvzL0+CXjqELYxs8NbL1lSVvaZdMhP+vlshX2dVLCvqTKSBoAVwPNlH6pKiH4PWCPpFElDZB21W9q22QJcpsw5wEhpf6iZHY56yZKysluAy9Pzy4G7c8svkTQs6RSywarvpv29JOmcdDbAZW1lWvu6CPhqWX8oVDicT/0C1wL3kp1acGtEPCLp6rT+ZrLRrfOBHWSnJVzZbb+kw/pFZDHVdzHVFRZXfRdTXaGP6ttLlhSVTbu+AbhT0lXAk8DFqcwjku4kG3yaAK5JI/MAf8KBU5y+lB4AnwI+mwahnicL61LqErJmZlbCZxKZmfXAIWpm1oM5D1FJ50n6iaQd6WqC9vWSdGNav13SGXNdpyIV6vqeVMftkr4l6bSFqGeuPqX1zW33BkmTki6az/q11aFrXSW9WdKDkh6R9PX5rmNbXbp9F1ZI+oKkH6b6VhkHmBOSbpX0rKSO513309/YkhQRc/Yg6wD+GfAqYAj4IXBq2zbnk3XqCjgH+M5c1qnHuv4usCo9X7dQda1a39x2XyXrsL+oX+sKrCQbAPi19Pr4fv7dAv8d+HB6/jKyQYihBarvm4AzgIcL1vfF39hSfcx1S3TqMq+IaACtS7Xypi7ziohtQOsyr/nWta4R8a1INzYguxLiJBZOld8twJ8Cn+PAuXMLoUpd3w18PiKeBIiIfq9vAEenU2SOIgvRBbnLd0R8g/JzGfvlb2xJmusQLbqEa6bbzIeZ1uMqDpwWsRC61lfSicA7gJvnsV6dVPndvgZYJelrkh6QdNm81W66KvX9GPCbZCdpP0R2F6B+nfC9X/7GlqS5vk3SrF0yOg8q10PSH5CF6BvntEblqtT3o8AHImIyazAtmCp1HSC7/di5ZOfufVvStoj46VxXroMq9X078CDwFuDfAvdJ+mZEFN9VZ+H0y9/YkjTXIbqYLhmtVA9JrwduAdZFxO55qlsnVeq7Frg9Behq4HxJExHxj/NTxSlVvwe7ImIUGJX0DeA0YCFCtEp9rwRuiKzTcYeknwOvBb47P1WckX75G1ua5rLDlSykHwdO4UAH/evatvn3HNzp/d2F6ByuWNdfI7uS4ncXujO7Sn3btv8MCzewVOV3+5vA1rTtcuBh4Lf6uL43AX+Znp8A/BJYvYDfh5MpHljqi7+xpfqY05ZozN0lowtV1z8nuy3WJ1LrbiIW6I4+FevbF6rUNSIelfRlYDvZTNm3RMSC3Cqx4u/2r4DPSHqILJw+EBELcos8SbeR3U9ztaSdwF+QZnPsp7+xpcqXfZqZ9cBXLJmZ9cAhambWA4eomVkPHKJmZj1wiJqZ9cAhambWA4eomVkP/j+qFKltCrvVCgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "B.assign(1000)\n",
    "solve_and_time()\n",
    "fig = plot(h)\n",
    "plt.colorbar(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Dealing with boundary layers\n",
    "\n",
    "#### Use a fine mesh\n",
    "* Computationally expensive\n",
    "* Offset by using multiple CPUs\n",
    "* FEniCS can be run in parallel using ```mpirun```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Using an adaptive mesh\n",
    "* FEniCS can do a goal-based adaptive mesh refinement\n",
    "* Implemented in the  ```AdaptiveNonlinearVariationalSolver``` class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Stokes flow in a channel with a cylinder\n",
    "<center>\n",
    "<img width=\"500\" src=\"flow.png\">\n",
    "</center>\n",
    "\n",
    "\n",
    "The incompressible Stokes equations are given by\n",
    "\\begin{align}\n",
    "\\nabla \\cdot \\mathbf{T} &= 0, \\\\\n",
    "\\mathbf{T} &= \\mu (\\nabla \\vec{u} + \\nabla \\vec{u}^T) - p \\mathbf{I}, \\\\\n",
    "\\nabla\\cdot \\vec{u} &= 0.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The aim is to compute the force exerted on the cylinder by the fluid\n",
    "\\begin{align}\n",
    "F = \\int_{Cylinder} \\vec{e}_z \\cdot (\\mathbf{T} \\vec{n})\\, ds\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Weak form\n",
    "\n",
    "The weak form for this problem is\n",
    "\\begin{align}\n",
    "\\int_{\\Omega} \\mathbf{T}:\\nabla \\vec{v}\\,dx = 0, \\\\\n",
    "\\int_{\\Omega} (\\nabla \\cdot \\vec{u}) q\\,dx = 0,\n",
    "\\end{align}\n",
    "where $\\vec{v}$ and $q$ are test functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Meshing\n",
    "\n",
    "* We will create a mesh for this problem using Gmsh\n",
    "* This mesh will be loaded into Python and used by FEnICS\n",
    "* The labels of boundaries will also be imported\n",
    "* Makes it easy to apply boundary conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Summary\n",
    "\n",
    "* FEniCS has many features that make it well suited for research problems\n",
    "* Relatively easy to implement complex models and perform complex calculations\n",
    "* Functionality can be extended through external packages (mshr, meshio, multiphenics, FEniCSx)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.8.10"
  },
  "rise": {
   "scroll": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}