3D7 - Heat Equation - very large mesh.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3",
      "language": "python"
    },
    "language_info": {
      "mimetype": "text/x-python",
      "nbconvert_exporter": "python",
      "name": "python",
      "file_extension": ".py",
      "version": "3.5.4",
      "pygments_lexer": "ipython3",
      "codemirror_mode": {
        "version": 3,
        "name": "ipython"
      }
    },
    "colab": {
      "name": "3D7 - Heat Equation - very large mesh.ipynb",
      "provenance": [],
      "include_colab_link": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/chrisrichardson/FEMexamples/blob/main/3D7%20-%20Heat%20Equation%20-%20very%20large%20mesh.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HfbyMIaJSuFx"
      },
      "source": [
        "# Heat Equation - with a larger mesh\n",
        "### Sparse linear algebra\n",
        "\n",
        "First import the required libraries for basic algebra and plotting. For larger problems, we use the scipy\n",
        "sparse package, only storing non-zero entries in the LHS matrix."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "Rm4W6yXFSuFz"
      },
      "source": [
        "from numpy import *\n",
        "from scipy.sparse import csr_matrix, lil_matrix\n",
        "from scipy.sparse.linalg import spsolve\n",
        "set_printoptions(suppress=True)\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "imL6xqS2SuF8"
      },
      "source": [
        "### Mesh\n",
        "\n",
        "The mesh is made the same way as before."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "rAzLybLHSuF9"
      },
      "source": [
        "\n",
        "# Make a rectangular mesh of triangles, nx by ny \n",
        "nx = 250\n",
        "ny = 250\n",
        "\n",
        "c = 0\n",
        "geometry = zeros((nx*ny, 2), dtype='float')\n",
        "for i in range(nx):\n",
        "    for j in range(ny):\n",
        "        geometry[c] = [float(i)/(nx-1), float(j)/(ny-1)]\n",
        "        c += 1\n",
        "\n",
        "ntri = (nx - 1)*(ny - 1)*2\n",
        "topology = zeros((ntri, 3), dtype='int')\n",
        "\n",
        "c = 0\n",
        "for i in range(nx - 1):\n",
        "    for j in range(ny - 1):\n",
        "        ij = j + i*ny\n",
        "        topology[c] = [ij, ij+ny, ij+ny+1]\n",
        "        topology[c + 1] = [ij+1, ij, ij+ny+1]\n",
        "        c += 2\n",
        "        \n",
        "mesh = (geometry, topology)"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZqZhkoQRSuGC"
      },
      "source": [
        "Here is a simple plotting routine to display the data values"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "CsJZ5TAZSuGD"
      },
      "source": [
        "def plot(mesh, data):\n",
        "    geom, topo = mesh\n",
        "    x = geom[:,0]\n",
        "    y = geom[:,1]\n",
        "\n",
        "    plt.gca(aspect='equal')\n",
        "    plt.tricontourf(x, y, topo, data, 40)\n",
        "        \n",
        "    xmax = x.max()\n",
        "    xmin = x.min()\n",
        "    ymax = y.max()\n",
        "    ymin = y.min()\n",
        "    dx = 0.1*(xmax - xmin)\n",
        "    dy = 0.1*(ymax - ymin)\n",
        "    plt.xlim(xmin-dx, xmax+dx)\n",
        "    plt.ylim(ymin-dy, ymax+dy)\n",
        "    return\n"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WiEMgjZ5SuGH"
      },
      "source": [
        "### Element matrix\n",
        "The \"element stiffness matrix\", $K$ depends on the conductivity tensor $D$ and the triangle geometry.\n",
        "\n",
        "$$ K_e = \\int_{\\Omega_e} B_e^T D B_e d\\Omega $$\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "Dn_luRBKSuGI"
      },
      "source": [
        "def Kmat(D, p, geometry):\n",
        "    ''' Calculate K from the D matrix and the points '''\n",
        "    # D contains the conductivity tensor\n",
        "    # p contains the indices of the three points\n",
        "    # making up the triangle\n",
        "    assert(len(p) == 3)\n",
        "    x0, y0 = geometry[p[0]]\n",
        "    x1, y1 = geometry[p[1]]\n",
        "    x2, y2 = geometry[p[2]]\n",
        "\n",
        "    # Element area Ae\n",
        "    Ae = 0.5*abs((x0 - x1)*(y2 - y1) - (y0 - y1)*(x2 - x1))\n",
        "\n",
        "    # 'B' Matrix - representing the 'gradient' operator\n",
        "    B = array([[y1 - y2, y2 - y0, y0 - y1],\n",
        "               [x2 - x1, x0 - x2, x1 - x0]])/(2*Ae)\n",
        "\n",
        "    K = Ae*matmul(B.transpose(), matmul(D, B))\n",
        "    return K"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-8TMJiHQSuGN"
      },
      "source": [
        "### Matrix assembly (LHS)\n",
        "The process of combining the local element matrices into a larger global matrix is called \"assembly\". It is really just a case of relabelling the indices with their global values, and adding the local matrices together. Now, we use the `scipy.sparse.lil_matrix` which only stores the non-zero entries."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "UutfeMhPSuGN"
      },
      "source": [
        "def assemble_matrix(mesh, Dmat):\n",
        "    geom, topo = mesh\n",
        "    np = len(geom)\n",
        "    Kglobal = lil_matrix((np, np))\n",
        "    \n",
        "    for tri in topo:\n",
        "        K = Kmat(Dmat, tri, geom)\n",
        "        for i, idx in enumerate(tri):\n",
        "            for j, jdx in enumerate(tri):\n",
        "                Kglobal[idx, jdx] += K[i, j]\n",
        "    return Kglobal"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "i064Egt_SuGS"
      },
      "source": [
        "Dmat = array([[15.0, 0.0],\n",
        "              [0.0, 5.0]])"
      ],
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "ASagR8n1SuGY"
      },
      "source": [
        "Kglobal = assemble_matrix(mesh, Dmat)"
      ],
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1SrxV4IKSuGe"
      },
      "source": [
        "### RHS assembly\n",
        "\n",
        "As before, we perform RHS assembly on the vector $f$, combining the local vector entries from each element. Each element has a constant source term $s$.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "fmUyrUI3SuGe"
      },
      "source": [
        "def fvec(s, p, geom):\n",
        "    assert(len(p) == 3)\n",
        "    x0, y0 = geometry[p[0]]\n",
        "    x1, y1 = geometry[p[1]]\n",
        "    x2, y2 = geometry[p[2]]\n",
        "\n",
        "    # Element area Ae\n",
        "    Ae = 0.5*abs((x0 - x1)*(y2 - y1) - (y0 - y1)*(x2 - x1))\n",
        "    return s*Ae/3.0*ones(3)\n",
        "\n",
        "def assemble_vector(mesh, s):\n",
        "    geom, topo = mesh\n",
        "    fg = zeros(len(geom))\n",
        "    for tri in topo:\n",
        "        f = fvec(s, tri, geom)\n",
        "        for i, idx in enumerate(tri):\n",
        "            fg[idx] += f[i]\n",
        "    return fg\n",
        "\n",
        "fglobal = assemble_vector(mesh, 3.0)"
      ],
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "ABII0pL7SuGj"
      },
      "source": [
        "def set_bc(K, f, row, val):\n",
        "    # Set BC\n",
        "    K[row,:] = 0.0\n",
        "    K[row, row] = 1.0\n",
        "    f[row] = val\n",
        "\n",
        "# Set zero BC on bottom edge\n",
        "for i in range(nx):\n",
        "    set_bc(Kglobal, fglobal, i*ny, 0.0)\n",
        "\n",
        "# Set u=10 BC on left edge\n",
        "for j in range(ny):\n",
        "    y = (j/(ny+1))\n",
        "    set_bc(Kglobal, fglobal, j, 0.04*sin(2*pi*y))\n"
      ],
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i9WkRimpSuGo"
      },
      "source": [
        "Finally, we can refer the $K.u = f$ problem to a sparse linear algebra solver which will use an LU method to get the answer. The solver requires the matrix in CSR format. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "T2ooql38SuGo",
        "outputId": "60151362-b8e1-4913-ca66-8485e6b2d75f",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 51
        }
      },
      "source": [
        "u = spsolve(Kglobal.tocsr(), fglobal)\n",
        "print(u.max(), u.min())\n",
        "print('Number of DOFS=', len(u))"
      ],
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "0.07354042572686259 -0.039999216713052996\n",
            "Number of DOFS= 62500\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "lSjzfqo8SuGt",
        "outputId": "9783aae3-c9a6-46bf-e196-2760fe416ae5",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 700
        }
      },
      "source": [
        "plt.figure(figsize=(12,12))\n",
        "plot(mesh, u)"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAArEAAAKrCAYAAADml9XQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dW6xm6X3X+d9Th67uuBO3Dz0C2Q6OhCORYUYT1Mow4oJIgZGTC/uCw9gogjAZLI3GiCEZpCBQQOFmAA1okMzBDFEACULIBWoJI18MiSIhErmjaDLYKKhlhrgNUppud0dtd1dX1V5zset17dq9D+t933V4nmd9PpKlrqq39l6uXbX3t/71X88qwzAEAABacmPtCwAAgH2JWAAAmiNiAQBojogFAKA5IhYAgObcWusdv//97x8+/OEPr/XuAQCo3C//8i//52EYnr3ox1aL2A9/+MN54YUX1nr3AABUrpTyHy77MesEAAA0R8QCANAcEQsAQHNELAAAzRGxAAA0R8QCANAcEQsAQHNELAAAzRGxAAA0R8QCANAcEQsAQHNELAAAzRGxAAA0R8QCANAcEQsAQHNELAAAzRGxAAA0R8QCANAcEQsAQHNELAAAzRGxAAA0R8QCANAcEQsAQHNELAAAzRGxAAA0R8QCANAcEQsAQHNELAAAzRGxAAA0R8QCANAcEQsAQHNELAAAzRGxAAA0R8QCANAcEQsAQHNELAAAzRGxAAA0R8QCANAcEQsAQHNELAAAzRGxAAA0R8QCANAcEQsAQHNELAAAzRGxAAA0R8QCANAcEQsAQHNELAAAzRGxAAA0R8QCANCcayO2lPKTpZTfKKX8m0t+vJRS/kYp5cVSyq+WUn7X9JcJAACPjJnE/lSSj17x49+f5CMP//epJH/r+MsCAIDL3bruBcMw/EIp5cNXvOTjSf7BMAxDkl8spTxTSvmtwzD8p4musSn/zaf/+mPfvvn2fO/rwRP1vL2TA37ugzsHvJ/bB/ycO8N+rz/g/8vJEyd7vX64s9/rk6Tc3v/n3Lpzf++fc/v2g71/zlN37u39c5LkqduH/bynn7h70M9Lkm+7ffjPffr2Wwf/3CR55vabR/38JHn3rePfRpK859bXJ3k7O8/c/Makby9J3nvzjcnf5lnvm/ntn/fMjcN/703pvTf2+5xYo3ffOOALyAY89Vv//dqXsKgpdmI/kOQrZ7790sPve4dSyqdKKS+UUl54+eWXJ3jXdZszYAEAtmzRG7uGYfjsMAzPDcPw3LPPPrvku16cgAUAmM+16wQjfDXJh858+4MPv2+TxCsAwPymmMQ+n+SPPjyl4HcneX2r+7BLE8zzuOHXFQCqd+0ktpTyj5N8b5L3l1JeSvIXktxOkmEY/naSzyX5gSQvJvlGkj8+18XSl5t397+568a9w27uAgD6MuZ0gk9e8+NDkv9lsitiLzffnv6UAgCA2nliVwesFQAAWyNi+SYxvJxy1x89ADiGr6SdWDtA3QwFACxpiiO2OOPm3cefhPLgTlnufduPvdCNu2Xvp3YBAHUTsTNbM2oPIYQBgBaI2IXNHbWtReghx2wBANiJXdnNu8M7wvbot3nkfuqhP99eLACwFBFbialjdu0bveZ0497aV/BON972RwkAluQrb2WmjNljQrbnCAYA2idiKzVVzC4do4esFNy8O/11HMtqBADUTcRO7NZbQ269NeFawMT7snu9byFXleGeP64AsOOr4kx2MTtF0B4bsmL09KxYAKAfInYBU8TssesFh4bsIT9viX+Kr/HmLgBgOSJ2QVPF7ME/t+KJbI17sYz35t3ba18CABsjYldwbMwuHbK1TmNbV+764wcAh/JVdEXHxOyaN3xxsRrPir1/10P5AOhTfV91N2jpkF1qGrv3+9hzpWDfvdh9b+7a6jT53r2ba18CAFxLxE7s5ttDbr69f1weOpU99IavJUJ2qxEIAMxPxM7kmJg96P0tFLIAADUQsTM7JGaPmcrObe5p7NwrBQBAH0TsQpaayu4bslubxrb+0ANP7QKAU74iLujQqeze72fmkN3SbmzL196CN+8ddr7sG2/fmfhKAGiNiF3BEiG7r5omsi0/+GDfY7acFQuwv9dPGv5CwWR8BV3JvlPZfUN27v3Ymqax9mIBYHtE7MRuvnWSm2+djH/9niG7T8zWtlaw19ue8S/Zre/FAgAidjZzhWyy31S2phu97Jeuw1O7AOiRiJ3RPlPZmkJ2H7Xs0s65UiC+2Xnj3pNrXwIAD4nYBewTsnPtye4TsrWEaas3eO17cxfb8dq9p9a+hG/62v13rX0JAEfx1XYhc+7Jjn67M4XsPq+tZapZ017svicULHFW7L17N/f+OW/ePey4LAA4hIhdUA0h2zOnFADAdojYhc25JzvqbTY2ja1lpaCWCTK8fr+elQSANYnYid28+2Dc6yYO2RrWCqAFv3nP074AeiBiZ9BbyI5+mxVMY/dZKahpLxYA2I+IncnNuw9GxeyaITuWaez+nFAAx3v1wdNrXwJQMV9pZ1ZzyM6xVtDiSQU12PeEgn154AEAvRGxC+glZMeaY3Jbw0pBTdG9xDFbAFAzXwkX0kPIzhGnNYUh63jznvNlWccr1hWgaSJ2QWuE7FrWnsbWwF4sAMzHV9mFLR2yLUxj5+DBB8c75KldS3rjbUdlAWyZiJ3YjbdGROqEITvGWiE79nVrrxS0uhcLcIhXTxwvSB9E7AxuvPXg2pidKmTXPHprDa2tFOxj7hMKlvDmXfutwDJeP+n4CwKjtP9Vs2K1hewYPU1jrRQAQL9E7MxqCtk1j92q3VxP75rz5q59j9lyViwAPRGxC2gxZMdYaxq79kqBvVgAWJ+IXcgUIbukXqaxVgoAoE8itiLXhewWprGj3699fgDYNBG7oCmO31oyZNeYxq75T/Vz7cXuo4cTCgBgCb5iTuzG3av//XpMyF5nyjNkr31fI0K2lQcgTGmf2G75yV21P/AAWvLaiQd0wJTa/epasWNDdoobvZZeKxhjrRu8et6L3feEAgDoha+AM1kiZK+z5FrBFqexLar1mK037y37kITfvGciBtA6ETujuUN26QchTKH22K1hL7Z1vT+16417T659CVV67cG3rH0JwMaI2JnVELLXqXEa28JKwVw3obm5CwCu56vlAq4L2bmZxq6v5Zu7AGr1+onzFrfMV9YKmMZebK3jtqwUAED9ROxC1l4rqO20gilj1ykF2/XG227QAtgqETuxG29eXkpLnCF7lSnWCpxU8Lg1H86w45ittrx276m1LwGgC776zeCYkL1KLWsFS6ohErfOAw8AqJGIXcFVIVvDWsF1pprGrrFSMNYce7H73Nw15wkFtZ4VCwD7ELEzuWoae+3PnXmt4Dq1TWOnZC8WAPogYmdU61pBa9NYKwUAwHkidkVzrhUcq8Vp7NQrBWOJbABYnoid2VprBbVMY2s0dqXAebGH6/3Rs0A9PPBgu0TsAmpdK7jOFNNYKwWPm+vJXT0cs/XmPeELS3n1xF/SaV/7X/kacWjIznmTV4/T2LVWCuYw5wkFANA6XyUb18s0Fhjv9fsemAAgYqf25luX/lCv09gpLL1SMPVebGvrDs6KBaB1InYOV4Tsoa4K2SWe5HXl269spYDpeWoXALURsQtb6yavK3/uNdPY1lYKetqLhTl97f671r4EgIOJ2LlUtlbQwjR2K6cUzHVCAQBsia+mc5phreAqtU9ja9PCI2j3OaGgh2O2lvSb9+6sfQkAHMFXvZX0OI299v03OkXt9eauXrzxthiFrfPAg20SsXNbeBp7jGOnsTWtFNiLBYC+idgVzTGNPWalAACgFSJ2CRVNY60UPK6FvdhWvXnXY2QBmI+IXVlt09hWVgpaV8MJBR54AEDL1v9KuhWmsZNaei/WzV2c9ca9J9e+BIDNE7FTe2v/cmptGnusLUxagTa88uDpRd/faydO04CpiNg5XBayFU1jj7HESgGP1HJW7BKPnn3znj1a4DCO2doeEVuJOR5He9U0tvaVgiWntW7ugnq9uvCkdEtePRm3JgW1ErFzWWgae8zDD67SwkqB82IBYLtEbEWumsbO4Zhp7BZXCsbe3AUAzE/EruGAaewWH37Q6g1gY08oqOGYLQBola+ic5r4pIKlzb1SAABwKBG7lomnsYdac6Vgqb3YUW+nnr87AHAgJxRsi4id2wHT2EP0ulIwBTd3AUB/RGyFrBRw3j5nxQLAFvjKuIQJj9s69Aavq8y5UnCsVm/uasX9u7fWvgQAOIiIrdTU09i1VgqW2ItdkmO2AKAOInZtDTyKtvaVgqlu7urFnI+eBYBa+Go3seEbb178AxPe4DXHKQW9G3Nz15QnFDgrFmAdTijYDl9BK3bISsGhpxTUvBcLAHCeiK1BAysFxzj2EbSt7c3yyJt3b699CQB0SsTOoMeVgjn3YkUqwDpePXGzKu0SsZWbeqXgKsesFLAN9+7dXPsSACCJiK3HQisFvT69a8kTChyzBVA3N3dtg4idSY8rBVdxcxcAsCQRyyJquLlrzDFbAEAbRGwDltyLvcoWbu6a8qxYAGA+InZGe68UVLAX6+auepW7/rgCwM6or4qllI+WUn6tlPJiKeXHLvjxby+l/Fwp5VdKKb9aSvmB6S+Vi9S0F8vjPA4XuMhrJ3fWvoRNcHNX/66N2FLKzSSfSfL9Sb4rySdLKd917mV/PsnPDMPw3Uk+keRvTn2h1M/NXQDAUsZMYr8nyYvDMHx5GIa3k/x0ko+fe82Q5Nse/ve7k/zH6S6RpJ692JqZfALAdoyJ2A8k+cqZb7/08PvO+otJfrCU8lKSzyX5k5NcXQdq3Ys91Jw3d23VjbftugLAvqb66vnJJD81DMMHk/xAkn9YSnnH2y6lfKqU8kIp5YWXX355onfNIda4uevYY7YAYB/2Yvs2JmK/muRDZ779wYffd9YPJ/mZJBmG4V8neTLJ+8+/oWEYPjsMw3PDMDz37LPPHnbFvEMvN3ddd8xWLcdwAfTk1RNPIaRNYyL2C0k+Ukr5jlLKEzm9cev5c6/59STflySllN+R04g1agUAYBbXRuwwDPeTfDrJ55P825yeQvDFUspPlFI+9vBlP5rkT5RS/p8k/zjJDw3D4N+OH7p0L3ZPh9zcxeM8tQsA+nBrzIuGYfhcTm/YOvt9P37mv7+U5PdMe2kb8Nbd5MkLzgt8863kqScneRc33nqQkydvTvK2xrj11pD7T/b/T1M37pac3PH3NIDavX5yN+++4WzeHrktGgCA5ojYDbvqhAIAgJqJ2ImdVH7O6xR6Pyv2htVjgK44aqtPIrYTUx+zNddZsQAAUxCxAAA0R8QuxDFbpzy1CwCYgohd21uX7OlsYLcWgDps4ald9mL7I2IB2LRXHjy99iUABxCxG3DjLUdpAQB9EbEAe3r6tnUfaJGVgr6IWGA2T91p+0ZEAOolYmewhQceAHC4107urH0J0DwRC4x2+7b9aqBtVgr6IWIBAGiOiAUANnFWLH0RsQCNefetaZ4A+J5bX5/k7UBrrBT0QcQCANAcEQvnnNwe8Zo7w3Tv74mTUa8b7ox7HeN8222TGICWiVigGk/ddq4ssAwrBe0TsQAANEfEAk17+gnTFOAwprFtE7EAC3rm9jQnCwBsnYhlUQ/urHsO4YMFn/R48sRy7+tQt+7cX/sSgIps8axY09h2iVjoTLntFAMA+idiAQBojojdgJMnb659CWzQU3cclwW0wUpBm0TsDG489eTalwCTu337wdqXUIWnb7+19iUAEBELAKt47WTBO025lmlse0QsF3rw5OW/NR48sc7dqw8auNt/X2MfOQsAPE7E1uqSlYSTp24vfCH7uf/k4YF7bKROcaTVyYhf3pM7w/HvaE/DHbELLGOLx2ztmMa2RcSu7Un/nARJ8tTt5W4E+7bbvlABtE7EduLkzv4T2gd3nFowlxYedNADj5wFpmYa2w4Ru5DyLU+tfQmr29LTutayz4MOPK0LgJaJWICGvPvWm2tfAnTPNLYNInZitZ0R60EHrMGDDoDWCdn6iVj24nitbfKgg2k8c9sUlTZs+YQC2iFia7Ty8VpXnRF7Fcdr7fn+Rp4R63gtgHWYxtZNxK5pxeO1nEwA++vpkbPvufX1tS8BmiBk6yViFzD3yQSHHK+1NUueTOB4rf0teUYsAH0QsSxi7eO1tmBrx2sdekasBx3M45mb31j7EmA2prF1ErEdczIBQN1eO9nAAdadELL1EbG1qfimrlpPJvDP9/Pa92QCx2tBH5xQQO1E7FomuqmrlsfNrnkywRScTADAdUxj6yJiZ+Zxs+vvw7qpC4CpCNl6iNhO2YcFgHkI2TqI2JpU9sjas9bah2Wclk8m2MrxWlM8revdtzzxC2ohZNcnYtew5z7slDd1XbUPW+OTutzUxUUOPV7rqPfZ0YMOeKdXHjy99iVUyc1dVxOy6xKxM/KQgza4qetyvZ1M4IxYYGpCdj0itkM17cO6qQuA3gnZdYjYpV22SuB8WABolpBdnoidSY2rBK2dD2sfdpyWb+oC6n9ql73Y8YTsskQss1l7lWCMWvdht2IrJxMA2yFklyNil7TnKsEhatqHXVvL+7Bu6rrcGicT9OY9t76+9iVA14TsMkTsDKZaJbhsH3bqVYJD92GPWSW4Tg2PoqUvx5xMcOzxWlOcEQu0RcjOT8SyimMj1T7sKfuwwBLsxR7m9ZO7YnZGInYpG1slWHsfdswqgX1YWuJpXdAuITsPETuxrawS8EiP+7BLcFMXsCWmstMTsRzE0Vp9q/WmLoDWidnpiNglLPCAA6sE+1l6lWAO++zD9sLJBLAee7HTErPHE7ENsUowzpJHa4219j7s1m/qOuZkAlhC7Q88YD5i9nC31r4Ajrf0FNYqwSP2Yft37PFaW/DMzW+sfQnQvLMh++4b/lIzhoid2wKrBIe4agp71NvtZJVgy/bdhz3Elm7qckYssC9BO46IbcTUqwRX/ryVHnAwhalWCWo/Wqum82Hd1AXb8epJyXtv1H2/QG+uWzXYRe7rJ3czzflI7RCxc9rYDV3XsUpwNasEV3NT1zQ8chb6suV9Wjd2NaCVG7qsElCjHm7q8qADgHcSsXPp8Alda97QNep9VLhKMIc5j9ayD/s4N3XBOzlqi1qI2EpM+YSuq7R6Q1erqwRj92HnWiWwDwtAr0TsHFaewrqh62JWCdrV6j6skwnYh7NiYT8itgJLTWHXssQqwVRqXyWY0xKrBEAfrBRQAxE7tQWmsFfp8Yau3lcJ9tH60VpL78P2cFMXABcTsStr4VitHm7oGrNKsNYU1tFacLn33nxjsff1yoOnF3tfwPFE7BIOmMI6VutUbVNYtqOWkwkcr0WtrBSwNhG7okN2YVuawtZiyhu61jyVoPWjtQ7V6k1dAMxLxM5tw7uwSzyhaytnw+7LPuz6nEzAIZxQAOOJ2JWYwrZnzRu62J+but7JI2eZmpUC1iRi51TBubCtTmHHXcP1r+nlhq6trhIAwGVE7ApaOBd27SmsG7oOV+MqwaHW2oet5aYuaIFpLGsRsXO5ZAo79ZFaprCXq/mGrq3Z2j7sljxz8xtrXwKwUSK2Eqawjyw5hbVK0OcqgX1YgP6J2DmYwh784+Ou4frXbHUKO/cqwZJaPlprqpMJnBG7TS2eUGClgDWI2KlV8GCDY5jCzq+WKewhHK0FQC1E7EKWPFLLFPb499Oifaewva4SHMtNXXAY01iWJmIXUMsawXVMYS94XUOrBJyyDwuwDSJ2RUvfzHXVFPa6gG1lCls7qwTjtLwPWysPOhjnlQdPr30JTTONZUkidma1TGGPWSNYwlRT2CkfbtDSFNYqQT08bpYptHhzFyxNxM7oqoBt6WYuU9hptDyF3Yqa9mGdTECrTGNZiohdwZZu5hrDFLZ+rZxKYB8WYDtE7ExqWSO4ztw3c5nCzq+3VQL7sJz1qh3VJpnGsgQRO4Op1wiOUfvNXFuZwra8SnDIFBY4nr1YuJqIXdChawRbvpnLubBXq3UK2+IDDqbYh3VTFzxiGsvcROzEpn6owVprBDXczDWlqaew+2h5Cru0Y1YJ7MMCbIuIXdmhAXudHm7mqnkKW8MNXftOYQ9hlWAdTiagF6axzEnELmCOPdhj1ghauJlrSqawp6wSbJcHHbTLXixcTsTOzBrBO5nCcpE1VwlqOh+W9Xlq17RMY5mLiJ3RGgG7lTWC0e+v0ynsEjd0WSXYj5u64HJCljmI2I5saY3AFLZOVgm25Zmb31j7EoANE7Ezqe04rbnXCMZYeo2glXNh99XjDV0ecACX62Uv1jSWqYnYGWxtD3bsa7hYLzd0raW3fVgnE9AzIcuUROzEWtuDnUKNN3OZwi7PKgEASxKxC1nrPNie1ghGv78Vb+baV41TWKsE+3NTF3PrZaUgMY1lOqMitpTy0VLKr5VSXiyl/Nglr/nDpZQvlVK+WEr5R9Ne5nbVvgfrZq7HbXUKuxZP6YI2CVmmcOu6F5RSbib5TJLfn+SlJF8opTw/DMOXzrzmI0n+bJLfMwzD10op/8VcF9yiFvdgx9jCGsGceprCtrpKUNs+7Bw86OAwrzx4Ou+7+cbalwFcYUwlfU+SF4dh+PIwDG8n+ekkHz/3mj+R5DPDMHwtSYZh+I1pL7Ndre7B1vRQg7HGBuxeb9MUdlY9rBLAUnpaKUhMYznemIj9QJKvnPn2Sw+/76zvTPKdpZR/VUr5xVLKRy96Q6WUT5VSXiilvPDyyy8fdsUNWStga9iDndLUawSmsNNaawpbyyrBlPuwTiZga4Qsx5jqxq5bST6S5HuTfDLJ3y2lPHP+RcMwfHYYhueGYXju2Wefnehd16nnG7m2skYw1xR27qdzHcoTuoA1CFkONSZiv5rkQ2e+/cGH33fWS0meH4bh3jAM/z7Jv8tp1G7SMQG75o1cY9S4RjCHWtYIDtHKubBrrxJsYR+W/vS2UrAjZDnEmIj9QpKPlFK+o5TyRJJPJHn+3Gv+WU6nsCmlvD+n6wVfnvA6mzFXwF5nioBdag92LDdz1T+F3foqATAdIcu+rq2mYRjuJ/l0ks8n+bdJfmYYhi+WUn6ilPKxhy/7fJJXSilfSvJzSf7MMAyvzHXRPar9Rq5J39cKawT7aHkK24q1p7BTcT5s31558PTalwBc4dojtpJkGIbPJfncue/78TP/PST5kYf/4wKt38jV+hrBVqawjtUaxyoBLXvt5E6eudHHXwTPe/Wk5L03ph9O0CdP7FqAgN1dy/WvSfq8mYvp9LpKMMfJBM6IpUXWChhLxM6s5oAdo/U92L3e94wBawr7SC+rBLCmXm/w2hGyjCFiZ1R7wC55I9dae7A1rBEAsD8hy3VE7ExaD9gptbAHawrbhilWCabah3VT1za4uWtdQpariNgZ9BCw9mCnMfeTuZLlz4W1SsDSXhWSF+p9pWBHyHIZETuxNQN2jC0E7Jzmvpmr9nNhmYfHzcLVXj0pYpZ3ELELWSJgazqJYKyxAbvX27RGMJs1p7A1rRJALbYyjd0RspwlYhewxYCdeg+2tTWCmpnCHs8+7LbYi62LkGVHxM5MwF7x9jpeIzCFfZxdWJjP1qaxiZDllIidUU8BO9aaAdvqGsGSWp3CWiUAzrMni4idSW8BW/uNXC2vEfQ+hWUdnta1LVucxu4I2e0SsTMQsFe8rYZu5EpMYadQww1dU5p6H9bJBG2wF1s3IbtNInZiNQTslNY6iWALe7CJKexSrBKwBVuexibWC7ZIxC5oqYCt9UausawRTK/VKSzAvoTsdojYhSy1QlBrwLa2B1vrGsHWprC1rRJA7bY+jd0xld0GEbuAWnZgx75OwPa1RrD1KeyUqwTOh11PDY+ebWUvVsg+Imb7dmvtC+jZdfGaCNhvvm6mgN3X3I+VTZZdIzjG2lPYLXBTFyxjF7LvvbH+/RZMxyR2JlME7IM7ZRMBu499A3buPdja1whaPRc2sUoAxzCNvZipbF9E7AymCtjr38/Y67n+NWsG7FwnEfSyRrCGY6ewPa4SQNLOSkEiZC9jxaAf1gkm1toJBEk7AbvFPdjEFLYW9mGhL2dD1ppBm0xiFyRgz7220YBd0hoB28sUdss8rYsd09hxTGfbJGIXImDPvbaSG7kO0fMaQU+sEjCXllYKEiG7DzHbFhG7AAF77rUzBqw1gvansFtZJXAywf5qOGaLbdjFrKCtm53YGU35CFkBO+b11ggALvLayZ08c6O+vyC2wO5svUxiZzLl9FXAjnn9/AHbwhqBKewpqwTMrbWVgsRawRRMaOsiYmew5PrAPq8TsKdqD1hTWGAuQnY6gnZ91gkmVmPA7hOFawfsvgTsI61PYWvlaC0u88qDp/O+m2+sfRlU4HzIWjtYhohdUM37r0kdATvnUVqHWOqJXMdq+UzYnS2tErip63CvPng67xWOR7EfO7+LprPCdnoidiE9rA8kbQfsUjdytbZGYAoL2yNkl3fV2oHAPYyIXYCAHfH6mc+CtUZQtxqP1YIxWl4pELL12GevVvA+ImJnNOX6wD6v6z1ge9uDPcaxAdvbFHbqVQL7sPROyLbnquD9LQteRw2cTjCTqY/PGnsDl4B93FIBe4xW1wimYgpL61o8bussJxbQKhE7sQd3StXrA4mAvcqhAbvFNYLaprBshyd3TU/I0iIRuwIBK2CTdQO2linslFo4lSBxMkGvWp/GJkKW9tiJXdDU8Zqsvz6Q9BGwh1pjD7YGU01ha14laHUf9j23vr72JdAwO7K0xCR2IXNMX1sL2H2u+dHPWSZgW9qDNYWFOvUwjU1MZGmHiJ3Z2Juydq8do5b1gTkfJXv6c5YJS2sE+6lxCtvKKkFvXnvwLau9b3ux8xKytEDEzmifeG0tYPexVMD2vgcL1K+XaWwiZKmfiJ3BXNPXfdYH9nmErIDdz1p7sKawy2l1H5Y6CFlYhoidWK/T10TAJscFbMvHadXMKsF2WSlYhpClViJ2BXNMX5O2A/bkiRMBOyNT2Do4Xms7eprGJqchK2apjSO2FjTH0VnJfvGazLs+kMx/AsHOVgK2hjWCWs0xhae6eWMAACAASURBVLVKAJdzBBc1MYldSA0BO/f+a9JnwK6pljUCU1hqVfNKQW/T2B0TWWphEjuzGuI1qW994PTnLLM+kBwXsC3vwfY8hYUWvPLg6bzv5htrX8bkdiFrKsuaTGJnsu+NWwL2elsK2ClMFbBTTWGn5oauOqx5VizrM5VlTSJ2BnNOX+deHxCwj7S8BzuVKQO2hVUC+7BtqnmlIOl3rWDHTV+sRcRObCvT10TAXqaGgLVGAHXpPWQTU1mWJ2JXUNP0NTlsfaDmI7SSdgN2CjWuEUw9hW1xlcDxWvOqfRqbbCdkxSxLEbELqnH6Wuv6QNJewE6hljUCgGMIWZYgYhdwSLzWNn1NBOx1elojqHkKOxf7sCxhC9PYHVNZ5iZiZ7Zv+NU6fRWwVxOwy2lxlaB3tZxQ0MJKQbKtkE3ELPMRsTPZ8vT19OcJ2LF6XiFoZQoLS9tayCZilumJ2InNHa/JMtPX5PDp66E3cG0xYKdiCgunWpnGJtsM2UTMMh0Ru5JD43Wp6WsL6wNJ+wFb0xrB1FqawtqHZS1bDdlEzHI8EbuCpeK15/WBRMBOrfYpLHWrZS82aWsam2w7ZBMxy+FE7IKWWB1Ilp++Ctj9TRWwW5rCtrxKsOQZsV+7/67F3hdMScyyr1trX8AW7BuuyXLxmmxn+poI2MuYwtKbVx88nffefGPtyxjtlQdP530NXe+czobsMzd8buJyInZGS8Vr0sbJA4mA3ULAtrQLm9iHpR67tQIx+8guaMUsF7FOMIOlbtpK2tl9TQQsh2t5lWBLatqLTdrbjd3Z+o7sRXarBtYNOEvETmzJ1YEld19Pf66APYYpLCxPyPZH0LJjnWBFS+69Ju3FayJg3/F2Kg7YuZjCslX2ZK9nf3bbTGJXcMzqwEHvr9Hpq4A993YqPYlgxxSWpL6VgqTdaWxiIrsPE9rtMYld0JI3bSWHh+vpz93u9DWpL2Cn1soUdm5u6tqO1k4rOMtEdn8mtNsgYhewdLwmAvYYNQbsVqewVgnglJA93PnJrKjth4id0ZbiNRGwZ9UcsKawzO21B9+SZ25+Y+3LeIeWp7GJI7imImr7IWJnIF73c2y8JgJ2rDkC1i5sfb52/115z62vr30ZVWo9ZBNT2aldtEMrbNsgYifWyokDj36+gE2cA1ujuVcJ7MPOq9ZpbC+E7LyEbRtE7Iq2HK9JnwFrCgv162Eam1gvWNpVpx4I3HWI2BWI1zriNRGwsFW9hGxiKluD6471ErnzELELOiZeEwG7I2DrMOcU1qkEfah9paC3kE1MZWt1yNm1wvd6InYB4nWaeE0E7CFMYeFyPYVsYirbEw9tuJ4nds3k5M7wzf8d/DaOeNLW6c8/6SZgb925L2APMFfAtj6F3cJNXV+7/661LyFJnU/wOq/lJ3pd5JUHT3vSF5tgEjuxtaeup2/j+GisJV6T6aavSb0BCzA1Kwb0ziS2EsdOXU/fxjST1x4D9vbtB1UHrCksPTONXZepLL0yiV1ZT5PXpL54TaY/A3bLATs3N3Sxpt72Y88ylaVHJrErqWXymgjYfbQQsHMyhW1LLXuxSRvT2KTviWxiX5a+iNgF7cK1pnidanVAwB7w9mYK2FansPRPyNZDzNID6wQLmGJl4PTtTBOKtU1ek2njNRGwc5h7CrvUKsEWTibgeD2vFpxlzYCWidgZ9RqvSb0BW3u8JgKWZX3t/rvynltfX/syvqn2ByCctZvIilmok4idWG3hmmwjXpM2AhZoz1amssnjJxkIWmpnJ7YyU+27JtPtvCbT7r0m2w1YU9iLOZVgW1rZjT1rC3uy59mbpXYithJbitep1wcErH/q52o1nVKwI2TbsYtZQUttROyKduFaY7wm064OJPVPX5P2AnZudmGZk5Btj5ilJnZiVzDlvmsy7c5rUn+8JgJ2p4cprFUCWrOlG74uY3eWGpjELmSuqevUk9eaVweSedYHEgF7EVPY6bx+/6m1L6HKlYKkzWnsztansjvWDViLSeyMpp64JtNPXZPpJ69JO9PXRMDC2lo6duu8LZ1cMIYJLUsSsTMQr9NqLV6Tdndgd5aawi69SrD1Bx3UdmbsWa2HbLLt9YKLCFrmJmInVuv5rme1Eq+JgL2IKSy9ajlkE1PZqwha5iBiKyRe54vXRMBexy5s32qexiZ9hGxiKnuV87uzopZDidhKtBSuOwL23NvuIGCBaYjZ8UQthxKxK5orXBPxet7cj5BtfQd2Z8kp7NaO1nr9/lN59606dnJNY5djxWB/opaxROzC5gzXpM14TQTsdUxh2ZreQjYxlT2UqOUyInYBLYdr0m68JgJ2H3Zht6X2aWzSV8gmYnYqF51HK2y3ScTOZO5wTdqO10TAjtHrBHZrqwS1aiVkk4hZrnTZgxbEbd9E7MRan7om4nXU++hkB3bHFHYZNe3Ftqa3qWwiZpcgbvsmYhvQQ7gm88dr0lfA9jqFpS4tTGN3egzZRMyu4apH5ArcdojYionX8ZaI16TPgDWFRcjWQczW4arATURuTURsZZYI12SZeE0E7CFMYPtmpeB4PYdsImZrd13knid65yNiK7BUuCZ9xWsiYI+1xhTWTV11amkam/R5w9d5YrYP+0Yv44nYFSwZrcly4Zr0F69Jfzdxsb5ap7GthWzS/1Q2eRSziaCFs26MeVEp5aOllF8rpbxYSvmxK173B0opQynluekusQ/l9sk3/7eUW3fuLzp57S1gn7p9b9GA3cIUdm2v3Xtq7Uuo3tfuv2vtS9jbbiq7Ba8+ePqxqIUtu3YSW0q5meQzSX5/kpeSfKGU8vwwDF8697pvTfKnkvzSHBfamqWnrTtLTl2T5SavSd/TV3uwcJwtrBecZToL4yax35PkxWEYvjwMw9tJfjrJxy943V9K8peTbHrhbelpa/Jo4rr02oCAncYaAbvFKWxtXr9f71S4xWnszpamsjums2zVmIj9QJKvnPn2Sw+/75tKKb8ryYeGYfjnV72hUsqnSikvlFJeePnll/e+WB5ZI1yTdeJVwPbDTV3taD1ktxyzgpatGLUTe5VSyo0kfy3Jj1732mEYPjsMw3PDMDz37LPPHvuuN2ftcO05XpPt3MBlCluPmqexSdshm2xzKrsjaNmCMacTfDXJh858+4MPv2/nW5P8ziQ/X0pJkt+S5PlSyseGYXhhqgvdqqWD9awlo/WspeM1WSdgtzaFpU0tnlhw1tZ2ZS9if5ZejYnYLyT5SCnlO3Iar59I8kd2PzgMw+tJ3r/7dinl55P8bwL2cGuGayJel7BWwJrC1qfW47bOaj1kEzG7I2jpybUROwzD/VLKp5N8PsnNJD85DMMXSyk/keSFYRien/sie7d2tCbrhWuyTrwm2wtYTr1276k8c7uuaBSyy9nCubJjnV81ELW0ZtTDDoZh+FySz537vh+/5LXfe/xl9a2GaN3ZYrwm2wxYU1iO1VPIJqay55nS0hpP7FpATdGarBuuO1sMWLhKC9PYpJ+QTcTsVUxpaYGInVhtwbpTQ7gm245XawT0YndqgZjdDlFLjURsx2oJ12TdeE0ErFWCx9W4F5u0M43d6Wkqm4jZfVx0dJewZWkitjM1hWuyfrwmArYWHnQwjpBdn5g9jGktSxOxHagtXBPxulNDwJrCtkfI1kHMHse0lrmJ2EbVGK5JHfGaCFiuV+tKwU6LIZv0syd71tknfwna41z2BDFxyyFEbCNqjdaknnDdqSFga2EKy9J6ncrumM7OQ9xyCBFbqZqjdUe8Xs4Utg2msfPoeSq7Yzq7jMviNhG4iNgqtBCsZ9UWr4mAbcEb9550c9cBWg3ZpP+p7I7p7DoELiJ2Ya0F606N4ZrUFa9JXQFrlWCc2qexSfshm/Q9ld0xna3HVYG7I3TbJ2Jn1Gqw7tQarjsCli1pOWSTbcVsImhbMCZ0dwRvnUTsxFoP10S8HkLAtq2FaWzSfsgm21kxOEvQtm+f4L2MEJ6eiCVJ/eG6I2DHsUrQr15CNtnOVPYsQbtdU4QwjxOxG9ZKuCZ1xit9aWUam/QRssm2YzYRtHAsEbshLUXrTu3xWuMUlsO1FrJJxGwnBC3sT8R2rsVw3RGwcL1eprKJmN05G7SJqIXLiNjOtBytO7XHa1J3wNa8D9vCWbEtTWN3egrZRMyeJ2rhYiK2cT1E604L8ZrUHbBMo9WQTfpYL9gRsxcTtXBKxDakp2A9q5V4TQTslrQYskl/U9lEzF5H1LJVIrZCvcbqeS3FayJgp9LCSkHregzZ5FHMJoL2KqKWrRCxK9lKqF6ktXhlu1qdxiZ9rhecZTo73vmoTYQtfRCxM9hyoF6l5XhtZQpb801drWo5ZJPtxGwiaPchbOmBiJ2YgH2nluM1aSdgW9LaSkHrIZv0u2JwlunscYQtrRGxzKb1eE0ELI/0ErJJv1PZHdPZ6VwUtom4pQ4ilkn1EK47AnZerU1jkz5CNtlOzCaCdi6XxW0icFmOiGUSPcVrImC5XC8hm2wrZhNBuxSBy1JELAfrLVxZXovT2OQ0ZJN0FbNbCdkdQbsOgcuURCx72UK4msIuq9WQTUxleyFo63BV4O4IXc4SsVxrC+G6I2DZV08hmzyK2UTQJqK2NmNCNxG7WyFieYctRetZAnY9LU9jk/7WC3a2PJ3dMaVt09jYPUv4tkfEkmS74bojYNfXesgm/U1ld8TsKVPavh0SvpcRxMsQsRu09WA9T8DWQ8jWTcw+TtRymSmDmMuJ2M4J1m35zXt3PHq2Ar2uF+xsfW/2MqIWliViOyBUD2cKW58eprE7vcdsYjp7lfNRmwhbmJKIrZAoXYaArVdPIZv0vWKwYzo7jmktTEfETkyAtkHA1q/HkE36nsruCNrxTGvhcCIWOtPTXmxvIZtsK2YTQXsIYQvjiFg2xxS2LW/cezJJxGwHBO3hLgrbRNyybSKWTdlKwPY0jd3pcSqbbDNmE0E7FXHLlolYoBm9hmyy3ZhNBO0cxC1bIGLZjK1MYXd6nMYm/a4X7Gw5ZpPHgzYRtVO7LG4TgUt7RCybsLWA3YKep7LJo5hNthu0iSntkq4K3ETkUh8RCx3rdRq70/tUdmfr09kdU9p1iVxqI2Lp3tansL2HbNL/VHZHzD5O1NblushNhC7TErGwAVsJ2aT/qWxi1eAyorZ+Y0J3R/ByHRFL17Y+hT1rCyGbbCtmE0F7FVHbtn2Cd0f4bouIhQ3ZSsgm21kxOEvQXu181CbCtjeHhO9ZIrgtIhY2Zmshm2xnKnuWoB3HtJazjo3gQwjnw4lY2KAthWyy7ZhNBO0+TGtZ2hrh3AsRS7fsw15tayGbPIrZRNAmgnasi8I2EbewNhELG/ab9+4kyeZiNjGdTQTtsUxtYV0iFtjkVHZHzJ46G7SJqD2UqS0sR8QCSbY9lU2sGpxnSjuty+I2EbhwKBFLt954+4692ANseSq7Yzr7OFPaeQlcOIyIndjTT9zNG2/fWfsy4Chbn8rumM5eTNQuR+DC5UTsDM5O/wQtLROzjwjay4nadVwVuInIpX8idma7oBWz67BSMA0x+zhBe7XzUZsI2zWIXHonYhciZtcjZKcjZt9J0I5jWluf6yJ3R+xSKxG7MKsG6xCy0xKzFxO045nWtmNs7CaCl2WJ2BW5CWxZQnZ6YvZygnZ/wrZ9+wRvIno5johdmTWDZQnZeexiNhG0FzkbtImo3Yew7du+0XueCN42EVsJMbscITsv09nrmdIe56KwTcTtFh0bweeJ4raI2MpYMVjG7tdYzM7HdHYcU9rpiFuONXUUMy8RWyFT2eWI2WUI2vFE7fTELfRJxFbMVHY5YnY5gnY/onY+4hbaJmLhDDG7LEG7P1E7v8viNhG4UBMRWznT2HWc/TUXtMsQtIc5H7WJsJ3TVYGbiFxYkohtgJBdl6Bd3tmgTUTtvkxr1yNyYTkiFvYgaNdhSnsc09p6XBe5idCFsUQsHEjQrsOUdhrCtl5CF8YRsTCB8+seonY5onY6F4VtIm5rNCZ0E7FL30QszMCUdj2idnqmtu0aG7s7opeWiFiYmSntukTtPExt+7Rv9CbCl/WIWFiYqF2XqJ3XZXGbCNxeHRK+Z4lgDiViYWWidl3nozYRtnMxveUix0bwWYJ4W0RsA5wRuy2idn3Cdlmmt0xlyiA+hIheloiFyl30lxhhuzxhuw6BS0vWjuitEbGVM4XlIsK2DsJ2XVcFbiJyoXciFjphDaEOwrYeIhf6JmIrZgrLMUxr63FR2Cbidm3XRW4idKFmIrZSApY5XPb7StyuQ9zWb0zoJmIX1iBiKyRgWZqpbV3EbXvELixPxFZGwFILYVsfcdu+sbGbCF64joitiICldtYR6nRZ3CYCt2X7BO+O8GVLRGwlBCwtE7f1Erjbckj4JuKXNonYCghYeiVu6yZw2Tk0fndEMGsQsSsTsGyRuK3fVYGbiFwed2wEX0QYcx0RuxLxCu8kbtshcpnbHGF8FdHcHhG7AgEL+7nqz4zArZPIpTVLRzPHE7ELE7AwLYHbpusid0fsApcRsQsSsLAsgdu+MbErdGGbROwCxCvUR+D2Y+xUNxG80BMROzMBC+0RuP0SvNAPETsT8Qp9uu7Ptsjtxz7Bm4heWJqInZh4hW0Tudu1b/QmwheOIWJpwpv3bq99CY956va9tS+BRo35i67Q3Y5DwndHALN1IpZV1Bal+5rq+sUwFzHNZYxjAnhHCNMyEcvkWg/UJR3yayV8Mc1lKlOE8HnCmKWIWA4mVtex76+76N2msfv5YpepzRHGlxHM2zYqYkspH03yfya5meT/Gobhfz/34z+S5H9Kcj/Jy0n+x2EY/sPE18pKxGrb9vn4Cd7tEbu0bMlgXsvYUN/Cr8V510ZsKeVmks8k+f1JXkryhVLK88MwfOnMy34lyXPDMHyjlPI/J/krSf6HOS6Y+YnW7Rr7sRe727PPySuCF6azxTgda8wk9nuSvDgMw5eTpJTy00k+nuSbETsMw8+def0vJvnBKS+S+QhWDmG6y1UEL7CEMRH7gSRfOfPtl5L8t1e8/oeT/ItjLor5iFaWZrrLVfY9W1v0AjuT3thVSvnBJM8l+b2X/PinknwqSb792799ynfNJUQrrTDdZYxDHigjfKFPYyL2q0k+dObbH3z4fY8ppfy+JH8uye8dhuHCzxjDMHw2yWeT5Lnnnhv2vlquJVrZAtNd9nHokxTFL9RtTMR+IclHSinfkdN4/USSP3L2BaWU707yd5J8dBiG35j8KrmScIWLme5yjGMeIy6AYX7XRuwwDPdLKZ9O8vmcHrH1k8MwfLGU8hNJXhiG4fkkfzXJ00n+aSklSX59GIaPzXjdmyZaYXqmu0zpmAA+SwzD5UbtxA7D8Lkknzv3fT9+5r9/38TXxTnCFepgusuSporhHVFMTzyxq2LCFdrm6WrUZuooPk8ksyQRWxHRCtsmemnd3JE8lpjeBhG7MuEKHEr0wsVqiWnmJWJXIFyBNRzyuUf4ArUSsQsRrkCLDv3cJX6BuYnYGQnXebx516/rZZ66IxyowzGf/wQwMIaInZhwHU+MTm/qX1NRzBoEMDCGiGUWArUPU3wchTBLmmqQIIahfiKWvQlU9nHs7xcRzBqm/Fc1QQzzELFcSKhSCxFM6+ZaMxPHbJ2I3TixSu+O+T0ugKnZkvdgCGZqJGI3QqzC/gQwnKrxpmVhjYjtkGCF9QlgmFeNYc2yRGzDxCr0yR4wwPVEbCMEKzCWCAa2QMRWSrQCa5nq848YBuYkYisgWIEeTfm5TRAD54nYFYhWgP14pDJwnoidmWCt1717N9e+hNncvv1g7UuAqs35uVkgwzJE7MRE6/R6js25rPlrJqDZuiW/DghmtkzEshgxug1zfpwFMjyulsGJmGYNIpajCFOWNMfvN2EMx6slpuck1OsjYrmUQGULpv59LoqhT1sI9daI2A0TqTA902KAZYjYjolU6INpMcA7idjGCVVgX3N93hDHwJJEbCPEKlC7uT9PiWTgLBFbEaEKcLklPkcKZWiHiF2BWAWo01qfn8Uz7E/EzkisAjBGrV8vxDU1E7ETq/UT0Rbdv+u395Ju3bm/9iUAE9vi1zTh3g5f5amC4GxfDR9DIQ0ca4vh3qr1v+rQlRpChu1a8vefYAZYl+LgSqIULiaYAdalUDZImEJblvozK5aBlqiZzghU4FBLf/4QzcAxFE9DBCrQkzU/pwloaJ8qqohIBVhGbZ9vRTXsr64/xZ2r7ZMmAHVo/euDCGcNbf+pqVDrn4gAYF++9rEGv+vgEsO9G2tfQnXK7ZO1LwEAkohYKick69Lqx0N8A/RHxDK5VkOHftX0e1JQA0xDxPIONX3Bh97U8udLTAOtE7Gdq+ULJlCXmj43CGrgECK2QTV98QE4Vs2f0wQ21EvEVqLmT+IAW9XD52YhTq9E7Mx6+AQIQLt8HVrG2n9Z2OLHWcRObIu/iQBg63z9X55fcQAAmmMSS1fKXX8v25rhjn0/gC0SscxOWDKnln9/CXCAw4lY3qHlKICWtPZnTXQDNRGxnWntiyLQjlY/v4hv6JOIrVSrXywAatP751ORzlaJ2Jn1/skTgHX5OsNWidiJ+WQCADA/xQUAQHNELAAAzbFOQHNuvO3vXszr5Ak3ygDUTsRyECFJz7bw+1uoA60TsR3bwhdi4DBb+Pwg1KFvIrYiW/iiArAUn1PHE/y0SMROzCdNAFrjaxct8rsWAIDmiFgAAJojYgEAaI6dWDbjxttrXwFLOnli7SsAYE4iloMIQmrn9+j+hD/QEhHbMF+kgSn5nLI/4Q/rEbET80UAYDt8zof1uLELAIDmiFgAAJojYgEAaI6IBQCgOW7sohs37pa1L4GOnNwZ1r4EAK4gYjdK8MHV/Bmpj79YAGeJ2In5wgcwD59fmUsPf0Ha4p8PEQsAbNoWA7AHbuwCAKA5IhYAgOaIWAAAmiNiAQBojhu7GO3GvbWvAKjJye21rwDYMhE7MaEHbIXPd8CarBMAANAcEQsAQHNELAAAzRGxAAA0R8QCANAcpxNwsJt3174CYAkP7qx9BQDvJGInJuyA3vi8BtTIOgEAAM0RsQAANEfEAgDQHBELAEBzRCwAAM1xOsHEbry99hUA0KKTJ9a+AmiLiAWAChiCwH6sEwAA0JxREVtK+Wgp5ddKKS+WUn7sgh+/U0r5Jw9//JdKKR+e+kIBAGDn2ogtpdxM8pkk35/ku5J8spTyXede9sNJvjYMw29P8teT/OWpLxQAAHbGTGK/J8mLwzB8eRiGt5P8dJKPn3vNx5P8/Yf//bNJvq+UUqa7TAAAeGRMxH4gyVfOfPulh9934WuGYbif5PUk7zv/hkopnyqlvFBKeeHll18+7IoBANi8RU8nGIbhs0k+myTPPffcsOT7Xsr/+9f+9NqXAADQvTGT2K8m+dCZb3/w4fdd+JpSyq0k707yyhQXCAAA542J2C8k+Ugp5TtKKU8k+USS58+95vkkf+zhf//BJP9yGIYuJ60AAKzv2nWCYRjul1I+neTzSW4m+clhGL5YSvmJJC8Mw/B8kr+X5B+WUl5M8mpOQxcAAGYxaid2GIbPJfncue/78TP//VaSPzTtpQEAwMU8sQsAgOaIWAAAmiNiAQBojogFAKA5IhYAgOaIWAAAmiNiAQBojogFAKA5IhYAgOaIWAAAmiNiAQBojogFAKA5IhYAgOaIWAAAmiNiAQBojogFAKA5IhYAgOaIWAAAmiNiAQBojogFAKA5IhYAgOaIWAAAmiNiAQBojogFAKA5IhYAgOaIWAAAmiNiAQBojogFAKA5IhYAgOaIWAAAmiNiAQBojogFAKA5IhYAgOaIWAAAmiNiAQBojogFAKA5IhYAgOaIWAAAmiNiAQBojogFAKA5IhYAgOaIWAAAmiNiAQBojogFAKA5IhYAgOaIWAAAmiNiAQBojogFAKA5IhYAgOaUYRjWecelvJzkP6zyzuf1/iT/ee2LYC8+Zu3xMWuPj1lbfLza0+vH7LcNw/DsRT+wWsT2qpTywjAMz619HYznY9YeH7P2+Ji1xcerPVv8mFknAACgOSIWAIDmiNjpfXbtC2BvPmbt8TFrj49ZW3y82rO5j5mdWAAAmmMSCwBAc0QsAADNEbEHKqV8tJTya6WUF0spP3bBj98ppfyThz/+S6WUDy9/lZw14mP2I6WUL5VSfrWU8n+XUn7bGtfJI9d9zM687g+UUoZSyqaOl6nNmI9XKeUPP/xz9sVSyj9a+hp53IjPi99eSvm5UsqvPPzc+ANrXCenSik/WUr5jVLKv7nkx0sp5W88/Hj+ainldy19jUsSsQcopdxM8pkk35/ku5J8spTyXede9sNJvjYMw29P8teT/OVlr5KzRn7MfiXJc8Mw/NdJfjbJX1n2Kjlr5McspZRvTfKnkvzSslfIWWM+XqWUjyT5s0l+zzAM/2WS/3XxC+WbRv4Z+/NJfmYYhu9O8okkf3PZq+Scn0ry0St+/PuTfOTh/z6V5G8tcE2rEbGH+Z4kLw7D8OVhGN5O8tNJPn7uNR9P8vcf/vfPJvm+UkpZ8Bp53LUfs2EYfm4Yhm88/OYvJvngwtfI48b8OUuSv5TTvyS+teTF8Q5jPl5/IslnhmH4WpIMw/AbC18jjxvzMRuSfNvD/353kv+44PVxzjAMv5Dk1Ste8vEk/2A49YtJniml/NZlrm55IvYwH0jylTPffunh9134mmEY7id5Pcn7Frk6LjLmY3bWDyf5F7NeEde59mP28J/KPjQMwz9f8sK40Jg/Y9+Z5DtLKf+qlPKLpZSrJkrMb8zH7C8m+cFSyktJPpfkTy5zaRxo3691Tbu19gVAbUopP5jkuSS/d+1r4XKllBtJ/lqSH1r5UhjvVk7/mfN7c/ovHb9QSvmvhmF4bdWr4iqfTPJTwzD8H6WU/y7JPyylbieYRwAAAddJREFU/M5hGE7WvjAwiT3MV5N86My3P/jw+y58TSnlVk7/GeaVRa6Oi4z5mKWU8vuS/LkkHxuG4e5C18bFrvuYfWuS35nk50sp/1+S353keTd3rWbMn7GXkjw/DMO9YRj+fZJ/l9OoZR1jPmY/nORnkmQYhn+d5Mkk71/k6jjEqK91vRCxh/lCko+UUr6jlPJETpfdnz/3mueT/LGH//0Hk/zLwZMl1nTtx6yU8t1J/k5OA9au3vqu/JgNw/D6MAzvH4bhw8MwfDine8wfG4bhhXUud/PGfF78ZzmdwqaU8v6crhd8ecmL5DFjPma/nuT7kqSU8jtyGrEvL3qV7OP5JH/04SkFvzvJ68Mw/Ke1L2ou1gkOMAzD/VLKp5N8PsnNJD85DMMXSyk/keSFYRieT/L3cvrPLi/mdAn7E+tdMSM/Zn81ydNJ/unDe/B+fRiGj6120Rs38mNGJUZ+vD6f5L8vpXwpyYMkf2YYBv9CtZKRH7MfTfJ3Syl/Oqc3ef2Qgcx6Sin/OKd/EXz/wz3lv5DkdpIMw/C3c7q3/ANJXkzyjSR/fJ0rXYbHzgIA0BzrBAAANEfEAgDQHBELAEBzRCwAAM0RsQAANEfEAgDQHBELAEBz/n9W3m8Mg/ytxwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 864x864 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": false,
        "id": "Kw6PDrRjSuGy"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}