FEM 2D Euler Beam Quadratic Load.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian_quadrature(n):\n",
    "        # n: number of points\n",
    "        # w: weight \n",
    "        # x: point of integration\n",
    "        \n",
    "        if n == 1:\n",
    "            w = np.array([2])\n",
    "            x = np.array([0])\n",
    "        elif n == 2:\n",
    "            w = np.array([1, 1])\n",
    "            x = np.array([-1/np.sqrt(3), 1/np.sqrt(3)])\n",
    "        elif n == 3:\n",
    "            w = np.array([5/9, 8/9, 5/9])\n",
    "            x = np.array([-np.sqrt(3/5), 0, np.sqrt(3/5)])\n",
    "        elif n == 4:\n",
    "            w = np.array([ (18-np.sqrt(30))/36, (18+np.sqrt(30))/36, (18+np.sqrt(30))/36, (18-np.sqrt(30))/36 ])\n",
    "            x = np.array([ -np.sqrt(3/7 + 2/7*np.sqrt(6/5)), -np.sqrt(3/7 - 2/7*np.sqrt(6/5)), np.sqrt(3/7 - 2/7*np.sqrt(6/5)), np.sqrt(3/7 + 2/7*np.sqrt(6/5)) ])\n",
    "        \n",
    "        return x, w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Input data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 1                             # Beam length    \n",
    "q = 100                           # load\n",
    "nel = 512                         # number of elements\n",
    "nnos = nel + 1                    # number of nodes   \n",
    "alldof = np.linspace(1, 1,2*nnos) # all degrees of freedom\n",
    "kg = np.zeros((2*nnos,2*nnos))    # global stiffness matrix pre-allocation\n",
    "I = 100                           # inertia\n",
    "E = 210e9                         # young's modulus   \n",
    "coord = np.zeros((nnos, 3))       # coordinate matrix pre-allocation\n",
    "inci = np.zeros((nel, 6))         # incidence matrix pre-allocation\n",
    "f = np.zeros((2*nnos, 1))         # external load vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coordinate matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(0, nnos):\n",
    "    coord[i,0] = i + 1      # node number\n",
    "    coord[i,1] = i*L/nel    # node position\n",
    "    coord[i,2] = 0\n",
    "l = coord[1,1] - coord[0,1] # element length"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Incidence matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(0, nel):   \n",
    "    inci[i,0] = i + 1          # element number\n",
    "    inci[i,1] = i + 1          # first node\n",
    "    inci[i,2] = i + 2          # second node\n",
    "    inci[i,3] = coord[i,1]     # first node coordinate\n",
    "    inci[i,4] = coord[i+1,1]   # second node coordinate\n",
    "    inci[i,5] = I              # inertia of beam section"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boundary conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#   bc=[node | degree of freedom | value]\n",
    "#\n",
    "#   Degree of freedom 1 --> y\n",
    "#   Degree of freedom 2 --> oz\n",
    "\n",
    "bc = np.array([[1,1,0],[1,2,0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mask stiffness matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = np.zeros((2*nnos,2*nnos))\n",
    "for i in range(0, np.size(bc,0)):\n",
    "    if bc[i,1] == 1:\n",
    "        mask[2*bc[i,0]-2,2*bc[i,0]-2] = 1\n",
    "    elif bc[i,1] == 2:\n",
    "        mask[2*bc[i,0]-1,2*bc[i,0]-1] = 1\n",
    "mask = ma.masked_equal(mask, 1)\n",
    "mask = ma.mask_rowcols(mask)\n",
    "mask = (mask==False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mask load vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "maskv = np.zeros(2*nnos)\n",
    "for i in range(0, np.size(bc,0)):\n",
    "    if bc[i,1] == 1:\n",
    "        maskv[2*bc[i,0]-2] = 1\n",
    "    elif bc[i,1] == 2:\n",
    "        maskv[2*bc[i,0]-1] = 1\n",
    "maskv = ma.masked_equal(maskv, 1)\n",
    "maskv = (maskv==False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Load vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#   F = [node | degree of freedom | value]\n",
    "#\n",
    "#   Degree of freedom 1 --> Fy\n",
    "#   Degree of freedom 2 --> Mz\n",
    "def Q(x, q, L):\n",
    "    Q = q + 3/4*q/L**2*x**2 - 7/4*q/L*x\n",
    "    return Q\n",
    "\n",
    "def phi(x, L):\n",
    "    phi = np.array([[2*x**3/L**3 - 3*x**2/L**2 + 1], [x - 2*x**2/L + x**3/L**2], [3*x**2/L**2 - 2*x**3/L**3], [x**3/L**2 - x**2/L]])\n",
    "    return phi\n",
    "\n",
    "\n",
    "n = 3\n",
    "x,w = gaussian_quadrature(n)\n",
    "fg = np.zeros((2*nnos,1))\n",
    "\n",
    "for i in range(nel):\n",
    "    node1 = inci[i,1]  # first node element\n",
    "    node2 = inci[i,2]  # second node element\n",
    "\n",
    "    loc = [2*node1-2,2*node1-1,2*node2-2,2*node2-1]\n",
    "\n",
    "    coord1 = coord[int(node1) - 1, 1] # node 1 position\n",
    "    coord2 = coord[int(node2) - 1, 1] # node 2 position\n",
    "\n",
    "    fe = np.zeros((4,1))   # local load vector \n",
    "\n",
    "    Q_ = 0\n",
    "    for j in range(n):\n",
    "\n",
    "        epsilon = (coord2 - coord1)/2*x[j] + (coord2 + coord1)/2     \n",
    "        fe = fe + ((coord2 - coord1)/2)*Q(epsilon  , q, L)*phi(epsilon, L)*w[j] \n",
    "    \n",
    "    fg[[int(loc[0]),int(loc[1]),int(loc[2]),int(loc[3])]] = fg[[int(loc[0]),int(loc[1]),int(loc[2]),int(loc[3])]] + fe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Global matrix assembly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(nel):\n",
    "    node1 = inci[i,1] # first node element\n",
    "    node2 = inci[i,2] # second node element\n",
    "    \n",
    "    # local stiffness matrix\n",
    "    inertia = inci[i,5]\n",
    "    ke = E*inertia/l**3*np.array([[12, 6*l, -12, 6*l], [6*l, 4*l**2, -6*l, 2*l**2], [-12, -6*l, 12, -6*l],  [6*l, 2*l**2, -6*l, 4*l**2]])\n",
    "    \n",
    "    # localization vector\n",
    "    loc = [2*node1-2,2*node1-1,2*node2-2,2*node2-1]\n",
    "    \n",
    "    # global stiffness matrix \n",
    "    kg[[[int(loc[0])],[int(loc[1])],[int(loc[2])], [int(loc[3])]], [int(loc[0]),int(loc[1]),int(loc[2]),int(loc[3])]] = kg[[[int(loc[0])],[int(loc[1])],[int(loc[2])], [int(loc[3])]], [int(loc[0]),int(loc[1]),int(loc[2]),int(loc[3])]] +  ke\n",
    "    \n",
    "kg_aux = kg[mask.data]\n",
    "kg_aux = np.reshape(kg_aux, (2*nnos-np.size(bc,0), 2*nnos-np.size(bc,0)))\n",
    "f_aux  = fg[maskv.data]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Displacement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = np.zeros((2*nnos, 1))\n",
    "u[maskv.data] = np.linalg.solve(kg_aux, f_aux)\n",
    "u_ = np.reshape(u,(nnos,2))\n",
    "displacement = u_[:,0]\n",
    "rotation = u_[:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Post-processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsoAAAGmCAYAAABlSKL8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXRV1fnG8e8hAQI3zGGSMAchEAHhMmkRFBGCosG2CFYRWpsWKrZi01ZBVBSqpeJAG63WQukQSlsBQYPiAKICIqgYSCBhTJhDIJBAQob9+wP1R+HIEG6yz02ez1pdy8DN2e/qfe7hzc1793aMMYiIiIiIyP+qZrsAEREREREvUqMsIiIiIuJCjbKIiIiIiAs1yiIiIiIiLtQoi4iIiIi4UKMsIiIiIuLCaqPsOM5fHMc56DhOSoCut8xxnKOO4yw9689fdRznC8dxNjqO8x/HccIDsZ6IiIiIVF6231GeCwwN4PVmAne7/PkDxphuxpiuwG7gvgCuKSIiIiKVkNVG2RjzAZBz5p85jtP+q3eG1zuOs8pxnE6XcL13geMuf37sq2s7QC1Ap6yIiIiIyHnZfkfZzcvARGNMT+CXQGIgLuo4zhxgP9AJmB2Ia4qIiIhI5RVqu4AzfTU7fA3w79Nv/gJQ86u/ux2Y5vJte4wxQy50bWPMOMdxQjjdJN8BzAlI0SIiIiJSKXmqUeb0O9xHjTHdz/4LY8xrwGuXc3FjTInjOP8CElCjLCIiIiLn4anRi69miXc4jvN9OD1T7DhOt8u55lfXiPr6v4HhQNplFysiIiIilZpjjL3PtTmOkwQMBCKAA8CjwHvAi0BzoDow3xjjNnLhdr1VnJ5BDgcOAz8ClgOrgLqAA3wBjP/6A34iIiIiIm6sNsoiIiIiIl7lqdELERERERGvsPZhvoiICNOmTRtby4uIiIhIFbF+/fpsY0zjS/0+a41yaWkpn376qa3lxaNGjRrF/PnzbZchHqNciBvlQtwoF+LGcZxdZfk+a6MXhYWFtpYWD8vIyLBdgniQciFulAtxo1xIIGlGWURERETEhbVGuX379raWFg9buHCh7RLEg5QLcaNciBvlQgLJWqOck5Nja2nxsKSkJNsliAcpF+JGuRA3yoUEkrVG+ciRI7aWFg9bsGCB7RLEg5QLcaNciBvlQgJJM8oiIiIiIi6sNcrNmze3tbR42NSpU22XIB6kXIgb5ULcKBcSSNYa5WrV9Ga2nCssLMx2CeJByoW4US7EjXIhgeQYY6ws7PP5TH5+vpW1xbv8fr8OopFzKBfiRrkQN8qFuHEcZ70xxn+p33fBt3Udx/mL4zgHHcdJ+Za/dxzHecFxnAzHcTY6jtPjUosQEREREfGai5l/mAsMPc/fxwIdvvpfPPDixSxcr169i3mYVDHDhw+3XYJ4kHIhbpQLcaNcSCBd1OiF4zhtgKXGmBiXv/sTsMIYk/TV11uAgcaYfee7Zvfu3c3nn39elpqlEjt06BCNGze2XYZ4jHIhbpQLcaNciJuyjl6EBmDtFkDmGV9nffVn5zTKjuPEc/pdZ0JCQvD7T9c7ceJEoqOjmTBhAgD9+/dn5syZ9O3bFwCfz8fKlSsZO3YsKSmnJ0CSkpJITk5m3rx5ACQkJNC8eXMmTZoEwODBg5k8eTIDBw4EICIigmXLljFq1KhvzoFfuHAhSUlJ3+y5OHXqVMLCwnj44YeB0z+VTpgwgdjYWAAiIyNZtGgRcXFxZGVlAZCcnExiYiJLliwBYMaMGRQUFDBt2jQARo4cyejRoxkxYgQAUVFRzJ8/n6FDh5KdnQ3AihUrmD59OsuXLwdg1qxZ7Nu3j5kzZwIwZswYYmNjGT16NAAxMTHMnTuXAQMG8PWc95o1a0hISGDVqlUAJCYmkpqayuzZswGIj4+nX79+jBs3DoBevXrx4osv0qdPH0pKSggJCWHt2rWMHz+edevWATBnzhxWr17Nyy+/XGHPU3p6Orm5uXqePP48VfTrafTo0URHR+t58vjzVNGvJ2MM1113nZ4njz9PFf16Sk1NJT4+Xs+Tx58nqNjXU1kF4h3lN4DfGmM+/Orrd4FfGWPWn++a+jCfuNGHMMSNciFulAtxo1yIm3L7MN9FyAJanvF1JLD3Qt9UvXr1ACwtlU1kZKTtEsSDlAtxo1yIG+VCzvZlVm6ZvzcQ7yjfDNwHDAP6AC8YY3pf6Jp+v9/oJz4RERERKQ+5J4qY+XYa/1i7m51P3VJu28MlAauBjo7jZDmO8yPHcX7qOM5Pv3rIm8B2IAN4BZhwMQt/PTcicqa4uDjbJYgHKRfiRrkQN8qFGGN4bUMWg2at4J9rdzP2mjZlvtYFP8xnjBl9gb83wM8udeGioqJL/RapAr7+cIPImZQLcaNciBvlomrLOHicKYtSWLM9h+4t6zN3XG9iWtTjsTJeLxC7XoiIiIiIWHPyVAmz30vnlVXbqV0jlBkjrmJUr5ZUq+Zc1nWtHWGtfZTFjfa/FDfKhbhRLsSNclH1vJt6gKmLN7Hn6Em+2yOSh4Z1IiK85v88xuauF2Vy8OBBW0uLhyUmJtouQTxIuRA3yoW4US6qjj1HTxI/71N+9NdPqV0jhH/F9+WZkd3OaZIvh7VGOTe37Ft1SOX19YbrImdSLsSNciFulIvKr6iklJdWbuPGZ1ayKj2bXw/txBv396dPu0YBX0szyiIiIiISFD7ZkcOURV+y9UAegzs35dHhnYlsULvc1rPWKLdo0cLW0uJhM2bMsF2CeJByIW6UC3GjXFROh/MKmfFmGv/dkEWL+rX48xg/N3ZuWu7rWmuUS0tLbS0tHlZQUGC7BPEg5ULcKBfiRrmoXEpLDfPXZfL0sjTyC4sZP7A9E2+IonaNimlhre164fP5TH5+vpW1xbv8fj86sVHOplyIG+VC3CgXlcemvblMWZTCZ7uP0qdtQ56Mi6FD0zplulZZd73QjLKIiIiIeMbxgiJmLd/KXz/eSYPaNZg1shsjrm6B41zenshlYa1RbtCgga2lxcNGjhxpuwTxIOVC3CgX4ka5CF7GGN74ch9PLN3MweOF3Nm7Fb8a0ol6tatbq8na6EXXrl3Nxo0brawt3pWZmUnLli1tlyEeo1yIG+VC3CgXwWlndj6PLE5hVXo2MS3q8mTcVXRvWT9g1w+6A0e2bdtma2nxsBEjRtguQTxIuRA3yoW4US6CS0FRCc8u38pNz33A57uP8tjwziz+2XcC2iRfDs0oi4iIiEiF+2DrIaYuTmHn4RMM73YFj9wcTZO6YbbL+h/WGuWaNQN3vKBUHlFRUbZLEA9SLsSNciFulAvvO3CsgGlLN/PGxn20jfDxtx/1pn+HxrbLcmVtRtnv9xtt3yIiIiJSNRSXlDJv9S5mLd/KqZJS7rs+ivjr2hFWPaTc1w66GeX09HRbS4uHDR061HYJ4kHKhbhRLsSNcuFNG3Yf4dY/fMS0pZvp2boByx+4jvsHdaiQJvlyWBu9KC4utrW0eFh2drbtEsSDlAtxo1yIG+XCW46eOMXTy7Ywf91umtYJI/EHPYiNaWZlT+Sy0If5RERERCSgjDH8d8MefvtmKkdPFvGja9vyi8FXEl4zuFpPazPKPXr0MBs2bLCytnhXXl4e4eHhtssQj1EuxI1yIW6UC/u2HjjOlIUpfLIzhx6t6vNk3FV0vqKu1ZqCbkZ53759tpYWD5s+fbrtEsSDlAtxo1yIG+XCnhOnivltcirDnl/F1oPHeer2q/jPT6+x3iRfDmuN8rFjx2wtLR62fPly2yWIBykX4ka5EDfKhR1vb9rP4Fkf8KeV27m9Rwvee3Ago3q3olq14JhF/jbBNSgiIiIiIp6RmXOCx5ds4p3Ug3RsWod//7Qfvdo0tF1WwFhrlHUOu7iZNWuW7RLEg5QLcaNciBvlomKcKi7llVXbmf1eOtUch4eHdWLctW2pHmJtWKFcWGuUi4qKbC0tHqbZdXGjXIgb5ULcKBflb/W2wzyyOIWMg3kM7dKMqcM7c0X9WrbLKhfWdr3w+XwmPz/fytriXX6/H53YKGdTLsSNciFulIvyk51XyIw3Unntsz20bFiLabfGcH2nJrbLuihl3fVCM8oiIiIi8q1KSg1Jn+zmd8vSOFlUwn3XR/Gz66OoVcPbp+oFgrVGuVGjRraWFg8bM2aM7RLEg5QLcaNciBvlIrBS9uQyeVEKX2QepV+7RjwRF0NUk6qzT7W1RrlevXq2lhYPi42NtV2CeJByIW6UC3GjXATGsYIiZr29lXmrd9LQV5Pn7ujObd2vCJqjpwNFM8riKZotEzfKhbhRLsSNcnF5jDEs2biPJ5ZuJjuvkLv7tubBmzpSr1Z126VdFs0oi4iIiEiZbT+Ux9TFm/gwI5urWtTj1Xv8dI2sb7ssq6w1yrVqVc5tROTyxMTE2C5BPEi5EDfKhbhRLi5dQVEJie9n8NLK7dSsXo0nbuvCnX1aExLkp+oFgrXRC7/fb/SrERERERF73t9ykEcXb2J3zgniul/BwzdH06ROmO2yAq6soxfWjk/ZsmWLraXFwwYMGGC7BPEg5ULcKBfiRrm4OPtyTzL+7+sZN2cdoSEO/7y3D8+NurpSNsmXw9roRWlpqa2lxcP0AU9xo1yIG+VC3CgX51dcUsrcj3fy7PKtFJcaEoZ05N7+bakZWvn3RC4LfZhPREREpApYvyuHyQtTSNt/nBs6NeHxW7vQsmFt22V5mrUZ5Z49e5r169dbWVu8q7i4mNBQ/fwm/0u5EDfKhbhRLs51JP8UTyWn8a9PM2leL4xHh3dhSJemVWpP5KCbUc7KyrK1tHhYQkKC7RLEg5QLcaNciBvl4v+VlhoWrMvkhmdW8J8NWcRf1453Jg1gaEyzKtUkXw5rP3Ll5eXZWlo8bNWqVbZLEA9SLsSNciFulIvT0vYfY8rCFD7ddQR/6wY8OSKGTs3q2i4r6Oh3EyIiIiKVRH5hMc+/m86rH+6gblgov/teV77XI5Jq2hO5TKw1yq1atbK1tHhYYmKi7RLEg5QLcaNciJuqmgtjDG9t2s/jSzazL7eAUb1a8uuhnWjgq2G7tKBmbUa5oKDA1tLiYampqbZLEA9SLsSNciFuqmIudh8+wQ/nruOnf99AvVrV+e/4fjz13a5qkgPA2q4XPp/PaK9DOZvf70cnNsrZlAtxo1yIm6qUi8LiEl75YDuz38sgtJrDA4OvZOw1bQgNsfY+qGeVddcLzSiLiIiIBJmPM7KZsjiF7YfyGXZVMx65pTPN69WyXValY61RjoiIsLW0eFh8fLztEsSDlAtxo1yIm8qei4PHC5j+RiqLP99Lq4a1mTuuFwM7NrFdVqVlrVEODw+3tbR4WL9+/WyXIB6kXIgb5ULcVNZclJQa/rF2FzPf2kJhUSn3D+rAhIHtCauuo6fLk2aUxVOq0myZXDzlQtwoF+KmMuZiY9ZRJi9M4cs9uXwnKoJpt3WhXWO94XgpNKMsIiIiUonknizi929t4e9rd9E4vCazR1/NLV2b61S9CmStUfb5fLaWFg/r1auX7RLEg5QLcaNciJvKkAtjDIs+38P0N1LJyT/FPf3aMOmmK6kbVt12aVWOtdELv99vKtuvRkREREQuR8bBPB5ZlMLq7Yfp1rI+0+NiiGlRz3ZZQa+soxfWNtqrihuCy4X16dPHdgniQcqFuFEuxE2w5uLkqRJmvpVG7PMfsGlvLk/GxfDa+GvUJFumGWXxlJKSEtsliAcpF+JGuRA3wZiL99IOMHXxJrKOnOT2q1vw0LBoGtepabssQY2yeExIiLa5kXMpF+JGuRA3wZSLvUdP8viSTby16QBRTcKZH9+Xvu0a2S5LzqAZZREREZEKVFRSyl8+3MHz76ZTagz3D+rAvd9pR41QHT1dXoJuRnn37t22lhYPGz9+vO0SxIOUC3GjXIgbr+di3c4cbnnhQ36bnMY17Rux/IEBTBgYpSbZo6yNXuiwEXGzbt062yWIBykX4ka5EDdezcXhvEKeSk7j3+uzaFG/Fq+M8TO4c1PbZckFaEZZREREpJyUlhr+9WkmTy9LI6+gmJ8OaM/9g6KoXUMtWDCw9iy1adPG1tLiYXPmzLFdgniQciFulAtx46VcbN57jCmLvmTD7qP0btuQJ+NiuLJpHdtlySWwNhCTl5dna2nxsNWrV9suQTxIuRA3yoW48UIu8gqLeWLpZob/4UN2HT7BM9/vxr/i+6pJDkLWdr3w+XxGc8pyNr/fj3ZDkbMpF+JGuRA3NnNhjOHNL/czbekmDh4vZHTvVvxqSEfq165hpR75f2Xd9UIDMiIiIiKXadfhfB5ZvIkPth6ic/O6vHRXT65u1cB2WXKZrDXKTZo0sbW0eNjEiRNtlyAepFyIG+VC3FR0LgqLS3hpxXb+uCKDGiHVeHR4Z+7u25rQEG33VhlcVKPsOM5Q4HkgBPizMeaps/6+FfBXoP5Xj/mNMebN810zLCysTAVL5RYdHW27BPEg5ULcKBfipiJz8WF6No8sTmFHdj63dG3OI7d0pmld9TeVyQVnlB3HCQG2AoOBLGAdMNoYs/mMx7wMfGaMedFxnM7Am8aYNue7rmaUxY1mDsWNciFulAtxUxG5OHCsgCffSGXJF3tp06g2026L4borG5frmnJ5ynNGuTeQYYzZ/tVC84HbgM1nPMYAdb/673rA3kstRERERMTLiktK+duaXTzz9lZOlZTyixs78NMB7QmrHmK7NCknF9MotwAyz/g6C+hz1mMeA952HGci4ANudLuQ4zjxQDycHr3w+0839hMnTiQ6OpoJEyYA0L9/f2bOnEnfvn0B8Pl8rFy5krFjx5KSkgJAUlISycnJzJs3D4CEhASaN2/OpEmTABg8eDCTJ09m4MCBAERERLBs2TJGjRpFRkYGAAsXLiQpKYkFCxYAMHXqVMLCwnj44YcBGD58OBMmTCA2NhaAyMhIFi1aRFxcHFlZWQAkJyeTmJjIkiVLAJgxYwYFBQVMmzYNgJEjRzJ69GhGjBgBQFRUFPPnz2fo0KFkZ2cDsGLFCqZPn87y5csBmDVrFvv27WPmzJkAjBkzhtjYWEaPHg1ATEwMc+fOZcCAAd+ccLhmzRoSEhJYtWoVAImJiaSmpjJ79mwA4uPj6devH+PGjQOgV69evPjii/Tp04eSkhJCQkJYu3Yt48eP/+ZUozlz5rB69WpefvnlCnuevv7/RM+Tt5+nin49ZWZm4vf79Tx5/Hmq6NfTtddeywMPPKDnyePPU0W/njIzM3nggQcC/jz5WnUhfOC9bNp7jOrZ6YRveZPv/uRVXnj2GT1PQfB6KquLGb34PjDEGHPvV1/fDfQ2xkw84zGTvrrWM47j9ANeBWKMMaXfdt2ePXua9evXl7lwqZyKi4sJDdVmLPK/lAtxo1yIm0DnIvdEEU+/lUbSJ7tpUqcmU2/pwrCrmuE4TsDWkPJX1tGLi/lIZhbQ8oyvIzl3tOJHwAIAY8xqIAyION9F09LSLr5KqTK+/qlS5EzKhbhRLsRNoHJhjOG/67O44ZkVzP9kN+Ouacs7kwZwc9fmapKrkIv5kWsd0MFxnLbAHmAUcOdZj9kNDALmOo4TzelG+VAgCxURERGpCOkHjjNlUQprd+Rwdav6zPtRb7pcUc92WWLBBRtlY0yx4zj3AW9xeuu3vxhjNjmOMw341BjzOvAg8IrjOA9w+oN9Y80FZjqqVdP+gnIun89nuwTxIOVC3CgX4uZycnHiVDEvvJvBn1dtx1czlN/efhV3+FtSrZreQa6qrB1h7ff7jbb1ERERES9YvvkAj72+iT1HT/K9npE8FNuJRuE1bZclAVKeM8rlYufOnbaWFg8bO3as7RLEg5QLcaNciJtLzUXWkRPc+9dP+fG8T/HVDGHBT/rx++93U5MsgMUjrE+ePGlrafGwr7eDETmTciFulAtxc7G5OFVcyqsf7uD5d7fi4PBQbCd++J22VNfR03IG7asjIiIiVcqa7Yd5ZFEK6QfzuKlzUx69tQst6teyXZZ4kLVGuV27draWFg9LSkqyXYJ4kHIhbpQLcXO+XGTnFTLjzVRe27CHyAa1ePUeP4Oim1ZgdRJsrP1+ITc319bS4mHJycm2SxAPUi7EjXIhbtxyUVpq+MfaXQx6ZiVLvtjLz65vz/IHBqhJlguytuuFz+czXx+bKPI1v9+PdkORsykX4ka5EDdn5yJlTy6TF6XwReZR+rZryJNxMUQ1qWOxQrGhrLteaEZZREREKp3jBUU88/ZW5q3eSUNfDZ69oxtx3VvoVD25JNYa5WbNmtlaWjwsISHBdgniQcqFuFEuxM0vf5nAki/28sTSzRzKK+SuPq355U0dqVe7uu3SJAhZa5SrV1dg5VzNmze3XYJ4kHIhbpQLOduO7Hz+nlWPlM8/46oW9XhljJ9uLevbLkuCmGaUxVM0cyhulAtxo1zI1wqKSkhcsY2XVmyjqPAEj93ek7v6tiZER0/LV4LuZD4RERGRy7Vy6yGGPPcBL7ybTuxVzWjw8WzuuaaNmmQJCGujF3Xr1rW1tHjY4MGDbZcgHqRciBvlomrbn1vAtKWbePPL/bSL8PGPe/twbVQED33Rz3ZpUolYG73o0aOH2bBhg5W1xbvy8vIIDw+3XYZ4jHIhbpSLqqm4pJS5H+/k2eVbKS413Hd9FPED2lEzNARQLsRd0I1ebNmyxdbS4mEDBw60XYJ4kHIhbpSLqmf9riMM/8NHPPlGKr3bNmT5AwOYOKjDN00yKBcSWNpHWURERDzt6IlTPL0sjaRPMmleL4yX7urBkC7NtCeylDtrjXJoqHp0OVdERITtEsSDlAtxo1xUfsYY/rM+i98mp5F7sogf92/Lz2+8kvCa395DKBcSSNZmlP1+v9G2PiIiIuJmy/7jTFn0Jet2HqFn6wY8GRdDdHNtBCBlE3Qzytu3b7e1tHjYqFGjbJcgHqRciBvlonLKLyzmt2+mcvMLq8g4mMfvvtuVf/+k30U3ycqFBJK1+YfCwkJbS4uHZWRk2C5BPEi5EDfKReVijOGtTQeYtmQTe3MLuMPfkl/HdqKhr8YlXUe5kEDSoLCIiIhYlZlzgkdf38R7aQfp1KwOL4y+Gn+bhrbLErE3o9y1a1ezceNGK2uLd2VmZtKyZUvbZYjHKBfiRrkIfqeKS3ll1XZmv5dONcdh0uArueeaNlQPKftkqHIhboJuRjknJ8fW0uJhSUlJtksQD1IuxI1yEdw+3pZN7PMfMPOtLVzfsQnvPjiAe/u3u6wmGZQLCSxrjfKRI0dsLS0etmDBAtsliAcpF+JGuQhOh44X8ov5n3HnK2spKjHMGdeLF+/qSfN6tQJyfeVCAkkzyiIiIlLuSkoN/1y7i9+9tYXColLuvyGKCddHEVY95MLfLGKJtUa5efPmtpYWD5s6dartEsSDlAtxo1wEjy+zcpm86Es2ZuVybVQjpt0WQ/vG4eWylnIhgWStUa5WzdrUh3hYWFiY7RLEg5QLcaNceN+xgiKeeWsLf1uzi0bhNXl+VHdu7XZFuR49rVxIIFnb9cLn85n8/Hwra4t3+f1+dGKjnE25EDfKhXcZY3j9i708sTSVnPxC7u7bmgeHdKRuWPVyX1u5EDdl3fVCM8oiIiISMNsO5TF1cQofZRyma2Q95oztxVWR9WyXJVIm1hrlevX0opFzDR8+3HYJ4kHKhbhRLryloKiEP76fwZ9Wbqdm9Wo8ERfDnb1bEVKt/MYs3CgXEkjWRi+6d+9uPv/8cytri3cdOnSIxo0b2y5DPEa5EDfKhXe8n3aQqa+nkJlzkhFXt+DhYdE0rlPTSi3KhbgJugNH0tPTbS0tHhYbG2u7BPEg5ULcKBf27T16kp/+bT3j5q6jRkg1/vnjPjx7R3drTTIoFxJYmlEWERGRS1JUUsrcj3by7DtbKTWGhCEd+XH/dtQI1Y5WUrlYa5SrVy//T75K8ImMjLRdgniQciFulAs7Pt2Zw5RFKaTtP86gTk147NYutGxY23ZZ31AuJJCszSj7/X6j7VtERESCQ07+KZ5KTmXBp1lcUS+MR2/twk2dm5brnsgigRJ0M8oZGRm2lhYPi4uLs12CeJByIW6Ui4pRWmr417rd3PDMCl7bsIefDGjHOw8OYEiXZp5skpULCSRroxdFRUW2lhYPy8rKsl2CeJByIW6Ui/KXuu8YUxalsH7XEXq3acgTcTF0bFbHdlnnpVxIIOnDfCIiIvI/8gqLeW75VuZ8vJN6tarz++9347s9WnjyHWSR8qR9lMVTtP+luFEuxI1yEXjGGJal7OfxJZvZf6yA0b1b8euhHalfu4bt0i6aciFugm5G+eDBg7aWFg9LTEy0XYJ4kHIhbpSLwNp9+ATj5q5j/D820MBXg9cmXMNvb78qqJpkUC4ksKw1yrm5ubaWFg9bsmSJ7RLEg5QLcaNcBEZhcQkvvJvO4GdXsm5HDo/c0pkl911Lj1YNbJdWJsqFBJJmlEVERKqojzKyeWRRCtuz87n5quY8cktnmtULs12WiGdYa5RbtGhha2nxsBkzZtguQTxIuRA3ykXZHTxWwJNvpPL6F3tp3ag2f/1hbwZcWTnmepULCSRrjXJpaamtpcXDCgoKbJcgHqRciBvl4tKVlBr+tnonz7y9lcLiUn4+qAPjB7YnrHqI7dICRrmQQLK264XP5zP5+flW1hbv8vv96MRGOZtyIW6Ui0vzReZRJi/6kpQ9x+jfIYJpt8XQNsJnu6yAUy7ETVl3vdCMsoiISCWWe6KImW+n8Y+1u2kcXpM/3Hk1N1/VXHsii1wEa41ygwbB+WlaKV8jR460XYJ4kHIhbpSL8zPGsPCzPcx4M5Wc/FOMvaYNkwZfSZ2w6rZLK1fKhQSStdGLrl27mo0bN1pZW7wrMzOTli1b2i5DPEa5EDfKxbfLOHicKYtSWLM9h+4t6/NkXAwxLerZLqtCKBfiJugOHNm2bZutpcXDRowYYbsE8SDlQtwoF+c6eaqE3y1LI/b5VaTuO86MEVfx2vhrqkyTDMqFBJZmlEVERCqBd1MPMHXxJvYcPcl3e0Ty0LBORITXtF2WSFCz1ijXrKkXr5wrKirKdgniQcqFuFEuTttz9CSPvx4mjyIAACAASURBVL6JtzcfoEOTcP4V35c+7RrZLssa5UICydqMst/vN9q+RUREpGyKSkp59cMdPP9OOgD3D+rAj77Tlhqh1qYqRTwr6GaU09PTbS0tHjZ06FDbJYgHKRfipirn4pMdOdz8wiqeSk7jOx0iWD7pOsYPbK8mmaqdCwk8a6MXxcXFtpYWD8vOzrZdgniQciFuqmIuDucVMuPNNP67IYsW9Wvx5zF+buzc1HZZnlIVcyHlRx/mExER8bjSUsP8dZk8vSyN/MJixg9sz8QboqhdQ/+Mi5QnazPKPXr0MBs2bLCytnhXXl4e4eHhtssQj1EuxE1VycWmvblMWZTCZ7uP0qdtQ56Mi6FD0zq2y/KsqpILuTRBN6O8b98+W0uLh02fPt12CeJByoW4qey5OF5QxONLNjF89ofsPnyCWSO7MT++r5rkC6jsuZCKZa1RPnbsmK2lxcOWL19uuwTxIOVC3FTWXBhjWLpxLzfOWsncj3cyuncr3ntwILf3iMRxHNvleV5lzYXYoeEmERERj9iZnc8ji1NYlZ5NTIu6/OluP91b1rddlkiVZa1R1jns4mbWrFm2SxAPUi7ETWXKRUFRCS+u2MaLK7dRM6Qajw3vzN392hBSTe8gX6rKlAuxz1qjXFRUZGtp8TDNrosb5ULcVJZcfLD1EFMXp7Dz8Alu7XYFU26OpkndMNtlBa3KkgvxBmu7Xvh8PpOfn29lbfEuv9+PTmyUsykX4ibYc3HgWAHTlm7mjY37aBfhY9ptMXynQ4TtsoJesOdCykdZd724qHeUHccZCjwPhAB/NsY85fKYkcBjgAG+MMbceanFiIiIVHbFJaXMW72LWcu3cqqklEmDr+QnA9pRMzTEdmkicpYLNsqO44QAfwQGA1nAOsdxXjfGbD7jMR2Ah4BrjTFHHMdpcqHrNmrUqOxVS6U1ZswY2yWIBykX4iYYc7Fh9xGmLExh875jDLiyMdNu60LrRj7bZVUqwZgL8a6LeUe5N5BhjNkO4DjOfOA2YPMZj/kx8EdjzBEAY8zBC120Xr16l16tVHqxsbG2SxAPUi7ETTDl4uiJUzy9bAvz1+2maZ0wEn/Qg9iYZtrurRwEUy7E+y6mUW4BZJ7xdRbQ56zHXAngOM5HnB7PeMwYs+zsCzmOEw/EA4SEhOD3nx4VmThxItHR0UyYMAGA/v37M3PmTPr27QuAz+dj5cqVjB07lpSUFACSkpJITk5m3rx5ACQkJNC8eXMmTZoEwODBg5k8eTIDBw4EICIigmXLljFq1CgyMjIAWLhwIUlJSSxYsACAqVOnEhYWxsMPPwzA8OHDmTBhwjcvusjISBYtWkRcXBxZWVkAJCcnk5iYyJIlSwCYMWMGBQUFTJs2DYCRI0cyevRoRowYAUBUVBTz589n6NCh35xHv2LFCqZPn/7N3o+zZs1i3759zJw5Ezj903FsbCyjR48GICYmhrlz5zJgwAC+nvNes2YNCQkJrFq1CoDExERSU1OZPXs2APHx8fTr149x48YB0KtXL1588UX69OlDSUkJISEhrF27lvHjx7Nu3ToA5syZw+rVq3n55Zcr7HlKT08nNzdXz5PHn6eKfj2NHj2a6OhoPU8ef54q+vVkjOG6667z9PNkgEE/fIjk/bU4XlBMWOZavtPa4bq2fejVq1eVeJ4q+vWUmppKfHy8Xk8ef56gYu97ZXXBD/M5jvN9YIgx5t6vvr4b6G2MmXjGY5YCRcBIIBJYBcQYY45+23X1YT5xow9hiBvlQtx4PRdbDxxnysIUPtmZQ49W9Xky7io6X1HXdlmVntdzIXaU54f5soAzNz2OBPa6PGaNMaYI2OE4zhagA7Du2y5aq1atSyxVqoKYmBjbJYgHKRfixqu5OHGqmOffTefVVTsIDwvlqduvYqS/JdW0J3KF8GouJDhdzDvKocBWYBCwh9PN753GmE1nPGYoMNoYc4/jOBHAZ0B3Y8zhb7uu3+83+olPREQqk7c37efxJZvZc/QkI/2R/CY2moa+GrbLEqnyyvqOcrULPcAYUwzcB7wFpAILjDGbHMeZ5jjOrV897C3gsOM4m4H3gYTzNckAW7ZsudRapQoYMGCA7RLEg5QLceOlXGTmnODev64j/m/rCa8Zyr9/2o/ffa+bmmQLvJQLCX4XtY+yMeZN4M2z/mzqGf9tgElf/e+ilJaWXuxDpQrR3Lq4US7EjRdycaq4lFdWbWf2e+lUcxweHtaJcde2pXrIBd+HknLihVxI5WHtCGsREZFgtnrbYR5ZnELGwTyGdmnG1OGduaK+Pn8jUplYO8K6Z8+eZv369VbWFu8qLi4mNFQ/v8n/Ui7Eja1cZOcVMuONVF77bA8tG9Zi2q0xXN/pgudsSQXR/ULclNuMcnn5ev9AkTMlJCTYLkE8SLkQNxWdi5JSw9/X7OKG369gyca93Hd9FG//YoCaZI/R/UICydqPXHl5ebaWFg/7eqNzkTMpF+KmInORsieXyYtS+CLzKP3aNeKJuBiimoRX2Ppy8XS/kEDS7yZERES+xbGCIma9vZV5q3fS0FeT5+7ozm3dr9DR0yJVhLVGuVWrVraWFg9LTEy0XYJ4kHIhbsozF8YYlmzcxxNLN5OdV8jdfVvz4E0dqVerermtKYGh+4UEkrUZ5YKCAltLi4elpqbaLkE8SLkQN+WVi+2H8rj71U+4P+kzmtUNY/HPrmXabTFqkoOE7hcSSNZ2vfD5fEZ7HcrZ/H4/OrFRzqZciJtA56KgqITE9zN4aeV2alavxq+GdOTOPq0J0dHTQUX3C3FT1l0vNKMsIiJV3vtbDvLo4k3szjlBXPcrePjmaJrUCbNdlohYZq1RjoiIsLW0eFh8fLztEsSDlAtxE4hc7Ms9ybQlm0lO2U+7xj7+eW8fronSv0/BTPcLCSRrjXJ4uLbVkXP169fPdgniQcqFuLmcXBSXlDL34508u3wrxaWGhCEdubd/W2qGhgSwQrFB9wsJJM0oi6dotkzcKBfipqy5WL8rh8kLU0jbf5wbOjXh8Vu70LJh7XKoUGzQ/ULcaEZZRETkPHLyT/F0chr/+jST5vXCeOmungzp0lR7IovIt7LWKPt8PltLi4f16tXLdgniQcqFuLnYXJSWGhZ8mslTy9LIKyjmJ9e14/5BHfDV1HtFlZHuFxJI1kYv/H6/0a9GRESkPG3ee4wpi75kw+6j9G7TkCfiYujYrI7tskSkgpV19MLagSPaEFzc9OnTx3YJ4kHKhbg5Xy7yCot5Yulmhv/hQ3YdPsEz3+/Gv37SV01yFaD7hQSSfu8knlJSUmK7BPEg5ULcuOXCGMMbX54+evrg8ULu7N2KhCEdqV+7hoUKxQbdLySQ1CiLp4SEaGsmOZdyIW7OzsWO7HymLk5hVXo2MS3q8qe7/XRvWd9SdWKL7hcSSJpRFhGRoFZQVELiim28tGIbNUOr8cshHbmrr46eFpH/F3Qzyrt377a1tHjY+PHjbZcgHqRciJvx48fz/paDDHnuA154N53Yq5rx7i8HcM81bdQkV2G6X0ggWRu90GEj4mbdunW2SxAPUi7kbHuPnuTdgrYkz1mno6flf+h+IYGkGWUREQkaRSWlzPloB8+9k86piA4kDOnIj/u3o0aotV+QikglZq1RbtOmja2lxcPmzJljuwTxIOVCANbtzGHKwhS2HDjOoE5NGNUxlMH9omyXJR6j+4UEkrUfwfPy8mwtLR62evVq2yWIBykXVdvhvEJ++e8v+P5Lq8krLOaVMX5eHduLHV/qV+xyLt0vJJCs7Xrh8/mM5pTlbH6/H+2GImdTLqqm0lLD/HWZPL0sjfzCYn58XTsm3hBF7RqnfxmqXIgb5ULclHXXC80oi4iI56TsyWXKohQ+zzxK33YNeeK2GDo01al6IlKxrDXKTZo0sbW0eNjEiRNtlyAepFxUHccKipj19lbmrd5JQ18Nnr2jG3HdW+A45273plyIG+VCAslaoxwWFmZrafGw6Oho2yWIBykXlZ8xhte/2MuTb6SSnVfI3X1b8+BNHalXq/q3fo9yIW6UCwkkzSiLp2i2TNwoF5XbtkN5TF2cwkcZh+kaWY8n42LoGnnho6eVC3GjXIgbzSiLiEhQOXmqhD++n8GfPthGWPUQnoiL4c7erXSqnoh4hrVGOTw83NbS4mH9+/e3XYJ4kHJR+byXdoCpizeRdeQkt1/dgoeGRdO4Ts1LuoZyIW6UCwkka6MXPXv2NOvXr7eytnhXcXExoaH6RYf8L+Wi8thz9CSPv76JtzcfIKpJOE/GxdC3XaMyXUu5EDfKhbgp6+iFtQNH0tLSbC0tHta3b1/bJYgHKRfB71RxKS+u2MaNz6xkVXo2vx7aiTfv71/mJhmUC3GnXEgg6UcuEREpV2u2H+aRRSmkH8zjps5NmTq8M5ENatsuS0Tkgqw1ytWqWXszWzzM5/PZLkE8SLkITtl5hcx4M5XXNuwhskEtXr3Hz6DopgG7vnIhbpQLCSRrM8p+v99o+xYRkcqnpNTwz092M3NZGieLSvjJde352fVR1KoRYrs0Eamigm5GeefOnbaWFg8bO3as7RLEg5SL4PFlVi63J37EI4tSiGlRj+SfX8cvh3QslyZZuRA3yoUEkrXRi5MnT9paWjwsJSXFdgniQcqF9+WeLOKZt7fwtzW7iAivyfOjunNrtytcj54OFOVC3CgXEkj6MJ+IiJSZMYbFn58+ejonv5B7+rVh0k1XUjfs24+eFhEJFtYa5Xbt2tlaWjwsKSnJdgniQcqFN2UcPM6URSms2Z5Dt5b1mTuuFzEt6lXY+sqFuFEuJJCszSjn5ubaWlo8LDk52XYJ4kHKhbecPFXC08vSiH1+Fan7jjNjxFUsHH9NhTbJoFyIO+VCAsnarhc+n8/k5+dbWVu8y+/3o91Q5GzKhXcs33yAx17fxJ6jJ/lez0h+E9uJiPBLO3o6UJQLcaNciJuy7nqhGWUREbmgzJwTPL5kE++kHqRj0zos+Ek/erdtaLssEZFyZa1Rbtasma2lxcMSEhJslyAepFzYc6q4lFdWbWf2e+lUcxweHtaJcde2pXqI/UOjlAtxo1xIIFlrlKtX1yei5VzNmze3XYJ4kHJhx8cZ2TyyOIVth/IZ2qUZU4d35or6tWyX9Q3lQtwoFxJImlEWT9FsmbhRLirWgWMFPPlGKku+2EurhrV5/NYuXN+pie2yzqFciBvlQtxoRllERC5LcUkpcz/eyXPvpHOqpJSfD+rA+IHtCauuo6dFpGqy1ijXrVvX1tLiYYMHD7ZdgniQclH+1u3M4ZFFKaTtP87Ajo15/NYutG7ks13WeSkX4ka5kECyNnrRo0cPs2HDBitri3fl5eURHh5uuwzxGOWi/GTnFfLbN9P474YsWtSvxdThnbmpc9NyPXo6UJQLcaNciJuyjl5Y+9jyli1bbC0tHjZw4EDbJYgHKReBV1Jq+Nvqndzw+xW8/sUeJgxsz/JJ1zGkS7OgaJJBuRB3yoUEkmaURUSqmM8zj/LIohS+3JPLtVGNePzWGKKa6B04EZGzWWuUQ0PVo8u5IiIibJcgHqRcBMaR/FP87q0tzF+3myZ1ajJ79NXc0rV50LyDfDblQtwoFxJI1maU/X6/0fYtIiLlr7TUsODTTJ5elsaxgmLGXdOGXwy+kvCaesNCRKqGoJtR3r59u62lxcNGjRpluwTxIOWi7FL25PLdlz7mN699SYcmdXjj/u8w5ZbOlaJJVi7EjXIhgWTtTllYWGhrafGwjIwM2yWIBykXly73ZBGz3t7C39bsoqGvBs98vxu392gRtGMWbpQLcaNcSCAF/1sKIiLyDWMMCz/bw4w308jJL+Suvq158KaO1KtV3XZpIiJBx1qj3L59e1tLi4ctXLjQdgniQcrFxdmy/ziPLE7hkx05dG9Zn7njehHTop7tssqNciFulAsJJGszyjk5ObaWFg9LSkqyXYJ4kHJxfnmFxUx/YzPDXljF1gPHeer2q3ht/DWVukkG5ULcKRcSSNYa5SNHjthaWjxswYIFtksQD1Iu3BljWLpxL4OeWcErq3Yw0h/J+w8OZFTvVlSrVnlmkb+NciFulAsJJM0oi4gEoW2H8nh08SY+zMimyxV1eemunlzdqoHtskREKhVrjXLz5s1tLS0eNnXqVNsliAcpF//v5KkSZr+XziurthNWPYRpt3XhB31aE1IF3kE+m3IhbpQLCaSLapQdxxkKPA+EAH82xjz1LY/7HvBvoJcx5ryniVSrZm3qQzwsLCzMdgniQcrF6TGLtzcfYNqSzew5epLbe7TgodhoGtepabs0a5QLcaNcSCBd8GQ+x3FCgK3AYCALWAeMNsZsPutxdYA3gBrAfRdqlH0+n8nPz7+M0qUy8vv96MRGOVtVz8Xuwyd49PUU3t9yiI5N6/BEXAy92za0XZZ1VT0X4k65EDdlPZnvYt5R7g1kGGO2f7XQfOA2YPNZj3sC+B3wy0stQkREzlVQVMJLK7eRuGIb1as5TLk5mnuuaUP1EP1GTkSkIlxMo9wCyDzj6yygz5kPcBznaqClMWap4zjf2ig7jhMPxAPUqlULv/90Yz9x4kSio6OZMGECAP3792fmzJn07dsXAJ/Px8qVKxk7diwpKSnA6e1fkpOTmTdvHgAJCQk0b96cSZMmATB48GAmT57MwIEDAYiIiGDZsmWMGjXqm1N7Fi5cSFJS0jefkJ06dSphYWE8/PDDAAwfPpwJEyYQGxsLQGRkJIsWLSIuLo6srCwAkpOTSUxMZMmSJQDMmDGDgoICpk2bBsDIkSMZPXo0I0aMACAqKor58+czdOhQsrOzAVixYgXTp09n+fLlAMyaNYt9+/Yxc+ZMAMaMGUNsbCyjR48GICYmhrlz5zJgwAC+fld+zZo1JCQksGrVKgASExNJTU1l9uzZAMTHx9OvXz/GjRsHQK9evXjxxRfp06cPJSUlhISEsHbtWsaPH8+6desAmDNnDqtXr+bll1+usOfp6NGjAHqePP48VfTrae/evfj9/ir1PFVrEUOt79zDrsMnqLH/S2qkv8WAu1/mxT/+wbPPU0W/noYNG8YDDzyg15PHn6eKfj3t3buXBx54QM+Tx58nqNjXU1ldzOjF94Ehxph7v/r6bqC3MWbiV19XA94DxhpjdjqOswL45YVGL7p3724+//zzMhculdOhQ4do3Lix7TLEY6pSLvYcPcm0JZt4a9MB2jX28cRtMVwbFWG7LE+qSrmQi6dciJuyjl5czO/vsoCWZ3wdCew94+s6QAywwnGcnUBf4HXHcc5bTHp6+qVVKlXC1z91i5ypKuSisLiExBUZ3PjMSj7Yms2vhnZk2c+vU5N8HlUhF3LplAsJpIsZvVgHdHAcpy2wBxgF3Pn1XxpjcoFv7uQX+46yiIictnLrIR5/fRPbs/O5qXNTpg7vTGSD2rbLEhGp8i7YKBtjih3HuQ94i9Pbw/3FGLPJcZxpwKfGmNfLsnD16tXL8m1SyUVGRtouQTyosuYi68gJnli6mbc2HaBthI+543oxsGMT22UFjcqaC7k8yoUE0gVnlMuL3+832r5FRKqigqISXvlgO39ckYGDw303RHFv/7bUDA2xXZqISKVUnjPK5eLrTyKKnCkuLs52CeJBlSkX76cdZMhzH/DM8q1c37EJ7zw4gJ9dH6UmuQwqUy4kcJQLCSRrR1gXFRXZWlo87OvtckTOVBlykZlzgseXbOad1NO7WfztR73p30GfzL8clSEXEnjKhQSStUZZRKQq+PrQkBdXbCOkmsNvYjvxw2vbUiNUh4aIiHidtRll7aMsbrT/pbgJ1ly8s/kAjy/dRGbOSW7p2pzJN0fTvF4t22VVGsGaCylfyoW4CboZ5YMHD9paWjwsMTHRdgniQcGWi12H8/nh3HXcO+9TaoaG8M97+/CHO3uoSQ6wYMuFVAzlQgLJWqOcm5tra2nxsK+P8BQ5U7Dk4uSpEma9vYXBz37A2u2HmTwsmuSf9+caHRpSLoIlF1KxlAsJJM0oi4hcJmMMb28+wLQlm9lz9CS3db+Ch4dF07RumO3SRETkMlhrlFu0aGFrafGwGTNm2C5BPMjLudiRnc9jr29i5dZDdGxah/nxfenbrpHtsqoEL+dC7FEuJJCsNcqlpaW2lhYPKygosF2CeJAXc3HiVDF/fD+DVz7YQY3QajxyS2fG9GtN9RDtZlFRvJgLsU+5kECytuuFz+cz+fn5VtYW7/L7/ejERjmbl3JhjGFZyn6eWLqZvbkF3H51C34zrBNN6mjMoqJ5KRfiHcqFuCnrrheaURYRuUgZB/N4fMkmVqVn06lZHZ4ffTW92jS0XZaIiJQTa41ygwYNbC0tHjZy5EjbJYgH2c5FfmExL7yXzl8+3EFY9RAeG96Zu/q2JlRjFlbZzoV4k3IhgWRt9KJr165m48aNVtYW78rMzKRly5a2yxCPsZULYwxLN+5j+hup7D9WwPd6RvLroZ1oXKdmhdci59L9QtwoF+Im6A4c2bZtm62lxcNGjBhhuwTxIBu5SD9wnB/8eS0Tkz6joa8G/x3fj99/v5uaZA/R/ULcKBcSSJpRFhE5Q+7JIp5/J52/rt6Jr0YIT9zWhTv7tCakmmO7NBERqWDWGuWaNfWujJwrKirKdgniQRWRi9JSw3/WZ/G7t9I4nH+KUb1a8cubrqRRuO5VXqX7hbhRLiSQrM0o+/1+o+1bRMQLPs88yqOLU/giK5cererz+K0xXBVZz3ZZIiISIEE3o5yenm5rafGwoUOH2i5BPKi8cnHoeCEJ//6CuD9+xN7cAmaN7MZ/x1+jJjlI6H4hbpQLCSRroxfFxcW2lhYPy87Otl2CeFCgc1FUUspfP97J8++kU1Bcwk+ua8fEQR0Ir6mPbQQT3S/EjXIhgaR/FUSkSvkwPZvHlmwi42Ae113ZmEeHd6Z943DbZYmIiAdZm1Hu0aOH2bBhg5W1xbvy8vIID1fTIv8rELnIzDnB9DdSWbZpP60a1uaRWzpzY3QTHEe7WQQr3S/EjXIhboJuRnnfvn22lhYPmz59uu0SxIMuJxcFRSU8u3wrN85aycqth/jlTVfy9gPXMbhzUzXJQU73C3GjXEggWWuUjx07Zmtp8bDly5fbLkE8qCy5MMaQ/OU+Bj2zkuffTWdw56a8++AA7ruhA2HVQ8qhSqloul+IG+VCAkkzyiJS6Ww9cJzHl2zio4zDdGpWh/nxfenbrpHtskREJMhYa5R1Dru4mTVrlu0SxIMuNhdnn6r3+K1d+EGfVoSGWPvlmZQj3S/EjXIhgWStUS4qKrK1tHiYZtfFzYVy4XaqXsKQjjT01aigCsUG3S/EjXIhgWRt1wufz2fy8/OtrC3e5ff70YmNcrbz5eLMU/V6tm7A47d2IaaFDgypCnS/EDfKhbgp664XmlEWkaB08HgBM5dt4d/rs2hSpybP3tGNuO4ttJOFiIgEjLVGuVEjfbBGzjVmzBjbJYgHnZmLwuIS5ny0kz+8l0GhTtWr0nS/EDfKhQSStX9Z6tXTr0blXLGxsbZLEA+KjY3FGMPyzQeY/mYquw6f4MboJky+uTNtI3y2yxNLdL8QN8qFBJJmlMVTNFsmbroPGEbM3Y/yYUY2UU3CmXpLZ667srHtssQy3S/EjXIhbjSjLCKVztETp3h2+VaO9vkpG7OO8tjwzvygb2uqa7s3ERGpANYa5Vq1atlaWjwsJibGdgniAcUlpfxj7W5mLd/K8YIirji5g6WPjtd2b/I/dL8QN8qFBJK10Qu/32/0qxEROduH6dlMW7qJrQfyuKZ9I6YO70ynZnVtlyUiIkGsrKMX1n5/uWXLFltLi4cNGDDAdgliyc7sfO7966fc9epaCopK+dPdPfnHvX3o1KyuciGulAtxo1xIIFkbvSgtLbW1tHiYPuBZ9RwvKOIP72fwlw93UCOkGr8a2pEfXtuWsOoh3zxGuRA3yoW4US4kkPRhPhGx4v+Pnd5Cdl4h3+sZya+GdKRJ3TDbpYmIiAAWZ5R79uxp1q9fb2Vt8a7i4mJCQ/XzW2X36c4cHl+ymS/35HJ1q/o8NrwL3VrW/9bHKxfiRrkQN8qFuAm6GeWsrCxbS4uHJSQk2C5BytHeoye5P+kzvvfSag4dL+S5O7rz2vhrztskg3Ih7pQLcaNcSCBZ+5ErLy/P1tLiYatWrbJdgpSDE6eKefmD7by0chvGwP03RPHTge2pXePibkHKhbhRLsSNciGBpN9NiEi5KS01LPxsDzPf2sL+YwUMu6oZD8VG07JhbduliYiIXJC1RrlVq1a2lhYPS0xMtF2CBMgnO3J4YunpOeSukfWYfefV9GrTsEzXUi7EjXIhbpQLCSRrM8oFBQW2lhYPS01NtV2CXKbdh08w/u/rGfmn03PIs0Z2Y9GEa8vcJINyIe6UC3GjXEggWdv1wufzGe11KGfz+/3oxMbgdKygiD+8l8Hcj3YSUs1h/MD2/Lh/O2rVCLnwN1+AciFulAtxo1yIm7LueqEZZRG5LMUlpSR9sptn30nnyIlTfLdHJAlDOtJU+yGLiEiQs9YoR0RE2FpaPCw+Pt52CXIJ3t9ykOlvpJJxMI8+bRvyyC2diWlRL+DrKBfiRrkQN8qFBJK1Rjk8PNzW0uJh/fr1s12CXIStB47z5BupfLD1EG0a1eZPd/fkps5NcRynXNZTLsSNciFulAsJJM0oi6dotszbsvMKeXb5VpI+2U14zVDuH9SBMf3aUCO0fD8XrFyIG+VC3CgX4kYzyiJSbgqLS5jz0U7++F4GJ4pKGNOvDT8f1IEGvhq2SxMRESk31hpln89na2nxsF69etkuQc5gjCE5ZT+/TU4lM+ckgzo14aFh0UQ1qdjRKeVC3CgX4ka5kECyNnrh9/uNfjUi4l3rd+Uw48001u86QqdmdZhyc2e+00EfwhURkeBT1tELBDV9NAAAHVpJREFUaweOaENwcdOnTx/bJVR5O7LzGf/39Xz3xdVk5pzgqduv4o37+1ttkpULcaNciBvlQgJJM8riKSUlJbZLqLJy8k/xwrvp/H3NLmqEVuOBG6/kx9e1pXYN+7cJ5ULcKBfiRrmQQLL/L6DIGUJCLv8UN7k0BUUl/OWjHbz4/jZOFJVwR6+W/OLGDjSp450DQ5QLcaNciBvlQgJJM8oiVVRpqWHhZ3t45u0t7M0t4MboJvwmthNRTerYLk1ERCSggm5Geffu3baWFg8bP3687RKqhA/Ts7ll9oc8+O8viKhTk/nxffnzPb082yQrF+JGuRA3yoUEkrXRCx02Im7WrVtnu4RKLW3/MX77Zhortx4iskEtnh/VneFdr6BatfI5US9QlAtxo1yIG+VCAkkzyiJVwP7cAmYt38J/1mcRXjOUycOiGXNNa2qGapZPRETk21hrlNu0aWNrafGwOXPm2C6hUskrLOZPK7fxyqrtlJbCD69ty303RFG/dnCdqKdciBvlQtwoFxJI1maU8/LybC0tHrZ69WrbJVQKp4pLmbd6JwNnvs/s9zIY3LkZ7z44gCm3dA66JhmUC3GnXIgb5UICydquFz6fz2hOWc7m9/vRbihlV1pqWLJxL8+8vZXdOSfo07YhDw2LpnvL+rZLuyzKhbhRLsSNciFuyrrrxUWNXjiOMxR4HggB/myMeeqsv58E3AsUA4eAHxpjdl1qMSJSNsYYPkjP5nfL0ti09xjRzesyZ1wvBl7ZGMfx9gf1REREvOqCjbLjOCHAH4HBQBawznGc140xm8942GeA3xhzwnGc8cDvgDvOd90mTZqUvWqptCZOnGi7hKDzReZRnl6WxsfbDtOyYS2eu6M7t3bz/k4Wl0K5EDfKhbhRLiSQLuYd5d5AhjFmO4DjOPOB24BvGmVjzPtnPH4NcNeFLhoW5p1Tv8Q7oqOjbZcQNLYfyuP3b2/hzS/308hXg8eGd+bOPq2pEWrtowflRrkQN8qFuFEuJJAuplFuAWSe8XUW0Oc8j/8RkOz2F47jxAPxcPqISb//9KjIxIkTiY6OZsKECQD079+fmTNn0rdvXwB8Ph8rV65k7NixpKSkAJCUlERycjLz5s0DICEhgebNmzNp0iQABg8ezOTJkxk4cCAAERERLFu2jFGjRpGRkQHAwoULSUpKYsGCBQBMnTqVsLAwHn74YQCGDx/OhAkTiI2NBSAyMpJFixYRFxdHVlYWAMnJySQmJrJkyRIAZsyYQUFBAdOmTQNg5MiRjB49mhEjRgAQFRXF/PnzGTp0KNnZ2QCsWLGC6dOns3z5cgD+r707j6+6vvM9/vpmgYQEEiAsISSyBUWRsgQBWxarKNjr0hZZasXtlrkwbac6ZWZaZ5wZZ9Q63LGdaw0ttyraxy2OtQ8V+wC5tCouZUeUHcIeSAhLCEkgZDnf+eMEJ8Wv5BDPOd9fkvfz8eCR7Zfz/ejvfZJPzvmc7+/pp5+mpKSEBQsWADB79mymTp3KrFmzABg6dCiLFy9m4sSJn+5HvWbNGubPn8/7778PQGFhITt27OCZZ54BYM6cOYwbN477778fgNGjR7Nw4ULGjBlDQ0MDiYmJrF27lrlz5366B+ULL7zA6tWrWbRoUdzO0549e6ioqNB5usR5WrpyFdP+8XkOJ/fFhOqZNaIHfSp38vO/+hE/j9N5ivf9adasWQwZMqRVnacg3J/a+s89ay0TJkzQeQr4eYr3/WnHjh3MmTNH5yng5wnie39qqWZfzGeMuQu4xVr7Pxs/vge4zlr7mec2jDHfBr4LTLTWnr/U7erFfOKiF2F8vopzdfxi1V5e+HA/DSHL3WOu4LtfHURWekffpcWcciEuyoW4KBfiEssX8xUDuU0+7gscdRRwE/AIETTJAOnp6ZHWKO3I+PHjfZcQODV1Dbz4pwMUvruXMzV13PGlPjw8+UryunfyXVrcKBfiolyIi3Ih0RTJI8pJwG7gRuAIsB74lrV2W5NjRgCvAlOstXsiWXjUqFF248aNLa1b2qj6+nqSknTBSID6hhC/21TMz/6wh5KKGiYO7sHfTLmSa/pk+C4t7pQLcVEuxEW5EJeWPqLc7Kt+rLX1hMcpVgA7gFestduMMY8ZY25vPGwBkA781hiz2RiztLnb3blz5+XWKu3AhTml9iwUsryx+QiTf/oef/u7LfTsksKS74zlxQeua5dNMigX4qZciItyIdEU0Z9c1tplwLKLPvdok/dvinJdIu2OtZaV24/x9Mrd7Cyt5MpenVl0zygmX91LeyGLiIh44O25iYSEtreFlXxxaWlpvkuIO2stHxSd4H///918fPg0/bPS+I+Zw7ltWNvaC/mLaI+5kOYpF+KiXEg0ebuEdUFBgdWrUqW923DgFAtW7GLt/lPkZKby/RsH8c2RfUlK1B+SIiIi0RKzGeVYOXDggK+lJcDuu+8+3yXExdYjFdz3wjqm/WI1e49X88+3X8PbP5zIjNF5apId2ksu5PIoF+KiXEg0eRu9OHfunK+lJcAubDDeVu0+VslPV+5m+dZSMjsl83dTr+Lecf1I7ZDou7RAa+u5kJZRLsRFuZBo0v4pInFw8GQ1P/vDHl7ffIS0Dkn81Y35PDi+P11Skn2XJiIiIp/DW6M8YMAAX0tLgC1ZssR3CVF1+NRZnn2niFc3FpOUaJgzfgB/MXEg3dI6+C6tVWlruZDoUC7ERbmQaPLWKFdUVPhaWgJs+fLl5Ofn+y7jCysuDzfIv91QTIIx3D0mj7+8YRA9u6T4Lq1Vaiu5kOhSLsRFuZBo8rbrRVpamq2urvaytgRXQUEBrXk3lHCDvJdXNx7GYJh5XS5zJw0kOyPVd2mtWmvPhcSGciEuyoW4tHTXC80oi0TBkdPnGh9BbmyQR+cx7wY1yCIiIq2Zt0a5d+/evpaWAJs/f77vEi7L0cYG+ZUNhwGYMTqXeZMG0SdTDXI0tbZcSHwoF+KiXEg0eWuUk5P1an/5rOzsbN8lRKSk4hyF7+zlP9cfxmKZMTqXuZMGkaMGOSZaSy4kvpQLcVEuJJo0oyyBEvTZspKKcyx8dy8vrws3yNMLcpl3gxrkWAt6LsQP5UJclAtx0YyySAwdOnmWhavC27wB3FWQy7xJA+nbtZPnykRERCRWvDXKXbp08bW0BNjkyZN9l/BnisoqKXxnL298fJTEBMOs6/KYM2GAGuQ4C1ouJBiUC3FRLiSavI1ejBw50m7atMnL2hJcVVVVpKen+y6DbUcrePadIpZvLSUlKZFvj83jO+MHaB9kT4KSCwkW5UJclAtxaenohWaUJVB8z5ZtOlTOs28X8cedZXTumMS91/fjga/015X0PPOdCwkm5UJclAtx0YyySAtZa1mz7xQ/f2cPHxadpGunZH5482DuGdePjFTtziIiItJeeWuUk5LUo8tnZWVlxW0tay2rdh/n528XseFgOVnpHXnk1iF8a0weaR2VzyCJZy6k9VAuxEW5kGjyNnpRUFBg9dSI+FDfEGLZ1lJ+uWov246eoU9GCv9r0kCmF+SSkpzouzwRERGJspaOXiTEophI7Nu3z9fSEmAzZ86M2W2fq23gpdUHuOHf3+X7Sz7iXF0DT33zWt6dfwOzx/VTkxxgscyFtF7KhbgoFxJN3p5fPn/+vK+lJcCKioqifpunz9by0uqDLP7TAU5V1zIiL5O//9rVTB7Si4QEE/X1JPpikQtp/ZQLcVEuJJo0iClt1pHT53ju/f28vP4QZ2sb+OpVPfmLCQO4rn83jFGDLCIiIpfmrVEeOHCgr6UlwF577bUvfBu7Siv55aq9LP34KAC3f6kPcyYO4KreushNaxWNXEjbo1yIi3Ih0eRtRvnUqVO+lpYAW7JkSYu+z1rLn/ae4IHF67nlZ+/x1rZSZo/rx6q/uYGnZwxXk9zKtTQX0rYpF+KiXEg06YIjEiiXu1F8bX2I339ylF+9v5/tJWfontaBe6/vx+xxV5DZSRcJaSt0AQFxUS7ERbkQF11wRNqV8upafrPuEC/+6QBllefJ75nOT75xLXeOyNHuFSIiIhIV3hrl7OxsX0tLgD366KOX/Pre41U8/8F+frepmJq6EBMG92DBXf2ZkJ+lF+i1Yc3lQton5UJclAuJJm+NckKCt/FoCbCUlJTPfC48f3yS5z7Yz9s7y+iQlMA3RuTwwFf6M7hXZw9VSry5ciGiXIiLciHRpBllCZSms2U1dQ28+fFRnvtgPztLK8lK78A9Y/tx99g8stI7eq5U4kkzh+KiXIiLciEumlGWNuPI6XP8vzUHeXn9YU5V13Jlr87827Rh3P6lPpo/FhERkbjx1ihnZGT4WloCyFrL6r0n6XTzDxj/1NsATL66F/eO68e4gd01f9zO3Xbbbb5LkABSLsRFuZBo8jZ6MXz4cLt582Yva0twVJ2v57VNxby4+iBFZVVkpibxrTFXcPfYK8jJTPVdngTE8ePH6dGjh+8yJGCUC3FRLsSlpaMXmlEWL/Yer+LXqw/yu43FVJ6vZ1jfDO4d148n5k5j47o1vsuTgNHMobgoF+KiXIiLZpQl8GrrQ6zcfozfrDvIh0UnSU40/I9hfZg97gqG52ZijOHJUL3vMkVEREQAj41ycnKyr6Ulzg6erGbJusO8uvEwJ6pqyclM5Yc3D2bG6Dx6dP7z3Sv69u3rqUoJMuVCXJQLcVEuJJq8jV4UFBRYPTXSdl149HjJukN8UHSCxATDjVf1ZNaYPCbk9yAxQS/OExERkfho6eiFt6t+FBUV+VpaYujgyWqeemsn1//kj/zlbzax/0Q1fz15MB/+7VdZNLuAG67seckm+c4774xjtdJaKBfiolyIi3Ih0eRt9KKurs7X0hJlNXUNrNhWym83FPNB0QkSDNw4pBffui6PCYMv79Hj4uLiGFYqrZVyIS7KhbgoFxJNejGftIi1lo8On+bVjcW8+fFRKmvqyclM5aGbBjNjdC69M3QJUREREWndtI+yXJZjZ2p47aMjvLqxmKKyKlKSE7h1aDbTCvoytn93Er7g7LH2vxQX5UJclAtxUS7EpdXNKJeVlflaWi7T+foGlm0p4f4X1jHuyT/yk+U7yUxN5qlvXsv6R27i6RnDuX5g1hdukgEKCwujULG0NcqFuCgX4qJcSDTpgiPiFApZNh4q543NR/j9JyWcPltH7y4pfHNUDtNG5dI/Ky0m62qjeHFRLsRFuRAX5UJcdMERiYpdpZW8vvkISzcf5cjpc6QkJ3Dz1b2ZNqovXx6UpW3dREREpN3w1ijn5OT4WloucuT0OZZuPsobm4+ws7SSxATD+PwsfnjLYG6+ujdpHeMXkyeeeCJua0nroVyIi3IhLsqFRJO3RjkUCvlaWoCyyhpWbDvGm5uPsu7AKQBG5mXy2B3XcOu12WSld2zmFmKjpqbGy7oSbMqFuCgX4qJcSDRpRrkdKTtTw/KtpSzbUsK6A6ewFgb1TOfO4X24/Us55HXv5LtEzZaJk3IhLsqFuCgX4qIZZXEqrahh+dYSlm0pYcPBcqyFwb3S+f5X8/nasGwG9+rsu0QRERGRQPLWKHft2tXX0m3egRPV/GHHMZZvLWXjwXIArurdmR/cOJivDevNoJ7BbY6nT5/uuwQJIOVCXJQLcVEuJJq8jV4MGzbMfvLJJ17WbmsaQpbNh0/zhx3HWLn9GEVlVUC4Of7atdncOiybgT3SPVcZmcOHD5Obm+u7DAkY5UJclAtxUS7EpaWjF5pRbqXO1TbwQdEJVm4v5e2dZZyoqiUxwTCmfzcmX92Lm4b0Ireb/5njy6XZMnFRLsRFuRAX5UJcNKPcDhw4Uc17e46zatdxPig6wfn6EJ07JjHpqp7cNKQnkwb3JKNTsu8yRURERNoEb41yx45+th9rTarO17N670ne232c9/Yc5+DJswDkdevEzNG5TL66N9f170aHJG9XIo+6QYMG+S5BAki5EBflQlyUC4kmb6MXBQUFVk+N/LlQyLKj9Azv7T7Bqt1lbDxYTl2DJTU5kesHdmfC4B5MHNyDfjG6fLSIiIhIW9TqRi/27Nnja+nACIUsO0srWbPvJGv2nWTt/lNUnKsDwi/Ee+Ar/ZmY34NR/brSMSnRc7XxMWXKFN566y3fZUjAKBfiolyIi3Ih0eStUa6vr/e1tDehkGXXsT9vjE+fDTfGed06ccs1vRg7oDtfHpRFry4pnqv148SJE75LkABSLsRFuRAX5UKiSS/mi6Hq8/V8XHyajw6dZtPBcjYdKqe8sTHO7ZbK5CHhxnjswO7kZKZ6rlZEREREmvI2ozxy5Ei7adMmL2vHgrWWQ6fOsulQOZsOnmbToXJ2llbSEAr//x3YI42ReV0ZM6A7Y/p3a5Vbt8VDVVUV6emtY89niR/lQlyUC3FRLsSl1c0ol5SU+Fr6CwuFLAdOVrPlSAXbjp5hS3EFW49WUFkTHidJ65DI8LxM5k0ayMi8rozIyySzUwfPVbcOjz/+OE8++aTvMiRglAtxUS7ERbmQaNIFR5pxpqaOPccq2X2sil2llWwvOcP2o2eoOh9uijskJTCkd2euyclgaJ8MhudmcmXvziQmGM+Vt07aKF5clAtxUS7ERbkQl1b3iHKQWGspP1vHgZPV7C2rYndjY7z7WCUlFTWfHpeanMhV2Z35xsgchvbJYGhOBvm90klObDv7GIuIiIhIWESNsjFmCvAfQCLwK2vtTy76ekfgJWAUcBKYYa09cKnbjPd12OsbQpRVnufI6XMcOFHNwZNnOXDyv99eGJuA8KPE+T3TGTugO/m90rmyV2cG9+pMTmYqCXqkOKaefvpp3yVIACkX4qJciItyIdHUbKNsjEkEngUmA8XAemPMUmvt9iaHPQiUW2sHGWNmAk8BMy51u3V1dS2vuona+hDlZ2s5VV1LeXUtJ6prKa04R0lFDaUVNZRU1FBScY7jlecJNZkySUww5HZNJa97GiPyMrmiexr9uneif1YaV3RP0+iEJ615dl1iR7kQF+VCXJQLiaZmZ5SNMeOAf7LW3tL48Y8ArLVPNjlmReMxq40xSUAp0MNe4sY7delqV2zcTX3IEgpZ6kOWhsa3NbUNVNfWc7a2gerz4bdna+upbvy4/Gwd5Y2NceV5937MaR0Syc5MJTsjhd5dUsJvM1Lpk5lC/6w0+mSmamQigDRbJi7KhbgoF+KiXIhLLGeUc4DDTT4uBsZ83jHW2npjTAXQHfjcXb8b0rK457l1zS6emGDo1CGRtA5JdOoYfpvZKZn+3TvRNa0D3Tp1CL9t/Nc9rQO9M1LonJIcwX+aiIiIiIhbJI2yawbh4keKIzkGY8wcYA5AaudMMtb/CqxlxvS7GDSgP088/q9gQ4wZNYJ//ed/YMqNN4BtID0tjVWrVnHfffexdetWSoAlS5awfPlyXnzpJQDmz59PenY2sx5+GIDJkyfzyCOPMGnSJACysrJ46623mDlzJkVFRQC89tprLFmyhFdeeQWARx99lJSUFH784x8DcNtttzFv3jymTp0KQN++fXn99de58847KS4uBmD58uUUFhby5ptvAvDEE09QU1PDY489BsD06dOZNWsWX//61wEYNGgQL7/8MlOmTPn06kHvvvsujz/+OCtXrgTC81UlJSUsWLAAgNmzZzN16lRmzZoFwNChQ1m8eDETJ07kws4ha9asYf78+bz//vsAFBYWsmPHDp555hkA5syZw7hx47j//vsBGD16NAsXLmTMmDE0NDSQmJjI2rVrmTt3LuvXrwfghRdeYPXq1SxatAiA733vewwZMoR58+YBMH78eBYsWMDYsWMBSLvoPNHkPL3U5DxlZ2fz8Oecpwv/PTpPwT5P8b4/lZWVUVBQoPMU8PMU7/vT3XffzUMPPaTzFPDzFO/7U1lZGQ899JDOU8DPE8T3/tRS3kYvrr32Wrtly5YWFy5t0549e8jPz/ddhgSMciEuyoW4KBfi0tLRi0iGdNcD+caY/saYDsBMYOlFxywF7m18fxrw9qWaZIB9+/Zdbq3SDlz4a1ekKeVCXJQLcVEuJJqaHb1onDn+LrCC8PZwz1trtxljHgM2WGuXAs8BvzbGFAGnCDfTIiIiIiKtVkT7KFtrlwHLLvrco03erwHuupyFU1NTL+dwaSeGDh3quwQJIOVCXJQLcVEuJJq8XcK6oKDAavsWEREREYm1WM4ox8SuXbt8LS0BNnHiRN8lSAApF+KiXIiLciHR5K1RDoVCvpaWALuwTY1IU8qFuCgX4qJcSDTp0nQiIiIiIg7eZpRHjRplN27c6GVtCa76+nqSkiJ6jam0I8qFuCgX4qJciEurm1G+cEUakabmz5/vuwQJIOVCXJQLcVEuJJq8NcpVVVW+lpYAu3DpTJGmlAtxUS7ERbmQaNKMsoiIiIiIg7dGOS8vz9fSEmCFhYW+S5AAUi7ERbkQF+VCoslbo1xTU+NraQmwHTt2+C5BAki5EBflQlyUC4kmb41yWVmZr6UlwJ555hnfJUgAKRfiolyIi3Ih0aQZZRERERERB2/7KBtjKgFdx1oulgWc8F2EBI5yIS7KhbgoF+JypbW28+V+k88duXe1ZONnaduMMRuUC7mYciEuyoW4KBfiYozZ0JLv0+iFiIiIiIiDGmUREREREQefjfIij2tLcCkX4qJciItyIS7Khbi0KBfeXswnIiIiIhJkGr0QEREREXFQoywiIiIi4hDzRtkYM8UYs8sYU2SM+TvH1zsaY/6z8etrjTH9Yl2T+BdBLh42xmw3xnxijPmjMeYKH3VKfDWXiybHTTPGWGOMtoBqByLJhTFmeuPPjG3GmN/Eu0aJvwh+j+QZY94xxnzU+LvkVh91SvwYY543xpQZY7Z+zteNMeb/NGbmE2PMyOZuM6aNsjEmEXgWmApcDcwyxlx90WEPAuXW2kHAT4GnYlmT+BdhLj4CCqy1w4BXgX+Lb5USbxHmAmNMZ+D7wNr4Vig+RJILY0w+8CPgy9baa4AfxL1QiasIf178PfCKtXYEMBMojG+V4sFiYMolvj4VyG/8NwdY2NwNxvoR5euAImvtPmttLfAycMdFx9wBvNj4/qvAjcYYE+O6xK9mc2Gtfcdae7bxwzVA3zjXKPEXyc8LgH8h/IdTTTyLE28iycV3gGetteUA1tqyONco8RdJLizQpfH9DOBoHOsTD6y17wGnLnHIHcBLNmwNkGmMyb7Ubca6Uc4BDjf5uLjxc85jrLX1QAXQPcZ1iV+R5KKpB4HlMa1IgqDZXBhjRgC51trfx7Mw8SqSnxeDgcHGmA+NMWuMMZd6REnahkhy8U/At40xxcAy4HvxKU0C7HL7j5hfwtr1yPDF+9FFcoy0LRGfc2PMt4ECYGJMK5IguGQujDEJhMez7otXQRIIkfy8SCL8VOokws8+vW+MGWqtPR3j2sSfSHIxC1hsrf13Y8w44NeNuQjFvjwJqMvuOWP9iHIxkNvk47589qmPT48xxiQRfnrkUg+bS+sXSS4wxtwEPALcbq09H6faxJ/mctEZGAq8a4w5AIwFluoFfW1epL9H3rDW1llr9wO7CDfO0nZFkosHgVcArLWrgRQgKy7VSVBF1H80FetGeT2Qb4zpb4zpQHiYfulFxywF7m18fxrwttVVUNq6ZnPR+BT7Lwk3yZo3bB8umQtrbYW1Nsta289a24/w7Prt1toNfsqVOInk98jrwA0AxpgswqMY++JapcRbJLk4BNwIYIwZQrhRPh7XKiVolgKzG3e/GAtUWGtLLvUNMR29sNbWG2O+C6wAEoHnrbXbjDGPARustUuB5wg/HVJE+JHkmbGsSfyLMBcLgHTgt42v7Txkrb3dW9EScxHmQtqZCHOxArjZGLMdaADmW2tP+qtaYi3CXPw18H+NMQ8Rfnr9Pj0Q17YZY5YQHsHKapxN/0cgGcBa+wvCs+q3AkXAWeD+Zm9TmRERERER+SxdmU9ERERExEGNsoiIiIiIgxplEREREREHNcoiIiIiIg5qlEVEREREHNQoi4iIiIg4qFEWEREREXH4Lw7N+vIpbNb+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "coord[:,2] = displacement\n",
    "fig, ax = plt.subplots(figsize = (12, 7))\n",
    "\n",
    "plt.plot(coord[:,1],coord[:,2])\n",
    "plt.xlim([0,1])\n",
    "plt.grid(b=True, which='major', color='k', linestyle='--')\n",
    "fig.savefig('fem_lista_2_Q4.png',dpi=300)                              # save figure as png (dpi = number of pixels)"
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}