fenics_intro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Introduction to FEniCS"
   ]
  },
  {
   "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",
    "Poisson's equation shows up in a huge range of applications, so let's look\n",
    "at how it can be solved in FEniCS\n",
    "\n",
    "Specifically, 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": [
    "# FEniCS Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from fenics import *\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": [
    "# Looking into the code\n",
    "\n",
    "The FEniCS packages are imported using the command"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fenics import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This provides access to libraries for...\n",
    "* Defining the weak form, differential operators, and mathematical functions (UFL)\n",
    "* Generating low-level C++ code from UFL code (FFC, UFC)\n",
    "* Mesh generation, assembly of discretised problem, numerical linear algebra (dolfin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Mesh generation\n",
    "\n",
    "The mesh for the unit square is generated by the code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 20\n",
    "mesh = UnitSquareMesh(N, N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$N$ represents the number of elements in each direction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The mesh can be visualised using the `plot` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f09340eaca0>,\n",
       " <matplotlib.lines.Line2D at 0x7f09340fb0a0>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD4CAYAAADMz1tMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOy9aXAUV772+aT2DQntEkIbQiubJDC7bcBsxsbGZvEGyO3udjh87e7p6I7eonvue2M+vHfinQ/vnZh750bPjRmDsc1uQGxmsYUwAoywwBZCSIiS0AKSoEBSlaRacz6UT3Iy85zMFEiWCuUT4QhKeixVHuWvMs//PPk/giiKMGXK1PhSwGi/AVOmTP38MsE3ZWocygTflKlxKBN8U6bGoUzwTZkahwoarV+ckJAgZmVljdavN2VqXOjy5cv3RFFMVH591MDPyspCdXX1aP16U6bGhQRBaGF93bzVN2VqHMoE35SpcSgTfFOmxqFM8E2ZGocywTdlahxKF3xBEP5fQRC6BEGo5XxfEATh/xQE4aYgCD8IglA6/G/TlClTwykjV/xPAKzW+P6LAHJ/+u99AP/3k78tU6ZMjaR01/FFUawUBCFLw/IqgO2i7/neC4IgTBQEIVUUxTtP+uZ27dqF5uZmzJ07V9NXWVkJAHjuueee2OfxeHDu3Dld371791BXV4eAgAAsXryY67tx4wY6OzuRnJyM/Px8ru/cuXPweDwoKipCQkKC7jEsXLgQQUH8P99wjolRnyiKOHv2rK6vv79fynBo+drb29HU1ITIyEjMnj2b67t69Sp6enqQmZmJzMxM3WMoLi5GdHS0ru/ZZ5+FIAi6vtEY44KCArzxxhuaP09LgpHn8X8C/7AoitMZ3zsM4F9FUfz2p9enAfxJFEVVOkcQhPfhuytARkbG7JYWZrZA0r/8y7/oH4EpU+NU//zP/6zrEQThsiiKc5RfH47kHusjkflpIoriPwD8AwDmzJmj+4lTVFSE7u5ufPjhh8zv//DDD/jyyy8BABkZGfjFL37B9H377bc4ffo0AN+n/auvvsp6bzh69Kh0JVqyZAmef/55lc/j8WDv3r2or68HAKxbtw6zZs1S+RwOBz777DO0trYCAMrKysCKKNtsNmzbtg337t0DAHz88ceIi4tT+e7du4dt27bBZrMBAP72t78hMDBQ5WttbcWOHTvgdDoB8E+O+vp67NmzB16vF/Hx8fjoo4+Yvurqahw5cgQAkJeXh7feekvlEUUR33zzjXS1nzdvHlavVs8OvV4vysvLceXKFQDAqlWrMH/+fJXP7XZj586daGpqAgBs2rQJhYWFKt/AwAA+/fRT3Lnju7n89a9/jUmTJql8PT092LZtGx48eAAA+P3vf4+oqCiVr7OzE9u3b0d/fz8A/thZLBZ8/vnncLvdCA0NxZ///Gemr7a2Fvv374coikhLS8OvfvUrpq+qqgonT54EAMycOROvvfaayiOKIr766itcvHgRgO9u5Ek0HFX9NgDp1OvJADqG4edqikCfnZ2NpKQkBAcHM30E+hkzZiAyMhIBAepDpqEnJyLrFo+Gft68edz3RqBva2vDM888w/UR6Ht6elBayq+JEuhFUcSMGTO4PgJ9VFQUcnNzmSc38Aj61NRUTJ48GWFhYUwfgT4vLw9xcXHMDxoa+tLSUgQHBzPHjoaeBTsRDb3WGBPou7q6MGeO6oImiUDf39+P4uJiro9AHxQUhKKiIuaxAo+gj42NxZQpUxAZGcn0EegzMjKQkpKCkJAQpo9AX1RUhAkTJnDPTwI9GROtKYgRDQf4hwBs/am6Px9Az3DM77VEQ//WW28Zgn7dunXMwaKhX7hwIZYvX878WTT0q1ev5p6UNPQbNmxAUVER00dD/8477yAjI4Ppo6EvKyvjzv9p6MvKyjBhwgSmj4Z+8+bNhqDfuHGjIehffvllXeh5d1KAHPpXXnkFJSUlTB8N/RtvvIHc3Fymj4Z+y5YtzLsBQA59WVkZYmNjmT4a+rKyMkPQv/3224agX79+vSHoV61axfxZQ5WR5bwvAJwHkC8IQpsgCL8UBOEDQRA++MlyFMAtADcB/D8A2Pflw6THhV7vSk+g17vSjwT0vIKUEvrERNVDVgDU0PMKV48LPauI6G/Qp6WlMX1K6FnTLGDsQP+kV3oiI1V99aRO/n0RwD8Ny7vR0c8BvbLYaUJvQv+0QQ/4UXLPvNKrZUKvlgm9MY3a8/hD0dMC/ZEjR3Shv3//Pg4dOqQLfVtbGz7//POfFXoAYx76vXv3GoJ+//794xZ6wA/A7+7uNgR9U1MTmpqaNKEHgJqaGgDQhB54FKbQgh4Ajh07BkEQNKEHgH379iE4OFgTegD4/PPPERkZqQk9AOzYsQPR0dGa0NtsNkPQt7e3o729XRf669evA4Am9KIo4rvvvgPAXxIlIkusWtADQHl5OQIDAzWhB4Ddu3cjNDRUE3pAPnY86D0ejyHorVarIegtFgssFosu9FevXgWAEYUeMBjgGQnNmTNH1OvAQwd4pk2bxoWerAsDvnVQvUEFgFmzZjEHtb+/Hw0NDdJr3hIQ/TtDQkK40NO+2NhYLvS0Ly0tjQs97ZsyZQoXetpXUFDAhZ72TZ8+nXt7b2TsPB4PfvzxR+k1b+w6OjrQ1dWl66PfW2RkJBd62peYmMiFnvZlZGRwoad9ubm5XOhpX1FRERf64Tw/vV4vfvjhBwBPFuDxC/AjIiK40Pf09Ej/1loH7e3tlV5HR0dzT1wSkAGAmJgY5u90OBwYHBzU9fX19cHr9er66GMYbl9AQAB3aY/28cbY6Nh5vV709fXpvje32w273a7rGxgYkEJIWr6RHDujvpCQEISHh+v6oqKiuMuijzPGo53cGzEVFhbi7t27+M1vfsP8Pp3Iy8rKQllZmcqjTOSVlJTglVdeUfmUibylS5cyM9PKRB4vudfX14ft27dLiTxeck+ZyOMl95SJPF5yj07kAcDf//53lQeQJ/ISEhLwT/+kXphRJvJ4+XBlIm/+/PnM9WZlIo+X3FMm8njJvZ6eHnzyySd4+PAhAH5yT5nI4yX36EQewAeLTuSFhYXhT3/6E9NHJ/LS09Px3nvvqTzKRN6sWbOwbt06lc/r9WLfvn2oq6sDoJ/719OYruoLgsCdb9KFvEmTJnE/SelCXlRUlG4hb+XKldz3QxfyeEEf4BH0PT09WLJkCddHV+8XLVrE9dHVe62UGl3ImzZtGje5RxfysrKyEBoaqvIoq/cJCQm6hbznn3+ee7tLQ681xnQh74UXXuD6CPQDAwOaENDVe63EIF3IKy4u5ib36EJefn4+IiIimD66kDd58mTuFIou5Gld6Qn0WmM3FI1p8Hkay4k8GvqRWLIzmshjwQyM3JLd888/z/2QG4nqPYH+aV+yU0K/YMEC5u8cqvwO/LG8ZDfS0I/VdXoT+kfyB+gBPwPfhF4tE3q1TOj15Tfgm9CrZUKvlgm9MY3pqj6REegBjHno9+3bZyiRd+TIEV3ob9y4gX379v2s0IuiOGrQnzhxwhD0p06d0oW+ubkZBw8eHLfQA34Afnd3t6ErPVki0oOehCT0EnlVVVWGEnmnT582lMg7cuSIoUTel19+ibi4OE3oAWDPnj1IS0vThJ7EhPWg7+joQEdHhy70ZKlTD3qydKqXyKusrDSUyPvqq68MJfIOHTqkm8gDfCnKxMRETeg9Ho8h6K1WqyHom5ub0dzcrAm9x+ORwjkjCT3gJwEeAMjMzNQcVKKsrCzmoDqdTrS3t0uvs7Ozmb/TarXKQhc8n8Vikf4dEBDAhZ72RUVFcaGnffHx8VzoaV9qaioXetqXnp7Ohd7I2Hk8Hty+fVt6zRuTwcFBae1dy0e/N6O+kJAQLvS0LyYmRvP2nigpKYkLPe1LS0vTjOESPen56XK50NbWJr3mjQl9Hj+1AR6iyZMnQxRFeDwe2ddFUZSCNIDvBKeTckQej0cGfUZGhupnAb7bexp6no9OT5H3x/KR8A5RXFwc00cfA+CLp7J8NHxhYWEIDg7W9SUnJ0MQBN2xmzx5MnPsvF6v7ITkjYnb7ZZBz/MNDAzIXvN8ygRdSkoK09fZ2Sl7HRMTozsmgG/89Hzx8fEIDAzU9aWmpv4s56fT6cTdu3dVX38cjWnwi4qK0NXVhV/+8peq75FCHhnYnJwcbN68WeUjc3qi0tJSrF27VuUjc3oiXnKPzOmJeMk9Usgj4iX3lIk8XnKvvr4ebW1t0onzhz/8gRkyqa6ulk7KsLAwfPDBByoPmdOTsUtNTWWOMSnkEfALCwuxadMmlY/M6Yl4yT0ypyfiJfdIIY+Il9zr7OyUjTEvuadM5PGSe7W1tWhtbZV6MvD6EFZVVUljHBMTg/fff1/lIXN6Msa8ZCmZUhDweT0hnU4nPv/8c6lnxFOd3APYvcWU1XveLaCykMdLstGFvPXr13PfC13IY/1xiOjq/Ysvvsj10dV7rZQaXb1fuHAh10cKebm5uZg1a5ZuIa+kpAQ5OTnMW1Rl9Z43RVEW8ni3xXQh7/XXX+ceA129Z31AExHog4KCNNtR0dV7rScF6ULe3Llzuck9upDH682nLOSlp6dzp1D79+9HXV0dVqxYwX0ugEB/+/ZtzbEbisY8+Er5UyJv69atSEpKYvpGKpGXm5uLTZs26fbIKykpwdq1a/12yY5A/+677w77kh3vgTAa+tdff90Q9FqFPBp63ge6Evrp01Ud7h9LfgW+P63TDwX64VqnJ9AbudKb0Pv0OEt2/g494Efgm9Cb0JvQD5/8AnwTehN6E/rh1Ziu6gPDD73T6Rx26E+dOmUI+m+//dZQIu/UqVOGoD9+/PiYhH5wcHDYoS8vLzcEfUVFhSHoT5w4MW6hB8Y4+KIo4t69e7h3754m9C6XS1oO0YK+v78ftbW1hhJ51dXVhhJ53377LSIjIzWhB4Cvv/7aUCLvxIkTuok8ALrQA4DdbteF3uPx4M6dO7hz544m9A6HA42NjQC0E3lOpxNXrlwxlMi7cOGCoUTemTNnMGHCBE3oAeDUqVO6iTzAlwTMzMzUhN7j8RiC/sGDB7rQu1wuKeOgBX1/fz/q6uogCMKIQg/4UXIvNjaWOahut1vWtoh3YthsNlkrJ57ParXKXg+nLygoiAs97YuIiOBCT/uio6O5V3qyTxzAHzuv1yt1r9E6BofDIWuXxfM9fPhQFlAZjTEODQ3VbIxJFBUVpRnDJYqJidGEnog3xh6PRxZG4h1Df3+/rJ0bz2e32+FwOACMg+Te9OnTuYNKWhEB4O4r53A4ZH9Mnq+vr8+QT5nI4/mUiTxWCAWAlIEnysnJYfroRpZZWVnMJUBRFFFbWyu9njZtGned/tq1a9Jr3jG43W6pu66Wb2BgwNDY9fT0GPIpE3k8nzL+m5eXx/TRYxIaGsqNxNJjXFBQwO1DSP+8Jz0/nU4nbty4oeuz2+2qD8PH1ZgGn+yWywrVKBN5OTk5zHADKeSRxBMvuUfm9CSiyUvuKRN5vORea2ur7I/JS+7V19fLAOQl96qrq2Un5ebNm1VXIjKnJ4qMjMSGDRtUP8vr9eLQoUPS67S0NObYKRN5hYWFTB8p5JGx4yX3lIk8XnJPmcjjJfcsFovsQ4mX3KutrZWB+tFHHzHDXFVVVdIYC4LA7C9I5vREsbGx3PNz//790uvs7Gzm2CkTebzknt1ul4qbbrf76U/usaQs5PHmh8rqPS+5pyzk8URX7998802uj67es7Y8JqKr91oJP/rRWt4fXFnIKykp4d7eHzp0CFevXsXzzz+PqVOnMn+espDHS+4pq/e822e6kMeKVhPR1XtWPJiIrt6zmqcS0dX7FStWcH10IW/BggWGHq3VutLThTzehqjKQh4vuUegf/DgAd5++23uMQxFfgf+07pkZ7RH3pMm8mjox9uSnZFE3pM+Tz/c1Xsl9LwpylDlV+A/rdAPVxMNE/pHGo0mGv4CPeBH4JvQm9Cb0A+f/AJ8E3rtDSxN6B/JhN6YDFX1BUFYDeDfAAQC+C9RFP9V8f0YADsAZPz0M/8PURT/v+F4g8OdyLPZbIagv3//vmHojx49agj6S5cuGYL+/PnzutB//fXXqKqq+lmhH4lEXmVlpSHoa2pqDEF/4cIFXejPnz+PqqqqcQs9YAB8QRACAfw7gBUA2gBcEgThkCiKdZTtnwDUiaK4VhCERAA3BEH4TBRFJ+NHGpbH44HVaoXVatWEvq+vD729vbqJPJvNhoaGBkOJvB9++MFQIu/y5cuGEnkXL140lMgzAj0AQ9D39/frQj84OIjOzk50dnZqQv/w4UN0d3frJvKcTieuX79uKJF35coVQ4m8S5cuGUrknT9/XjeRB0AXesA3fkYSeT09Pejp6dGE3mazwWq16ibyyM8KCgrShf5Jg3e6yT1BEBYA+G+iKK766fVffvrF/53y/AVAOnwfAFkATgLIE0VR3WfoJw01ucdrjKBsUfQ0+0RRlCXjAgICuNDTf9exdAxjxScIAvdKb2SMx8IxjHRyLw0AHUFrA6C89P5fAA4B6AAwAcAbLOgFQXgfwPsAuGubLPH2lXM6nbh06ZKuz2azybYg5vnu378vS9HxfG1tbWhpadH13bhxQ5by4/kuX74s/UFjY2O5dyznzp2T/j1nzhzunndVVVW6v9PtdksbNWr5BgcHcfnyZV1fb2+vLGDE83V1dUmZfy1fc3OzrA8dz3ft2jVZ7JjnO3/+vPTvrKws7p0IPca8n+X1emU/70nPT7vdLttKm+d78OCBLAn4JDICvvrjDlDeJqwCcAXAMgA5AE4KgnBWFMVe2f8kiv8A8A/Ad8XX+8UkucfqsqNM5OXk5DB9pJAXHBwMl8uF0tJSpo8U8sLDwzEwMMBN7rW2tuLSpUvS7+Ul9+rr62UnBy+5V11dLctob968WXXbq0zkAb6HkZRXBWUiLyoqinmsykReWloa06dM5BUWFjJ9ZE4fGhoKh8OB+fPnM32kkBcSEgKn08lN7lksFtmHEi+5V1tbK/uQ4yX3qqqqZFfx9evXq8JcykReQEAA8xiUiby4uDimT5nIy87OZvqUibzi4mKmz2q14pNPPkFERAT6+/t/luReG3y38UST4buy0/oFgP2iTzcBWAAUPNE705Cyes/79FZW73nJPWX1nidl9Z4nunrPin0S0dV7Xn855WYXvLm6spBXWlrK9CkLebzknrJ6z0vuKQt5vLk1gT44OFhz7Ojq/ZYtW7g+unqv1SeRrt6vXr2a6VEW8hYtWmSokKeVvacLeby7W2Uhj5fcI9B7PB7N83MoMgL+JQC5giBkC4IQAuBN+G7rad0G8AIACIKQDCAfwK1heYcKPa1LdnqJPL1trfxhyY5AP9aX7Fgareq9Evrk5GSmb6jSBV8URTeAjwB8BeA6gN2iKF4TBOEDQRBI7+b/DcBCQRB+BHAawJ9EUbzH/omPr6cV+ifdy86EXq2nYZ1+pKAHDK7ji6J4FMBRxdf+k/p3B4CVw/auGDKhN6E3oR8+jenHconGE/QAhrSBpQm9XEagBzCuoQf8AHyj0A93Iq+trQ3nzp0z1CPvxo0bhqC/du2aIehra2s1oQd8G0T++OOPfg/96dOnDUF//fp1Q4m82tpaQ9DX1taOW+iBMQ5+V1cXent70dvbizVr1iAjIwM2m03lo9d7N27ciPj4eKbPZrNJX9+6dSsiIiJUPrLs09jYiJCQEGzYsAEBAQEqH2mndP36dcTExOC1116D2+1W+UjHlGvXrmHSpElYs2aNbPmOiKz319bWIj8/H0uWLJG1uyIie6f9+OOPmDNnDubMmcM81ubmZgwODqKpqQnLli1Dbm4u03fz5k3p36+88gpSU1OZvu7ubunfb731FmJiYpg+p9MJi8WCgIAAvPnmmwgJCVH5yP55DQ0NiIiIwPr16yGKospH1ufr6uoQHx+PV199FU6nU9ZCDfB9mAO+scvMzMSqVavQ39/PPYba2lpMnz4dixcvZo4xOZ/q6uqwcOFCzJw5k3msDQ0NcDgcsNlsWLVqFbKyspg+ep+99evXIzExkekjyT3At6wbGRnJ9AGQtZt7HPlNzz1TpkzJ9dT33FuzZg3z6zabDZWVlbq++/fvywIhPF9bW5u0P7mW78aNG2hqatL1Xb58WeodFxERwb0l/+abb6SrYGFhIfM2UBRFHDt2THq9YsUKZmMJt9uNEydO6L63wcFBfP3117q+np4eWZqN5+vs7JQl/Hi+5uZmWfqM57t27ZosHcnznT9/Xmp6mZycjNmzZzN9X331lZSOfOaZZ5i5BK/Xi+PHj+v+TpfLhZMnT+r67HY7zpw5o+uzWq24cOGCrq+jo0OW8HsSjWnwSXLvmWeeUX3PZrNJ80iXy4WcnBym7969e6isrERkZCTsdjtKS0uZvtbWVpw6dQoTJ07Ew4cPsXTpUqavvr4eFotFSlDxknvV1dXo7OyUUmobN25UJfdI9X5gYEBKbi1fvlw1JyZ97wFIKbp58+ap1v7pRF5AQAAiIiKYx6BM5KWlpTF9PT092LZtm5TIKywsZPo6OztRUVGB6Oho9Pb2Yv78+UyfxWLBiRMnJN+qVauYvtraWty+fRthYWEYHBzkJveqqqrw4MEDaexeeeUVVXKPVO89Ho90vM8995wqzOX1erFv3z4Aj3LyrPdGJ/IAX8Sa5VMm8rKzs5k+q9WKs2fPSudTcXEx09fe3i7VQx48eDA+e+4R6Ekhj1dkoqv3ZWVl3OTe4yTytPrGKav3LCmX7F566SWmj97sYsmSJdw/uLKQV1xczPQpC3m85B6Bvr+/H1u2bOEm9zo7O7F9+3YEBQWhrKyMW3gj1fu4uDjN9Jmyes+TsnrPktFwDoG+rq4OK1euZMaIgccv5PGSe8pCHi+5197ejk8//RQRERGa5+dQ5HfgK6E3smRXVlbGPXHH+jo9DT1vm+eRqN7T0GtV72no9ZbsCPRjKZGnhH7BggVM32hV75XQ8z4chiq/At+EXi0T+scP54xX6AE/At+EXi0TehP6x9WYLu4RjSfoSSLPCPR79uwxoadkQm9cYx58o9DTiTwt6Nva2lBbW2sokXfmzBlD0Le1tekm8ioqKnD79m1Dibzbt29rQg8Au3btgsViGZPQNzc3o7q62hD0HR0dhhJ5ra2tuom8EydO4Pbt25rQA8D+/fvR0tKiCb3H4xk16K9fvz6i0ANjHHyy3jswMICVK1ciODgYHR3KVgDy5N66devgcrmYPjq5t2LFCtlrIrLW29HRAUEQsGzZMuZ+ZSRB19bWhqioKCxatAhdXV0qH+noc/v2baSkpKC0tFTaOZUW6RB0+/Zt5OXlITc3l3kMpPGExWLBnDlzkJyczPR9//33AIA7d+5gyZIliIyMZPro5N6aNWsgCALTRyf31q9fj8HBQabP6XRKY7N8+XJZGo2IjHl7ezuCg4OxZMkS1X6E5BgB3x3axIkTsWDBAuln0yLnye3bt5GRkYEZM2Ywx5jkDFpaWjB9+nRkZmYyj4HkFlpaWrBgwQLExcUxfaTjEGkWExoayvTRyb1XXnkFHo+H6aPH6dVXX4XdbmcmC8l7exKZyT1TpvxUT31y76233mJ+/d69e7IEFc/X2tqKb7/9VtdXX1+PmpoaXV91dbXUNy4sLIy5Px6Z05PkXm5uLubMUY2/VMgj2fJnn30WkydPVvnInJ6IN60YGBjAgQMHdI+hp6cHR48e1fV1dnbKEn48n8VikaXPeL7a2lpZbz6e7/z582hubgYAJCQkMPe9I3N6ktybNWsW8yEuMqd3u90AgJUrVyI+Pl7lczqdUohH673Z7XZZizOez2q1ytp58Xzt7e2yBCrP19jYCHKxfNIpwJgGnyT3WFsf37t3D+Xl5YiMjITL5UJ6ejrT19raiu+++w5xcXGwWq0oLS1l+urr63H16lWkpqbizp07WLp0KdNHoE9MTER3dzdWr16t8tHQk9+7cOFCVXKPhp74iouLVXNnOpFHklv5+fmq5B6dyIuMjITX62UeQ09PD44dO4bQ0FAEBAQgLi6O6evs7MTBgwcRHR0Nu92OvLw8ps9isaC6uloak/nz5zN9ZNfa5ORkdHZ2YtWqVUxfVVUVmpubER8fj/v372PZsmXMMSbQk7GbO3euKrlHQ098M2bMUIW56EReTEwM+vr6mO+NTuSFhYUhJCSE6bNarTh8+DAiIiLg8XgwadIkpq+9vR0XL16U/q7FxcVMX2NjI2pqaqTzk5UWHYr8ZjmP1uNW73nJvZFM5L388stM32gm8ow+WjuU6j15tJZXoBvJRN7KleweMI+byJs2bRrTN1rV+8bGRuzatQtJSUmafQiHIr8D/2lYshvNdfqRhN5csht56MPDw5m+ocqvwDehV8uE3oT+ceQ34JvQq2VCb0L/uBrTxT0iE3q1BgcH8dlnn5nQUzIKvcvlwt69e8ct9IAfgG8U+ra2NkPQ37hxA1euXDEEPans6vXI6+3t1U3klZeX4+HDh4YSeVarVRN6ANi+fTvu37+vCb3NZhtW6Jubm9HY2GgI+osXLxpK5PX29hrqkdfb26ubyNu3bx8ePnyoCT0AfPbZZ7BarZrQezyeUYP+xx9/HFHogTEOPklk2Ww2LF26FPfu3WMmvOjk3ty5c9He3i77GhGd0ps9e7aUDKNF1nr7+voAANOnT5ft9UbU1tYGwNf7LDw8HDk5ObJ994hI55+HDx8iOTkZSUlJuH79uspH1oWtVityc3MRFhbG9JGttLq6ulBSUgK32830keTew4cPsWjRIql3oVIkuUeagJCdc5Wik3tz586VpdFo0T3xSkpKZJ2KiMj7IDscFxUV4caNGypfQ0OD5IuJiUFGRgZzjEnO4OHDh0hPT8fEiROZY0JyC/fv30dRURECAwOZPpLcI01gBgcHmT6SR+jr68Pzzz8v7eysFD1WS5YsQUdHh25yb86cOVKOgaXLly9j6dKl3O/ryUzumTLlp3rqk3sffPAB8+utra04cuSIrq++vh4VFRW6vurqaikZNWHCBLzzzjsqD5nTk6tRfn4+85OXbHZBsuWLFy9m3l6UZtYAACAASURBVFaSOT25G1m3bh1SUlJUPlLII12Af/WrXzGnHz09Pfjiiy90j7WzsxNffvmlrs9iscjSZzxfbW2tLB3J850/f156LiElJQXr1q1TecicnlzxiouLmWvwpAU26bK7fPlyZn6BFPJId+M333wTEydOVPnsdjs+/fRT3WOwWq3YvXu3rq+9vV1qmabla2hokKUjeb6amhrpDjInJ4fpMaoxDT5J7rHmTa2trTh58iTi4uLgdruRmJjI9NXX16OyshJpaWno6urCjBkzmD4CfU5ODpqamqQHYGjR0GdnZ8NisaCwsFDlo6HPyspCc3MzcnJyVD4aeuJLT09XzbEHBgZw8OBBCIKA9PR0tLa2IiUlRZXc6+npwc6dOxEaGorY2FjYbDbmsXZ2duL48eOYMGECgoODER4ezvSRvveJiYmw2+3IzMxk+mpra3Hu3DlkZmZKT8exfFVVVbh69SqmTJmCW7duYdasWcwxJtCTMcnLy1P5aOiJLzs7W+Uj0DscDmRmZqKlpQVpaWmqMJfdbse+ffsQFBSE5ORk3L17l3kMVqsVx44dQ3h4OCZMmAC32830tbe348SJE5g4cSJEUURcXBzT19DQID0FSqYgLB+BnpyfvHqNUfnNch4tZfV+woQJTJ+yes/aTx54VL3Pzc3Fpk2bmB66el9SUqKZyKO3tTJavR/uRB5ru2jg0WYXQUFBePfdd4e9es/q/AsMb4885WYXRpfsjCbyWNuZA48KeW63G2VlZbqFvPDwcLz77rvcQl5DQwN2796tW8irqanBoUOHkJOTgzfffJPpGar8DvyRWrIj0Ost2ZWUlGDt2rV+vZfdSEE/HvayG23oeatLQ5VfgW9Cr5YJvVom9PryG/BN6NUyoVfLhN6YDIEvCMJqQRBuCIJwUxCEP3M8SwRBuCIIwjVBEM6wPI+r8QT94OCgCb1CJvTDX4PX/YmCIAQC+HcAKwC0AbgkCMIhURTrKM9EAP8BYLUoircFQWBvQfsY6u7uNgR9U1MTLBaLLvQk2KIFPQBp+U8LeuBR8EYLegDYuXMnPB6PoUSe2+3WhB4Atm3bBofDoQk96VeoB317e7sh6K9fv476+npd6ElDDiM98gAY6pEniqIm9ADw+eefw+12a0IP+MbY6XRqQu/xeAxBb7VaDUFvsVhw+/ZtXeivXLmCK1eujCj0gLHlvLkAboqieAsABEHYCeBVAHWU520A+0VRvA0Aoiiqm889hkhyz+l0Ij8/X7bnGi2S0vN6vZg6dSp3fzE6uZeRkQFWgIgk94gSExNl++4R3bp1S/p3WFgYwsLCZB1oiE6dOgXAdyIlJSXB4XAwfWSt3OFwYMqUKbh//760Pk2LJPd6e3sxffp0tLa2orW1VeUjH3ADAwN45plnpNyBUnRyLz8/X9YdhxZJ7omiiClTpkg/Xyl6J9u0tDR89913Kg/ZBZcoJiaGOcZkvV8URURHRyMgIEBz7FwuF9LS0mCz2TR9/f39yM/P56YUSXKvr68PxcXFsFgszJQnGSuHw4GSkhJmug94lNzzeDzIzc2VjkspOrmXlZXFPD+JKisrRza5JwjCBviu5L/66fUWAPNEUfyI8vxPAMEApgGYAODfRFHczvhZ7wN4HwAyMjJm6zUMNJN7pkzxNdLJPdb9l/LTIgjAbAAvAAgHcF4QhAuiKMouM6Io/gPAPwBfZFfvF5M2Q3/605+Y36+vr8fBgwd9b1IQ8Mc//pHpq66uxunTpwEAkyZNYnYxIXP6S5cuAfClxVjryOQpO5IZX7ZsGXOTQzKnJ3cj69evZ67Bk0Ie6Rv3y1/+EgkJCSofmdOT9Nkf/vAHVYAHeDSnJx/ovLGzWCyy9BnPV1tbK6UjIyMj8dFHHzF9VVVVOHv2LABfqmzDhg0qD5nTkyvevHnzmNMjMqcnd1WrV69mtpoic3py1X7rrbeY+9SROT15/uLDDz9kZj/InN7lcgHgj0l7ezt27Nghveb5GhoaZOlInq+mpkaa9iQlJeEXv/gF01dRUSHdGbH6Nw5FRsBvA5BOvZ4MQPmEQRuAe6Io2gHYBUGoBDALAPv+0qBiY2PhdruZ8/X6+nqUl5cjLS0NDocDMTExTB+BPi8vT0q8KX009MXFxbhy5QpiY2NVPhr6WbNm4erVq4iOjlb53G439u7di/b2dsycORM//PADoqKiVL6BgQHs2bNHum2vra1FRESEykcSeaIoIj8/Hzdu3EBYWJgK/M7OTuzatQtRUVFISEhAd3c3c0wsFgv279+PxMREBAUFISAggOmrra3F0aNHkZmZiZ6eHqSmpjJ9BPqioiI0NDQgMTGROcYEejLGEydOVPk8Hg/Ky8ulZB9vjJ1OJ/bv3y+lMX/88UfmGNvtduzatQsDAwMoKipCXV0dwsPDVT6r1YqdO3ciODgYWVlZuHXrFvNY29vbsWfPHkycOBHR0dGw2WxMX0NDAw4dOoTU1FS43W7mewMeQZ+Tk4O7d+9i0qRJTB+BnoxdRESEyjMUGanqXwKQKwhCtiAIIQDeBHBI4TkI4FlBEIIEQYgAMA8Ae8IzDBpqIo88Wsu6Qg4lkUeepx9KIs/okt1wJ/JiY2OZvpFO5Bmp3vN65D1uIs9o9V4vkUcerU1KYtemR6t6X1FRgTNnzqC4uBhr165leoYqXfBFUXQD+AjAV/DBvFsUxWuCIHwgCMIHP3muAzgO4AcA3wH4L1EUa4flHSo0Uk009JbsaOjNdXqfnqYlO6PP04829LwVkqHK0FqBKIpHARxVfO0/Fa//B4D/MSzviqOfA3plsdOE3oT+aYMe8KPknnmlV8uEXi0TemMa04/lEo036A8dOmQI+h07dpjQU7Lb7Thw4IAh6Pfu3TtuoQf8APzu7m5D0Dc1NaGpqUm3Rx7ZIksvkXfmjC91rJfIO3z4MADoJvJ2796NwMBA3UTe9u3bERoaqgk98YWEhGhCb7PZDEHf3t5uCPrr16/j+vXruok8suSkl8gj4Sa9RN6BAwcgCIJuIu/zzz9HUFCQJvSAL/UYFhamCb3H4zEEPWnKoQc9CQHpQU/CZyMJPTDGwSdJPa/Xi8TERAlGpej+etHR0dKaPS1RFGXJvZCQEGntlJbD4VC9pjvQENXWPqpdhoaGoquri+mjE2QTJ07ErVu3ZKk/li85OVnabornI40lSO5AKZKsc7vdSElJkXXHoUWSe6IoIj4+XkoGKkX33IuKipLtWUgkiqIsuQeAOcb0zwJ8CTm9sYuMjOT2UqR9CQkJaGhoYCYVic/r9SI5OZmb8KTHeMqUKcwUIPAouefxeJCcnCzb/44W3XMvNjaWeX4C8uReSEgId4yBnyG5N1Iaas893lWIPtECAwOZS3ZKX3BwMPMqJIqiFN7Q+p1ut1tqgWX0vT0tPrLur+fj/Syv1yuLRY/lYx1uX0BAAPdK/zjn51Pbc6+goAAdHR343e9+x/w+mdMDQHp6Ot577z2Vh57TA75bqFdffVXlo5+yA4ClS5cy97Oj5/SAr0ceK1VG5vRkn/aysjLmOjIp5JH8+scff8y8dSeFvIGBAQDA3/72N+aHHJnTE7h4JweZ04uiiIkTJ+K3v/0t00fm9IDvwSbWvnf0nB7w3eKvXr1a5aPn9ACwatUqZi89Mqcnke5NmzahsLBQ5SOFPHIH8etf/5qZcyCFPJLc+/3vf8/cR5HM6cldH2/syJze4/EgICAAf//735m+77//Xuq5l5qaivfff5/pI3N6AJgxYwZef/11lUcURZSXl0tT1WeffZb5s4xqTFf1AwICDIVzUlJSmJ+4ys0uIiMjmVcrGnqtAaWh1/LRhTytuStdvZ87dy7XR6APDg7mtugC5IW8goIC7iahdCEvIyODO/enC3nx8fHcYildyAsJCdEt5PE2CAXkhbzFixdzfXT1nrcZJiCv3s+ePZvro/vez5gxg3vnSBfypk6dymzaCTyCPicnh5vGA+SFvAkTJnBDZgR6ct7x6iZGNabB52moiTy9HW7GQhMNvURecHAwysrKDCfyeJHOka7es0RD7w/bWvEe/TZavaehf9LqvRL6J5nX0/I78Id7W6uxAL3eOj2B3t+X7PwB+iddshtp6J/0Sk/kV+Cb0KtlQq+WCb2+/AZ8E3q1TOjVMqE3pjFd1ScyAj2AMQ/9wYMHDUFfXl5uCPpdu3b9rNADGPPQnzhxwhD0X3755biFHvAD8Lu7uw1BT5J7WtCLoiiFNvQSeWT5Ty+R99VXXxlK5B04cMBQIm/37t2Ijo7WhB4AvvjiC8THx2tCb7PZDEFPgjF60JPWUnpbVZNQkd6utd98842hRN6RI0cMJfL27duHiIgITegBX/9D8oHJg97j8RiC3mq1GoKeJPf0oCfLySMJPeBHAZ68vDzuEgvd66ygoIB7QtI7srLWhQEfLHQPO56P/p1BQUFc6GlfTEwMt3pP+5KTk7nQ0z6tpTjaN3XqVO7z9LQvPz9fF3qAPyZut1u2szDPd/v2bdjtdl0f/TvDw8O5z9PTPt5WVUrfpEmTuNDTvuzsbO5SHO3Lzc3lQj+c56fH45FSiU9tgIcoLi5Oak2lFB3/TEhIYDaoFEVRtr12YmIic7ttt9st+z08H2l/RRQbG8v00RFMwJfuYvmUEVav12vINzAwIAV6eL7IyEjV+2D54uLimFs8K328MfF6vbKx5/lcLpcMep6P9gC+mDDLp3zPgYGBhsbO5XLp+oKDg2Gz2WRRb5YvOjpa1UCU5XvS89Pj8XD/RkPVmAa/sLAQnZ2d+Pjjj1XfI4U8MrBTpkxh9tIjT9mRgSwtLWV2MSFzegI+L7lH5vQkBcZL7pFCHoli8pJ7nZ2d2L59O/r7+wHwk3vKRB4vuUcn8gIDA/GHP/xB5QHkibzk5GTmDq2kkEfGuLCwkLm3oNfrxb59+6STev78+cw1fTKnJ5DwkntkTk/GhJfcUybyeMk9ZSKPl9xrbGzErl274PF4AAB//etfVR7g0VN2gO+DlZcspRN5GRkZzF56ZE5Pzs9Zs2YxdxB2u93Ys2ePBL5WCMqIxnRVXxAEQ+GcSZMm6T5au2TJEm6SjS7k8VpvAfJC3osvvsj10dX75cuXc30E+qCgIM0/JIE+Li5OM+FHF/JmzJjBnZfShbwpU6YYSuQlJiYyfxaBnhTytDLrpJD30ksvcY+BLuTxAkGAvHq/bNkyro+u3i9atIjrI9AnJSVh9uzZ3Gkl/WhtYWEhN1lKF/ImT55sKJEXHR3NPI8J9A0NDZpjNxSNafBZGsqSHQ39cPfIG44lOwK9kSW7uLg4bN26dVR65PEAVELvr0t2NPRbtmzhzumH83l6o9V7JfRP2l2XyK/AN6EfO+v0Tyv0P0cTjdGGHvAj8E3oTehN6IdPfgG+Cb0JvQn98GpMV/UB49CTQR3L0FdWVupC39zcjGPHjhmC/uDBg+MG+mPHjhmC/vTp04agP3LkyLiFHhjj4JP1TSNXetLOSgt6l8uFH374AYB+Iu/ChQuGEnnffPONoUTe8ePHDSXyysvLkZSUpAk94EsCZmZmakJvs9kMQd/R0YGOjg7dGC4JmOhB//333xtK5H377beGEnmnTp0ylMg7evSobiIPgLTDjRb0Ho/HEPRWq1UXeq/Xi9u3b6O5uVkTepfLJbVbG0noAT9K7mkt2ZEuNwC48A0ODsrCEzxfV1eXrL0Rz6fs/WbEFx4errlVNdHEiRM1G2MSJSYmasZwiVJTUzWhJ+KNsdvtlu0qyztWu90uC7IM59gFBQVpXumJoqKiuNDTvri4OC70tC85OZkLvZExVp6fvDF2OByy0I6R8/ipT+7l5OToXukBcLehcrvdsj8SzzcwMCCDXmtbK1o8n3ILZt4fU7kFM2vTTOBRY0zAd+LyTnDal5WVxV2nJ+3DAO0xNjJ2TqdTBj3Pp0zk8XzKhBovrqv8EElJSWH66DEBwP0Apn2TJ0/mLu3RvuzsbO46vZExNnp+Ki9eT6IxDX5RURG6urqwefNm1fdIIY8oJycH77zzjspH5vREvOQemdOT3XR4yT0ypyfiJffIo7VEvOSexWJBc3Oz9JqX3KutrUVTU5PUZfXDDz9knmxVVVXSSRkeHo6ysjKVh8zpyUk5adIk7hjv27dPes1L7pE5PRk7XnKPzOmJeMk9Usgj4iX3SCGPiJfca2xshMVikRJ5vOReTU2NDOhf/vKXKg/gm9MT38SJE7F161aVh8zpibKysphj7Ha7ZbsW83pCDg4OYseOHQgICIDX6326k3sAu7eYsnrPu5IqC3m85J6ykMcTXchbv34910c/T6+1ySFdvV+xYgXXR1fvtfrQ0YW84uJi7lWILuRpXenpQh4vuacs5PGmHnQhj/XhQURX71977TWuj67er1mzhuujq/daCT+6kDdv3jxuco8u5E2bNs1QIS89PZ17pd+9ezcaGxvx0ksvce/gCPR37tzRHLuhaMyDr5Q/LNkNpYnGaCfyxlr1fjwl8mjoeYU8JfT5+flM31DlV+CPd+if9iW78bROP5rQA34Evgm9Cb0J/fDJL8A3oTehN6EfXhmq6guCsBrAvwEIBPBfoij+K8f3DIALAN4QRXHvcLzB4YZ+cHBw2KH/6quvDEFfVVVlCPqKigpd6M+fP4/Tp0//rNC7XK5Rg/7gwYOGoD979qwh6L/++utxCz1gAHxBEAIB/DuAFQDaAFwSBOGQKIp1DN//DkC9++FjinR1uX//vib09PqmFvQ2mw11dXWGEnnff/+9oUTe+fPnDSXyKisrDSXyvvnmG91EHgBd6AEfgHrQu1wu3LlzB3fu3NGEvr+/H93d3bqJPKfTiR9//NFQIu/SpUuGEnnnzp0zlMirqKjQTeQBMAS9x+MxlMh78OABHjx4oAm9w+GQMh1a0Pf19aGnpwcBAQEjCj1gILknCMICAP9NFMVVP73+CwCIovjfFb7/BYALwDMADutd8Yea3OMtxTmdTtmGgzyfsn3SaPgEQdBsjEkUEhKiGcMlioiI4F6F6KBMZGQkN4ZLt+7SWu4ka+BavrEwxoGBgVzoaV9YWBgX+uEcY7fbLWvVNpzHOtLJvTQArdTrNgDzFD88DcBrAJbBBz7vTbwP4H3A14rIqEpLS5lf93g8UldSLd/g4KC0UaOWr7e3Vxbg4Pm6urpkfwCejw7mAODeiVy7dk36tyAI3Ksp2foa8IVptNo96703r9cr2yaa53O5XNJ20Fq+/v5+1NfX6/oePHhgaOw6OjoM+ejmngCYYSoAsmONi4vjJgHpsSsuLuZCTzav1HpvRs9Ph8MhOwd4vr6+PtXxPq6MgM/q76u8TfifAP4kiqJHqx2wKIr/APAPwHfF1/vFRUVF6O7u1uyRR5STk6OZyAsMDITH4+Em98icPjg4GC6Xi5vcUybyeMk9i8UiPRAE8JN7tbW1spPoo48+Yk4ZqqqqZK/Xr1+vCpmQOT1RZGQk81hJ33uitLQ0pk+ZyCssLGT6yJw+KCgIbrebm9wjc3qSPuMl99rb22Ug8JJ7DQ0NMrB4yb2amhoZ0L/4xS+YV9SKigrp34IgMBN0ykRebGws9/ykE3nZ2dlMnzKRR6YWSvX29mLbtm0ICQmB0+n8WZJ7bQDSqdeTAXQoPHMA7BQEoRnABgD/IQiCumPgMElZyOPNwZXVe97tk7KQxxOBPigoiBm/JKJ3uNmwYQPXR1fvtXqp0dV7Xp1DWcgrLS3V3exixYoVmtl7upDHS+4pC3m8KYqykMcTqd6Hh4fjrbfe4vrozS5YzSmJ6Or9ypUruT66kLdw4UJDhbzp06cbKuTx7m7pQt7GjRu59QsCvc1m0zzvhiIj4F8CkCsIQrYgCCEA3gRwiDaIopgtimKWKIpZAPYC+FAUxQPD8g4VGu0lu6CgILz77rvDvq3Vz9EjTwk9bwvv0V6yCw8Px7vvvvvEO9w87pKdVr+Hkajeb9y4EQUFBUyfEvr09HSmb6jSBV8URTeAj+Cr1l8HsFsUxWuCIHwgCIK6J/MI6mmF/ufey86E/pFGe8luNKAHDK7ji6J4FMBRxdf+k+N998nfllom9Cb0JvTDJ79I7pnQD88Glib0jzSeoQfG+PP4wOgm8i5cuGAI+mvXrhmC/vvvvzeUyKuurtaF/vjx46iurvZ76CsqKgxB/+OPPxqCvrq6Whf6M2fOoLq6etxCD4xx8AcHB6Vk1EsvvYSZM2fKgiREnZ2d0tZSb7zxBqZMmcL0kX3QgoODsXnzZqSkpKh8Xq8XAFBXV4fIyEhs2bIFMTExKh/5fVevXkV8fDw2b96MsLAwru/y5cuYPHky3nzzTWlpkeWrrq5GQUEB1q1bB1EUVT7SIai6uhpz5szBihUrpPdMi+yrV1dXhxdeeAHz5s1jjonVapW63axbtw6FhYVMH9lGKyAgAG+//TYyMjKYPqfTiYaGBoSGhmLr1q1ISEhQ+cjr2tpaREdHY8uWLYiKiuKOyZUrV5CUlIR33nkHISEhmmOXnZ2NjRs3QhAEzbGbMWMGXnrpJeYYOxwOeL1e1NTUYOHChXj++eeZY9zX14fe3l709vZi9erVKCkpYY5JV1eXFOLZuHEjcnNzmb6enh709PRIq0aTJk1i+uhjflz5Tc89U6ZMyfXU99xbunQp8+uDg4M4f/68rq+np0cW4OD5urq6ZMERnq+5uVnWJ4/nu3btGrq6unR958+fl64IaWlpyMvLY/q++eYb6d+LFy9mLgF6vV5po0at3+lyufDtt9/q+ux2O7777jtd34MHD2TpOJ6vo6NDth00z9fY2Ii2tjZd35UrV6SNTkNDQ7nTmTNnzkhX7cLCQmZvPlEUZSEe3u/0eDyorKzU9Rk9P/v6+kBfBHm+7u5uqQvvk2pMg0+Se1q71pLb5pycHM0eeaGhoXA4HCgtLeUm8i5cuIDIyEjY7XZucs9iseDs2bNSSo2X3KutrZVtkcxL7lVVVcmy3K+//rqqnqBM5AG+NuLK5J4ykRcVFcU8BmUiLy0tjelTJvIKCwuZPjKnj4iIQH9/P+bPn8/0tbe3o6qqCuHh4RgYGOAm9xoaGlBZWSml2XjJvZqaGtm25lu3bmUm9yoqKmS36mvWrFGFuZSJvICAAOYxKBN5cXFxTJ8ykZednc30KRN5xcXFTF9XVxe+++47REVFwWazPf0991hSVu95BTplIY+X3FNW73lSVu95oqv3Wukzunr/yiuvMD3KJTutqxBdyOPlvZWFPK1uuHQhj5fcUxbytNp908/T86Ss3vNEV+83btzI9dHVe94Ox8pC3qJFiwwV8mbMmMH8ecpCHi+5pyzk8YqbXV1d2L59OwRB0By7ocjvwB9PS3ZPWyJP73n6sZLIY2m0qvdK6Hmt14cqvwJ/vEPvr0t2YwF6f1yyGynoAT8C34TehN6Efvg0pot7ROMJesBYIs/r9eLAgQMm9AqZ0BvTmAd/tKBvbm7G2bNnDUHf1NRkKJHX0NBgCPqGhgZN6AFg//79qK+vH7PQ19TUGOqR19LSYgj6xsZGQ4m8hoYGTegB4PDhw2hoaNCE3uPxPLXQA2Mc/ObmZvT398Nms2H16tWIj49X7akGyPdPW79+PcLDw5k+ktwDIG23pfSRpJTFYkFQUBDWrl0Lh8MBh8Mh85H1+Zs3b2LChAl48cUXVe2T6PfW0NCA1NRULF26VLbHHFFLS4vky8vLw9y5c2VLVUQNDQ0AgPr6esyZMwcFBQXMY62pqYEoimhpacHSpUsxadIkpo/uOLR27VrExMQwffTS5KZNmxAcHMz0OZ1OtLb6Gja98cYb8Hg8Kl9/fz8A4NatWwgNDcXLL78sJQ1p3b17F4DvAyIuLg4rV65Eb2+v6neS9f6GhgZptyHWGJPsRUNDA6ZPn47i4mLmGJNuTY2NjViwYAGmTJnCPFbSmai/vx8rV65EUlIS03f79m3p36+99hoiIyOZPpLcA4C3334bAQEBTB8A2QabjyMzuWfKlJ/qqU/u8Tqs9PT0yNJsPF9nZ6csQcXzNTc3y9JnPN+1a9dkvc94vvPnz0vdVRMTE7Fo0SKm76uvvpKudiUlJcjMzFR5vF4vDh161P/k5ZdfZt7yulwuHDlyRPe92e12nDx5UtdntVplKTWer6OjQ5bw4/lu3rwpS5/xfFeuXJF6FgYEBHBzDmfOnJGu2tnZ2cwwlSiKOHLkiJRvX7x4MfMW2uPxyEI8vPfmcDhw7NgxXV9fXx9Onz6t6+vu7sa5c+d0fbdv35YSqLy9/YxqTINPknusP2ZPTw+2bduG0NBQOJ1OTJkyhbtr7YkTJxAdHY3e3l6UlpZye+TV1tYiPj4e9+/fx9KlS7mJvJs3b0o/j5fcq6qqQmdnp5RmW7NmjSq5R6r3AwMDUppt8eLFqroDvWst8bE2xaQTeSEhIQgODma+NzqRJ4oiUlJSmD6r1Yqvv/5aOobCwkKmr729HceOHUNsbCwePHiA+fPnM32NjY24fv265Fu1ahXTV1NTg+bmZimltmHDBmZyr6KiAg8ePJDe3/Lly1XJPVEUcfjwYbjdbmns5s2bpwpzud1u7NmzB4CvA6/L5WK+NzqRFxgYiAkTJjB9dCLP5XIhKyuL6evq6sLJkyelYy0uLmb6WlpacOTIESQkJODevXvci4hR+c1yHi0CfX9/P7Zs2cKMaQI+6Ldv346goCCUlZVxk3uPk8h7++23uT5l9Z4lo+Ec5WYXrIgroC7kTZs2jekbSiFv27ZtcLvd2Lp1Kze5p6ze84qbyh1ueFJW73lSVu9ZItB///33uuGcPXv2oKGhAS+99BJmz57N9I10IY9XBG1pacFnn32GmJiY8ZvcU0KvVb2noR9LS3aPC/3Pua0VDb2/bmv1ONCPteq9EnrexWuo8ivwn1box9pedk8j9Lx1+vEIPeBH4JvQq2VCb0L/uBrTxT0iR56/7QAAHFhJREFUE3q1XC4Xdu7caUJPySj0Ho9nXEMP+AH4ww19c3MzfvjhB0PQ37t3z1Air7OzUzeRd+LECdy9e9dQIk9vA0sA+Oyzz9De3v6zQt/R0WEY+kuXLhmCvqury1Air7Oz01Ai7+7du5rQA8Du3bvR0dGhCb3H4xk16K9duzai0ANjHHySoHI6nVi+fDlcLpdqPzpAntxbu3at1AdNKTpZ9/zzz6O7u1uWSAMeJffI1xcuXIiODuXGQZC+RpbsZs+eLUtoEZFcwN27d5GSkoL8/HwppUeLrOPeuXMHeXl5SE1NZR4ryS20tbVh9uzZiIqKYvrIei9pZCIIAtNHknsDAwNYvXo1BgYGmD56nF588UWpF6JSTqdT2rn42WeflXIMtMjfoaurC4GBgZg/f76s2w5RU1MTAN8YT5w4ETNmzGCOMQmC3b17FxkZGcjOzmaOMemu09HRgenTpyMhIYF5rORv0d7ejnnz5iEsLIzpI8k9q9WKZcuWwePxMH30e37ppZdkCVJaJLUH+M7Pe/fucRN6N27c4BYsjchM7pky5ad66pN7vLXLzs5OHD9+XNdnsVhk6TOer7a2FpcvX9b1kQduAF8ib82aNSoPmdOTK96sWbNQXFys8pE5PcmvL1++nDmdcTqd+OKLL6TX77zzDvPW2G63Y+/eRzuU847BarXKUmo8X3t7O06dOqXra2xslG3syfPV1NRIm4kGBgZy94I7c+aMdPXMzMzEkiVLVB4ypyd59nnz5jFvyUkhj2ynvnbtWuZUcHBwELt27dI9ht7eXnz55Ze6vq6uLlnCj+draWmR9frj+erq6nDp0iUA4E7HjGpMg0+Se6xedZ2dnaisrER0dDRcLhcmTZrE9FksFlRVVSExMRH37t1DSUkJd9fa77//Hunp6WhtbcXSpUu5PfLIAzd37tzBokWLuIm8zs5OpKSk4O7duyguLlb5aOiJr7CwUHVS0om8xMREdHV1ITs7W5Xcs9vtOHr0KIKCghAdHQ2n08k8BqvViv379yMiIgKBgYGIjo5m+trb26UnFO12O3Jycpi+xsZGXLx4URqT+fPnM30E+rS0NLS3t2P58uVMX0VFBZqbm6WfN2/ePOYYl5eXw2q1SmM3c+ZMVZiLhj45ORmdnZ3Iy8tTzZ3pRF5cXBwePHjAfG+9vb0oLy9HSEgIIiIiEBAQwPR1dXWhsrISUVFR8Hq9SE5OZvpaWlpw7tw5JCQk4MGDB5gxYwbTV1dXh+rqaun8zM/PV3mGIr9ZzqM10tV7nh4nkcfboXWkE3m8/d+H+3l6ZfWeVwQd7kQe/Tw9624AePxEHg+q0are19XVYe/evZg8ebL0VOmTyu/AH09LdmO9icZoLdkZbaLxNCzZKaEPDQ1l+oYqvwLfhF4tE3oT+seR34BvQq+WCb0J/ePKEPiCIKwWBOGGIAg3BUH4M+P77wiC8MNP/1UJgqB+rvAJNJ6gd7lcJvQKDTf0DodjXEMPGKjqC4IQCODfAawA0AbgkiAIh0RRrKNsFgDPi6L4QBCEFwH8A8C84XiD3d3dhhN5LS0thqCvqakxlMgbHBw01CNvcHBQN5FHqvdGEnk2m00TegDYvn07enp6NKG32WzDDn1DQ4Mh6B0Oh6FE3uDgoKFE3sDAgG4ib8+ePbDZbJrQA8Cnn36Kvr4+Teg9Ho8h6K1W67BDf/Xq1RGFHjC2nDcXwE1RFG8BgCAIOwG8CkACXxTFKsp/AcDk4XhzJLnX39+P5557Dq2trVI/N1p0cq+4uFjWR44WnZaaNm0arl+/rvKQLi1kW6upU6dKCS1aJI01ODiI8PBwpKSkSOvTtMj+dP39/UhOTkZERASuXr2q8h04cACAr2vL1KlT4fF4mD6S3Lt//z6Ki4vR29vL9JHkXl9fHxYsWIC7d+9KPexokbFyOBxYuHAhM3kGyJN7M2fOlHIMSpG1csC3Rx29FyERSVUODg5CEARkZGQwx5j8fQYGBhATE4O4uDjmGJMsh81mQ3p6OoKDg5ljcvToUQC+hFxhYSEcDgfTR5J7VqsVs2fPlu0oTIu8Z7vdjsWLF6O9vV12LhKRc8Xr9aK4uFhKJCpFJ/emT5+O+vp6pg/wnVcjmtwTBGEDgNWiKP7qp9dbAMwTRfEjjv8PAAqIX/G99wG8DwAZGRmzWbFKWmZyz5QpvkY6uce6r2J+WgiCsBTALwEsZn1fFMV/wDcNwJw5c3SzwmQDy9/85jfM71ssFln6jOerra3F119/DcDXv+2jj5ifWaiqqpJy35MnT8brr7/OOgapBTYAlJaWYvFi9eGSzS7u3LkDAFixYgWzfRQp5JGr4BtvvMG8JSdzenJH8uGHHzJvoa1WK3bs2CG95o1Je3u71M5Ly9fY2ChLn/F8NTU1OHv2LADfZp3vvfce03fmzBnpKpubm8vcz44k8khX3IULFzJv3Unfe5JnX7t2LXPaQwp5JB3Jm86QOT25GPKOtaurCzt37pRe83wtLS04ePCgrq+urk6WjuT5Ll68iIsXLwLwhdueREbAbwNAT3AmA1A9tSIIwkwA/wXgRVEU7z/Ru/pJmZmZ6O7uRmxsrOp7FosFx44dQ1JSEpxOJ+Lj45m+2tpafPPNN8jMzMTdu3cxbdo0po9AX1BQgPr6euTm5qp8NPT5+fm4ceMGMjIyVD4a+ry8PDQ0NGDSpEkqn9PpxMGDB9HX14fc3Fw0NjYiKSlJ5bPb7VJcNzs7GxaLBXFxcarkntVqxcGDBxEREYG4uDg8fPiQeazt7e04fPgwYmNjERQUhJCQEKavsbERJ06cQGpqKmw2GyZPnsz0EehzcnJgsVgwffp0pq+iogJXr16VxnjKlCnMMS4vL4fFYpHGmPV7aejJGKekpKh8g4OD2LdvHwYHBzF16lTcvHkTCQkJqjl2b28vDh06hODgYKSmpqKtrY15DF1dXSgvL0dUVBQmTJgAh8PB9LW0tODo0aNISEiA1+tFTEwM01dXV4fTp08jPT0d3d3dKCgoYPoI9GTsnrTfvpGq/iUAuYIgZAuCEALgTQCHaIMgCBkA9gPYIooie/I3jCLV+7i4OGzdutVw9Z61nzwwvIk8eluroSzZDXciLykpiel73Op9WFgY0/e41XuWhpLIoze7MJLI27RpE/Ly8pg+ZfV+8mR2iWq0qvcXL17E8ePHUVBQgA0bNjA9Q5Uu+KIougF8BOArANcB7BZF8ZogCB8IgvDBT7b/FUA8gP8QBOGKIAjaj909gR4X+p97Lztznf6RRntbq02bNo35GK5R6J+0rTaRoYd0RFE8CuCo4mv/Sf37VwBUxbzh1s8BvdfrlXlM6E3onzboAT9K7plXerVM6NUyoTemMf1YLtF4g37nzp2GoP/0009N6CkZhb6vrw979+4dt9ADfgB+d3e3Ieibmppw69YtXehJsEUvkUcaI+gl8shyol4i74svvoAgCIYSeQA0oQcgLTtpQW+z2QxB397ebgj669ev4/r167rQX7hwAQB0E3knTpwAAN1E3v79+wFAN5G3Y8cOCIKgCT0AfPLJJwCgCb3H4zEEvdVqNQS9xWJBc3OzLvSkVdtIQg+McfBJcs/tdiMrK0vWHYcWSUuJoojJkydLJ55SdHIvMTFRStXRcrlcstfh4eHS+jQtOrlGtlyiu/wQ0Xv7JSYmqvaiU/rcbjeys7O5KUXis9vtKCoqwo0bN3Djxg2Vj3zAORwOTJs2jZlQAx4l9zweDzIzM6UOL0rRyb1JkybJuu3QopN7MTExzDFW9uoLDAxkjjHdmi06Ohr9/f2aYyeKIlJTU9HZ2cns9Ud8TqcTubm5sFgsUlaAFknu2e12zJw5UzoPlSLJPZfLhezsbGnMlSLJPVEUkZ6eLq3FK0Un95KSkmT76RGRjEFlZaXZc8+UqfGop7bnXk5ODtra2vDHP/6R+f3a2lqp91l8fDw+/PBDpq+qqkratZS3Vk/m9OSKt3jxYuYnKpnTkxz1yy+/jJKSEpWPzOnJVXvz5s3MW3dSyCOdaT/88EPEx8erfKSQZ7fbAQB//etfmbeBZE5P7lz+/ve/M0bEd8eye/duiKKI4OBg/PnPqocuAfjm9IcPHwYApKen491332X6KioqpKt2SUkJXn75ZZWHJPLI7ewLL7yAhQsXqnxkTk8y7bycAz2nB4D33nuP2a+QFPIePnwIAPjd737HvCUnc3qyczFv7MicnnRk5vnq6uqkdGRMTIxmIo9Me/Ly8vDGG2+oPKIo4tSpU9Ld7FO9aWZoaCiio6MREBCg+q+urg4HDhxAZmYmkpKSMHHiRKbvwoULOH36NIqKihAREYGwsDCVRxAEnDx5EpcuXcIzzzwDAAgODlb5RFHEgQMHUF9fL801g4KCVD63242dO3dKLbAB3+2s0jcwMIBPP/0UPT09UiNOlu/hw4fYvn07RFGUAGAd6507d/DZZ58hMjISU6dORVRUFNN38+ZN7N27FykpKUhLS0NSUhLTd/XqVRw+fBg5OTmIjY3l/rzKykqcPXsWxcXFCA4ORmhoKHOMjxw5gitXrkhjzBo7r9eLvXv3oqmpSXOMyQdrZ2en9MHLGjtS5+jv78fMmTO5Y3fv3j3s2LEDgYGBKCgoYP6sgIAAtLa24osvvkBsbCyysrIQFxfH9NXX12P//v1IT09HcnIy13fp0iWcOHECBQUFiIqKkvr4Kcfu66+/xoULF6QxedK5/5gGn6fHTeQ9yQaWdPV+xYoVzCsV8PhLdsOVyAsPD8e7776L6Ohopo9c6Ueqes8q0LGq9yw9TvV+KM/TG03kse64AHn1fuvWrZgwYQLTN5zVe3Klr6qqwpw5c5gdnR9Hfgf+aC/ZjQT0w7VOT6DnVe9HGvqfe8lutJpojCb0vJWPocqvwDehV8uEXi0Ten35Dfgm9GqZ0KtlQm9MY7qqT/S0QH/06FFD0B86dMgQ9J9//rkJPaXe3l7s27dPF/ru7m7s27dv3EIP+AH43d3dhqBvampCU1OTbiKvpqYGgH4ijwRFtKAHfO2cjCTy9u7di6CgIN1E3meffYaIiAhN6AFfSi0yMlITepvNZgh60jJKD3qS3NODngRU9BJ5ZIlVL5F36NAhBAQEaEIPALt27UJISIgm9IAvHRkVFaUJvcfjMQS91Wo1BD0JC+lBT4JWIwk94EcBnhkzZnCr93Riqri4mHtCEugB31oza1AHBgZkvfhKS0uZv7OmpkZKUYWEhHChp99bXFwct3pP+9LS0rjQ076cnBwu9LSvqKiI+zw97Zs5cyYXeiNj7PV6pTV6gD/Gd+/ele1AzBtj+ndGRkZyY7i0LykpiVu9p32ZmZnc6j3tY223xfJNnz5dNyYO+PZQ5EFv5Pykx/hJAjx+AX5oaCgXejqGS9ZA9Xy8P6TX65XaM2n5XC4XHA6Hrs9ut4MeX55PuWXyaPiMjnFkZCR3yY6Ei7R+p8fjkQIyWr7BwUGpzZjRYxgtX0BAACIiInR94eHh3PX3xzk/n9rkXkFBAe7evYvf/va3zO+TOT3g+wRnpcroOT3g+yR95ZVXVD56Tg8AS5cuxXPPPafykTk9aRS6bt06zJql3kaAFPJIxr2srIy7geUnn3yCvr4+AMDHH3/MbCFOCnnkA+dvf/sb8yQic3qSKuOdHGROD/hSj7w+hGRODwD5+fnMfe/oOT3gm0atXr1a5aPn9ACwatUq5p6BZE5PnsHYtGkTs18hKeSRDri//vWvVZtmAo8KeeSD6fe//z0TLjKnJ6lH3tiROb0oiggJCcFf/vIXpo/M6QHfXdyvfqVuWUHP6QHfHcG6detUPq/XiwMHDkjPB7DOzaFoTFf1AwICDIVzUlNTmbeoykJeVFSUbiHvhRde4L4fupC3bNkyro+u3mv9gejqvdaTffSjtbzbYkBeyCsqKuJeOehCXlZWFncaQBfyEhISDBXyQkJCdAt5WmNMF/J4rbcAefWe1eyUiK7ez5vH3+qBLuQVFxdzr8x0IU9rGkAX8tLS0pjTAGUhLzo6mnt7T6DXGruhaEyDz5M/JfJ+riYayur9cPfIG8uJPKNLdkYSeVqbsYxG9V4JvdaH3FDkd+CP9SW7sQD9075k5++dc0YbesDPwDehV8uEXi0Ten35Dfgm9GqZ0KtlQm9MY7qqT2QEegCGd60dLej3799vCPojR44Y6pG3d+/ecQP9qVOnDEF/6tQpQ9AfPHhw3EIP+AH43d3dhqAnyT096En4QS+Rd+7cOUOJvJMnTxpK5JWXlxtK5O3fv1/a8ZcHPeBLqaWmpmpCb7PZDEFPknt60JPmI1rQezweqWWWXiLvzJkzhhJ5x48fN5TIO3jwoG4iD/ClKBMSEjSh93g8hqC3Wq2GoG9ubkZzc7Mu9GRT0JGEHvCTAA/g2zpKC3qinJwcpsftdoPepJPne/jwodQNR8tH/87AwEBuIo/2RUVFcaGnfQkJCVzoaV9aWhq3ek/7MjMzudDTvilTpnCr97du3ZJe88bE6XTK+gTyfLdu3ZIFm4yMcWhoKDeRR/smTpzIrd7TvuTkZC70tC89PV0zJk6UlZXFXQI0cn56PB7ZTsU8n8vlknr4PbUBHqLU1FRZE0da9LbEaWlp0vbWtLxer9SeScvncrlk0PN8dEINAFJSUpg+5dbKMTExTB8dXwV8JznLRx8rAVnPFx8fD7fbLUvBsXypqamyNCKRKIqy98cbE4/HI9uGm+dzOBwy6Hk+EmgiSkhIYProJqCAL1moNyaAb/z0fNHR0fB6vbq+pKQkuFwuVaNWpW/SpElPfH52dXWpvv44GtPgFxYWoqurC++//77qe6SQRwY2JycHmzdvVvm8Xi/27dsnDWxpaSlz7zYypyfiJffInJ6Il9wjhTwiXnJPmcjjJfdI33uSyPvzn//MvMLU1NRIYxISEqKZyCO+1NRU7hgfPnxYAr+wsBCbNm1S+dxuN/bs2SOBP3/+fGZugszpif7/9s4vtoq0DOO/h60NEKxbOGgMtrKVFahGyCq6boJBDdmFXhCSvTCaLG5MNhI1Xla9UBMvXO/UqItks9l4s3ujMUgX5I/UkizQVdKlVLIEl2hqTfQosWFJ4HD6ejFnyvT0fDNf2/Nn2vP9koZzmJczz/mYpzPzfu98r6tyL07kxbgq9+JEXnxScFXuVVfkuSr3khV5EK3NV4tLly7Njl1XVxeHDx+eFxPf08dxmzdv5tChQ/Pi4nv6+PjcuXMnBw4cmBdXKpVmG6fCCq/ck+SVva+1wCI8MH3c7MJ1aZdM5B08eNCpJ5nIczV+hLnZ+1qlqzHJ7H3aUsnJZhdpFX7JRN6OHTucl6jJRF5fX59zjI8fP87ly5fZvXs3GzdurPlZsemvX7/OwMCAc5/JRF6tktSYZPZ+YGDAGZfM3u/du9cZl8zep5klmcjbtWuX87I9mcjr7+93VowmE3k9PT1eiTzX7V1s+ps3b6Yenwsh18avhW9FXrXpV0pbK9+KPNeB65u9T5o+rSIvafqV0tbK9curFdn7atPHC4YulWVlfN95+pVq+mZN2VWb3jVltxJNn6cpu0aZHpaR8YPpg+mD6evHsjB+MH0wfTB9ffHK6kt6CvgJ8BDwopk9X7Vdle37gTvAl82sdiOxBVJv05dKpbqb/uzZs16mP3/+vJfpT58+7WX6oaGhppq+XC63zPRDQ0Nepj937lym6a9du8aJEyfa1vTgYXxJDwE/B/YCk8Abko6ZWbKT4D7g0crPp4AXKn8uCTOjWCxSLBYzTR83P8wy/fj4uFdF3ujoqFdF3sjIiFdF3pkzZ7wq8k6ePJlZkQd4mf727dtepp+ammJqairT9HGj0DTT37t3j7GxMVatWpXZtfbChQteFXnDw8NeFXmnTp3KrMiDaJ3Enp6eVNOXy2Uv09+6dcvL9JOTk0B6RV6pVGJiYgKgoaYHj8o9SZ8Gvm9mT1befxvAzH6YiPklMGxmr1TevwXsMbN/1vhIYOGVe4VCwTmoyaIb17RTqVSa7Z2WFjc9PT2nkMUVV1044hPX2dnpNH0ybt26dU7TJ+O6u7udpk/GbdiwwWn6YrE4+941xuVyeU4xkuu73r17l+np6cy4xYydJKfpk3Fr1qxxmj4Z19XV5TR9Mm79+vXOGZJknGvsqsc4bVo02UXYFRfvd/Xq1QwODjpjYpZSubcJSPZrnmT+2bxWzCZgjvElPQc8B9Db2+ux64jt27c7Vxs1s1nj9/f3Oz9jZmZm1vhpcV1dXbMllmlxd+7cma3gS4tLHhxbtmzJjFu7dm3q2MRxfX19zqm9ZNy2bducpc7A7EGZNcax8bPGODZ+WlxHR8dswYrv2NUq4KkVl3Z1FscVCgVnK7Jk3NatW1N71MVxaWMHD8Y47bvGVw9ZcfF+XSXivvgYv9Y3qr5M8InBzI4CRyE642ft2KcWORAILByfrP4kkLwB+wAwtYiYQCCQE3yM/wbwqKRHJHUCXwCOVcUcA55RxOPA/9Lu7wOBQGvJvNQ3s/uSvg78nmg67yUzm5D01cr2I8BrRFN5N4im855tnORAILBUvObxzew1InMn/+5I4rUBX6uvtEAg0CiWReVeIBCoL8H4gUAbEowfCLQhwfiBQBvSssU2Jf0b+FtmIBSAYmZU68izvjxrg3zry7M28Nf3QTObV//bMuP7IulPtWqN80Ke9eVZG+RbX561wdL1hUv9QKANCcYPBNqQ5WD8o60WkEGe9eVZG+RbX561wRL15f4ePxAI1J/lcMYPBAJ1Jhg/EGhDcmN8SU9JekvSDUnfqrFdkn5a2X5F0mM50valiqYrkl6XNL+nVgv1JeJ2SSpLejpP2iTtkTQmaULSH5ulzUefpPdI+p2kNyv6mvbkqaSXJP1L0lXH9sV7wsxa/kP0uO9fgT6gE3gT6K+K2Q+cIFrt53HgUo60PQF0V17va5Y2X32JuD8QPWX5dF60AQ8DfwF6K+/fm6exA74D/KjyeiPwX6CzSfo+AzwGXHVsX7Qn8nLG/yRww8zeNrN7wKtAdefAA8CvLOIi8LCk9+dBm5m9bmbxSokXiVYgahY+YwfwDeDXQH3ardZP2xeB35jZ3wHMLG/6DHh3ZQn5dUTGn996uAGY2Uhlfy4W7Ym8GN+1WOdCYxrBQvf7FaLfws0iU5+kTcBB4AjNxWfsPgx0SxqW9GdJzzRNnZ++nwHbiZaSGwe+aWYzzZGXyaI9kZc22XVb0LMBeO9X0meJjF974fTG4KPvx8CgmZXTVoNtAD7aOoCPA58H1gAXJF00s+uNFoefvieBMeBzwIeA05LOm9n0vH/ZfBbtibwYP88LenrtV9LHgBeBfWb2n+rtDcRH3yeAVyumLwD7Jd03s9/mQNskUDSzd4B3JI0AO4BmGN9H37PA8xbdVN+QdBPYBow2QV8Wi/dEsxIpGUmMDuBt4BEeJFk+UhUzwNxExmiOtPUSrTf4RB7Hrir+ZZqX3PMZu+3A2UrsWuAq8NEc6XuBqKEMwPuAfwCFJv7/bsad3Fu0J3JxxrccL+jpqe27wAbgF5Wz6n1r0pNdnvpago82M7sm6SRwBZgh6s1Yc/qqFfqAHwAvSxonMtigmTXlcV1JrwB7gIKkSeB7wLsS2hbtiVCyGwi0IXnJ6gcCgSYSjB8ItCHB+IFAGxKMHwi0IcH4gUAbEowfCLQhwfiBQBvyf9OLaeF/omOLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Defining the elements\n",
    "\n",
    "* We need to define the elements that will be used to construct the solution\n",
    "* We'll use P1 elements (continuous polynomials of order 1 on triangles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "P1 = FiniteElement('CG', 'triangle', 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Defining the discrete function space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Now we create the function space that contains the basis functions on the mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vh = FunctionSpace(mesh, P1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This function space is used to create the solution $u$ and test function $v$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = Function(Vh)\n",
    "v = TestFunction(Vh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Defining Dirichlet boundary conditions\n",
    "\n",
    "Typically done in two steps:\n",
    "1. Specifying the boundaries where a Dirichlet condition is to be imposed\n",
    "2. Specifying the value of the solution at these boundaries\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In our case, we do these two steps using the code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "A Python function `boundaries` is used to define the boundary $\\partial \\Omega$\n",
    "\n",
    "`DirichletBC` is used to specify that functions in $V_h$ take the value of 0.0 on $\\partial \\Omega$\n",
    "\n",
    "The boundary conditions are stored a variable called `bc` - this will be passed to the solver (which is when they are imposed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Defining the PDE via 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",
    "\n",
    "The weak form is coded up very naturally (the key benefit of FEniCS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "f = 1.0\n",
    "F = -dot(grad(u), grad(v)) * dx + f * v * dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quantity `dx` represents integration over the domain $\\Omega$ defined by the mesh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Solving the problem\n",
    "\n",
    "We can now solve the problem by simply calling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "solve(F == 0, u, bc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* This uses a Newton-like method to solve the problem $F(u) = 0$\n",
    "* Jacobian automatically computed!\n",
    "* The boundary conditions are passed as the third argument"
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}