005_1D_FEM_Code_Complete.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as nla\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import splrep, splev\n",
    "from scipy.integrate import quad\n",
    "from sympy import init_printing, Symbol, Matrix, simplify, lambdify, Array, diff, symbols, integrate, exp, sqrt\n",
    "from numpy.polynomial.legendre import leggauss  #Gauss-Legendre Quadrature for 1D\n",
    "from jupyterthemes import jtplot\n",
    "jtplot.style()\n",
    "init_printing()\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rcParams['font.weight'] = 700\n",
    "plt.rcParams['xtick.top']='True'\n",
    "plt.rcParams['xtick.direction']='in'\n",
    "plt.rcParams['ytick.right']='True'\n",
    "plt.rcParams['ytick.direction']='in'\n",
    "plt.rcParams['ytick.labelsize']=22\n",
    "plt.rcParams['xtick.labelsize']=22\n",
    "plt.rcParams['xtick.minor.visible']=True\n",
    "plt.rcParams['ytick.minor.visible']=True\n",
    "plt.rcParams['xtick.major.size']=6\n",
    "plt.rcParams['xtick.minor.size']=3\n",
    "plt.rcParams['ytick.major.size']=6\n",
    "plt.rcParams['ytick.minor.size']=3\n",
    "plt.rcParams['lines.markersize']=np.sqrt(36)\n",
    "plt.rc('text',usetex=True)\n",
    "plt.rcParams['text.latex.preamble']=[r\"\\usepackage{amsmath,amsfonts,amssymb,bm}\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computing the exact solution\n",
    "x, r, xi, Ao, l, rhO, g, E, m, L = symbols(\"x, r, xi, Ao, l, rhO, g, E, m, L\", real=True)\n",
    "area = lambda x: Ao*exp(-x/l)\n",
    "bf = lambda x: -rhO * area(x) * g\n",
    "innerIntegral = lambda r: simplify(integrate(bf(xi), (xi, L, r)))\n",
    "usol = simplify(-integrate(1./E/area(r)*innerIntegral(r), (r,0,x)) -m*g/E*integrate(1/area(r), (r,0,x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class geometry():\n",
    "    def __init__(self, ne, deg):\n",
    "        self.L = 324.\n",
    "        self.g = 9.81\n",
    "        self.E = 5.e8\n",
    "        self.m = 2.e4\n",
    "        self.l = 90.\n",
    "        self.Ao = 120.**2\n",
    "        self.rho = 1500.\n",
    "        self.nNNodes=int(ne*deg+1)   #Isoparametric deg-th order lagrange elements\n",
    "        self.ndim=1\n",
    "        self.nLNodes=ne+1       #Linear elements\n",
    "        self.nQNodes=2*ne+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsoAAAAVCAYAAABfYciTAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAQzUlEQVR4Ae2di5EdNRaG1y4HwNoRsGRgmwiwM4AlAkwGdhEBBRkAEfDIAIjAizOADGwmA/b/eiRZrVa3fnXfeUtVmlZL562jo9OvO/f++eeff40yLDAsMCwwLHB6C9y7d++bQPWtjo8Ub1+1uLRwNP6BaPxX9UvRe1LSOzoOvaM0juLfFRkcPYEZZVhgWODqLPDg6lgPzsMCwwLDArfXAkoWSZLfKpn9Fi11/kz1V50/X9O6haPxx8J9FvBJmGfl6DjEjtI4in9XZHD0BGaUYYFhgau1wEiUr9b+g/uwwLDA7bXAC6n2SVRPCfJvSiJ/Vf2P2n/F/uK4iSO8N4J/IxqfFnjT6dFxiBylcRT/rsjg6AnMKMMCwwJXa4GRKF+t/Qf3YYFhgVtoAZJhqcUd37NCPc65K7xIlPfgFLTH6S20gPziZVDrIx15SvGlKr71jS5KFn6k/lGGBYYFTmiB+yekNUh1WIBNUfWFKgFvlEuygOxNkjIrzIEqic21KTdFzlMb7BbpHdf1u8JGnD8s+uLpHpyIO4630AJaDyTJvygh5vUd/OY7tXnP/anqIpapb5QVC9yE2HITZFwx7+Hu66z7ye4oS8n8o5XpyrfnatfFd+AEQ9ITP5ohoLA5vZI8PLZclIxmHPtxDTYCCIcgxcc0XN2non42u69Up493dOT8Z8H9pmNewP+OKpy8P7bPhPPveHIdj5nd0PVK5zzaJ5MpdpVz+XuYo+gLMTlZfBQVCQieufpZ9Ynm5Cz2l0cHLpOvZTNLzqCL42+luPFd1JoPW+tnD2/DRq7eloxRafFlrcXyUI0vtuYSwCDrwj6RyM4jvKPPuST24Li0B9z1tgD7QLxrjN/8GsTFf2MMu1QNtC6uzV4fFc9kil1l3Kffii2RQFj/zbgPvAObybgV+20ZRQ9/6I79QdZFXFO/HVN7eQeeLVtaunfK2R33s/l8byMtNt5JO1RF+A/VTyMdtZnAP1V5F69JW3AWvgMHT1US08RX5yxsft7jWdEPLHKmfrWZzD9zuFo74M34AKfCFX/iHfqgmewT+pAJvRljMvOKTDP4kuZVn0s+a87W5HTxO+CsuRQ9bEtlsjgyDx+UctKnGucGXYHfDQd9FdtmgnXltPyt1C/IA4+ZD+u8Z/1YvEXTsmUm0+b8dMoYY9HLaAPhc+HzN3RiX+2o8YV9anC1PmirMukzHqGvurZ7cATLO8qrceroODodpXEU/67I4OgZYJiUx7Svqoq/HcNqMrr4HXBW3A/2a8ZU8e2JVT2wlt3EvyljtKtgrfgb4eMx8Ngd94Mtm7x7bBloNnUXTWt/Cryh1x33M1mSjQ4vOAnCxyeLgK0+EhC+8N7k4eJ3wJHc1BIaNsa/c3l0jvPOJlznXLFvyq1xHocxCcmQwbjY4kXOI/SzoP7I+3U+4xvH1I8jVMcizKmO4sPX8wt5W/TBUb02c468KtZcCm42Zy1dA23mGyYLv8rxNb4Kp7Eumwm+KWeguZg/9S/8LZcz06nmw9b62ctbeKs2CnI5elsyRnriOVv3oZ91vrrOgpwL+4DrVviqzhIbneNHs+Q5p+fiCG4kyufr/sbbIfhjSw8u7gR67n9bPhRhTn0Uz64YVvJ38V046KtYcT/ANmNLLrNob8YqF7ZTH0vGQLM79gedFnFN/T0x1c5zoo0C3809FBki/NrRlTPAdcd9+Kow7zMb3af3YPlM+LVHQK/Vz88hsWlvFRffhSP5+7vCl1cf0ruoGicwEXziKxpqTpHouYy19fNN4JyFCkpeeP0AOZ0SH6GVsHygMXudowQ44Tlz05qfGjt3Lmq49Ln4FtzeuVwT7oL6LV06eff4WyIte235sLV+RGwX7yTEsYYrI1xY5/HRdc6VmEXAX5SGfRbwGx3fawxZpyK6tN9ofU/y6Jx4xA2FvGzi5ICjfbstIN/gOxYutiifq+b77GXtERPz8OdoDHPxLTjZZtcenit0SW1Ln05ZuuOv7HWKuI+Y3bw7ddsCd2N/d9yH6ZqNTpEoI/i7imYkk5S0UZyfLv66+C4cCfFf2owi/5JhTAx5t+dsA67Ei+efC4fNrFZIfqffSpXBIx/g2Azzd2W4NfALA3kRDnBf533XtO3OxZr4Lr4Lt3cu1+S7iH5Xlx7etr8VRLd82F0/e3kXouw6tWTM1mAtPr2FcwiMpRBb9ilhV8+1xrkIf8S6VuUuBRvmJxkCT4/4oJfjVFo4wEJPwPg8bb5xgPZUjo5D5CiNo/h3RQZDT/yWnxTkgo71Rvul6lXtE0djmIvvwt2EuM80u/oA65Y98XcrrlkxNQi3h7erVwuuKafWR8y9euM+vKs2etCSams8E2gL7OHaoIvvwsFHGw2bUa1wNcX4mzDIOb9HypGrdTZOrpS4/c9kLIpg2ZBmCW8OBJ5gSIC5muGuNhslNHmVY5EYqz8VwbJZ8pgWnFnRGBMf7359rPYXqsAjN+V1i/452PG/QZYWoUudcwnTNZfSAVvGxcT88Ag++kVLt+7xvTZrySmZu/1NNFs+bK2fPbxdwxl6uzKeiRZsa/74KMgzG2vZJ9dBsPgQm/WUdOuIL72Sbc50nIrai/WcjeFziw92Gzjcjd6ieWgc2cT/EI2j+HdFhpaesiN+lPt6dV+CzkWX4OstNrO1lAO7+C5coN0V98ER/Wsf+x0Z5RtdsV80TxL3sWEvb3Dc0tJdvPP1kJPFF1J+Jzqc1vyxGvcB3rLRoUQ5EyRtDDAMJWbzbCZrJSrSwnfhqnxkAIxIYjltMDrPZXoq46eNR2MkuHxRPEts1Qc+d6Cnx6ZVRupkIgXLFT+OyZHNkGOrALMGl17HEG3o/qBKcvxK5yTu4M3k1flFlUNzIaFcfAtO+vfOJfA/yXaTzwmfef1TR165uaiNyNKlmDBLzh5/C7o2fbiQI951TesnjvfwjjjG0dK7pCPdHqtvIaP6WBfc0SkL8BRwptJjH8FCk7X3mewwXWQFGUicUzyZCI8/wwI33wJ7YliutYtvwWmtESdisfZwAe+KLZHJzqOlT0bbltGNvz1xLZNjNe4D4/LO6RltW/ecVoi7Zey34z60Wja6nzO8oHbM4PeSd/G34HhR/XtN7vSvZDNBuINbvkbxk8Z/KBYiKPxUSAmbkTpvCo+7yRTuLpF4PVYlEatt1hpKk8SX8ItETXhcAZcJNDxiYsxCLMfVdaVlay4cwVz8HM6aS9mYhPgsCqE2Fz7YnaTnKkuuC4HIkrPT3ywfrhihun46eVfILrtcvZeY08cotTXO0xcCYVp/arMmow/EC3rALPsIn6DM40cuYKckGWQVeNA/yrDAXbTALIbtMICLn8NZcR9ZDsSWHap0oSR9emTsiL9WXKtIXI37wHXwrpCtd/XoXlCoydkT9yG3aaMHUpgs/ndVjm6Jd1HyTabEjVdS8bFkOc65i+/CLXhIPxKg3zQJ6eMHteMj2drd4T8ET3L6FDxVnILzZiIV4Pit3cjrufpIapnIn9X+EN5qlwX4fMPNx/8nnFxOPjRMHwRpbO1RRKIhvsieEoU0EO7uajzKmw1NPGq0d89FIO7iW3DYU/JDOrdR1GMxl3GgOILLu+X8GkGNTgHefWrpYlCdySl58UvL3wJs04dLGYS3WD/A9PAuae44n+ld4q/JCFzwjw/V5D1hEmQ2pdehsjan+Q76uPaJcHyMxxOeuNHxitXiYlfjowwLXLoF5JtH9vZS3qMxzMW34E4U99FxM7aURthxbunToLuQMcSrZuzvjGtJDOFV4z4ALu9E7FhjoXtObk1ON+5Dy7HRAwgK9knO3G1nzlpLsmPftBHVaLr4LlzJIxiAzbKWCCIXuq+V6ZGsaHDkZ8FW9cgIcGeXTTkV4f0iGtxd5udGSFbjneAEo0basPNO2sJ/U/RBo+vR7or+OAh8SQ7LO+0Fy/engo2JaZzf94PvL7ZWbeXiu3CBuTuX3O17KNpr/l7TKddvV7tTF+bFldPyt04fTjoKj0R8bf1YvBMxo9Ghd6LWkHGCw/5qzGKA8GKy+5faPWscmtM6Fl173UyCjD/DApdogeD3a7GuS5LeGFYSd/FduEDfivvA7oktpQ57znv06ZSxGX9Fj9zBzV2SesLbivvANXkLppbnJB5lo1P3Cb0lJ7YX4Grch4hoWLH/wcTx2B/uoExJZUEm3lFu3WFx8V24SQwZgCTwIxkrGSoYhY2fBUatyT3hh3HawHwsXO4K54W7U3x1Tj+0vlbFKZmcWYGf4HCcaJM0rn6SM3iUCXGCiQ3BskFTkhMGfHRa8D0HvZC/XXNRkcDFd+HcuXwqWWpX+NO8yIbNOajo4na5ukCvKWeYd9ffLB+W/ukCTPRX14/kw4Yub/RxS1PvnNCWjKy5HLbShhdPmrjwo91c47l9BN+iX2E5uoYFbrQFemJYTVEX34Vz4z6ydMWWmvAH+lx9LBk7Yv9J4z4xtYN3r7ks3SPRA7EfPlPcD7QsG5FgHapixtVH7YeduWMz+ycbNV4uvgsHDxWSWN4fnOmmPh6TssHrMN0VUnMBg9wLfSpw/L7l7AeywVOt/jMB9XOXcDGmPpJfBFr88wP1kUSDN/3nQB1Jymcy63yhZynr2rlwP1VN/7lmDa7sF861mnPJM9mwIudsLgVXtZX6mbfVfzKjMfyGOZp8p+QTz7fgNGbbTLCWnILr9rcoK8eAP/Ph0O+sn128xXPVlhqz9HZlDHD4OLKmuVOb4IgBeL9xFiPyc42DV7MPT4fWZK3+172c7miv23zY5vraRj5vx7DaPLr4HXBW3EcWlbX1yhqvxn71r8aqUr8tWI1ZdhOcLaNgd8XfYIu1uNaM+xn+IpcJY2t5zqYtO3Vvyil6u+N+pmOK/aubBMBulVBsHGmDUJsEj8lIG1How2MXybP6mvhB+CacaLEJAkeSlFeSzNl/k9M5k5qcU+0od9JlzQaCRZfZAtM5C5f3YtOmHORmoVQTUvUzodBKckSe2RiOgWxl4ge/pqyRXnkM9KtylbDluXCduUDmC59zZFNpzqVg8I3ZBYnOWcD46mzOcn0DDExWYYIMrWDQtFmgY8kpebr9rdALnUofttbPXt7CW7WRxly9LRmzOcHuae7UxlcW/9Uqt03AXdgn9LOeZ3EM+qqsz+oGUtIe5/7+Mmx1fWwl/27GsLAWWDuzNRLWThO/E64Z9wM9K7bkvib5V2NVDhfob8KadrNlFL3dsV+4i7imvp6Y2s1b9Fv2sXR35Qz8dsX9MJ8zG90LnaK7v4Tb8V+JQvxw72O1vxbt2aNswSH47MM6uHbgsxlt8hEtElUSy1rhI7jZO1uC530b6FIeqi7knkbCH8GzGT5VjTx4rEKyMb2vqHH6kTF/xM+VCXCLIngcBJlrP0mHXMiH3bhw+VbwJNZ80AfOO/Wl1zB03lUCLTb27ncthducC4QR3IXPeVRavJpzKRjsHV8xYL6Zp9lv32b0uLhCT+abI4/68OnXuc1E04WzbCb68d2pppzi3eVvgfaqD4uevX56eHfYqDk/PTIGffGLWJiD1fUIkOiv2icSEQyB/5Hq29iX+0TsG8dhgdtkAfm9FcMEd63iPnMgmZqxJcBZ8bwHtsNuloyBd1fslwyrcU1jdtzv4S26PbZs6t4jp2C74n7Qq2qjkyTKMBjl5llAjkTSvStRvnnaDomHBYYFhgWGBYYFhgWGBfoscL8PfEDfMgucSR/qKMMCwwLDAsMCwwLDAsMCwwKFBf4P2GrgHZWpbuYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle - 0.0026487 x + 0.00651105842815294 e^{0.0111111111111111 x} - 0.00651105842815294$"
      ],
      "text/plain": [
       "                                    0.0111111111111111⋅x                      \n",
       "-0.0026487⋅x + 0.00651105842815294⋅ℯ                     - 0.00651105842815294"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geom=geometry(2, 2)\n",
    "usol = usol.subs([(Ao, geom.Ao), (rhO, geom.rho), (L, geom.L), (g, geom.g),\n",
    "                  (m, geom.m), (E, geom.E), (l, geom.l)])\n",
    "usol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALQAAAAPCAYAAABNwOWEAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAG8ElEQVRoBe2a/3EVNxDH/TIu4MWpIKYDMB04HZihAqADMlSQgQ4gFYDpwHQAdgd2B+DXgfP96GnFSqe727Pvn2SyM3qSVvtbqx939zZ3d3cH/8N/OwKbzeax5vnKeyncVv0j4W88/t/e3viElpNvs0PfVT9SebvE4Qi/aI4l91XWQ1Dpo+dLxlWVk2n4j6Itk6NxZLxRwebfVOif9+QtoZWMAuJ7rM4ryTS705jw2P5nJjxR/YO+ty+PpUr0i+M7pntErqHbGN1qgLhY3GgDT2Trbt/c/0rfe9c/UvtFS+PGD7J958INZBld1O9VdMvYA4rgUuXM9XH6WuXYcFN1hF80yHzv5aj/UgUDim7GBSQL+k+NXm0Cd239TFfJyzjoKnlLaRsd2MEiSbHKsrCvxZGwGF9sNh7h7hVf8Q10O5nRGCGDgm3U2Lk1Odkfm+/XhhcNC5nFUOUAvCrEmOTHL+RW8pyMWb/hVcGuB+u2ZCapqkTJTuL4hRk3VosmxI/BKjhfkk1tnAF36eXTV6mSVf0LitGpjd6X1rdaOGS28sK0JodagM0Eu01eJnQwicKRALeNjFB8PM+UbqOTntkYZTmV7cbva8nCn8ruzEvMq3lo+GxOe7EI+b2m7l8kDHimYsdRQuSfr6pPdRSQIFMQ5UfHLpckT8GhX4H0nQnB7mDHeRoX7R8UR8y1CN0RWEKb5MkObKjsdYpO1b7txIar01Z4dk+DaHyM3o7yMd2MR2NUZM40kNe7TzNnJOZ9IOr3arotoZmcHx2LLdkYn4IQv5Lxi8qv1CYsTwxdf3fjTrwTnek38rZm92DBXahs3SAni5fH0BJaE/VcNnywTlPjw82Ejd6eUHwa+VO6IY3GqBE77LrY9XLgOxyiYXEvhVm/19Z96AROGXs0NvgQfvHiMMnHA5dPHIJ3pXHq5yoElR2Wo7MsBtqi+Sw8K5zdkh0dOq4l4AssoYVJsjhK20Xh5Y2dDGnipY+dDTk+sQt/06jiO6c784ZiZHokk13WbCFGXCOSjap3Goe0sgOEgAdtoDe2H+n8Rv1eW/ehM3TXsctWrAWiQ1IcDfPLWSaDZH6qQlC/qSRoAnEih8u1Q2MkLU/dJVnVfiYci4IEpEYe9QCitJLHdYEToncED+QaQnz4BW+xWW1LhFB8IrpF4+djNkayAfpP8ifZkHVcq+YKZxsEMWVOWsAnwF+h9pjp3yV+r6bbrhzTpv1cpXN0Y+O2ytO4gnil8k6FXe6jyqWCyy7rgXenftdm7JPK335CHR+7DpPDBDBZg8lZQNueGBIZAh6sPuBbiPonkY/PEt2hGMkeErcsKLVZqMTKn0AvMMfHTW1iaXy2uUG2Fpjfq+kmoacMtVWW7lEjXjyIX8FldRK0cwWQJ2ULYG935Kl+q3KiQvA5Rpks3v1yn+WB0a4CSR50QJQ20/mJ3guY+RUfPDwjVO+qhQvHJ6pbOsIxmjCb+B5LZ9p5s8zfhePEe63CKccYLwaA3nzsR/q/Yb/X1H2IMBmPSSRKC4YbdWYJv/Sw4nlVeNUo+qY+OyqFBEefTZqaA0iTICxBZxIKSPZn6WG35lWbyWN8llZ82MWiGvUXQS2IL70FEF+bzPgaiq9k4NMS3aEYSS4Pw3wRfNLanfs2x8lW4SofxG+Le1FMon6bTdCvofswC+T4sSQxHdS2Q9s9y4/59iy/AkPg2GHZLXnTMZWwBK9nD+zATZbnd/T9iH4lm3EWRrJ/AS06n4qeq4MHFiK7GXhOgnJHFo6r0iPhSiIIl2zHjixkNj6i49RZons2Rlk3cnu7ZYqNbGw3l8xWKvg5eabmqxA3jYjfDUvVXa5bhiKBHab3Up3VWX2ggL4tUX50qAxe8mc8hpCgqtLOqmatS3jsKXbSVqm+YhmP8OxMZWwJrcmwOvP27CbR+WxfxUQ4HlCTL4wJ7h1f8Y7FjNNH4ge62xgN7Ms2Idd/pDrLurzdLEwc4K5e+Wh9jc19WCnz5XiqvJKM1XQXIyWUI7r9gofTxRm12WVxcJDkwkX4mdjqy576OIPMFk9ClslQG93Y421kUtO92oJFLUBP+YyacWFaLyvzIrRMfsYx2fjM5PjCLt776jobn1bvmG6jk55IjLCz+tKnPklILH3ygsNGj0N+NS+m22qNjyZ0tn/W7yxjFd3lz0k6JkmYNyr2AMgrtb9kVHUkiQ7Fg4efBfwklj24qZmuFiTu4Fojmdx7sQvgiOzZwy6J3f5Yrd5XayyB5IVpYRA9icqxBx+AjST2O42xkAzPmAfe4lR3VtGH4mtCpnQbDbXoIjEiqe2adKQ2seJBeqe6QJZlfeztxhEC0bJwoSE+1FyByJWvxEd1AtExFskr/DC4t+5/ADg5YtSIoR3ZAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle 0.0368842046591091$"
      ],
      "text/plain": [
       "0.0368842046591091"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sqrt(integrate(usol.diff(x)**2, (x, 0, geom.L)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAFmCAYAAAA4UA4EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXQc1YEu8O+WFkutlmXJ7mpbttVqebeWsMS2TFi8sXkPCYtNwiQsw5sT3gDvPfIGyDvAyZickHMCzGSWxGxhsYFsYLaE2MIsgyV5EgKSvFstCWxLLXmR1d3W5rrvj7aEhbaW1F23qvv7nZOD011V+uh7xOeqvnVLXHTJNRJEREQWo6kOQERENBAWFBERWRILioiILIkFRURElsSCIiIiS2JBERGRJbGgiIjIkpIj3TA9PQ3FRUWYMcMLh8MBAGhq8qOqugb19Q0xC0hERIlJDHejrseTh5s33ojiokL4/c1o8vvh89XD7XZB13W4dRfa2gLYvuM9vPnWO2blJiKiODdkQd1+2/dQ4M3Ha9veRGXlfw96EF13Ye2aVSgqnI9Njz6G5paWGEQlIqJEMmRB5ed7UFdXH/HB0tPT4XRmoLmZBUVERGMz5CSJ88tp4YKvD7rdsqVXAADOnDnDciIioqiIeBbfnX9/K+7/p/+DtLS03tdckybh0U2PYOOGG2MSjoiIEtewkyR6OBwO3LzhBixatBAvbXkZkye7sXbNKuwo24mnnn4uxjGJiCjRRFxQPdasXomNG24AADz+xC9QuXvwyRNERESjlTRl+syHI9kwPT0Nd/3gTlx7zVXY/PRzOHDgAO64/ftISUlGdc2e2KYkIqKEE/F3UM889Z+QErj19n9AWdlObHvjbdx97334Wkkxnvj5Y7HMOCaZWTkAhOoYNGKCY2dLHDd7sua4RXwGtX//Abzx5tvo7u7ufS0YDGFH2U6kpKTgwIGDMYo4NuPS0tHZ0a46Bo2Y4NjZEsfNnqw5bhEvdTTUZbw33nw7KmFiwZM3DcFAFgAJn69OdRwiIorQkJf4Ht30CFavugbp6WlDbQYAKCychwf+6T4UFs6LWjgiIkpcw87iu/2276F00UJUVdfA7/fj8GEfmvx+OJ0Z8Obnw+3WUVJchGAwiF899eyIVp4wQ2ZWDtpaTwIY0WRFUk4gMyubY2c7HDd7sua4RTTNPD09HYsXL0JJcSG8+fnQdRdCoRCa/M2oqqrBrvIKyxVTDxaUXVnzF4aGw3GzJ2uOW0TfQZ05cwZlZTtRVrYzxnGIiIjChvwO6u5//MGAfyYiIoq1Ic+gSooLcdutfwefrw6lixag6tyisAMpe+/9qIcjIqLENWRBPfHkv2HNmpVwu3UAQGnpwoE3lCwoIiKKriELqqq6BlXVNQCAB+6/D4/+5GemhCIiIop4qSOWExERmSnigiIiIjJTxEsd2VVySorqCERE8S1pHJDiBFIygJQMiHP/BASQ4gBEMqAlQYhkyJP7IVs+i+iwcV1QKampWPftW/DSM/+Crs4O1XGIiOxLSwXSXRAOF5CWA6RlQ6TlAONygPbjEM5cSKML6Aqe+18IsuMUIDTA6ALOtkMa3ZBdgYh/ZFwXFBERjUJyOuCYDJExBcgI/1OkT4Rs+wJISgXaT0C2n4RsrYNsPwF0nAK6AuEiimaMqB6NiIjsJ3U8usfPhsiZBDE+DzDOAmnZQLARMngM8siHMIKNQHsLIA3TYrGgiIgSTdI4IKsAItMDkT0LIi0bXV2ngVN1kEd3QbY1AO3HVadkQRERJYR0F0TWDIicuUDmNKD7DOSpw5Cfvw8j0IDMNM2ei8VGi8PhwIrlS9DU1AynMwNut44tW18d8THWr1uNQCAAp9MJAHjt9TcRCoViEZmIyL4yciFy5kHkzIFInwSj+VPI1sOQ9X8CgsfO21CEL+lZjKkFde89d2HzU8/C728GAOi6Cw8+8ENsevSxiPbXdRfuuP37fY6xaOECrF+3esRFR0QUlxxuiElFQNokiOxZwOk6yMZKGCf2hScy2IhpBbV82RL4m/y9xQKgT8lUVO4e9hj33n0Xtmx9tc8xLlm8CE1+f9TzEhHZRsp4CFdJuJjSXeFSav4URu02oPuM6nSjZlpBlZYuRHl5Zb/Xq6qqsWLF0mELatHCBeHtz60N2OPxJ38RvZBERHYhksKX7lwXAg430NkK2fw3yJYaoKtNdbqoMK2gSoqLsG3bW/1eb2pqxrq13mH3X7FiKWprfbGIRkRkH2mTINwXQUwqAQDIls8g6/8MnIm/K0mmFJTD4QAABAL9r3+GQiE4nRnDHqPA60VVVTWKiwqh6y4EAkG43S5s37GTEySIKM4JIHs2tMkLgfH54VLyvQ15cj8gz6oOFzOmFFQkBeRwOIYsGqczA06nE4FAoPcyn8PhwE8efQT3P/BQBCUlRhKZLIVjZ08ctzFLdkDoF0K4LwaSxkE2fwrpextoP3HeRtH+nM0Yt8imstviPihddwEA3LoOX1197+uhUAjVVTXDzuJ75OEfIUnr+6F/8OEufPDRrtgEpqjJzLLe1FcaHsdtbIzUCejKuRAydQJkUhqST/wVya0HIGQXMA7AuNh8vmaNW1vrieE3gk0KKhAIAgAO19b2e6/WV4eNG24csqAeevif0dXZGbN8FBuZWdnnbhwkO+G4jUHmdGhTFgPZs4HWwzDqdwCn69Ftxo+24LiZUlA9BTOUoS7R9bzX1NTc771AIBjRJUQr3R1NkTj/jJdjZx8ct1GZMBPa1MuBjCmQx6shP/ulyZMerDluphRUKBRCIBCE/pVLdAB6JzwMp7GpacAiiqyciIisRkBMnAeRexkwbgJkYwXkwd8AnfExRTwaTLvEV+vzDVgmbreO6q/c2zSQOl893G59wPcam5rGnI+IyBwCYlIxxNRLgRQnZGMlZGOFrW+ojRXTHvleXl6JkuKifq8XFxXi410Vw+7/8a4KFBcV9nu9pLgIFRXDr0JBRKSWgJhYBK3kf0BMWQTZUgXjkychv9jJchqEaQW1o2wndN3VOyMPALz5HgSCwX6rSDz5xM/wk02P9HmtonI3AsFg74oSPfvnez1ch4+IrC1nLrSSOyG8KyGPV8Oo+TXkkQ+Bs3zS91BMncX3400/xfp1q9HUFP7yz+3WB1woNhQMIRDs/73Upkcfw8YNN2DGjC9Xnrj7nvtiF5iIaCyyCiBcF0Bkz4I8VgF5rBw42646lW2YWlChUCiis537H3xo0Pd4tkRElpcxBVreciAzD7JxN4xP/oWX8UbBFvdBERHZQuoEiLxlEBPnQ7ZUQf7t34DOVtWpbIsFRUQ0VknjIKZeDjF5IWTz32BU/QoIxd/irWZjQRERjZbQIPSLIaZdAXS2wdi/FWjtv+INjQ4LiohoNLK80PKvBZLSIBt2QDb/DVZahSEesKCIiEYiLQea5yogqwDyi/chG3cDBtf6jAUWFBFRJLRUiKmXQUwpBU4dhPHpfwAd1lpcNd6woIiIhiEmFkK4vw6kZMDY/zLQelh1pITAgiIiGkxaDrT8lUDmtPD3TP6/ANJQnSphsKCIiL5KJENMvRQi9xLg5EEYn/470HladaqEw4IiIjrfeC+0gtUAAOPAq8CpQ4oDJS4WFBERAKRkQHiuhnDmhleBOPIRIM14li0NhgVFRAlPuC6A8FwJhJph7HsZaG9RHYmQIAXl8eThbHcXfL461VGIyErScqB5VwMZkyEbtkP6/6o6EZ0nIQqKiKgvAZG7GCJnPtBxKjwJoiugOhR9RUIUVH19A7o6+WAwIgKQrkObsRZIHQ/D9xZwcr/qRDSIhCgoIiIILbwSRO6lkMerIPe+yIcHWhwLiojin2MyRO43IDKnwTjwCqeO2wQLiojil9Agci+FmHoZZPOn4fXzuLCrbbCgiCg+pbugzVgHpDp51mRTLCgiijMCwv11CM+VkMdr+F2TjbGgiCh+jJsAbcZ6IH0ijP2vcNVxm2NBEVFcEK6vQeRfA7TWwfjsP4GuoOpINEYsKCKyt2QHtIJVQNYMyLo/QTZ/ojoRRQkLiojsK2tGeCJEx0kYn/2ST7iNMywoIrIfkRRe3NWhQzZWQh79LwBSdSqKMhYUEdlLugvarOsAkQzj0B+A4FHViShGWFBEZBt9po/7/sibbuMcC4qIrC85HVrBWmC8B/LwNsjjNaoTkQlYUERkbePzoc38ZvixGFW/BDpaVScik7CgiMiiBMT0JRATZkP6/wr5xQfgRIjEwoIiIutJHQ9t5nVAWjaMQ78HTterTkQKsKCIyFomzA4/UDBwJHxvU3dIdSJShAVFRBYhIDxXQrgXQH5eBnlsl+pApBgLiojUG5cFbda3IbuCMPY8BwSOqE5EFpAQBeXx5OFsdxd8vjrVUYjoqybMhjZzHXC6AfLw63w0BvVKiIIiIgsSGsT0ZRCTF/GSHg0oIQqqvr4BXZ0dqmMQUY/UTGgzvwWMmwBjz/NA4HPViciCEqKgiMhCxnuh5V8FdAZgVP2Ks/RoUCwoIjKNyL0UYtoVkJ/vgDxWAd54S0NhQRFR7CWNCz+KPXM6jP0v81HsFBEWFBHFlsMNbfb1QHc7jOrNXEuPIsaCIqKYEZNKIApWQTZ/Bln3R0CeVR2JbIQFRUTRJ5Igpl4GkVsK6XsbsvlT1YnIhlhQRBRVRnIGtPnrASFgVD8LhBpVRyKbYkERUfRkTke793og0Bhehbz7jOpEZGMsKCKKivDj2K9G8om/oePQOwAM1ZHI5lhQRDQ2IgnCuwpi4nwYh/6A1LON6OD9TRQFLCgiGr3U8RCeqyAcbhjVTwNnWoCsbNWpKE5oqgMQkU05p0Mruh0wDBjVTwFnmlUnojjDMygiGjGhXwiRfy3ksXLIz98DlyyiWGBBEVHkhAbhuRrCdQHk4dchj9eoTkRxjAVFRJFJdoSXLBo3AUYN72+i2GNBEdHwHG5oc24EOk7DqNrMR2SQKVhQRDS0CXOgzVwHebz63Hp6vL+JzMGCIqJBiWlXQOReCqP2DaDlM9VxKMGwoIioPy0ZYsY6iKwCGPteAk7XqU5ECSghCsrjycPZ7i74fHWqoxBZX2omtNk3AkkpMKqeAjpOqk5ECSohCoqIIpSRG54MEWqCsfd3wNkO1YkogSVEQdXXN6Crk79oREMREwshZqyFbPoLZP2fwZtvSbWEKCgiGpqYtgQi9xvnHi74ieo4RABYUESJTUuGyF8J4ZwCY++LQFu96kREvVhQRIkqxQltzk3hyRD7tgIdp1QnIuqDBUWUiBzucDmdaYGx9wVOhiBLYkERJZrs2dBmXgfZ8hlXhiBLY0ERJRAxpRRi+nLIhu2QjRWq4xANiQVFlBAEhGcFRPZcGAdeBU4dVB2IaFgsKKJ4l5QKbda3gXQdxr4twBm/6kREEWFBEcWz1PHQ5mwAYMCofhroalOdiChiLCiieJUxJTxTL3gMxsHfAUaX6kREI8KCIopHPTP1/J9A1r8LLltEdsSCIoo3rguheVdCNvwZsrFSdRqiUWNBEcUNAeG5CsJVAuPAK8CpQ6oDEY2JqQXlcDiwYvkSNDU1w+nMgNutY8vWV0d9vOXLlsDvb0ZVdU30QhLZkZYMbeZ1gHNaeE294DHViYjGTDPzh917z10or9iNisrd2FG2E9t3vIcHH/jhqI7lcDjw93fcGuWERDaU7IA2/++AtIkwap5mOVHcMK2gli9bAn+TH35/c+9rPX9etHDBiI+3uHQhAoFgtOIR2VPaRGhFtwFnu2DUPAt0tKpORBQ1phVUaelC1A7wyPWqqmqsWLF0RMcqLirkZT0i5zRohX8H2fY5jH0vAWfbVSciiirTCqqkuKjP2VOPpqZmFHi9IzqWrrsGPBZRwsiZB23+LTCOfAx5+DVAnlWdiCjqTJkk4XA4AACBQKDfe6FQCE5nRsTHWr5sCXaU7YxSMiL7EZMXQeStgKx7B/D/VXUcopgxpaAiKSCHw4FQKDTkNmM7cxKj3I/U49j1EJ4rIfSLYBz4zbkFX6382Vg5Gw3OjHGL7MZxW90HVVxUOKqzp0ce/hGStL4f+gcf7sIHH+2KUjKKlcysbNURLEGKJHTmXomzjlyMa3gNSbIFsPBnw3GzJ7PGra31RETb2aagFi1cgF3lo7sr/qGH/xldnZ1RTkSxlpmVjbbWk6pjqJecDm32DUBKBoyqZxDqsPZnwnGzJyuOmykFFcl08KEu7/V8hzXcJcChcS0yezn/jDeBxy41C8JzFSAEjJpngO4zqhMNg+NmT9YcN1MKKhQKIRAIQtd1+Orq+7yn665hC2z9utXIcDgwY0bf2X5OZwbWrl2F4uJCVFXVcOo5xReHG9rcjUDgCIw9LwCyW3UiIlOZdomv1ucbcLKE262jephiGWw5pHVrV2PbtrdYTBR/xudDm30j5PEqSN87sNLfaonMYtp9UOXllSgpLur3enFRIT7eVWFWDCLLExMLoc3dCHnsY0jf22A5UaIyraB2lO2Errug667e17z5HgSCQVRU7u6z7ZNP/Aw/2fRIRMft+X6KKC7oF0PMWA/pexvyyIeq0xApZeosvh9v+inWr1uNpiY/gPDlvU2PPtZvu1AwhEBw8O+lNm64AV5vfvjPG29ASXEhXn/jLa4uQbYm8q4MPypj/1agtVZ1HCLlTC2oUCgU0eM17n/woSHfH8sjOogsR2gQM9ZBZBXA2LcVCB5VnYjIEmxzHxRRXNJSoc2+HkjLCa9G3h7ZDYxEiYAFRaRKSga0ORsAoYXLqav/WpVEiYwFRaTCuAnQ5n0H6Dgdfjz72Q7ViYgshwVFZDaHG9rsGyBPHYKsf5ePyiAaBAuKyEzjPeEbcFuqIOv+CN7jRDQ4FhSRWbLnQpt1HeSRjyCPfKA6DZHlsaCITCD0iyDyr4Ws/yNk019UxyGyBRYUUYyJqZdBTL0MxqHfAyf2qo5DZBssKKJYmnoZxJRLYOzbApyuU52GyFZYUESxcP7qEDW/Bs40qk5EZDssKKJo01LCq0OkT+LqEERjwIIiiqakNGhzNwBJ42BUPwt0talORGRbpj1uQyWPJ6939XOimEnJhFb4PQACRs1zLCeiMeIZFFE0pOVAm3sz0H4cxoHfAEaX6kREtpcQBVVf34CuTq51RjHimAytYBVk4Ajk4dcAaahORBQXEqKgiGImMw/anJsg/Z9ANvxZdRqiuMKCIhqtCbOgzf425NFdkF/sVJ2GKO6woIhGQUwsgpixDrJhO2Rjheo4RHGJBUU0QsL9dQjP1ZC+NyGbP1UdhyhusaCIRqB3Xb2DvwVO7lcdhyiusaCIIiSmLYFwfY3r6hGZhAVFNCwBUbAaInsOjH0vAcFjqgMRJQQWFNFQRBK0md8EMqfD2PNr4Eyz6kRECYMFRTQYLQXa7BuAtJzwoq8dp1QnIkooLCiigSSNgzZ3Y3jR15pnga6A6kRECYcFRfRVyQ6I/GsAoYUv63WfUZ2IKCGxoIjOlzoe2rzvAJ1tMPY8z0VfiRRiQRH1GJcNbf53gVBzeEVy2a06EVFCY0ERAUC6C9q870C2NUAe+gNXJCeyABYUUUYutLkbIU/ug6x9C4BUnYiIwIKiRJfpgTb33OMy6t9VnYaIzsOCosQ1YSa02ddDHv0Y8ov3Vachoq9gQVFiyp4HbfoSyIb3IBvLVachogFoqgMQmU1MKoE26zrIY7tYTkQWlhBnUB5PHs52d8Hnq1MdhRQT7oshPNdAHn4d8ni16jhENISEKCgiABBTFkNMX8pnORHZREIUVH19A7o6O1THIIXEtCsgplwCY/8rQOth1XGIKAIJUVCU2ITnKgjXheFnObU1qI5DRBFiQVFcE9OXQ2TPhbH3BSB4VHUcIhoBFhTFKQExYx1EVgGMPS8C7X7VgYhohFhQFH+EBm3mdUDmtPDjMtqPq05ERKPAgqL4IpKhzb4eSJ8Eo+Y5PgWXyMZYUBQ/tBQIz9VAWnb4zKnztOpERDQGLCiKD0mp0OZsBJLHwaj5NdAdVJ2IiMaIBUX2l5QGbd7NAET4Kbh8RDtRXGBBkb0lO8KPaDe6Yex7ETjLG7KJ4gULiuwrJQPavFuAriCM/S8DRqfqREQURSwosqfU8dDmfRfytC/8oEGjW3UiIooyPm6D7GfcBGiF3wPOtEDW/YnlRBSneAZF9pI2Edq87wKBL2Ac+j0gDdWJiChGWFBkH+mu8GW91lrIw68DkKoTEVEMsaDIHhyToc37DuTJfZC1b4HlRBT/WFBkfRm50PKuhGypgqz/k+o0RGQSTpIga3NOC1/WCx1jORElGJ5BkWWddeRCm7YasrES8vMy1XGIyGQsKLKm8V50TF8DefRjyC/eV52GiBRIiEt8Hk8evN58xSkoYhNmQpt7E1Ja/hvyiw9UpyEiRXgGRdaSPRti8iLIz99DSmg/2lXnISJlEqKg6usb0NXJRUQtL2cutJnfgqx/F7Lpv4GsbNWJiEihhCgosj4xsRBixnrIuncg/X8FIFRHIiLFWFCknJhUAlGwBtL3JmTzp6rjEJFFsKBIKeG6AMK7CrJ2G2RLleo4RGQhLChSZ1IJRN4KyMN/gDy+R3UaIrIYFhQpIdxfh/BcDePAb4BTB1THISILYkGR6cTkRRB5y2EcZDkR0eBYUGQqMWUxxPSl586cDqqOQ0QWxoIi04jcSyGmXQZj/ytA62HVcYjI4lhQZAox9XKI3G/A2PcycNqnOg4R2QALimJO5F4GkXsJjP1bgNP1quMQkU2woCimxPRlEBNmwti3BWhrUB2HiGwkIVYzJzVE3goI9wIYvndYTkQ0YjyDopgQeVdC6BfC2PciEDiiOg4R2RALiqJOeK6CcF0AY++LQPCo6jhEZFOmFpTD4cCK5UvQ1NQMpzMDbreOLVtfjXh/b74HxcWFcDqd8Hrz4W/y46WtryIUCsUsM42M8FwN4SqBsfcFIHhMdRwisjFTC+ree+7C5qeehd/fDADQdRcefOCH2PToY8Pu6833oKDAi21vvN372to1K/Hs0/+J/3n3/+49Jqkj8q8Fxk2AsecFINSoOg4R2ZxpkySWL1sCf5O/T5H0/HnRwgXD7l9cXIgdZTv7vLbtjbdRW+vDHbd/P5pRaRSEdyXExELIz8tYTkQUFaYVVGnpQtT66vq9XlVVjRUrlg67//LlS1FcVNh//+oalBQXRSMijZLwrobImQdj7/NAqEl1HCKKE6YVVElx0YCX4ZqamlHg9UZ0DF13RTsWjZEoWAORM+fcZT2/6jhEFEdM+Q7K4XAAAAKBQL/3QqEQnM6MYY9x9z33Dfi6W9dRWxvJ0jl8hHh0CYjpSyEmzAqX05lmxO4z5tjZE8fNnswYNxnRVqYUVCQF5HA4Rjwbz+FwoLR0IX7++L8Oud0jD/8ISVrfD/2DD3fhg492jejnUZiEQOeUZTibMQ3jGl5DUmo3kJodk5+VmRWb41Jscdzsyaxxa2s9EdF2tr4P6uYNN6C8vBIVlbuH3O6hh/8ZXZ2dJqWKdwJixloIx3QYNS8g1N4Ss5+UmZWNttaTMTs+xQbHzZ6sOG62LajiokIUFHhx/4MPRbhHZKeUNJRz5ZTlhbHn10D78Zj+rC9x7OyD42ZP1hw3UwoqEAgOu81ILu85HA6sXbsKP97007HEohEREDPXQ4zPh7Hn+RiXExGRSbP4QqEQAoEgdF3v956uuyIqsPPdecet2PzUs1xBwjQsJyIyn2nTzGt9vgEnS7jdOqqrayI+zsYNN+Clra/0mbLuzfdEJSMNRITX1hvvMeGyHhHRl0wrqPLyygFvqC0uKsTHuyoiOsbyZUuwa1dFv/upiov738BL0SAgZn4TInsWjJpfA+2RzbwhIooG0wpqR9lO6Lqrz8223nwPAsFgv1l4Tz7xM/xk0yN9XisuKoTbrffu5833oLioEIsWLoB7gEuHNFbnyikzD8bel4AOa83uIaL4Z+osvh9v+inWr1uNpqbwigNutz7gQrGhYAiBYN/vpX704P8FAKxbu7rf9tu3l8UgbSITELOug3BOC1/WYzkRkQKmFlQoFIro8RoDTR2/ccMtsYhE/ZxfTs+znIhIGT7ync7DciIi62BB0TkComA1REomy4mILIEFReg9c8oqgHH4dZYTEVkCCyrh8bIeEVkTCyqhsZyIyLpYUAmL5URE1saCSkgsJyKyPhZUwhEQ+ddApE1iORGRpbGgEkrP2npzYBx4leVERJbGgkoY58ppfB6MvTxzIiLrY0ElhJ7nOeWde54TVyUnIutLiILyePLg9eYrTqHKVx82yHIiIntIiIJKXAJixjqWExHZkqmrmatSX9+Ars4O1TFMdq6csrx8TDsR2RLPoOKSgChYBSQ7WE5EZFssqDgkClZDTJgNWf8uy4mIbIsFFWdEwWqI7NnhqeTtLarjEBGNGgsqjgjv6vBNuHueB86wnIjI3lhQcUJ4V0HkzIGx5wXgTLPqOEREY8aCigPCuxIiZx6MvS8AZ/yq4xARRQULyubE1CsgcuaHyynEciKi+JEQ90HFK+G5GmJS8bnvnFhORBRfeAZlU8JzFYSrBMbeF1lORBSXWFA2JPKuhHBdEC6nUKPqOEREMcGCshmRtxxCvzBcTsFjquMQEcUMC8pGxLSlEPrXYex7CQgeVR2HiCimWFA2IaYtgRjvCZdT4IjqOEREMceCsgEx9XKIKaUwPi8DAl+ojkNEZAoWlMWJqZdB5F4CY99WoK1BdRwiItOwoCxM5H4DYuqlMPa/DLTVq45DRGQqFpRFiSmLIaZdAWP/K8DpOtVxiIhMx4KyIOFeCDFxfricWmtVxyEiUoIFZTHCvQDCswLGF+8DrYdVxyEiUiYhCsrjyYPXm684xfCE+2IIz1UwDvwWOHVIdRwiIqUSoqDsQOgXQniugXHwt8CpA6rjEBEplxCrmdfXN6Crs0N1jEEJ1wUQ+SthHPo9cHK/6jhERJbAMyjFxKQSCO8qyMOvASf2qo5DRGQZLCiVsudB5F8Nefh1yOM1qtMQEVkKC0oRkTMf2qzrYPjegTxerToOEZHlsKBUyJ4LMfObkL63AZYTEdGAWFBmy54Nbda3IOv+CNn8ieo0RESWxYIy04SZ0Ey6ew4AAAqgSURBVGZ9G7L+XUj/X1SnISKyNBaUWbIKICaXQjbsgGzarToNEZHlsaDMMD4f2pwbgdZayMYK1WmIiGyBBRVrmXnQ5twE+cWHkMc+Vp2GiMg2WFCx5JwObe5GyGO7II9+pDoNEZGtsKBiJWMKtNnXQzZWQH7xvuo0RES2kxBr8ZkuYwq0ed+FPLoL8uiHqtMQEdkSz6CizeGGNu87kM1/YzkREY0BCyqa0l3hcmqphqx/V3UaIiJbY0FFS9qk8GW9E/sh695RnYaIyPZYUNGQlgNt/nchWw9B+t5UnYaIKC6woMZq3ARoM78FeboO8vAbqtMQEcUNFtRYpGZBm38LZPAo5KHXAEjViYiI4gYLarRSM6HNvwUINp77zonlREQUTSyo0UhxQpt3C3CmBcbB3wHSUJ2IiCjuJERBeTx58Hrzo3OwlAxo878LdJyCceBVQJ6NznGJiKiPhCioqElOh8i/FugMwDjwCsuJiCiGEmKpo/r6BnR1doztIElp0OZ9B+jugLF/K2B0RyccERENiGdQkUgaB23ezYBxFsaBlwGjS3UiIqK4x4IajpYKbe5GAALGvpeAs52qExERJQQW1FC0FGhzNwBaCoy9LwJnx3iZkIiIIsaCGoxIhjbnJiA5/Vw5tatORESUUBJiksSIiSSImesh5VnIvS8A3SHViYiIEg7PoL5KaNBmfRsiY3J4bb2uoOpEREQJiQXVh4A28zogww1jz/NAV5vqQERECYsF1UtAzFwPOKeFy6nztOpAREQJjQV1jihYAzE+H8be54GOU6rjEBElPBYUAOFdCaRkwNj7AtB+QnUcIiICCwrCczVEznzIz3cAZ1pUxyEionMSuqBE3goI19fCK0SE/KrjEBHReRK2oMS0JRDui8PlFDymOg4REX2FqTfqOhwOrFi+BE1NzXA6M+B269iy9VXT9u8hpl4GMaUUxr4tQODIiPcnIqLYM/UM6t577kJ5xW5UVO7GjrKd2L7jPTz4wA9N2x8AoF8EkXspjP0vA20NI/w3IDNdfuli1RFoFDhu9mTFcTOtoJYvWwJ/kx9+f3Pvaz1/XrRwQcz3BwDhXgAt/5rw85xO10WcndS4/DLr/cLQ8Dhu9mTFcTOtoEpLF6LWV9fv9aqqaqxYsTTm+wv9IgjPVTAO/pblRERkA6YVVElxUZ+znx5NTc0o8Hpjur/IKYLIvxbGod8DJw9EnHkwXm8+vN78MR/HrOPG8tixOm6seL32+hzsdtxY8Xrt9TnY7bix4vWOLa8pBeVwOAAAgUCg33uhUAhOZ0ZM9xf5V0Mefh04sTfCxEREpJops/iGKxAgXEKh0MCPtRjr/slHdsBoOwikjhv2OJEQWvhjS4nS8WJ93FgeO5aZQ6F2pKSmRvWYdvsc7HZcgONmx+MCZo+bRFfn8E8nj+vnQaWkhD/s61ddDuBytWFoxOr9Hbj51n9UHYNGiONmT2aP20vPPDlsScV1QYWCAbz12lacOtH/uyuyPuf4CQic5sK9dsNxsydzx81CZ1CBwPAP/Rvs8txY928/Ezr3Qchhj0FWItDd1cWxsx2Omz1Zc9xMmSQRCoUQCASh63q/93TdNWwBjXV/IiKyH9Ommdf6fANOdnC7dVRX18R8fyIishfTCqq8vBIlxUX9Xi8uKsTHuypivj8REdmLaZMkdpTtxIrlS6Hrrt4bbr35HgSCQVRU7u6z7ZNP/AyhYAj3P/jQqPY/f1HZSS4dWeOd2LL1lRj/G9JwvPkerF+3Bo8/+YsB3++7GLAT0/Py8Nxzzw+z3egXDabhefM9KC4uhNPphNebD3+THy9tfbXfd76R/s5x7Myh6y6sWP7lCjtuXcdLW1/pt9iB1cctacr0mQ/H/Kecs6u8EiuvvRq67kKBNx9z587Bv//Hr/ptd8VllyIQDOLDD/9rVPv/8L57se2Nt3DgwEE0+o/j84YG3HP3D/odj8xxx23fw8UXXYDSxYuQnT0BO8p2Drjd+ePm89WjtS2IH/zDbf3Gre92dWjy+3HP3XdxfKPMm+9BQYEXf/zTn1FVXYMPP/wvTJuWi//3o3/CBx9+hGDwy5L6ckwOodHfMujvHMcu9nTdhdJFC/Db372GquoaVFXXQEqJ/3vf/0J5RaWtxs3U1cxDoRC2bH0VO8p2YkfZzkEb+P4HH8KmRx8b1f7RWFSWomvz089h89PPoby8ctBtBhq3lpbjAPqOG8fXPMXFhf3+MrHtjbdRW+vDHbd/v/e1SMeEY2eO88+celRU7obTmdHnPTuMW9w9sHCsi8qSGoOPW02fceP4mmf58qUoLirs93pVdU2f74MjHROOnXkWly7q91ogEETGuWXjAHuMW9wV1FgXpSU1Bh83f59x4/iaS9ddw24T6Zhw7MyxZeurfb6/B8LfITmdGSiv+PL7ejuMW1ytJDHWRWVJjUjHjeNrrrvvuW/A1926jtpaHwCOnV3cecet2L69DFXnbsmxy7jFVUGNdVFZUiPSceP4qudwOFBauhA/f/xfAXDsrKxnBuaMggJ8vKuiz2xnu4xbXBXUQNpaT6iOQKMiw2M3fYrqIHSemzfcgPLyyn63dnxJ8nfOInx19fDV1UPXXVi3ZhWczoxBZ9Baddzi7jsoIoqN4qJCFBR4B72PjazJ72/G5qefw8YNN2LtmpWq44xIXBXUWBelJTUiHTeOrzoOhwNr167Cjzf9tM/rHDv72FH2Hm7eeBMA+4xbXBUUF5W1p0jHjeOrzp133IrNTz3b7z9GHDtrcTgcePCBH8Kb7+n3Xs9EB1132Wbc4qqgAC4qa1eRjhvH13wbN9zQb5mc8/8DyLGzDrfuQklx0YCF4nQ6AXx59mSHcYu7guKisvYU6bhxfM21fNkS7NpV0e8+mOLiL2/g5dhZh6+uHq9ve3PASSzFRYX4rKq69yzYDuMWdwW1o2wndN3V5wbDwRaVJXNlZDgGfS/SceP4mqe4qBBud/hv4t58T3jaclEhFi1cAPd5f0Pn2FnL4cM+LF+2pM9rxUWF0HUdm596tvc1O4ybuOiSa6zz+MQocTgcWL9uNZqa/ADAFZMVu+O278HpdKKoqBBOZwY+q6qGv8mP7Tveg6+uvne7SMeN42uOV7b2X0m+x/btZdj89HO9/59jZy3efA8WL16EQCAAp9MJt67jl5ufGXAVeiuPW1wWFBER2V/cXeIjIqL4wIIiIiJLYkEREZElsaCIiMiSWFBERGRJLCgiIrIkFhQREVkSC4qIiCyJBUVERJbEgiIiIktiQRERUUwVFxXila3Pw+EILxjtcDjw9Ob/wKKFC4bcjwVFREQxVVVdg/LyStx5x60AgJs33IBan2/Y1dBZUEREFHO/3PwMSksXYu2alSgtXYTHn/jFsPuwoIiIKOZCoRB+tfkZ3LzxJmzZ+kq/R38MhAVFRESm6Hloac+DMIfDgiIiopjTdRdu3ngTfv74v2Ld2tV9ntA7GBYUERHF3L1334XXt72Jisrd2L69DPfefdew+7CgiIgoptauWQlHhqP3MfGbn34Ouq5j+bIlQ+7HR74TEZEl8QyKiIgsiQVFRESWxIIiIiJLYkEREZElsaCIiMiSWFBERGRJLCgiIrIkFhQREVkSC4qIiCzp/wMsB5GVNKkRWQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 460.8x403.2 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<sympy.plotting.plot.Plot at 0x22aae4fc588>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import plot\n",
    "plot(-usol, (x, 0, geom.L))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "uplot = lambdify(x, usol, \"numpy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GPXi():\n",
    "    def __init__(self, ordr):\n",
    "        self.xi=leggauss(ordr)[0]\n",
    "        self.wght=leggauss(ordr)[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "class basis(): # defined on the canonical (master) element (1D : [-1,1] ) \n",
    "    def __init__(self,deg,basis_type):\n",
    "        deg = int(deg)\n",
    "        if basis_type == 'L':          # 1D Lagrange basis of degree `deg`\n",
    "            z=Symbol('z')\n",
    "            Xi=np.linspace(-1,1,deg+1)\n",
    "            def lag_basis(k):\n",
    "                n = 1.\n",
    "                for i in range(len(Xi)):\n",
    "                    if k != i:\n",
    "                        n *= (z-Xi[i])/(Xi[k]-Xi[i])\n",
    "                return n\n",
    "            N = Array([simplify(lag_basis(m)) for m in range(deg+1)])            \n",
    "            dfN = diff(N,z)+1.e-25*N\n",
    "            self.Ns=lambdify(z, N, \"numpy\")\n",
    "            self.dN=lambdify(z, dfN, \"numpy\")\n",
    "            self.enrich = 1  #enriched dofs (not needed for lagrange elements)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ N_1 = \\frac{\\xi(\\xi - 1)}{2}\\\\\n",
    "N_2 = {1-\\xi^2}\\\\\n",
    "N_3 = \\frac{\\xi(\\xi + 1)}{2}\\\\\n",
    "x = \\sum_{k=1}^3 N_k(\\xi) x_k \\\\\n",
    "Je = \\frac{dx}{d\\xi} = \\sum_{k=1}^3 \\frac{d N_k}{d\\xi}(\\xi)x_k \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x22ab06a9408>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAFYCAYAAACI+iGlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeXBc133g++/degNAcG2AJAigwZ0EqY27SYoiKcmWJUpeJI/sScZxYs9SeRnnvXp5lThVtiuRXbFrauz3ZpLUyB45mUiyFY9sS6ZsS9xESuImcQNIgiCJBkAsRIMkSKDR6PWe90eTNCliaQDdffvePp8qlGXg3r4HF83763PO7/yO8uCGTwokSZIkKQ9UqxsgSZIkFQ8ZdCRJkqS8kUFHkiRJyhsZdCRJkqS80a1uQDb4SkpJJOJWN0OSJEkCDMNFZDA87M9sH3QWL1nC+i1PWd0MSZIk6Q6v/a9/GDbw2D7opFJJAA7ue5NUKklbW/s4X0Hh29/6a775rb8Fxp89XlNTDTCB62bvNaw+f7L3MBttsPv58h5m59+S1f+erT4/O68xsXt467qaprN+y1Mjjj7ZPujckkolSSWTJOLjHWZT0FTl5nnjf5OmkumgN/7rZu81rD5/svcwG22w+/nyHmbn35LV/56tPj87rzGxe3jrumOxfdBpa2tn483/ndybdWKCwVbLX8Pq87PB6t/B6vOzwerfwerzs8Hq36EQnieTva7hcrFxlONk9pokSZKUNzLoSJIkSXkjgw6w/8BBq5tge/IeTp68h9kh7+Pk5fIeyqAD7H9PvkknS97DyZP3MDvkfZy8XN7DjBMJArU1PPP0U/zXH/63cV3A5/OxfdsWenp6KS0toaLCzyuvvjbh4yRJkiT7GjPofPWPvwyAv8JPaUnJuC/w51//U1780UuEQr3p1/HP4ht/9Re88J3vTeg4SZIkyb7GHF578cc/4cUf/4RDh46M+8W3bd1CqCd0O5AAt/977ZrV4z5OkiRJsreczumsW7eGlmFyxhsaGtm+/ZFxHydJkiTZW06DzsoV9Xf1Xm7p6emlLhAY93GSJEmSveWsIoHP5wMgHL634FskEqG0tGRcx41Nufk1TrqXlMcPN/omdr70MfIeZkJRBW6vwOURqBpouqCkbAiXL0kqpWCmIB5ViA0pCFPe04mR920ianWVQQQDwMTu4ejn5CzoZBIsfD5fxsdFIpFRj/n2t/4aTb37l91/4OCYqX+JaSuJT7+f0mgvygRrXklpZeXTrG5CgRF4SpOUz4xTUp7AV56gZEoCd0kKl8e868hUUkGYQyiaQNPufpV4VCU6qBO5oTPYbzB43aD/ipvooIZ8sA5PvhcnRkPwVT3GrlSKDzO4h5s3rmfzpvV3fS9lCprb+0c8x/a112755rf+dmK118IH0WatJaz54VpT9htWJMrKpzFwo8/qZlhMUDbdxF+dxD8vxbTKFG6vIDKg0H9VI9yncrlVI9KvE42oxCIK8Wi6VwPKHfcw3ftxeQUer8DtM/FNEZROS1I2LcacBSbeUkF0UKGvRyN0SSPUrhPuU5FBSL4XJ2ONW4MSgxNCy+ge7tz5Fjt3vnXX9wyXiy995c9GPMcxQSddDXUCPZVUHP36GUTlWsxrZ7PequJw54OuuHqLiiKYMTfF3AVJZtcl8ZYKrodUQpc0Wk+7uXZZIxYZ/9SpmYJoWCEaVhhu6tVTajK9IsX02SlqlyW47+EYkQGF7hadzvM6V7s1EMUYgIr3vZgNW7w670WTJI1b93Ei93D0c3IWdMLhwTGPGWvIbLzHTZTed4rE/D+Akjkw2JXTa0nOMGVGitrlCeYuSmK4BJdbdRrfdxNq14gP5b7QRzSs0hVW6bpoAOD2mfir08HvE58ZIj6k0NGs03baYKBPG+PVJAkWGiqVmsI/DCXByN11chZ0IpEI4fAgfr+fYGvbXT/z+2fdDkqZHpdLamIAca0JZfY6xIXXc349yZ5UTTBvcZLa+jjTKkxC7Rqn3nVzuVUnlbC2VxGLqFxqUrnUZKC7BLPrklQvTbD9DyJc7dIINhp0NOsyKUEa0SNegw9jKQYElOXwOjkdXmsJBodNFKio8NPYeHrcx+WS6D6EuvzLiPZdEB95EkwqPi6PoLY+zvz7EgC0Nhoc/Y1BZKAwSxcm4wqXmgwuNRmUlJvULE+wYlOM5RtiXDxhEGx0kYzL4CP93ixVYZmh8nfXozm/Vk7/1Rw6dISVK+rv+f6K+uV8cPDwuI/LqXAnhLtQKmUFBCnNcAuWbYjx+JfDzFuS5MxBN7/7SQlnD7sLNuB83OANlTMfuPndSyU0HXFRW5/gk18Js2RtDN0l5zyktC1eneaESXcq9++JjP/llJT4Rv35D3/wfb77wrfv+t7uPfvw+2fh98+6/b1AbQ3hwUEOHzk67uNyTXQfQvE/BGoOBzSlgqcbgsVr0sGmsjbJh2972P0vPtrOGJgpe/YQUkmF1kYX7/yvEo7v9lC1MMnjXw6z6KEYmi6DTzHzKbDGo7NvKJGX62VU8LO0tJT6+uWUlpbwjb/6C0I9IXbt3nvXHExkMEJ48N75l7954e945ukn6ekJAekhs+GKeGZ6XC6Ja2dRaraj+B9AXB5/rTnJ7gTVS5Ms3xAjEYfjezx0ntdxVBqyUOg8b9B5QWfe4iRL18YIrExw+n03Hc0O+12ljHzCo9OXEpxNmGMfnAVjBp0Xf/yTjF7oL7/xzWG/H4lEMtqiINPjcksgug+jVK5FXD6KTLksHtMqU9z3cJSScpOmw25aGgxnT7qL9LxP53md+ffFuX9rlLqVKiff9XCjV2a7FQsd2Ow1eGswnrennT0GpvNI9B4H3QvTl1jdFCkPdJdg5cNRNn8+Ql+Pxjv/XMLFky5nB5w7mCmF88fcvPPPJYT7VLY8F6F+Y1QOuRWJB90aKnA0lsrbNR20ODRLUnFE6Bjq7PVysajDVdYmuX9rlERM4cDPfVy7XLyf8GMRlWO7vbQ1JXlwa5Q585Oc2Osh1C4fEU72iNfgwFCCZB6v6Zh3VE1NNalkkuAwWySMl7h8BKVyLZRWQbhj8o2TCopuCFZsjjFvSYJzR1w0f1Q8PZuxXO3U2f1KCYtXx1n/1BDBRoPT77tJJeX9cZolhsosTeG9aHZCTiBQC4Cmjx5WHBN0sirej7h2Jt3bOf+vVrdGyqLps5OseixKKqmw72c++q8Ub+9mJGZK4ewhN5eDOqseG8I/L8WHb3u4HpL3ykke8RocjSYJ53kk1TFzOm1t7Vnp5dwiug7C9MXgltVqnUGwaFWMTZ8douuizt6fyoAzlr4ejT2vltDbofHwsxEWPBBHJtc4w2xNYbGhsjdLvRyAYLCVYLCVtrb2UY+TPZ2RRC5DfxvK7LWI1t9a3RppElwek1WPR5nqNzn4ay+hNvm2z1QqqXByn4dQu8aDj0aZMSfFsXc8JGRFA1vb6jU4kzAJ5WEx6Mc5pqeTC2b3QZRZ96ez2SRbmlaR4pHnI+iGYO+rPhlwJqi7xWDvqyV4S00eeX6Q8pn5y3aSsqtcVXjIrbE7kp/FoB8ng85orl+A2HWUilVWt0SagHlLEmz6XISuCzoHXvcxFJZv98mI9Kvs/7mPULvO5mcjzFlgzUNLmpyHPTodSZOLyfwsBv04+a9wDKL7IErFalDk+L9tKIL6jVEe2BrlxB4PDQc8MjstS8yUwom9Hhrfc7P68ShL1saQ8zz24VHSFQj2DOUzSfpuMuiMQVxpBATKrJVWN0XKgGYI1j85xLzFSQ687qO9SdbRy4Vgg4v3f+mlbmWCNZ+Komoy8NjBerfOoBCcjFs3PCqDzlhEKr1uZ/b6sY+VLOX2mWz6bARvmWDfz3z0FfFiz3y40qmz72c+ymem2PiZIVweGXgKmQo87NXZO5S0tG8qg04GRM9H4JoC0xZZ3RRpBKXTUjz8bIREXGH/z+X8Tb5E+lXe/VcfINj8+Qi+KdbME0hje9Ct4VYUDmcxTXoi5L/MTKSiN0vjbLC6JdIwpvpTPPz5CFe7NT74lVduUJZn8ajKe7/w0X9N5eFnI0yZITPbCtFWr8GBaJK4xe2QQSdDovsQlM5Nl8aRCsbMuUk2fjZCe5PBR2/LhAGrmCmFI7/x0B3U2fS5CNMqZOApJEsMFb+msD9Pe+aMxjFBp6am+nbtn5yI9yOuNqLOkb2dQlFRm2TD00NcOOai4YAbuReMxYTCiT1u2k4bbPxMhJlV1g7jSL+3zWtwOMclbwKBWgKBWmpqqkc9zjFBJx9E10GYthg8M6xuStGbXZdg3aeHOP2Bm6YjMuAUDoXG9900f+Riw44h/NUy8Fhtnq6y0FDZa2Ga9J0cszy7ra2dRDzHo5VDIbh+AWXOekTLr3N7LWlEs+vSabqn9rsJNrisbo50D4VzR92YKYV1Tw5x6NdeuUWChbZ5dU7EU1wxc5uzdqv2peFysXGU42RPZ5zM7g9QZq4Eo9TqphQlGXDs4/wxF2cPuVn3pOzxWGWmqnC/S2N3gfRyQAad8etvg0gPyuy1Vrek6FQGkjLg2MydgWfWvMJ58BWLR7w65xMmlywqeTMcGXQmwOx6H8X/EGhuq5tSNGZVJVnzqSEaDsiAYzfnj7loOuJi3ZNDTK+UWW35UqrAWo/OngLIWLuTDDoTce0cJMIoFQ9Z3ZKiMK0yxbonh2g64qLllAw4dtT8oZuWky7W74gwRVaozouHvQahlOBsonB6OSCDzgQJRNdBlMp1shBojk2ZkWLDjggtp1w0fygDjp2d/sBFR7PBxmeGKJ1aWA9Cp/EosMmj845F2xeMRgadCRJXTpEuBHqf1U1xLF+ZyYanh+g8b3D6AxcyLdruFE7uc9PTrrHhmQhunww8ubLBY31hz5HIoDNRIoXoPoQyZwPyYZh9Lo9gw9ND9PVonNgn1+E4h8KxXR7C11U2PD2E7pJFQrNNBx65uX1BIYZ1GXQmQYQ+At2HMmOp1U1xFFUTrHsyQjwKR3/rASEDjpMIU+HITi8IWPvEEIoqA082rXZrKAVQ2HMkMuhMRiqO6DmKMucTVrfEQQSrH4/i8sChX/swUzLgOFEyofDBG15Kyk0e3B5FbgSXHQqwzWewbyhBYYYcBwWdnNdeG4HoPgzemVA+P+/XdqLln4gxfU6KD37lJR6VAcfJYhGVD37lo7I2yeLVVtc+dob7XBplisL7FvRyAgFZey0/khFE6Djq3NEKP0iZqFkeZ/59CQ7/2ktkQL41i0H4usqRt7wsWRNn7sLCy7Sym0d96e0Lhgq44+iYf9ltbe23a//km+g6mN7yoGyeJdd3gplVSe7fEuPYLg/X5I6fRaW3Q+fkPjcPPRplmlw8OmFLDZUKTWGfRYtBg8FWgsFW2traRz3OMUHHUvEbiKsNqHNkb2ciSspN1j4xxLmj6XUcUvFpPe2i5ZTBuieH8JYWYs5V4XvUZ3Awx9sXZIMMOlkiOj+AqQvAV2F1U2xFNwTrnhyi95JO0xG5+LOYNb7v5npIZe2nh1C1An9yFpiArhLQVfYUUGHPkcigky3RK3CtSWayjYvgoUejCAEf7fIg1+IUOaHw4W+9GC7BA1tlRtt4POoz+DCWoi/H2xdkgww6WWR2vYcyYxl4plvdFFtYvDrOzLlJDv3aSyohA44EibjCoZ1eZs9PMv8+mViQiTmawjJDZVeBFfYciQw62TTYDTeCsreTgYqaJEvWxDnyWy+Rfvk2lH5v4JrGR297qN8YY+bcwh8ustpjPoNT8RQ9qcLv5YAMOllndh5Ib/LmmmJ1UwqWr8xk1eNDnDnkoveS3FFSuld3i0HzRy5WfzIqa7SNwq+lN2l7uwALe45EBp1sG2iHcOfNmmzSx6maYM0TQ1zt0jn/kUwckEZ29rCL/qsqaz4VRVHs8Sk+37Z7Dc4mTDps0ssBGXRywuw8gOJ/AIwSq5tScFZsiuHyCD56WyYOSGMQCkd/56Gk3GTZelmx4OOmqwqr3fbq5YAMOrlx4yJEelFmr7e6JQWlalGCmmUJDr/lJRGXAUcaW3xI5chvPCx4IM7sOns9XHNtm1enJWESLKCtqDMhg06OmJ0H0juLah6rm1IQSspN7t8a5dR+Nzd6ZcUBKXPXunVOf+Dmwe1RuXD0pimqwnqPzu9skrF2J8cEHasKfo6o7xzEbqDMXmt1SyynaoLVnxqip02ntVFWHJDG78Jxg2uXNVY9Lud3ALZ6dTqSJs0FtBV1ICALflpOdB5AqVwLmtvqplhq+Yb0PM6J3XIeR5oohY/e8VA61WTJmuKe3ylVYKNH53c2m8u5xTFBx8qCnyMRV89AIoxSsdrqplimojZJ3coER38r53GkyYkPqXz4Ow+LVsWZWVW863ce8Rr0pASnC6iXA7LgZ4EQiM73UGavA7X4hpXcPpOHHo1y9pCLPlk5WsqC3g6d88dcrHosiuEuvmE2nwKbbNzLARl0ck5cbYRUDKVildVNyTPBg9ui9F9VaT4m1+NI2XP2sIuhsML9RVifbYvX4JopaIjbdwuInAcdv39Wri9R2ISJ6HovnT6tFs/q+8CKBNPnpNLrcYQcVpOyR5gKH73tpbI2ybwlxTPM5lFgs0fn7UjC1qE2o6egz+dj+7Yt9PT0UlpaQkWFn1defS2jC2zf9giBQC27du0lEokQDofvOSbY2gbAd1/4Nrt27+XgoSMArKhfzsoVy3n51deIRCKZ/UYFSPSeRJm7GcX/IOLyEaubk3Ol01LUb4xxfLeHobDsTEvZF76u0nDAzX0PR7naqRXFTrObPToDpuC4jXs5kGHQ+fOv/ykv/uglQqFeIN17+cZf/QUvfOd7Y55b4fezckU9K1fUD/vzXbv28OKPfwJAXV2Ar9UF+NpXvwLA5Z4eXvjO92wdcICbvZ33UeZuQvR8BMLeb5rRKKpg1WNRult0uSGblFOtjQaVgSQPPRblwOteR/eoPUo6geD1wbitezmQQdDZtnULoZ7Q7YAD3P7vtWtWc/jI0VHPD4fD/NEf/4dhA8eOp57gjTffuv3/X37lpxw6fJQKv59wOHy7B+QEInQcZe7GdG+nZ/R7ZmeLVsXxlAje/6VcFCvlmsLxXR62/dsI8+9LcPGEc+cON3l0BoXgo5j9P7CO2Sddt24NLcOkIjc0NLJ9+yNjXqAnFBo24KyoX04weG9QCYV6aWg87aiAA4BIITrfR5n7CVCcmclVPivFktVxju32kIg591OnVDhiQyon9rpZviFG6dTCSiHOFrcCW70Gb0cSOOE3HDPorFxRf1cv55aenl7qAoExL3BnT+ZOgUANDY2nM2iic4jQMUBJFwN1GEUVPPToEO1nDUJtxZMwIVmv64JBd4vOQ48NgQOrFWzy6AwJwYcO6OXAGEHH5/MBDDv5H4lEKC2dWBXljw+rfVygtoYV9csn9NoFTaTScztznNfbqbvvOoZL0PBecVdfkKxxcp8b3xTBwgedVa3AhbN6OTBG0MkkqNwKTJkK1NbQ03Nvzwlgfl0da9esZjASoaHxNF98/jm2bd2S4SsrE/zKxmtk/iV6joOi3uzt5P56+fia6k9RU9/Psd1eknHV8vbY8+tOVrfFfl/xqMaJPR6WrotRUn4r8Fjfrsl+bfIaRIXgSMzM43XvlI3XuFvex0Geefop/usP/9uwP9u1e+9dQ26vvPoaP/zB92/P84zm29/6azT17l92/4GD7H/vYEbtKiufltFx2ZC4doJk1SY8sVYUYe/PL4oiWPVYN13nS4n2z6Cs3OoW2Vs+34dOE74KobZelm64ymB/pe2z2QwE2/QYb6d0SsrH9+E+GzJ5L27euJ7Nm+7ewiVlCprb+0c8J69BZ6yFosMFlsaGdI/nL7/xzVHP/ea3/pZEfGJd67LyaQzc6JvQuRPS/x7q9PsZdNemU6htbNGqGLoryfmPKvN7Dx0o7+9DBzq2W+XRP0gyq6aHlpP2Hurd5tWJqRr7rw/kfWgt0/fizp1vsXPn3VMlhsvFl77yZyOeM2rQCYcHx7zoeNbQbN/2CBdbWjI+HqAl2Mr27VszOFIwsZIYd34aytMkpEjcntsRoeO2XbdTNi3FkrUxjrzlJZW4NVLrvInc/LDgfehA8SGVc0ems3zDFS636LZdNOoiHXTeGIxj5v39MNn34ujnjPoXSVcQGMTv99/zM79/VkZB6U5r164ecT7nhz/4PmvXFE81ZtFzjN/P7diR4IFtUbov6lwOykWgUuHoCfq40qFz/zb71mbb7NVvzuXY8wPpaMb8GNASDA6bUFBR4adxHCnPPp+PyoqKEXtGkcEIwdbWYa/T0hLM+Dq2IZLpCtRzNmLHTLa6lQnKppmcfNfeQxiSEymc2OthemWKahvWZru1Lud3DspYu9OYQefQoSPDlrBZUb+cDw4ezvhCFTfnc4ZLvwY4eOjwsOuB1q5dza7dezO+jp2k1+2A4n/Q4paMj6fEZNmGGA3veYgP2XP4QnK2obDKmQ/c1G+K4fLY69G9+ea6nKMO7OVABkFn9559+P2z7koCCNTWEB4cvKcEzg9/8H2++8K3h32dW0N0gyP0dBoaTt+THr3jqScIhXrZvWffWM20J5FKV6Ce+wlQ7LOgcuXmGNdDGu1n7dNmqfi0NBhE+hXqPxGzuikZc3ovBzLMXvubF/6OZ55+kp6eEJAe8hqu2GdkMEJ4cPh5nlAoRDg8OOI80K2yN198/jkASnw+ekKhjIqK2pkIHUeZ8wmUCntUoK6oTVJZl2TPKz7GyseXJEsJheN7PGz5QoS2piRXOwv/Q9LDN2usOaX6wHAy+itEIpGMtjIYLa052NrGH3/1P456frC1zXk118YiUojOAyhVW9LDbWbhjkFruuC+LVGaP3QR7rPfPJRUfG70alw8YfDAIzH2vKphpgr3g5L3Zi/n54Nxx/ZyQO4cWhBE7wkwEygVhZ29t2RNHJGC5g+dW81Xcp6zh91oRuGXyHnEa3DDdEYl6dHIoFMIhIno3J+uyaYW5gO9bHqKBQ/EObHXU9CfFiXp41IJhVPvulm8Oo5vSmH2IXwKbPHo/Mbmu4JmwjFBp6ammkCg1uJWTJzoPQXJIZTKNVY3ZRiC+x6O0XVRp7ej8MfFJenjult0ejs0Vm4uzKSCbV6Dq6bgpI13BQ0EagkEaqmpqR71OMcEHfsTN3s760ErrLUvcxcmmVaRkhWkJRtTOPWuB39Nkorawpo3LVPSi0HfKoJeDjgo6LS1tRMcZrM5OxFXGiEeRpm9fuyD80Q3BCs2xWg64iIadszbRSpCgzdUzn/kYuXDUVStcB7v230GPUmTBhv3cgCCwVaCwVba2tpHPU4+RQqKwOzYhzJ7Lej5ryo7nMVrYiTjChccvBWwVDyaP3ShKLDwocJIKpiqKmz06OyMJKxuSt7IoFNorp2F6LV0UoHFSqelWHB/gpP73AhTJg9I9pdKKjTsd7N4VRxfmfVJBY97ddqTJmcT1rclX2TQKUDmpb0olavAKLW0HSs2xegOyuQByVm6W3SudGks32htUsFMVWGdR+fXg8XTywEZdArT9QsweBll7ibLmlBRk2RWVYpGmTwgOU66tzOnLsmMudYlFXzKZ9CcMLmYLJ5eDsigU7DMS3vShUDdU/N+bUUVrNgc5fwxF5F++RaRnGfgmkZLg5FOoVbyn1RQqSk85NaKai7nFvlEKVT9bdDfhlK1Oe+Xnn9fAt2A5o9k8oDkXE2H3XhLBbXL8v/gf8Jn0BBP0V5kvRyQQaegmZf2oMxcCd6Zebumy2uyZE2M0++7SSVk8oDkXImYwtlDLpZtiGO48tfbqdZVVro03irCXg7IoFPYBrug7xzqvEfydsll6+IM9KlcOieTByTna200iA4qLF6dv6SCJ30GH8ZSdKcKZ61QPsmgU+DMS3th2mIomZPza5VNT1GzPEHDfg9y2wKpGAih0Piem7r7E3mpy7bQUFlgqEXbywEHBR27114b0dAVxJVTqNVbc36p+o3p+mrXLsttC6TiEWrX6b2ksXxD7ns7T/kMPogmuWY6r5cTCMjaa44hOt6FshqYEsjZNfzVSWbNS3H6fZkiLRWfxvfczFmQZHpl7krRrHBpzNZVflfEvRxwUNBxQu21EcVuIHo+zF1vRxHUb4zRctKQKdJSURq4ptF22mDFpijkoOymAnzaZ7BvKMmA8zo5gKy95jii6z3wzoLpS7L+2jVLE3hLTM4dlb0cqXidPeSibIbJ3IXZXzC62q1RrirsGSruXg7IoGMfiUFE90HUeVvJ5iS/pguWrovTdMRNIiaTB6TiFRtKV6FetiGGomavO6KTXpfzzlCCIYf2csZDBh0bEd0HQfei+B/I2mvOvz+OmYJgo5G115Qku7pw3IWuQ6A+ez2STV4dBdg/VFj7+FhFBh07ScURnQdQqh4GdfLraFwek0UPxTlz0C23oJYk0lWozx52sWRNHD0LC0Y9CjzmNXgrkkCGnDQZdGxG9HwEZhKlcu2kX2vx6jiDN1Q6muVCUEm6pe2MQTwGCx6Y/J4727wG/abgSMzeG7Rlkww6diNSiEt70/vt6N4Jv4yvzCSwMkHj+27kQlBJ+j1hKpz5wM3CB+K4fRNfMDpFVXjEq/NmkWxDnSkZdGxIXG2EWB/KnI0Tfo2l62Nc7dTovSR7OZL0cV0XdfqvqixZM/Hezqd8BpeSJo0234Y622TQsSmzfTdK5Wpwl4/73LLpKeYtSnL6A5kiLUnDU2j8wE3t8omVx6nQFNa5Nd4osg3aMiGDjl3daLm59cH4i4EuWx+nq0XnekiWu5GkkVzt1Ont0Fi6dvzlcXbc3LogWIRbF4xFBh0bM9t3o8ysB19lxudMq0hRGUhy9qDcK0eSxnLmAzfzFicpm575EFmdrrLMpfHrIi93MxLHBB3HFvwcTeQy4kojas32jE9Ztj7GpSadgT7Zy5GksVzv1ehq0Vm2PvO5nadL0kU9Q0W2dUEgIAt+FgXRsRfKqqF8/pjHzqpKMnNuiqbDci5HkjJ19pCL2YEk0yrG7u3c59KYo6v8VvZyRuSYoOPogp+jid1AXD56s7czWuqzYNmGGMFGg8iAY/7skpRzA9c02pt0lkBWU5EAACAASURBVK0ffW5HA3aUGOyOJBxb1HM0suBnERGdB8AoQ5l134jHVNSmmDLD5NxROZcjSePVdMTNzLkpZswdua7ARo+OC9gjy92MSgYdJ0hFEZ37UeY9AupwNdQES9fFCDYYxCLyTy5J4xXpV2k7Y7BsXZzhtj7wKvBJn8HOSILJ1zFwNvkEcgjR8yGYCZTZ6+/52ey6JGXTTJo/kr0cSZqoc0ddTKtMMbPq3rmdx7wG103BYVnuZkwy6DiFMDHbd6HM2QBG6Z0/YMnaOC0nXcSH5J9bkiZqKKzS2nhvb2e6qvCwV+dXg3FZ7iYD8inkJNeaINKDUrXl9rfmzE9SUm5y/pjcukCSJqv5QxdT/Sn81b/v0TxVYnA+YdKUkAtBMyGDjsOYbW+j+O8Hn5/0XE6ciydcxKPyTy1JkxUdVAk2GCxdFwMEtbrK/S6NXw7KmZxMySeR04Q7EVdPo9Y8xtyFSbylJheOy7kcScqW5o9cTJlhUlGd4rMlBoeiSbqLbCHoZMig40CifTeUzWXxOrh40iW3oZakLIpF0nM722uhUlN5Sy4EHRcZdJwo3k+l8ktKykwunpTVByQp24LHXKxsLeF9V6woF4JOhmM2U6mpqSaVTBZnVYJ7CBbXHaelex2JKW4YOmp1gyTJUdabLpIJuLE+DL/wIjdC5HbtS00fPazIno4DVQZSlJanuHDERKl6GDSP1U2SJMcoU+BRr8EbgwmmzBl+3Y40MscEnaKtvXYPwZI1MVoaDGKdpyHej1K12epGSZJjPFnioitlcrjfpO20MandRZ1E1l4rUhU1Kcqmm1w45gIEZtvvUCpWg2eG1U2TJNur0hTWuDVev7kjaPNHLmbMTjFjjqy3likZdBxFsHhNjNZGg9it6gP9bdDXjFrzmLVNkyQH+Fypiw9jKdpv7gg6NKDS3mSweLXs7WRKBh0HmTk3xTS/yfljd6/LMdvfgfIATF1gUcskyf4ecGlU6SpvfixFuvlDF/55Kab65dxOJmTQcZDFq+O0nTWIDn7szxq7jug+hFrzKCjyTy5J42WQ3ivn7UiCfvPuHOnBGyod53UWr5K9nUxk9ATy+XzseOoJ1q5ZzbatW/ji88+N6yLffeHbbNu6BZ/Ph8/nY+2a1Xz1j7+Mz+fL6nWK2bSKFLOqUpwfoZK06HwPNE96fkeSpHHZ6k2nAe8dYa+c5g9d6Wru02VvZywZrdP586//KS/+6CVCoV4A/P5ZfOOv/oIXvvO9jC5SVxfga3UBvvbVrwBwuaeHF77zPSKRSFavU8wWrYrT0awzeGOEzxFmHNG+G6X2k4grDZCMDH+cJEl3maYqPOoz+OeBOCOlC/Rf1bjcqrFoVZyP3vbmtX12M2bQ2bZ1C6Ge0O1AANz+77VrVnP4yNgLD19+5accOnyUCr+fcDhMsLUtJ9cpVmXTU8yuS7LnFd+ox4krp1AqVqHM24oI/jpPrZMke9tRYhBMmJyKj96LOfehm82fj3D2kEmkXw5jj2TMO7Nu3Rpahln/0tDQyPbtj2R8oVCol4bG08MGnGxepxgtWhXnclCn/6o25rFm62/TVahLZuehZZJkbwtuVpF+PYMq0n2XNa52aix6SM7tjGbMoLNyRf1dvY9benp6qQsEstaQfF3HaXxlJlWLkpz7MMNK0oNdiN6TqLWfzG3DJMnmVNIp0u+No4r0uQ9dVC9N4PbJvXVGMmrQuTXRHw6H7/lZJBKhtLRkXBcL1Nawon55jq6jTPArG69h3deCB+Nc69bou6xnfI64tBe8fpSZK7PUDnvfw8L4kvew0O7jBo9OuarwViSZ8Tm9l3QGrqnMvy9RAPfCyns4slHndDJ52Pt8vnsSAj5ufl0da9esJtjaSrC1jS8+/xw9PSF279mXtet8+1t/jabe/cvuP3CQ/e8dHPO1AcrKp2V0XCEx3Clqlw9wat8sysrHN3mZuHqURO2jeFM9KGZ2SrPb8R4WGnkPs2Oy99GL4Ek9xtspHX2Kl7JxnHvprJsl66/Sdb6CVMK+czuZ3MPNG9ezedP6u76XMgXN7f0jnpOXKtO7du+lofH07f//yquv8cMffP/2PE82fPNbf0siPrGx1LLyaQzc6MtKO/Jpydoo4esKrWeGgOj4Tu5/F7VsMYNlKxHtuybdFrvew0Ii72F2ZOM+PlFi0IvK3hsDjHfngvApQeB+mDkvxIXj9txaJNN7uHPnW+zc+dZd3zNcLr70lT8b8Zy8hOHhAktjw+ksr8MRE/zKxmvk/0vTTeavjN+xLmecryFS6aSCyjXgnTHJ9tjzHhbWl7yHhXIf5+kK6zwaPx+MIyZwvhBw4ZiLBQ/EUTWzAO5J/u/haEYNOuHw4KgnA2MOrY2kJdhKXV0g59dxqtrlCZIJhY7zk+is9gcRfU0yqUCSblKAZ0sMjsRStCUnngzQdsZAUWHeYlkI9ONGDTqRSIRweBC/33/Pz/z+WRkFix/+4PusXTP6KvhsXKeYKKpgwQNxzh93IczRJ+3GItregdIqmL40S62TJPta49bwaypvZpAiPRozpXDxhMHCB+OM9cm/2Iw5vNYSDA470V9R4acxg/mYyGCEYGvrsOe3tASzdp1iMndBEt0laDtjTP7F4v2IzgPpKtRqFl5PkmzKq8COEhdvRRKEsxAngg0uvKUmlQFZGudOYwadQ4eOsHJF/T3fX1G/nA8OHh7zAgcPHR52/c3atavZtXtv1q5TPAQLHowTbHCRSkyul3P7FbsPgZmQm71JRe1Jn8ENU/BeNDtDYomYQusZgwUPyMWidxoz6Ozesw+/fxZ+/6zb3wvU1hAeHLynNM0Pf/B9vvvCt+/6XkPDabZt3XLX93Y89QShUO/tlOnxXqeYzaxKMWWGycWTWeyViBRm8DcolevAOzN7rytJNjFPV9ng0XktHCebyzovnnAxc47c9uBOGc1C/80Lf8czTz9JT08ISA95DVeEMzIYITx49/zLrbI3tzLVSnw+ekKhYc/P9DrFbOEDcS41GcQiWU48vJ1U8ATm2X/O7mtLUgFTgOduJg+0TiJ5YDiRfpXOCzoLH4xz9LeyEChkGHQikQivvPramMf95Te+Oez3g61tI9Zcm8h1ilXZ9BQVtSl2v5yb3H/R+jbK/f8JZUY94mpjTq4hSYVmg0dnpqbyj/1DOXn9C8ddPPxshNNlJpEB+y4WzRZ5B2xkwQMJelo1Bq6NXdhzQhIDiI53UWoeA82ei9okaTxKlfRczpuRBIM5SjLr69G42q0x/345twMy6NiG22syb0mC88czLOw5QeLyEUiEUeY9ktPrSFIheLrERShlcjBLyQMjOX/MRe3yBIZbpk87JujU1FQTCNRa3IrcCaxMMHBN5UpHjno5twgTM/gWSsUqKJmT22tJkoUWGCqr3Bo/C8dzvpLmclAjOqhQuzw7dQ4LUSBQSyBQS01N9ajHOSboOJmqCepWJLhw3MVYFVyzItyBCB1HDTyRn+tJUp5pwHMlLt6NJunKcNuCyVG4cNJF3X1xFKW4ezuOCTptbe0Eh9kEzgnmLU4iBHROpuTNOIn23eAuT/d4JMlhtnp13Ar8JpK/nkf7WQPdEMxZ4MzSOMFgK8FgK21t7aMe55ig41yC+ffHaTllYKby2OtIRRFtb6PM2wpGaf6uK0k5NkNVeNxn8L8HE8Ty2OlIJRRaTxtFn1Agg06BmzUvRelUk2Bj/kvUiCsNMNglC4JKjvJcqYvmhMmpeP4XbLacdDGtwmRaRfEuFpVBp8AtuD/OpXMG8SFr/lRmcCdMWwRTF1lyfUnKpofcGnWGyr+GreltDIVVui7oRd3bkUGngJVONamoTXHhuIWFOKPXEJ3voQY+JQuCSrbmU+AzJS52DiboM62bzL9wwsXchUm8pdmtfmAXMugUsLr74vReyuFi0AyJrvfBjMu1O5Kt7ShxccMU7M/xmpyx9F3WuN6jUrfSuenTo5FBp0DpLkHN0gQXT+R2MWhGRAqzZWd6l9GS2Va3RpLGrU5XWevW+OlALKsFPSfq4kkXtcvjqFrxpU/LoFOgapYliA4pXG6ztpdz20A7InQCte5JUOTbRrIPHXi+1MX+aJJLeVmTM7bOCzqplMK8xcXX25FPj4IkqFsZp+WkC0ThLM4U7bvAKEWZvc7qpkhSxh7zGegK7BwsnAe8MBWCDQbz70tQbDuLyqBTgCpqU3h8gvZs7AyaTakoZutvUKoeBs90q1sjSWOarSls9+r8LByn0PLFWhsNSqebzJhTXOnTjgk6Tqq9Nv++OO1NBol44fRybrvWBNcvoAaetLolkjQqhfSw2kexFE2JQpjJuVtsSKWzWWf+/YXTA5uMQEDWXrOl0mkp/NUpWrK5M2iWmcHfQEklyqwHrG6KJI3oYY/ODE3lF4OF1sf5vYsnXcypK6706fwV88qxtrZ2EvHCfXNlqm5lgtAljYG+AkkgGE4ijGjfhVLzKOL6eUgMjn2OJOXRDFXh0yUGrw7EiRTwlMn1kEbfzfTp0x/Yew+rW7UvDZeLjaMcJ3s6BUR3CaqXJtIJBAVOhI7BYDdq4NNWN0WSPkbwb0oNmuMpjllQ6ma8Lp50UbM8UTTp0zLoFJB5SxLEowqXWwu4l3MHs+XXUF6HMmO51U2RpNtWKSmqdZXXCihbbTRdF3WEgKqFzqw+/XEy6BQMQd3KBMEGo6DSpEcV60Nc2oNS+zhC81jdGkliqgqf0pL8cjDBDQtL3YyHmVJobTQIrLT/9EAmZNApEDOrUpSUm7QVWpr0GMTlIxDtI16x2eqmSBLPlbroECoHY4U/rHan1kaDaX6TqX57tXsiZNApEHUrE3Q265ZVk544gdnyJqmyOpi22OrGSEVstVtjoaHyy5T98qOGwirdLTp1RdDbsdsTzpE8pSaz65K0nCr8BIJhDV3BuHIkvb217rW6NVIRmqLA50pcvDGYoM+mj7WWUwZVi5K4PPYYFpwoe/51HKZ2eYIbvSp9Pfb9c+hXj0P8Borc8E2ywBdKXXQmTd6L2nd4qrdDI9KvUrPMHgkQE2Xfp5xDKKogUJ+42cuxSQLBMBQE5sU3UKYvlcNsUl6tcmsscmm8HI7bvIqZQsspg8CKOE6uxyaDjsVm1yVRNUHnefuNQ99j6AqiY1967Y4cZpPy4M5htWs2yVYbTXuTgdsn8NfYt8c2FscEHbvWXgusSNB+1iCVtG8v506i6+DNYbZPWd0UqQj8mzL3zWE1Z6xxScYVOs4ZBOrtN8QWCMjaawWvdKrJrKoUwUZ7pUmPTmBe+CXK9MUoM5ZZ3RjJwda5NRboKq/YfljtbsEGg9kB59Zjc8CYTpoda6/V1se50qERLuQ6axMRvYpo34NS+wSiv03WZpOybrqq8NkSF68Pxh0xrHan670a13tVapYnaDpsn3pssvZagVM1Qc2yhMN6Ob8nLh+GoRBq3VNWN0VyGAX4YqmLC0mTQzZbBJqpYINB7fIEiuqsgAoy6Fhm7sIkZkqh66JjOpv3MC++AVNqUGbdb3VTJAfZ5NGZq6v8dCBmdVNypqPZQDMElbXOmKu6kww6FgmsiNN2xkCYzkggGFbsOqL1dyi1j4N7qtWtkRygQlPYUWLwWjhOv/M6AbelkgqXzhoEVtgvoWAsMuhYYMqMFNMrTVodOrR2J9F7Am60os5/GjuvQ5KspwF/WObmZCzFcRtsWTBZwUYDf3W6JqOTyKBjgdr6BD1tGpGB4rj9Zsub4J2JMnu91U2RbOxxn0GpAj8v4J1As2ngmsbVLo2a5c7q7RTHU6+AaLqgekmiKHo5tyUjmC1voszbAr4Kq1sj2VCtrvKoV+flcJwhBw+rfVzraYOaZc5KKJBBJ8/mLkySSipcbnVuAsGw+poRV06hLvgMKEX2u0uT4lbgD8pc7I8maU44a6hpLJ3ndVRVMDvgnIQCGXTyrLa+CBIIRiBafweKhlK9zeqmSDbyuRIXCQFv2mQn0GwyUwqXmgxqbVihYCQy6ORR2fSbCQSni2ho7U5mAvPiL1EqVsHUBVa3RrKB+10aD7k1/mkghnM+649P6+l0QoFvijN6eTLo5FFtfYJQe7p8edEKdyI696PW7QDdZ3VrpAI2VVX4Qmm6mGd3yjlzGuPVf1Xj2mXnbHngmKdfoRf8VLUiTCAYgeh8D6LXUOfvsLopUoFSgC+VumhPmux3SDHPyWhtdBV8QkEgIAt+FpS5C5OYpkJ3UE6ig8C8+Asoq0apWG11Y6QCtM2brjrw8kDMUcU8J6rzvI7ukAoFjnkCFnrBz9rlCdqLNIFgWLEbiOCvUeY/gxhog0jI6hZJBaJGV3nCZ/Cj/pijqw6MRyqpcOlcuh5bd0thjpbIgp8FpHSqycy5KdrOFOabxSri6hnElQbUBZ8D1TGff6RJ8Cjw78pcvBdNcqbI0qPH0nraoKImhafE3vdFBp08qFmW4EqnRvi6vN0fJ1p/C4qCUvO41U2RCsBzJS6iAt4owvTosdzo1bhxVaV6qb3vjXwK5piiCKqXJmgr1jTpsZgJzPOvo8y6D6Yvtbo1koXWuDVWuDV+UsTp0WNpO23cLItj33HHjMY0fD4f27dtoaenl9LSEioq/Lzy6msZXyRQW8OKFcspLS0lEKgl1BPi5VdfIxKJ3HXcd1/4Nrt27+XgoSMArKhfzsoVy4c91i4qalNohqDzghw+GlHkMqLtHdS6pzAHuyF23eoWSXlWqSk8W+riX8NxQkWcHj2WS+cM6jfGmFmV4kqHPZ8pGbX6z7/+p7z4o5cIhXoB8Ptn8Y2/+gte+M73xjw3UFtDXV2AN9586/b3djz1BC/9+B/5P/7z/3X7NQHq6gJ8rS7A1776FQAu9/Twwne+Z9uAA1C7PE7HOYNUUiYQjEb0HEUpD6Au/Bzm6ZdA2HvcWsqcAXz5ZvXoIw7dlC1bEjGFrgs6tcsSzg0627ZuIdQTuis43PrvtWtWc/jI0VHPX7Fi+V0BB+CNN99i/bq1fPVP/uiuwPXyKz/l0OGjVPj9hMNhgq1t4/ldCo7bZ1JRm+Ld1+yz5ayVzItvoK78Gkr1dkTb21Y3R8qTz5a60IDXwoWbfVpIWs8YbNgxhPGuIBGz34fZMed01q1bQ8vNVLg7NTQ0sn37I2NeYNu2R1hRv/ze8xtPs3JF/T3fD4V6aWg8bfuAA1CzNMHANZXrITl1lpFUND2/U7Eapi2yujVSHjzo0ljj1nhpIIYMOZm50qExNKgwb7E9EwrGfBquXFF/Vy/nlp6eXuoCgYwu4vfPGn/LbE9Qs+xWAoH9Po1YJtyBuLQnvembu9zq1kg55NcUni9z8fpggi45jzMOyh0JBfYzatDx+dK1scLh8D0/i0QilJaWjHmB//z1/5vde/bd8/0Kv5+WluCw5wRqa4btHdnJjNkpfFMEl87Zc9zVSqL7IPS3oy78PCia1c2RcsAAvlLmpiGW4n1Z5mbcLjUZlM80KZ9pvzmwUYNOJkHlVmAaD5/Px7p1a/jlr3591/fn19Wxds1qBiMRGhpP88Xnn2Pb1i0Zvqoywa9svMa9X9XLknQHdeJRLauvW5hf2b+H5sU3QPehVD9aAL+fPe9hIX89W+pCBX4aTsj7OIGvobBGqF2jelkyB6+fjXs4Mks+hn/p+ec4dOjIPUkIu3bvpaHx9O3//8qrr/HDH3z/9jzPaL79rb9GU+/+ZfcfOMj+9w5m1Kay8mkZtn5sqm5StWiAhndnUlZePJWUs3kPAVLd7xCr+Ryu1FX0gYtZfe1Cle17WIgeVJI8qCX5h6QLd7mPXKTZFMN9DLW6WLz2Gm0NU3NSXiuTe7h543o2b7p7G/qUKWhu7x/xnLwHnRX1y6mrC/CX3/jmPT8bLrA0NqR7PMMdf6dvfutvJ1x7rax8GgM3+iZ07nCql8RJxKDtbBQhYll73UKW7XsIwI0+FN4mVr2VoStBiF7N7usXmJzcwwIzV1N4aqqb18IJLsaGcnKNYriPAIONgsVrTUqn99J1MbuLzzO9hzt3vsXOnXdnJxsuF1/6yp+NeM6oQSccHhzzouNZQ+Pz+dix49P8zQt/l/E5LcFWtm/fmsGRgomt0r3zE0J2JjNrlie41GQgRPZes7Bl/x7efrWeo1BWhbro85iNPwbTnpOnY8vdPSwUXgX+eIqbo9EkR2K5msdx/n28xUxBR7NBzbI4XRez2X+Y7D0c/ZxR53QikQjh8CB+v/+en/n9szIKSnf691/9Ci/+6KVhA9UPf/B91q6xf5n7kvKbxT3PygSCbBEt6bk/pe4pi1siTZQC/EGZm0FT8L9lXbWsaTuTLgLq9tlnMfWYKdMtweCwCQUVFX4ax5hnudMXn3+Ol1/92V3p14Hamtv/HRmMEGxtHfY6I2W5FaLqpQmudauE+2TWVdaYCczm11CmLkSpXGN1a6QJeNSrE9BV/udAXNZVy6LrIZWBPpXqJfYJ5GMGnUOHjgy7iHNF/XI+OHg4o4ts27qFgwcP37PeZ8WK36dFHzx0788B1q5dza7dezO6jvVuFveUWxhkX/Qa5sVfpbPZykbfmVAqLEsNlU/5DP5pIEaf6ewhr/xTaDtjULMsiV2GE8cMOrv37MPvn3XXAs9AbQ3hwcF7ss9++IPv890Xvn3X91bUL6eiwn/7vFtrcNauWU3FHcN2DQ2n70mP3vHUE4RCvcOu8ylEs6pSuL2CjvMy6OREXxOi+2B6/Y5RZnVrpAzMUBX+sMzNW5EETXJ/nJy41KRTMtVkWoU97m9GEw9/88Lf8czTT9LTk97dsaLCP2yxz8hghPDg3fM8f/2N/weAp3c8ec/xu3btuf3ft8refPH55wAo8fnoCYUyKipaKKqXJui6qJOMZz99UUoTl/ailFSiLnoW88w/gbDf4rhi4QL+ZIqb84kU7wzJQbVciQ2phNo0qpcm6Osp/GH9jIJOJBLJaCuD4dKav/D8H2bcmGBrm21rrumGYM6CJIff8lrdFIcTmOdfR13xJyi1n0QEd1rdIGkEz5elC3m+LAt55lzbWYMHtkZpOODGTBX2h15ZiTJL5ixIkogphNoL/5OG7aWi6cSCmStQ/A9Z3RppGFu9OssMjRf7Y8TsMdVga5eD6f5DZW3h9yhl0MmS6qWJdJ01UdifMhwjEkonFtR+Espqxjxcyp+lhspTNxMHemXiQF6YKYWOZsMWW1k7JujU1FQTCNRacm1fmcmsqhTtZ2UCQV5dO4voeh910efBJStSFwK/pvDvytzsjCQ4IxMH8qq9Kb1mx+W15r4HArUEArXU1IyeXeqYoGOleUsS9PWoDFyTQ2v5Jjr2wUAH6uIvgCqDvpU8Cny1zM2ZeIpdMnEg7/ouqwz2K8xbXNj33jFBp62tneAwm83lnqB6SUL2cixkXvwFKGp6Dx7JEirwR2VuYkLwqkwcsIhC+1nrhtiCwVaCwVba2tpHPc4xQccq0ytNfFMEHc2y7I1lUnHMcz+FKbUoVVusbk1ReqbEYI6m8OJAnMKfVXCuW/vsTJlRuEsJZNCZpHlLEvS06cSj8lZaKnY9ndE25xMoM+6toCHlzgaPzgaPzosDcW7IxAFLDYVVrnRozCvgsjjySTkJiiqoWpSgvUn2cgrCQDsiuBNl/g4onWt1a4rCQkPl8yUGLw/EaU/KxIFC0N5kpOd1lML8ACCDziRU1iZB+X2OvGQ90XsCcfko6qIvgFtmtOWSX1P44zI37wwlOR4v3OGcYtN1UcflEcyqKsy/iQw6kzBvSZLO80bBrwAuNqJ9F4Q7UBc/D1ou9qWUShX4D1PSmWq/iRTuUE4xSsYVui7qBVt5WgadCTLcgspAkktyaK0ACcwLvwAzibrwWVDk2zybdNI11W6YgldkplpButRkMHt+Ek0vvCE2+a9xguYuTBANK1ztkmtzCpKZSGe0eWegBJ6wujWOoQBfLHVRpir8qD8m98YpUKF2jVRSYXZd4f2FZNCZoOolSS6dM7h7a1epoCTCmE2vokxfhjJno9WtcYQnfQZLXRr/eCPGYOF9iJZuEkKh45xekFlsMuhMgG+KyYw5KS41yQWhBW8ohNn8ryhVD6PMXGF1a2xto0fnYa/O/+iXNdXs4NI5g4rqwtvK2jETEjU11aSSybxUJZi3OMG1yyrh6zJm20J/ENHyJkrdDkQ8DP322f68UNS7ND5XYvCTgThBmRptC9dDKgPXVaoWJbl4wpXz6wUCtQBo+uhhRT41x00wb/GtoTXJLsSVU4jO/aiLngNfhdXNsZVaXeXLZS5+OZjgpEyNthGFS0068xYX1hCbY3o6bW3tJOK5z6SZOsukdKpJ53nH3LqiIToPgKsMdckXMU+/BLHrVjep4Pk1hX8/xc2BoSTvRgtvUloaXUezwfINcUqnmjkfmbk1ymS4XIw2gyp7OuNUtThB6JJGLCJvnR2J4G/SVamXfAl0n9XNKWjlqsJ/urkW5w25FseWIv0qV7tVqgqotyOfnOOhCKoWyaE1exOYF16HRBh1yfNyO4QReBX4j1Pc9KQEL4fjyLQB+7p0zrg5xFYYf0UZdMZh1twULo+gu0UOrdmaSKXX8Ch6ulyOItda3ckF/PspbhJC8OP+GDJtwN46z+v4pgimVRTGX1IGnXGYtzhBd4tOMi7X5theKobZ9DJ4pqIu+CxyvVWaBnxlihuvovAP/TFkvQH7iw+phNq0gkkokEEnQ6ommLNADq05SiKMefZfoKwKpe5Jq1tjOQX4t2UuKjSFv++PESmM0RgpCy6dM5i7KIlSAJWnZdDJUGVtEiGgp00OxThK7Drm2X9Bmb4EpfpRq1tjqWdLDBYaGv/9Rkzui+Mw3UEd3RDMmmd9yrsMOhmatzhdUVqYchjGcYZ6MZteQal4EKXqYatbY4lnSgwedOv8/Y0oV2TAcZxUQqH7YmGs2ZFBJwO6S1BRm+SS3JLaucKdmE0/RZm9AWX2Bqtbqa/M3wAAG41JREFUk1dP+Aw2uHX+vj9GV0oGHKe61JyuPK1q1v6NZdDJwJz5SeJRhaudcmjN0Qba0ltez9uCUrHK6tbkxXavzlavzj/2x+TOnw4XatcQKagMWLvIVwadDFQtStDRrCMznIrAjYuY53+OUvM4iv9Bq1uTU494dD7lM/hRf4wWGXAcT5gKnRcMqhZZG3QcM16Uq4KfLq/JrHkpzhyUO1AWjb5mzAuvp1OphYnoPWF1i7Jus0fnqZJ0wGlKyIBTLDrO62zYMYTuEllf+hEI1AKy4OekzV2YJNKvcD0kb1VRuXYWcfEXKIFPo8xcaXVrsmqTR+eZEoP/ORDnjAw4ReVKp0Y8qjDHws3dHNPTyVXBz3mLknTIzdqKkrh6BlBR5j+d/v9XTlnboCzY6NH5zM0tChplxejiIxQ6z+vMXZSgPcv7gcmCn1ngLU1v1tYhK0oXLXG1EXHxVyh1T6HMut/q5kzKZo/OZ0sMXhqIc0oGnKLV0Wzgr07h8lrTy5VP01FULUpwo1dl4JrMWitm4mojYKLM/wwoCiJ03OomjdsWj86Om0NqsodT3Pp6VCL9CnMXJAk25H5zt4+TQWcUVYuTN7PWpGInrp5BCJFOLlA0RM+HVjcpY9u8Ok/czFKTczgSKHQ0p7PYrAg6cnhtBKVTTabOMulolrXWpJuuncU8/68oNY/ZZgHpEz6DT/kMXpQBR7pDR7POzLkpPKX5f0/IoDOCuYsSXLusEhmQt0i6Q18z5rmfoVQ9jFK1xerWjOozJQZbPDp/f0OmRUt3G7im0X9VZe6C/GexySfqCKoWJmUvRxrejYuYTS+jVK5FqXnM6tbcQwG+UOpijVvnv8mFn9IIOs7rVC3Kfy02GXSGUTY9xZQZJp0X5HyONIKBdsyz/wtl5sqbKdWFkVKvA18uc1Hv0vj/bkRlaRtpRJ3NBtMrTXxl+X2PyKAzjKpFSa50akTD8vZIoxjswjz9E5QpAdRFz4Fi7YcUF/C1KW7m6So/uB6VxTulUYWvq1zvVZm7ML+9HflUvYegamGCTrk2R8pE9Arm6ZfAOwN16ZdAs6ZcUqkCf1rupkxV+MH1KFfl9gRSBtILRfM7r+OYoFNTU3279s9klM8yKSkXcmhNylz8Bubpn4Cqoy7/Mhhleb38dFXh61M9JIH/90aUfhlvpAx1NBtM85uUlE9+iC0QqCUQqKWmpnrU4xwTdLKlamF6aC0WkbdGGodkBPPMP0OsH7X+K+CdlZfLVmkK/+dUD11Jk7+/EWNIBhxpHCL9Kn09+R1ic8yTta2tPQsVpgVzFyZk2RtpYswEZvPPEDeC6R5PWU1OL7fUUPmzcg8nYkleGohjbcF6ya5uLRSdrGCwlWCwlba29lGPc0zQyYZpFSbeMkGXHFqTJkqYiJY3EJePoC79EsrMFTm5zAa3xlenuPndUIKfDyaQHRxpojrP65TPNCmblp/ySPLpeoe5CxL0dmjEozIWS5MjOt6FaB9K3Q5wT0V0HsjK6yrAUz6DzV6dfx6Ic0LWUZMmaSiscrVbZc7CJOeO5L7OpHy63iaYszAps9akrBFXTqUXkc5ejzL/GVAm9w/aBfxRmYt1Hp3/fiMmA46UNZ3nDaoW5meAVgadm6ZVmHhLBd0XZdCRsqi/FbPxxyilc1GX/TswSib0MlNVhf881UOFpvJfrkcJykWfUhZ1XdCZMiM/Q2wZPWF9Ph/bt22hp6eX0tISKir8vPLqaxlfJNPzJ3udyZBDa1LORK9iNv4YddGzqPV/gnnupxDpyfj0Gl3lT6a46Uia/NNAjKicwJGy7O4httx+8M7oCfvnX/9TDh0+yuEjR9m9Zx+7du/lG3/1FxlfJNPzJ3udiRPMWZigSw6tSbmSimI2vYzoa0Zd/kcoM5ZldNpat8aflbv5KJbkf/TLgCPlTr6G2MYMOtu2biHUEyIU6r39vVv/vXbN6jEvkOn5k73OZEyZGcdbKuiSQ2tSLgkT0fobRNvbKPOfQanezkg121Tg8yUGz5a6eDUc55cyQ03KsXwNsY0ZdNatW0PLMOtfGhoa2b79kTEvkOn5k73OZPhr/v/27jw4jupO4Pj39UzPjGYkS8LyyBjb0sictmQSQ7BNHAI2UGAweANJNuySg1TYJZUNbG2FHGQroRKHAKkKTqWWTRGybEggoXYhEMwGYnMGLNvEgGUwPpB8IVuSbcnS3D3db/9o2dix0TWXZvT7VKnskuZNPz319O8dv34dl6k1UTC6eyPOO79G1bVgnH0DeCuO+/mkwS1tmn0eVvYleT0lCQMi/46dYsunYa+yc1uajxt9HNHV1UNTJDLsAUZaPtvjjJ2mvjFG53Z5jIEooOhenLZfgseH0XIzhKYBMMtr8I3aCmwNP+lLskc27RQF5E6x5Xd3giHnk4LBIADRaPSEn8XjcSorh87EGWn5bI/jUoxle/npYYfJKTU4tTY+tqcvbdKGI2ZF3RHPzCUYsz/Pwvf/lyu9CVYnMjwTz6ClLbMk7TdanTtM5l6UIlgDA4dhbG04dJkhg85ILvbBYJB4PJ5V+WyPA3Dn97+Lxzj+l335lbW8/Je1Q77vgtoY5z8XpKvCZJ3fg5yoY1dVXVvsKpSkir4NfCrxFhFS/HfNfHY471JppotdrZIm5+LoaSBzykc52P8mldMrifcNX+aiRQu56BMLj/ue7Wi27e7/0DJls3L+ve//ECs9+g/q/xyGztoqliuLmU6KR6Np2TRxDKqqaxk43FvsapSciNfgC1Um/akEPzPD7AueCVXNONsfh9i+YlevJMm5OAbeCoxZ10LVdF55ag+V9oERteGqVc+watUzx33P9Pn4h5u+/uGHyrqy44Ye/BotRZvjZdtAlC9N8nF7jZ9fD6Tl5rtROXZ0KBF7JBRwWYWXK4MmLyYyPB3PEKyO4mz+FarhMow5X0TveR69r7XYVS0xci6O2qRGjNOXQ+owzqZfQHoAjo4Ux9KGQ5cZMuhEo7Fh336oKa9sy4/ldWN1wNH8tC/JspDJ16v9PJfI8GzcQkKPyLVaQ3FjlY+wx+CB/hTvWA5HL5badtOqD3dgNF2Nqj4d570/gHXieqcQWVEGasZi1NT56H2vufsF6mPOxTwZMnstHo8TjcYIh8Mn/CwcnjJsUBlp+WyPkysZ4ImYxX/2p7jQ7+HWaj91hqzxiNw5z+/hmzUBUhp+3JsYDDgn0fuu2+tUYMz9Z6g9q7AVFeUtUIcx5ybU5Nk4Wx5G73lhMODk37Ap0+0dHSdd6K+vD7N589vDHmCk5bM9Ti5ttRx+3Jek39F8szbAokAZzUKKoggpd7POz1X6WBW3+EV/iuhwMxfWAM6W36A7/4JxxnWoWdcW7XHYonyoqQsw5t6MTh50OzYDQz//JteGDTqtreuZ29J8wvdbmufw2tp1wx5gpOWzPU6uxTQ8OJDmsWiaq4MmX53kp0ZGPWIMmn0evl1bQY2huLs3ySvJ0d18p/e14rQ9gKqYgjH3FqielaeairLmr8GY/QXUaYtwdjyB3vEE2KmCV2PYoLPm+RcJh6cQDn/w+N1IYwPRWIx16zcc99qV993LXSvuHFP50RynkDakbH7clwTgOzUBLgx4JalajEilgs9X+ripysdLCYuVh1P0OGNc3E704Lz9K3TPGxhn/T2qaRl4ArmtsChTCjV1vjtNm0nibLofDm0pWm1GNG/0gxV3s/zaq+nq6gbcKa8VP7rnhNfFY3GisRPXX0ZafqSvK7Q+R/Mf/SkW+D0sD/k4z+fh0WiaA2O9gIiyN8/n4fpKHz22wz19SfbnYmcB7aD3voQ+9C7GrGtQ596C0/EM9G7N/r1FeQrUYcxaBoFT0O1/RB8s7FLFyYwo6MTj8RE9YuDbd3wvq/IjfV2xtKZs3rGSfDpk8q3aAM/FLdYkMsjOWOKIyYbi05U+ZpkGq2IWLyUzuU/cjXe5z+g5dSHGGddB3w6cnX+C9IffkCcmGOVFnbYINe3j6ENb0Ft/D5n8ZgCPlKyQj1K/o3lwIM1cn4frQybn+738PprmPbmvZ0LzAJdUeLkiaLLNcrirN8mhfI6EtYPufBV9aAtGZCnGubeg97yI3r8euT9lgpsUwYgsBWXgbPs99O0odo2OI0FnjDalbbZaNkuDJl+r9rMxZfNk3KJfptwmnLNNg+tCPvwKHh5I81YhHyOdPISz5Teoyc2ohstRUz6Cs/P/Cp6RJMYB3ySMhsuh9iz0vrXo918GpzCPoB6Nsgk6DQ0zsTMZOk7yeIR8SWn3vp51yQzXV/r4bk2AZxMWLyUyjL8/tci1UwzF34VM5vg8vJDI8FzCIlWkPoc+uBndtx01/ZMY59zoTqns/rN7d7kob8qLmrYQNW0RDOx206CTBwpejUikEQCPd+iwUjZBp5g6bc3PDqeY5/OwPGTy8YCXJ2NWYXu8omACCi6vMPlkhZftg/d0dY+HRxDYKfSu59Ddb2A0XoE692tuj7fzVXDyu129KA41uRk1czEAzo7HSyKppGyCzq5du8e04WcubUzbbE7bLK7w8o9VPi7OODwZs9gp6z1lwQNcGHD3S+t3NL/sT7Hlw3YUKKZED86Wh6H2TIyZl6HCH3XXe3reRNZ7ykTVTIyGyyBQh+58Bb1vHejidnKPzDKZPh+Lhnhd2QSd8SIN/CmR4bWUzVVBk9uq/bSlbVbFrdykzYqCU8A8v4ergiY+pVgVt1ibzIz/ffl6t+H07UDVn4+auQQ1bSHOnufh0LvFrpkYq2A9xozFUN2E7t6IfvfRcZOVNlISdPKk39E8Gk3zQkJxVdDkmzUBXk/ZPBu35P6eEtLi83Bl0GSyoViTsHgxkaGknnSjHfT+9eieN90U61nLYVoPzt6Xxl1WkxhCYDJq+kWoyXPQB99Bb7ofkoeKXasxkaCTZ/ttN8W6wWuwNGhyR22ADSmb5yT4jGstPg9XBE3CHsUriQyrExbxUv5z2Wn3xtL9G1CnLcI489MQ75bgM94dE2zo3e4+4jy+v9i1yooEnQLZlXG4vz9Fo9fgysHg80bKZnXColOm3cYFA3cabUmFSZ1H8XIiw/0Ja/iNOUtJJu4mG3S+ipr28cHg04PT+ergtFs5/bIlLDQNY9qFcMrZbrDZ/GDZPNRPgk6B7RwMPjO9BpdWeLm9JsAWy+GFhMW28bgoPQH4Fcz3e7mkwktAKV5OZng5YREr5+uvFfsg+Eydj9G0DGYsRu9rRR94a1ze3zEhVM/CmLYQJjW602hlMLL5WxJ0imR3xuFXA2nCHsWSCpN/muSn29a8lLB4PWXLfT4FUGsoLgp4WRjwktCalxIZXkuW2JpNtqyY+4TSzldR4fNQp30CNeMSd5G6a4Pc51MIhhdVNxd16nzw1aAPvIVuXwWp8nzktgSdIuu23YSDP8ZgUcDL1SEf14agNZnh1WRG1n1yTAHnmAaLKkxmmwY7Mw6/i6bZlLbHfzZaPtkp9+mR+1tRp5yDOnUBatqF0LsVp+uvcLi92DUsP4E6VP15qCnngpNG79+A7t4ImUSxa5ZXEnTGiah2U61XJzLM9XlYVOHlkooA2yyHdakMm1I2cnvf2J1iKC7we5gf8FJpKF5PZvhJLM1eWU87nnbQB992dyMOTUPVn49x1mch1Y/uecudepPRz9gZphvUwx+Bqgbo78Bpfwp6txXsyZ3FJkFnnMng3mS6MW0z1aNYGPDyqZCPz4RgYyrDX1M272UcWe4dgYCCc30ezvd7OcM02JNxWJPI8HoqQ1IacHixTnT7U+hdz6GmzEVNORc142I43I4+0IY+tBWcCTUZOUYKJjWi6ppRk2e7mYQHNqHbny7ZtOdslE3QKcbea/m239Y8EbN4KmYx2+fhAr+Hr1b7iTqajWmbN1M2uyQAHcevYLbpYZ7fw2yfh5ij+WvK5vFYmn0yqhkbO+ne67N/PQTr3eAz8zJU5Gp03zZ3VNT3nmy1cxwFVTPcUc3k2eANoHu34Wx/fDBFvfzOxUikEZC918qCDbSlbdrSNhUK5vo8zPN7+WS1l6ij2TT4sx2WMyETECYZijmmQYvfy9mmQULDW6kM9x9Oyagw1+Jdbtbbrj8f7b0bkavBY0Lfe+jerei+HWCd+DDHsmd4YVIEVXsmqvYs8AbgcAd692oZFR6jbILOeNh7rRASGtalbNal3AA0x+dhrs/DTVV+DAXbLId30zbbLLtst90xgYhpcJbp4Ryfh+legwO2Q1vaZk3col0CTQFo6O9A93egeRomNaBOOdvd5XrWtehoJ7pvB7q/Awb2Fn1fsLypCKOqI6iaWTCpEZwM+vB76F3Povu2g13+16QjZO+1CSCh4fWUzespGw/QZBrMNj0sDHi5vtJHv6PZbtm0Ww4dlk2nrUsyQ6tCQaPXIGJ6mOU1iJgGNtBuOaxPZnjIssfHLs8Tlob+nej+neidf4KKOlTNGajqWahpC92XDOxB9+9GR/dCdG9pXoyVAcGpqKrp7tTZpEbwBiG+H324Hf3+qxDdM2ESAsZKgk6ZsIHtlsN2y+HJuEWlgtNND2eYBgsCXq4LmVjA3ozD3ozDnoxDZ8ahy9bjKiuuSsFUr8F0r8EMj/tv2KNIaffG2u2Wu3nqrowjjwkfrxIH0IkD6H1rQXmgcjqquhFVNdMNQoYPEj3o2D6I7UfH90O8Z3xtXGmYbvAM1kNoKip0KgSnguFx6zywB6f9aejfBXay2LUtKRJ0ylRUw5tpmzfTNmDhVzDTazBj8OvyoEmdoQA46Gi6bM0B2+GArTloOxx2NH2OJqZzu+TpAaoNRbWhqDUUdR5FncdgikdR7zGoNBRprXl/MDg+n7DYnXHYZ2uZMitF2oaBXeiBXYN/PwXBMKpyGoROdTO6gotRhom2Ym7ASh6CZC+kDqFT/W6KtjWQ+xGEGQKzCvyTUP5aCNS6/1bUoQK1aO1A8hA6ts99KN6e5yHaKQkTWZKgM0Gk9AcjoSNMoN6jmOo1CHsM6gxFo99gssekajAgZbQmNhh8YlqTdDQpDSnA0u50naPBNCwyQS+Gcvcw8yuFX4FPKYIKQoYipBShY963z9EcHAx2m9MOqzMW+2xNryMBpnxpNxkh3gW88UEg8tdAcAqqog78tajqCATmocwqlOFBa+3eNJmJD/6bRDtpd5rOsdyApB3QmnSgAlWTdt/X8ILHhzJ87sK+t2LwK4gyPG6N0lFI9aGTh9CxTjjQhpPogcSB8l2LKiIJOhOYBey1NXttG/5mssqLmxVWbSgqDUVIQdBQBI4JJiGlUIDHAFNp0h6Fo913SmtN1IG0dohriA8Gr4E8jaBEKdPuli+pXnTvtiPf+cCREYkZQpnBDwKHYbpfHr+73qIMFAba9IJOuUHIsdytfuxesFPuhqeZBFjx/I2gxJAk6IiTygCHHM2hEW3Do6iqrmRgIIaEEpFzVuxoCvZwZ5dG4a+uJX24dwSvFsVgFLsCQgghJg4JOkIIIQpGgo4QQoiCKZs1nXLce00IIUpFJNIIDL/3mox0hBBCFEzZjHQmyt5rQggxHo107zUZ6QAXLVpY7CqUPGnD7Ekb5oa0Y/by2YYSdICLPiEnabakDbMnbZgb0o7Zy2cbStARQghRMCUfdBoaZhb1+JFI49GsjWK9R7HL50K2dSj18rmQbR1KvXwuZFuHYpfP1XvkU8kHHSGEEKWjbLLXTNM3xpIK29GYPh9j2avpSE66W35ssn2PYpfPtg1zUYdSLy9tmJvPUrE/z8Uun5v3yK4Nh7sWq3kXXlHSu+IFQ5V85sZbil0NIYQQx3js4fuJx6InfL/kgw64gcey5B4dIYQYD0zTd9KAA2UyvfZhv5wQQojCG+pGfUkkEEIIUTASdIQQQhSMBB0hhBAFI0FHCCFEwXhOnXH694tdCSHKVTAY5MorLqOmpoazzzqThQsuoG3z2wUrXy6ybYe7VtwJQFdXN6ZpMu+jH2HpFZfz7tZtWJaVr2qPO5HGBr70hRtpXbd+VOVyeR6WRfaaEOPVv972NR745X/R3d0DQDg8hTu+czsrfnRPQcqXi2zboakpws1NEW7+yk0A7O/qYsWP7iEej+etzuPJV778RQDC9WEqQ6FRl8/leTjhg06ksYHl1y7jpyt/PqpywWCQS5dcTFdXD5WVIerrwzzy6GN5quX4lG0b3LXiTlaveYG1rW6vq6V5DnNb5vDbRx8ri4vBksUX093VffSDChz9//wLPsa69RvyWr5c5KIdfvvI72hdt4H6cJhoNErHzl35qu649MCDDwFuW1665JJRlc31eThhg854ivylSnqfQ1uw4AJaW0+cxmhr28yll14y7Ic12/LlIlft0N3dc9yFU4xMrs/DCZtI8MCDD/HAgw+dtDGHM1zknwhy0Qa/feR3/Mut/8YPV9zNt77979x62zfK6qIwt6X5pL9PV1cPTZFI3suXC2mH4sp1+0/YoJONBQsuoH3w0azHOhL5J4JctUF3dw9tm98uu+mOYDAIQDR64m4Z8XicysqhR9fZli8XuW6HSGMDLc1zclK3iSAf56EEnTGQnpe0wXBG8mE88oHOR/lykat2mNXUxPwLPkYsHqdt89vc8LnPsGTxxdlXsMzl4zyUoDNK0gOV3qcoPavXvMC69RuOdpQeefQxrrnmKjnvikCCzihJD1R6n6L0nOyeks1t7jknCkuCjiiacu59RqOxYV8zVJZetuXLRT7bob1jJ01NMhU8lHy0vwQdUTTl3PuMx+NEozHC4fAJPwuHpwz7Yc62fLnIRTusvO/eCZNVmmv5OA9L8j6dlffdO6rXx2Nxvn3H93Jy7HLpgWbThvnufV566eIxlR1v2js6TjoVWV8fZvMIthDJtny5yLYd4rE4HTt3nrR8e3tHLqpY1nJ9HpZk0Ln1tm8U7djHRv6/TfMtpR5oNm2YizZYed+9PPLIY2V9g2Nr63rmtjSz5vkXj/t+S/Mc/vDk03kvXy6ybYe1retOmmk5f/7HeOqpVbmqZtnK9Xko02tjID1Q6X2OxJrnXyQcnkI4POXo9yKNDURjsROC7cr77j26KeVYypezbNuxre3tExJUrlm2lO7unhMupOUuFBo6wacQ52FJjnSKTXqg0vscqR+suJvl115NV1c34AbVk20TFI/FicZOHCGOtHy5y6Ydj4zGj6wVhoJBurq7J1Q7fuXLX6SyspLm5jlUVoa44zu3093Vzeo1Lxw3W1GI81DNu/AKPaaSZeKaZUtZuGD+h675rLzv3pOuCd214k5+uvLnRy+ckcYGbrjhsxPqRB5pG5ysDSONDTQ1RY4LWtcsW0pLS/OEakMhJpoJO9IZT5G/VEnvUwgxWhN+pCOEEKJwJJFACCFEwUjQEUIIUTASdIQQQhSMBB0hhBAFI0FHCCFEwUjQEUIIUTASdIQQQhSMBB0hhBAF8/+FedjLKCW6iAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 460.8x403.2 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xpl_shp = np.linspace(-1,1,101)\n",
    "nBasisPlot = basis(2, \"L\")\n",
    "N = nBasisPlot.Ns(xpl_shp)\n",
    "plt.figure()\n",
    "plt.plot(xpl_shp, N[0])\n",
    "plt.plot(xpl_shp, N[1])\n",
    "plt.plot(xpl_shp, N[2])\n",
    "# plt.plot(xpl_shp, N[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 18.25766979,  81.        , 143.74233021])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xi = leggauss(3)[0]\n",
    "(162+0)/2.+xi*(162)/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loc_mat(nodes, ElFlag, Ndnode):   #Computes the element quantities when supplied with gloabl nodes\n",
    "    xi=GP.xi;W=GP.wght\n",
    "    if etypes=='L':\n",
    "        x = np.array(nB.Ns(xi)).T @ nodes\n",
    "        print(\"Gauss Points: \",x)\n",
    "        Je = np.array(nB.dN(xi)).T @ nodes\n",
    "        print(\"Element: Jacobian = \", Je)\n",
    "        a1 = np.array(B.dN(xi)).reshape(-1,1,W.size)  \n",
    "        a2 = a1.reshape(1,len(a1),-1).copy()\n",
    "        a1 *= geom.E/Je*A(x)*W       # multiply by weights, jacobian, A, in bitwise fashion \n",
    "        mat = np.tensordot(a1, a2, axes=([1,2],[0,2])) #+ np.tensordot(b1,b2,axes=([1,2],[0,2]))     #could potentially use einsum to clean it up (but works for now!)\n",
    "        elemF = (np.array(B.Ns(xi))*W*Je*body_force(x)).sum(axis=1).flatten()\n",
    "        return mat, elemF\n",
    "\n",
    "def fsample(x, dof, nodes, sample_type): #giving out the displacement, strain and stress at x: GP and dof value of disp. at gauss point\n",
    "    if sample_type=='spline':\n",
    "        if int(El[1]) >5:\n",
    "            k0=5\n",
    "        else:\n",
    "            k0 = int(El[1])\n",
    "        dspl=splrep(nodes,dof,k=k0)\n",
    "        return (x,splev(x,dspl))\n",
    "    elif sample_type=='vandermonde':\n",
    "        Xinew = np.linspace(-1, 1-1.e-13, 20)\n",
    "        N = np.array(B.Ns(Xinew)).reshape(int(El[-1])+1, Xinew.size,-1)\n",
    "        N = np.einsum('jik',N)\n",
    "        Nx = np.array(nB.Ns(Xinew)).reshape(int(MapX[-1])+1,Xinew.size,-1)\n",
    "        Nx = np.einsum('jik', Nx)\n",
    "        dofarang = np.vstack((np.arange(k,k+dof.size-1,int(El[-1])) for k in range(int(El[-1])+1))) #arange dofs for multiplication with n\n",
    "        dofarang=dofarang.reshape(int(El[-1])+1,1,-1)\n",
    "        xarang = np.vstack((np.arange(k,k+nodes.size-1,int(MapX[-1])) for k in range(int(MapX[-1])+1)))\n",
    "        xarang = xarang.reshape(int(MapX[-1])+1, 1, -1)\n",
    "        dofs=dof[dofarang]\n",
    "        nodesx=nodes[xarang]\n",
    "        soln = np.einsum('ikj,kmj->imj',N,dofs)\n",
    "        xpl = np.einsum('ikj,kmj->imj',Nx,nodesx)\n",
    "        return (xpl.T.ravel(),soln.T.ravel()) # return the solution evaluated at the given points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def A(x):\n",
    "    return geom.Ao*np.exp(-x/geom.l)\n",
    "\n",
    "def body_force(x):\n",
    "    return -geom.rho*A(x)*geom.g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "El = \"L1\"\n",
    "MapX = \"L1\"\n",
    "Nel = 2 #number of elements\n",
    "Np = 3\n",
    "etypes = El[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\users\\bshri\\appdata\\local\\programs\\python\\python37\\lib\\site-packages\\ipykernel_launcher.py:8: FutureWarning: arrays to stack must be passed as a \"sequence\" type such as list or tuple. Support for non-sequence iterables such as generators is deprecated as of NumPy 1.16 and will raise an error in the future.\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "if El[0]=='L':                   #1D Mapping function for the physical coordinates\n",
    "    nB=basis(int(MapX[-1]), MapX[0])\n",
    "    geom=geometry(Nel, int(MapX[-1])) \n",
    "    probsize=geometry(Nel, int(El[-1]))\n",
    "    B=basis(int(El[-1]), El[0])     #Basis for FE fields (isoparametric) \n",
    "    GP=GPXi(Np) \n",
    "    nodes = np.linspace(0.,geom.L,geom.nNNodes)\n",
    "    elems = np.vstack((np.arange(k,k+nodes.size-1, int(MapX[1])) for k in range(int(MapX[1])+1))).T\n",
    "    globK = 0*np.eye(probsize.nNNodes)   # This changes because now we no longer have iso-p map\n",
    "    globF = np.zeros(probsize.nNNodes)    \n",
    "    prescribed_dof=np.array([[0,0]])\n",
    "\n",
    "dof=np.inf*np.ones(globK.shape[0])\n",
    "prescribed_forc = np.array([[-1, -geom.m*geom.g]])\n",
    "dof[prescribed_dof[:,0]] = prescribed_dof[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1],\n",
       "       [1, 2]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "elems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gauss Points:  [ 18.25766979  81.         143.74233021]\n",
      "Element: Jacobian =  [81. 81. 81.]\n",
      "Gauss Points:  [180.25766979 243.         305.74233021]\n",
      "Element: Jacobian =  [81. 81. 81.]\n",
      "Kglobal = \n",
      " [[ 2.06095915e+10 -2.06095915e+10  0.00000000e+00]\n",
      " [-2.06095915e+10  2.40163341e+10 -3.40674256e+09]\n",
      " [ 0.00000000e+00 -3.40674256e+09  3.40674256e+09]]\n",
      "fGlobal = \n",
      " [-1.02262169e+10 -7.38220848e+09 -9.41048743e+08]\n"
     ]
    }
   ],
   "source": [
    "for k in range(elems[:,0].size): #Assembly of Global System\n",
    "    elflag = 'nn'\n",
    "    ndnode = 1\n",
    "    elnodes=elems[k]\n",
    "    if etypes=='L':\n",
    "        elemord=int(El[-1])+1     # basically selecting the corresponding dof\n",
    "        globdof = np.array([k*(elemord-1)+i for i in range(elemord)], int)\n",
    "#         print(globdof)\n",
    "    if etypes != '2D':                                                         # 1D lagrange or legendre polynomials \n",
    "        nodexy=nodes[elnodes]\n",
    "    else:\n",
    "        nodexy=nodexy[:,elnodes]\n",
    "    kel, fel = loc_mat(nodexy,elflag,ndnode)\n",
    "    globK[np.ix_(globdof,globdof)] += kel                                      #this process would change in vectorization \n",
    "    globF[globdof] += fel\n",
    "\n",
    "globF[int(prescribed_forc[:,0])] += prescribed_forc[:,1]    \n",
    "fdof=dof==np.inf                                                               # free dofs\n",
    "nfdof=np.invert(fdof)   \n",
    "print(\"Kglobal = \\n\",globK)\n",
    "print(\"fGlobal = \\n\",globF)\n",
    "if etypes=='L':  \n",
    "    AA=globK[np.ix_(fdof,fdof)]\n",
    "    bb=globF[fdof]-globK[np.ix_(fdof,nfdof)] @ dof[nfdof]                                                                           \n",
    "    dofme=nla.solve(AA, bb)\n",
    "dof[fdof]=dofme.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xplot=np.linspace(0.,geom.L,10*geom.nNNodes-1)  #sampling points for the exact solution\n",
    "fittype='vandermonde'\n",
    "plot_it  = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_it:\n",
    "    if etypes=='L':\n",
    "        if fittype == 'spline':\n",
    "            xplot, uapp = fsample(xplot, dof[0:len(dof):int(El[1])], nodes[0:len(nodes):1], fittype)\n",
    "        elif fittype == 'vandermonde':\n",
    "            xplot, uapp = fsample(xplot, dof, nodes, fittype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uexact = list(map(uplot, xplot))\n",
    "uexact = -1*np.array(uexact)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "plt.plot(xplot, uexact, \"b-\", lw=2, label=r\"$\\widetilde{{u}}(x)$\")\n",
    "plt.plot(xplot, -uapp,'--r', lw=2 , label=r'$u^h_{0:1d}(x)$'.format(int(El[-1])), markevery=4)\n",
    "plt.legend(loc=0, fontsize=22, fancybox=False, edgecolor=\"white\", shadow=True)\n",
    "plt.ylabel(r\"$u^h_{0:1d}(x)/\\widetilde{{u}}(x)$\".format(int(El[-1])), fontsize=24)\n",
    "plt.xlabel(r\"$x$\", fontsize=24)\n",
    "plt.grid(which=\"major\", color=\"k\")\n",
    "plt.savefig(\"Displacement_deg_{}_nel_{}.eps\".format(int(El[-1]), Nel))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}