{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving a Laplace problem with Dirichlet boundary conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, we will solve a simple Laplace problem inside the unit sphere with Dirichlet boundary conditions. Let $\\Omega$ be the unit sphere with boundary $\\Gamma$. Let $\\nu$ be the outward pointing normal on $\\Gamma$. The PDE and boundary conditions are given by\n",
    "\n",
    "\\begin{align}\n",
    "\\Delta u &= 0&&\\text{in }\\Omega,\\\\\n",
    "u &= f&&\\text{on }\\Gamma.\n",
    "\\end{align}\n",
    "\n",
    "The boundary data is a source $\\hat{u}$ located at the point $(0.9,0,0)$.\n",
    "$$\n",
    "\\hat{u}(\\mathbf x)=\\frac{1}{4\\pi\\sqrt{(x-0.9)^2+y^2+z^2}}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, we will use a direct integral equation of the first kind. Let\n",
    "$$\n",
    "g(\\mathbf x,\\mathbf y) = \\frac{1}{4\\pi |\\mathbf x-\\mathbf y|}\n",
    "$$\n",
    "be the Green's function for Laplace in three dimensions. From Green's representation theorem, it follows that every harmonic function $u$ in $\\Omega$ satisfies\n",
    "\n",
    "$$\n",
    "u(\\mathbf x) = \\int_{\\Gamma} g(\\mathbf x,\\mathbf y)\\frac{\\partial u(\\mathbf y)}{\\partial \\nu(\\mathbf{y})}\\mathrm{d}\\mathbf y-\\int_{\\Gamma}\\frac{\\partial g(\\mathbf x,\\mathbf y)}{\\partial \\nu(\\mathbf{y})}u(\\mathbf y)\\mathrm{d}\\mathbf y,~\\mathbf x\\in\\Omega\\setminus\\Gamma\n",
    "$$\n",
    "\n",
    "or equivalantly\n",
    "\n",
    "$$\n",
    "u(\\mathbf x) = \\left[\\mathcal{V}\\frac{\\partial u(\\mathbf y)}{\\partial \\nu(\\mathbf{y})}\\right] (\\mathbf{x}) - \\left[\\mathcal{K}u\\right] (\\mathbf{x}),~\\mathbf x\\in\\Omega\\setminus\\Gamma,\n",
    "$$\n",
    "\n",
    "where $\\mathcal{V}$ and $\\mathcal{K}$ are the single and double layer potential operators.\n",
    "\n",
    "Taking the limit $\\mathbf x\\rightarrow \\Gamma$ we obtain the boundary integral equation\n",
    "\n",
    "$$\n",
    "\\left[\\mathsf{V}\\frac{\\partial u}{\\partial n}\\right] (\\mathbf x)=\\left[(\\tfrac12\\mathsf{Id}+\\mathsf{K})u\\right] (\\mathbf x),~\\mathbf x\\in\\Gamma.\n",
    "$$\n",
    "\n",
    "Here, $\\mathsf{V}$ and $\\mathsf{K}$ are the single and double layer boundary operators."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now demonstrate how to solve this problem with Bempp. We begin by importing Bempp and NumPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bempp.api\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next define a mesh or grid. For this problem, we will use the built-in function `sphere` that defines a simple spherical grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = bempp.api.shapes.sphere(h=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define the spaces. For this example we will use two spaces: the space of continuous, piecewise linear functions; and the space of piecewise constant functions. The space of piecewise constant functions has the right smoothness for the unknown Neumann data. We will use continuous, piecewise linear functions to represent the known Dirichlet data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp0_space = bempp.api.function_space(grid, \"DP\", 0)\n",
    "p1_space = bempp.api.function_space(grid, \"P\", 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now define the operators. We need the identity, single layer, and double layer boundary operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "identity = bempp.api.operators.boundary.sparse.identity(\n",
    "    p1_space, p1_space, dp0_space)\n",
    "dlp = bempp.api.operators.boundary.laplace.double_layer(\n",
    "    p1_space, p1_space, dp0_space)\n",
    "slp = bempp.api.operators.boundary.laplace.single_layer(\n",
    "    dp0_space, p1_space, dp0_space)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define the GridFunction object on the sphere grid that represents the Dirichlet data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "@bempp.api.real_callable\n",
    "def dirichlet_data(x, n, domain_index, result):\n",
    "    result[0] = 1./(4 * np.pi * ((x[0] - .9)**2 + x[1]**2 + x[2]**2)**(0.5))\n",
    "    \n",
    "dirichlet_fun = bempp.api.GridFunction(p1_space, fun=dirichlet_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next assemble the right-hand side of the boundary integral equation, given by $$(\\tfrac12\\mathsf{Id}+\\mathsf{K})u.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "rhs = (.5 * identity + dlp) * dirichlet_fun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now solve the linear system using a conjugate gradient (CG) method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "neumann_fun, info = bempp.api.linalg.cg(slp, rhs, tol=1E-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now want to provide a simple plot of the solution in the $(x,y)$ plane for $z=0$. We first define points at which to plot the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_grid_points = 150\n",
    "plot_grid = np.mgrid[-1:1:n_grid_points*1j, -1:1:n_grid_points*1j]\n",
    "points = np.vstack((plot_grid[0].ravel(),\n",
    "                    plot_grid[1].ravel(),\n",
    "                    np.zeros(plot_grid[0].size)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable `points` now contains in its columns the coordinates of the evaluation points. We can now use Green's representation theorem to evaluate the solution on these points. Note in particular the last line of the following code. It is a direct implementation of Green's representation theorem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "slp_pot = bempp.api.operators.potential.laplace.single_layer(\n",
    "    dp0_space, points)\n",
    "dlp_pot = bempp.api.operators.potential.laplace.double_layer(\n",
    "    p1_space, points)\n",
    "u_evaluated = slp_pot * neumann_fun - dlp_pot * dirichlet_fun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now plot the 2D slice of the solution. For a full three dimensional visualization, Bempp can export the data to Gmsh. Since the solution decays quickly, we use a logarithmic plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f4835f5a310>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUUAAAEICAYAAADIsubvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO29ebwkVX3w/f1V32VWGAaGGYZFXFDBBaIjStQAEQnwJBJNNBqjxJiH+D4hm3FBUUPcgqJiEk2U5PGNeaNREyUQRdmioiYaEEFAFkGHOAwMM8Aw+723u37vH1XVfar6VNWp7urtzvl+Pv251afOVtW3f/3bzilRVTwej8cTEYx6Ah6PxzNOeKHo8Xg8Bl4oejwej4EXih6Px2PghaLH4/EYeKHo8Xg8Bl4oeoaOiGwUkdNq7lNF5Ek9tn2hiNxV53w8k4sXimOEiPymiNwoIrtE5AER+aqIvGDU8zIRkVNEZNOo59EPWQGqqt9S1aeMck6e8cELxTFBRN4IfBR4P7AWOAr4G+DsUc7L49nf8EJxDBCRA4F3A7+vql9S1d2quqCq/66qb47rzIrIR0Vkc/z6qIjMxudOEZFNIvIWEXko1jJ/VUTOEpG7ReQREXm7Md6FIvKvIvJ5EdkpIjeJyPHG+ZQmJSL/ICLvFZHlwFeB9bE2u0tE1otIICLni8i9IvKwiHxBRFYb7V8jIvfF5y4ouRdniciP4nndLyJvMs79bxG5J76eK0RkfU4f3xCR3zXe/7aIfDs+vj4uviWe/29ktV8ROTbuY7uI3C4iL8nci4+LyFfiOX5PRJ5YdE2eycILxfHgJGAJcFlBnQuA5wEnAMcDJwLvMM6vi/s4HHgX8HfAbwHPBl4IvEtEnmDUPxv4F2A18Fng30RkumiSqrobOBPYrKor4tdm4A+BXwVOBtYDjwIfBxCR44C/BV4TnzsYOKJgmP8L/J6qrgSeDvxH3M8vAn8BvAI4DLgP+FzRfHOu4Rfiw+Pj+X/ePB/fg38HrgYOBf4A+IyImOb1q4A/Bw4C7gHeV3UenvHFC8Xx4GBgm6o2C+q8Gni3qj6kqluJvpSvMc4vAO9T1QUiYXEI8JequlNVbwduB55p1P++qv5rXP8jRAL1eT3O//eAC1R1k6rOARcCvy4iU8CvA19W1evjc+8EwoK+FoDjROQAVX1UVW8yrv9TqnpT3M/bgJNE5Oge55zH84AVwEWqOq+q/wF8mUgQJnxJVf87/rw+Q/RD5VkkeKE4HjwMHBILkTzWE2lHCffFZe0+VLUVH++N/24xzu8l+rIn/Cw5UNUQ2JTprwqPAy6Lzc3twB1Ai8g3uj4z1m6i683j14CzgPtE5JsiclJcnrp+Vd0V93N4j3POYz3ws/ieJNyXGedB43gP6fvqmXC8UBwP/gvYR2SC5rGZSPgkHBWX9cqRyYGIBEQmbdLfHmCZUXedcWzbVulnwJmqusp4LVHV+4EHMmMtI9KMrajqDap6NpHp+m/AF+JTqeuP/ZsHA/dbutldMP8yNgNHxvck4aiccTyLEC8UxwBVfYzID/jxOECyTESmReRMEflgXO2fgXeIyBoROSSu/099DPtsEXlZrJ3+MTAHfDc+dzPwmyLSEJEziHyFCVuAg+PgUMIngPeJyOMA4jkmUfN/BX5ZRF4gIjNEASXr/52IzIjIq0XkwNis30GkcULk93ydiJwQB5jeD3xPVTdauroZeFl8H58EvD5zfgvwhO5mAHyPSKi+Jf4MTgF+hR78l57JxAvFMUFVPwK8kSh4spVI+zqPSFsCeC9wI/BD4FbgprisVy4HfoMoKPIa4GWxIAL4IyJBsJ3Il5fMAVW9k0hA/yQ2l9cDfwlcAVwtIjuJhOtz4/q3A79PJNQeiMcrynN8DbBRRHYAbyAKFqGq1xH5I78Y9/NE4JU5fVwCzBMJv08T+f1MLgQ+Hc//FeYJVZ0HXkIUUNpGlBb12vi6PfsB4jeZ3f8QkQuBJ6nqb416Lh7PuOE1RY/H4zGoRSiKyKfipOHbcs6LiPxVnHj7QxF5lnHuDBG5Kz53fh3z8Xg8nl6pxXwWkV8AdgH/qKpPt5w/iygJ9iwiX9NfqupzRaQB3A28mMjPdAPwKlX9Ud+T8ng8nh6oRVNU1euBRwqqnE0kMFVVvwusEpHDiFZl3KOqP4kd3J/Dr/X1eDwjpChZuE4Ox0jgJdIKD88pf66tAxE5FzgXYPny5c9+6lOfOpiZejwevv/9729T1TVV2/3Sqcv14Uda5RWB7/9w7ipVPaPy5AbMsISiWMq0oLy7UPVS4FKADRs26I033ljf7DyFvPjn38PuI5YytTda5LHjqCmCeEHiwjJhyaMaH8PCiugjXfKw0lwaf7yJPRJ/so19yr410bnpXdDYG52YWy1MxWtxwmlYsakV9xuwctM8X7/Gu5yHhYjcV16rm4cfafHfVx3lVLdx2I8P6WWMQTMsobgJY1UDndUTMznlHo9nAlEgLFzaPv4MSyheAZwnIp8jMo8fU9UHRGQrcIyIPJ5oGdUrgd8c0pw8Mac9P8oB33b8sramtrBCmNkZaXDN45fTmhGm9kYqn4QQNqJ6jTmluTTuSGBqT3TYmgWJg3iqscYYa4qtJcLU7s74rSUS99XpV1qwZ230prkM5lbPcsL/+QgA8wcK03H75jI4+NYo5/z6r7yl73vh6Q9FWVA383lcqUUoisg/A6cQbWqwCfgzYBpAVT8BXEkUeb6HaF3t6+JzTRE5D7gKaBDtgnJ7HXPydHPKGR8AQANh55HRR9+aFWaOiZYJB81ImAEEC0prJm6o0JhX1AjLSSzgVKJXpzw+IaBiCEOh4yyRdPt2uab9KYmADOZpzxWgMd+ZpzRh+5OiHc+e+s5LaMR1V9zf0Va+909/WnhfPPXiNUVAVV9Vcl6JlnrZzl1JJDQ9Hs+EoyitCV8lNyzz2TMCTntBxyxuzEF4dKRRNeZpm7LBvNIytpYV8/85GwYz3qstRIZ53lQfHdtXOJ/MXzQKykD6uvasDdra5NPeegmtGTj4R5FZ950vvgnP4AjtsdKJwQvFRcTpz7mQ7U9dCcDcQQGzT4icfY19oEHkp4PouC1UQuw5AOQLvq7ylDCzCMPscZY8YejQPmt+J4SNzvUikcDccXRkjz/5PZcwuz06teqeJt+6/M0Fk/NUQYGWF4oej8fTwWuKnpHygpdezO51kQY0fexKwqk4GLFPaU131CvRjhaWa/q6aGyFbXK0xAx5fbuO6VTPos0mgaPGvs7xI0+Z4ikXXgLA0oeUm//mjfkDe0pRYMH7FD3D4pdW/jYA+15wLHvWRo608OBGOzobTklbYBSZxVEF+3GRrzBfmJU4GG3CuIow7LV9zrTEWDagUxDEu0juO1g47u2RgFyyTVlzQ2RjX/WDd9s78nShqDefPR6Pp41Ca7JloheKk8AZx76N3U9ezWO/Ez2aeWaHtqPEYqSE9aLl9R0AKWjj0r4wil21fVkdy1ipHMuw837vocK9r1oFwPF/cAkHblzg+n/3yeFlRCtaJhsvFMeY5/7WhwFoPW8NQVOZ3mUkRlfopyx9xlrPuU15tLlIKKWo0L6SMCwqz5mDCkztjt7MHwBbnj3NsRdEpvUd7/uTgg72d4SW6z/PmOKFosfjqY0o0DJ4oSgiq4HPA0cDG4FXqOqjlnobgZ1ED0BrquqGsr69UBwzTj4zWoq384hpgpnOMjltSMdUNnPzqNm0tWDX1HrIR7T16TqvGrTEylH3ZCliECeJxxHrY973EZY9EFW85WNeazSJ8hSHoimeD1ynqhfFO/afD7w1p+6pqrrNtWMvFMeAF77kYgB2r20wvTr6SBrzpL7MzonUeedc/0+dhJJbZwMVanWZzFXmlWx9Ni/MHRQdP+kDH2HplqjCrR/xAhIgHIKmSLQZ9Snx8aeBb5AvFCvhheKIOPMJ0SYFu56+jrkDjd1npvv4hxpUAMNxzKEJtToEYh99qxBvdxKtmpmPn4D99Ddfwqp7oo0mv33Z/rlKpqKmeIiImBujXhrvm+rCWlV9ACDecevQgildLSIKfNKlfy8UPR5PbShCy/0pJ9uKfHwici2wznLqggpTer6qbo6F5jUicmf8+JRcvFAcAS98ycXsPflwAAIjqSso24auF03ObNerEtrW9CpEmh3LK7UdhcmcN35S1JJ2+fwBsOXEaHXR086/hNsv2j/N6brMZ1U9Le+ciGwRkcNiLfEw4KGcPjbHfx8SkcuIngtVKBT9c5+HzAteejF7D24gYZIXJ+1XNs8m5UssM4cd/I+Ffdvap8qH4idyox+TueQ6C+tYxknVE40+Q4W5VdpeHbM/oQjz2nB69ckVwDnx8TnA5dkKIrJcRFYmx8DpgPUxzCZeKHo8ntqIkrcDp1efXAS8WER+TPSI5IsARGS9iCT7s64Fvi0itwD/DXxFVb9W1rE3n4fAyWd+gF3rI8/81LKAoKl2DcxRGRvU6hJnRmH+9hI4cqnjOnbuOU2XJ/n1LWguhSd+KHqEwuwjwo/ev3+Y08NIyVHVh4EXWco3E+3yj6r+BDi+at9eKA6QX3zRXwCwb91MZ/t9R6uht9Ul1du4tG8L7lFFjC1C1Sb8ygRiz/5Da/uMMLT1NaXtNW8Ly5VnvDEypxdz6o6q0NLJNkBrmb2InCEid4nIPXEiZfb8m0Xk5vh1m4i04ox0RGSjiNwan/PPLfV4JpwQcXqNK30LRRFpAB8HzgSOA14lIseZdVT1YlU9QVVPAN4GfFNVHzGqnBqfL12CMymc+NoPs2v9DLvWz4BA0NTIbI5pB1egyzzLIzc4UEZpEKViNLtXnAI6BQEQOnXMPnPnn+k3VW5pr0L3mFW1ROOlQjugRgB71yp71ypP/OBHLJNdHESBlimn17hSx8xOBO6J7Xfix5ieDfwop/6rgH+uYdyx48wnvZmHTjkMiB5NmJD9whZFckdhNqewmc3Zfuv0BxqUmbeV+uzJJC5rrz20Md4bZY//aLTZx6E3Lq6nDSaBlkmmjtkfDvzMeL8pLutCRJYBZwBfNIqTjPPvi8i5NczH4/GMkJaK02tcqUNTtF1d3s5WvwJ8J2M6O2WcxwLzXICjjjqq3znXyi/93LsA2P2MQ5naG2sTrj83g9QG+9DoKj2Nr6bgShdVx6swl2FpiKlzRpNgPjqx/cnCz70hMqd/8InJfxRCxRUtY0kds98EHGm8PwLYnFP3lWRMZzPjHEgyzrtQ1UtVdYOqblizZk3fk66L05/7bractIotJ61ifmVg9/UNUPC5fFFHTq9+RDrnC9tYyI1eVxGIok4CsdAnmfSXnWBcb2FFyCPPavHIs1oc/4eLI9k71MDpNa7UMbMbgGNE5PEiMkMk+K7IVhKRA4GTMTLPe80493g840m0IUTg9BpX+jafVbUpIucBVxHFFz6lqreLyBvi85+Iq74UuFpVdxvN1wKXSWSqTQGfdck4HwdOOf0iABaOWMbMzpxopAVrzl8JTkv9KravZDY7t7WXd1HR3C5sk6nTS39FJnO6rGQuWQ2yYDyzXWN3JCD2HKYc+454d+/3TmYuoyIs9L+Eb6TUEhdX1SuBKzNln8i8/wfgHzJlPWWcj5pffNFfsPvIWQCCBU1t/tq3CezIOJnNlf2IWRzbuKxWKerTaWVNWbpN3rlsXy7zyFRuLVEWVkb/TMdecMlEPvZAFZ+8vb9xyukXsfOIWaSlSEvTuWkFpPISi6jgNysjb0ODXrXEXvMka/MJ2ugSaqSFlLid60sgmn3Fr1R5/Erfv9hfmTkvoSChsHd9i6f8+ST6GN0St8c5eXt8Myg9Hs/EoUy+puiFoiOnP+dCAJpHLk9t6GDirNW5mpUubYb0g9tPCo5rX5X8iGX3sErbsnXMruZ30bnMOHY/ZtIJSFNoLlWe/qZIW7ztQ5NjSo9zEMWFyZ79kPilE97J1mcfwNZnH8DCsviWte2kHHoRVv0KuBJhlZ8Sk7PksIKJnBorb7ycemVmda9mtvV6Xcex1M+dM3aBmE7V0YJzln4FWstCdj5jjp3PmOMZfzIZprQihOr2Gle8pujxeGojesTpZIuVyZ79gEkeLrXnGeuY3p1RCyv80NW1Z+KoflyrBGusOM6/l+BK13vbWIUaoZZHkrMasxSn26TGSDTEkv665hYHYWRP9BXds14nZOWLDOsRpwPDC8UCtsbPUZnap+3np9i+2MNKw0lR0cfXT15i3/NwbFPaR54fzqGvxGTtvM9ZqVI0hmEul7fJEYaO82rXidO9WstCHjmx7CE+o0dhrFeruOCFosfjqZVJ1xQnW6QPkOe++sM0FpTGgqJBSdAgh75zE7Nm1pBxynMsa5NnOlcJeljGKO2rUIO0uEJs42cCIKnYWvZ8u42mtD1rQMV4pdqYeYvJGEH00oZCM4BmwOMuvdh6X8YBVRnK2mcRebmI3C4ioYgUPSa1cANsG15TzPCik98PgBw+WxhcTuglpaYXKq8aKTHbBmU22yj6ManFj2jrK1MnN/WmyJQ1y7P9FpnMtr4c23TPxzivnfInvzfyL979jvHyL0aBlqEs87sNeBnwybwKxgbYLybauOYGEblCVfP2egW8UPR4PLUynGe0qOodAFJsiVXdABvw5nOKk8/6IHvWzbBn3QwYFo0zzmZ19TbDpjTP0bGNFVfzO2/8QtO4u8yqJWbN1DzNrqTMZjIjWM3lrrzFzDK/7vo55xXm1rSYW9NqPy1wXIgCLc55ioeIyI3Gq+5Npp03wDbxmqLBnjVTSBj9cwd5gT5nwVdd2jlFjntlBKaze1pM+ZC9RKKd/JQ2YZjtP3lfxWTO/O2KMOfOucCcTggUpqJ6rUPHLyJdYUXLtqLnMonItcA6y6kLVPVyS3lXF5ayUjXHC0WPx1MbyYqWWvpSPa3PLqpsgN3GC0XgF375gwBMLwsIp9w+0Eqfe11an4smWaeG6dhXX2Zz0Rhl713aANYE7by5ZDVES5+2ZXud8nR/pQGV1PVr7mec6kOAMDqpCwFH/3X0EKyNfzAeD8AaowdXtTfABu4n2gD7N8sajc3sR8XJZ36AuVUN5lY1CBvp/8IyH1qbvC+ZYxvncerE4lOrNJeSa3Y2nR3GdBKWeYKrxBzt8i8m5dn6kvEHYiunq9zmb8z6JlXU7j/EaB/kbFM3FcKBC3DgQls4jhJVWAgDp1c/iMhLRWQTcBLwFRG5Ki5fLyJXRnPRJpBsgH0H8AVVvb2s7/1WUzzzyD8CYN/JR6U2ic3FUWj14kvMG6eWNByjvNe5OednVrl/VfvKE6x5mpmtbkk9JMfhZKtfoB12yvO1wy7fofm3XZ7ne0xrlJrI3IPm2xtH3HrJaHbVicznoUSfLyN6plO2fDNwlvG+awPsMvZboejxeAbDpK9o2W+F4t6nrQdAzOCdszboOEiRv2wM6cU3WBjhLTGbrWW2vjN/B+lHtJV3+QWLNE1L3e7+ihK287RLh/JQ2PnE0Uajk5ScSaYWPbdsKY2InCIij4nIzfHrXa5tB8Ev/uJfMHdgg7kDG9XyEOtiBMKylzScvGV+XeQJlrwxHK65Nj+ipa2THzHPlM0KRKN+mQ+x0HeYGaPzMvIZ4yV/ZrkI7RfQ8S/+zYcYDcNZ5jdI+tYUKyyl+Zaq/nKPbT0ez4Qwzs9fcaEO87mnpTQ1tO2Z3etn2r/2QQurd703E7nPf4aSQEmvAZZe51BE6f0pMZ3zzM8qfVjb5NUtMnnLxsumyph1bSZzcj5rMud9NlnN0ixPmd+d/sUyNxElaMSLD9bsZRRE0efJfsRpHTqs61Kak0TkFhH5qog8rWJbROTcZDnQ1q1be57si05+P405JWhC0CT6RlQ1/1wY9x/LjM+sUMhZ/GJd58gxeV38egZFy/zy3ltTZczx8uYllt/DPFM4NU6mXt6YWZPZNv8uk9myY05c3jaVg3QdEUXiY1VBVVjYN8XRHxu+Cb0YHkdQh1C0XV32f+0m4HGqejzw18C/VWgbFapeqqobVHXDmjVrep6sx+MZLJP+iNM6hGLpUhpV3aGqu+LjK4FpETnEpW3d7DxqFm30/oH08gPXt0nbJ9aHUpW26e1cbgCkKmV95Ghq1kRqMhqebayie2PRENPn7Zs7aJ7mmtsmP9Aikm4jgUYvo43pvQmmQ5as351zQYOj4oYQY0kdPsXSpTQisg7YoqoqIicSCeOHge1lbeviBS+LNuYMlgS9f+FN+vxMK/sEHer0Inz7WqJn66OX+gXCyK1eTgqB4R7IllfxIxamEnX5/XLmbDOZ220sPsW4LM932Dkmdf3t6qKEofDEz70PgHtfeUH3RQyIcY4su9C3UFTVpogkS2kawKdU9XYReUN8/hPArwP/j4g0gb3AK1VVAWvbfufk8XhGg6rQ3N+FItiX0sTCMDn+GPAx17Z1c+ZRf0zz+bGVPoC8xL6jzsOmDm24onboMmaelu4aSbaZzanxSzRSa7DG0l6z2p11zto9Zp5GaNEasxpgElDpHHfam9phdv5T01Ey91MuvIS7LhzO0r9xNo1d2C9WtOx78jqCZvTPog0ZrIk8Rv8PfT9atcI5l/qV0nBqMJ1zhWzROLb2QmF/Zl9qE3zJPB1M5qy5LBbhawpMwazTaZ+UhfFuOvNHz+dMvl4Ww4qWRS0UX9x4BQA7Xvc8pvfE/0RhplINwsKpTZ7vr1eG5UssOe+25K7CXBz6zhW2vWiJqfbFfsROf3Y/ombLk7rZObhoh4lAtPgO8wQhRrnE7RuNSFM87KAdnPD70S7dN398sM918ULR4/F4YurcZHZULGqhGD7/eACmd6cTcPumzs/cReNz0V4qjtVXXQftsNR0rqJRWutWjDhD9+NJy0iZxWDX+jImc2F/biYzdLTEMt9hohF2H6c1yId3L2PXzzUdLrp/xjkH0YVFLRT3rp0FIGgpGkQfVJ0+s6jNBPwDDMCfWJcy0HN6jtnewddXOL5t95ter6/Ij5jbJmMyt48754p8h+njtIAU48fjgEN39XhR7qhCs88NZEfNohaKHo9n+Ey6+TzZIr2AU0+7iHBaCKelrSVWopfPddT/CylT3HEyZWZyPyZ3FTPVQRPNXeNs9JO3jju7LVjpem/SdfMeR9q1FVi7bdlqFXOlCunAiqEltk8ba6JNLTAIlCAI45car5BAlEb8mgpCZqebzE43eUKc0D0I/NrnMWb3+unomzCK/RIZ8NK+PD+kI1WX6RWeswiznvyJtnKbsCkYx7Vtul5mn0Sb4MwTwtmx2v2ZArB4c4fuXMSMQGxv9tBpYwo/EXNPRe28gMB43wiUQKLX4w59pOSm9EeyKUXZqx9E5OUicruIhCJS9JjUjSJya7yP640ufXvz2ePx1MqQAi23AS8DPulQ91RV3eba8aITimc8+S0ATD/jUMLk6mymZIXPrU5N3+bIr7zeuRfKgg0FDMLScc5NLJuDgzaojvUKSWmDVSLOZvuCfMRMecpkpqMhgj2gArQ1QbNO4jkyAy7b9y4Z2GNRVYfjU1TVOwBkAIHORScUdx8bbSsWDurKxtcV0s2A5lpHwnZpPWufxWZsYSS6q15OxFnsddLncqLKeatWkiLpriepNt1pN0mbxFROznXKO3pZIOk6icmcvE+OZxotZp7c+56kxQgt9+jzIRmT9lJVvbTmCSlwtUQ35pMu/S86oejxeEZLBX/hNlUt8gdeC6yznLpAVS93HOP5qrpZRA4FrhGRO1X1+qIGi04o7j042gq9MafRowYcmPj8xJSZ7TgvV7PPtb6FfnMQe8pNNNrmbQ82ULLaoPneqGPd3AGzXAtN5iAI28MlGmAUce70E4jSsNRbOrXAoTP7arrgNHWufVbV02roY3P89yERuYzoESiFQnFRRZ/PePoFLNneYsn2VnsDiASnFIxe6ddf1QO9XI9TxNYVyzW7zqc2f2I8D9f0mlQ0WTplbik1kdlsTcExI8yZMa2bxNJpb5YL6fqpaHIq+hy2hw/ESMkx2jSCkIZRbyp+3whCmhrw4O4DeHD3ATzu//1AyY2riEZ+RZfXoBGR5SKyMjkGTicK0BSyqISix+MZPcN4HIGIvFRENgEnAV8Rkavi8vUikmxFuBb4tojcAvw38BVV/VpZ34vKfH7k51bTmO/zJ6hPjW8Q+Ym1Lj3L69vhnJMm2O8155rSDp9rDRp7T9aELbjSpTVmgitGuRlcyS7fs65rxm4yN4IwFYluZKLRU7EpHaBMz8wBcNJT7+3hgvPRaoGW3sdRvQy4zFK+GTgrPv4JcHzVvheVUGzMa8dBU4N+PpBUHMj4AO3lAxMudeM6TlUfpeV8lR+HdhpO/Lffz7Iw4mySEZB2f2Hn31RS5ba1y522gSHwEp9iIHT8hkadRCA22vWUhnSOZ2KH+57mTOV7UcYwTONBsiiE4ot//j0A6NHLuvdLtDFGMZI6qW1pX5X6uXNxrNijFptKnXGllzSckjnmBlMw/IWZPrpXspiCMKnT8SUm79PpNel6QEozbMTL/KakIzCnDAEZxOXLphY49rILAbjjpRcWXKg7/a5WGTW16LkicoaI3CUi94jI+ZbzrxaRH8av/xSR441zlZfheDye8SQKogx+md8g6VtTFJEG8HHgxUSPLL1BRK5Q1R8Z1X4KnKyqj4rImcClwHON85WW4eTRmFfCRr+95DCun2Ed86rax5Cizp3xHNYem38HRZHZ3DWXrNaYjk7n+RHJlJt+xLZGaCRyJxoh0GUuT0nHxzgVhG2tMZBOmz3NaaYaLuaVO+O82YMLdZjPJwL3xE5NRORzwNlAWyiq6n8a9b9L9Hzn2th11LJkoE6hoyk5iTmKtW8CMWxygyn5TZwFr9lXTya2rb0l3SY1ToGJnC2DlO8w9XhSS2AlaWIGUbL+xYbpa0yEoIQpgTklIVOxH9E0q2cbTZ63fmPJTanGpPsU6zCfDwd+ZrzfFJfl8Xrgq8b7ZBnO90Xk3LxGInKuiNwoIjdu3TqoJUoej6cfFCEMA6fXuFKHpmj7DbX+VojIqURC8QVGsdMynHjN4qUAGzZsSPXfXBJNoTHnlrlRK+MW5XVsU3UTiFrXOxfhEgTq6V5oWnO0RaZd/3nyNGCUuyMAACAASURBVMdsAIV0xNms1jlPSgPMHtuCK3krVRqGiZwNtEwFrXagZUpCpmOtcWljgdXTu92u25EJVxRr0RQ3AUca748ANmcricgzgb8HzlbVh5NycxkOUd7RiVUGP+Op5zP7WIvZx1pIaE+XH5bpWPc4A12FM0xyBUh+Exeh7boBRN8UjWOucLGVJ2/b5nRnf8TkOL0fYtp0Nne7ye6XKKRXqpj7JyYCMRKGrbZPMRGIyXFLhfv3reL+fat46pf+vP97tQgCLXUIxRuAY0Tk8SIyA7wSuMKsICJHAV8CXqOqdxvlPS3D8Xg8Y4w6vsaUvs1nVW2KyHnAVUAD+JSq3i4ib4jPfwJ4F3Aw8Dfx/mfNeHeMtcBlcdkU8FmXZTgmjz57Tdc6556p88fL6Mtp2/u65tBP+xquv5f12Lk4m7NxX9WGLuwrMqsLesyZd1eydoHZnKqXOU7M4vRWYGaeYXfQpSFhykTOM5nTxy0OmIo2h/iNY27Kv94KjLMW6EItyduqeiVwZabsE8bx7wK/a2nX0zIck+m9Ia2ZIOmwn66GQy8R0RJq2xmnrL6lfb/pOLW3zfMVpt6T/hwsfj9733lL+DLnrXW6k7UhKxy7BaTpRwwyArNdTnf0OYgFZOJTTEzmznEkFGeDJrNB9OjTHc0lRVfvhAJh6IWix+PxRCgT7wifeKE4vyKgMVe/hjjhn+t40kN+YqU6g8IlIt5+n9YWS83kbMQ5Pm8GTfLOJceJWR1phq32cUBnOV/WZE60w2lpMR23OWJpPQ+0mgSDrYiJF4qz20NacUrOUHZxGRcqXEdteygOsg+Xfob42TknbJvFOQIwaWLWy5rN2TaJ6Zwyk83Um1gQNozVKZEQTCdrT+eYzNPGcbJRxKb51TkXXBEvFD0ejydhvNNtXJhYoXjmk94MQPPZawc/2GR/xpWw/T87P1bAlT7yE9t1BvGZlPVblKuYVycnwNIuypRln8aXR1FwJSlviLaX9k0ZZvK0YT4vCRaYDRYAOGx6Oydf96bcMZ3xmuJo2PukQwBozQiNhbhQNW2XLDIm7gd4SInVvXwHs5Hp8nFKNoCw1CnyJ+al4XSZ3TkR59RjTels+tARlnGitrGHYiIwp4O0gFwWzAOwZmoHpxz6Y6DkISZFKKiPPo+G6R3RBzm9ojG4nXEmmT61sUHT98441jbkptpkl/Y599fLOeN8Xm6iWWbTCBN/Yq4gTB135y9OSxjvhhMLwoymOB1rirPBAg2iOlubB7B1fmXJhbkwBv9gfTC+q7I9Hs9kMoQVLSJysYjcGe/RepmIrMqpV7jXq42JFYo7H7eMnY9bhgakHpI2aEaybdiE/PD2vF9iQX8j02orrLSx7ZjdVS+TsG3zGZppONk1zokWaJrL0Ik4J+VBW0NUpiVsa4jZ1xJpskSa8fsm09JkVWMPrzj4e7zi4O/1cMMMhrPM7xrg6ar6TOBu4G3ZCsZer2cCxwGvEpHjyjqeWPM5WdqnDekpMWocTMixoYd7sVjvX+WlfRV+jW1mclEwJYttaZ+5YWx6b8Ww61zKfJbIfJ6RJkviQMvKYC8L2qcvakjJ26p6tfH2u8CvW6qV7vVqY2I1RY/HM55UeO7zIckeqfErdz/VEn6H9B6tCVX3egUmWFNsLo1+jQaxmsWJAf8Y9vWozToYliY4Tgn3/Sa5lyRkd87ZE7tta51N7dAkMN6ntEaSzSE0E3hptZO0zaDLEllgiUSa4oy0eNzUDocLLcE9+rwt3hjGiohcC6yznLpAVS+P61wANIHP2LqwlJUKjIkVikseiT7U5tIA9fpuOQPMLRwKgxx/BD9AvbimzdzEdlnWx2isbmkQtoVfQzrHgYQpATkTlwPctXBwfPRg9QnG1OXfV9XTCscROQf4ZeBFqlYfmtNer1m8OPF4PPXhGmTpP/p8BvBW4CWquienWulerzYmVlNsLovk+aJy+C+ma8kyyddWpvrknC/LUYzKi/s2AyVmmU1rbKRyFsM48BK/NzaHMM3nSIOMgi4rg32sD3YWzqecoaUMfAyYJXqECcB3VfUNIrIe+HtVPStvr9eyjidWKM6viITi9J6O13ZULCrBPAoc799A7nPdpnPBuTzBaFvaZ+6KU0bKZLYkeEMk/BqJv5GwnbA9LU3Dv9jiydPLHUctYAhfR1V9Uk75ZuAs433XXq9lTKxQ9Hg8Y0q9j5EeOhMrFJdui37dWkvEa2qe/QJJmdBuksdcC520aZjHcUAGYEEDvj83398kh5SnOEhqCbSULaWRiL+Kz/9QRJ7l2jaP1hJp76M40SyCS5hIhrkMqgaCnLnmlZsmsv280kAJCNvH0xKyvjHP+kZ/gjG5tWWvcaVvoei4lOZM4Jj4dS7wtxXaejyeSWI4y/wGRh2aYnspjarOA8lSGpOzgX/UiO8Cq0TkMMe2VsKGEDa86eyZPKos7XPBNI2zNIycxHb9HA1yWkKWBQ2WBfv3tlN1CEWXpTR5dZyX4YjIuclyoK1bt/Y9aY/HMxj2e/MZt6U0eXWcl+Go6qWqukFVN6xZs4agpQQtHeub6/HYqHu7/lADwpxlXS0NaGXOhTlf+wUN2BO22BO2rOedUKJlfi6vMaWO6LPLUpq8OjMObT0ezyQx4YpKHZqiy1KaK4DXxlHo5wGPqeoDjm2tNPYpjX0Tfvdh4v+BJpaRbtZYnTBnrnnlLQJaBV/vFkILISRoH7cQNrdm2Nya6Wuuk24+960p5i2lEZE3xOc/QZRRfhZwD7AHeF1RW5dx9x4SOYPHYUWLxzMMTLO7Yy4X5yuGsYfKNLFb5jHSFp4NlGfO9icQo4n238UoqSV527aUJhaGybECv+/a1uPxTDBeKI6GmV3RL+Q4WECJKTAOc5lIki9Ryf0byH3OC/f12qbgWhKDJrsGOjGBRaVdR1VydwEPVVL9J+1DhFZyrJIyrVsamcmQNq0XdKq92/Y+bXD3wu6cC3Nj3E1jFyZWKE7tiYTiotpPsZcv6KQwydemUvxNT4RPpo55yXlP31WVwrxFU2CmBCEd4ZeUt1TaZnFkLks78hxK59yCNtqCsKUBCxqJgZ3hEra2kqf5PZB/vWWMcWTZhYkVih6PZzzxmuKI2Lc6+qUb2eMIJo1+NbVRa3qDHN/RfO+7jdk8x5QuItEOJRVwkXZ58j6p2yJoa4TT2kodtwytcT4uXyILPGX64d4uyGTCv5ITKxSn90R3Ppzq7Wl+fdPnl6IM0R58Z3UKjmEJQZdxxmEuFeeZmMtJxNg0kbMmsxrlYfy/LCrRKz4XZP4Z2v5CMQShBoRxv63Yp5icW9AG09pqH3f8iNMEGrmilugCP2keGI/Q4+MIvE/R4/F4MnihOBrCqfgXODR+cStoEz5ibNCD1rtY719imlojv7b7VBaEMas6aI1FJFpfgKS1w1i7jAIr6bzETj1pa4cNDdvBlUCVIIw0xZ0s5eDGLqe5FOG41ePYMrFCccX/RM+q2bN+KeEQN/UQVbSXR7H1w6j9eY7kmfw9uQIYseAtNKXTgtD0D6YEn1nPOI6ee9wtIMPYZDbrJMLZPBeKkPzLh6aAxIw+hyxoQBALwkCVIJ5oQ0P2tYViSCMWA9tbU3z9sWPjnu92vVMjQUQuBn4FmAfuBV6nqtst9TYCO4EW0Cx6pGrCYklm8Xg848Jw9lO8Bni6qj6TSIK/raDuqap6gotAhAnWFFtLo6kvLBUaC3GhX+7XoUDT6VVzq5NSLbCCWdppY/QPdB79pPF70hVc+yvay6kkMGNWs0Wc87XGTrAlaR8YprBpFie73jRDJQiiPiItUVO5iaamGMRupzmU6UYUgFk3tZ01M30+zW9IgRZVvdp4+13g1+vqe2KF4ux9jwCwd826tDDsJddhQpg4P96AI/TJGMlHXeW72L6XprgsmmeekE5laHeb1e25Gb5DNVauiGhbKKpxbNZNzqVXvhhCUUwB2fEpNlUJYqEYaEgzjAUkjbbwbGjInjBa77y1eQBf3/LkgpvgiPsHcYiI3Gi8v1RVL+1hxN8BPl8wm6sluqGfdOl/YoXiV++5GIAXvuTi9rNaBiYsJsSnVwc2wdulWQ4w59FF8Isa37u6l/wV9Zsn5M0VLdk6qXPdK1+ygZauJX95S/2MnMWOTzGgGWuAQaC0VGjGDveG0U+gmtpHIslZ3LawkutPuzie/oes4zrhLhS3FZm0InItsM5y6gJVvTyucwHQBD6T083zVXWziBxK9IzoO1X1+qJJTaxQ9Hg844dQX/RZVU8rHEvkHOCXgRfFm87Y+tgc/31IRC4jegTK4haKcwcG1XyKi0Xrq2CaFvoQ67gfdZnJZXMZhjmeDKGC5pnVahQ6mMtJEzHqpSPW3W2SRG4xNMf2NEIjPhrQ9hVG65ujWk0NIm0wrhqEjU5Y1SgnhEYjkmJHzDxSfFNcGJJPUUTOAN4KnKyqe3LqLAcCVd0ZH58OvLus74kXijO7Qlqz9QfRxyEYsejINT9zyrN1RvV5FKbnZM5lTOn2pWXSddq+wow7PJuCY9s1J1RJBV1aGSEZnYCg0Qm0NE2TOSB1nOQv3rn3sLI74cZw4p0fA2aJTGKA76rqG0RkPfD3qnoWsBa4LD4/BXxWVb9W1vHEC0WPxzNmDCf6/KSc8s1EG1qjqj8Bjq/a98QLxYVlAUFzglJxYs2izkiyxOpEaVJ5VfOzK2jQ3d5Vo+5H867UNpOWExVlE3TSdVJnyiLQtnpm1KcT1s5f+ZJJ5O6Y3N3aZHYtNECQiUQn0Wdr0CXHZE4dx8wGzYKLd2fS1z5PfPL2Qd/fSmNeaczX8EnU+WGaX86yftWobx73Om6v7Wu4/qqbjBbWVXGThvE112Jda+dlCqLcetliw1+oKql6qY85U69jTkfHSmQaq/FqhQGtMGjnKSbnk+OWBjTD+KUBzbCRer8QNlgIG6njuXCKbfMr2Da/gn+/9+l13MHUPRxg8vbA6EsoishqEblGRH4c/z3IUudIEfm6iNwhIreLyB8Z5y4UkftF5Ob4dVY/8/F4PCNGo+izy2tc6VdTPB+4TlWPAa6L32dpAn+qqscCzwN+X0SOM85fEi/BOSF+XkslvnbnRcwd2GDuwAYaSGSLZMzIYanzdY+zGLZ2B+xaQYm2UHbd7XszjPtTNE5GG+wqT95qUhZrhcZxojV2HdPRHEMlpTUmQzZj7THRINXQHjvaYqOjPcZaYnLcEOWw2cc4bPYx7nzZn9V7vyZUU+zXp3g2cEp8/GngG0Rh8jbxo0wfiI93isgdwOHAj/ocu81U8qhTM5Q3rBUtw4qK9jJOgQ+xzKeZ9eNZE7izfddxL4r8nllnYKV+OytXjDUsqa7U9Ak69tXtW0yn6rR9hIZ9b94mNdrYFmYlCpWoECbzV2mn5TSCsLO3Yhhk/IYhzaTA0MzM1B2AhXAldTLpP+T9aoprY6GXCL9DiyqLyNHAzwHfM4rPE5EfisinbOa30fZcEblRRG7cunVrn9P2eDwDY7FrikVLbaoMJCIrgC8Cf6yqO+LivwXeQ3SL3gN8mGgdYxfxmsVLATZs2JC6pdZtxPKeFJSdV1YDcrmWUWwflhq/lzlXbzMw8jTKAk3TOcqddKXpsl6i1532cSJ3rvbaHZVWjdqZGmU7YTvSGzvzNSLRHa0xuhlJbmIrDGgEYbueqUGm5pokbxMlbyc7cZu5iYGE7Yj1fNjgxw+vybsb1RlzgedCqVAsWmojIltE5DBVfUBEDgMeyqk3TSQQP6OqXzL63mLU+Tvgy1Umn6U1Ix0Hbt075hSZdaNk0Carw5iu6UV59freMacfs7oK5n2yrGNOzyU9Z1VDFJruRjWkJZLZCKKT5I2R8N0KoRHbeK0wgFhYCnRWu8Sm81QsCZuGaR0gTBHtjLNsaoH5hfoy8wRvPl8BnBMfnwNcnq0gUTr5/wXuUNWPZM6ZKfQvBW7rcz4ej2fEJEGwste40u9PxEXAF0Tk9cD/AC8HyCy1eT7wGuBWEbk5bvf2ONL8QRE5geiHcyPwe71M4pr/fCcAJ/3Gh9wCLMMKjgwZZ7O+6PprCqA4m+slJjPk9GPEOdxN43Q6t9FNpy/rmfw5pqtltcPEhI4rdpnWpG6UooaB0x18Sa3Si08EdHIcW+Z0wiAyvdvaoRrmc2ffxR0LS7jr196Vc3E9MsYCz4W+hKKqPgy8yFJuLrX5Njn/Uqr6mn7Gz9KalVofeVrvqhOjL0MQpIRHKixJf4J73KLi2XplJrulX3OxSNmYQh9+RVt/5gYR5vhZ0z4jJJM2qbXPZk/GxNJzjn2KhmmdSu0zzORUJDo78fhcaAhIVWVfPObqWeteCv2xPwvFcWP1TY+w68mrAGhN9/gN6NN32CX8+ujL7LOv/noMYGR/FJx+JPq95rz2ZX5Fs20f46c2n3W18bLaYlKY+hHofIBq+A4xNnrI+hdTPkXjQzQ3nMXwKRIGnX7jsqTeVNDZZDYQZS72I97108PgJLfLdGLMTWMXFpVQ9Hg8Y8CEC8WJX/ts8rXb3se+VQ32rWqgjWyqwgB/wUaQhtDL9ZS2qdKf5Zpd55NXr2p5Mg+ne2HWM1ahtFfF2D7DzLn2Wuhs/fYKFVv76Jy5osXsO7XShXT91AoWY7VLGAbt4UMVwjCIXkabZJVLUs9c+TIlIeuW72Dd8h3c97q3UjeTvsxv0WmKSx6NUg1aM2bOSLEt1Yvv0HlnmmFgWFy17Zhj8wH2EGyJ5lIyZsVyp8cRKKndZ4YWWEs7C40ThjmunYCO6V9M+x2jwEsYRu+DlPqihIkpnFnRkpjiIsQRmSB+3xln98IMj80t6ftS8/Dms8fj8SSMwGqqm0UnFFfcHi0B3PmMQwkHcXXD1Dr6paZAT5Za10Hn1bOV5yVMm/Nqv8kfUgxNLX1k1Mk9Q2eFS3YsLdAO6cRQzJUuaq5vyQRdEm0RIAylS1uMyo3gCml/WCvsWAwi2u57oRXw6E9zV9T2jxeK48XX7v4gACe+9sME8bNbrN+hCgKjziVytkhybsQ6U6evORQIqaqbQ9RB1dUtuXNwiDhn03N6+pHI+k+TSHA2ylzYrtPGzF9sC1/tlGt8wWatPFM66Tc0hhMVAkMQNgLaN/eAJXP84Lw3FV9vjwg537e6xxF5D9GGNCHRSrrfTh5Slal3BvCXQIMod/qisr4XVaDF4/GMHgnV6dUnF6vqM1X1BKLlwV0Z6CLSAD4OnAkcB7wqs22hlUUrFJdvXoh/tkYzfurXsu5fTtNy66Hv0gh03vm8CG2m72z/XeMV9W8bz1LfKZLu4t8qiESnotWFfRjjmJFotb3HiC5n+o7LknKlE00224ShtF/d+zHad+5uhZ1x/ufB1SU3pQ+0wqufYTqbygAsz+nxROAeVf2Jqs4DnyPSLgtZtELx69eeT9BUgmaPv0q9fGij9qWkhKXjZHoVkC59uba31LMJvFIhpfnCUsxxCurZ5pVN40n3Y0nRwVbf9j5J1SEWiGIc0xGM0BaM2XSdSDjGKTmGsAzjTWdb8asZBswtTDG3MMVPX/32kgvvj2GtfRaR94nIz4BXY9EUifZt/ZnxflNcVsiiFYoej2dEuGuKhyR7pMavc81uRORaEbnN8jobQFUvUNUjgc8A51lmUuaFtrLoAi0myx6YA2D34UsG9sS5Ue+t6ITSd0Ape66udeFOkeyk3DJWYcTZ/Pcvuq68SHQ/ARkBazQ6GTTVpjvoksS3s4GXdtd5EbnURNJbkW1/aEUPF1SdClrgNlXdkHeyaNvCDJ8FvgL8WaZ8E3Ck8f4IoCsYk2VRa4rynZuR79zM3Epp+xdriYzVaSanTN6cMbRzvq/5VzGHy8zqkn6d/Iq2Mte5mKtDLP3mmdBZM7qQjLlt9w+SXuniQsaMTvkXM6Z0l4+RMlPaNKk79+egZXtZedc0K++adpxkH7hrij0jIscYb18C3GmpdgNwjIg8XkRmgFcSbXdYyKLWFD0ez5BRhrWE7yIReQpRSs59wBuA1LaFqtoUkfOAq4hScj6lqreXdbyoheI14b8AcOqLL2JuVXSp2TXRVUxLqGYu5uUf1pL3l+nPNk4vc8sdyzWHMd+Sc5tLUXmiMWfqlm4pVnJ9qZ1x4srZMdL9Zc3ipHtLUnfXGsf2icz7xHhOm9JJYre2x42qty/J3D1IzBsd7bAzNxdphj97dBn3XfwnORdUH8PKU1TVX8spb29bGL+/Eqj0lNBFbT4nLLl7C+GUEE5JNWHUi5pfg2lQF6LaiUI7zqnIRK9qvmfrW9vnmdXOpnX/kejS9kphf7mmdMoFYo8+l5nSGkr71TGtO+0Tc9qWtgMgoogoS346m3OhA6CTI1T8GlP2C6H41fsuobGgNBZ0IHmLKeEzCVT1fbmeK+rXYUznXXJyBFQqz9A2fsm5lN+woH2Xf9E6Z0uqTo4gtAlLjQWi2aYj/LD6GkMVNAzQMGhfQ7PZoNlscNeFg9cSO9c+nJScQbGozWePxzNkxshS6pW+hKKIrAY+DxxN9IyVV6jqo5Z6G4GdRI+RaCZheNf2dfCdf43Wep742g/TmO+Uu6697aLt1OltPkX+Rlu56ZdL+fEcysvm4dSm4Jz1npXUd07DMctz6+WlunSqZ8fL21bM1b9oq5eukJlzqiwnVccsU4l9h8kcOm1SvkalfTEitNc6qwpTgfKTV1Z6EnEtjPNeiS70az6fD1ynqscA18Xv8zhVVU/I5CVVaV8LK382R9DUXCupjL6X1Y3gVzTlV3Qcv+g6qy4T7MlUKuvDNGUz9SRvDr2Y1hjj5NUrSNWxjpPbpruP9MoXw5ROfI2WVTCoEC4E7Htwec4FDZZJ32S2X6F4NvDp+PjTwK8Oub3H4xknFPb3QMtaVX0AIP57aE49Ba4Wke9nlvK4tkdEzk2WA23durXnCV/39bezsEwIp+jst1imQdWg6YwdxjU7ba6QadN1jnxtLve9BecIta1NnmaWNy+1WPimBpoxJ7LRaJumatMWuxK7U8d5Uen8yHRRACaJNk/NNtn4fwazPVgZiz7QIiLXAussp6o4K56vqptF5FDgGhG5U1Wvr9AeVb0UuBRgw4YNfd3SZQ8uMHdQfOliT9NxziVUDJ9en0v+Un0Z4yflyVUX+B57yVnsmoNDm1Lfq3EtSf3SZYKZNtm5VPFd2vyWKZeeWe7qX7TNzWiSegJgVnCbOYxd/mKbAzJb0RyJtK/RmGcreRzFnhlGxhgLPBdKhWLR+kMR2SIih6nqAyJyGNFmj7Y+Nsd/HxKRy4i29LkecGrv8XgmA2G8tUAX+jWfrwDOiY/PAS7PVhCR5SKyMjkGTgduc20/CL5x9fnMPtZi9rFW/ytLeqHI5BwQqTzKPFOuq40xvzLzuor7weGanZO8s21sETRrPUu52W/XeB1ztitwlmdyF5nStnEKAi0o0YK2MF3e5abbPg3bp7nv3DczEtRtg9kaNpkdGP3mKV4EfEFEXg/8D/ByILX+EFgLXCaRjTIFfFZVv1bUfhgsvXUTAHtPPbpyJKySaQ3dZu2Y0kt6TqFZ7GD+FqXndC3fM63PspQeNdJW8q1Qp/L2eMbNyTXFzRqi1h8AwXwWS2Z864UVTE4FgtiUbigr7210DzhsxlfeOdGXUFTVh4EXWcrb6w9V9SfA8VXaezyeyWXSzef9dkXLV+//awBOPuuD7Fsd/bpKi2LtiILzyel+nwddFDSxaHFVErl7nVu7r7J74Khdtt9X6Ss3CELqvtg0uJ6CLnTXtylxKSU0q82WBF46XVmiNu3J2SSMEc0Wbc9ZHpnh1kuGt5zPigJjbBq7sF+sfS7im1e+hdlHW8w+2iJopT9M59SBMp9aSZuRpCiYc+5lLiXXXJSwnptw7dKXpT8X/2KuHzFbnvj/sv2qMU9TOqbK6Sov8jF27nmBrzFpH0r3OYBmADumYMcUP/3DP2UssF1nzrWPI/utpmhy/VfeAsBzfvsjSCwYg1Zxm0o+wjwNriqmdlWkUdrK+xyvaN7OaUC2c0X3Jlvf0t5pyaBKWrKaX0ib7zI7rEVL1dSkO/WzWmO3FzI7AXO8zNZjtutI2k/F9RrKxnPfau1vVHjz2ePxeAzGObLsgheKBsseajK/MvYoiKA250KR38ygl0RuJ/9gr6T6M+ZmXo/jtbm2yfWJ2s5RPHalKHX8Jus3rJQwbvQlkF6VlhozPYGUUpsXMc/1FRb4Gs1rCwXmozdT28cg2mwy5qaxC/u9T9Hkm1e+hWUPzrPswXk0SLuDnHD8Zxj1BhEu2OZYyfdXdF0lfeT5EPN8d1118/yLqfdGvcSPmO3X9PtZylJ9Z/yGXRtIaGf+qTYFPjernzEW4rNbppjdMsW9b3kj44QQ/ei6vMYVryl6PJ56GcIOOCLyHqINZUKilXC/naycy9TbiGXbwiK8UMxw3TejB4U/99Uf7nbwW+glkbsXchOr88xsW3nK5IzeqEh6bq7zrNCmKKWnl4TxrvFtfWXqZIMjVVJ6ykzsIlPauqy5oE2W1CNOjTnc/c4Rp94UMCQt8GJVfSeAiPwh8C7ih1dZOFVVt7l27M3nHL73mT+lNS20piXa/800r5zNZEczwTTTsuOM0MrIu+YyE9o2f6tJXGCWV07RqWRiZyRvrsltlMcySbJtSkzptMnc/Uqb0tI1ZrK0T1oCDYWGjm4JnwsF7oCuVz/DqO4w3i7vv8cOXlP0eDw1Umld8yEicqPx/tJ4NywnROR9wGuBx4BTcycUbVuowCdd+vdCsYA137wfgF3PWMfcAZFS3aWF0NFcqmzXVWcOoTUSaym3m9aDiUQXuhWKxsm7lrxzuPXVZRonSd2og4lL6j6bHoeu/4V2caZThyh5VztDNW3sCjjg1tiwez3jjbv5vK3Ix1e0baGqXq6qFwAXiMjbgPOAP7PUrbxtoReKBXz1Jx8G4JdOeCe7nncQADM7M/4fX37HhAAADttJREFUB8GRK3gK23QLsVFtKmH3SRpzKbsXjvOvxb+oxmFOX+lxyjeOaOdnm/LK6E+LxrdI3PyEHNM8F7ShsKIJwPK7Z/nBJ8bXj9hGqe1RA0XbFmb4LPAVLEKxYNvCXLxP0ePx1MsQHkcgIscYb18C3GmpU7RtYS5eU3Tgqpvfw+nPuRCAPUcuZ2Fp4BadrarR9dImp32pppnS+hJNKT8S3UuUvSh528n8T4oLxi7SLq3X6zqOpb4a71NzTppk3AJdZrElyTs1VGaejd0BSzdGD7H/4UcnQEtMqC3kUchFIvIUolDUfcSR5wrbFubihaIjV99wIQCnnPEB5lY22s5kMxJaSXC4+N5c2vQrSHvBZS4F83JuUySUSwRp4T3OnrOl6hQJ0q4L6vzpUoBs45i9ZPuNbTdtKNMPB9x28QQJwxgJB5+oqKq/llPutG1hEV4oejye+lCGkrw9SLxPsSLf+NpbWblpDm0I2ojNzeRVQN85iz3glGeYKldruXU+DtfcNV6Fa+tlyWDe8j2Xc9algHTXzT3XUfzsuYyavX9CdpmfhKCBooGy9IEGd757ArVE3Jb4jfMyPy8Ue+A/rnsbK+6fZ8X980zt02iddGAxgx0FRy/CL1dAVRBWdeEyl0IhVyasbOPktbeN0SX87PWyydfW6yr6gTAFXMG1pRK6DRpzwpItDZZsaXDHeydPILaZ8Oc+92U+i8hq4PPA0cBG4BWq+mimzlPiOglPAN6lqh8VkQuB/w0kD3J+u6pe2c+chsXXrzkfgNOf+24eevYKIErXKU870e6ARglFy+RcqJJSUxh0KfEj9uLrrORfLPPDGsKo637Z+jPPxce25Xdd98+cP/nzt6btWGitiOzNlRunuOWvJlgYJoyxwHOhX03xfOA6VT0GuC5+n0JV71LVE1T1BODZwB7gMqPKJcn5SRGIHo8nB6Xz1MGy15jSr1A8G/h0fPxp4FdL6r8IuFdV7+tz3LHh6u+9i7X/tZ21/7WdmV1huSlXgLMfzaXNOP1YF5mcBeUJtfkkDfO161zm2GZK57Wxrm/O9pedYFxvelfA6psarL6psTi0RKLos8trXOk3+rxWVR8AiB9of2hJ/VcC/5wpO09EXgvcCPxp1vyeBK76wbsBOPNJb+ahUw4DoDGvboKpXzPTkaomr3UVjpP5bS/vyczO9Ou0AW/RXCw5hNZ5ZtvHta0rX3LbdGPeznA2qnHgrcp//39j8myVWhhvf6ELpZqiiFwrIrdZXmdXGUhEZogyz//FKP5b4InACcADwIcL2p8rIjeKyI1bt27Nq+bxeEaJsvgDLUXrD0Vki4gcFmuJhxFt9pjHmcBNqrrF6Lt9LCJ/B3y5YB6XApcCbNiwYSzv6Ffvubh9/JxzPtIOWmS1maJHFfS9uUSOBuNMUdClkqZZf9DFSTulu701UKXGYUmgp9M+rTHS6ab7WjLzFAVtdBr99I8Wk3aYYXwtYyf69SleAZwTH58DXF5Q91VkTOdYkCa8FId1iZPCDZ9+Iys2z7Ni8zwA4ZQQTpl5cEauls0PaaHXPMFKPr1BkjemQ3mXv4/u8rLxrG3KfILW9gX5jJl+tRELQ4WlDwpLH5Sxe4RA3ezveYoXAS8WkR8DL47fIyLrRaQdSRaRZfH5L2Xaf1BEbhWRHxLth7Y4PM0ez/7MYjefi1DVh4kiytny9vrD+P0e4GBLvdf0M/648x/XvQ2AX/hfH2TP2uhWS4voKYEOJnFt5qcjpUGXsuBCWXm/5neJyVyUv1gaBMr2m+2jq31nonnrpaUlMB0VzD4q3PqR/eA3XxVak20/+7XPQ+D6r7ylffz8X/sQcwcE9j3n6hZ+2S9+H/6+3HlWEXB1CVXbfGxtcubskjDeVccg7bo0TGm0/VhcDWB6N9z9zsVtKlsZYy3QBS8UPR5PvUy4UPRrn4fMd774JpY+3Go74FOO54x2k7ehQxc9BE2qBjrGyjHuEHQpWyPtul7bev+zmmM7SCMkO0LMPir86P37gbmcRYFQ3V5jitcUR8C3rngzZz4hSsnY/bS17DsoytUIGxC0Chq6mIl57XpN0THau25G255XprzOR7RWNr8t11M0ZuH4ZldTiixEFWd2wKp7ow/w25eN8RP3BoqCep+ix+PxRCg+0OLpjeShWAAveFmU9L13dYNgIVJHwinpKwCSG4ktqufaPqevsn4rRX8HoWlm510yZsocJ92HNOODGZjeGR3e9qH90Fy2MU6ulh7wQnEM+PaXOqbWL/yvDwKwa/1Ux5TWfCFVJLyco7cmDtHfohU5LuP3KuBqE6q26yUj+Gzzit+HU8qybdGbW96zH0aXyxiiUBSRNwEXA2tUdZvl/BnAXwINome3XFTWpxeKY4aZvvPc34q0yda0IKESTkdfxPbKCygWFllcBIkFm0/Pxb9YRE8Crk6hakONrgwhmOSWBtHiJO76My8I89GhCUUROZJoUcj/5JxvAB+P62wCbhCRK1T1R0X9+uizx+OpDwXC0O3VP5cAb8EaAgOiZzzfo6o/UdV54HNE2x0W4jXFMeZ7/xRFqM849m3seupqdh4RfVzTu9KrKMp8f67+wd5W0Vi2GMscW6PHNiq071fTbA+Z42sUYGFZ9Gb5JuHAjQtc/+9vweOAu6Z4iIjcaLy/NN74pRQReQlwv6reIvmunMOBnxnvNwHPLevbC8UJ4Gt3/AUAZxx8LgB7T3wiu9dNA7RXUADRg4/KTGDIFT5duAQ6Ctq4tC8NCFVpXyBUze7MuZvtNaC90mjpFuXIa3YB0UbCHle0SvR5m6puyDspItcC6yynLgDeTvRw+yJKvMh2vFD0eDz1oaA15SnmbVsoIs8AHg8kWuIRwE0icqKqPmhU3QQcabw/AthcNq4XihPE1x7utixe+JKL2XV4lPw9vVvb25OFjShA0KYgJcU16JIf2TWCLjZKoselgaJe2+eY6yodjVBa0FoWHS/bovzgb30QpW8GvFpFVW8F2rv8i8hGYIMl+nwDcIyIPB64n2jn/98s698LxQnnW1d00nlOf86FbH/qSgAWVgTMbo+++a1pqbSbTFedMnPZlqpT0KanjWmrzs3SvrEQHc6tgtnt0fGqe5p86/L9dfXJgBhhnqKIrCdKvTlLVZsich5wFVFKzqdU9fayPrxQ9Hg89aFaV2S5wpB6tHGc3bbwSqDSU0K9UFxEXH3Dhan3p73gvQBse+YyGvMQxp92MA8aH2tgmNkZjSvPlO7SNF3zFx0DKu0xs+avxRwWNTznQadO0ILWTFwnhHAaVt0TXeh3vvimnIl4asGvaPGMK9d++x2p96ec8QEANBB2HBV99OG0MLMv+iduzXaEXTZBPGrYeV+e8F3haYAJZVFpcz5G1404qXp+ZaffZZvDdvvF9bS8cUfRVtGuJuOPF4oej6c+lLHeFswFLxT3I77xtbday097fmxmH7+MqX1R2cJyYWaHxsegU0Jjb2waT9EVuIHo+6CBxGXaHZU2k6YDo23SvpHWVIN404XmssgEntodz+cAYWpPdG5+Jay+M4qgXP9ln1w9FuzPW4eJyMuBC4FjgRNV9cacetZF2SKyGvg8cDSwEXiFqj7az5w81bn2O+8oPH/a89/LnsOX0IjN7MceP8XU3kRgCrPbO+b3wgFRmyXbOj49YmGX+C4bc8q+QyLpN70DGvNR+7lVwnQs+JpLhZUPRA20ISzdusDXrzm/luv1DA4FdMI1xX7XPt8GvAy4Pq+CsSj7TOA44FUiclx8+nzgOlU9Brgufu/xeCYV1UhTdHmNKf0+ze8OgIK1h2Asyo7rJouyfxT/PSWu92ngG4DdxvOMjDJN0uMx8YGWcooWZa9V1QcAVPUBETk02zhBRM4Fzo3fzonIbYOY7Ig5BOjaE26RsFivbbFe11N6abSTR6+6Vv/1EMfqY3nfSoVi0aJsVb3cYYyeFmV3NYh2z7g0ntONRQvJJ5XFel2weK9tMV9XL+1U9Yy65zJsSoVi3qLsChQtyt4iIofFWuJhwEN9juXxeDx9MYxNZtuLskVkhmhR9hXxuSuAc+LjcwAXzdPj8XgGRl9CUUReKiKbgJOAr4jIVXH5ehG5EkBVm0CyKPsO4AvGouyLgBeLyI+JtgwvfX5CjNNGlBPIYr0uWLzX5q9rkSE64esUPR6Pp078M1o8Ho/HwAtFj8fjMZgIoSgiLxeR20UkFJGiZzqcISJ3icg9IjL2q2NEZLWIXCMiP47/HpRTb6OI3CoiN/eaKjEMyu6/RPxVfP6HIvKsUcyzFxyu7RQReSz+jG4WkYl4sIuIfEpEHsrL+53kz6xnVHXsX0Rrq59CtOJlQ06dBnAv8ARgBrgFOG7Ucy+5rg8C58fH5wMfyKm3EThk1PMtuZbS+0+0+edXiXJXnwd8b9TzrvHaTgG+POq59nBtvwA8C7gt5/xEfmb9vCZCU1TVO1T1rpJqPT3jdcScTbS8kfjvr45wLv3icv/PBv5RI74LrIrzU8edSfzfckJVrwceKagyqZ9Zz0yEUHTEtpzw8BHNxZXUMkeMh/FkUOBqEfl+vNxxHHG5/5P4GYH7vE8SkVtE5Ksi8rThTG3gTOpn1jNjs5/iuCwnrJuSZ9e68nxV3RyvDb9GRO6Mf+HHCZf7P5afkQMu874JeJyq7hKRs4B/A44Z+MwGz6R+Zj0zNkJRB7uccGQUXZeIOC1z1OhhPKjqQyJyGZE5N25C0eX+j+Vn5EDpvFV1h3F8pYj8jYgcot2P3Zw0JvUz65nFZD4XLSccV0qXOYrIchFZmRwDpxPtYzluuNz/K4DXxhHN5wGPJe6DMaf02kRkncR76InIiUTfrYeHPtP6mdTPrHdGHelxeQEvJfrFmgO2AFfF5euBK416ZwF3E0UKLxj1vB2u62CizXV/HP9dnb0uoojnLfHr9nG+Ltv9B94AvCE+FqINh+8FbiUnk2AcXw7Xdl78+dwCfBf4+VHP2fG6/hl4AFiIv2OvXyyfWa8vv8zP4/F4DBaT+ezxeDx944Wix+PxGHih6PF4PAZeKHo8Ho+BF4oej8dj4IWix+PxGHih6PF4PAb/P5i36Ch6YDR8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The next command ensures that plots are shown within the IPython notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Filter out solution values that are associated with points outside the unit circle.\n",
    "u_evaluated = u_evaluated.reshape((n_grid_points,n_grid_points))\n",
    "radius = np.sqrt(plot_grid[0]**2 + plot_grid[1]**2)\n",
    "u_evaluated[radius>1] = np.nan\n",
    "\n",
    "# Plot the image\n",
    "import matplotlib\n",
    "matplotlib.rcParams['figure.figsize'] = (5.0, 4.0)\n",
    "\n",
    "from matplotlib import pylab as plt\n",
    "\n",
    "plt.imshow(np.log(np.abs(u_evaluated.T)), extent=(-1,1,-1,1))\n",
    "plt.title('Computed solution')\n",
    "plt.colorbar()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}