FEMfromScratch_2DLaplacePb.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "2DLaplacian_FEM.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i9jVyxtEkpHh"
      },
      "source": [
        "# FEM from scratch: Laplace problem"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2e1y-9bVR-nR"
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FqCaUpqFim8z"
      },
      "source": [
        "We implement the $\\mathbb{P}^1$ FEM for solving the Laplace problem with Dirichlet and Neumann boundary conditions. On the unit square $\\Omega = {[0,1]}^2$ we consider \n",
        "\\begin{equation*}\n",
        "  \\begin{cases}\n",
        "    -\\Delta u = f, \\qquad & \\text{in }\\Omega, \\\\\n",
        "    u = g_D, \\qquad & \\text{on }\\Gamma_D, \\\\\n",
        "    \\boldsymbol{\\nabla}{u}\\cdot\\boldsymbol{n} = g_N & \\text{on }\\Gamma_N.\n",
        "  \\end{cases}\n",
        "\\end{equation*}\n",
        "\n",
        "We test the solver using the simple analytical solution given by \n",
        "$$\n",
        "u(x,y) = (1-x^2)(1-y^2),\n",
        "$$\n",
        "from which we infer $f = 4 - 2(x^2+y^2)$. .\n",
        "We consider three different configurations of boundary conditions:\n",
        "\n",
        "1.   **Problem 1**: full Dirichlet conditions, i.e.  $\\Gamma_D = \\partial\\Omega$\n",
        "2.   **Problem 2**: mixed homogeneous conditions, i.e. $\\Gamma_D := \\{ (x,y)\\in\\partial\\Omega\\ | x=1 \\} \\cup \\{ (x,y)\\in\\partial\\Omega\\ | y=1 \\}$ and $\\Gamma_N :=\\partial\\Omega\\setminus\\Gamma_D$ (which, indeed, yields $g_D=0$ and $g_N=0$). \n",
        "3.   **Problem 3**: mixed non-homogeneous conditions, i.e.\n",
        "$\\Gamma_D := \\{ (x,y)\\in\\partial\\Omega\\ | x=0 \\} \\cup \\{ (x,y)\\in\\partial\\Omega\\ | x=1 \\}$ and $\\Gamma_N :=\\partial\\Omega\\setminus\\Gamma_D$. \n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KsYDSdp5K0H9"
      },
      "source": [
        "### Step 1: Mesh generation and reordering of boundary nodes"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jx8ChM2nA9YI"
      },
      "source": [
        "def left_mesh(a, b, nx, ny):\n",
        "  '''\n",
        "  construct a left-type uniform triangular mesh of [0,a] x [0,b] starting from a \n",
        "  cartesian grid with nx cells in the x-direction and ny cells in the y-direction\n",
        "  '''\n",
        "  P = np.zeros(((nx+1)*(ny+1),2))\n",
        "  T = np.zeros((2*nx*ny,3), dtype=np.int)\n",
        "\n",
        "  x = np.linspace(0, a, nx+1)\n",
        "  y = np.linspace(0, b, ny+1)\n",
        "\n",
        "  for i in range(0, nx+1):\n",
        "      for j in range(0, ny+1):\n",
        "          # create node matrix P (ordered from bottom-left to top-right)\n",
        "          P[i*(ny+1)+j,0] = x[i]\n",
        "          P[i*(ny+1)+j,1] = y[j]\n",
        "\n",
        "          # create correspondence matrix T with local backward orientation\n",
        "          if (i<nx) & (j<ny):\n",
        "              T[2*(i*ny+j),0] = i*(ny+1) + j\n",
        "              T[2*(i*ny+j),1] = (i+1)*(ny+1) + j\n",
        "              T[2*(i*ny+j),2] = i*(ny+1) + j+1\n",
        "              \n",
        "              T[2*(i*ny+j)+1,0] = i*(ny+1) + j+1\n",
        "              T[2*(i*ny+j)+1,1] = (i+1)*(ny+1) + j\n",
        "              T[2*(i*ny+j)+1,2] = (i+1)*(ny+1) + j+1\n",
        "\n",
        "  return P, T"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CalT2Et5-96j"
      },
      "source": [
        "def Dirichlet(a, b, Point, tol, pb = 0):\n",
        "  '''\n",
        "  condition defining the Dirichlet boundary\n",
        "  '''\n",
        "  edge_0 = Point[0] < tol\n",
        "  edge_1 = Point[1] < tol \n",
        "  edge_2 = Point[0] > a-tol\n",
        "  edge_3 = Point[1] > b-tol\n",
        "  \n",
        "  if pb == 0:    # full Dirichlet problem\n",
        "    diric_bc = (edge_0 or edge_1) or (edge_2 or edge_3)\n",
        "  elif pb == 1:  # mixed homogeneous problem\n",
        "    diric_bc = edge_2 or edge_3\n",
        "  elif pb == 2:  # mixed non-homogeneous problem\n",
        "    diric_bc = edge_0 or edge_2\n",
        "  else:\n",
        "    print('wrong pb argument:') \n",
        "    print('  0= full Dirichelt, 1= mixed homogeneous, 2= mixed non-homogeneous')\n",
        "\n",
        "  return diric_bc"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "r--FAF3T_qdt"
      },
      "source": [
        "def Neumann(a, b, Point, tol, pb = 0):\n",
        "  '''\n",
        "  condition defining the Neumann boundary\n",
        "  '''\n",
        "  edge_0 = Point[0] < tol\n",
        "  edge_1 = Point[1] < tol \n",
        "  edge_2 = Point[0] > a-tol\n",
        "  edge_3 = Point[1] > b-tol\n",
        "\n",
        "  if pb == 0:    # full Dirichlet problem\n",
        "    neuma_bc = False\n",
        "  elif pb == 1:  # mixed homogeneous problem\n",
        "    neuma_bc = edge_0 or edge_1\n",
        "  elif pb == 2:  # mixed non-homogeneous problem\n",
        "    neuma_bc = edge_1 or edge_3\n",
        "  else:\n",
        "    print('wrong pb argument:') \n",
        "    print('  0= full Dirichelt, 1= mixed homogeneous, 2= mixed non-homogeneous')\n",
        "\n",
        "  return neuma_bc"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "b_4dk1nCL1bT"
      },
      "source": [
        "def ordered_mesh(a, b, nx, ny, pb = 0):\n",
        "  '''\n",
        "  construct a left-type uniform triangular mesh of [0,a] x [0,b] starting from a \n",
        "  cartesian grid with nx cells in the x-direction and ny cells in the y-direction\n",
        "  P,T are modified in order to put internal nodes first and Dirichlet nodes last\n",
        "  '''\n",
        "  P, T = left_mesh(a, b, nx, ny)\n",
        "  P_new = np.zeros(P.shape)\n",
        "  T_new = np.zeros(T.shape, dtype = int)\n",
        "  \n",
        "  n_bndry = 2 * (nx + ny) \n",
        "  n_inter = 0\n",
        "  n_diric = 0\n",
        "  n_neumn = 0\n",
        "\n",
        "  tol = min(a/(2*nx), b/(2*ny))\n",
        "  v = np.zeros(P[:,0].size, dtype = int)     # permutation map\n",
        "\n",
        "  for i in range(P[:,0].size):\n",
        "\n",
        "    if Dirichlet(a, b, P[i,:], tol, pb):\n",
        "      # node i belongs to Dirichlet boundary\n",
        "      v[i] = P[:,0].size - n_diric - 1\n",
        "      n_diric += 1\n",
        "    elif Neumann(a, b, P[i,:], tol, pb):\n",
        "      # node i belongs to Neumann boundary\n",
        "      v[i] = P[:,0].size - n_bndry + n_neumn \n",
        "      n_neumn += 1\n",
        "    else: \n",
        "      # node i is internal \n",
        "      v[i] = n_inter\n",
        "      n_inter += 1\n",
        "    \n",
        "    # update nodes matrix P \n",
        "    P_new[v[i],:] = P[i,:]\n",
        "\n",
        "  # update elements matrix T\n",
        "  T_new[:,:] = v[T[:,:]]\n",
        "\n",
        "  return P_new, T_new, n_diric, n_neumn"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZP5sTx4UDe0o"
      },
      "source": [
        "def plot_mesh(P, E):\n",
        "  for i in range(E[:,0].size):\n",
        "    x = [P[E[i, E[0,:].size-1], 0], P[E[i, 0], 0]]\n",
        "    y = [P[E[i, E[0,:].size-1], 1], P[E[i, 0], 1]]\n",
        "    plt.plot(x, y, \"k\")\n",
        "    for j in range(E[0,:].size-1):\n",
        "      x = [P[E[i, j], 0], P[E[i, j+1], 0]]\n",
        "      y = [P[E[i, j], 1], P[E[i, j+1], 1]]\n",
        "      plt.plot(x, y, \"k\")\n",
        "        \n",
        "  plt.show()"
      ],
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GtvGwAfYGNnB",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 591
        },
        "outputId": "e4cc3b78-4003-4951-f752-ec602e3d044a"
      },
      "source": [
        "P, T, n_diric, n_neumn = ordered_mesh(1, 1, 4, 2, 2)\n",
        "plot_mesh(P,T)\n",
        "print(P)"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:7: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.\n",
            "Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations\n",
            "  import sys\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deVjUVfs/8PcZRkRARUIeFxSXzFwyU0xZBBQRQQJ3RY0kV8yFNHdDIrc090Rzw8glVDQJcQdkVTPtccn8uT3mWuSKC7Kd3x8wNPlV54MOnDPD/bour4tlYt4DejecmbnfjHMOQgghhk8lOgAhhBD9oIFOCCFGggY6IYQYCRrohBBiJGigE0KIkVCLumIbGxter149UVdPCCEG6Zdffvmbc179eZ8TNtDr1auHY8eOibp6QggxSIyxKy/6HB25EEKIkaCBTgghRoIGOiGEGAka6IQQYiRooBNCiJHQOdAZY+sYY38xxk6/4POMMbaUMXaBMXaSMdZK/zEJIYToouQe+noAXV7yeW8AjYr+DAOw4vVjEUIIKSmdz0PnnCczxuq95CL+AKJ44R7ew4wxK8ZYTc75TT1l/JeGDRvi6tWrcHR0BGOsNK7CqKSkpKCgoADt27eHSkUnbLqkpqYiPz8fTk5OqFChgug40svJyUFGRgZMTEzg4uIiOo5BuHDhAqpXr44TJ07o/Wvr4194bQBXtd6/VvSx/4MxNowxdowxdiwzM/OVruz69evIzc3FmTNnXum/L28KCgoAAL/++qvgJIYhPz8fAHD8+HFQV4Bud+/eBVD4faPvlzIPHz7Eq84/Xcr0laKc81UAVgGAg4PDK/3033//fRw+fBi3b99G3759ERwcrNeMxsbKygpZWVnIysqCs7MzZs2aJTqS1KysrPDw4UNkZ2ejXr16+O6770RHkl6FChWQl5cHU1NT7N+/X3Qc6bm7u5fa19bHPfTrAOpovW9X9LFSoVKp0KpVK6jVaowaNYruqStgYWGBSpUqYfbs2UhKShIdR3rm5uawtrZGVFQUNm7cKDqO9CpVqgTGGA4cOID58+eLjlOu6WOgxwIILHq2SzsA90vr/FzDzMwMGzduLD4bzsnJKc2rM3gqlQp79uwBYww+Pj64d++e6EhSU6lUSE1NhYmJCQYNGoTLly+LjiQ1lUoFS0tLVKxYEZMmTaIdTQIpedriZgAZABozxq4xxgYzxkYwxkYUXSQewCUAFwCsBjCy1NJq6dOnD4KCgnD37l107NixLK7SoLm6uiI0NBRPnjxBu3btRMeRXpMmTRAREYG8vDw4OjoiLy9PdCSpqVQqbN++HZxzdOzYEY8fPxYdqVzSOdA55wGc85qc8wqcczvO+VrO+UrO+cqiz3PO+Sec84ac83c452X2v+d169bhrbfeQlpaGsLCwsrqag1WWFgYnJyccO7cOQwZMkR0HOkNGzYMPXr0wJ9//glfX1/RcaTn4+ODkJAQZGVl0TNeBDH457FlZGSgUqVKCA8PR3Jysug40ktMTES1atWwdu1abNmyRXQc6W3duhV16tTB3r17sXDhQtFxpLdo0SK8++67OHHiBMaNGyc6Trlj8APd2toa8fHx4JyjS5cudD6sg6mpKVJSUqBSqTBgwABcufLC1coEhUcJhw8fhqmpKT777DMcP35cdCTppaamwtLSEosWLcLu3btFxylXDH6gA4VPA5o6dSqePHkCR0dH0XGk16xZM3zzzTfF58Oa56qT56tVqxa2bNkCzjnc3d3pfFgHS0tLJCQkgDGG7t2746+//hIdqdwwioEOALNmzULbtm3x+++/Y9iwYaLjSC84OBjdunXDzZs38cEHH4iOIz1/f3+MHDkSWVlZcHV1FR1Hem3atMHs2bPx9OlTtG3blu40lBGjGegAkJSUBCsrK6xevRoxMTGi40gvJiYGdnZ2iI+Px7Jly0THkd7y5cvxzjvv4JdffsGECRNEx5He5MmT4eHhgf/9738YMGCA6DjlglENdDMzs+Lz4YCAAPzxxx+iI0lNpVIhIyMDFSpUQEhICK0HUCA1NRUWFhb4+uuvsXfvXtFxpLdnzx5Ur14dP/zwA73qtgwY1UAHgObNm2PJkiXIzc2l82EF7OzsEB0djYKCAri7uyM7O1t0JKlVqVIF+/fvB2MM3bp1w99//y06ktTUanXx8q7Bgwfj/PnzoiMZNaMb6AAwatQo+Pr64saNG+jevbvoONLr3r07RowYgfv378PNzU10HOk5Ojriyy+/RHZ2Np0PK9CwYUOsW7cO+fn5cHZ2phdplSKjHOgAsHPnTtSqVQuxsbGIiIgQHUd6K1asQNOmTXH06FFMnTpVdBzpTZs2DW5ubrh06RICAwNFx5FeYGAgAgICkJmZCS8vL9FxjJbRDnTt8+HRo0fj9OnnFi4RLRkZGTA3N8ecOXOQkJAgOo709u3bhzfeeAMbN27Ehg0bRMeR3oYNG1C/fn0kJCRg7ty5ouMYJaMd6ABQt25dbNq0qXiJF50Pv1yVKlWwd+9eMMbQtWtX3LlzR3QkqZmamiItLQ0mJiYICgrCxYsXRUeSmuZFWhUrVsTUqVNx5MgR0ZGMjlEPdADo1asXhgwZgnv37tESLwVcXFwwY8YMZGdn0xIvBRo3boxvv/0WeXl5cHJyovNhHWxtbbFjxw5wztGpUyc8fPhQdCSjYvQDHQBWr16Nt99+GxkZGZg2bZroONKbMWMGXFxccP78eQwaNEh0HOkNHjwYvXr1wl9//QUfHx/RcaTn7e2NcePG4eHDh7TES8/KxUAH/lniRSUPyhw8eBDW1tb47rvvsHnzZtFxpBcdHY26deti//79VPKgwIIFC/Dee+/hv//9L8aOHSs6jtEoNwPdysqKSh5KwNTUtLjkITAwkJZ46aBSqXDkyBGYmppSyYNCqampqFy5MpYuXYq4uDjRcYxCuRnoQGHJw/Tp06nkQSHtkod27drR+bAONWrUQExMDJU8KGRubl68xKtnz564deuW6EgGr1wNdAAIDw+Ho6MjlTwopCl5uHXrFi3xUsDX1xdjxoyhkgeFHBwc8NVXXyEnJwft2rWjF2m9pnI30AEgISGBSh5KQFPysGfPHixevFh0HOktWbIELVq0oJIHhSZMmABPT09cuXIF/fr1Ex3HoJXLga69xGvgwIG0xEsH7ZKH8ePH0xIvBdLS0qjkoQTi4+Nha2uLrVu3Yu3ataLjGKxyOdCBf0oecnNz6Vc9BTQlDwUFBXBzc6PzYR0sLS1x4MABKnlQSK1W4/Dhw1Cr1Rg+fDjOnTsnOpJBKrcDHSgsefD398fNmzfh5+cnOo70NCUPDx48oCVeCrRt27a45IHuNOhWv359REZGFi/xysnJER3J4JTrgQ4A27dvR+3atbFr1y4qeVBg+fLlaN68OY4dO4ZJkyaJjiO9yZMno0OHDrh8+TIGDhwoOo70Bg4ciAEDBuD27du0xOsVlPuBrjkf1pQ8nDx5UnQk6aWlpcHCwgLz5s3D/v37RceR3r59+2BjY4PNmzcjKipKdBzpRUVFoUGDBkhKSsKsWbNExzEo5X6gA4UlD5s3b0ZBQQFcXV1piZcO2ku8/Pz8qORBB7VajfT0dJiYmODjjz+mkgcdNC/SMjMzw+eff46MjAzRkQwGDfQiPXv2xPDhw6nkQSFnZ2eEh4cXL/Gi8+GXa9SoEVavXk0lDwrZ2Njgxx9/BOccnp6etMRLIRroWlauXIkmTZpQyYNC06dPh6urKy5evEhLvBQICgpC3759kZmZiS5duoiOIz0vLy9MmDABjx49gpOTk+g4BoEG+jPS09Op5KEE9u/fjzfeeAPff/89Nm7cKDqO9DZt2oR69erh4MGDVPKgwLx589C6dWucOnUKn3zyieg40qOB/gwrKysqeSgB7ZKHQYMG4fLly6IjSU1zPkwlD8olJyejSpUqiIiIQGxsrOg4UqOB/hwuLi4IDQ2lkgeFtEseaImXblTyUDLm5uZITEwEYwy9e/fGjRs3REeSFg30FwgLC4OzszPOnz+PoKAg0XGkp13y0LVrV9FxpOft7Y1PP/2USh4UatWqFb7++mta4qUDDfSX0CzxWr9+PaKjo0XHkZ6m5GHfvn1YuHCh6DjSW7hwIVq2bEklDwqNGzcOXl5euHr1Kvr06SM6jpQUDXTGWBfG2DnG2AXG2OTnfL4uYyyRMXaCMXaSMWYUPVya82HNEi8qeXg57ZKHzz77jEoeFEhJSYGlpSWWLl2K+Ph40XGkFxcXh//85z+IiYnBmjVrRMeRjs6BzhgzAbAcgDeApgACGGNNn7nYdABbOOfvAegHIELfQUVp0qQJVqxYQefDClHJQ8lYWloWlzxo9s6TF1Or1cjIyIBarcaIESNw9uxZ0ZGkouQe+vsALnDOL3HOcwD8AMD/mctwAFWK3q4KwKgetRg2bBi6d+9OJQ8K+fr6YvTo0cjKykL79u1Fx5FemzZtMHfuXFripVD9+vURFRWF/Px8uLi40BIvLUoGem0AV7Xev1b0MW1hAAYyxq4BiAcw+nlfiDE2jDF2jDF2LDMz8xXiirNt2zYqeSiBpUuX4p133sHx48cxfvx40XGkN3HiRHh4eODKlSvo37+/6DjSCwgIQGBgIO7cuQNPT0/RcaShrwdFAwCs55zbAfAB8D1j7P98bc75Ks65A+fcoXr16nq66rJBJQ8ll56eDgsLCyxcuBB79+4VHUd6e/bsQfXq1REdHY3IyEjRcaQXGRmJN998E8nJyZg5c6boOFJQMtCvA6ij9b5d0ce0DQawBQA45xkAzADY6COgTKjkoWQsLS2xf/9+MMbg7+9PJQ86aM6HTUxMMHToUFripYNKpUJGRgbMzMwQGhqKtLQ00ZGEUzLQfwbQiDFWnzFmisIHPZ99udYfADwAgDHWBIUD3bDOVBTy9/dHcHAwlTwo5OjoiFmzZuHp06dwdHSk82EdGjZsiHXr1iE/Px+Ojo50PqyDjY0N4uLiwDmHl5cXHjx4IDqSUDoHOuc8D8AoAHsBnEXhs1nOMMbCGWOamp/xAIYyxv4LYDOAQZxzXlqhRYuIiECzZs2o5EGhKVOmoEOHDrh06RI+/PBD0XGkFxgYiP79+1PJg0IeHh6YOHEiLfECAM65kD+tW7fmr8LNzY27ubm90n+rT/fv3+cWFhYcAN+3b5/oOC9UtWpVXrVqVdExeG5uLrexseEA+HfffSc6zgvJ8v3Kz8/nDRo04AD4zJkzRcd5KVm+Zw4ODhwADw4OFh3lpV53hgE4xl8wV+mVoq+ISh5K5tmSh4sXL4qOJDXtJV5U8qBMSkoKqlSpghUrVmDnzp2i4whBA/01ODs7IywsDNnZ2XQ+rIB2yYOjoyO9SEsHGxsb7Ny5k0oeFDIzM8OhQ4egUqnQp0+fcrnEiwb6awoNDUX79u1x4cIFWuKlgHbJg7e3t+g40vPy8sK4cePofFihli1bYsGCBcjJyUHbtm3L3Z0sGuh6cODAAVhbWyMqKopKHhTYtGkT7O3tceDAAcybN090HOktWLAArVq1wqlTpzBmzBjRcaQXEhICb29vXLt2Db169RIdp0zRQNcDU1PT4vNhKnnQTaVS4ejRo6hYsSImT56Mn3/+WXQk6aWkpKBy5cpYtmwZlTwoEBcXhxo1amDHjh349ttvRccpMzTQ9aRx48ZYuXIlLfFSyNbWFtu3by9e4kXnwy9nbm5evMSrd+/etMRLB80ru9VqNUaOHIkzZ86IjlQmaKDr0ZAhQ9CzZ08qeVDIx8cHISEhePjwIS3xUsDBwQHz588vt+fDJWVvb48NGzagoKAA7du3Lxcv0qKBrmdbtmyhkocSWLRoEd599138+uuv+PTTT0XHkd748ePh6emJP/74A3379hUdR3p9+/bFoEGDcPfuXXh4eIiOU+pooOsZlTyUXGpqKiwtLbF48WIqeVAgPj4etra22LZtG5U8KBAZGYlGjRohNTUVX3zxheg4pYoGeimoUaMGtm7dSiUPCj1b8kBLvF5OrVYXnw9TyYMyhw8fRqVKlfDFF18gNTVVdJxSQwO9lPj5+WHUqFHIysqCq6ur6DjS0y55oPNh3erXr4/169dTyYNC1tbW/1ride/ePdGRSgUN9FK0bNkyvPPOO/jll18wYcIE0XGkpyl5+N///ocBAwaIjiO9AQMG4MMPP6SSB4U6duyIKVOm4PHjx0b7Ii0a6KUsNTUVFhYW+Prrr6nkQQFNycMPP/yA7777TnQc6a1fvx4NGzakkgeFZs+ejffffx9nz57FiBEjRMfROxropaxKlSrFJQ/dunWjJV46aJc8DB48mEoedNA835pKHpQ7dOgQqlatim+//RYxMTGi4+gVDfQy4OjoiC+//BLZ2dl0PqyAdsmDk5MTvUhLBxsbG8TGxlLJg0JmZmZITk6GSqVCQEAArl27JjqS3tBALyPTpk2Dm5sblTwoFBgYiICAAPz999/o3Lmz6DjS8/T0xIQJE/Do0SM4OzuLjiO9Fi1aYPHixcjNzUW7du2M5k4WDfQytG/fPrzxxhvYtGkToqKiRMeR3oYNG1C/fn0kJiZizpw5ouNIb968eWjdujVOnz6NTz75RHQc6Y0ePRpdu3bF9evX0aNHD9Fx9IIGehkyNTVFWloalTwopDkfrlixIqZNm4YjR46IjiS95ORkVKlSBREREeW25KEkYmNjUbNmTezcuRMrVqwQHee10UAvY40bN6aShxKwtbXFjh07wDmHh4cHLfHSwdzcHImJieW65KEkNHcaKlSogFGjRhn8Ei8a6AIEBQWhd+/eVPKgkLe3d3HJA50P69aqVavikgdjOh8uLXXr1jWaJV400AX54Ycfikse5s+fLzqO9BYsWID33nsPJ0+epJIHBUJCQuDl5YWrV6+id+/eouNIr0+fPhg8eDDu3r2LDh06iI7zymigC6L5Vc/U1BSTJk2iJV4KpKamFpc8xMXFiY4jPU3Jw/bt27Fq1SrRcaS3Zs0aNG7cGOnp6QgNDRUd55XQQBeoRo0aiImJoSVeCmmXPPTs2ZNKHnTQXuIVHBxs8OfDZUGzxGvmzJlISkoSHafEaKAL5uvrizFjxiArKwsuLi6i40hPu+SBzod1s7e3R1RUlFGcD5cFKyur4hXOPj4+BrfEiwa6BJYsWYIWLVrgxIkTGDdunOg40tOUPFy5cgX9+vUTHUd6AQEB+Oijj8pNycPrcnd3x7Rp0/DkyRM4OjqKjlMiNNAlkZaWBktLSyxatIhKHhTQlDxs3bqVSh4UWL9+fbkpedCHL7/8Eu3atcPvv/+OoUOHio6jGA10SVhaWuLAgQNU8qDQsyUP586dEx1JepolXsZe8qAviYmJsLKywpo1a7Bt2zbRcRShgS6Rtm3bYvbs2VTyoFD9+vURGRmJ/Px8ODs70/mwDs+WPNASr5czMzNDSkoKVCoV+vfvjz/++EN0JJ1ooEtm8uTJ6NChA5U8KDRw4EAMGDAAt2/fpiVeCnh4eGDy5Ml4/PixwZ0Pi9C8eXMsWbIEubm5cHR0lP5OFg10Ce3btw82NjZU8qBQVFQUGjRogEOHDlHJgwJz5sxBmzZt8NtvvyE4OFh0HOmNGjUKvr6+uHHjBvz9/UXHeSlFA50x1oUxdo4xdoExNvkFl+nDGPuNMXaGMbZJvzHLF7VajfT0dCp5UEilUuHIkSPFJQ8ZGRmiI0kvOTkZVatWxcqVK7Fjxw7RcaS3c+dO1KpVC3Fxcfjmm29Ex3khnQOdMWYCYDkAbwBNAQQwxpo+c5lGAKYAcOacNwMQUgpZy5VGjRph7dq1xefDtMTr5WxsbPDjjz+Ccw5PT086H9bBzMwMSUlJUKlU6Nu3r1GVPJQGzZ2GChUqYOzYsTh58qToSM+l5B76+wAucM4vcc5zAPwA4NnfO4YCWM45vwsAnHN6ioYefPTRR+jXrx8yMzPh5eUlOo70vLy8ikse6EVaurVs2dIoSx5Ki52dHTZv3oyCggK4ubkhOztbdKT/Q8lArw3gqtb714o+pu0tAG8xxtIYY4cZY12e94UYY8MYY8cYY8cyMzNfLXE5s3HjRtSrVw8JCQmYO3eu6DjS05Q8nDp1ikoeFBg9ejR8fHxw/fp19OzZU3Qc6fXs2RNDhw7FvXv3pFzipa8HRdUAGgFwBxAAYDVjzOrZC3HOV3HOHTjnDtWrV9fTVRs3za96FStWxNSpU6nkQQHtkofY2FjRcaT3008/oWbNmvjxxx+NouShtK1atQpvv/02Dh8+jGnTpomO8y9KBvp1AHW03rcr+pi2awBiOee5nPPLAP4fCgc80QPtkodOnTpRyYMO2iUPvXv3ppIHHVQqFTIyMqBWq42i5KEsZGRkoFKlSpg9ezYSEhJExymmZKD/DKARY6w+Y8wUQD8Az97t+RGF987BGLNB4RHMJT3mLPe8vb3x6aef4uHDh1TyoACVPJSMvb09Nm7cSEu8FLKyssKePXvAGIOvry/u3LkjOhIABQOdc54HYBSAvQDOAtjCOT/DGAtnjPkVXWwvgNuMsd8AJAKYwDm/XVqhy6uFCxeiZcuWVPKgEJU8lEyfPn0QFBRk8CUPZcXV1RWhoaFyLfHinAv507p1a/4q3NzcuJub2yv9t8YgKyuLW1pacgD8p59+0nn5qlWr8qpVq5ZBMjnl5uby//znPxwA//bbb3Vevrx/vzjn/K233uIA+Oeff67o8uX9e+bk5MQB8KCgIEWXf90ZBuAYf8FcpVeKGhhLS0sqeSgBtVpdfD48cuRInD17VnQk6WnOh2fOnInk5GTRcaSXmJiIatWqITIyEtHR0UKz0EA3QG3atMHcuXPpfFih+vXrIyoqCvn5+XBxcaHzYR2sra0RHx8Pzjm6dOlicCUPZc3U1LR4idfAgQNx5coVYVlooBuoiRMnwsPDg0oeFAoICEBgYCDu3LmDTp06iY4jPXd3d0ydOlWu82GJNWvWDN988w3y8vKELvGigW7A9uzZg+rVq2Pr1q1Yu3at6DjSi4yMxJtvvomUlBSEh4eLjiO9WbNmoW3btvj9998xbNgw0XGkFxwcDH9/f9y8eRMffPCBkAw00A2Y5nzYxMQEw4cPp5IHHTTPtzYzM0NYWBjS0tJER5JeUlISrKyssHr1asTExIiOI73t27ejdu3aiI+Px7Jly8r8+mmgG7iGDRti3bp1VPKgkI2NTXHJQ+fOnWmJlw5mZmY4dOgQVCoVAgICDKLkQSSVSoXDhw+jQoUKCAkJwa+//lq211+m10ZKRWBgIPr374/bt2/TEi8FtEsenJycRMeRXosWLQyq5EE0Ozs7REdHo6CgAO7u7mW6xIsGupH4/vvv0aBBAyQlJWHWrFmi40hPU/Jw5swZKnlQwJBKHmTQvXt3jBgxAvfv34erq2uZXS8NdCOhvcTr888/p5IHBajkoWQMpeRBFitWrEDTpk3x888/Y8qUKWVynTTQjYiNjQ127txJJQ8KUclDyWgeVJa95EEmGRkZMDc3x9y5c3Hw4MFSvz4a6EbGy8sL48aNo5IHhVq2bImFCxcWlzyQl6tbty42bdokdcmDTKpUqYK9e/eW2RIvGuhGaMGCBWjVqhVOnTqFJ0+eiI4jvbFjx8Lb2xvXr1/Ho0ePRMeRXq9evTBkyBBpSx5k4+LighkzZiA7O7vU7zTQQDdSKSkpqFy5MnJycpCbmys6jvTi4uJQo0YN5OXl0VM/FVi9enVxyQPdS9dtxowZcHFxwfnz5/H777+X2vWoS+0rl5KUlBQUFBSgatWqYIyJjiM1zdPLHj9+jMqVK8PExERwIrlpvl9PnjyBhYUFKlSoIDiR3DTfr6dPn+Lp06ewsvo/JWVES0FBARhj+PPPP1FaFZwGew+dWnt0U6v/+f81HSXoplL988/h8ePHApMYBpVKhYoVKxa/X7jZlbyISqWCubk5gNL7XhncPfT27dvj2LFjePToEYKCgrBo0SLRkaRmZWWFhw8fIj8/H56enti6davoSFKzsrLCo0ePkJeXh+bNmyM1NVV0JOlVrFgROTk5qFmzZqkeJxgLFxeXf93Z0ieDvIfesmVLWFpaYvHixYiPjxcdR3rm5uawtbXFtm3bsGbNGtFxpGdhYYFGjRohLS0NYWFhouNIr1KlSjAxMcG5c+cwZMgQ0XGkV1rDHDDQga5Wq4tLHnr06EElDzpo9kuo1WqMGDGCSh4UOHz4MCpVqoTw8HC6l66Aubk5qlWrhrVr12LLli2i45RbBjnQgX9KHp4+fUolDwrUr18f69evp5IHhaytrYuXeHl5eVHJgw4qlaq45GHAgAFCSx7KM4Md6EBhyUOnTp1w5coV9O/fX3Qc6Q0YMKC45MHT01N0HOl17NgRU6ZMwePHj6nkQQFZSh7KM4Me6ACwe/du2NraIjo6GpGRkaLjSE9T8pCcnIyZM2eKjiO92bNnU8lDCQQHB6Nbt25CSx7KM4Mf6Gq1Gunp6VCr1Rg6dCjOnz8vOpLUtEseQkNDqeRBASp5KJmYmBjY2dkJK3kozwx+oAOFJQ9r165Ffn4+nJyckJeXJzqS1LRLHry8vGiJlw5U8lAy2ku8RJQ8lGdGMdCBf0oe/v77b3Tu3Fl0HOl5eHhg4sSJePToEZydnUXHkR6VPJSMyJKH8sxoBjrwT8lDYmIi5syZIzqO9L766is4ODjg9OnT+OSTT0THkd6oUaPQtWtX3LhxA927dxcdR3raJQ9ubm6i45QLRjXQtUsepk2bhiNHjoiOJL2UlBRUqVIFERER2Llzp+g40ouNjUXNmjURGxuLiIgI0XGkpyl5OHr0KKZOnSo6jtEzqoEO/LvkwcPDg3a+6KB9PtynTx/cuHFDdCSpaZcAjx49GqdPnxYdSXqakoc5c+YgISFBdByjZnQDHfh3yQOdD+vWsmVLLFiwADk5OfQiLQXq1q2LDRs2oKCgAO3bt6fzYR20Sx66du1a6iUP5ZlRDnTgn5KHkydPYsyYMaLjSC8kJARdunTB1atX0bt3b9FxpNenTx8MHjwY9+7dQ8eOHUXHkV5ZljyUZ0Y70IF/Sh6WLVuGuLg40XGkt2vXLtSoUQPbt2/HqlWrRMeR3skiTcAAABsQSURBVJo1a9C4cWNkZGTg888/Fx1HetolD0FBQaLjGCWjHujm5ubFS7x69uxJS7x00F7iFRwcjDNnzoiOJD3NEq9Zs2YhKSlJdBzpHTx4ENbW1li/fj2io6NFxzE6igY6Y6wLY+wcY+wCY2zySy7XkzHGGWMO+ov4ehwcHDB//nzk5OSgbdu2dD6sg729/b/Oh2mJ18tZWVlhz549AAAfHx9a4qWDqakpUlNToVKpMHDgQFripWc6BzpjzATAcgDeAJoCCGCMNX3O5SoDGAtAuucKjh8/Hp6envjjjz/Qt29f0XGk17dvXwwaNAh3796Fh4eH6DjSc3V1xfTp0/HkyRM6H1agSZMmWLFiBfLy8tCuXTt6ZbceKbmH/j6AC5zzS5zzHAA/APB/zuW+BPAVACkf8o+Pj6eShxKIjIxEo0aNkJqaii+++EJ0HOmFh4fD0dGRSh4UGjZsWHGXAS3x0h8lA702gKta718r+lgxxlgrAHU457te9oUYY8MYY8cYY8dKqyT1RdRq9b9KHs6dO1em12+INOfDX3zxBZU8KJCQkAArKysqeVBo69atqFOnDvbs2YPFixeLjmMUXvtBUcaYCsBCAON1XZZzvopz7sA5d6hevfrrXnWJaZc8ODs70/mwDs+WPNASr5czMzMrLnkYOHAgLfHSQfMgvKmpKcaPH09LvPRAyUC/DqCO1vt2RR/TqAygOYAkxtj/ALQDECvTA6PaBgwYgA8//BC3b9+mJV4KUMlDyTRv3hzLli1Dbm4uvUhLgVq1amHLli0oKCiAm5sbHj9+LDqSQVMy0H8G0IgxVp8xZgqgH4BYzSc55/c55zac83qc83oADgPw45wfK5XEerB+/Xo0bNgQhw4dopIHBWbPno33338fv/32G4KDg0XHkd7IkSPh5+eHmzdvws/PT3Qc6fn7+2PkyJF48OABLfF6TToHOuc8D8AoAHsBnAWwhXN+hjEWzhgzyL+tml/1NCUPGRkZoiNJ79ChQ6hatSpWrlyJHTt2iI4jvR07dqB27drYtWsXlTwosHz5cjRv3hzHjh3DpEmTRMcxWIrO0Dnn8ZzztzjnDTnns4o+Fso5j33OZd1lvneuYWNjg9jYWHDO4enpSefDOpiZmSE5ORkqlQp9+/bFtWvXREeSmvYSr5CQEJw8eVJ0JOmlpaXBwsIC8+bNw/79+0XHMUhG/UpRXTw9PTFhwgQ8evQILi4uouNIr0WLFli8eDGdDytkZ2eHzZs3o6CgAK6urrTESwftJV5+fn74+++/RUcyOOV6oAPAvHnz0Lp1a5w6dYpKHhQYPXo0fHx8cP36dfTo0UN0HOn17NkTQ4cOxf379+Hu7i46jvScnZ0RHh6O7OxsaoZ6BeV+oANAcnIylTyUwE8//YSaNWti586dWLFiheg40lu1ahWaNGmCI0eOYNq0aaLjSG/69OlwdXXFhQsXaIlXCdFAR+ESr8TERCp5UEj7fHjUqFG0xEuB9PR0mJubY/bs2VTyoMD+/fvxxhtvICoqChs3bhQdx2DQQC/SqlUrKnkogWdLHuhFWi9nZWVVfD7s6+tLJQ86mJqaIi0tDSYmJhg0aBAuX74sOpJBoIGuJSQkBF5eXlTyoJCm5OHu3bvo0KGD6DjSc3FxQWhoKC3xUqhx48ZYuXIlLfEqARroz4iLi6OShxLQlDykp6cjNDRUdBzphYWFwdnZmUoeFBoyZAh69eqFv/76C127dhUdR3o00J+hvcRr5MiROHv2rOhI0tMs8Zo5cyaSk5NFx5FeQkICqlWrRiUPCkVHR6Nu3brYt28fFi5cKDqO1GigP4e9vT2ioqKQn58PFxcXOh/WwcrKCvHx8QCALl26UMmDDqampv9a4kUlDy+nUqlw5MgRmJqa4rPPPsOxY9K/blEYGugvEBAQgI8++gh37txBp06dRMeRnru7O6ZNm4YnT57QEi8FmjVrhoiIiOLzYXoQ/uVq1KiBmJgYcM7RsWNHWuL1AjTQX2L9+vVo1KgRUlJSEB4eLjqO9L788ku0a9cOv//+O4YNGyY6jvSGDx+Obt264datW/D19RUdR3q+vr4YPXo0srKy4OrqKjqOlGig66BZ4hUWFoa0tDTRcaSXmJgIKysrrF69GjExMaLjSC8mJgZ2dnbYvXs3lixZIjqO9JYuXYp33nkHv/zyCyZMmCA6jnRooOugXfLQuXNnWuKlg3bJQ0BAAJU86KBSqZCRkYEKFSpg3LhxVPKgQHp6OiwsLPD1119j7969ouNIhQa6Ah4eHpg8eTKVPCjUvHlzLFmyBLm5ubSPQwE7OztER0ejoKAA7u7utMRLB0tLS+zfvx+MMXTr1o2WeGmhga7QnDlz0KZNGyp5UGjUqFHw8/PDjRs34O//vE5xoq179+4IDg7G/fv36XxYAUdHR8yaNQvZ2dlo27Yt3WkoQgO9BJKTk6nkoQQ0JQ9xcXH45ptvRMeRXkREBJo1a4aff/4ZU6ZMER1HelOmTIG7uzsuXbqEwMBA0XGkQAO9BMzMzJCUlEQlDwppL/EaO3YslTwooDkfnjt3Lg4ePCg6jvQ0S7w2btyIDRs2iI4jHA30EmrZsiUWLlxIJQ8KaZc8uLm50fmwDtolD76+vnQ+rINarS5e4hUUFISLFy+KjiQUDfRXMHbsWHh7e+P69evo2bOn6DjS05Q83Lt3j5Z4KeDs7IywsDAqeVCocePGWL16NfLy8uDo6Fiul3jRQH9FcXFxqFmzJn788UcqeVBg1apVePvtt3H48GEqeVAgNDQU7du3p5IHhYKCgtC3b19kZmbC29tbdBxhaKC/Is3zh9VqNZU8KJSRkYFKlSph9uzZSEpKEh1HegcOHIC1tTWVPCi0adMm2Nvb48CBA5g/f77oOELQQH8N9vb22LhxI5U8KGRlZYU9e/aAMQYfHx9a4qWDqakpUlNTqeRBIZVKhaNHj6JixYqYNGlSuVziRQP9NfXp0wdBQUG4e/cuOnbsKDqO9FxdXankoQSaNGlSvMSrvJ8PK2Fra4vt27eX2yVeNND1YN26dXjrrbeQlpaGsLAw0XGkFxYWBicnJ5w7dw4ff/yx6DjSGzZsGHr06IE///yTlngp4OPjg5CQEGRlZcHFxUV0nDJFA11PNOfD4eHhVPKgQGJiIqpVq4bIyEgqeVBg69atqFOnDvbu3UslDwosWrQI7777Lk6cOIFx48aJjlNmaKDribW1NeLj48E5p5IHBajkoWQ0L9LSlDwcP35cdCTppaamwtLSEosWLcLu3btFxykTNND1yN3dHVOnTqWSB4W0Sx7o+da61apVC1u3bgXnHB06dCh358MlZWlpiYSEBDDG0L17d/z111+iI5U6Guh6NmvWLLRt25ZKHhTSlDzcvHkTH3zwgeg40vPz88OoUaPw4MEDWuKlQJs2bTB79mw8ffq0XLyymwZ6KUhKSqKShxLQlDzEx8dj2bJlouNIb9myZcUlDxMnThQdR3qTJ09Gx44dcfnyZQwcOFB0nFJFA70UmJmZ4dChQ1TyoJB2yUNISAiVPCiQmpoKCwsLzJ8/H/v37xcdR3p79+5F9erVsXnzZkRFRYmOU2oUDXTGWBfG2DnG2AXG2OTnfH4cY+w3xthJxthBxpi9/qMalhYtWlDJQwlQyUPJVKlSpbjkwc/Pj5Z46aC9xOvjjz/G+fPnRUcqFToHOmPMBMByAN4AmgIIYIw1feZiJwA4cM5bANgGYJ6+gxqiUaNGoWvXrrhx4wa6d+8uOo70unfvjhEjRuD+/ftwc3MTHUd6jo6O+PLLL6nkQaFGjRph7dq1yM/Ph7Ozs1G+SEvJPfT3AVzgnF/inOcA+AHAvypoOOeJnHPNQ+6HAdjpN6bhio2NRa1atRAbG4uIiAjRcaS3YsUKNG3aFEePHsXUqVNFx5HetGnT4ObmRiUPCn300UcICAhAZmYmvLy8RMfROyUDvTaAq1rvXyv62IsMBvDcJ30yxoYxxo4xxo5lZmYqT2nAtM+HR48ejdOnT4uOJL2MjAyYm5tjzpw5SEhIEB1Hevv27aOShxLYsGED6tevj4SEBMydO1d0HL3S64OijLGBABwAPHfVGed8FefcgXPuUL16dX1etdTq1q2LTZs2FS/xovPhl9MueejatSvu3LkjOpLUTE1NqeShBDQv0qpYsSKmTp2KI0eOiI6kN0oG+nUAdbTetyv62L8wxjoBmAbAj3P+VD/xjEevXr0wZMgQ3Lt3j5Z4KeDi4oIZM2YgOzublngp0LhxY3z77bfIy8uDk5OTUZ4P65OtrS127NgBzjk6deqEhw8fio6kF0oG+s8AGjHG6jPGTAH0AxCrfQHG2HsAvkXhMDf+l2O9otWrV+Ptt99GRkYGlTwoMGPGDLi4uOD8+fMYNGiQ6DjSGzx4MHr16oW//voLPj4+ouNIz9vbG+PGjcPDhw+NZomXzoHOOc8DMArAXgBnAWzhnJ9hjIUzxvyKLjYfgCWArYyxXxljsS/4cuUelTyUzMGDB2FtbY3vvvsOmzdvFh1HetHR0bC3t8f+/fvLbclDSSxYsADvvfce/vvf/2Ls2LGi47w2RWfonPN4zvlbnPOGnPNZRR8L5ZzHFr3diXP+H855y6I/fi//iuUXlTyUjHbJQ2BgIC3x0kF7iVd5LXkoqdTUVFSuXBlLly5FfHy86DivhV4pKoCrqyumT59OJQ8KaZc8tGvXjs6HdahRowZiYmLKbclDSZmbmxcv8erRowdu3bolOtIro4EuSHh4OBwdHXHu3DkMGTJEdBzpaUoebt26RUu8FPD19cWYMWOQlZWF9u3bi44jPQcHB3z11VcGv8SLBrpACQkJsLKywtq1a7FlyxbRcaSnKXnYs2cPFi9eLDqO9JYsWYIWLVrg+PHjGD9+vOg40pswYQI6deqEK1euoH///qLjvBIa6AKZmZn9q+SBlni9nPb58Pjx42mJlwJpaWmwtLTEwoULsXfvXtFxpLd7927Y2toiOjoakZGRouOUGA10wZo3b45ly5YhNzfXoH/VKyu1atXCli1bUFBQADc3Nzof1sHS0hIHDhwAYwz+/v7louThdajVaqSnp0OtVmPo0KE4d+6c6EglQgNdAiNHjoSfnx9u3rwJPz96gpAu/v7+GDlyJB48eEBLvBRo27ZtuSp5eF0NGzb81xKvnJwc0ZEUo4EuiR07dqB27drYtWsXlTwosHz5cjRv3hzHjh3DpEmTRMeR3uTJk9GhQ4dyUfKgD4GBgRgwYABu375tUEu8aKBLQnM+rCl5OHnypOhI0ktLS4OFhQXmzZtHJQ8K7Nu3DzY2NkZf8qAvUVFRaNCgAZKSkjBr1izRcRShgS4ROzs7bN68GQUFBXB1daUlXjpoL/GikgfdNOfDxl7yoC8qlQpHjhyBmZkZPv/8c2RkZIiOpBMNdMn07NkTw4cPp5IHhZydnREeHl68xIvOh1+uPJQ86JONjQ1+/PFHcM7h6ekp/RIvGugSWrlyJZo0aUIlDwpNnz4dbm5uuHjxIi3xUuCjjz5Cv379kJmZiS5duoiOIz0vLy9MmDABjx49gpOTk+g4L0UDXVLp6elU8lACmpKH77//Hhs3bhQdR3obN25EvXr1cPDgQcybR42RusybNw+tW7fGqVOnMGbMGNFxXogGuqSsrKyo5KEEtEseBg0ahMuXL4uOJDXN+XDFihUxefJk/Pzzz6IjSS85ORlVqlTBsmXLEBsr50JZGugSc3FxQWhoKJU8KKRd8kBLvHTTLnno2LGj9OfDopmbmyMxMRGMMfTu3VvKJV400CUXFhYGZ2dnnD9/HkFBQaLjSE+75KFr166i40jP29sbn376KR4+fEhLvBRo1aoVvv76a+Tk5KBt27bSPQhPA90AJCQkoFq1ali/fj2io6NFx5FedHQ06tati3379mHhwoWi40hv4cKFaNmyJX799Vd8+umnouNIb9y4cejcuTP++OMP9O3bV3Scf6GBbgBMTU3/tcSLSh5eTnM+bGpqis8++4xKHhRISUmBpaUlFi9ebPAlD2Vh165dsLW1xbZt27BmzRrRcYrRQDcQzZo1+1fJg2y/6smGSh5KxtLS0mhKHsqCWq3G4cOHoVarMWLECJw9e1Z0JAA00A3K8OHD0a1bN9y6dQu+vr6i40jP19cXo0ePRlZWFlxdXUXHkV6bNm0wd+5cWuKlUP369bF+/Xrk5+fDxcVFiiVeNNANTExMDOzs7LB7924qeVBg6dKleOedd/DLL79QyYMCEydOhIeHh0GXPJSlAQMGIDAwEHfu3IGnp6foODTQDY32+TCVPCiTnp4OCwsLKnlQaM+ePahevbrBljyUtcjISLz55ptITk7GzJkzhWahgW6AqOShZCwtLXHw4EEqeVBIrVYjIyMDJiYmGDp0KC3x0kGlUiEjIwNmZmYIDQ1FWlqauCzCrpm8Fn9/fwQHB1PJg0LaJQ+Ojo50PqxDw4YNsW7dOuTn58PJyYlepKWDjY0N4uLiwDmHl5cXHjx4ICQHDXQDFhERgWbNmlHJg0KakodLly7hww8/FB1HeoGBgejfvz/+/vtvdO7cWXQc6Xl4eGDixIl49OgRnJ2dhWSggW7gNOfDVPKgjKbkYdOmTVTyoMD333+PBg0aIDExEXPmzBEdR3pfffUVHBwccPr0aXzyySdlfv000A0clTyUzLMlDxcvXhQdSWraS7ymTZuGI0eOiI4kvZSUFFSpUgURERHYuXNnmV43DXQj4OzsjLCwMGRnZ9P5sAKNGjXC6tWrkZ+fD0dHRzof1sHGxgY7d+4E5xweHh60xEsHMzMzHDp0CCqVCn369MGNGzfK7LppoBuJ0NBQtG/fHhcuXKAlXgoEBQWhb9++yMzMhLe3t+g40vPy8sK4ceOEng8bkpYtW2LBggXIyckp0xdp0UA3IgcOHIC1tTWioqKo5EGBTZs2wd7eHgcOHMD8+fNFx5HeggUL0KpVK5w8eVLqkgdZhISEoEuXLrh69Sp69+5dJtdJA92ImJqaIjU1lUoeFFKpVDh69CgqVqyISZMm0RIvBVJSUlC5cmUsW7YMcXFxouNIb9euXahRowa2b9+OVatWlfr10UA3Mk2aNCle4uXo6Cg6jvRsbW2xffv24iVe5OXMzc2Ll3j17NmTlnjpoFKpipd4jRw5stSXeCka6IyxLoyxc4yxC4yxyc/5fEXGWHTR548wxurpOyhRbtiwYejRowf+/PNPPHr0SHQc6fn4+CAkJARZWVnIysoSHUd6Dg4OmD9/fvH5MHk5e3t7REVFFS/x4pyX2nXpHOiMMRMAywF4A2gKIIAx1vSZiw0GcJdz/iaARQC+0ndQUjJbt25FnTp1kJeXh6dPn4qOI71Fixbh3XffRUFBAZ48eSI6jvTGjx8PT09PXLlyhe40KBAQEIBBgwbhzp07pbp/Sa3gMu8DuMA5vwQAjLEfAPgD+E3rMv4Awore3gbgG8YY46Xwv6ILFy7g4cOHcHd31/eXNjp2dna4evUqsrOz0bp1a1SuXFl0JKlZWloCAHJyctCiRQtYW1sLTiS3goICVKhQAbm5ubh//z79m1SgUqVKePDgAY4fP14qX1/JkUttAFe13r9W9LHnXoZzngfgPoA3nv1CjLFhjLFjjLFjmZmZrxS4evXqxf/wyMuZmprCxMQEAMAYE5xGfmq1uvj7pVLRw0u6qFQqvPfeewBQ/H0jL9eyZUuoVCpUrVq1VL6+knvoesM5XwVgFQA4ODi80r33EydO6DUTIYQYCyV3Q64DqKP1vl3Rx557GcaYGkBVALf1EZAQQogySgb6zwAaMcbqM8ZMAfQDEPvMZWIBfFT0di8ACaVxfk4IIeTFdB65cM7zGGOjAOwFYAJgHef8DGMsHMAxznksgLUAvmeMXQBwB4VDnxBCSBlSdIbOOY8HEP/Mx0K13s4GUDavbSWEEPJc9FA+IYQYCRrohBBiJGigE0KIkaCBTgghRoKJenYhYywTwJVX/M9tAJS3rjW6zeUD3eby4XVusz3nvPrzPiFsoL8OxtgxzrmD6BxliW5z+UC3uXwordtMRy6EEGIkaKATQoiRMNSBXvpdTvKh21w+0G0uH0rlNhvkGTohhJD/y1DvoRNCCHkGDXRCCDESUg/08lhOreA2j2OM/cYYO8kYO8gYsxeRU5903Waty/VkjHHGmME/xU3JbWaM9Sn6WZ9hjG0q64z6puDvdl3GWCJj7ETR328fETn1hTG2jjH2F2Ps9As+zxhjS4u+HycZY61e+0o551L+QeGq3osAGgAwBfBfAE2fucxIACuL3u4HIFp07jK4zR0AmBe9HVwebnPR5SoDSAZwGICD6Nxl8HNuBOAEgGpF79uKzl0Gt3kVgOCit5sC+J/o3K95m10BtAJw+gWf9wGwGwAD0A7Akde9TpnvoReXU3POcwBoyqm1+QP4rujtbQA8mGGXZ+q8zZzzRM7546J3D6OwQcqQKfk5A8CXAL4CkF2W4UqJkts8FMByzvldAOCc/1XGGfVNyW3mAKoUvV0VwI0yzKd3nPNkFPZDvIg/gChe6DAAK8ZYzde5TpkHut7KqQ2IktusbTAK/w9vyHTe5qJfRetwzneVZbBSpOTn/BaAtxhjaYyxw4yxLmWWrnQouc1hAAYyxq6hsH9hdNlEE6ak/951KtOSaKI/jLGBABwAuInOUpoYYyoACwEMEhylrKlReOzijsLfwpIZY+9wzu8JTVW6AgCs55wvYIw5orAFrTnnvEB0MEMh8z308lhOreQ2gzHWCcA0AH6c86dllK206LrNlQE0B5DEGPsfCs8aYw38gVElP+drAGI557mc88sA/h8KB7yhUnKbBwPYAgCc8wwAZihcYmWsFP17LwmZB3p5LKfWeZsZY+8B+BaFw9zQz1UBHbeZc36fc27DOa/HOa+HwscN/Djnx8TE1Qslf7d/ROG9czDGbFB4BHOpLEPqmZLb/AcADwBgjDVB4UDPLNOUZSsWQGDRs13aAbjPOb/5Wl9R9CPBOh4l9kHhPZOLAKYVfSwchf+ggcIf+FYAFwAcBdBAdOYyuM0HAPwJ4NeiP7GiM5f2bX7mskkw8Ge5KPw5MxQeNf0G4BSAfqIzl8FtbgogDYXPgPkVQGfRmV/z9m4GcBNALgp/4xoMYASAEVo/4+VF349T+vh7TS/9J4QQIyHzkQshhJASoIFOCCFGggY6IYQYCRrohBBiJGigE0KIkaCBTgghRoIGOiGEGIn/DzYaPXOuXsn+AAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[[0.25 0.5 ]\n",
            " [0.5  0.5 ]\n",
            " [0.75 0.5 ]\n",
            " [0.25 0.  ]\n",
            " [0.25 1.  ]\n",
            " [0.5  0.  ]\n",
            " [0.5  1.  ]\n",
            " [0.75 0.  ]\n",
            " [0.75 1.  ]\n",
            " [1.   1.  ]\n",
            " [1.   0.5 ]\n",
            " [1.   0.  ]\n",
            " [0.   1.  ]\n",
            " [0.   0.5 ]\n",
            " [0.   0.  ]]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "b_NIHQIdK6Kn"
      },
      "source": [
        "### Assembling of stiffness matrix and right-hand side"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lDvEVoQP9ThD"
      },
      "source": [
        "Let $\\phi_i$ be the hat-function associated to the $i$-th node in the mesh $\\mathcal{T}_h$. We denote with $\\mathcal{E}_h^N$ the set of boundary edges belonging to the Neumann boundary. The linear system associated to the FEM approximation reads $M \\boldsymbol{u} = \\boldsymbol{F}$, where\n",
        "\n",
        "*   The stiffness matrix $M$ is given by:\n",
        "$$\n",
        "M_{ij} = \\sum_{T\\in \\mathcal{T}_h}\\int_T \\boldsymbol{\\nabla} \\phi_i \\cdot \\boldsymbol{\\nabla} \\phi_j\n",
        "$$\n",
        "*   the right-hand side vector $\\boldsymbol{F}$ is given by:\n",
        "$$\n",
        "F_i = \\sum_{T\\in \\mathcal{T}_h}\\int_T f\\phi_i +\n",
        "\\sum_{E\\in \\mathcal{E}_h^N}\\int_E g_N\\phi_i - \n",
        "\\sum_{P_j\\in \\mathcal{V}_h^D} g_D(P_j)\\sum_{T\\in \\mathcal{T}_h}\\int_T \\boldsymbol{\\nabla} \\phi_i \\cdot \\boldsymbol{\\nabla} \\phi_j.\n",
        "$$\n",
        "\n",
        "In **problem 1** and **problem 2** the contribution to $\\boldsymbol{F}$ associated to the Neumann BC vanishes. \n",
        "Moreover, we can deal with non-homogenoeous Dirichlet conditions using the usual trick. Thus, we start by focusing on the assembling of $M$ and the term corresponding to the volumetric load function. We will explain later how to take into account the BCs."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ERG2lYz2JLec"
      },
      "source": [
        "Let the nodes $P_i$ and $P_j$ belong to the boundary of element $T$ and let $e_i$ and $e_j$ be the edges of $T$ opposed to the vertices $P_i$ and $P_j$, respectively. We observe that \n",
        "$$\n",
        "(\\boldsymbol{\\nabla}\\phi_i)_{|T} = \\frac{R_{\\pi/2} \\boldsymbol{e}_i}{2\\ |T|}, \\quad\\text{ and }\\;\\;\n",
        "(\\boldsymbol{\\nabla}\\phi_i\\cdot\\boldsymbol{\\nabla}\\phi_j)_{|T} =\n",
        "\\frac{\\boldsymbol{e}_i \\cdot \\boldsymbol{e}_j}{4\\ |T|^2}.\n",
        "$$\n",
        "Thus, $\\int_T \\boldsymbol{\\nabla}\\phi_i\\cdot\\boldsymbol{\\nabla}\\phi_j = \\frac{\\boldsymbol{e}_i \\cdot \\boldsymbol{e}_j}{4\\ |T|} $.\n",
        " "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iLAnWmyOL0PW"
      },
      "source": [
        "The term associated to the load function $\\boldsymbol{f}$ is computing using a simple quadrature rule (evaluation at the center of each element)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "j74S1MzRSOHR"
      },
      "source": [
        "def assemble(P, T, f):\n",
        "  '''\n",
        "  construct the Grad*Grad stiffness matrix and rhs vector \n",
        "  corresponding to mesh (P,T) and voulmetric load function f\n",
        "  '''\n",
        "  n_nodes = P[:,0].size\n",
        "  n_elems = T[:,0].size\n",
        "  M = np.zeros((n_nodes, n_nodes))\n",
        "  b = np.zeros(n_nodes)\n",
        "\n",
        "  for k in range(n_elems):\n",
        "    # compute edges of element k\n",
        "    e = np.zeros((3,2))\n",
        "    e[2,:] = P[T[k,1],:] - P[T[k,0],:]\n",
        "    e[0,:] = P[T[k,2],:] - P[T[k,1],:]\n",
        "    e[1,:] = P[T[k,0],:] - P[T[k,2],:] \n",
        "    area_T = np.linalg.det(e[0:2,:])/2.0\n",
        "    bar_T = np.array((P[T[k,0],:] + P[T[k,1],:] + P[T[k,2],:])/3.0)\n",
        " \n",
        "    # each local contribution is stored in the corresponding global position\n",
        "    for iloc in range(3):\n",
        "      b[T[k,iloc]] += area_T * f(bar_T[0], bar_T[1]) / 3.0;\n",
        "      for jloc in range(3):\n",
        "        M[T[k,iloc], T[k,jloc]] += np.dot(e[iloc,:], e[jloc,:])/(4*area_T)\n",
        "\n",
        "  return M, b"
      ],
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lminviDUjWvF"
      },
      "source": [
        "### Main program"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hRM7sQSaTSpv",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 529
        },
        "outputId": "91964bc4-c812-46ef-cb8b-2c119b02411f"
      },
      "source": [
        "pb = 2\n",
        "n = 20\n",
        "f = lambda x, y: 4 - 2* (x**2+y**2)\n",
        "u = lambda x, y: (1-x**2)*(1-y**2)\n",
        "def grad_u(x,y): \n",
        "  grad = np.array([-2*x*(1-y**2), -2*y*(1-x**2)])\n",
        "  return grad\n",
        "\n",
        "# Step 1: meshing\n",
        "P, T, n_diric, n_neumn = ordered_mesh(1, 1, n, n, pb)\n",
        "\n",
        "# Step 2: definitions of discrete space (dofs)\n",
        "n_nodes = P[:,0].size\n",
        "n_dofs = n_nodes - n_diric\n",
        "uh = np.zeros(n_nodes)\n",
        "\n",
        "# Step 3: assembling matrices and rhs vectors\n",
        "M, b = assemble(P, T, f)\n",
        "# accounting for Neumann BCs\n",
        "for j in range(n_dofs-n_neumn, n_dofs):\n",
        "  b[j] += neumann_rhs(1, 1, P[j,:], grad_u, 1/(2*n))\n",
        "# accounting for Dirichlet BCs\n",
        "uh[n_dofs:n_nodes] = u(P[n_dofs:n_nodes,0], P[n_dofs:n_nodes,1])\n",
        "b[0:n_dofs] += -np.dot(M[0:n_dofs, n_dofs:n_nodes], uh[n_dofs:n_nodes])\n",
        "\n",
        "# Step 4: solving linear system\n",
        "A = M[0:n_dofs, 0:n_dofs]\n",
        "uh[0:n_dofs] = np.linalg.solve(A, b[0:n_dofs])\n",
        "\n",
        "# Step 5a: compute errors\n",
        "err_h1 = 0\n",
        "err_l2 = 0\n",
        "for k in range(T[:,0].size):\n",
        "  # compute edges of element k\n",
        "  e = np.zeros((3,2))\n",
        "  e[2,:] = P[T[k,1],:] - P[T[k,0],:]\n",
        "  e[0,:] = P[T[k,2],:] - P[T[k,1],:]\n",
        "  e[1,:] = P[T[k,0],:] - P[T[k,2],:] \n",
        "  area_T = np.linalg.det(e[0:2,:])/2.0\n",
        "  bar_T = np.array((P[T[k,0],:] + P[T[k,1],:] + P[T[k,2],:])/3.0)\n",
        "\n",
        "  gradu_bar = grad_u(bar_T[0], bar_T[1])\n",
        "  grad_uh = np.array([sum(-e[:,1]*uh[T[k,:]]), sum(e[:,0]*uh[T[k,:]])])/(2*area_T)\n",
        "  err_h1 += area_T * np.dot(gradu_bar - grad_uh, gradu_bar - grad_uh)\n",
        "\n",
        "  u_bar = u(bar_T[0], bar_T[1])\n",
        "  uh_bar = (uh[T[k,0]] + uh[T[k,1]] + uh[T[k,2]])/3.0\n",
        "  err_l2 += area_T * (uh_bar - u_bar)**2\n",
        "\n",
        "err_h1 = np.sqrt(err_h1)\n",
        "err_l2 = np.sqrt(err_l2)\n",
        "print('n={} eL2={:.2e} eH1={:.2e}'.format(n, err_l2, err_h1))\n",
        "\n",
        "# Step 5b: plot the results and compare with true solution\n",
        "Q, T = left_mesh(1,1,4*n,4*n)\n",
        "u_ex = u(Q[:,0], Q[:,1])\n",
        "fig, axes = plt.subplots(1,2, figsize=(12,6))\n",
        "\n",
        "axes[0].tricontourf(P[:,0], P[:,1], uh , 20)\n",
        "axes[0].set_title(\"FEM approximation\");\n",
        "\n",
        "axes[1].tricontourf(Q[:,0], Q[:,1], u_ex , 20)\n",
        "axes[1].set_title(\"exact solution\");\n",
        "\n",
        "fig.tight_layout()\n",
        "plt.show()"
      ],
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:7: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.\n",
            "Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations\n",
            "  import sys\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "n=20 eL2=9.43e-04 eH1=4.29e-02\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 864x432 with 2 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAGoCAYAAABbkkSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5Rl51nf+d9TfZEbVV+qJTwGSUj2RJ6gcTJchIGESZgYiDBgkRXGsQnjGAgeIIYMeMiYS4AxMBPCCjBOHECOHWMSbAOZkAbEMouLFwvH9lgJxMEmBiFsJIEjdFfL7b6+88c55T5VdS57n/Pu/T7P+34/a/Wyuqu6ztZxq779vPvde1tKSQAAAACAzW2VPgAAAAAAqAUDFgAAAABkwoAFAAAAAJkwYAEAAABAJgxYAAAAAJAJAxYAAAAAZMKABQRiZmfN7DkjvdaPm9k/HOO1AAAxmdnLzey3Nvj9tAbVYcCCS2b2ITM7Nx0odn98spndYmZp36+fNbO/Nf19b5p+/M59X+9Hpr/+8iL/QpmklLZTSvfl/rrzAplS+vqU0vflfi0AwPimffz+wsdAa9CEw6UPAFjiy1JKvzr7C2Z2y/QfT6WULi34fb8v6WWS/t309xyW9GJJfzjMYXZnZibJUkpXSh8LAAAA8uMMFmr0C5I+z8x2pj+/Q9L7JH1k0W8ws+eb2bvM7HEz+1Mz+2dmdnTm48nMvtnM7jOzh83sh8xsa/qxl5vZO6e/5wkz+y9m9oKZ3/sOM/sBM3unpI9Keo6Z/SUze+/0899rZn9p+rl/y8z+yMxOTH/+xWb2ETP7xJnj+HPTf36Tmf1zM/vl6Vm8d5rZs8zsR83sselxfPrMcbzazP7QzJ4ysw+Y2d+Y/vqnSvpxSZ87/TqPz3z975/5/V9nZvea2aNmdsbMPnnf+/P1ZvYH0/fwddNhEgAwx3RXxr8xsz+bft//5umvnzazB8zsy6Y/355+733Z9OdfYma/bWZPmtn9Zva9+77u55nZv59+L75/2qhXSPrbkv7B9Pv8L8w5Hpvu9nho+rX/s5k9b/qxk2b25umxftjMvmu3gfu+xu4uk8Mzv/YOM/u7tAYtYcBCjT6mydmrl0x//jJJb17xey5L+hZJ10v6XEkvkPSN+z7nb0i6XdJnSLpT0tfMfOyzNTlDdr2k75H0/5rZ6ZmP/y+SXiHpuKSnJP2SpNdKuk7SD0v6JTO7LqX0Nkn/XtJrzew6SW+Q9HdTSn+24LhfLOm7pq97XtK7JP3H6c9/bvq1d/2hpP9R0klJ/6ekf2Vmn5RS+j1JXy/pXdMtiKf2v4iZ/TVJ//f09T5J0oclvXXfp32ppM+S9Benn/fXFxwzADRtOpz8gqT/JOkGTZrzv5nZX08pPapJX15vZs+U9COSfieltNuxpzXp2ilJXyLpG8zsy6df92ZJvyzpn0r6REmfNv29d0n615L+8fT7/JfNOawvkvRXJD1Xk068WNIj04/90+mvPUfSX52+/lf3+XemNWgJAxY8+/npCtXjZvbz+z728MzHHp+ujM16s6SXmdkpTWKw//fvkVL6Dymld6eULqWUPiTpJ6a/b9YPppQeTSn9saQflfTSmY89JOlHU0oXp0PSBzUJ3643pZTeP93W+EWS/iCl9FPT13uLpP8iaTd4f0/SX5P0Dkm/kFL6xSWH/m+nx/4xSf9W0sdSSm9OKV2W9DZJHz+DlVL62ZTSn6SUrkyP8Q8kPX/Z+zLjb0t6Y0rpP6aUzkv6dk1WIW+Z+Zx/lFJ6fPr+/IYmYQcAHPRZkj4xpfSalNKF6bW1r9d0YTCl9CuSflbSr0l6oaT/dfc3ppTekVL6z9Pv5e+T9BZd7dVXSvrVlNJbpj16JKX0Ox2P6aImi4B/XpOt7L+XUvpTMzs0Pa5vTyk9NW3kP9Fk4TA3WoMqMGDBsy9PKZ2a/vjyfR+7fuZjp6YrYx+XUvotTVbvvlPSL6aUzi17ITN7rpn94nQ73pOS/i9NzgLNun/mnz8s6ZNnfv5gSikt+fjs7/3k6cdnfViTVUyllB7XJKzP0yRiy/zXmX8+N+fn27s/MbOXmdnv7A6l06+//99xkT3HnFI6q8nK5g0znzO7BfOjs68NANjjZkmfPLtQKOk7JP03M59zlybfp9+UUto9kyQz+2wz+43pdr0nNDkrtPu9/Cateb1xSunXJf0zSa+T9JCZ3WWT7erXSzqivd36eLMyozWoAgMWavavJL1Kq7cHStKPaXIW6daU0glNQrd/X/dNM//8KZL+ZObnN+zbB77/47PD159oEtdZnyLpQUkys0/TZHvIWzTZRrix6baR10t6paTrplszfldX/x3Tot8775jN7FpNtjc+mOP4AKAx90v6o30LhcdTSi+UpOlZo7s06dc32vTa26mflnRG0k0ppZOaXNdkM1/3v13wmqu+zyul9NqU0mdKuk2TrYLfJulhTc5uzXbr483a5+np/37CzK89q8cx0BpUgQELNXutpC+U9JsdPve4pCclnTWzPy/pG+Z8zreZ2Y6Z3STp72uyBW/XMyV9s5kdMbP/WdKnSrp7wWvdLem5ZvaVZnbYJreYv03SL5rZMzQZDL9Dk/3tN5jZ/mvB1nGtJmH7M0kys6/WZGV013+VdKPN3Nhjn7dI+moz+zQzu0aTM3zvmW4VAQD08/9JesrM/g8zO2Zmh8zseWb2WdOPf4cm37O/RtIPSXrzdOiSJr16NKX0MTN7vibbAnf9a0lfYGYvnvbluuminTT5Pr/wOYpm9lnTs2NHNBmUPibpynTL+c9I+gEzOz5dsPtWTVq1x/R64QclfdX03+lrtHfgozVoAgMWonrc9j4H61v3f8L0eqlf27d1b5H/XZNIPaXJmZ63zfmcfyfpP0j6HU1uUvGGmY+9R9Ktmqz0/YCkr5jd0rHvuB7R5CLdV2my9eEfSPrSlNLDmlzce39K6cem+8+/StL3m9mtHf4dFkopfUCT7Ybv0iRwf0HSO2c+5dclvV/SR8zs4Tm//1cl/UNJ/0bSn2oSzJfs/zwAwGrToeVLNbl+6I80ace/kHTSzD5TkwHmZdPP+0FNhq1XT3/7N0p6jZk9Jem7NRl+dr/uH2tyzdarJD2qSa/+h+mH3yDptgXXNUvSCU3695gm2/Qe0WS4k6Rv0mTouk/Sb2lyFu2NC/71vk6TM1+PSPrvNblx0y5agyZYt797Am0zs6TJ9sF753zs5Zrc6e/zRj8wAAAAuMIZLAAAAADIZOWAZWZvtMlD5353wcfNzF47fSjc+8zsM/IfJgAAi9EqAIAXXc5gvUnSHUs+/sWaXHtyqyYPUv2xzQ8L8CWlZPO2B04/9ia2BwLFvUm0CgDgwMoBK6X0m5pcKLnInZLenCbeLemUmX1SrgMEAGAVWgUA8OJwhq9xg/Y+RPWB6a/96f5PNLNXaLJyqEM6/JnXbp3I8PIAgGievPLowymlTxzxJWkVAKCXdVuVY8DqLKV0lyYPztPJQ9elz3nGl4z58gAAJ37loz/14dLHsAitAgBI67cqx10EH5R008zPbxRP3AYA+EKrAACjyDFgnZH0sukdmj5H0hMppQNbLgAAKIhWAQBGsXKLoJm9RdLnS7rezB6Q9D2SjkhSSunHJd2tyVPD75X0UUlfPdTBAgAwD60CAHixcsBKKb10xceTpL+X7YgAAOiJVgEAvMixRRAAAAAAIAYsAAAAAMiGAQsAAAAAMhn1OVgAAIRy6LC2Tu8M+hJXHn1s0K8PABhXuQFrhGihDP6yAADd0cJy6BWAIXAGC9lF/8sCwQWANnjvFT0CYmLAAvbJFVzCCADYxNgDIN0C8mDAAgYyVhgJIgAgh9zdok9oFQMWEFyOIBJBAEBu9AmtYsAC0DmChA4AMKZ1hjRahdIYsAB0RugAAN71aRWNwhAYsAAMqkvoCBwAoAR2cGAIDFgAilsVOMIGACiJTqGPcgPW4S2lnePFXn4d9thTpQ8BaBJhQzEBWkWbgPKWdYpGtYczWD14j2wfBBk1WRQ2ooYW0CbANxrVHgasRpUIMuHE2IgaEEupYZE+oYR5jaJPdWDAwmiGCidhRF+zUSNmAIboE23COhi66sCAhfD6hpHoYdZuzAgYgJxyDW00Cwxd8TBgoTmbRo/Y1Wl/wIgXAA/WbRatqhvN8o0BC+iJM2ZtYMUQQGS0qi0MXL4wYAEDY/WxHgQMQK36tIo++ce1xmUxYAFOrYodgSuPs1wAWtRlGKNRfjBsjY8BCwhqUeCIWlmc5QIAhjCvGLbGUWzASltbunziWOfPP/TkuQGPBqgHg5cvxCy2vq1aB31Dq5YNYTRrePRpOGHOYA0dOA+ILIbUdX89URsOt4THPJH7RrcwFIavcdGnvMIMWC0oFVkCiVlc+zW8rdM7RAxVGKtbdAqz2KkxHM5q5cGAheyBJIR1mxc2otYfq4VAdzk7RaPqxeCVF51aHwMWsls3hEQvrv1RI2bdETBgXJsMa3QqJhYGN8Oui/4YsOAGg1k9WEXsj4AB/q3TKRrlE0NXPywG9sOAhfC4G2UcnOlajoAB9aFRcdCo1ehUNwxYaEqX0BG48RCz+QgY0KauwxidGgdnuRajU8sxYAH7ELhyZmNGxNg2CGA+FgvLoVN70an5GLCANa0KHHHbDGe3JogXgHXQqOHRqQnOZh3EgAUMZFncCFt/La8aMmQByI0BLL/WBy5adVWxASsdMl04eaTUy2d19ImLpQ8BwcwLGzHrrsWIEa4yamkVnUJfiwYwWtVdiwuDnM2a4AxWBmPFl0DWjZitr5WIMWRhXXQKubA7Yz2tdGpX671iwAokZyCJYBzErJ/aI9Z6tOAbnWobuzO62e1UjY2a1XKvGLAatU4EiZ0/nPVartZhq+VooR2bDGv0yg86tVitjZrVaq8YsNBZ39gRuHJYRTyothXDVqMFdEGv/NvfKRpV77DVYq8YsDAYAucLQ9dETYNWi9EChtCnV7RqGDTqqpo6tau1XjFgwY1VgSNq+c0GrbWQ1RKw1qIFlNZlGKNXebQ+dNXSqV0t9YoBC2EQtWG1ul2jhm0ZLUULiIBeDafFVtU0aLXSKwYsVIWzYPm0eHYrcsRaiRZQC3qVR0utityoWS30igELTVkWNGK2WEsBk+JGrIVoAa2gV/21cnYraqNm1d4rBixgitXEbloJmBQzYrVHa2zpsOnCyUMLP370icsjHg0wsahXdGqv2nsVsVGtKDZgrYpWJAS2DawmztfC2a20czxUwBiyxhOtY/Sqbgxey9Xaq6iDVs2t4gxWBrkDSwDjmRe1FoNW8x2fogWs5nBhfUMMhDTLPwavg2octqJ1Sqq3VQxYDuUKINEri6FrYjdiNQUsSrxqDRd82aRZdKosOjVR27AVqVNSna1iwKpY3+gRuuG1HLOaAhZplbDGcKEe6wxntGpYLXdKqqdVkTpVIwYsfFyf0BG4fFqMWU0BI17AuGjV+FrslFTHDowonaptMZABC2vpGjjitp79Mas5ZNGHrQirhLWFC+iqS6vo1HroVBwROiXV1SoGLAxqVdwIWzethCzyaqH3VcKawgXkxBCWRytnuegUumDAQlE8X2Y9tQ9cUVcLiRdQJ1q1ntlW0SkfvJ/NqmUxkAELbhG07moeuKKtFnoesmoJF+AJreqGTvniuVU1YMBCSARtuRpDFilgnsPFkAWMh1YtRqfK89qqGjrFgIXqLApayzGraZvG5RPHQsTLa7gA+ECr9qqtU1KMQctrq6IPWcUGrCtb0vnjW6VevrNrnrpS+hCQCTGbqGHVMEq8CFd8UVo1i27FNq9VLXcqYqN20ap2cQZrhSHCSvx8aX3w2g1ZxIhFOJtFuDC2nN2iVz603Kkahq0Ig5b3m19Ew4BVwCbxI3bj2R+02kMWNWJRwuUtWpzFQhfr9opWjaO1s11RO7WLXvUTuVMMWMF0jR1xy6+lgSviWS3v4fIUrV2R4wXf+g5mNCufVloVediiV/VjwKoUcRteCxGLGDDP2waJFjAfi4fDaeEsV8RWSfSqi6gLgQxYkNQtboRtudoHrkhntTyvDnqJ1q6o8UKbaFUeNfcqUqskelUrBix0tipsRG2vWlcNI8XLa7iIFjAcWtVfjQNXtLNa9GqxiAuBDFjIZlHUiNlVNUUsUrw8hstDtIAWLRvA6NVETa2S4i0MemoV1sOAhcERs8VmIxY5YFHi5W3Q8jJkRVwdBIZAr+arrVWS717RqoOidYoBC0URs6tqWDEkXv15CBeA1dilMVFDq6QYC4OezmbRqn4YsOBW6zGLvmIYYdjyEi/CBcRFq+poFZ1ajVZ112nAMrM7JP0/kg5J+hcppX+07+OfIuknJZ2afs6rU0p3L/ua6ZB0cXutY+7tyNlxXgfjmBez2kNGwIbj5WwW4dpc9FZJ9KomtCpWq+hUNyVbFWmb4MoBy8wOSXqdpC+U9ICk95rZmZTSB2Y+7bsk/UxK6cfM7DZJd0u6ZYDjXcvQcSSI5bUUssjbM7wHrHS8GLLWV0OrpGF7RavKa7VVdCoPD52SaFUXXc5gPV/SvSml+yTJzN4q6U5Js9FKkk5M//mkpD/JeZDebRJEgjecVkIWMWIXTh4hXguUCleklcEFaNUKtMqnFloVtVOSv0HLQ6ewWpcB6wZJ98/8/AFJn73vc75X0q+Y2TdJulbSF8z7Qmb2CkmvkKQjx3f6HmuV+gaPyG1mf8iIWDnEC5nRqgH1aRWd2lzNrYrUKclnqzx0isXA5VY/Er2bl0p6U0rpRkkvlPRTZnbga6eU7kop3Z5Suv3QJ1yb6aXbcnF79Q90d/741oEftbhw8tDchx17c+HkkT03xPBgd797KWnneNHXrxitGkGXTtGqfmrvFK3q7/KJY7TKsS5nsB6UdNPMz2+c/tqsr5V0hySllN5lZs+QdL2kh3IcJPrpEi5WGBerbeUwymqht1XC0iuE7HHvjVYFQ6vWV1unpKut8twpyd8W99KtwnxdlkHeK+lWM3u2mR2V9BJJZ/Z9zh9LeoEkmdmnSnqGpD/LeaDIi1XF7mpaOYywUuhplbD06uDYtk6H3g5HqyrEWbBuauyU51Z56pRUtlUlzmJFaNXKM1gppUtm9kpJb9fktrZvTCm938xeI+melNIZSa+S9Hoz+xZNLiJ+eUopDXngGM6yaLGauHflMOqqYYSzWl5WCUuuDnIWqzta1SZ6NV8tZ7i8n9XytPOCVvnS6TlY0+eE3L3v17575p8/IOkv5z00eETM9qpp2PIYMC/xYgtGDLQKsxb1ilbFa5XnTkksCOKgTgMW0EXrMSNgw/EQr1LhYmUQyItWxT275Xn3ResLgmO3yvvdBBmwMLgWY1ZDwIjXQS2sDnqPFjCUea2quVO7oi4Oem4VQxYYsFBMSzGLGDDiNV+JcBEtoIzWFghpVR4sCKLYgJW2hrsDUK3f+FrQwtAVLWDE6yCGrHYM2ap5avt+Vyta5Y/XVrU2ZNGqiSrPYA0Zw9q+gUZQc8giBezCyUOuwiWVjRergxhC7n7V8r0ygv3/39X03kdrleRn0GpxyBqL5y3tVQ5YQ1o3fjV9o/WgxpDtBsxzvLyFSyp7NmvscLEyiL7WaVYN3089qHVxMMqw5alXre26oFUMWKPpGrkavvmWUNPAFSFensK1q9QqYY2rg55XBTG8PkNZ5O+1JdTUKomFwb7YddEOBixnCFsetUTMe7w8hUtqY8hiZRCedGlW1O+/Y6itVZLvXnloVektg2NpvVUMWIGtClvUb9RDmH2vIr4v3uPladBqYcgCIqFV3UVvleR7YdBLq+hUPl53XDBgVWxZ1KJ+484h+oqh52GLeI0TrzFWBr1GC/WhVfPRquF4aFWp67LYcTEOBqxGtfZsj2UiR8xrwLzEiyELiI1WXVVDqzx1SvKxbbBEq2o8k+UNAxb2YCUx7hYNjwErPWi1std9KJzFglcMXjFbRacWHUO9Q1ari4Fbqz8FmLi4Pf9HzSL+u54/vrXnzJYHuwEr89pHRn/NyyeOjfI6aef4KK8DRNFip6R4/650at7rj9+qWmyd3il9CAf4+tONkFoKWaR/R28Ba23IAuBHS0NXpH9Hj51qqVUsBg7Hz59qVKOFkEX6d/MUMMIVj8eVQSAHWuXHbqdoFa2qBddgYRS1PtFeirMP/vzxLRf73kvudx97n/sYe9xb3d8+mi3p0nZa+WmHz9oIB4Oh0aryaBWtqkG5AatjtMZAGMuIfEekRXb/nbz+u3i6wLjU3ZtqDBfKG7NnNGtctGp8tIqbNEXHGSxtFkZCl09NEfO+UuglXoQL6G/dZtGrPGjVeDy1Shr/bNaYrYq+GOjtrrcMWBvqEzri1k8tEfMcMA/xIlyba23rBdZDr4ZBq4bnoVVSmUXBmhYEW2oVA9aIusaNsM1XQ8S8bsvwEK/awxV9dRBt4bqz9dGq4dCqYdGpfBiwHFoVNqI2ETliXlcKS8er5nBJw8ZryJVBb1sv4AOt6sbr9/suGLTmq7lV7LjIgwEroGVRazloUSPmMWAl41VzuICW0KqDoi4Meu1ryTsOltjeTqviYMCqzKKgtRYzrzFYxuMxl4pXzeFiCwZAq3Z5/L6/irdFwdbOZo3RKjq1OR9PdMPgLm2nuT9aEOVBi7M8HW/Jhz+O/bDHsR/wmFvaOV76EICN0Co/3/tX8XasJR9SXPLBxEMZ8gHEQ7Vq6/TOIF93HZzBaty8cNW8ghhttdDTSmGpVUJWBwHQKr/o1MSYrWKroH+cwcIBraweRlot9HScJVYJOZMFYL9WznZFaZWnY2xh18UYnYp4FssLBix0UnvECFg/pYasmrZhDBWu2qMFrEKryvJyfKW2t9c2ZGE9DFhYW61Dl5c4LOPh+GqPF+Hay9PedqCP2lvloQfzeDm2mjs1ea1hWzXkWayacQ0Wstofrsh75L3fTtfLvvcSdxoca6/70Pvch7oWq5XnjIwhHZIubc//8334LGuUUdXUKslPD+bxcGwlrs0q8cgR+FFswFoWrTEQxnHUFDEPkZjHw3HVPGQBi4zZMJo1rFpa5fnmGC22isXA5YZYDNw6vaMrjz6W9Wuuo9kzWJuEkdCtr4aIeQ1Y6XjVukIYNVyoyzrNolXrq6lVnjollT8uhiyModkBaxN9QkfglosesdKhmKf0MdUYr4jhYpsguraKTq0226qonZL8taqVBcEahiwWA/thwBrYqsARtr2iRqz0UDNP6XiVeA5JVIQLJXUZxGjVVZEXBr21qvTxjNkqtrXPV+tiIANWYcvC1nrQIg5bpWOxX8njqS1cEc9i5eZlbzvGR6sWo1Wbo1U5X4OzWB60/V3RuUvbV+b+aFG02+t6uT3trlLHM+at3Me4Le6Qt8Md4la4PBMLY1jUqhZ7Ras2U+pYxryVe/RWoRvOYAXU+i2DI60Wetv7Xmrb4FgrhGzBAHxpuVe0aj2lzmbVdiZrKEOcxapxm2D93+Ea0uLqYaTVQi8rhSXPZo1h6NXBaGexAI9a6xWtWu84xlbTrouWz2Jtnd4pfQicwapdS6uHUVYLvex9L3E2q5YzWVyPBeQ3r1e0qhwPrar9bFbkM1lYrr7vXOik9pXDCKuFHlYJSxzDmHvdwXVYiK32M120qvsxjK2GVg11FovdFqvF/9ODLAhYOS3Ga4xwRd1+QbiA5WrtVYRWlX59WrXO14+xVbC2xUAGLMxV66qh93iVDBjhAhBRjb3y2qrSndo9hjHRqvlYDFyOAQud1RQwzyuFJQM29mtHD1eLZ7E8XDwMrFJLr7y2qvSgxZDV92vHOItVk3I3udhK0val1Z93lvtweLU/WlEvRt4Nl7eLjUvdUn3s1x7jYuIhLySOcMOLGm+BO5qurVqFlhU12ytalU/JG2GM/drRWzUEHjy8mP/v+DnCtgjByyr6wOXxzk6l41VTuICicrWMbm2sllZ56ZTEgmAELAaOq+3v1H2DR9h6iRwxbwErFa+awhXtLBYrg5irT7doVidRz255WxRkQTCPaK3CfHz37aNr2IjaXBEj5mnQKhWvMV+XcAEZsQ2/t6gLg7SKIasUr4uBW6d3dOXRx4q9Pt9Zh7AqagQt3LBFvMaLF+G6Kme4atp6gYyW9YpW0aoNlNh5wYJgl6/LYuAY+O5ZAkHbI1LAiBfhWv51CRcqQqv2oFX9sSCIvmpZDPT9HaJF25fm/2hElNvqerltbolb5Y71emPcFhfAmmgVreqhxC3da2jVULduz33bds+PFimFv8FE0WDMIgSs1XgRrmVfl3ChYY11SqJVXbEguJ6IDyEGA1Z8jYTMe8BajNdYr0e48kk7x0sfAlrU0NDluVOSn1bV/HpRsBg4LAas2jQQMs8BI17DYLsgUJnKW8Wi4HIlFgSHxmJgPjUsBvK3lhZUGjLPAfMQr5pfL6chwsXKILCGCjsl0aplGLK6i9AqXMWA1arKIka8DqptyIoWrpptnd4pfQhoQYWLg55bVQpDVj1YDLyq3T8FuKqigHldKWTIyqPlcOVSw9YLNKySXnntVMlWjdWryEMWZ7H6KbkYyN9WcFAF8ZL8BaxUvGrb6x4lXGwTBEYQvFceFwVLD1o1vQ7axYCF5SpYLSReEzUNWQBwQPBeeWxVCbUMWVEWA3PLtRgYfbdFsUexbx26omPb5wf52ufOXjPI14UORutssT9Cve2G6/BZH+sKl7aTDp+1UV9zrKfbD/1a549v6Zqn8v9F5MLJQzr6xOWMX++Ijj5xMdvXw/iGbNV+tCuz2V7RqrXsDlm1tmro16FV7YrzHaeHdWNI3NYQMGDe4lVruIZ+rSjhyunyiWM69OS50oeBgfRtF83qgVZtpFSrpOF7FXXIgm8xvsuMpGvciNoCwQLmJV4lVgjHCtfua4010HnkcWUw7RyXPfZU6cPAhro0i17NQavWPI56z2ZF7BRnsXzz/53FoVVRI2gKFTBP8aoxXENq8SwW0Ae9WoFWrXEcZRYFIw9ZrZ3FYrcFA9YgFgWt2ZAFCZiHeBGu/lobsggXcqJXM2hVz+MYd1GQVh3ktVNgwBoVIVOIgHmIF+Hqx/vqIFsvEE3zvaJVHY+BVvVBq/qJvJ3d50aIpM0AACAASURBVHeNxswLWRMRcx6w0vEa+2zWGNdlRd+SiKu2Tu/oyqOPlT4MjKzJwYtWrXj98VsVecjKzetZrNZ3W5S/NQ3mOrZ9fu6Pajl+bknpZ5OM/SySoZ8NMtTXH+KZIzmfN5L7wcObiv6MEfhCq8prqVVjPCsrUqty8taqqHz/v4wDqo+Y44dElowXQ1Y33sOVS64HOQJDqX6B0Hmryr12Gq1XDFlXeX/4cElbp3fKvG6XTzKzO8zsg2Z2r5m9esHnvNjMPmBm7zezn857mFik6og5jleJgI0ZLmmceEVQ81msmtAp/2jVeFo5m0Wn4NXKzcRmdkjS6yR9oaQHJL3XzM6klD4w8zm3Svp2SX85pfSYmT1zqAPGavvDFX5vvNP975e2rxTZ8z7mRcVD7kMf6mt7v4gY+dGpmKptlbNOSWWuzxqrVUNfLxWlVTmvxcp1s4sc12FFvdFFl//ini/p3pTSfSmlC5LeKunOfZ/zdZJel1J6TJJSSg/lPUxsoqpVQ2crhS2czRpyhZDtF8iETlWgmlY53D5YslVjiLqtHfXq8reQGyTdP/PzB6a/Nuu5kp5rZu80s3eb2R3zvpCZvcLM7jGzey49+dH1jhgbqyJiDuNV5nXjx6u1cLFNcBDZOiXRKi/Cd0qiVRpvQTDikOV5MZBWbSbX/7OHJd0q6fMlvVTS683s1P5PSindlVK6PaV0++ETn5DppbGp0BFzFC9WCH19bc/hyoEbXfTWqVMSrfIo/MIgrapiyBpCKzdnak2X/1cflHTTzM9vnP7arAcknUkpXUwp/ZGk39ckZAgmbMCcxWv812TImodwLVfRrdrpVGNo1eZqbVW0TuXGYqAPXa7EfK+kW83s2ZoE6yWSvnLf5/y8JiuC/9LMrtdkK8Z9y77ooa2kU9f2u/Dt8afb/D+plNlwhbn42MlFxiUuLK7lguIW5LqAGB83SKek9Vo1D/0aDq1aX62tinaDJq83Z6JV61v5X3ZK6ZKZvVLS2yUdkvTGlNL7zew1ku5JKZ2ZfuyLzOwDki5L+raU0iO5D3bTyBG49YULmKN41RauIXkPV867NCEfT51aZN1+0a1+aNV6amxVtCErJ1pVXqf/olNKd0u6e9+vfffMPydJ3zr94VafwBG1xUIFzEG8xl4hJFwHeV0dRD61dGq/Lt2iV/PRqn5oVVl0arGIt2r387AGZ1ZFjaBNhAmYk3gRrvJfe1O5VgZzbL3I8YwRxEevVqNV3dGqcl/X444LWrUeBqw1LQtaqzHbDRjxWoxwlcPqIFq1qFe0ynmrGhqyJA3aq0hDFurAgDWAeTFrKWTEa7kxt2FEHrI8h6um/e2bbL3YOr2jK48+lvmIMJbWBy/3Z7UcLAhK9WwZ9NyUWR7PYqE/BqyRtDh0Ea/lxlohjDxk5ebtLBZ3aIJHLfeKVh1U286L3KL0bxO0qj8GrIJaihjxmo8ha/yvCaC/VnrlemGwgVZNXmu4XkXZceFtMTC6ErsteAqnM6euPXfgR01cPxyy0MMfx3rY41gPI84t94Mdcz182MvDHFt9iCPKq7lVkmjVPpe2r1TRqwgPC86pllalneOZjmQcnMEKYH+4alk1dHlWq9AK4Vh73dnj7g9bL1CLWs9y0aq9ath5EWHHBWexYuMMVkC1rRq6XCXcvlRklXCM1cGhz2QNsTrIWSwgphpb5apX7LxYW4RW5ZKjVRdOHslwJO3gDFYFajnD5XLve4G7DY6xOhjxTBZnx4DYammV5OysVsGdF9HPZHnHWay4OINVoRpWDF2tEhY4m1XDmSzvajqLxXVYiIhWZVbprovJ6wzTq5bOYnnQUqsYsCoXfTthy/Eaa8giXOWV3noR7eJh1IdWZVLpguCQvLfK02Jg6VZFwoDVmKgBazVeY921qdUhK1e4AORFqzZUYMgaulXRrsdC2/jbRcMiBqzleA3/GnGGLG9YGQSGQ6vWVOHZrEhb21kMbBv/j0FSvIAVD5dUXbgiaWFo66ulve1oF61aQ2WtYjFwfR6uGd5EpO3sxe4ieGTrsp517VPZv+5Hno7z5nu1Gy7vd3hycRenke/eNPRdm4a6W5PnOwBylyYsM1SrJHq1KVrVw8h3xKVV+b+el1Zt+uzGyyeO6dCTMRZINlHdbdq7hpCwrTa7Qug5YK3Fi3D5dOHkIR194nLpw0AgfQY3mrUYreqIBUFgNNUNWF11CRtBuypCwFzEi3CNhpVBtGRVs+jVBK3qgFYt5fksVg4sBo6j2QGri0VBaz1k3rdlHNs+38QKIeHCfq1svcBB9OogWrVERUPWELy2isXAOBiw1kDIJjzHq5UVwtbD5e0sFiuD8IZe0aqFKlkQbG3HBWJgwMpoXshaiJjnLRnF40W4ULmt0zu68uhjpQ8DPbXYK1q1AK2ay+tWQRYDY4i19B3Qs659as+P2nm9fW6xW+WOdHvcIW+LO8QtcXPeDre2W+vyPCyUsr9XNTeLVu1Dq+aqrS9elHisyNbpnXFfb9RXA/EqiHD54jFcOR7mGP05I8Cs2pvltVVFelVBq7w/iDhX9zw8eJjFwOXYIlhY7ds0vO19L7YNI/gWDO9bBdnfvrm0c1z2WH1/gUZe+5tVS6+8tUoqdBOM4K0aQq19ibpNMEqrYvzpbkyNK4beVgmLrQ6OsEI41OpgK9svPKwMbqLE1gu0q7Ze0SqFP5PlvVUeu4f8Yv9NohE1BcxTvGrehtFiuDxFa9Ntgmy9QFS0Kr8irRppQXAo3rcK5hB9MbB2/L8TUA0B8xav0QUeslpAuIDN0aq8amxVpE55WsTLhWuGh8PfIoKLHjAv8aoxXNIw8eIsFoC+aFUexc5mDSjSjotcamnVJrstat/OzoBVkcgB8xCvGsM1FM/hAuAbrdpcba2KMmR5G4zYbeEX/89UKmrAPMSLcHX9uj7D5eU2uFyHBawWuVWl0aquX9fngqCXYY1tgsNgwGpAxICVjtfoZ7OChgs+rbv1Iu3UccttxBStVV4WBGnVuLwMRrs4i+UT/680JlLAvMRrNAHDVftZLABl0Kp+amrVEGjVcNhtMR8DVsOIVzeEazmv2y9yKL1NEECsVpVUS6ta2yqI/iLsthj+cd0LHNm6rBuPPZ71az5w7lTWr9eK2XB95Gm/f2hPXXtOjz9d5q4zu+E6d/aa4V9s+5J0dpj/NC9tX9Hhs77XVS5uS0fO+vk6aBut8iNCq3aHrJKtGqVTEq1y1Krzx7d0zVPrD6YXTh7S0Scub3YQa7h84pgOPVn+esYhFBuwhrAqgkRttd2AEa/5RotXsHBd2k46fNayfk2gVl0HNpq1GK1ajCFr2dekVRhHVQPWKl2iRtAmIsSLIWt93sPlaWVwU5usDF44eURHn7iY+Yjy2zq9oyuPPlb6MKqzrFm0aoJWzceui3HU1Crk1dSA1cW8oLUcMs/xYsjCGDbdegEMgeFrL++tkio/mxWoVbWexSq5TTDKYuCYYvzXUNiikLUUMa/xIlzr4yxW3Wre247lWl4o9HydVsmzWaNtGRxAK61CXRiwNtBixDwPWlWHq9EhC0Ae+3tVe6skn70qtSg4Sqsa3iqYA4NaXRiwMmslYoTrKlYHh+PlLFbUOzQBi7S0QOi1V1W2qtGtgl6GI7a0++Hzb1UVufHY43t+1Mbj80lKPItklOePBHpGFs8byYeHOGIsLfTKE1rVz1DPx/Kk9IOHSzy78fKJMtfSD40Ba2S1BszboFXi4cSEaxilgwO0ilYNj1b1k7tVORcDW28Vi4F7MWAVVlvAPMZrTIRr9uv5Oou1afzOH9/s22WJlcF1pB0/26jgR22Lg7SKVl39er5atalNW1WzrdM7473WaK+ElYjXMAhXfK2vDALe1NIrWgVv6F033hcDGbCcqmW10Eu8CFc3rAwC6ItW5VNdqxo7i1XDcBRlt4V3DFhBRA+Yl3CNGS/ClVeOcNUQvz5qvXgYftGqzdGqbry2CpAYsEKKGq8WVwijhiunms5ildrbzsXDiCjqsEWrBtBQqzwsBtKq8hiwAiNe6yNcy3ldGYx+BoqtF2gVrVoPrVrOa6vQT427LRiwKhE1XiURruVyhsvTWazoQxoQWcSFQVrVDk9nsRAbA1ZlosWr9Aoh4RqPpyFrE9wCF8iDVnVXTaucLwZ6UnKbILstNsffFCoWadgqHa6x4kW4NudhfzuAfGhVN7RqsVp3XCAuBqxGRIhXKyuEhCu+EmexuHgYLaBVq9GqGFgMbBsDVmMirBQSLqwSOTpjb71Y5+Jh7w9wRP28d0qiVd6wGOgDi4ETh0u98DVbl3TzsYezfb0Pn7s+29dqxY3HHtcD506VPoy5dsP1kafH/4veqWvP6fGnh7+jzbHt8zp39pr8X3j7knQ273/al7av6PDZutZjLm5LR86WPgp4R6vK2h2yaNVBY7VqMAO0ypvSnTl/fEvXPLXe4Hnh5CEdfeJy5iNqRzV/spcFkKAtFiFehGsNjsN1aTvp8FkrfRhAEV2GNZp1EK2ab4xWDbYYKGVvVc7FQC+tKj2kYT0+/waWGcPXap7jRbh88HYWq3R0NlkZBJahWYvRqoNo1V60KqbLJ47p0JPj3S1zaE0MWMssClmrEfMaL8K1Bs5iAdWZ16wWe+W5VdL4WwZp1TBqaFXNi4Fp57jssbLPqVvE559oB1qPmMd4Ea41ON5+kcOmK4OsLKIGLffKY6ukMouCtOoqb61Ce/jT18PNxx7e86MFHu/kVOLOTWM+5LEVLd+lad07CXJ3JnTVWq9o1QStyi9Hqza9822pO+fywOH1MWBtoKWAEa7hwzXY7XAzP2/E48OHSynxPCxgHa30ytugRat6oFXVYDGQASur2gPmMVxjxytsuDLLFa6WVwaBklrolRe0qocGHkBcAouB4+MdH1Ct8fI4aI0pZLiIFoAlauwVrQrYqsxYDEQpDFgjIV7Dqi1cg3C6/aLla7HGcvlE4Ge6YXS19YpWBcOCICrAgFVAjfHyoMRe96FEWBn0pOTK4LpbL7h4GBHU1CtalV+EVrEYuBlatR4GrMJqiZeXFcIxwxVy+4XTs1gA/KuhV7Qqvwit8oJtgu1gwHKklniVRrhWcBiuVlcG1+H17kxbp3dKHwJGFL1XtCov763iLNb4N7pYp1U1bWdnwHIqcrw8rBDWFC50w8ogUAatWh+tAurEgBVA1GHLQ7jGiteQ4WJl0D9ugQvQqnXVMmS10qpNsZjXhk5/KzCzO8zsg2Z2r5m9esnn/U0zS2Z2e75DxKxo8SodLmm8eIULF4rzfPFw2jle+hB6o1V+RGxVyV4xZMVRejGQAS2GlQOWmR2S9DpJXyzpNkkvNbPb5nzecUl/X9J7ch8kDooUr9LhkuoYsrJjZRAVoVU+RWqVVHZRsJYhK7sKW1UKd70dT5d3+vmS7k0p3ZdSuiDprZLunPN53yfpByV9LOPxYYVI8WplyBqK9+0XOZReGURotMoxWtVN9E5JbZzF2hSLgfn03W0x1g2ZugxYN0i6f+bnD0x/7ePM7DMk3ZRS+qVlX8jMXmFm95jZPWcfvdD7YLFYlHh5OJs1tFa3X9SwMhglel7vJFgYrQqAVq1Ww46L7K1ydhaLxUCssvGV2Wa2JemHJb1q1eemlO5KKd2eUrp9+/TRTV8ac0SKVwk1hCu7ys5ilRqSvN/ooqbb366DVvlCq5ajVVgmymJgy7r8jeBBSTfN/PzG6a/tOi7peZLeYWYfkvQ5ks5w8XBZEeJFuNbDWay6sLc9G1oVEK1ajFbtw2JgOC3vtugyYL1X0q1m9mwzOyrpJZLO7H4wpfRESun6lNItKaVbJL1b0otSSvcMcsToxXu8ag/XUGoPV0ktRK9StCowWjVf9FZ51fJioPfdFrU4vOoTUkqXzOyVkt4u6ZCkN6aU3m9mr5F0T0rpzPKvMN81dlHPOfqQJOm+C89c50ugh91wffjc9YWP5KAbjz2uB86dGv11n3XtU/rI08PeivrUtef0+NNtb81ax6XtpMNnbe3ff3FbOnI24wHBvTFatQot25z3VkkavVeRW3Vs+7zOnb0m3xfcviSdXflXVwzgwslDOvrE5cFf5/KJYzr0ZPytq53+lKaU7pZ0975f++4Fn/v5fQ+ia7yWIWzdeI0X4erPa7gubV/R4bOskGF8Q7dqlT4to1nLeW2VVGZRkFbNcNSqkouBLCT6Vs0ywKKwEbH5bj72MOGaGiNcaMP541u65ql2t56gu2XDGN26ilZdFXnIAlpTzYC1CIPXYl5XCEttGRwSK4PrfA1WBpe5cPKIjj5xsfRhYADzutVys2jVVVEXBLO3KhN2XGAozf6pes7Rh+b+aJHHC4vHvqB4jAuJw9ypCYA79MrnjTBK3Pxi6F6FuHW7oxszRbybIDe6GB7v8D6tBoxwxR6yssoUrpbv0tQXt2rHOlrulSel7jA4pCFaxWJgftz11i8GrBVai5e3QavGIWsItYZr05VBoCUtDVy0ilaxGDi+dRYDW30WFgNWD63Fy4vawtXSWazSNlndY2UQ0bXQK1o1nBCtcoLFQOzHgLWB2uPlaYWQcK3m8SxWiyuD7G2HRzX3ilYNx32rKjmLxYJeffibQCY1n90iXCiJlcF8Lp/g9suod9iiVWjRusMZi4HD4t0dSG0BazVcQ2JlEN6lnXi3g0Y/NbbKQ69o1XIeW7UpFgMxiwFrBLXEq8VwNb/9ogJRtl5wJ0GUVtOwRavQF4uByIkBa0S1xItwNaTxlUHvw1mrd2fC8GhVHrW0yv1ioJNWbcJ7b8ZUw3Z2BqxCosfLw9kswhVHayuD7G1HLaIvDNKqfFpo1aZa2CbocbeFx+3s/C2gsMjhksqvEBKu+VgZBJBb5F611KpIvG1pj7oYyI0u/OGddSLyKiHhwhg2WRlk6wWQD61az1itirQYmBWLgW61uJ2dAcuhiPEiXJtzfRYrg6grgwDmo1X91dCq3Ly1alMsBkJiwHItWrwIV8VYGRyFx73twCq0qp/orar9LBaLgciBASuASPFqJVxDqf0sVkSsKALdRGtVyV7Rqr1oFWrDgBVIlHi1EC7OYq1n05XBSFsvuHgYrYrSKqn8ouDQmmxVo1gM9IW/AQREuHxocvsFAHQUZdAq1Sq2Cu6V7SwWi4GDYzv7agxYQRGu5aJvv8iJrReYVcMDHBELrVqMVqE0dlsMg3c1OMK1WOStgm7PYnGzCwBriLAoSKv6o1XAfAxYFSBci0UOV06ezmKx9QJol/de1dyqCGpqFdrGgFURz9GS6r8mKzdWBv3wPJi1+ABHxOe5V7W2qrmzWIWVWAz03KrWMGBVhtXBgziLNeFpZbAV7G0HFvPcK1pVjqebXSCf1hYDqX+lvEZLqjdcQ6h1ZZCtF4txdya0xmuvaFV3tAq5Rb8h0+FSL3zULumWI498/OcfunhdqUOp1m607rvwzMJHctDNxx7Wh89dP+pr3njscT1w7tRgX/9Z1z6ljzx9fLCvn8Ox7fM6d/aazb/Q9iXpbLFvH7q0nXT4rBV7fbRjf6tWoWXr8dorWlVGtlahk/PHt3TNUwyTOZX7G9I+fQI2i5itRrjGM0S4Tl17To8/HXslx5OL29KRs6WPIo60c1z2WLdtRVund3Tl0ccGPiLfVrWMZi33nKMP0aqgXLaKxUAUEn6L4C1HHpn7Awd53IYx9haMqNsvahRt6wUXDyMHmrUarRq+VRGuxfKiVKu83+iC7ezLhR+wFiFg83m8qJhwrZZzfzsXEAP+0Ky9aBWt4sZMiKzaAWseAnYV4eJMVnSb3AIXiIBe0SpkwGIgCmhqwJqn5aGLcA2HlcFuom29WAe3akcuLXZK8teqMXEWy0erNsFiYJso/xwtBczbNowxhyzOYmXAyqALrT1fBO2d3Wq1UxKt8iDaNcMojwFrCeJVvyHDxcpgfbh4GF610CtPC4LsuFjO3XOxgi4GstsiLt7NjmoPl+RnyKopXBgWWy+Ag2rvVYutavksFouB/XHX2/IYsHqqfZWQcOXFyuBqbL0AhlFzr2hVXty2fVgsBg4v7fh6eDYD1gYI17A4kzW8llcGWeFDS2rsFa3yjcVA7NfS9cIMWBnUuEroJVxjYWUQQAtoVWyttqrlxUDPol4vvHV6Z/jXGPwVGlNTvDyEq5btF7m5WxlEKJdPHCt9CCispoVBWuUXrdocuy1iYsAaCOHKp4btF9WvDBa8QxN724HN1NArWpWH51ZlEfRugoiHAWtghCsWVgbLiLS3ndVE1Cp6rzy0aqwhK1KrcvGwTXCTVrEY2BYGrJEQrs2wMojIeL4IIoncq9KtqkHuVtW0GBjJOouB3loVeTu7r3eyAVGjJZUPFyuDw4m+TXBdnI0CFovaK1oFoDQGrAJYHVxf9DNZrAyiq6h3Z0JdovaKVm3Ga6taXQxEPAxYBREuv1gZHB972wG/IvaKVmEIJVrFbot4GLAciBYtqWy4om+/qHplEMW09ABHlBOtVy20aihcNwysjwHLCVYH+4kermqx9QKoXrRe1d6qKGexWAxESxiwnCFcvkQ5i4VxsV0DiHU2i1ZhDxYDs+F64fkYsJwiXKtxFuuqWlYGSzwPa6xhydvtb4EcIi0K0qr+WAycL8qzG1kMLIfiO0a4Vou8/aLacBVaGeRGF0A5UVpVSuRW5VTLYiDW18r1wgxYAUQZtGrfggEAWIxO1YnFwLy830mQ3RZ58C4GEiFeJbAyOMHKIIDSIiwIsuMCwNAYsIIhXHWpdmVwA1H2tgNYjFZhERYD0QIGrIAI10GsDDrTwB2auHgYWI5WHRT1hhcsBh7EYiCWYcAKynu4SogarpxyrQxGxY0uAF9o1UFDt4rFQKA8BqzAPIer1u0XQ4TL48pga1svOBu1XNo5XvoQEJjn67JqbdUQcrbKzWJgsBtdIA4GrOC8Rkti+wWwKR7giJp47VWNrWrlLFZri4HrYAGxDAasCniNFrqrcmVwA7Xubef2t2id115xJgsRMCzFQe0rQbSuYmXQiQZudBHZ5RPHSh8CGuW1V2OL2CoWA/eqdTEQm2PAqojXaLEyGBNbLwAMxWOvaFWjWAzcGNvZD2LAqozHaJXQ+spgy7h4GIjBY6/GHrIitsobFgPzYzv75oq9g9dYqVeuH9FCDVsvAA+uMek5h9PKH1iPx15hObYJ5sFiYN0Ol3zxrlG67xLTWF+3HHlEH7p4XenD2OM5Rx/SfReeOdrr3XzsYX343PWjvR7m2L4knV3v28yl7Ss6fNb3KtrFbenI2dJHcdWFk0d09ImLpQ+jSct6RsOW89YrWgX4cfnEMR16Mt4g7vtvL1OsHq6HlcFht1+0sE2wpa0X3J0JQ6FXq3nrVU27LlrYJlgSN7rAPCEGrC6I13xEq10tb70APGPgms9br8YU7RmO1S0GBrrRBYuBMVQzYM0iXHt5i1ZNQxYrgwA2xcB1lade0apxsBg4PIay8VU5YM0iWhOeojU2VgYBRNJ6s1rtVbRWAeu6cPJI6UMYXKcBy8zuMLMPmtm9ZvbqOR//VjP7gJm9z8x+zcxuzn+om2t9ldBTtGpaGaxdxK0X3J2pPbV0albLvfKCVi3GYiCw2MoBy8wOSXqdpC+WdJukl5rZbfs+7bcl3Z5S+ouSfk7SP859oENoMV6tDlnRbnaRS/StFzVePMzzRfKruVO7aFU5tGp4OVoV8aZMLAbmlXaOlz6Ej+tS+udLujeldF9K6YKkt0q6c/YTUkq/kVL66PSn75Z0Y97DHF5Lw5ancGE+Vgax68LJQ6UPIYImOiXRKmAwG+y2qHExEJvpMmDdIOn+mZ8/MP21Rb5W0i/P+4CZvcLM7jGzex5+xO8fxhbi5SVcrAzCEy4EDitbpyRa5Qmt8o3FQGC+rHtVzOyrJN0u6YfmfTyldFdK6faU0u3XX+d/mwzhArAKQ1ksqzolxWxVzb2iVfmwGAiMo0s5HpR008zPb5z+2h5m9gWSvlPSi1JK8TbCLlB7uDxgZXBYre5tR1Oa7tQuWjUsbngxrOjXDNeG64U30+Xde6+kW83s2WZ2VNJLJJ2Z/QQz+3RJP6FJtKr8DlRruFgZ9KuqrReBHuKIkOjUVK2Lgq21KtJioKdWtbIYyM4J/1YOWCmlS5JeKentkn5P0s+klN5vZq8xsxdNP+2HJG1L+lkz+x0zO7Pgy4VWY7QkH+GqYWWQrRfDWPfiYe7O1A46dVCNgxatyoNWLcFi4Nq4IdNeh7t8Ukrpbkl37/u175755y/IfFxu7QbrvktW+EiwrpuPPawPn7u+9GEAyIhOzfecw6mqXt1y5BF96OJ1pQ9jFLQKOV3clo6cLX0U7WCD5ZpYGcyvhpVBr9jbDrSrxrNZJdEqzMNuC8xiwNoA0cJ+ubdeeNrbDiC2WnrlYUEwOq+tYjEQtWDAyoBo5TPWymCkC4i9aOXiYaBm9CoPWuUXrYIHDFiZ1HI2q3S0UDEuHgZcqKVXAOAVA1ZmRGtz0fe3c4em2Lj9LVoRvVelFwRpVaVYDEQGDFgDiL46WDpaY4my9YK97etfPIzlLp84VvoQUFjkVrWCVgHxMGANKHK4Sg9Z0VcGgVIunDxS+hAQDK1aH60CMA8D1sAihwvrq3XrRbSLh8e6/W3fbYXnj/OtF75EblXpIWsMUc5ieVGqVWPutmA7u29UfgRRw1U6WmOsDBItAJiIvr0d/dW6GBgNz8KqDwPWSIgWNlHN3nYuHgbci9irkguCbBO8imc3AhMMWCMiWj5xFgsA9orYK9Sj+GIgJLGdfRO8cyOLGC1WBtfD1gvkcuHkodKHALhXe6uGWgykVXNUutuC67bGw4BVQMQhC9jVwsXDAGgVAKyLAauQaOGqfWUwAva2AxgbrQIwhNofKcKAVVC0cNWM67AAYD5a1U3kbYI5eVkMjPZYEdSFAQudsTLYn8e97Vw8BOkIYAAACiJJREFUDKCvSEMWrerPY6uiYjs7JAas4iJFC4honeeLcCHw+rZO75Q+BAyEXq3GlvZ8WAxEZAxYDkSKVqmVQaJVkUrvzgS0IEqvaj6LFWGbINA6BiwnokSrZhGi5WVvOwAAKMvbbouhHily+cSxQb7ukBiwAPTGxcNAu6IsCLLjopwqFgNH3m2xznZ2+MWA5QjRWi5qtLh4GEBtovQK3dEqIB8GLGeIFgAgggi9qvVarAhb2j1gtwVKYcDCWohWbFHvzuT59rfe9sIDQHRRWwUwYDkUYVWwlKjbBAGgRvRqPloFtI0By6kI0ar1LNYQcu5tr+LiYQDV8N4rWgVgbAxYjnmPFgAAKGOILe0sBubheTs7xsGAhXCG3nrRynVYRfGwYZ0/Huvbb9o5XvoQ4BgLggexTRBoV6zCN8h7tNh60S7uzgQgCloFYEwMWAF4H7IQU0t3Z+IBjsDwaBWqU+FuC+54Ow4GLGyMlcFueIgjAJRDq9rEbguUwIAVBCuDe7V+HVbLFw+PhVU+oD9atVfrrcqhpd0WHkW7XtgL3rVAPIeLlUEAAOJjtwWwOQYsAACQjefFQIyL3Rb9cL1wPRiwAABAE9htAWAMDFjBsDKI1vEAR8A/WnUVz8MC2sOAhbC4eBgAALSMGzL5xICFbNh6sRoXD2MMF04eKX0IgNuzWLW1isVAwB8GrIC8RgvxcPtbAMB+LAYCm2HAAoLycHcmHuAIAACwFwMWgDK2L5U+AgADY8fFBDe62EzE3RbckKltDFhBEa0JooVWXDh5qPQhANWo7Tosz9htgRYxYCErogUAmMWCIEJjtwXWwIAFLMHdmQAA3tEqwBcGrMBYFQQAAAB8YcACRsbtb+vFAx+B+VgQBNASBixkx3VYAADvxm4VN2UC2sGAFRyrgkQL3Vza5r8VAOiK3RbA+hiwAAAAACATBiwgMA/PFwGALthxAaAVDFgVIFrDqv32t6euPVf6EAAAAKrBgIVBcKMLAIB3tArAEBiwAAAAgqt9twUQCQMWAAAYBVva28T1wmgNAxYAAABc43phRMKAhSrwLKxyjm2fL30IABACrQLawIBVCbZdAAAAxHZpm7/P1YABCyjgxmOPlz4EAACArC6cPFT6EFxgwAIQzqXtK6UPAQAAYC4GLABY4OJ26SMAAADRMGABAIDRcM0wSuCGTBgTAxYAAGjWLUceKX0IACrDgAUAAIADuCFTndj+PjwGLADlbF8qfQRFnT/Ot2AAAGpD3TEYtl0AAACgNQxYFeHCYQAAAKAsBiwAAAAAyIQBC+jg5mMPlz4EAACAalw4eaT0IQym04BlZneY2QfN7F4ze/Wcj19jZm+bfvw9ZnZL7gMFAGAZWgUA8GDlgGVmhyS9TtIXS7pN0kvN7LZ9n/a1kh5LKf05ST8i6QdzHygAAIvQKgCAF13OYD1f0r0ppftSShckvVXSnfs+505JPzn955+T9AIzs3yHCQDAUrQKAODC4Q6fc4Ok+2d+/oCkz170OSmlS2b2hKTrJO25cMXMXiHpFdOfnt+54YHfXeegG3a99r2n/t2/+lOGFfA9K473rD/es/7+u8xfj1b5EfC/B1oVEO9Zf7xn/a3Vqi4DVjYppbsk3SVJZnZPSun2MV8/Ot6z/njP+uM964/3rD8zu6f0MSxCqzbDe9Yf71l/vGf98Z71t26rumwRfFDSTTM/v3H6a3M/x8wOSzopiafMAgDGQqsAAC50GbDeK+lWM3u2mR2V9BJJZ/Z9zhlJf2f6z18h6ddTSjz1FgAwFloFAHBh5RbB6T71V0p6u6RDkt6YUnq/mb1G0j0ppTOS3iDpp8zsXkmPahK2Ve7a4LhbxXvWH+9Zf7xn/fGe9Zf1PaNVrvCe9cd71h/vWX+8Z/2t9Z4Zi3cAAAAAkEenBw0DAAAAAFZjwAIAAACATAYfsMzsDjP7oJnda2avnvPxa8zsbdOPv8fMbhn6mLzr8J59q5l9wMzeZ2a/ZmY3lzhOT1a9ZzOf9zfNLJlZ87cp7fKemdmLp3/W3m9mPz32MXrT4b/NTzGz3zCz357+9/nCEsfphZm90cweMrO5z5GyiddO38/3mdlnjH2M0+OgUz3Rqf7oVH90qj861d8grUopDfZDkwuN/1DScyQdlfSfJN2273O+UdKPT//5JZLeNuQxef/R8T37nyR9wvSfv4H3bPV7Nv2845J+U9K7Jd1e+ri9v2eSbpX025J2pj9/ZunjDvCe3SXpG6b/fJukD5U+7sLv2V+R9BmSfnfBx18o6ZclmaTPkfQep/+/0qn+7xmd6vmeTT+PTvV4z+jUWu8ZnTr4vmVv1dBnsJ4v6d6U0n0ppQuS3irpzn2fc6ekn5z+889JeoGZ2cDH5dnK9yyl9BsppY9Of/puTZ730rIuf84k6fsk/aCkj415cE51ec++TtLrUkqPSVJK6aGRj9GbLu9ZknRi+s8nJf3JiMfnTkrpNzW5W98id0p6c5p4t6RTZvZJ4xzdx9Gp/uhUf3SqPzrVH51awxCtGnrAukHS/TM/f2D6a3M/J6V0SdITkq4b+Lg86/KezfpaTabqlq18z6anc29KKf3SmAfmWJc/Z8+V9Fwze6eZvdvM7hjt6Hzq8p59r6SvMrMHJN0t6ZvGObSw+n6/K3UMdGovOtUfneqPTvVHp4bRu1Urn4MFv8zsqyTdLumvlj4Wz8xsS9IPS3p54UOJ5rAm2y8+X5PV5980s7+QUnq86FH59lJJb0op/RMz+1xNnrn0vJTSldIHBpRAp7qhU2ujU/3RqREMfQbrQUk3zfz8xumvzf0cMzusyenKRwY+Ls+6vGcysy+Q9J2SXpRSOj/SsXm16j07Lul5kt5hZh/SZP/smcYvIO7y5+wBSWdSShdTSn8k6fc1CVmrurxnXyvpZyQppfQuSc+QdP0oRxdTp+93Do6BTu1Fp/qjU/3Rqf7o1DB6t2roAeu9km41s2eb2VFNLg4+s+9zzkj6O9N//gpJv56mV5Q1auV7ZmafLuknNIlW6/uNpRXvWUrpiZTS9SmlW1JKt2hyPcCLUkr3lDlcF7r8t/nzmqwKysyu12Qrxn1jHqQzXd6zP5b0Akkys0/VJFx/NupRxnJG0sumd2j6HElPpJT+dORjoFP90an+6FR/dKo/OjWM3q0adItgSumSmb1S0ts1ubPJG1NK7zez10i6J6V0RtIbNDk9ea8mF5i9ZMhj8q7je/ZDkrYl/ez0Ous/Tim9qNhBF9bxPcOMju/Z2yV9kZl9QNJlSd+WUmp21b7je/YqSa83s2/R5ELil7f8F3Eze4smf/m5frrf/3skHZGklNKPa7L//4WS7pX0UUlfPfYx0qn+6FR/dKo/OtUfnVrPEK2yxt9TAAAAAMhm8AcNAwAAAEArGLAAAAAAIBMGLAAAAADIhAELAAAAADJhwAIAAACATBiwAAAAACATBiwAAAAAyOT/B5S5QkkrcHRyAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pr1Zyw-2Aaw9"
      },
      "source": [
        "def neumann_rhs(a, b, Point, grad_u, tol):\n",
        "  '''\n",
        "  returns the average flux in the two Neumann edges sharing 'Point' \n",
        "  '''\n",
        "  # on the domain boundary we define the outward-pointing normal vectors\n",
        "  # corresponding to the two edges sharing 'Point' \n",
        "  if (Point[0] < tol) and (Point[1] < b-tol):\n",
        "     normal_l = np.array([-1, 0])\n",
        "     if Point[1] < tol:\n",
        "       normal_r = np.array([0, -1])\n",
        "     else:\n",
        "       normal_r = normal_l\n",
        "  elif Point[1] < tol:  \n",
        "     normal_l = np.array([0, -1])\n",
        "     if Point[0] > a-tol:\n",
        "       normal_r = np.array([1, 0])\n",
        "     else:  \n",
        "       normal_r = normal_l\n",
        "  elif Point[0] > a-tol: \n",
        "     normal_l = np.array([1, 0])\n",
        "     if Point[1] > b-tol:\n",
        "       normal_r = np.array([0, 1])\n",
        "     else:\n",
        "       normal_r = normal_l\n",
        "  else:\n",
        "     normal_l = np.array([0, 1])\n",
        "     if Point[0] < tol:\n",
        "       normal_r = np.array([-1, 0])\n",
        "     else:\n",
        "       normal_r = normal_l\n",
        "  \n",
        "  tangent_l = np.array([normal_l[1], -normal_l[0]]) * tol \n",
        "  tangent_r = np.array([-normal_r[1], normal_r[0]]) * tol \n",
        "  flux_l = np.dot(grad_u((Point+tangent_l)[0], (Point+tangent_l)[1]), normal_l)\n",
        "  flux_r = np.dot(grad_u((Point+tangent_r)[0], (Point-tangent_r)[1]), normal_r)\n",
        "  flux = tol*(flux_l + flux_r)\n",
        "\n",
        "  return flux"
      ],
      "execution_count": 9,
      "outputs": []
    }
  ]
}