Part2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ffc29f77",
   "metadata": {},
   "source": [
    "## First let's do some symbolic derivation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "de1dec6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Math,display\n",
    "import sympy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3eebb9da",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}x_{0} & y_{0}\\\\x_{1} & y_{1}\\\\x_{2} & y_{2}\\\\x_{3} & y_{3}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[x0, y0],\n",
       "[x1, y1],\n",
       "[x2, y2],\n",
       "[x3, y3]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xs=x0,x1,x2,x3=sympy.symbols('x:4')\n",
    "ys=y0,y1,y2,y3=sympy.symbols('y:4')\n",
    "xy=sympy.Matrix([xs,ys]).T\n",
    "display(xy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f73c0051",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\left(1 - \\eta\\right) \\left(1 - \\xi\\right)}{4} & \\frac{\\left(1 - \\eta\\right) \\left(\\xi + 1\\right)}{4} & \\frac{\\left(\\eta + 1\\right) \\left(\\xi + 1\\right)}{4} & \\frac{\\left(1 - \\xi\\right) \\left(\\eta + 1\\right)}{4}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[(1 - eta)*(1 - xi)/4, (1 - eta)*(xi + 1)/4, (eta + 1)*(xi + 1)/4, (1 - xi)*(eta + 1)/4]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\eta}{4} - \\frac{1}{4} & \\frac{1}{4} - \\frac{\\eta}{4} & \\frac{\\eta}{4} + \\frac{1}{4} & - \\frac{\\eta}{4} - \\frac{1}{4}\\\\\\frac{\\xi}{4} - \\frac{1}{4} & - \\frac{\\xi}{4} - \\frac{1}{4} & \\frac{\\xi}{4} + \\frac{1}{4} & \\frac{1}{4} - \\frac{\\xi}{4}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[eta/4 - 1/4, 1/4 - eta/4, eta/4 + 1/4, -eta/4 - 1/4],\n",
       "[ xi/4 - 1/4, -xi/4 - 1/4,  xi/4 + 1/4,   1/4 - xi/4]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# corrdinate of parent domain\n",
    "xi,eta=sympy.symbols('xi eta')\n",
    "# shape functions\n",
    "N=sympy.Rational(1,4)*sympy.Matrix([(1-xi)*(1-eta),(1+xi)*(1-eta),(1+xi)*(1+eta),(1-xi)*(1+eta)]).T\n",
    "display(N)\n",
    "# partial(N1,N2,N3,N4)/partial(xi,eta)\n",
    "GN=sympy.Matrix([[sympy.diff(i,xi) for i in N],[sympy.diff(i,eta) for i in N]])\n",
    "display(GN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9a86cfa8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{x_{0} \\left(\\eta - 1\\right)}{4} - \\frac{x_{1} \\left(\\eta - 1\\right)}{4} + \\frac{x_{2} \\left(\\eta + 1\\right)}{4} - \\frac{x_{3} \\left(\\eta + 1\\right)}{4} & \\frac{y_{0} \\left(\\eta - 1\\right)}{4} - \\frac{y_{1} \\left(\\eta - 1\\right)}{4} + \\frac{y_{2} \\left(\\eta + 1\\right)}{4} - \\frac{y_{3} \\left(\\eta + 1\\right)}{4}\\\\\\frac{x_{0} \\left(\\xi - 1\\right)}{4} - \\frac{x_{1} \\left(\\xi + 1\\right)}{4} + \\frac{x_{2} \\left(\\xi + 1\\right)}{4} - \\frac{x_{3} \\left(\\xi - 1\\right)}{4} & \\frac{y_{0} \\left(\\xi - 1\\right)}{4} - \\frac{y_{1} \\left(\\xi + 1\\right)}{4} + \\frac{y_{2} \\left(\\xi + 1\\right)}{4} - \\frac{y_{3} \\left(\\xi - 1\\right)}{4}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[x0*(eta - 1)/4 - x1*(eta - 1)/4 + x2*(eta + 1)/4 - x3*(eta + 1)/4, y0*(eta - 1)/4 - y1*(eta - 1)/4 + y2*(eta + 1)/4 - y3*(eta + 1)/4],\n",
       "[    x0*(xi - 1)/4 - x1*(xi + 1)/4 + x2*(xi + 1)/4 - x3*(xi - 1)/4,     y0*(xi - 1)/4 - y1*(xi + 1)/4 + y2*(xi + 1)/4 - y3*(xi - 1)/4]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# J is partial(x,y)/partial(xi,eta)\n",
    "J=sympy.simplify(GN*xy)\n",
    "display(J)\n",
    "Jdet=sympy.simplify(J.det())\n",
    "Jinv=sympy.simplify(J**(-1))\n",
    "# partial(N1,N2,N3,N4)/partial(x,y)\n",
    "NablaN=sympy.simplify(Jinv*GN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ab8d0d2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "Ksym=[]\n",
    "for i in range(4):\n",
    "    Ksym.append([])\n",
    "    for j in range(0,i+1):\n",
    "        integrand=sympy.simplify(Jdet*(NablaN[:,i].T*NablaN[:,j])[0])\n",
    "        Ksym[-1].append(integrand)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85e17ca8",
   "metadata": {},
   "source": [
    "## Then I read a mesh exported from blender"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "24763630",
   "metadata": {},
   "outputs": [],
   "source": [
    "#%matplotlib notebook\n",
    "import meshio,numpy\n",
    "meshio.__version__\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cf826178",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<meshio mesh object>\n",
      "  Number of points: 36\n",
      "  Number of cells:\n",
      "    quad: 25\n",
      "  Point data: nx, ny, nz, s, t\n"
     ]
    }
   ],
   "source": [
    "sample=meshio.read(\"square5x5.ply\")\n",
    "#sample=meshio.read(\"square20x20.ply\")\n",
    "print(sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b02d9586",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1.         -1.        ]\n",
      " [-0.6        -1.        ]\n",
      " [-0.6        -0.6       ]\n",
      " [-1.         -0.6       ]\n",
      " [-0.19999999 -1.        ]\n",
      " [-0.19999999 -0.6       ]\n",
      " [ 0.20000005 -1.        ]\n",
      " [ 0.20000005 -0.6       ]\n",
      " [ 0.6        -1.        ]\n",
      " [ 0.6        -0.6       ]\n",
      " [ 1.         -1.        ]\n",
      " [ 1.         -0.6       ]\n",
      " [-0.6        -0.19999999]\n",
      " [-1.         -0.19999999]\n",
      " [-0.19999999 -0.19999999]\n",
      " [ 0.20000005 -0.19999999]\n",
      " [ 0.6        -0.19999999]\n",
      " [ 1.         -0.19999999]\n",
      " [-0.6         0.20000005]\n",
      " [-1.          0.20000005]\n",
      " [-0.19999999  0.20000005]\n",
      " [ 0.20000005  0.20000005]\n",
      " [ 0.6         0.20000005]\n",
      " [ 1.          0.20000005]\n",
      " [-0.6         0.6       ]\n",
      " [-1.          0.6       ]\n",
      " [-0.19999999  0.6       ]\n",
      " [ 0.20000005  0.6       ]\n",
      " [ 0.6         0.6       ]\n",
      " [ 1.          0.6       ]\n",
      " [-0.6         1.        ]\n",
      " [-1.          1.        ]\n",
      " [-0.19999999  1.        ]\n",
      " [ 0.20000005  1.        ]\n",
      " [ 0.6         1.        ]\n",
      " [ 1.          1.        ]]\n",
      "[[ 0  1  2  3]\n",
      " [ 1  4  5  2]\n",
      " [ 4  6  7  5]\n",
      " [ 6  8  9  7]\n",
      " [ 8 10 11  9]\n",
      " [ 3  2 12 13]\n",
      " [ 2  5 14 12]\n",
      " [ 5  7 15 14]\n",
      " [ 7  9 16 15]\n",
      " [ 9 11 17 16]\n",
      " [13 12 18 19]\n",
      " [12 14 20 18]\n",
      " [14 15 21 20]\n",
      " [15 16 22 21]\n",
      " [16 17 23 22]\n",
      " [19 18 24 25]\n",
      " [18 20 26 24]\n",
      " [20 21 27 26]\n",
      " [21 22 28 27]\n",
      " [22 23 29 28]\n",
      " [25 24 30 31]\n",
      " [24 26 32 30]\n",
      " [26 27 33 32]\n",
      " [27 28 34 33]\n",
      " [28 29 35 34]]\n"
     ]
    }
   ],
   "source": [
    "points=sample.__dict__['points'][:,0:2]\n",
    "cells=sample.cells_dict['quad']\n",
    "Npts=len(points)\n",
    "print(points)\n",
    "print(cells)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6a51fa23",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD4CAYAAAAQP7oXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0iklEQVR4nO3de3gUBZb38e8xQEIuQDDRBKMJFwUMII7t4KwKXnAE1lUGBUFFRGZZBNydcX1HVBAGnR15dXVU1NXH1XG8IYIiCorAyGV0UMMYMIBAEFBiDCC3cAmQ5Lx/dCVvE9NJJV3dBeF8nqeedF3zO9WdPl1VnW5RVYwxxhg3TvE7gDHGmBOHNQ1jjDGuWdMwxhjjmjUNY4wxrlnTMMYY41ozvwM0Rlpamubk5PgdwxhjTigrV67cqarpkWzjhGwaOTk55OXl+R3DGGNOKCKyNdJt2Ompetxyyy1kZmbSqlUrzjnnHF544QUAjhw5wg033EBOTg4iwpIlS/wN6qFwNQMsXryYLl26kJiYyOWXX87WrRE/Bo8L4WpesWIFV111FW3btiU9PZ3BgwdTXFzsc1pvhKt57dq1BAIBUlNTSU1NpW/fvqxdu9bntN6o67FdZerUqYgIixYt8iGht2rWC6QBiEiOiKiI7A8ZJrnaqKqecMMFF1ygsVJQUKBlZWWqqrpu3To9/fTTNS8vTw8fPqyPP/64Ll++XDMyMvTjjz+OWaZoC1fzjh07tFWrVjpz5kw9dOiQ3n333dqrVy+f03ojXM3z58/XmTNn6t69e/XAgQM6cuRIvfrqq31O641wNe/evVs3b96slZWVWl5erk888YR2797d57TeCFdzlcLCQu3WrZtmZmbqwoUL/YrpmZr1AkeBC4AcQIFm2sDnXzvSqEdubi7x8fEAiAgiwqZNm2jRogW/+c1vuOSSS4iLi/M5pbfC1fz222+Tm5vL4MGDSUhIYMqUKaxatYqvv/7a58SRC1dz//79GTx4MK1atSIxMZHx48fzySef+JzWG+FqbtOmTfURtKoSFxdHYWGhz2m9Ea7mKuPGjWPatGm0aNHCr4ieqlkvwUbRMZJtWtNwYezYsSQmJtKlSxcyMzMZMGCA35Girraa16xZw3nnnVe9TFJSEh07dmTNmjU+JvWOm/t52bJl5Obm+pAuOuqquU2bNiQkJHDnnXdy3333+ZjSW+Fqfuutt4iPj29yf9+h9RI80pgfMnuriGwTkZdEJM3N9qxpuPDMM89QWlrK8uXLGTRoUHXnbspqq3n//v20bt36mOVat25NaWmpTym9Vd/9vHr1aqZOncojjzziU0Lv1VXznj172Lt3L9OnT+f888/3MaW3aqu5tLSU++67jyeeeMLveJ4LrRfYAxwGdgIXAtkET1elAK+52Z4nTUNEXhSR7SJSEGa+iMiTIlIoIqtF5Gch80aIyEZnGOFFnmiIi4vjkksuYdu2bTz77LN+x4mJmjUnJyezb9++Y5bZt28fKSkpPiX0Xrj7ubCwkP79+/PEE09w6aWX+pjQe3U9tpOSkhgzZgy33nor27dv9ymh92rWPGXKFIYPH05TfSt/Vb1Ac+AOVd2vqnmqWq6qJcB44JciUu8fs1dHGn8G+tUxvz9wtjOMBp4FEJG2wGSgF/BzYLKIpHqUKSrKy8uPOQd6MqiqOTc3l1WrVlVPP3DgQPX0pib0ft66dSt9+/Zl0qRJDB8+3Odk0RPusV1ZWcnBgwcpKiryIVV0VdW8ePFinnzySTIyMsjIyOC7775jyJAhTJs2ze+IXhNqv6ZR9XHn9feEhl45DzcQvBpfEGbec8CwkPH1QCYwDHgu3HLhhsa+e2riO19phwnzdOI7X7lavqSkRN944w0tLS3V8vJy/fDDDzUxMVHfffddVVUtKyvTQ4cO6RlnnKELFizQQ4cOaWVlZaOyRYuXNW/fvl1btWqls2bN0kOHDunvfve74/LdU17WvG3bNu3QoYM+8sgjUU4dGS9r/uijj/Qf//iHlpeX6969e/XOO+/UzMxMPXToUJSraBgva965c6cWFxdXD1lZWTpz5kwtLS2NchXueVEvUAFcS/CFemenSZwKvAl8rG6e690s5GpDdTeN94FLQsYXAwHgbmBiyPRJwN1htjEayAPyzjrrrEbt9A4T5mn2Pe9rhwnzXC2/fft27d27t7Zu3VpTUlK0W7du+vzzz1fPz87OVqdDVw+bN29uVLZo8brmhQsXaufOnTUhIUH79Olz3NWr6m3NU6ZMUUCTkpKOGY43XtY8c+ZM7dy5syYlJWlaWpoOGDBAV61aFc34jeL1YztUdnb2cfeWWy/qBbZo8Pl0GLAZOAAUA38BMtTFc704G4iYiOQA76tqt1rmvQ88rKp/c8YXA/cAlwEJqvqQM30ScEhVH63rdwUCAW3Mf4RPmlPA6599y029zuLBgT+J2SRZzVZzU3Wy1exFvSKyUlUDkeSIVdN4Dliiqm844+sJNozLgMtU9d9qWy6cxjYNY4w5mXnRNGL1ltu5wK3Ou6guAvaqajGwgOAV+1TnAvgvnWnGGGOOQ558YKGIvEHwqCFNRLYRfEdUcwBV/R+C/0wyACgEDgIjnXm7RORB4AtnU1NVdZcXmYwxxnjPk6ahqsPqma/AuDDzXgRe9CKHMcaY6LL/CDfGGOOaNQ1jjDGuWdMwxhjjmjUNY4wxrlnTMMYY45o1DWOMMa5Z0zDGGOOaNQ1jjDGuWdMwxhjjmjUNY4wxrlnTMMYY45o1DWOMMa5Z0zDGGOOaNQ1jjDGuWdMwxhjjmjUNY4wxrnnSNESkn4isF5FCEZlQy/zHRSTfGTaIyJ6QeRUh8+Z6kccYY0x0RPzNfSISBzwNXAVsA74QkbmqurZqGVX9bcjydwLnh2zikKr2jDSHMcaY6PPiSOPnQKGqfqOqR4AZwHV1LD8MeMOD32uMMSbGvGgaZwDfhYxvc6b9hIhkA+2Bv4ZMThCRPBFZISIDw/0SERntLJe3Y8cOD2IbY4xpqFhfCB8KzFLVipBp2aoaAG4C/iQiHWtbUVWfV9WAqgbS09NjkdUYY0wNXjSNIuDMkPEsZ1pthlLj1JSqFjk/vwGWcOz1DmOMMccRL5rGF8DZItJeRFoQbAw/eReUiHQBUoG/h0xLFZF453YacDGwtua6xhhjjg8Rv3tKVctFZDywAIgDXlTVNSIyFchT1aoGMhSYoaoasnpX4DkRqSTYwB4OfdeVMcaY44sc+xx+YggEApqXl+d3DGOMOaGIyErnGnKj2X+EG2OMcc2aRh0OHz7MqFGjyM7OJiUlhZ49e/LBBx8AsGXLFkSE5OTk6uHBBx/0OXHk6qo51NSpUxERFi1a5ENKb9VX88GDBxk7dixpaWm0bt2a3r17+5jWG3XV/Nprrx3zuE5MTEREWLlypc+pI1Pf/Txz5ky6du1KSkoK5557LnPmzPEvrAdqqxdoVTVfRH7tfIrHfhH5UETaudqwqp5wwwUXXKCxsH//fp08ebJu3rxZKyoq9L333tPk5GTdvHmzbt68WQE9evRoTLLESl01VyksLNRu3bppZmamLly40L+wHqmv5ptvvllvvPFG3b59u5aXl2teXp6/gT3g5n6u8tJLL2mHDh20srIy9kE9VFfN27Zt0+bNm+v8+fO1srJS33//fW3ZsqWWlJT4HbvRaqsXqABygMuA7UAu0AJ4FliqLp5/fW8AjRli1TRq0717d501a1aTbRq1qaq5ytVXX63z5s3T7OzsJtE0alNV87p16zQlJUX37t3rd6Soq3k/V7nssst0ypQpPiSKvqqaV6xYoenp6cfMS0tL008//dSnZNEBHASuBx4FntaqRgDtAAU6aj3Pv3Z6qgFKSkrYsGEDubm51dOys7PJyspi5MiR7Ny508d00VGz5rfeeov4+HgGDBjgc7LoCa35888/Jzs7m8mTJ5OWlkb37t2ZPXu23xE9V9tjG2Dr1q0sW7aMW2+91adk0RNacyAQoGvXrsydO5eKigrmzJlDfHw8PXr08DumZ0pKSgASgDXOJAmZXXW7W70bqq+rHI+DH0caR44c0SuvvFJHjx6tqqqlpaX6xRdf6NGjR/WHH37Q66+/Xn/5y1/GPFc01ax537592qlTp+pTGE3xSKNmzX/4wx8U0MmTJ+vhw4d1yZIlmpSUpGvXrvU5qXdq1hxq6tSp2qdPn9iHirLaan7hhRc0KSlJ4+LitGXLlvr+++/7mNBbVfUCOzR4ZNEX2An0AFoCzwGVwDC101ORq6io0BtvvFH79++vR44cqXWZ4uJiBXTfvn0xzRYttdV811136e9///vqZZpa06it5scee0ybN29+zGnIa665Rv/0pz/5FdNT9T22O3XqpC+++KIPyaKntpoXLlyobdu21S+++EIrKir0888/14yMDP3yyy/9DeuB0HqBlfr/T0mNAzYCJcC9wF7gUq3n+ddOT9VDVRk1ahQlJSXMnj2b5s2b17qcSPDorrKyMpbxoiJczYsXL+bJJ58kIyODjIwMvvvuO4YMGcK0adN8Thy5cDXXdnqi6r4+0dX32P7kk0/4/vvvueGGG3xK6L1wNefn59O7d28CgQCnnHIKF154Ib169Trh3x1Ys16C1y2q5j2tqmer6unAbIL/7F3gaqMn2tDYI42J73ylHSbM04nvfOV6nX/7t3/TXr16aWlp6THTV6xYoV9//bVWVFTozp07dciQIXrZZZc1Klc0eVnzzp07tbi4uHrIysrSmTNn/mQ5v3lZ85EjR7Rjx446depUPXr0qP7tb3/T5ORkXbdundexI+JlzVX+9V//VYcPH+5VRM95WfOSJUv01FNPrT6y+Mc//qFt27bVBQsWeBk5Il7US/BTOiB4baMbwWsZZxH83L//UhfPv743gMYMjW0aHSbM0+x73tcOE+a5Wn7Lli0KaHx8vCYlJVUPr776qr7++uuak5OjiYmJmpGRocOHD9fi4uJG5YomL2uu6Xg9PeV1zQUFBXrRRRdpYmKidu3aVd9+++1oxm8Ur2s+dOiQtm7dWhctWhTN2BHxuuannnpKO3bsqMnJydq+fXt99NFHoxm/wbyol+Bbbm8G2gCrgQPAD8AfgTh18fwb8WdPnUhu6nUWr3/2LTf1OsvV8tnZ2VXn/mo1bNgwr6JFjdc1h9qyZUsEyaLH65pzc3P5+9//Hnb+8cDrmhMSEtizZ49H6aLD65rHjx/P+PHjvYrnOS/qFZEvVfU1Z7RRbw2zz54yxpiThH32lDHGmJiypmGMMcY1axrGGGNcs6ZhjDHGNU+ahoj0E5H1zsfsTqhl/m0iskNE8p3h1yHzRojIRmcY4UUeY4wx0RHxW25FJA54GrgK2AZ8ISJz9adf2/qmqo6vsW5bYDIQIPifiiuddXdHmssYY4z3vDjS+DlQqKrfqOoRYAZwnct1rwYWquoup1EsBPp5kMkYY0wUeNE0zgC+Cxnf5kyr6XoRWS0is0TkzAaui4iMFpE8EcnbsWOHB7GNMcY0VKwuhL8H5KhqD4JHEy83dAOq+ryqBlQ1kJ6e7nlAY4wx9fOiaRQBZ4aMZznTqqnqj6p62Bl9AbjA7brGGGOOH140jS+As0WkvYi0AIYCc0MXEJHMkNFrgXXO7QXAL0UkVURSgV8604wxxhyHIn73lKqWi8h4gk/2ccCLqrpGRKYS/BjeucC/i8i1QDmwC7jNWXeXiDxIsPEATFXVXZFmMsYYEx32gYXGGHOSsA8sNMYYE1PWNIwxxrhmTcMYY4xr1jSMMca4Zk3DGGOMa9Y0jDHGuGZNwxhjjGvWNIwxxrhmTcMYY4xr1jSMMca4Zk3DGGOMa9Y0jDHGuGZNwxhjjGvWNIwxxrhmTcMYY4xrnjQNEeknIutFpFBEJtQy/y4RWSsiq0VksYhkh8yrEJF8Z5hbc12/TZ8+nUAgQHx8PLfddtsx81544QU6depEcnIy/fr14/vvv/cnpIfqqnfmzJl07dqVlJQUzj33XObMmeNLRq8dPnyYUaNGkZ2dTUpKCj179uSDDz6onr948WK6dOlCYmIil19+OVu3bvUxrTfqqvnIkSPccMMN5OTkICIsWbLE37AeqavmFStWcNVVV9G2bVvS09MZPHgwxcXFPieOTG31Aq0ARORcEckTkd3OsEhEznWz3YibhojEAU8D/YFzgWG1/PIvgYCq9gBmAf83ZN4hVe3pDNdGmsdr7dq1Y+LEidx+++3HTF+yZAn33Xcf7777Lrt27aJ9+/YMGzbMp5TeCVdvUVERt9xyC4899hj79u3jkUce4aabbmL79u0+JfVOeXk5Z555JkuXLmXv3r089NBDDBkyhC1btrBz504GDRrEgw8+yK5duwgEAtx4441+R45YXTUDXHLJJbz66qtkZGT4G9RDddW8e/duRo8ezZYtW9i6dSspKSmMHDnS78gRqa1eoKOI5ADfAzcAbYE0gl/RPcPVhlU1ogH4BbAgZPxe4N46lj8f+CRkfH9Df+cFF1ygsXb//ffriBEjqsf/8z//U8eOHVs9XlRUpIAWFhbGPFs01Kx3xYoVmp6efswyaWlp+umnn8Y4WWx0795dZ82apc8995z+4he/qJ6+f/9+TUhI0HXr1vmYLjqqag51xhln6Mcff+xPoBiorWZV1ZUrV2pycrIPiaILOAhcr8c+JzcDxgEH1cXzrxenp84AvgsZ3+ZMC2cU8EHIeIJzmLRCRAaGW0lERjvL5e3YsSOiwF7RkK/KrbpdUFDgV5yoCgQCdO3alblz51JRUcGcOXOIj4+nR48efkfzXElJCRs2bCA3N5c1a9Zw3nnnVc9LSkqiY8eOrFmzxseE3gut+WRRV83Lli1rcvuipKQEIAGofvCKyB6gDHgK+C8322kWhWxhicgtQADoEzI5W1WLRKQD8FcR+UpVN9VcV1WfB56H4HeExyRwHfr168fQoUMZM2YMZ599NlOnTkVEOHjwoN/RoiIuLo5bb72Vm266ibKyMlq0aMFbb71FUlKS39E8dfToUW6++WZGjBhBly5d2L9/P+np6ccs07p1a0pLS31K6L2aNZ8M6qp59erVTJ06lXfffdendN6rqhf4UVW/rpquqm1EJAkYAbi6WOfFkUYRcGbIeJYz7Rgi0he4H7hWVQ9XTVfVIufnN8ASgqevjnt9+/bl97//Pddffz05OTnk5OSQkpJCVlaW39GiYtGiRfzud79jyZIlHDlyhKVLl/LrX/+a/Px8v6N5prKykuHDh9OiRQumT58OQHJyMvv27TtmuX379pGSkuJHRM/VVnNTV1fNhYWF9O/fnyeeeIJLL73Up4TeCq0X+LbmfFU9APwP8BcROa2+7XnRNL4AzhaR9iLSAhhK8KJKNRE5H3iOYMPYHjI9VUTindtpwMXAWg8yxcS4cePYuHEjJSUlXH/99ZSXl9OtWze/Y0VFfn4+vXv3JhAIcMopp3DhhRfSq1cvFi1a5Hc0T6gqo0aNoqSkhNmzZ9O8eXMAcnNzWbVqVfVyBw4cYNOmTU3i1EW4mpuyumreunUrffv2ZdKkSQwfPtzHlN6pWS8Q7izNKUAidV9aqF4w0lDlwHhgAbAOmKmqa0RkqohUvRvqESAZeKvGW2u7Ankisgr4GHhYVaPWNCbNKaDjvfOZNMf9dYfy8nLKysqoqKigoqKCsrKy6mkFBQWoKt9++y2jR4/mP/7jP0hNTY1W/EZpaM3h6r3wwgtZvnx59ZHFl19+yfLly4/LaxqNuZ/vuOMO1q1bx3vvvUfLli2rp//qV7+ioKCA2bNnU1ZWxtSpU+nRo8dxdxrHy5oh+HbNsrIyIPgW3LKysmOu4R0PvKy5qKiIK664gvHjxzNmzJhoxI2Yl/WKyFUicr6IxIlIK+AxYDfB5/C6ublafrwNjX33VIcJ8zT7nve1w4R5rteZPHmyEuzO1cPkyZN19+7d2r17d01MTNTTTz9dJ0yYoOXl5Y3KFU0NrTlcvaqqTz31lHbs2FGTk5O1ffv2+uijj0YxeeM1tOYtW7YooPHx8ZqUlFQ9vPrqq6qqunDhQu3cubMmJCRonz59dPPmzVFM3zhe15ydnf2Tx8HxVreXNU+ZMkWBY6YnJSVFuYKG8aJeoAK4GRgMfA3sB3YA84Ae6uL5V/Q4e/XgRiAQ0Ly8vAavN2lOAa9/9i039TqLBwc2zdNINVnNVnNTdbLV7EW9IrJSVQOR5DipmoYxxpzMvGga9tlTxhhjXLOmYYwxxjVrGsYYY1yzpmGMMcY1axrGGGNcs6ZhjDHGNWsaxhhjXLOmYYwxxjVrGsYYY1yzpmGMMcY1axrGGGNcs6ZhjDHGNWsaxhhjXLOmYYwxxjVPmoaI9BOR9SJSKCITapkfLyJvOvM/E5GckHn3OtPXi8jVXuQxxhgTHRE3DRGJA54G+gPnAsNE5Nwai40CdqtqJ+BxYJqz7rkEv1M8F+gHPONszxhjzHHIiyONnwOFqvqNqh4BZgDX1VjmOuBl5/Ys4EoREWf6DFU9rKqbgUJne8YYY45DXjSNM4DvQsa3OdNqXUZVy4G9wKku1wVAREaLSJ6I5O3YscOD2MYYYxrqhLkQrqrPq2pAVQPp6el+xzHGmJOSF02jCDgzZDzLmVbrMiLSDGgN/OhyXWOMMccJL5rGF8DZItJeRFoQvLA9t8Yyc4ERzu0bgL+qqjrThzrvrmoPnA187kEmY4wxUdAs0g2oarmIjAcWAHHAi6q6RkSmAnmqOhf4X+AVESkEdhFsLDjLzQTWAuXAOFWtiDSTMcaY6JDgC/4TSyAQ0Ly8PL9jGGPMCUVEVqpqIJJtnDAXwo0xxvjPmoYxxhjXrGkYY4xxzZqGMcYY16xpGGOMcc2ahjHGGNesadRj+vTpBAIB4uPjue2226qnr127lkAgQGpqKqmpqfTt25e1a9f6F9Qj4epdsWIFV111FW3btiU9PZ3BgwdTXFzsX1APhas51NSpUxERFi1aFNtwURKu5i1btiAiJCcnVw8PPvigf0E9VNf9fPDgQcaOHUtaWhqtW7emd+/e/oT0WLiaReRmEdkfMhwUERWRC+rbZsT/3NfUtWvXjokTJ7JgwQIOHTp0zPRZs2aRnZ1NZWUlTz/9NEOHDmX16tU+po1cuHp3797N6NGjufrqq2nWrBnjx49n5MiRfPjhhz6m9Ua4mqts2rSJt956i8zMTB/SRUd9Ne/Zs4dmzZrW00NdNY8ePZry8nLWrVtH27Ztyc/P9yekx8LVrKqvAa9VjYvIbcAk4B/1bbNpPSqiYNCgQQDk5eWxbdu26ult2rShTZs2AKgqcXFxFBYW+hHRU+Hq7d+//zHLjR8/nj59+sQ0W7SEq7nKuHHjmDZtGmPHjo11tKipr+amKFzNX3/9NXPnzmXbtm20atUKgAsuqPcF9wmhAffzCOAv6uK/ve30VITatGlDQkICd955J/fdd5/fcWJm2bJl5Obm+h0j6t566y3i4+MZMGCA31FiKjs7m6ysLEaOHMnOnTv9jhNVn3/+OdnZ2UyePJm0tDS6d+/O7Nmz/Y4VMyKSDfQG/uJmeWsaEdqzZw979+5l+vTpnH/++X7HiYnVq1czdepUHnnkEb+jRFVpaSn33XcfTzzxhN9RYiYtLY0vvviCrVu3snLlSkpLS7n55pv9jhVV27Zto6CggNatW/P9998zffp0RowYwbp16/yOFiu3AsudL8Krl52e8kBSUhJjxowhPT2ddevWcdppp/kdKWoKCwvp378/TzzxBJdeeqnfcaJqypQpDB8+nJycHL+jxExycjKBQPCjiU4//XSmT59OZmYmpaWlpKSk+JwuOlq2bEnz5s2ZOHEizZo1o0+fPlx++eV89NFHdO3a1e94sXAr8F9uF7YjDY9UVlZy8OBBioqa7teBbN26lb59+zJp0iSGDx/ud5yoW7x4MU8++SQZGRlkZGTw3XffMWTIEKZNm+Z3tJgJfitz8PHdVPXo0eMn06rqbupE5GKgHcGv4XblpGoak+YU0PHe+UyaU+B6nfLycsrKyqioqKCiooKysjLKy8tZuHAhX375JRUVFezbt4+77rqL1NTU4+6VSUNrDldvUVERV1xxBePHj2fMmDFRTh0Zr2pevHgxBQUF5Ofnk5+fT7t27XjuuecYN25clCtoOK9q/uyzz1i/fj2VlZX8+OOP/Pu//zuXXXYZrVu3jnIFDedVzb179+ass87ij3/8I+Xl5XzyySd8/PHHXH311VGuoGG8eP4CxPkivCojgNmqWup6o6p6wg0XXHCBNkaHCfM0+573tcOEea7XmTx5sgLHDJMnT9aZM2dq586dNSkpSdPS0nTAgAG6atWqRuWKpobWHK7eKVOmKKBJSUnHDMcjr2quKTs7WxcuXOhxWm94VfPrr7+uOTk5mpiYqBkZGTp8+HAtLi6OcvrG8fJ+Ligo0IsuukgTExO1a9eu+vbbb0cxeeN49fwFTNHgm6QSgD3AldqA59+T6vs0Js0p4PXPvuWmXmfx4MBuUUh2/LGaream6mSr2Yt6vfg+jYiahoi0Bd4EcoAtwBBV3V1jmZ7As0AroAL4g6q+6cz7M9AH2Ossfpuq5tf3e+1LmIwxpuGOhy9hmgAsVtWzgcXOeE0HgVtVNRfoB/xJRNqEzP8/qtrTGfIjzGOMMSaKIm0a1wEvO7dfBgbWXEBVN6jqRuf298B2ID3C32uMMcYHkTaN01W16lPrfgBOr2thEfk50ALYFDL5DyKyWkQeF5H4OtYdLSJ5IpK3Y8eOCGMbY4xpjHqbhogsEpGCWobrQpfT4MWRsBdIRCQTeAUYqapVb/q+F+gCXAi0Be4Jt76qPq+qAVUNpKfbgYoxxvih3v8IV9W+4eaJSImIZKpqsdMUtodZrhUwD7hfVVeEbLvqKOWwiLwE3N2g9MYYY2Iq0tNTcwn+cwjOz3drLiAiLYB3CH6C4qwa8zKdn0Lweoj7/1oxxhgTc5E2jYeBq0RkI9DXGUdEAiLygrPMEIKfoHibiOQ7Q09n3msi8hXwFZAGPBRhHmOMMVF0Uv1znzHGnMyOh//TMMYYcxKxpmGMMcY1axrGGGNcs6ZhjDHGNWsaxhhjXLOmYYwxxjVrGsYYY1yzpmGMMcY1axrGGGNcs6ZhjDHGNWsaxhhjXLOmYYwxxjVrGsYYY1yzpmGMMcY1axrGGGNci6hpiEhbEVkoIhudn6lhlqsI+QKmuSHT24vIZyJSKCJvOt/yZ4wx5jgV6ZHGBGCxqp4NLHbGa3NIVXs6w7Uh06cBj6tqJ2A3MCrCPMYYY6Io0qZxHfCyc/tlgt/z7YrzveBXAFXfG96g9Y0xxsRepE3jdFUtdm7/AJweZrkEEckTkRUiMtCZdiqwR1XLnfFtwBnhfpGIjHa2kbdjx44IYxtjjGmMZvUtICKLgIxaZt0fOqKqKiLhvnA8W1WLRKQD8FcR+QrY25Cgqvo88DwEvyO8IesaY4zxRr1HGqraV1W71TK8C5SISCaA83N7mG0UOT+/AZYA5wM/Am1EpKpxZQFFEVfksVtuuYXMzExatWrFOeecwwsvvOB3pKg6fPgwo0aNIjs7m5SUFHr27MkHH3zgd6you+yyy0hISCA5OZnk5GQ6d+7sd6SoqqqzaoiLi+POO+/0O1ZUrVu3jiuuuILWrVvTqVMn3nnnHb8jRd306dMJBALEx8dz2223VU8XkRYiMktEtoiIishlbrcZ6empucAI5/YI4N2aC4hIqojEO7fTgIuBtaqqwMfADXWt77d7772XLVu2sG/fPubOncvEiRNZuXKl37Gipry8nDPPPJOlS5eyd+9eHnroIYYMGcKWLVv8jhZ106dPZ//+/ezfv5/169f7HSeqqurcv38/P/zwAy1btmTw4MF+x4qa8vJyrrvuOq655hp27drF888/zy233MKGDRv8jhZV7dq1Y+LEidx+++21zf4bcAvBSwuuRdo0HgauEpGNQF9nHBEJiEjVS/KuQJ6IrCLYJB5W1bXOvHuAu0SkkOA1jv+NMI/ncnNziY+PB0BEEBE2bdrkc6roSUpKYsqUKeTk5HDKKadwzTXX0L59+ybdKE92s2fP5rTTTuPSSy/1O0rUfP3113z//ff89re/JS4ujiuuuIKLL76YV155xe9oUTVo0CAGDhzIqaeeesx0VT2iqn9S1b8BFQ3ZZkRNQ1V/VNUrVfVs5zTWLmd6nqr+2rn9qap2V9XznJ//G7L+N6r6c1XtpKqDVfVwJHmiZezYsSQmJtKlSxcyMzMZMGCA35FipqSkhA0bNpCbm+t3lKi79957SUtL4+KLL2bJkiV+x4mZl19+mVtvvZXgGxpPHqpKQUGB3zFOOPYf4S4888wzlJaWsnz5cgYNGlR95NHUHT16lJtvvpkRI0bQpUsXv+NE1bRp0/jmm28oKipi9OjR/Mu//EuTPqKssnXrVpYuXcqIESPqX/gE1rlzZ0477TQeeeQRjh49ykcffcTSpUs5ePCg39FOONY0XIqLi+OSSy5h27ZtPPvss37HibrKykqGDx9OixYtmD59ut9xoq5Xr16kpKQQHx/PiBEjuPjii5k/f77fsaLulVde4ZJLLqF9+/Z+R4mq5s2bM2fOHObNm0dGRgb//d//zZAhQ8jKyvI72gmn3rfcmmOVl5c3+VegqsqoUaMoKSlh/vz5NG/e3O9IMSciBN+r0bT95S9/YcKEcB/k0LT06NGDpUuXVo//0z/9U5M/woqGk+pIY9KcAjreO59Jc9ydx9y+fTszZsxg//79VFRUsGDBAt544w2uvPLKKCf1TkNrBrjjjjtYt24d7733Hi1btoxiuuhoaM179uxhwYIFlJWVUV5ezmuvvcayZcvo169flJN6pzH386effkpRUdEJ+66phta8evVqysrKOHjwII8++ijFxcXHvA31eNeY+7i8vJyysjIqKiqoqKiA4IdxNHNuxItIgrNoCxFJEBcXtk6qpvH6Z99Socrrn33rankR4dlnnyUrK4vU1FTuvvtu/vSnP3HttdfWv/JxoqE1b926leeee478/HwyMjKq38f/2muvRTmpdxpa89GjR5k4cSLp6emkpaXx1FNPMWfOHM4555woJ/VOQ2uG4AXwQYMGkZKSEsVk0dPQml955RUyMzM57bTTWLx4MQsXLjyhrk825j5+6KGHaNmyJQ8//DCvvvoqwM+Aic7s9cAhgp/EscC5nV3fNuVEPAQPBAKal5fX4PUmzSng9c++5aZeZ/HgwG5RSHb8sZqt5qbqZKvZi3pFZKWqBiLJcVI1DWOMOZl50TROqtNTxhhjImNNwxhjjGvWNIwxxrhmTcMYY4xr1jSMMca4Zk3DGGOMa9Y0jDHGuGZNwxhjjGvWNIwxxrgWUdMQkbYislBENjo/U2tZ5nIRyQ8ZykRkoDPvzyKyOWRez0jyGGOMia5IjzQmAItV9WxgsTN+DFX9WFV7qmpP4ArgIPBRyCL/p2q+quZHmMcYY0wURdo0rgNedm6/DAysZ/kbgA9U1b4uyxhjTkCRNo3TVbXYuf0DcHo9yw8F3qgx7Q8islpEHheRsJ9TLCKjRSRPRPJ27NgRQWRjjDGNVW/TEJFFIlJQy3Bd6HIa/LjcsB+ZKyKZQHeCn9te5V6gC3Ah0Ba4J9z6qvq8qgZUNZCenl5fbGOMMVFQ79e9qmrfcPNEpEREMlW12GkK2+vY1BDgHVU9GrLtqqOUwyLyEnC3y9zGGGN8EOnpqblA1ZfsjgDerWPZYdQ4NeU0GpyvGBwIuP8eQ2OMMTEXadN4GLhKRDYCfZ1xRCQgIi9ULSQiOcCZwNIa678mIl8BXwFpwEMR5jHGGBNF9Z6eqouq/ghcWcv0PODXIeNbCH4Pbc3lrojk9xtjjIkt+49wY4wxrlnTMMYY45o1DWOMMa5Z0zDGGOOaNQ1jjDGuWdMwxhjjmjUNY4wxrlnTMMYY45o1DWOMMa5Z0zDGGOOaNQ1jjDGuWdOox65du/jVr35FUlIS2dnZvP76635Hiqrp06cTCASIj4/ntttu8ztOTG3cuJGEhARuueUWv6PEzIwZM+jatStJSUl07NiR5cuX+x0pqrZs2cKAAQNITU0lIyOD8ePHU15e7nesqKnr71lErhSRr0XkoIh8LCLZbrZpTaMe48aNo0WLFpSUlPDaa69xxx13sGbNGr9jRU27du2YOHEit99+u99RYm7cuHFceOGFfseImYULF3LPPffw0ksvUVpayrJly+jQoYPfsaJq7NixnHbaaRQXF5Ofn8/SpUt55pln/I4VNeH+nkUkDXgbmETwC/DygDfdbNOaRh0OHDjA7NmzefDBB0lOTuaSSy7h2muv5ZVXXvE7WtQMGjSIgQMHcuqpp/odJaZmzJhBmzZtuPLKn3xoc5M1efJkHnjgAS666CJOOeUUzjjjDM444ycfRt2kbN68mSFDhpCQkEBGRgb9+vVr0i8C6/h7HgSsUdW3VLUMmAKcJyJd6tumNY06bNiwgWbNmnHOOedUTzvvvPOa9IPsZLRv3z4eeOABHnvsMb+jxExFRQV5eXns2LGDTp06kZWVxfjx4zl06JDf0aLqN7/5DTNmzODgwYMUFRXxwQcf0K9fP79j+SEXWFU1oqoHgE3O9DpZ06jD/v37adWq1THTWrduTWlpqU+JTDRMmjSJUaNGkZWV5XeUmCkpKeHo0aPMmjWL5cuXk5+fz5dffslDDzXt70Hr3bs3a9asoVWrVmRlZREIBBg4cKDfsfyQDOytMW0vkFLfihE1DREZLCJrRKRSRAJ1LNdPRNaLSKGITAiZ3l5EPnOmvykiLSLJ47Xk5GT27dt3zLR9+/aRklLvfjUniPz8fBYtWsRvf/tbv6PEVMuWLQG48847yczMJC0tjbvuuov58+f7nCx6Kisr6devH4MGDeLAgQPs3LmT3bt3c8899/gdzQ/7gVY1prUC6n1FHOmRRgHBc2PLwi0gInHA00B/4FxgmIic68yeBjyuqp2A3cCoCPN46pxzzqG8vJyNGzdWT1u1ahW5ufUewZkTxJIlS9iyZQtnnXUWGRkZPProo8yePZuf/exnfkeLqtTUVLKyshCR6mmht5uiXbt28e233zJ+/Hji4+M59dRTGTlyZJNulHVYA5xXNSIiSUBHZ3qdImoaqrpOVdfXs9jPgUJV/UZVjwAzgOsk+Ai9ApjlLPcyMDCSPPWZNKeAjvfOZ9KcAlfLJyUlMWjQIB544AEOHDjAJ598wrvvvsvw4cOjGdNTDa25vLycsrIyKioqqKiooKys7IR7S2JDah49ejSbNm0iPz+f/Px8xowZwz//8z+zYMGCGCT1TkPvZ4CRI0fy1FNPsX37dnbv3s3jjz/ONddcE8WU3mpozWlpabRv355nn32W8vJy9uzZw8svv0yPHj2inNQbjbmPa/49AyIizYB3gG4icr2IJAAPAKtV9ev6thmLaxpnAN+FjG9zpp0K7FHV8hrTayUio0UkT0TyduzY0aggr3/2LRWqvP7Zt67XeeaZZzh06BCnnXYaw4YN49lnnz2hjjQaWvNDDz1Ey5Ytefjhh3n11Vdp2bLlCXeeuyE1JyYmkpGRUT0kJyeTkJBAenp6DJJ6pzGP7UmTJnHhhRdyzjnn0LVrV84//3zuv//+KKb0VmNqfvvtt/nwww9JT0+nU6dONG/enMcffzyKKb3TmHpr/j0DPwMmquoO4HrgDwTP8vQChrrZpqhq3QuILAIyapl1v6q+6yyzBLhbVfNqWf8GoJ+q/toZH+4EnAKscE5NISJnAh+oarf6QgcCAc3L+8mvqtekOQW8/tm33NTrLB4cWO+vaRKsZqu5qTrZavaiXhFZqaphrz+72kZ9TcNlkCWEbxq/AKao6tXO+L3OrIeBHUCGqpbXXK4ujW0axhhzMvOiacTi9NQXwNnOO6VaEDwEmqvBbvUxcIOz3Ajg3RjkMcYY00iRvuX2VyKyDfgFME9EFjjT24nIfADnmsV4YAGwDpipqlVX6O8B7hKRQoLXOP43kjzGGGOiy5PTU7Fmp6eMMabhTpTTU8YYY5oIaxrGGGNcs6ZhjDHGNWsaxhhjXDshL4SLyA5gayNXTwN2ehjHK5arYSxXw1iuhmmqubJVNaKPOzghm0YkRCQv0ncPRIPlahjL1TCWq2EsV3h2esoYY4xr1jSMMca4djI2jef9DhCG5WoYy9UwlqthLFcYJ901DWOMMY13Mh5pGGOMaSRrGsYYY1xrkk1DRAaLyBoRqRSRsG9PE5F+IrJeRApFZELI9PYi8pkz/U3nI929yNVWRBaKyEbnZ2oty1wuIvkhQ5mIDHTm/VlENofM6xmrXM5yFSG/e27IdD/3V08R+btzf68WkRtD5nm6v8I9XkLmxzv1Fzr7Iydk3r3O9PUiUu93xnic6y4RWevsn8Uikh0yr9b7NEa5bhORHSG//9ch80Y49/tGERkR41yPh2TaICJ7QuZFZX+JyIsisl1Eav0uVwl60sm8WkR+FjIvavuqVqra5AagK9AZWAIEwiwTB2wCOgAtgFXAuc68mcBQ5/b/AHd4lOv/AhOc2xOAafUs3xbYBSQ6438GbojC/nKVC9gfZrpv+ws4Bzjbud0OKAbaeL2/6nq8hCwzFvgf5/ZQ4E3n9rnO8vFAe2c7cTHMdXnIY+iOqlx13acxynUbML2WddsC3zg/U53bqbHKVWP5O4EXY7C/ehP8KtaCMPMHAB8AAlwEfBbtfRVuaJJHGqq6TlXX17PYz4FCVf1GVY8AM4DrRESAK4BZznIvAwM9inadsz23272B4FfgHvTo94fT0FzV/N5fqrpBVTc6t78HtgPR+ILvWh8vdeSdBVzp7J/rgBmqelhVNwOFzvZikktVPw55DK0Asjz63RHlqsPVwEJV3aWqu4GFQD+fcg0D3vDod4elqssIvkAM5zrgLxq0AmgjIplEd1/Vqkk2DZfOAL4LGd/mTDsV2KPBL48Kne6F01W12Ln9A3B6PcsP5acP2D84h6ePi0h8jHMliEieiKyoOmXGcbS/ROTnBF89bgqZ7NX+Cvd4qXUZZ3/sJbh/3KwbzVyhRhF8xVqltvs0lrmud+6fWSJyZgPXjWYunNN47YG/hkyO1v6qT7jc0dxXtWoWzY1Hk4gsAjJqmXW/qvr2tbF15QodUVUVkbDvd3ZeRXQn+I2HVe4l+OTZguD7te8BpsYwV7aqFolIB+CvIvIVwSfGRvN4f70CjFDVSmdyo/dXUyQitwABoE/I5J/cp6q6qfYteO494A1VPSwi/0bwKO2KGP1uN4YCs1S1ImSan/vruHDCNg1V7RvhJoqAM0PGs5xpPxI89GvmvFqsmh5xLhEpEZFMVS12nuS217GpIcA7qno0ZNtVr7oPi8hLwN2xzKWqRc7Pb0RkCXA+MBuf95eItALmEXzBsCJk243eX7UI93ipbZltItIMaE3w8eRm3WjmQkT6EmzEfVT1cNX0MPepF0+C9eZS1R9DRl8geA2rat3Laqy7xINMrnKFGAqMC50Qxf1Vn3C5o7mvanUyn576Ajhbgu/8aUHwATJXg1eXPiZ4PQFgBODVkctcZ3tutvuTc6nOE2fVdYSBQK3vtIhGLhFJrTq9IyJpwMXAWr/3l3PfvUPwfO+sGvO83F+1Pl7qyHsD8Fdn/8wFhkrw3VXtgbOBzyPI0qBcInI+8BxwrapuD5le630aw1yZIaPXAuuc2wuAXzr5UoFfcuwRd1RzOdm6ELyw/PeQadHcX/WZC9zqvIvqImCv86IomvuqdtG8yu7XAPyK4Lm9w0AJsMCZ3g6YH7LcAGADwVcK94dM70Dwj7oQeAuI9yjXqcBiYCOwCGjrTA8AL4Qsl0PwFcQpNdb/K/AVwSe/V4HkWOUC/sn53aucn6OOh/0F3AIcBfJDhp7R2F+1PV4Inu661rmd4NRf6OyPDiHr3u+stx7o7/Hjvb5ci5y/g6r9M7e++zRGuf4IrHF+/8dAl5B1b3f2YyEwMpa5nPEpwMM11ova/iL4ArHYeSxvI3jtaQwwxpkvwNNO5q8IeVdoNPdVbYN9jIgxxhjXTubTU8YYYxrImoYxxhjXrGkYY4xxzZqGMcYY16xpGGOMcc2ahjHGGNesaRhjjHHt/wH/j7RcGiDGTQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "textdx=0.2*(2/numpy.sqrt(Npts))\n",
    "plt.scatter(points[:,0],points[:,1],s=4)\n",
    "for i,p in enumerate(points):\n",
    "    plt.text(p[0]+textdx,p[1],str(i),fontsize=12)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "35d2f81e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bd1st (len=20):\n",
      " {0: -3.0, 1: -2.600000023841858, 3: -2.200000047683716, 4: -2.199999988079071, 6: -1.7999999523162842, 8: -1.399999976158142, 10: -1.0, 11: -0.20000004768371582, 13: -1.399999976158142, 17: 0.6000000238418579, 19: -0.5999999046325684, 23: 1.4000000953674316, 25: 0.20000004768371582, 29: 2.200000047683716, 30: 1.399999976158142, 31: 1.0, 32: 1.800000011920929, 33: 2.200000047683716, 34: 2.600000023841858, 35: 3.0}\n",
      "\n",
      "bd2nd (len=0):\n",
      " {}\n"
     ]
    }
   ],
   "source": [
    "# define boundary conditions\n",
    "def bdfun(p):\n",
    "    return p[0]*1+p[1]*2\n",
    "\n",
    "counter=numpy.zeros(Npts,dtype=numpy.uint)\n",
    "for c in cells:\n",
    "    for i in c:\n",
    "        counter[i]+=1\n",
    "        \n",
    "bd1st={};bd2nd={}\n",
    "for i in range(Npts):\n",
    "    if counter[i]==4:\n",
    "        continue\n",
    "    if numpy.isclose(points[i][1],-1) or numpy.isclose(points[i][1],1) or\\\n",
    "       numpy.isclose(points[i][0],-1) or numpy.isclose(points[i][0],1):\n",
    "        bd1st[i]=bdfun(points[i])\n",
    "    else:\n",
    "        bd2nd[i]=None\n",
    "print(\"bd1st (len=%d):\\n\"%(len(bd1st)),bd1st)\n",
    "print(\"\")\n",
    "print(\"bd2nd (len=%d):\\n\"%(len(bd2nd)),bd2nd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "11ba0be5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{2: 0, 5: 1, 7: 2, 9: 3, 12: 4, 14: 5, 15: 6, 16: 7, 18: 8, 20: 9, 21: 10, 22: 11, 24: 12, 26: 13, 27: 14, 28: 15}\n"
     ]
    }
   ],
   "source": [
    "counter=0\n",
    "g2i={}\n",
    "for i in range(Npts):\n",
    "    if i not in bd1st:\n",
    "        g2i[i]=counter\n",
    "        counter+=1\n",
    "Nipts=len(g2i)\n",
    "print(g2i)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f610c1b",
   "metadata": {},
   "source": [
    "## Then begins the real FEM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "31316007",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm import tqdm\n",
    "import numpy,scipy.integrate,scipy.sparse.linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "e591aadd",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 25/25 [00:01<00:00, 15.73it/s]\n"
     ]
    }
   ],
   "source": [
    "def gen_K_b(kappa=1.0):\n",
    "    Kval=[];Kij=[]\n",
    "    b=numpy.zeros(Nipts)\n",
    "    for c in tqdm(cells):\n",
    "        gen_K_cell(c,Kval,Kij,b,kappa)\n",
    "    row=[i[0] for i in Kij]\n",
    "    col=[i[1] for i in Kij]\n",
    "    K=scipy.sparse.coo_matrix((Kval,(row,col)),shape=(Nipts,Nipts))\n",
    "    return K.tocsr(),b\n",
    "\n",
    "def gen_K_cell(cell,Kval,Kij,b,kappa):\n",
    "    xyval=[(xs[i],points[p][0]) for i,p in enumerate(cell)]+[(ys[i],points[p][1]) for i,p in enumerate(cell)]\n",
    "    for i in range(4):\n",
    "        if cell[i] not in bd1st:\n",
    "            f=sympy.lambdify([xi,eta],Ksym[i][i].subs(xyval))\n",
    "            k=kappa*scipy.integrate.dblquad(f,-1,1,-1,1)[0]\n",
    "            indi=g2i[cell[i]]\n",
    "            Kval.append(k)\n",
    "            Kij.append((indi,indi))\n",
    "\n",
    "        for j in range(i):\n",
    "            if cell[i] in bd1st and cell[j] in bd1st:\n",
    "                continue\n",
    "            \n",
    "            f=sympy.lambdify([xi,eta],Ksym[i][j].subs(xyval))\n",
    "            k=kappa*scipy.integrate.dblquad(f,-1,1,-1,1)[0]\n",
    "            \n",
    "            \n",
    "            if cell[i] in bd1st:\n",
    "                #print(\"detect boundary: %d->%d\"%(cell[j],cell[i]))\n",
    "                indj=g2i[cell[j]]\n",
    "                b[indj]-=bd1st[cell[i]]*k\n",
    "                del indj\n",
    "            elif cell[j] in bd1st:\n",
    "                #print(\"detect boundary: %d->%d\"%(cell[i],cell[j]))\n",
    "                b[indi]-=bd1st[cell[j]]*k\n",
    "            else:\n",
    "                indj=g2i[cell[j]]\n",
    "                Kval.append(k);Kval.append(k)\n",
    "                Kij.append((indi,indj))\n",
    "                Kij.append((indj,indi))\n",
    "        try:\n",
    "            del indi\n",
    "        except:\n",
    "            pass\n",
    "        try:\n",
    "            del indj\n",
    "        except:\n",
    "            pass\n",
    "\n",
    "K,b=gen_K_b()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "4058206e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K:\n",
      "   (0, 0)\t2.6666666666666674\n",
      "  (0, 1)\t-0.3333333333333333\n",
      "  (0, 4)\t-0.3333333333333333\n",
      "  (0, 5)\t-0.33333333333333337\n",
      "  (1, 0)\t-0.3333333333333333\n",
      "  (1, 1)\t2.6666666666666674\n",
      "  (1, 2)\t-0.3333333333333333\n",
      "  (1, 4)\t-0.33333333333333337\n",
      "  (1, 5)\t-0.3333333333333333\n",
      "  (1, 6)\t-0.33333333333333337\n",
      "  (2, 1)\t-0.3333333333333333\n",
      "  (2, 2)\t2.6666666666666674\n",
      "  (2, 3)\t-0.3333333333333333\n",
      "  (2, 5)\t-0.33333333333333337\n",
      "  (2, 6)\t-0.3333333333333333\n",
      "  (2, 7)\t-0.33333333333333337\n",
      "  (3, 2)\t-0.3333333333333333\n",
      "  (3, 3)\t2.6666666666666674\n",
      "  (3, 6)\t-0.33333333333333337\n",
      "  (3, 7)\t-0.3333333333333333\n",
      "  (4, 0)\t-0.3333333333333333\n",
      "  (4, 1)\t-0.33333333333333337\n",
      "  (4, 4)\t2.6666666666666674\n",
      "  (4, 5)\t-0.3333333333333333\n",
      "  (4, 8)\t-0.3333333333333333\n",
      "  :\t:\n",
      "  (11, 7)\t-0.3333333333333333\n",
      "  (11, 10)\t-0.3333333333333333\n",
      "  (11, 11)\t2.6666666666666674\n",
      "  (11, 14)\t-0.33333333333333337\n",
      "  (11, 15)\t-0.3333333333333333\n",
      "  (12, 8)\t-0.3333333333333333\n",
      "  (12, 9)\t-0.33333333333333337\n",
      "  (12, 12)\t2.6666666666666674\n",
      "  (12, 13)\t-0.3333333333333333\n",
      "  (13, 8)\t-0.33333333333333337\n",
      "  (13, 9)\t-0.3333333333333333\n",
      "  (13, 10)\t-0.33333333333333337\n",
      "  (13, 12)\t-0.3333333333333333\n",
      "  (13, 13)\t2.6666666666666674\n",
      "  (13, 14)\t-0.3333333333333333\n",
      "  (14, 9)\t-0.33333333333333337\n",
      "  (14, 10)\t-0.3333333333333333\n",
      "  (14, 11)\t-0.33333333333333337\n",
      "  (14, 13)\t-0.3333333333333333\n",
      "  (14, 14)\t2.6666666666666674\n",
      "  (14, 15)\t-0.3333333333333333\n",
      "  (15, 10)\t-0.33333333333333337\n",
      "  (15, 11)\t-0.3333333333333333\n",
      "  (15, 14)\t-0.3333333333333333\n",
      "  (15, 15)\t2.6666666666666674\n",
      "b:\n",
      " [-3.80000001 -2.19999999 -1.79999997 -1.26666665 -1.39999998  0.\n",
      "  0.          0.60000002 -0.59999994  0.          0.          1.40000006\n",
      "  1.26666671  1.80000001  2.20000003  3.80000007]\n"
     ]
    }
   ],
   "source": [
    "print(\"K:\\n\",K)\n",
    "print(\"b:\\n\",b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "18d75f5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "x=scipy.sparse.linalg.spsolve(K,b)\n",
    "#print(\"x:\",x)\n",
    "gx=numpy.zeros(Npts)\n",
    "for i in range(Npts):\n",
    "    if i in bd1st:\n",
    "        gx[i]=bd1st[i]\n",
    "    else:\n",
    "        gx[i]=x[g2i[i]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "95c2fd21",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOsAAADsCAYAAACVD1FtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABsQElEQVR4nO2dd3gU5drGfzM7u5uebHolNEGkiiAoxV6OR1E/sRfAcjwqdkUEK0WxKxZAURBUUI8FFcXu0SMooXdCJ6QnpCfb5/tjmM2W2c0CqbD3deXKJjOzM7M79/u871PuR5BlmRBCCKH9Q2zrCwghhBCCQ4isIYTQQRAiawghdBCEyBpCCB0EIbKGEEIHQYisIYTQQSA1sT0U1wkhhJaHEMxOIcsaQggdBCGyhhBCB0GIrCGE0EEQImsIIXQQhMgaQggdBCGyhhBCB0GIrCGE0EEQImsIIXQQhMgaQggdBCGyhhBCB0GIrCGE0EEQImsIIXQQhMgaQggdBCGyhhBCB0GIrCGE0EEQImsIIXQQhMgaQggdBCGyhhBCB0GIrCGE0EEQImsIIXQQhMgaQggdBE2pG4YQQptCbZzm/luWZXQ6HYIQlCjgMYMQWUNoNWgRD8DpdLr+tlgsWCwWoqKi8Nfh0GKxEBERgcFgaIWrbj8IkTWEI4I/4qmWz/219z4AgiBokrG6upqysjJ69OiBw+HAbrdjtVqxWq3YbDasVisOh4OMjAz0ev1xZV1DZA0hoMVzOp2IohiQeP7gTiSn0+kinEo6dwKqrxsaGrDb7Rw8eBCdTofBYECv17t+6/V6SktLCQsLIyIiAr1e39wfR7tFiKzHGNxJFGi6GQzxqqurycvLo3fv3j7EU0mnRTz33+p7i6LoIp07AcPDwz3IWFFRQU1NDd27d/d7f5GRkWzfvp2IiAhSUlKOG+saIms7RnMSzx3uD7csy9hsNtd00518tbW1VFZWsnbtWqxWq+u9BUHwIJj6OzIyEpPJ5PE/UTy8gIMsy02SLzIyEoPBwIYNGzj77LORpOPjMT4+7rKdIBgHiz/ilZeXY7PZSE1N9XhPb+I5HA4fC6dl9RwOh+s4dXrpbvXCw8PR6XQ4HA569eqFXq9Hp9O1wKfiCXXa3RSysrJYv349DQ0NREVFHRfWNUTWI0QwDhZvEgYDremm1Wrl4MGDmM1mTeLZ7XbXMZIk+RDPaDQSFRXl8b9gQh+VlZXU1dURFhYW5Kdy9AjGssqyTEREBIIgUF5eTnh4+HFhXY/9OwwCh+vZLCwsJDk52WUB/Hk2va2eSix/zhX1twqdTucinsViARQrqIYtVOJJktQiliXYAaY5EaxldTqdREREsH//fhITE4mMjDzmresxR9ajDSn4g/uDcODAAaKjoxEEoUniua/z3Nd46uvo6GiftZ/WQ1dYWEhDQwMZGRlH8ekcPlqbALIsB01Wg8GALMuUl5cTFhZ2zFvXdn13TTlY7HY7DocDvV7fIg4W79fqtLa+vp6tW7diNBo9iBcZGUlcXJzH/w7XwdKe0FaWNZi1sdPpRBAEOnfuzK5du0hISDjmrWu7JqvVanURRAtFRUWYzWa6dOkCHJ6DxW63u75YdweLe0ghJibGZ50HsGbNGnr16kV4eHjLfwiH4G+q3ZIIZv3YEucMZoBT91NnOBUVFce8dW3Xd6amn/kjXk1NDRaLhaKiItcxaiDdnXjuDhb1f0eTW9oWxDleoFrMYPZTSd2lSxf27NmDyWQ6pq1ruyZrXl4eZWVlPus8NXOlrq6O6upqevXq1apfkCiKAS1+S6AtHsC2sKyH42BS94uJiUGWZSorK49p69qu76pLly506tTJ73ZZlqmpqWn1B6qtLOvxMg0O1rK679elSxf27dtHXFwcUVFRLXmJbYaO6/0AV9C+taHmyobQ/AjWsnqvbWNjY3E4HFRVVXnEnY8ltGuyNjXCtsV0FJTraotp8PFiWQOR1el0IsuyJqlV62o2m1v6MtsE7Xoa3BTakqwhy9o8cDgcLseh6jQEKC0t1UyPVLOXEhMTfcgaFxfHnj17jlnPcIe+m7acBoccTNr7u3vtvb347n+rEEXR5Tw0GAyuuHlsbKym997pdLJy5UoaGho047GdO3dm27ZtOBwOOnfu3BwfQ7tBhybr8WZZW/ucauJJVVVVQOJ5FwW4k8+9Gse9JtXfVHfDhg2kpKQQHR3t97o6derE3r17fYoaQLGuDocDs9mMw+FoleKD1kK7JmtTo3rIwRQ8vK1eIAuoQv1s7Xa7BwnVErXmiFl7IxgHk8lkIjc3V/O7FwSB6OhoSkpKyMrKIjIyslmuqz2gXZO1KRzPDibvtd7hWj31daAa1KKiIurr6+natWur3WcwU29BEIiLi+PgwYNkZ2f7bJckCbPZTFVVFWFhYceMde3QZG2rTJXmnga7Wz1v0qn/q62txWw2U15eDniu9VrK6rXnqpvw8HAqKyuxWq0+wmmyLJORkUFeXh6xsbHHjHXt0GRtKzRl0d2tXqDpppbV8yaeWhhQV1dHRUUFvXr1ao1bBNpn6Ma9YCMxMZH9+/f7SMDIsozJZKKwsJCamppjxrqGyNoEtKxeVVUVsixTVVXlV2/In9Vz//twrJ7Vam2TmURbpBsGc05ZlomPj2f37t3YbDYP4TS1cqdz587s37+fmJgYIiIiWvKyWwXtmqwt8aCoVq+pEEMg2RMAg8FAcnLyMVMOp4W28ngHmxus0+nIyMjgwIEDrsorUL5jURRJSEhg7969VFdXYzQaO7x1bddkbQpqAXltbW1A8mlZPe8Qg7v6QlNWLy8vD1mWSUhIaLV7PV4ymA636iYtLY2cnByysrJcSRAq4QVBIDs7m7y8vGPCurZrsjqdTvbu3RvQ6tXX17Nz504PonlPN5vb6rWFN/h4QVMDhLv4gCAIiKJIeno6BQUFrqIPd8InJiayd+9eampqOrx1bddkFQQBSZICWr2//vqLvn37tuqXIIqih1ZSa+B4sawQ3PLHfbqcnp7OqlWryMjIQKfTeWxzt67R0dEd2rq260WWKIpkZmaSnJxMXFwcERERPuJgbZX619GSIo41uId4dDodKSkpFBYWau6blJREdXU1tbW1bZJE01xo12QNBm2RxdQRM5iOBG1lWZuCugxxX9pkZGRQUFCgOXALgkBWVhZ5eXkulciOiHY9DYamrVhbWdbjxZq3JlnV3jrV1dU+iSHuzsIePXr4kFWv15OQkEBxcbHme6ekpLB//36Kioro1KlTh6zI6XhX7IW2sKzHyzT4aO9RSzFS68ddrtVsNrNv3z4fh6HqKKyqqmLfvn2aVl9V6fen4ZyVlUVubi4JCQnExsYe1b21BTo8WdvCsh4v62R3Qsiy7CFQrmXxvAsBvLvAqeJ1qlay+uNuIZcvX07fvn39XpMkSezZs0czHmswGIiLi6O0tFTz2JSUFLZv347ZbCY6OrrDxcVDZD0CHAuWVav21PunuroaWZbZt28foBDFOzMrIiKiVbWS3S2k1nkyMzPJz8/XtLxqjD0/P79DeobbPVmbIkbIwaTAvS9OoJ+mak/de+IUFRURHh7e6l0A5swRefttHX//bUNraZmUlMTWrVux2Ww+SfxqLLW8vJzExESfY3U6HRUVFdTX1xMWFtahrGu7J2tTOFYdTN4lcNXV1dTU1JCbm+vxf/U63Nswuv/ExsZ6/P9w+uLodLpWfZh37oTHHuvD8uVKnu9NN0ksXGjHO4Su3uuBAwd8yvecTidhYWHs27ePhIQETeuanp7OgQMHiIiI6FDWtcOTtSNYVrVTQKB1npazxZ10qvdTVVxwTxBpKbTW7KGqCmbM0PHGGzokycTUqXYaGuCZZyTCw+HttxsJq16T2gHd27PrdDpdudyVlZWYTCaf86WlpbFq1aoOZ107PFnbqgDdbrdTV1fXos4Wd9TU1LBnzx6SkpJa6xZbPM7qcMD8+SJPPy1RWgo33GDn8stXc9FFJwMgSTBlioQgwJw5nhZWEARSU1MpLCwkKyvL9X91UMvOzmbnzp2aZFWta35+PpGRka3aBuVo0O7JGowc6dFa1mCcLTabzZViqPbT2bFjR5s6Wzoyfv9d4OGHJdavFzntNCdffmmnb18r69Y1av5OmuRAlmHqVOUxnTPHUw84PT2dNWvWkJGR4fqcVbJGRUWh0+morq4mJibG5/yqdc3IyMBoNHaI76ndk7Up6HQ6nzxdp9MZVHxPFYPWWu95O1vc2zHW19ezfft2BgwY0AZ33HpoCcu6Zw9MmiTxxRc6srJkFi60MXq0E0EAi8W34mbyZGUgnjpV4vPPRfbutWI0KtskSSIxMZHi4mLS0tIAzzTE7Oxs9u7dS79+/Vz3o0Kn07kKADqKde0QZA2kN1RZWelqTtWSzhZ3HE/houYia00NPP+8jpkzdeh08MQTdu6/34E7R7QkXfLzYdkyEZCprxe4+OJwli2rd23PzMxk/fr1pKam+qQhxsTE4HQ6qa2tJSoqyuf93QsAOoJ1bfdkzcvL88lo8S6Fa2hooGfPnq1WeXMsxFmDQXPco9MJH3wg8sQTEkVFAtdd52DqVDta0SB3S261wp136vjoIx1Op0DXrk5OPdXO4sV67r47jLFjlWMMBgMxMTGUlZWRlJTkQ8js7Gz27dtH7969fbbpdDrS0tJc1jUsLOyo77cl0e7J2qlTJ9LT0/1uV7uPt3aJXMiyNo3ly5V16erVIqee6uSTT2yceqr/e1DJ9OabIo89JtHQIBAVJfPKKzZuvFFpm9G5s5MZM4yUl/di8WIQReUZ2bJlC4mJiT6ENJlM7Nmzh/r6ekRR9HlOMjIyWLVqFenp6QEdfO0B7ffKgkQoN7jlcKT3uH8/3HijxNlnGygqEnjvPRu//RaYqAB//qnjkksG8OCD+kOW1U5JiZUbb2wcGCdNsjJhQgPLlqUxfnwYTqeidBgWFkZlZaXmVFq1rlrbdDqdy6vs7r1vj2j3lrUpHC+5wW2Bw3Uw1dXBSy/peOUVHbIMjz5q56GHHDSlBFpYCNdco+fvv5VspOHDbcyZc5DoaCsFBYpvQm2qLUkSEyZkU1JSwvz52axfL/D11w1kZ2eza9cu0tLSfAgZHx/Pnj17aGho0LScGRkZrF69ut1b1w5PVp1Od1xMSdvzNFiWYfFiZeqany9w5ZUOpk2z06mTkgxSX6/tja+vt/LMM1ksXZqCLAukptYzceI6hgwRMZsNOJ2KbyI8PNzlHMzNzaWqqoqRI6v57LMsNm6U6N8/kkPpy67prvc9dOrUicLCQs3SOEmSSElJYcOGDfTp00czTbE9oN2TtTXirIeL9liQ3RLQGhwcDgcWiwWbzYbFYiEnR2DKlATWrzdy4ol1vPnmLk466SAHDjg5cKAx+V+v12M0GtHr9YSHh/PFF2lMmxZDQ4NIZKTMCy/YuOyyKoqLI/xqI8uyTExMNlddFcnq1YqaYWysk6oqkbvvNjJ1aja7d+/U7IGTlJTErl27iIuL03zvzMxM9u3bh8ViCVpovLXR7snaFI6XKWlLDRBaaZDqT2lpKRUVFQCuz1jNxPrppxTeeiuToiIDSUl2Xn65guuusxMWlo3BcILfh335chgzRk9enogoyvzrX3ZeftmBJEFZmX+SOJ0wYYLErFkZOBwCCQkW5s+3M3Kkk+nTDTz/vBFBSOHaa7dqNlMWBIHExESqqqo0318dVMrLy0lISGiXnuEOT9a2ak7VFghmGhyIfFoFACr53H/UNEiz2UxKSgrJyckuEjU0wKuv6nj6acWrOmKEg88/dxAdHTghvqgIrrtOz/LlyqAzfLiTjz6ykZzseX9aZP3oI5F775WoqREIC5N54IEizjtvI0OHDgVg8mQrsgwvvGCksnIA99+/TbPdY2xsLMXFxT6i4CoMBoOroZXRaGx3M6gOT9Zj3bKq2Vi1tbVYLBYKCgpcDhd1KqqSTxAEnwIAlXxRUVEYjcYmc5DdUVJS4koWkGX47DORSZMk9u8XOOMMB8uXi/zxh8gffwhcdJH2QGK3wz336Jg/X4mXdurkZMECG4d45nOv7gRZv14h+K5dAoIA111n57XX6mlosLFxo5IYYzAYEAR47DHFk/vCC8lYrVYWLLASFubbAycmJsZHFNx9u5oRpXqY2xPaPVmDWbN2JLKq5Au27lQtmNbpdC6SukudtEb1zZo1Ag89JLF8uUi/fk7eecfGGWfILF8O555r4Mor9fz9t5U+fTyPe/ddkQkTJOrqBCIiZJ57zsZtt/n/rlTLWlkJ118v8fPPyoAyaJDMokU2MjNlrFaor5eJioryKJFTCVteXsl772WSkeGkuLjOox7W6XQSGxtLUVGRhyi4+/asrCzWrVtHSkpKu7Ou7Z6sTaGtP0z3PGTVygVDPq1srEDkM5vNbN68WdMitBRKS3U8/3wsH3+sJzER3nzTxtixTlf1y+mnK+Vrt94qccYZBrZvt5KYCCtXwg036Nm/X1mX3nyznZkzHZqF5O6w250891wyH3xgwOEQSE2Vee89K2efrWxXVwFOp5Po6GhKS0vJzs52fVbz5kl89JHiybXZRG6+2cD8+VbUSYTD4UCSJFfWkioK7g6DweASXmtv1rXDk7W5EaziQm1tLX/++SeiKGrmIXsLkx+tml5rhm7MZnj9dR3PPtsbm03kvvscTJzoQEtj7IYbnOTmOnj+eR19+xo44QSZlSuVAfS005wsWmRDwznrg8WLRcaP70RtrQ6jUWbaNBv3369thVULrCYz5Od34uabw11Oq2uvrcFotDB/fiL33Sfw6qsWRFH5bg0GA6mpqaxevdolCu4N1bqmpqa2K+t6zJNVluWgpp3uFThali82Ntbj77///pthw4a18d01L2QZliwRmThRYu9egZEjK3nuOQcnnxwV8LinnnIwf75ISYnIypUCmZlO5s+3MXx40+fcsEFZl+7cqaxLL7mkmg8/NOKl1uIBNbSi16dz5ZUCGzcqzq3TT7fzxBMb6dMnmYiISCyWPcyfr8xEXn21MSSj0+lISkqiqKhIU7LGYDAQHx/vsq5GtcynjdHuyeo9qrnXnqpOFqvVyo4dOzymooHK34xGow/5DrfpcGuPti19vg0blHXp77+LnHSSk6VLbSQn79TsLO6O998XefBBidpaAZABgWHD5CaJWlmppCT++KMyRx04UObll/eSmSljMGRqHqM2IrNanUydmsFnn5lwOgXS0228/76VIUOcbNliRRRFJEnHI4/UEB1dxttvJyIIcPfdjaGhzMxM1q5dq5nxBEq+8fr160lJSTnkxGp769ruyQqwevVqzGYzoE0+UMqhvKed7eEDbk60xDS4pASeflpi3jwRkwlee83GLbc4kSTYuNH/ILF6NVx/vZ69e5Wp54032nnlFQe9ehn4+GORnj11TJrkG1JzOmHSJEXCxW4XSEmRee89G+ecI7N3rw1RDGBSgXnzDEya1I2GBpGICJmnn65m0KA1DBw4EPAskcvKyuTKK9cSFTWcl182Ul6ewZQpdSQmeoqCp6Wl+aRWqrKmJSUl7ca6dgiy9u/fP6DlKy4uJjk5uU26dHfUAcFqhbfe0vHMMzrq6+GuuxxMmuTAXQVFa3A4eFDx1P76q0KIU0918vHHNg7VfpOTY6VXLwNTpujo0UNm9OjGdefHH4vcc49EVZWA0SgzZYqNBx5o3B6o3eNff8FNNxldTqurrqpg9mwJSRLYsiXcpbfkTla9Xk98vIk77jjAtm2d+OqrRJYti6ekpA5RbBQFT01N1cxa6tSpExs2bHD14W3r77r95VRpQFVo8Ie2yA9u7ZBRczmYZBmWLhUZOFDPxIkSp5/uZPVqGy+84ElU9/OCYhEffFBHVpaBX3/VkZ4O339v5fffG4kKkJoKP/2keGBvukkiJwc2bYJ+/fSMGSNRXQ1XX+2guNjqQVTl2nyTIoqK4Jxz9Jx5poH9+wVOP93Ob79tZ+rUUpd3uVOnTuzfv//QdXqSrq4um3PPTeXbbyVAxmoVefBBI04nLl9EaWmpJlnV5VJpaWm7qMjpEGRtCm2VH9zRyuS2bBG4+GI9V1yhR6eDr76y8uWXdnr21L4P9f4++kgkJcXAm29K6PXw/PM2du+2csYZ2uc55RRYsMCO0wkjRhgYNMhAbq7AySfLbNtm5f337XhHRNQQWENDA2VlZezfX8CNNzbQtauBP/8USU218MYba3j33Z2YTHYPYkVFRbmaaqukq62Fq68OY+TIePbti6BXLyuffLKW8eNrefddg4uwKtHtdrumZ7hTp07k5eXR0NDQ5t93h5gGN4XjoUzuaAaH8nJFw+idd0Sio+Gll+z8618ONDLuPLB9exjXXhvH7t2iK4No9mxHQE8tKFZ49WrXlQMyc+eWMHKk0nBqy5bGsje73e6y3jabjYiICBYujODll7NpaBAJD3cyZUoNt97qRBC6k5OTQ1JSkuaUVSGdk6lTw3n9dSN2u0BiopOXXqqkc+dNSJLEY4/Vo9freeUVZQ360ksQGRnJwYMHNR1NYWFhLiWKtl67dgiytkdV/o5gWUtK4PbbJVasEKmpgdtuc/L443YSEgIfV1kJN9wg8dNPSs8ZNYMoPV2puPFX8maxWPjxx1hefvlEamv16PVOoqKsVFQYmTEjlhEjajykV1X1fJWsH3+8j0mTupKfLyGKMrfcYue11xxIkkIQp9NJfHw8tbW1PoqFJpOJefOqeemlgdTW6jEYZB57zMyECTZAz8aN0iFFEZGnnlKmtCphp03LZsuWzUT6Kbzt1KkTGzdupKamhgEDBrTZ2rVDkLUpHE8yK8Hi++9FLr9cwukU6NfPyfz5dk46yfN6vVMfzWYr06bFsHBhPA6HgMlkZtKkrZxySiV5eZCfL3rkF7vXmu7eHcYtt8SSm6tY4dGjHcydq0x3hw93smpVGA880IXPPvOtiCkqUjzLf/55AgDDhjn58EPfZApZlsnIyHDJj6rYulXkxhvDyM1V8h0vu8zG7Nlm3MX2VWeRMjjAU08pyf+vvmoETFx3nV6zWgdwWVRVJ7qtrOsxQdaQg0mBLMts3uxg4kQ9P/2kJzHRRlmZRG6uTH39VtaubXDFoL2T/n/6KYnp09UMIieTJ9dy/vkbOPHEE4mKOsnvOaurFUeSokAIAwYoVtg9K/K332z06GFg6VKRyZN1TJ+uzIK8k/xTUy3MnVvHuef6l5YICwtDkiTq6+vR6xO45ZYwvv9eeYz79nUwfvx/ufLKgT4ZY7GxsTidTsxmM1FRUQgCPP20lYMHBd5918Affwxm5sz/+j1vUlKSS22irTzDxwRZj2UHk9pm0Ww2Y7PZOHDggIfMiZqXXFMj8f77nVmyJIOwMCf335/PTTdV8eabacyda+LOO0/it99qfRJANm1SMohyc5UMoquvdjBnjp2wMANr1vgKjKlwOuGJJ3S8+qoSL01Kkpk718YFF/h+JpIEq1ZZOeEEAy+9pIR0rFZ45BGJ+vrGJP+hQzc2mYQBEBYWwbRpkXz6adShczuZNcvM+ec7+PNPW8C837y8PHr16sXBgzBjhpEPP1SmzF27KgOYP1FwtSVHRUVFm61dOwRZ22PlzZGeU6vqRqsAQB0I3Ntu2O12HA4HERERxMbGHipfM7BgQRhTp0pUVMC4cU6efNJBcnIikMgbb8CKFU42bNDz4ovRPPGEMqhVVysZRN9/r1jEk0+WWbzYhjtX/MWRP/9c5M47JSorBQwGmaeftvHII4E/i/h4+OMPKwMHGrj9dgkQvNalsGGDdj2rOz77TOLee4dQWythMMg8/riZhx9uFHnX6/UUFRWRmZmpKY5WVVXPG28IPP98JNXVMHas7VA9bDEHD5o8RMHd4XA4iI+PZ//+/bzwwgu89dZbAa+zJdAhyNoU2tLBpNV6Q4t8WlPPw+l543Q6+fvvvz0sz88/C0yYILF5s8jIkU5efNFOv36+lu3nn2107mzg2Wd1XHCBg6++0vHaa4pFTE6WefddG+ed1/QsYcsWuPZaPdu3K1b4iiscvPuubxhGCyUlcPfd7u5nmT/+sHLKKZ736G9g3rQJrr02jB07lHP/4x9VTJq0g/79e3rsJwiCK683zT0ADKxcmcDcuSexe7eBM8+08+yzFnr3VgaZwkIHkZGRmM1mlyi4OxwOB+Hh4RQUFJCbm9v0DbcAjgmyNqdl1er4pkW+6upqKisr0el0LZJ7HAg7dwpMnKjjm290dO6sWMRLL1VaUGghLg4+/dTGJZfoOeMMA6BkED3+eGCLqFrW2lplXfrtt8oA0r+/zEcf2ejWrelrtdvh3nt1zJvXWHw+dKiDTz6RuPhiAzt2WFF5oZUUUVmpnPuHH9RzO5g8eR0jR2axaVOVpupDZmYmGzZscKn0b90qMmmSkZ9/HkC3bg6mTl3PHXdkYTA0Hufe0EoVBXeH0+lEkiTCw8PZtm1b0zfeAjgmyNqUZXWfegaqOfWeerqTLzIy0qPd4o4dO0hPT9fsUtYSEASBmhqRRx9V8mqNRpg61c7ddzuatGxbtsBDD6lftUB4uExhobXJ45xOmDrVyGuvGbDZlHXpO+/YuPDC4Nbq77wj+qxL1eLzyEgH8+ZJDBmiZ+NGm6uEzT1jyj2HODlZ5u23zZx5ppWNG83odDoyMjLIz893Sbiox6szlZ07K5k9O5X33tMTFQV33LGDqVNTKS01kp/vqRbhcDhc+cCqKLh771aHw4HRaKR79+5YLBZ+//13Ro4cGdTn0FzokGT1rrypqanBbrfT0NCgWfDtPfU0Go2H1W5RC63Z/Vxpjahj8uShVFXpuPFGJ08/bcdrlucDb4vYt69McbFMSYnArFmi33pRgC+/FLnttlOpqVEcME8+aePRR4ObvSh5vI3F5+7rUhWzZjnIzRX580+RSy6RWLrU7rJuH38scvfdEtXVygzgySdtPPywE7vdicPRmKyfmprKqlWr6NSpk2t2JYoiVit8+21PXnwxnIYGiZtvtjFpkpXdu/MwGFI9jlMdaOqxqmzp/v37OfHEE92+Awc6nY76+nr69OnDySefHNRn0ZzoEGQtLS1l165dHt3i3KeeanzMe+p5tAXfgdBaZP3vf5UWFBs2iPTpU8m330YwcGDg8zqd8PTTOl5+WedjEUtKoHt3A5MnS1xwgZWTvKIy27Yp69KtWxULd+mlNubNcxJMg/CSEuXYP//0L4rmjh9/tNGrl56ff9Zx330yp54azrhxka7a1kbPtLK/anXV6blal6quTx0OJytWJHLLLZHs2iVy6qkVzJhhY9CgcGRZZvdu5X10Oh0pKSkUFhaSmamU46lkBEhISGDv3r1YLBaX19fhcCCKIvX19SQkJBAdHd30B9LM6BBkjY+PJzo62m/ZW1FREXV1dST7eypaAGrHspbC7t1Ka8Qvv2xsjZievpaBA08PeNySJSL//rdERYWgaRGTk5W83WuukTj3XAP791uRJMUKjxkjsXRpoxWeNGk1F1/cS1MJ0B1a61J/omjuEEXIyVGcX7Nn65g9eyCg1LYuWuTpmYZGOVT3pPuMjAw2bNhAeXk6EydG8/vvCfTo4eDTT+s59dRa8vMPAH18EvUzMzNZs2YN6enpHlYZ8LCuJ5xwguucOp2Ouro6v5lOLY0OkcgvSVLAyptjKYOppgYee0zHgAEGfvhB5Ikn7GzYYOXKK/07kAByc2HgQD1XXy1RWQmXXeagqMiqOXW97DIn11zj4OBBgVGjJJ56SkdamoGlS3UkJMDnn9vIybGRnt7QpFPsnXdEkpMNvPuuRFgYvP66jdzcpokKygxg+nSlRE/JIYYpU2wsX65NVHVm5R5Sqq4O4403+jB8eCQbNkjcf/8eVqyo54ILHJhMcVgsFhoaGlyWUYUkScTHx1NSUgJ4WlZQkiAqKytd1Tbq9rq6Oo+1bGuiQ1jWph6YtgjdNPcA4XTCwoUiTz7ZdGtEd9TXw9ixEl9/rciF9umjWKVDBsEv3nvPwddfi/zyi8gvvyhW+Ikn7JoF41oIZl0aCB9+KHD//XrXunTYsBp+/TWKKVMk+vbdSmpqjY/fwWq1utT6rVaBOXP0PP+8kfr6SP7v/wp5+mmRmpoS9PrG9hdq1Yz7+lRFVlYWGzZsICUlxcfyCoJAZmYmeXl5dOvWzWMaHCLrUaCjW9bnnxf48kuJNWtEhgxpujUiNFql559X1qXx8TJz5ti55JKmP4fcXKURVF2dOgjKLF/uKyWqlRQRaF2qpPMFTvjYskXPU0/1Jj8/HEGACy8s5amnDhAdbWT2bJnXXz+JG2/sxYYNVSQne4a8ysrK2Lt3H//7XyK33hrJnj0iF1xgZ/p0C/X1uzAa0w/JyzQiMTGRPXv2kJqaqlmvGh0dTXl5uY/lBUhJSXE5olTL2tDQ0CzTYLPZzMiRI7FYLGzYsGEz8B9Zlp8MdMwxQda2yg0+WrLu3w9DhuipqBAxmWTmzbNx9dVOmnJKf/ONwO236ykvF9DrZSZPtjN5sqPJ42pr4eabG61w794yV1zhYMoUJea5e7fV4z2cTqert01dnYVHHolk8eIonE6BtDQLU6bsoGfPCnbutLNzp7bYnEqIujoDt98ex2+/KY/c4MFKrDYrKwZQvFyXX74cWe7JG29InHlmLFu22Dym/vv2mRg/PoJ160yceKKDzz+v59xzFctbXp7FgQMHCHdvpY5yTRkZGRQVFfnVWtq2bRs6nc7H8oqiSHp6Ovn5+c2+ZjUajfzyyy+H8pSFAcD/BEH4Tpblv/wdc0yQta1yg490gKirgxdfVFojKpct43DQJFHz88MYPFjPxo3KE3zxxQ7mz7cTFVh8UMMKO3n11RrOOacOi8XC//6XyC+/RHPxxXVMnbrNFXOura1l7dq1LF2aweuvd6ahQUd4uJPHHjvIzTfbMBjSMRq7BEz4cDrh4Yd1zJ6tw+EQSEtTtIDPOkv7Wl980cH27QI//qjjvPP0/PqrjaIieOopifffNxAXZ+TOOzcxbVq2x5Q7Pj6e3NxcTZ3f1NRUVq5cqRkTVyVjtbrPAaSlpbH6UHGuKIrNRlZBENyzpPSHfgKO/scMWTvCNNjpVPRxH3+8sTXi9Ol2Hn9c4uOPdYwfr+Ott3wHnfp6uOUWiS++ULw2J52kZC316KFsV62f1vTzxx8jmDatC9XVEpLk5Kab9nLrrQWEhRkoLVUs35w5pZx2WgS//GLimmv6c8MNIqIo8s47G3nhhYEa69ImRodD+OADkfvvb+xR8+yzNu65p+nvackSOwMGCKxYIdKnj56iIgGzGe66y8ZDDzWwYcMeIAP3x1dtPFVdXe3zfjqdjri4OBoaGjTPl52dzdq1azWLFtRmy/v370cQBOrr6zXlS48EDoeDU5R8yxLgTVmW/w60f4cga3t1MB3OOVeuVKQ+V64UGTjQycKFNk4/XSH7u+/aWbpUZN48Hffd5+CEExqTPp57Ts+rr0ZjswlERVl54on9DB9eSnm5heXLZde1uNeZGo1GSkpiueOOZLZtkw7p8TqYN89OVFQ6kO5zfb/+amfAAAN33RXBySdbue8+PX/+qSTuNhUv9cbatUp9qqIwIXP99XZmzWpaYUKFIMDEiQ5uvllg506RQYMczJvnIDvbgiwLHsn67oiJiaGoqMjHswuK5d2xY4fmOjw6OhpZlqmrq9OMn6anp7N7926cTicNDQ3N5mDS6XSsW7cOQRAygS8EQegjy/Imf/t3CLI2hfZsWfPz4bHHJBYt0pGS4uS116q59NIq7HYru3c3WsBHH41i8uReXHihgw8++JvVq5N45pmeVFTokSSZO+44yKhRqxg4cAAGQ6Lf/jb19TBunMRXXynrUm8r7A8nnADPPWfnoYckBg9W8odTUhr4+GMxqDAMaCsfLlpka9Kj7Y41a5QkkD//FOna1Unnzk6+/VZJerFaFaeXXq+noKCAjIwMD+LJskx0dLSmeLcoioSFhVFWVkZSUpLPefV6PXl5eZzknSWCEuaRJInCwkLq6+ubPc4qy3KlIAi/AhcCxzZZ26rqxm63U1tbq+n5rK62sXBhMh991AmHA66/fj9jxxZiMknU1Bg9VBaMRiO9ehn4+muZv/6K4JJLzqSuTjnPRRcp69KYmEiWL9cTq9XD4hCmTdPx3HPKutRkkpk9286llwY3iL39tsgTTyilawB9+jiYOfMvhg4NnIQByvT+oYd0zJmjrEvT02XmzfMvqKaFggJ47rkT+eEHpa/OG2/YGDeusa+OOjCqyf6qF9e9S7ksy8TFxVFQUEB6eroHkZ1OJ3FxceTl5WmSVfX0NjQ0+DipQCFzfn5+s3mDS0tL0ev1xMXFIQhCOHAe8FygY44JsjZXGEVN+Pe3/lN/AFeAvq6uzmP6GRMTy3ffRTNtWjQHDui4/HIHzzxjp0uXFCDF77nNZkhOlgGZujro1k3m009tPumAWvD2Dk+caOeJJ5r2DoNvc+PrrrPzzTc6Nm0S+d//Eji9Ca66K/KHhSnJ+uPHBz/LUfu9vviiDqs1hQcecDBhgnZfHbUkURRFMjMz2bFjhwdZnU4ner2emJgYHyKrifgGg4GqqiqfQU/NWsrLy6OH1zREjcGaTCbi4+ObhayFhYWMGTNGNTI5wCeyLH8T6JgOQdam1qz+tsuy7FJaCERC74R/7/Wft9q/IAiUlJRQXV1N9+7dXefzbo343ntWRo5sehB54QWRadMkLBYBQZCRZUWpvimi7tqlxEtV7/A//+ng/feb9g5D4ObGGzY4GDLEwNSpfbj5Zhtuz7wLOTlw442eivxvvhn8ulSW4dNPRSZPlsjLE7j0UgejR6/kyisHNnGcsuaMjIw8VL7XWHuqlrFlZmaSm5vrQ2SdTkenTp38FpgnJiayd+9e1L6v3sdmZGTwww8/cM011wR3kwHQr18/1q5dq/7ZJ9C+KjoEWd3hT2nBbDazceNGV/8b1dKq7efdSRgdHe1ByCNJ+HcP3RQVwRNPSCxcKJKUBG+9ZWPMmMYpnD/8/LPAuHF6SkoEJEnmwQftTJ3q4KSTFBJ9843AxRf7kl31Dn/5pbIu7dVLyVpyKxLxi2CaG/frB1OmKF7qs8/Ws2FDYwFFWRlcd53SEwdgyBAnixfbmqwAckdOjrIu/esvkf79ncyda2P4cDs5OYGFtNXP3L1fzYEDB1zVMar1jIyMRKfTUVNT43IYqUSOiYnB4XB4ZCKpA4Aak3Xv+6q+r06nIzw8nKqqKpYvX07fvn2Dv+FmQoch6+rVq7FYLH4D75IkkZ2djdFoRK/XH1a525FAEJRwwgsvKOtEi4WArRHdsW+fkgW0Zo1i1c4/38EHH9hRpX++/NLGyScbGDdOT2Gh1SOe+OyzIs8+K2G1KuvSWbPsXHZZ8OvSiROV+tLISJnnn1f62mjh4YedLFpUzZYtsdx3n46XX3Zw//065s5V1qUZGTLz51sZMSKoUwNw4AA8/rjqbJOZPdvGjTcqg5rd3nQrElmWPcgaHx/P7t27XZbQ3dOblZXl4TByP867BM59m3v5nDqIu2c3xcfHh4rPm8KgQYMCrkv37t2rKXTVEpBlWLYsgiefTCQ/X+Liix3MmOGge/fAU16rFW67TeKTTxSLeMIJMh9+aMN7RnbiiTBmjIP58yVuvVXH/PkOli0TGDNmGFVVevR6mUcesfPkk0e2Lr3tNqWJVFMTipdfXsfo0SOZPVshqd2uFK5Pm2bjrruCX5fW1cHLLysle0qShJ0JExy4R0m02ldowbs6Jj09nYKCAjp37uyxLS4ujl27drnK3AKR3H2bWifrLrqmToPV188++2zQ996c6DBkbS9obI1oolu3BpYutXLOOU2vS195ReTppyXMZoGYGJmZM+1cc43/B/6ttxx8+aWOxYt1/Pe/IoWFisVo9A43fa2B1qXBYOfOSIxGmdpaEbtdcT69/XbwyfreSSCjRzuYNs3OIWEHr32drhRO99ad7r/1esUb7m6B3S2hN5HVaXK3bt38buvatatPXDY9PZ3Vq1e7RNfct9tsNo/1bGvimCGr93qmufHTT/D55xLz5yutEWfMqOGss3bQv3/vgMf9+iuMG2egqEhAp5O55x47M2Y0bRGtVujWzcnq1QpRTzxRZsKEHK67bkCT1xrMujQQSkqUpIY//lCSIsLCnJjNIps2iUhSsFU5yro0J0fk5JOdzJ1bxymnNGCxWMjP93T2WSwWl99hxYoV6PV6l09B/R0ZGcmOHTsICwvz+I51Oh2JiYmUlJT4fP/Jycnk5OSQnZ3tQ0g1ST87O9vnOEmSXM2W09PTPabBbdk5sMOQtanwjHcBcXPBaoVrr9WxdKmi5D5+vNIaUaczs2+ff8uYn694anNylC/2nHMcfPihnbi4ps/53HMi06cr61K1QfFVVzno3Lm+yWMPZ13qDbsd7r9fx7vvKiRPTjazeLHA6acrBQfr14s89piOadMcLlE5bwu4b5/Mq6+m8MMPiSQkWHjkkd2cf34xYWEG9u1r9DGoPWTUv61WKzt37qR///4e37P62uFwkJqaSnFxsU8xfEZGBps3byY8PNzj+1entEVFRT7PhiiKpKSkUFBQQFxcnE+CSWZmJuvXryctLc01DW7rDgwdhqxNobmzmGQZvv1W5JFHdOzcKQIyoghTpyoCZTU12oOH1Qp33KHjo490yLJA165KK4hgJHu+/17g1lv1lJYq3uGHHlLWdpmZBqZP19G/v/+v60jXpSrmzhV45BGlbC483MnkyQcZNGgVSUkZbN5s5cUX7Vx88Sm8+KKO7Ox19O3rGV+22YzMm5fOO+8o3rUJEyw89JCTqKjOQGePpIbGz7jxtZqU4N6oChrDcqIokpSUxP79+32SGsLCwggLC8NisfgM1unp6axdu5aIiAifbWorjujoaJ9tajvIsrIyjzWr6jU+GuTl5XHTTTdRXFyMIAhs2bLlXlmWX2vquGOGrM2ZxbRlizKF+/lnkZ49nXz1lZXffxd58UWJ8eN1zJ3r0Bwc3nxTiRuazQJRUTKvvmrjhhuaHkD27FG8w+vWKQ/BhRc6WLCgcV36/PN27rtPz6RJ/fnnPz2P9V6XjhihDA7JyQoZrFabpgVUX69da2Tq1JMoLg5HEGQuvbSYyZMLiIgwcOBAY0f5zp0NfPqphVGjwnn00VPYu9dMeLiM0wkffqjjqaf0FBWJXHmljaeftpKVpVpEz+v1JqL6t9FoJCIigpqaGuI0ph9qx/vo6GjNhPysrCw2btzoQzo1S6i6utpnm6oWUV5erpm6mZWVxZYtW0hNTW1WyypJEi+99BIDBw6kpqaGmJiYuwRB+FGW5S0Bj2uWs7cDNIdlDdQa8dxzHcyapWPRIh0vvujAYGi0rH/8AWPHGsjPV9ald95p58UXm16Xms1w660Sn32meId79FBqPL2LwP/9bydvveVk+/ZoPvnExqWXmqmvt/Lgg2EsXhyJ0ymQnm5h2rRd9OhRwc6dVnbuVI51F5ZzX/9VVBi5+24TK1YoD+nppzv44AMrKSnRgCKcXVRUhMlkchHqrLOc/PvfNmbPNnDRRQamTLHw6KNG1q7VMXiwgw8/bGDIEFVO9PCtT1aWUpOqRVZQBuTY2Fjy8vJ81o5q/NRqtfqkC2ZlZZGTk6NJyMzMTNatW6eZghgeHk54eDh1dXVERUV5CKgdDdLS0lwC5IfiwFtRyoiODbIGU3lzpGS12WDOHB3Tpun8tkYURXjsMTuPPqrn1lslPvhApKRExxln6Pn7b+XazjhDWZdqZfx4wz1rKTZW5uWX67n4YqVL2YEDvlbwySd13HDDUP71Lx05OSW8804XzGa1vrSMMWOsGAxJGAwZrjizmp6nQsnogvvu07NggR6nUyAry8m8eWaGDFHFyJR91XrL6upqNX8VUAawpUslVq7UceGFEaSnO3n3XQtXXaUOTkc+RYyJiaG2tpbKyspDswLfbLO6ujr0ej3V1dU+KYNGo5Hi4mKf/4eHh7uyneLj4z22hYWFYTAY/HY2V9s9xsbGtoj+0t69ewFOBgKWx0EHImtTONIC9O+/F5kwQcf27SLnnOPkhRd8WyOquP9+J88/L/PNNyLXXx/Ot98OQpYFOndWSt4GD27cVyvN0WKx8NtvRh5/vBMVFXp0OifXXLOf227bQ1iYgby8xqwqNclf/fuUU/ScckoJq1al8Prr3Q/Vl9p46SUrkhSBLDdaE6fTqTlwzZ2rZ/Jkg8v59OyzFm65Rf3MfKcBGRkZFBYWuoq2a2rghRf0h8JIMhkZTtautRAoVdY95VNLYF3NOFOv1263s3XrVg9B9ejoaBITE13ZZhUVFeTl5fmQUqfTUVFRgd1u98lKUxPxvckKSm/X4uJizetXUxnVgvzmJGttbS1XXHEFwH2yLPsW4nrhuCXr9u0CEybo+P57Hd27O/nsMxsXXeRfQVB96IYMkVm2LIylSw2EhdmZOLGIiy8uxmq18tdfVtdD515jajQaKS2NYPz47mzerMTozjvPzsKFNmJikoFkTceL+ru4GG680cCqVcmoYgKfftrA+ec3Kth7/3bHihUiY8caOHBAdJH85ZdtTTqfTCYT27dvx2p1sGiRgaeeMlBSInDttTbuu6+ezp0tWCxWqqu1uxx4p3y6p3iqA5H6P3U9abfbycnJoWfPnn5nU0o/2N0+01JZlklJSdGsddXpdFitVs2qGlU5U6vHDShT1dLSUiIiIpqtPM5ms3HFFVdw/fXXs3r16s+DOeaYIWuw0+CKCkXiZPZsHRERMH26mbFjawErhYW+00/1odu4MYbp03tTUhKGSpjXXvubf/yjE3p9huuhc3dEyLKM2Qy3367ns88kV9bSggVm+vRRrlVrfFEfUodD4L77DCxYIB1al1oZN+4g06encOedYezebQ54r8XFcP31RlasUIgwfLiDhQstmkkR3skI6s/GjQn8+9+wc6eRk06q4skncznxxBosFj179viuhY+0w4EKSZIwmUx+607VfdLS0lwJD+7IyMhg7dq1PrWu0LgmPsFL+tHpdJKYmMj+/fs161n1ej319fXU1NQ0i2WVZZlbbrmFXr168cADDwR9XIcha7BrVn99berqLHzySSxz5mRQU6PjoosKGDduD8nJkJ/vG4BXH7iKCgM33BDm8rYOG+Zg3Dgbt95q5KWX+nLddY0PpGp9Vbzyip7p0w1YLErW0muvWbjySsehewn8IL/zjo5JkxqnrDNmWLj22lq2bNnBqlWJfP+9xOOP65k61eZzrN2udGz74APp0LrUzuuvl9G7dy319VZycz2noOrn6+6MKiqK5JVXsvjhhyhSUszMnVvL1VdLiGLgJJDmgFr+5o+sgiCQlJTEunXr6NKli8egoGY5eZfIQWNVjfc02eFwEBUVRUVFBWaz2UfHyeFwkJ6ezo4dO5rFsv75558sXLiQvn37MmDAANavX78OmCTL8reBjuswZLXZbFRXV2umoamOB0EQ2L17t0+Sf05OLM88k8i2bXpGjLDzwgsW+vdPQJaV9Yv3tFN1xDzwgIH339d2xDz7rIHc3Gj++quW004TPQaTn38Wue02A8XFIjqdzL332pg2zRaUAybwlDUch8PBrFnlnHRSMq++KnHZZXnExzf2+Fm0KIHZs7tgsegIC7Nzzz17GD26EoPBgM3WqDaoDkze4ulVVfD883reektCr4cnn7Ry+uk5DBrUG1EMordjMyAqKgqn0+lXo1cURfR6vd/Wjlq1rupxaWlpPs2W1Yoc1fK6lz2CQlaTycQff/zht7n04WD48OHeYaABwRzXYchaXV1NXl6eh/WLi4tz/V1QUOCqOVQ/iJ07BR58UM/SpRKdOzv54IMGRo1yIAiK9fGG+tDOm6fn0UcN1NUp3c9mzPB1xLzzjpWzzzbyr38Z2bRJebN9+5Rp59q1yj7nnmvn/fetTWYtORwO8vKsjBsXxcqVypr2lFPqmDFjF1FR9axbZ3VZbKvVyu7dm3j44S5MndqFW29N57vvCti6NYl7702goECHTifzr3/ZeOEFG5KUgRIVCAy7Hd5/X2LKFD3l5XDDDQ6efNJGWprMgQNKXxj3rmstDTV317sQXIVOpyMtLY1t27a5Wjuq0Kp1VZGWlsaaNWs8mi2r2U1JSUns27fPx/KqZK6rq2Pr1q0tcLfBocOQ9bfffuPcc88lLCxMMxtG/TAdDsch62Bg1iw9RiM8/bSVu+6yERYGguB/+pmTI3LTTQaXmt/YsTZee03bETNkiJPevW1s3qznm2+cfP65xCef6JBl6N5dZsGCBnr0MGOxWCgp0Ra+djgc2O3w2msn8t13aciyEi996aUihgxxYjAkYDCkubyggiDQ0NDA1q1bmTgxlS+/dLJxo5FBgzpTUaEmRdj54ANrUOEjFb/+KvLIIwY2bxYZNszBc89ZOfnkxs82JSWF1atX07lz51bLi01KSmL37t2a4megWMnw8HCMRqNmGMe71lWFmghRWlpKSoqi3KGeQ7W8+fn5Hk2r1e29e/dmyZIlfqVfWhodhqw//PADJpOJ0047zcPzqb6Oj49n374Cli3T88ILBjfrYD1UGO3/ISsrgxtuMPDHH8pDcdppDj780EKKmwqLd5Nli8XCjBkyl1zShauvNgICkZF27rlnG+eeW4LVqiM317PuNioqyuPvefOMPuvSm292ANprNVBihor2k5Vu3fRs3ChSUSHSqZOTDz6wcMopwWfZ7NghMGmSnm+/lcjOdvLhhxYuvdTh4xHX6/VERkZqkqKlIIoiycnJFBcXk57uq8YICvHS09M1wzhqGZzNZtPM+928eTPJyck+BSBqxU1WVpbrf2oiv9ls5qGHHmoTokIHIuvw4cNZsWIFw4YN09xuMBgZNWoodXV6TjvNwZdfeloHLTgcMg88IPHeewacToHUVBsvvlhA375VlJdbKSxUiKl6mXU6nWvavWaNiUmTOqEOAoMHW/jxRyt6fVegq/+TcuShFBU//NCTf/wjBrNZyVkG6NnTGTRRKypgxgw9s2dLhIfDlClW7rrLHrC5slo32lpkBcWzu3HjRr9kFUXRFcYxmz0942qta35+PpmZmR6EVXOJq6qqiIuL86iqUat4vNfCap+bzlr1fa2EDkPWESNGMH78eL/bBQEGDarlv/81MXNmHZ06mTl40L/o2XffpfDGGz1oaJAwGh3cc08BY8ZUH7J6nn1e3b/ovDxlXbp6tfLl9ulTz6ZN4ezaZUCvDxznPZxQihZ+/13kllsMFBREIIoy//63si494YQwfvxRx++/i4wc6T98ZbfDu+9KTJ+u5+BBGDvWweOPWz1mEP6gKt63ZBmiN9Tsourqar/CAu5hHO/rSk1NZfXq1Zp9brKysti3bx9xcXEeifrQWHHjvRaur69vM6sKHYSsTqeTqKgo8vLy+OKLL+jfv78ryK2S0Gazcfvt8OefZ/Dss7VMmJDn4RFW2yRs3hzGrbfGuIS+brjBxuuv2zAYEoAEv9dgtcKdd+pZvFiJl3brJrNwoZlevaxcdFEDK1Yk8MYbEuPH+3qulP6lele8VKkvtTJ4cHDpkfn5CslzcpQH7swz7dx771+cdVZfRFHPJ59YGDkyjBtuMLJ3b4NmTvKPP4o8+qiBrVtFRo50MGOGlf79g58yq4r37mu91oAqz9K7t3bISA3jrF271qcoXKfTkZCQQGlpqQ9ZY2JisNls1NfX+zSlUrOmDh48SIJbzml9fb1m0kRroUP0Z73hhhu46aabsFgsLFmyhKqqKpfrPjs7m969ezN48GAuuKA/555byvffp5GV1ZeePXvSpUuXQw1zE7n22kTOPTeWvXsFTj3VQW5uA3Pm2JpU5HvjDYm0tHAWLdITFQVz51rYsMFM//7KFztx4lZEUWbKFD3eeRnvvKMjLS2c+fP1hIfD669b2LrVHBRRFRkYPT17hpOTo6NLF5k//jCzdKmVE04wUVpaCsApp8hcdpmD8nKBhx7yrPXctk3g8suNXHZZGBYLLF5s4dtvLYdFVBXqVLg1YTKZqKmp8eh67w41jBMXF+cR41ahNqXScoypoRpvywqNOk3uaAmB78NBhyDrRx99xPfff8/48eMZNmwYJ598Munp6a528Uaj0aW4fuml+6irE1iwQJk0OJ0wYYKezp3D+e03ibQ0me++s/Drr5YmFfl+/12kW7cwHnnEgM2m9FopKGjg2ms9idalSyyXXFJLXZ3A1KkKWVasEDnxxDDuu8+I2Qy33KIcqziQmoY6QHz0kZ7ISJgzx8KmTWYGDlRIlpKS4pHPOm+elchImbfflti1S6kgeughPaeeGsZff4k884yVVavMXHKJrwMpWERGRmK327FYLEf2BkcAd50lf9DpdCQlJfmsW0GZSuv1es1U1MTERCorK3E4HD5kjoiIQJIkqqqqXP9rDrLefPPNJCcn08e7tCoIdAiyqhgxYgQrVqwIuM8ppwgMGWJh9myJDz8USU8P58039ej18NxzVnbuNAdc14Ey7TzjDCP/+IeRoiKBs86ys3dvA88/b9OcYiYmJvLAA9vR62VeeknirLOMnHuukbw8geHDHeza1cDMmcE5kP77X5Hu3RsHiDvvtFFY2MANN3g+bBEREa7pP4DBALNmWZBlGD48jP79w5kzR2LcODsbNjRw7712mqG6y6W80JpIS0ujsLDQp57U4XDQ0NBATU0N9fX16HQ6D3KpSE5O1qyBVQcCm82maXlV66pOkZuDrGPHjmXZsmVHdGyHWLOq6NevH1u2bAmog2MymRg5sowXXsjgX/8yIghw7bU23nqr6emu1Qp33aVn0SLJpfKwcKGFAQMCTxljYmJoaNiCyeSkpETHypW6I1qXXnedkVWrlAfjrLPsLFhgRaNIxIXk5GRKS0tJT09HliEyEnQ6qK4WOeEEJ8uWWejTp3mlSFJTU1mzZg2dOnVqkZirv64ITqeTVatWAbispCiKHo5AtVO5Vomc0ujZN5UwLS2NnTt3ajrOYmNjPWYRzdHuceTIkWpZ3GGjQ5FVkiRSU1Nd7ngtmEwmRozYzgsvpGMwwIYNDWRlNf3eb76p44knDJjNAtHRMi+/bOG664Kbsr7zjsSjj56J2axDDaXk5JiDUsb3dFwFP0CAQtbc3FwqKzOYONHAzz/ryM52kpHh4PvvrZqzgKOFXq8nPDxcVTgI6hin0+lTJOD9Wl1vanVFiIiIIDs7m6KiIpdzUWugUJQv1vpU4zidTqKjozVTCdUmylppi6Cs01VyNTQ0tKmDqUORFWDYsGH8+eefXH311Zrbw8PDiYw088QTVqZMMXLwoOCSGNFCYzhEPKTyYOO557Snu97wjpdeeWU5J5wQzTPPGLjjDgMLFwZWmH/tNSW9Tx0gXnnF7LMeDoSGhkheeCGbr78OIyZGmeb/61/2oFtYHCnU/FpV6EyLfO5TdG9hdqPR6NGUyz1DKxAKCgqw2Wx+Oyiog7l3NY4aTSgrK9PMiFJrXb1DNaDMmux2u8tz3JYOpg5H1hEjRvDxxx/7JSso9YfXXlvBSy+l8NZbeubM8SVNYaEy7Vy5sjEcsnBh4GmnCu946bBhDubPr+XAgY0MGnQqb7yh58svdZSUoBlD/fVXkVtvNVBUpAwQ48fbePbZ4AYIUKzx7NkSM2boqa3N4Prrq5g+3eChbHEkUMvkmmrOJcsy9fX11NXVYTQaPep23fsCeRcJHC0yMzPJz8/3sY4qBEEgOTnZpxrH4XAgSRIpKUqOs/usTJZldDodERERVFRU+BSnO51OIiMjXV7jI2m10lzocGQdMmQIEyZMaHLdarUe5IYbEpg3T2LKFFyBf6tVKR/76CMl5tmli5MFCywuL2sgaMVL58+3HqrE0ZOfL2C3W5kyRcd99xm59VYDX33VOFDk5SkDxJo1ykN0zjnKujQYeVJQFRd1PPqonl27RM47z8ETT1QhittISBjg5xjZ5cHVIp93byB3zV6VgGqapNqaRBAEcnNziYuLIznYjI5mQEpKCn///Tddu3bVTMxwr8YpKSkhNTUVaEzUT0lJYd26dR61ruq2Tp06sWvXLh+yquVzJSUlzVJxczTocGRVa03VVDEtmEwmtm3bxr//bWfOHD3vvSfx6KN2Zs3S8fjjBhoaFPXBl16y+HhZ/eHtt5X60oYG7zzeRiQmJlJeXs5tt6UxfbrML7/o2LMH0tLgzjsNXon+Spw2WGzcKDBxooHfftPRs6eT//ynnpEjFc2mzZtr2bNnDw6Hw8cCgqdomnuSiHuZ3OFmJaWlpbFr165WJataGeNORG94V+NAIyHVdpDuyQ5qQkRUVJSr+7n7VFe1prIsN8u9Xnvttfz222+UlZWRmZnJ008/zS233BLUsR2OrABDhw7l77//5oILLtDcHh4ejtlsZsAAJxdc4GDmTD1z50quaaeapnck69JAebxqcXNaWhovv2zhxhuNDB9upKFBdBWgv/KKxW/bDC3dpsJCBzNnJrFkSSJRUXbuvnsnF19cQFiYjl27GvWaGhoaSEpK8rCKLZkWGB0d7brW1mwnoQp6+yOrWo3j3ofVffqq1rqqZHVPiFCzpdwrddQ1bpcuXVi5cqVmAsXhYNGiRUd8bKuS9dNPP+Wpp55i69atrFy5kkGDBmnut2zZMu69914cDge33norEydO9Ng+YsQIli9f7pesoDxMNTU1ZGQYqK7WUV0NZ5yhTDuDKR87kjzeqKgoamtrkWXZ1UmuslKHKMrcfnsNEyaUY7db2bXL0xnjXiigWj4wsnhxMrNnJ9DQIHDbbWYmTbKTmJgFeLq3a2tr2blzp19lhZaCGnN1L+RuaYSHhyNJkkc7R294V+O4h2W8a129O8S5N6yCRrImJiYiiiJfffUVl19+eevcrPd9tebJ+vTpw+eff87tt9/udx+Hw8Fdd93Fjz/+SGZmJoMHD2bUqFEe2jjDhg3jlVdeCXguk8lERUUFr7wSw8KFEhkZMkuXWpvM3gkmj9fpdPoNQ+TnS5x8soMdO1SFA5nU1AZuuy2XurpGpQb38IT7SC3L8PXXOiZN0rNnj8g//uHgmWcs9OghA9ojelRUFGazWVPVryWRmprKunXrWpWs0Fir2qtXL83tajXOnj17XIOhd6K+WuvqrbafkZFBfn6+q9BeJXpDQwNdunQJaCBaGq1KVn8frjtWrlxJ9+7dXc1sr7nmGpYsWeJB1qSkJGpqagKKLptMJnJzc8nOzubll63cfbeR5ctFhg3TnoI6nU5mzRJ48slwGhpEwsOdTJhQwKhRinLhqlWNag3ewXij0YhOF8nTT3fl66/DXcJoCxeaGTfOyNat4ZSU9OaMMwKHZdatU9alf/yho1cvJ199Zeacc4IL5SQlJVFaWqoZK2wpqPceyMq1BBISEti5c2fAwUmv17vCON4JD+61rt7b1IZVnTp1cnV5MBgM1NXVERMT0+y6wYeDdrdmzc/PJ8stiyEzM5O///bUPxYEgYEDB7JmzRpOO+00zfdR13Fms5l//tPG448n8tJLdjIy9ntYQpvNxqZNsUyb1pvS0rBD8dIynniijMhIIwZDikeHdC0P9IsvSjzzjN61Lr333i1MnKiMzO++a+H008O54w4DW7ZoqxEWFcGUKQYWLNCRkACvvWZl7Fh70PWtoDxku3btalWyQmMqYGuSVRAE13mz/GS8uIdxYmNjPQjpXusaExPjYXV1Op1H0bs6Da6trW1TokILkPXcc8/VzB2dPn06l1566VG/v9PppLy8nOzsbD799FMOHjxI7969fUSjAcxmMxs2bCAyMpL/+z+YPz+ZsrIounZVYoGVlUbGjInir7+816URQNNTux9+EPnXv4yuRlIPPGDj6adtrF17ELM5jbCwMPr3h8GDneTkiHz5pejRpdxshtdfl3jxRT0WC9x9t51HHrEFHcpxR1RUlKuxU2tOhRMTE9m1a1er1rmCZw9VrQHUPYxTU1OjWeu6atUqIiMjNRtWuXeQUwvP2zIhAlqArD/99NNRHZ+RkUFeXp7r7wMHDpCR0Sj49fPPP/Piiy8SERFBfn4+J510EgMHDtRU7FOnQJ06deKRRwTefx8+/zyNKVNsR1VfumeP4nxav175ki+4wM78+VZXI6mEhATKy8td1z1vnoU+fcK5/34Dl11mRpbhiy90PPaYnn37RC65xM706Ta6dTu6PN7ExETKysr8ekpbAqIoupo7taaDS21SpZXIoEIN4xQVFfkQUnUaVVRU+Hiz3etZVcvqT2mxNdHuqm4GDx7Mjh072LNnD1arlcWLFzNq1CjX9vPOO4/vv/+ezz77jIaGBsaNG0daWhoJCQmu4L060qpOJoDMTJnLL3cwa5ZESsqR1ZeazYoyft++4axfr3SY+/vvBj7/3OrRiVwljYouXeDssx2UlIg89pjE+ecbufFGI9HRsHSpmcWLrUdNVPAtm2stqOmHrQ3VUeQPahhHEATq6up8tmdkZFBeXq45I1DDOCpZmyOJH5RIR8+ePenevTszZsw4rGNblaxffPEFmZmZrFixgn/+858uz1pBQQEXXXQRoLjd33jjDS644AJ69erFVVddpakSIIoiPXv2ZPv27X7PFxERQX19vSs5ICFBxmoVjqi+9PnnlfrSzz+XiImBBQssrFlj8en4Bkp4oL6+3qNDwDPP2ACZV14xsGOHyJtvWli+3MyZZzZfT1l1KtxcrS+DRXR0NGaz2W+BeEtBrYrRqmOVZRmLxUJDQwM6nU5zaRYWFoYkSZpNqdQQj9rztTnIqkY6vvvuO7Zs2cKiRYvYsiVg4zgPtKqD6fLLL9eMUaWnp/Ptt41i5BdddJGLvIEwbNgwVqxYodnyABRHQmRkpKtl34sv2vjwQ8WJ8+qrwSVFLFsmcvvtRsrK1AbHNp58MvCxgiAQFxdHZWUlYWHxzJwp8dJLSlF6ZKTM//7XgJ+ioaOCKr1SVlbWqtIrgiC4Yq7+HD7NAa0mV2FhYWzcuJHw8HCX30K9JjV1Ui0y14oexMbGUllZqXm+rKwstm7d6mr0fLTT4GAiHYHQ7rzBh4MRI0bw8ssvB0zXUqfCUVFRiCK8/rqNceOM/PCDyIUX+rdqu3Yp69KNGxVW/uMfdt57z3O6GwgJCYksXGjnrbfCOHBA5PLL7UydaqNLl5ZtdZ+SksKePXtalaygOGzWr19/RGT1lnjV+q3OUiRJ8igcMJlM7N27lx49ehAeHq5ZPGC32zlw4AD5+fkuoqjQ6/XY7XbNWleTyYTdbsfhcDSLgymYSEcgdGiy9u/fn82bNzeZ1L97927Xh3T55Q4mTXLy1lt6LrzQV56kvh5uu83AkiVKHu+JJyr1pUEOfoAiFj5hQiYrV+ro39/Ju++aGT68+aa7gRAVFeUSAWvNxHM1yV/NDPJudOX+W33tHbd2J6G7w1CV7fGH+vp6zGazX5lUVYN43bp1dO7c2eO91KZUWrWuqnUuLCykrq7Oby56a6FDk1Wv15OUlERhYaFfbVl1GqwSWq+H22+389RTBrZuFejVq9HSzZgh8dxzeqxWgbg4mTfftHiEWprCgQMCTzyh5+OPJVJSZB59NJf7708gMrL15CtbciqsCp37s35ms5nVq1djNBo9pqEq4cLCwoiJiTmsGtZgkJmZydatW/3er79qHFDIqiZZaA1wkiRRXl6O2Ww+6sLzpiIdTaFDkxWUdevy5csZPXq05nZBEIiIiHCtWwHGjbMzY4aeWbMkZs608d13Iv/+t7Iu1etlHnnExmOPBV9fWlcHr76q55VXpEMCbTYeeMBGRYWTgwfLiIxsuXWcFg53KuwtpaL1230a6m0F1QoeSZJYt24dQ4YMabU2G6A4EgP1VwXtahwIXOsKjckVUVFRR71mdY90ZGRksHjxYj766KOgj+/wZB05ciSffvqpX7KC57oVIDERrrnGwcKFEv/7n8j27Qor//lPZV0a7ADqdMLixTqefFJPQYHI6NHKurRTJ8Va63SJbN++vUWdLlqIioqirq7OlSShRUB3JQf39EmVhFFRUcTHx2vmLweCyWTSbLfY0lBlRb1726hwr8ZxFw1Xkx7S09P99nXNyMjgjz/+CNoR5A/ukQ6Hw8HNN9/sVw9Z8/ijOns7wJAhQ5g4cWLAdauaC+pOmiFDbMyfL7F9u8hJJyl9Xvw0LNPEX3+JTJigZ/VqHaec4mDBAjOnneY5ZY6IiMBisTTr+tEf+bzrWK1WK+vWrfNY+7XUNNQdaWlp7N+/v9XJmpiY2GS+sHs1jkoStepG7evqLewNynJr//79rF+/ngsvvPCorjPYSIfm9R/VmY8CBw8e5Oqrr2bv3r107tyZTz75BJPJ5LOfTqejb9++gCIN+dVXX3lsj4qKQpKkgE2TvNetADfdJHP//U4MBoHlyy3o9ZqH+mD/foHHH9fzn/9IpKU5eecdC9dc4/A7ZTaZTFRWVvo8AO7wNw317jinfh7uU1CDwUBcXJyHkoMoilRXV7Nv374j0qc9GsTExFBXV4fNZkMf7IfaDHAPH/kT0/OuxjEajR5VN961ru5ITU1l2bJlPPLIIy16H4HQZmSdMWMG55xzDhMnTmTGjBnMmDGD5557zme/8PBw1q1bF/C9hgwZwsqVKznvvPM0t6vrVm/3+8KFNq680shXX+m44orAiQS1tfDSS3pmzpQQBHj0URv332/Dnzdf9YZGRkaSl5enSUZvQTF3EkZGRhIfH+9W1XN4ljk6Opra2tpW9wqrpCkuLvZLmpZCRkaG36msCneFzK5du3rUs6rPhrr2dTqdrvfR6XS8+eabrXMjftBmZF2yZAm//fYbAGPGjOHMM8/UJGswUIvR/ZEVGtet7mS98EIHXbs6eestyS9ZnU748EMdTz5poLhY4KqrbEyaVENSkpnaWgsHD/qS0FtOpbKyEpPJ5CEopjpkWsoRo3qFy8vLW1V6BRQrtHHjxlYnq9pkO5Dkj3cYx7vW1X3t607k+vp6zZlfa6LNyOqu06qOxFowm80MGjQISZKYOHEil112mc8+w4cPZ+bMmQHPpwbP3R8gUYQ77rDz8MMGli+30adPgwfx/vrLwIsvprN9ewS9elXx2GM76Nu3jvp6A0VFjVbQXVJTS05l3bp1JCYmtnrVRnJyMvv37291soaFhTVrPu3hQBX6DkRWvV6PyWSipKTEw3qCb62rSuS2bkoFLUzWQOVy7nBviuyNffv2kZGRwe7duzn77LPp27evhyYsKA9lVVWVhxyHt6qf2Wzm4MGDbN++3RWst9ls9OypIyJiGC+8YGH69AMYDAZKS6N46aUufPddJOnpDubOreeaa/QIwpF5A1UL19oPbkxMTJtMhUFJIS0sLPQrG9pSiI2NJTc3N6A2lE6nIz093ZVX7v7sude6Jicnuz43s9ncpu0eoYXJGqhcTo1rqUXE/kZ/NWjctWtXzjzzTNauXetB1i+++MKlwn7BBRdwxRVXuLSd1Gmod2aMOhVVU9PGjpV5++0kJCma99+XeOMNCUmCxx+3cs89do62MiohIYGtW7e2uvyJIAiucr3Wtq5JSUns2bOHbt26tWrM1Z1sqjSLN9Qwjl6v10ziV2tdTSaTi6yqvnBbos1K5EaNGsX7778PwPvvv69ZmF5RUeHqNVJWVsaff/7pE+uKiori1FNPZfTo0ZxzzjncdtttDB06lKFDh3LKKafQt29fevToQefOnUlLS0OWZZ9Sun//247dDsOGhfHyy3pGj3awfr2ZiROPnqigOMlsNlurV8OAMiiWlJS0+nl1Op0rFNLaUGtYvRtZuUO1rlqVQmqt68GDBxFFEVmWA77X4eLTTz+ld+/eiKLo6t8TDNqMrBMnTuTHH3/khBNO4KeffnIpGK5atYpbb70VgK1btzJo0CD69+/PWWedxcSJE33Iet5553HxxRdz4403snHjRr+aTOBZ3+qObt1kUlJkDAb4/Xcz77xjJT29eRPu4+Pj2+TBjYmJoaamxqNcr7WgToVbGzqdjvj4eI+aYq19VOVDrRaWGRkZlJSU+MjBNAdU4cCRI0ce1nFt5mBKSEjg559/9vn/oEGDmDt3LgCnn346GzduDOr9unfvzp49ewLKi6jypFpYvtxMcjIt0swJlPstLi5udblQ96lwa587NjaWbdu2tbrUDCiOpu3bt7vuWU0UcXcgNjQ0oNfrKSgo8Jkyq7WuWg2ajxbBCAdqocNnMKkQRZEePXqQm5vrN+VMEATCwsJoaGjwcRa0tBJKXFwc27dvD5hp1VJITk7mwIEDbTJQqOoVh5OwfjjwLixwf11dXc2KFStcDky1sED1YURGRnLSSSe5VDC9B3mTSeku3xYOOi0cM2SFxmJ0f2SFxqlwa3v2RFH0KShoLcTGxrJ169ZWFzUDZf24adOmwyKrd4OspsrrvJtjqaG0uLg4amtr6dmzp99z2Ww2VxjHW7sqLCwMp9PpE58PBi0hHHhMkXXEiBHMnDmTcePG+d3HZDJx4MABvyV1LQm1dK21ySoIQpuImoHywAuCQH19PWFhYQHrWy0Wi0+DLK28ZjWjq6kZSlRUFH///XdAy+gexvEmq8PhID4+nvz8/MOuuDla4UAtHFNkPfnkk9m4cWPAqWZ0dDTV1dWtfGUKEhMT2bRpE507d271c6ekpJCfn9/sZPWWWtEiY319PTk5Oa6qF9UCqqG0xMREFymb0/KrneNUDWB/+6hhHPdqHFDIGhcXR1FRUZsnRMAxRla9Xu9y5ARqXGQ0GjVlPFoaRqMRh8PRJg6Xw50Ku3cr92cN/UmtqCV27gTMyclh8ODBrb5eVzWAA82kVOvqXo0DjU2rvv/+e011xCPFF198wd13301paSn//Oc/KSkp+V6W5Sb7chxTZAUl9XD58uX83//9n9991HVra6vXA67paGtrJKlTYXUa7m8a6i445m0F1eKCw61xBSWEpOZItyZUhQq1o5wW1DCOKn/r3ZTqjDPOYPHixc225tcQDgyqgc4xR9YRI0bwxRdfNEnWgoKCNiFrYmIi+fn5zUpWLb0j79eqRS8qKiIuLs5jPagWGajkbAnrl56eTkFBQZskw6v6wv7ICp7VOGoYRyWrJEmcdNJJAfvCtgbalKxNtXa0WCzcdNNNrF69moSEBD7++OMm13tDhw5l0qRJAdetMTExbNu2rblu47CgTkeDCeG4d6vz5xX155AJDw93eUXdS+z++usv+vbt2+pe4bi4OLZt29YmYRCTycSOHTsC5gurjZrXr1/vCuO4t84YOXJkmxIV2pCswbR2fPfddzGZTOzcuZPFixfzyCOP8PHHHwd8X0VyVKSmpsbDWeAOtfKiLdatoChIlJSUoNfrNaeh3nIr7utB9xK7I3HImEwmDh482OpKDqqWUSBnT0ueW7Xs/gZ79bN2D+O0p9YZ0IZkDUbweMmSJTz11FMAjB49mvHjxzdpkQRBYMiQIeTk5HDOOef43a+5162q2kMgr6jqkHE4HOzevdu1/nOX3mzpOle1gKK1yQrKVHjLli1tEjZLS0tj1apVZGdn+/1sVVG13NxcH7K2dVMqaEOyBiN47L6PJEnExsYGJcalOpkCkTU+Pj6odauW5tHROmSsVivr168PGKxvKajT0bZIkAgPD0eWZc0MspZGMM+PGsZRZYLU4vO6urpWdwhq4ZhzMIHiZJo1a5bf7bIsExYWRkVFBWVlZZok9Kd55C69eaQOGXXd1No6RaAMKuqsIpAuVEtBTe73VsZvDWRmZrJz586Ag727qJoaujnuLWswgsfqPpmZmdjtdqqqqpp8wPbt20dBQQHFxcXMmDGDM888k5iYGB+HjMFgwG63U1paSmRkJOHh4R7CYy0dB1WT69vCaZGcnExRUVGbkDU5OZlVq1bRpUuXVo+5RkdH43A4Alp2nU5HXFwce/fudalINDdZH374Yb7++msMBgPdunXjyy+/jJNlubKp49qsRK6p1o7gWfP6n//8h7PPPrvJL3jy5MksXLiQyMhIqquriYuLIyMjg169ejF48GBXrevAgQPJzMzEZDLRqVMnUlNTMZlMREZGtkrCgndbyNaEyWSiqqqqTcrmJEkiKiqKqqqqVj83NN0mEhrDOKqjrzmaUrnjvPPOY9OmTWzYsIEeiv7to8Ec12aW1Z/g8RNPPMGgQYMYNWoUt9xyCzfeeCPdu3cnPj6exYsXN/m+H3zwAQDvvPMOZWVlAZP64+PjKSoqahPrppbrtUUVjtrlri2nwgUFBW3SOyY5OdmlYOFvza6GcdSO7s1tWc8//3zX66FDhwIEpSzXpmtWLcHjKVOmuF6HhYXx6aefHtF7jxgxgkceeYR7773X7z4xMTEB+7u2JARBICYmJqDecUsiJSWlzabCJpOJ7du3t0nMVSWiu2Cf1j46nQ6dTkdpaWmLCr+99957AN8Fs+8x6WAC6NGjh2tk9DeCqtkpgYLlLQl1KtwWZFUJ0xZeYUEQSEpKoqSkpNWzyOx2O/Hx8a5B2tu56B7jVrvLHYmyYTAlctOnT1eXXB8G857HLFlFUaR79+7s3LlTXRdoQvWMtoVrPj4+nn379vmoNbYG3Bs+x8fHt/r509PT2bZtW7OQ1V3JUivTS8u773A4qK6uJjo62lV04B3jtlgsbN68+YiaUjVVIjd//ny++eYbfv75ZyIiIoLSEDpmyQqNxehNkbWkpKRNyKq2utDqyN0aUMvH2oKsEREROByOgFlk7jnP/kioOsnclSzVH/f6V2+nYXl5eZPZVJIkuXwLzVkit2zZMp5//nn++9//HtYgcEyTdcSIEbz11luMGTPG7z6xsbHs2LGjFa/KE2oIpy2yetpKakbN9lJTSePi4vwWoXsnmkRERGAymY64rYiK+Ph4duzYETDWrdPpyMzMJCcnh6qqqmbzCI8fPx6LxeLqILF+/frZsiz/u6njjmmyDhw4kA0bNgR8GHU6HaIotum6dc+ePW1CVrVRU0VFRbNYV63CA3cSqnFuNdtLr9dTXl5OTEyMqwg9mE7nzQFBEEhLSyM/P18zX1iWZZekiyRJzJ8/n0cfDSrC0iR27tzp/a8miQrHOFnVxOxAxejQ2OmttYWwobGXals4eqBRVzgQWVVRMn9EdG+w5U8Pyb3DnTs2btxIbGxsqznZ1O7tpaWlFBQUcM8993D77bdTUlJCUVERhYWFlJSUYLfbXcqQ6enpAaMKrYVjmqzQuG71Kvb1gMlkoqysrE3IKggCsbGxVFVVtWqtp/rQGo1GysrKKCgo8FkfaomSedfAqiQ80mm0GnM9WrKqTqaSkhIKCgooKipyka+4uJiioiKXUqEavklNTcVoNJKbm8tZZ51FWloa6enprv+3dvy7KbQbsjZV2zp//nwefvhhV0ri+PHjXWLggTBixAi+/vrrgGSNi4vTmpq0GtQQTnOQVcszqlWIDo2eUVEUXZ3X/HlGWwrx8fHk5ub6jbmqer/FxcUUFha6CFhUVOQiYVlZGU6nE0mSSEpKchEuLS2NoUOHkpaWRlpaGqmpqT4Dy9q1a/n999+54oorWvQ+mwNCE20BmleW3g8cDgc9evTwqG1dtGiRR7nc/PnzWbVqFW+88cZhvXdVVRXnn38+v/zyS8AHLycnhwEDBrR6Yj0osb/Vq1czZMgQv/toeUb9aSJpeUbd/9byjJaWlgbM9mpuyLKM2WymqKiImTNnEhERQUZGhssiFhcXu9Ix9Xo9KSkpLtKpFlD9URtItTdLeBgI6sLbhWUNprb1SBETE4Msy03q9aoxx9aW6nQ6na7pZn5+PrIs+5DR3TOq1XD5SDSR3GEymcjNzW0Wr7Asy9TX17usn/tUVP1dUVHhqnxKSUkhKiqK3377jQcffJBevXq5yJiUlIQoih2ZhM2KdkHWYGpbAT777DN+//13evTowSuvvOJxjD8IgsCpp55KTk4OZ511lt/9TCZTs+rqqp7RQN5R9frUMruSkhISExNbVJ5TC6IoNiloJssytbW1PlNR9x81OT88PNw17UxPTyctLY1+/fq5LGFCQoJPm8+xY8cyduzYNvHItyYef/xxlixZ4mrqPH/+/KAjAe2CrMHgkksu4dprr8VoNDJnzhzGjBnDL7/8EtSxajF6ILLGxcWxe/fuJt8rWM+oliSLu2fUuw62rq6OnTt3BjUANRfU7mjV1dVs2rSJ5557jgsvvNC1FlQtodofKDIy0mMqmpaWximnnOKyhCaTKWCv3UCYP39+M99d+8TDDz/M1KlTAZg5cyZTpkxh9uzZQR3bLsgaTG2re8L5rbfeyoQJE4J+/5EjR/L222/73a5OM51OJ6WlpR5NmJvyjLoXox+NZzQiIoL6+vpmCeGoJKysrNS0hOqUtL6+HlCWCikpKfz++++ccsopZGVlMWTIEJcljImJOWIShuAJd12wurq6w/pM2wVZ3WtbMzIyWLx4MR999JHHPmrjZYCvvvoqqE5csixTUVHhSi17/fXXGTJkCLGxsZqeUYfDQX5+vkuqs7V0kaBRwSFQrq5KwvLycr9T0aKiIsxmM6DMFtytYNeuXRkxYoTLEkZHR3vc0wMPPMBVV13VJrnKxxMmT57MggULiI2N5ddffw36uHbhDQb49ttvue+++1y1rZMnT/aobX300Uf56quvkCSJ+Ph4Zs2a1aT30ul0cuGFF5KSksLGjRsZOnQoY8aMoVOnTpqKEKWlpVRUVATMJW4JqNkyy5YtY+nSpYwaNcojPqi+tlqtLlKr4Ql1TZiWlkZGRgZpaWlERkaGrGAbItimVM8++yxms5mnn346qC+r3ZC1pTFnzhwqKyu5++67/e5jt9tZs2YNp556arOc0z1bprCwUNMSlpSUYLPZEEWRhIQE1q9fz7hx43xCFGlpaYSHh4dIeAxh//79XHTRRWzatKnjhG5aAyNGjGDy5MkByapa2aZ60agxz5KSEs1AfWFhoUdfTzVbRiXfwIED/WbLjB49mgcffLBNsqlCaHns2LGDE044AVCkdg8ntn3cWFan00m/fv3466+//DpwZFnmtddeIykpiW7duvnECYuLiyktLXVlyyQnJ/sN1KekpByRs6ktZF5CaD1cccUVbN++HVEUyc7OZvbs2WRkZISmwd644IILuO6665AkyaWAqBKyvLzc5Ym1WCwMHz7ctSZUp6EZGRnHQrZMCO0PIbJ644EHHmDfvn3069fPxxKq2TJVVVVcdtll/Pbbb219uSEcPwiR9UjRVsoNIRy3CJE1hBA6CIIia5uJfIcQQgiHhxBZQzjucPPNN5OcnEyfPn00t8uyzD333EP37t3p168fa9asaeUr1EaIrCEcdxg7dizLli3zu/27775jx44d7Nixg7fffps77rijFa/OP9otWTvq6BdC+8fIkSMDak4tWbKEm266CUEQGDp0qKsgoq3RbsnaUUe/EILHsmXL6NmzJ927d2fGjBk+2+fPn09SUhIDBgxgwIABzJ07t1WuS6u+Oj8/v1XOHQjtNt1w5MiR7N271+92f6Nfa7djCOHI4HA4uOuuuzykfEaNGuWjDnL11VcftpTPsYp2a1mbQnsd/UIIDu5SPgaDwSXl0x4QTH11W6DDkjWEjo1gB9vPPvuMfv36MXr0aA8CtSRGjRrFggULkGWZv/76i9jY2HYxY+uwZG2vo19HQ1PrRovFwtVXX0337t0ZMmRIwKVJc+OSSy5h7969bNiwgfPOOy9gG5TDwbXXXstpp53G9u3byczM5N1332X27NkueZWLLrqIrl270r17d2677TbeeuutZjnvUUNVH/Dz06bYs2eP3Lt3b81t33zzjXzhhRfKTqdTXrFihTx48OBWvrqOD7vdLnft2lXetWuXbLFY5H79+smbN2/22OfNN9+Ub7/9dlmWZXnRokXyVVdd1SznXr58uXz++ee7/n7mmWfkZ555JuC1xsTENMu52yGa4qEiPdTEDm2Ga665Rk5NTZUlSZIzMjLkuXPnyrNmzZJnzZoly7IsO51O+c4775S7du0q9+nTR87JyWnLy+2QCIYw559/vrx8+XJZlmXZZrPJCQkJstPpPOpz22w2uUuXLvLu3btdA8WmTZs89ikoKHC9/vzzz+UhQ4Yc9XnbKYIia7v1Bi9atCjgdkEQePPNN1vpaloHBw8e5Oqrr2bv3r107tyZTz75RFMaVKfT0bdvXwA6derEV199dUTnC0YC1n0fSZKIjY2lvLycxMTEIzqnCkmSeOONN7jgggtcUj69e/f2kPKZOXOmh5TP8aKA6A/tlqzHI2bMmME555zDxIkTmTFjBjNmzOC5557z2S88PJx169a1/gU2My666CIuuugij/9NmTLF9frZZ5/l2Wefbe3LarfosA6mYxFLlixxOVHGjBnDl19+2aLnC8ZJ576P3W6nqqrKQxY2hNZDiKxu+PTTT+nduzeiKLJq1Sq/+zXlQT1SFBcXu0IEqampFBcXa+5nNpsZNGgQQ4cOPSpCu0vAWq1WFi9ezKhRozz2GTVqFO+//z4A//nPfzj77LNDKhlthSYWtccVtmzZIm/btk0+44wz/DqsgvGgBsI555wj9+7d2+fnyy+/lGNjYz32jYuL03yPAwcOyLIsy7t27ZKzs7PlnTt3Bn1+byxdulQ+4YQT5K5du8rTpk2TZVmWH3/8cXnJkiWyLMtyQ0ODPHr0aLlbt27y4MGD5V27dh3xuULwi47tDW5LBCLr4YYcDgc9evRweUALCgrkHj16NHnMmDFj5E8//bRZzh9CmyEosoamwYeJlkxzdJ9yvv/++x6C0CoqKiqwWCwAlJWV8eeffzZLt70Q2j+OO29wsGrpbYGJEydy1VVX8e6775Kdnc0nn3wCwKpVq5g9ezZz585l69at3H777YiiiNPpZOLEiSGyHic47sj6008/HdXxLZnmmJCQwM8//+zz/0GDBrnKw04//XQ2btzYLOcLoWOhKcG04xKCIPwGPCTLso9LWBAECcgFzgHygRzgOlmWN7fqRYZw3CG0ZnWDIAiXC4JwADgNWCoIwveH/p8uCMK3ALIs24HxwPfAVuCTEFFDaA2ELGsIIXQQhCxrCCF0EITIGkIIHQQhsoYQQgdBiKwhhNBBECJrCCF0EITIGkIIHQQhsoYQQgdBiKwhhNBB8P9bHqFtWIyJ9AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot3d(gx):\n",
    "    fig=plt.figure()\n",
    "    ax=fig.add_subplot(projection='3d')\n",
    "    drawn=set()\n",
    "    for c in cells:\n",
    "        for i,j in ((0,1),(1,2),(2,3),(3,0)):\n",
    "            if (c[i],c[j]) in drawn:\n",
    "                continue\n",
    "            xs=[points[c[i]][0],points[c[j]][0]]\n",
    "            ys=[points[c[i]][1],points[c[j]][1]]\n",
    "            zs=[gx[c[i]],gx[c[j]]]\n",
    "            ax.plot(xs,ys,zs,'b')\n",
    "            drawn.add((c[i],c[j]))\n",
    "    ax.azim=-10\n",
    "    plt.show()\n",
    "    \n",
    "plot3d(gx)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}