FEM 2D Euler Beam.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": [
    "# Input data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 1                             # Beam length in meters    \n",
    "q = 100                           # load\n",
    "nel = 1024                        # 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": 3,
   "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": 4,
   "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",
    "    if inci[i,4] <= 0.5:       # inertia of beam section\n",
    "        inci[i,5] = 2*I\n",
    "    else:\n",
    "        inci[i,5] = I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boundary conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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],[nnos,1,0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mask stiffness matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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": 7,
   "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": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#   F = [node | degree of freedom | value]\n",
    "#\n",
    "#   Degree of freedom 1 --> Fy\n",
    "#   Degree of freedom 2 --> Mz\n",
    "j = 1\n",
    "k = 0\n",
    "for i in range(0,nel):\n",
    "    node1 = inci[i,1] # first node element\n",
    "    node2 = inci[i,2] # second node element\n",
    "    if coord[int(inci[i,1]),1] <= 0.5:\n",
    "        loc = [int(2*node1-2), int(2*node1-1)]\n",
    "        f[loc] = f[loc] - [[q*l/2],[q*l**2/12]]\n",
    "    if coord[int(inci[i,1]),1] <= 0.5:\n",
    "        loc = [int(2*node2-2), int(2*node2-1)]\n",
    "        f[loc] = f[loc] - [[q*l/2],[q*l**2/12]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Global matrix assembly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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,1), 2*nnos-np.size(bc,1)))\n",
    "f_aux  = f[maskv.data]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Displacement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAGmCAYAAACOfxn/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3RVVd7G8Wen90AIoYUOgQRISIigiAVFBwUsqAhWxoK9d1ERHRv2Oo46dgRkFFFQ7GMdhIQSEhKKhF4DJIT05J73D3FeZaJAcpN9y/ezVtZKyL3nPJoN6+Fwzm8bx3EEAAAA4P8F2A4AAAAAeBpKMgAAALAfSjIAAACwH0oyAAAAsB9KMgAAALAfSjIAAACwH2sl2RjzqjFmuzEm103Hm2eMKTbGzNnv1183xhQaY5bs++jvjvMBAADAd9m8kvy6pOFuPN6jks7/g+/d4jhO/30fS9x4TgAAAPggayXZcZxvJe367a8ZY7rvuyKcbYz5zhjT+xCO96WkUnfnBAAAgP/xtHuSX5J0jeM4AyTdLOkFNx33AWNMjjHmSWNMqJuOCQAAAB8VZDvAr4wxUZIGS5ppjPn1l0P3fW+0pPvqedsmx3H+coBD3yFpq6QQ/VLCb/uDYwEAAACSPKgk65er2sWO4/zPg3WO47wv6f2GHNRxnC37Pq0yxrymX65QAwAAAH/IY263cBxnj6RCY8xZkmR+kdbY4xpj2v16PEmnSXLLNA0AAAD4LuM4jp0TGzNN0rGS4iVtkzRJ0leS/i6pnaRgSdMdxzmoWyOMMd9J6i0pStJOSRc7jvOpMeYrSa0lGUlLJF3uOM5e9/7XAAAAwJdYK8kAAACAp/KY2y0AAAAAT2Hlwb2QkBAnNTXVxqnhwVatWqWePXvajgEPw7pAfVgXqA/rAvXJzs4uchyn9aG+z8rtFpGRkU5ZWVmznxeeLTMzU1lZWbZjwMOwLlAf1gXqw7pAfYwx2Y7jZB7q+7jdAgAAANgPJRkAAADYDyUZAAAA2A8lGQAAANgPJRkAAADYDyUZAAAA2A8lGQAAANgPJRkAAADYDyUZAAAA2A8lGQAAANgPJRkAAADYDyUZAAAA2I9bSrIxZrgxZoUxZrUx5nZ3HBMAAACwpdEl2RgTKOl5SSdJSpE0zhiT0tjjAgAAALa440ryQEmrHcdZ4zhOtaTpkk79szc4AUHasKtcO0qrVFFd54YIAAAAwO8tXLurwe8NcsP5O0ja8JuvN0oa9GdvqItK0FFTvv7v19GhQUqICVXb2DB1i49SzzZR6pEQpb4dYhUTFuyGiAAAAPAn367coQteXdDg97ujJJt6fs35nxcZM0HSBEkKDI9WVN4sOYHBcoLCVBMSpQ2h0VoXFqMfI1vLCQrbdxSXAst2KLh4g4KK1ypk588KqCl3Q2R4ovz8fGVmZtqOAQ/DukB9WBeoD+sCv3JMgIoHXa6AgIZXXXeU5I2SOv7m60RJm/d/keM4L0l6SZIiIyOd3I9eqfdgjuNo254qrdhWqqUbipW9ro0Wre+g0spMGSP179hCx/VK0Mi09uoaH+mG+PAUmZmZysrKsh0DHoZ1gfqwLlAf1gV+9eZ/1uqe2Xn6x/kDNPzZhh3DHSV5oaSexpiukjZJGivpnIYezBijtrFhahsbpmOSWkuSXC5HyzaV6OsV2/X1ih164ouVevzzlUpLjNUp/TtoVGo7JcSEueE/BQAAAN6spLxGT3y+UoO7t9KJKW0afJxGl2THcWqNMVdL+lRSoKRXHcfJa+xxfysgwCitYwuldWyh64claXNxhebkbNbsJZt1/5zlevDjfA3v01YXHNFZA7vGyZj67gABAACAr3vqy5XaU1Gju0emNKoTuuNKshzH+VjSx+441sFo3yJcE47urglHd9fq7aWasXCD3s3aqLnLtqhXm2iNP7KLRmd0UGhQYHNFAgAAgGWrt+/VW/9Zp7EDOym5XUyjjuX1O+71SIjWxBEpmn/H8ZpyRqqCAo3ueH+Zjpnyb732QyEj5gAAAPzEA3OXKzw4UDeekNToY3l9Sf5VeEigxhzWUXOuGaK3Lx6kTq0iNPmj5Tpqyld6+ds1qqyhLAMAAPiqf+97du3a43sqPiq00cfzmZL8K2OMhvSM17uXHaEZEw5XcrsYPfBxvo5//BvNWrxRLtf/TKcDAACAF6upc+lvc/PVpVWELhzcxS3H9LmS/FuDurXSWxcP0tRLBqllZLBumLFUo577Xj/+XGQ7GgAAANxk6vx1Wr19ryaOSFFIkHvqrU+X5F8d2SNeH141RE+P7a+Sihqd8/JPunbaYm3fU2k7GgAAABphd1m1nvxilYb0iNew5AS3HdcvSrL0yxi5U/t30Bc3HqPrh/XUvLytOu7xb/TP7wtVW+eyHQ8AAAAN8PSXq1RaWaO7Ria7dQyw35TkX4UFB+r6YUn67PqjNaBzS90/Z7lOe+EHFWzdYzsaAAAADsGqbaV6a/46nTOok3q3bdzIt/35XUn+VZf4SL3+18P0wrkZ2lJcqVHPfq/nv17NVWUAAAAv8be5+YoICdSNJ/Ry+7H9tiRLv0zCOLlfO312w9E6sU9bPfrpCo3++49ata3UdjQAAAD8ia8LtuublTt03fE9FRcZ4vbj+3VJ/lWrqFA9f06Gnj8nQxt3V2jEM9/rtR8K5TiMiwMAAPA0NXUu3T93ubrFR+qCI7o0yTkoyb8xIvWXq8pHJ8Vr8kfLdembWdpVVm07FgAAAH7jrf+s05odZZo4ItltI9/2R0neT3xUqF6+IFOTRqXo25VFOvnp7zR/zU7bsQAAACBpV1m1nvpipY7qGa/jertv5Nv+KMn1MMbor0d21ftXDlZ4SKDOeXm+nv5iFbv1AQAAWPbUFytVVl2nu0emuHXk2/4oyX+ib4dYfXTNEJ3av4Oe/GKlLn0zSyUVNbZjAQAA+KWV20o19af1OndQJyW1iW7Sc1GSDyAqNEhPjEnT5FP66JuVO3Ta8z9oJdMvAAAAmpXjOLp/znJFhgTqhmFJTX4+SvJBMMbowsFd9M6lh6u0slanPf+DPl62xXYsAAAAv/FVwXZ9t6pI1w9LUssmGPm2P0ryIRjYNU5zrhmiXm2jdeXURXr8sxXcpwwAANDEqmtd+tvcfHVvHanzj+jcLOekJB+itrFhmj7hcI3JTNSzX63WtdMXq7KmznYsAAAAn/Xmf9aqsKhMd41MUXBg89TXoGY5i48JDQrUI2ekqlvrKD38SYE2F1fopQsyFR8VajsaAACAT9m5t0pPf7lKxyS11tBeTTfybX9cSW4gY4wuP6a7XjwvQ8u37OGBPgAAgCbwxOcrVV5dp7tHJjfreSnJjTS8bzu9e9kRqqp16Yy//6if2HgEAADALQq27tG0Bet1/uGd1SOhaUe+7Y+S7AapiS30wVVHqk1MmM5/dYE+zdtqOxIAAIBX+3XkW0x4sK4f1rPZz09JdpMOLcI187Ij1Kd9jK54O1vTFqy3HQkAAMBrfb58m35YvVM3DEtSi4imH/m2P0qyG7WMDNHUSwbp6KTWuuP9ZXr2y1VyHEbEAQAAHIqq2jo98HG+eiRE6ZxBnaxkoCS7WURIkF6+IFOj0zvo8c9X6t4P85ilDAAAcAje+HGt1u0s193NOPJtf4yAawLBgQF67Kw0tYoK0cvfFaq0qlZTzkhVkKUfMgAAgLco2lulZ79creN6J+iYpNbWclCSm0hAgNHEESmKCQvW45+vVFWtS0+d3d/a34YAAAC8weOfrVRFTZ0mjmjekW/7oyQ3sWuO76nQ4AA9+HGBqmtdeu6cdIUGBdqOBQAA4HGWb96jGQvXa/zgrureOspqFi5rNoMJR3fX5FP66PPl2zThzWy2sQYAANiP4zia/FGeYsODdd3xzT/ybX+U5GZy4eAuenh0P327aocuen2hyqpqbUcCAADwGPNyt+qnwl266cReio0Ith2Hktycxg7spCfGpGn+mp26+I2FqqjmijIAAEBlzS8j33q3jda4gXZGvu2PktzMTk9P1BNj+uunwl269M0sbr0AAAB+75Xv1mjj7grdMypFgQHGdhxJlGQrTkvvoEfPTNMPPxfpsreyVVVLUQYAAP5pa0mlnv/6Zw3v01aDu8fbjvNflGRLzhyQqAdP76dvVu7QVVMXqbrWZTsSAABAs3tkXoHqHEd3nmx35Nv+KMkWjRvYSfef2kdf5G/XNdMWqaaOogwAAPzHovW7NWvxJl16VFd1ahVhO87vUJItO/+ILrpnZIo+zdum62csUS1FGQAA+AGXy9Hkj5YrITpUVx7bw3ac/8FmIh7goiFdVety6cGPCxQRHKgpZ6bKGM+4aR0AAKApzFq8SUs3FOvxs9IUGep5ldTzEvmpCUd3V1lVnZ7+cpViw4M1cUQyRRkAAPikvVW1emRegdI6ttDp6R1sx6kXJdmDXD+sp0oqavTK94VqGRmiq4Z63j89AAAANNYLX6/W9tIqvXj+AAV4yMi3/VGSPYgxRveMTFFJRY0e/XSFYsODdd7hnW3HAgAAcJv1O8v1yveFGp3eQRmdWtqO84coyR4mIMBoypmp2lNRo7tn5yomPFinpLW3HQsAAMAtHvw4X0EBRrcO7207yp9iuoUHCg4M0PPnZuiwLnG6ccYSfb1iu+1IAAAAjfbj6iLNy9uqq4b2UNvYMNtx/hQl2UOFBQfqlQsz1atttK54O1vZ63bbjgQAANBgtXUu3TdnuRJbhuviIV1txzmgRpVkY8xZxpg8Y4zLGJPprlD4RUxYsN64aKDaxoTpkjcWas2OvbYjAQAANMi0hRtUsLVUE09OVlhwoO04B9TYK8m5kkZL+tYNWVCP+KhQvf7XgQowRuNfW6gdpVW2IwEAABySkvIaPfHZCh3eLU7D+7a1HeegNKokO46T7zjOCneFQf26xEfqn+MP0/bSSl38xkKVV9fajgQAAHDQnvpypUoqanTPyD5esw9Es023MMZMkDRh3+fKzOTujEMVHJ+knLRxSrviGcXkTJdxfGsL6/z8fNYF/gfrAvVhXaA+rAvPVBvZWsWDrlDY5kW64JRJtuMcNOM4zp+/wJgvJNV3XXyi4ziz973m35Judhwn62BOGhkZ6ZSVlR1iVEjS1J/WaeKsXI0b2FEPnt7Pa/42djAyMzOVlXVQSwh+hHWB+rAuUB/WhedxHEcXvLpASzYU6983H6tWUaHNnsEYk+04ziH/7emAV5IdxxnWsEhoCucO6qwtxZV67uvVah8brmuO72k7EgAAQL2+Ktiu71YV6e6RKVYKcmOwmYgXuunEJG0uqdDjn69UuxbhOnNAou1IAAAAv1Nd69Lf5uare+tIXXCE9+0g3NgRcKcbYzZKOkLSXGPMp+6JhT9jjNHDo1M1pEe87ng/R/PX7LQdCQAA4Hfe+HGtCovKdPfIFAUHet/WHI2dbjHLcZxEx3FCHcdp4zjOX9wVDH8uJOiXXfk6t4rU5W9nq7CIe7wBAIBn2FFapWe+XKXjeifo2F4JtuM0iPfVevxXbHiwXr3wMAUYo4teX6ji8mrbkQAAAPT4ZytUUVOniSOSbUdpMEqyl+vUKkL/OH+ANu2u0OVvZ6u61rfGwgEAAO+Su6lEM7I2aPzgLureOsp2nAajJPuAw7rE6ZEz+2n+ml2664NlOtBYPwAAgKbgOI4mf5SnuIgQr5/AxXQLH3F6eqLW7CjTs1+tVvfWUbrsmO62IwEAAD8zd9kWLVy7Ww+N7qfY8GDbcRqFkuxDbhiWpDVFZXp4XoG6xEfqL328Y290AADg/Sqq6/TQxwVKaRejMZkdbcdpNG638CEBAUaPn5Wm1MQWun76EuVuKrEdCQAA+ImXvl2jTcUVmjQqRYEB3r8jMCXZx4QFB+rlCwaoRUSwJryZpaK9VbYjAQAAH7e5uEJ//2a1RvRrp0HdWtmO4xaUZB+UEB2ml87P1M6yal359iImXgAAgCb18CcFchzp9pN6247iNpRkH9UvMVZTzkzVgrW7dO9HebbjAAAAH5W1dpc+XLpZE47upo5xEbbjuA0P7vmwU/t3UP6WUr34zc9KaRej8w73vn3TAQCA56pzOZr0YZ7axYbpimN9a7IWV5J93C1/6aWhvVrr3g/z9NOanbbjAAAAHzJj4Qblbd6jO09OVkSIb117pST7uMAAo6fHpatTqwhdMXWRNu4utx0JAAD4gJLyGj36aYEGdo3TyNR2tuO4HSXZD8SEBevlCzJVU+fShDezVV5dazsSAADwck9+sVIlFTW6d1QfGeP9I9/2R0n2E91bR+mZcenK37pHt8zMYetqAADQYAVb9+it+et07qDOSmkfYztOk6Ak+5GhvRJ02/Demrtsi17498+24wAAAC/kOI4mf7hc0WFBuvGEJNtxmgwl2c9cdnQ3jUprr8c/W6HvVu2wHQcAAHiZT3K36j9rduqmE3upZWSI7ThNhpLsZ4wxeuSMfuqZEK1rpy3Whl08yAcAAA5ORXWdHpibr+R2MTpnYCfbcZoUJdkPRYQE6cXzB6i2ztEVU7NVWVNnOxIAAPACL37zszYVV+jeUSkKDPC9h/V+i5Lsp7rGR+rJs/srd9Me3f1BLg/yAQCAP7VhV7le/OZnjUprr0HdWtmO0+QoyX5sWEobXXtcD83M3qh3Fqy3HQcAAHiwBz/OV4AxuuOk3rajNAtKsp+7bliSjkn6ZUe+Ret3244DAAA80A+ri/RJ7lZdNbS72rcItx2nWVCS/VxggNHTY/urbWyYrnx7kXaUVtmOBAAAPEhNnUuTP8pTx7hwXXJUN9txmg0lGWoREaK/nztAu8urdc20Raqtc9mOBAAAPMTb89dp5ba9umtEisKCA23HaTaUZEiS+naI1UOj+2n+ml16ZF6B7TgAAMAD7NxbpSc+X6mjesbrxJQ2tuM0K0oy/mt0RqIuOKKzXv6uUHNzttiOAwAALHvssxWqqK7TpFEpMsa3R77tj5KM37lrRIrSO7XQbe/laM2OvbbjAAAAS5ZtLNH0hRt04eAu6pEQbTtOs6Mk43dCggL03DkZCgo0unLqIjYaAQDADzmOo3s/ylOryBBdN6yn7ThWUJLxPzq0CNeTZ/dXwdZS3TM713YcAADQzGYv2azsdbt16196KyYs2HYcKyjJqNfQXgm6emgPvZu1UTOzNtiOAwAAmsneqlo9+HG+UhNjdeaARNtxrKEk4w/dcEKSjujWSnfPzlXB1j224wAAgGbw/Nertb20Svee0kcBAf71sN5vUZLxhwIDjJ4e11/RYcG6cuoi7a2qtR0JAAA0ocKiMv3zu0KdkZGojE4tbcexipKMP5UQHaZnxqZrbVGZbn8vR47j2I4EAACayN/mLFdIUIBuG97LdhTrKMk4oCO6t9JNJ/bSnJwtemv+OttxAABAE/i6YLu+LNiua4/voYSYMNtxrKMk46BccUx3De3VWvfPWa6lG4ptxwEAAG5UXevSfXOWq1t8pMYP7mo7jkegJOOgBAQYPTGmvxKiw3Tl1EUqKa+xHQkAALjJP78vVGFRme4elaKQIOqhREnGIWgZGaLnz83Q9tJK3TRziVwu7k8GAMDbbSmp0LNfrdIJKW00tFeC7Tgeg5KMQ9K/YwtNPDlZX+Rv18vfrbEdBwAANNIDc/NV53J0z8gU21E8CiUZh+zCwV10Ut+2evTTFVq0frftOAAAoIF+/LlIc3K26Ipju6tjXITtOB6FkoxDZozRw2ekqm1smK55ZzH3JwMA4IVq6lyaNDtPHePCdfkx3W3H8TiUZDRIbHiwnjsnQ9v2VOo25icDAOB13vhxrVZt36t7RvZRWHCg7Tgeh5KMBuvfsYVuHd5L8/K2Mj8ZAAAvsn1PpZ76YpWG9mqtYck8rFcfSjIa5ZIh3TS0V2v9bU6+8jaX2I4DAAAOwkOfFKi61qVJo/rIGGM7jkdqVEk2xjxqjCkwxuQYY2YZY1q4Kxi8Q0CA0eNj+qtlZLCufmex9lbV2o4EAAD+xILCXZq1eJMmHN1NXeIjbcfxWI29kvy5pL6O46RKWinpjsZHgreJiwzR02PTtW5nme6atYz7kwEA8FC1dS7dMztXHVqE66qhPWzH8WiNKsmO43zmOM6vlw7nS0psfCR4o8O7tdJ1xyfpgyWbNTN7o+04AACgHm/PX6eCraW6a0SywkN4WO/PGHdd9TPGfCRphuM4b//B9ydImrDv8wEZGRluOS88hyOjPRkXqCY2US0WvKSgsh2H9P78/HwlJyc3UTp4K9YF6sO6QH1YF3/OFRKp3Udco6A9mxSz+C35y53I2dnZ2Y7jZB7q+w5Yko0xX0hqW8+3JjqOM3vfayZKypQ02jmI1h0ZGemUlZUdalZ4ge17KnXS098pPipUs68+8pBGymRmZiorK6sJ08EbsS5QH9YF6sO6+HO3zFyqD5Zs0ifXHa0eCVG24zQbY0yDSvIBb7dwHGeY4zh96/n4tSBfKGmkpHMPpiDDtyXEhOmJs/trxbZSTf5oue04AABAUva63ZqZvVEXDenqVwW5MRo73WK4pNskneI4Trl7IsHbHZPUWpcf013TFqzXR0s3244DAIBfq3M5mvRhrtrEhOra43rajuM1Gjvd4jlJ0ZI+N8YsMca86IZM8AE3nZikjE4tdMf7y7RuJ7fWAABgy7QF65W7aY8mjkhRZGiQ7Theo7HTLXo4jtPRcZz++z4ud1cweLfgwAA9My5dAUa6+p3Fqq512Y4EAIDf2V1Wrcc+W6HDu8VpVGo723G8CjvuockktozQo2eladmmEj322QrbcQAA8DtTPl2h0spaTT6lLzvrHSJKMprUX/q01bmDOumlb9fo25WHNhIOAAA0XM7GYk1fuF7jB3dRr7bRtuN4HUoymtxdI1LUMyFKN767VEV7q2zHAQDA57lcju6enadWkaG6fhgP6zUEJRlNLjwkUM+ek649lTW6eeZSuVxMCgQAoCnNzN6gpRuKdefJvRUdFmw7jleiJKNZ9G4bo7tGJOvfK3botR/X2o4DAIDPKimv0SPzVuiwLi11enoH23G8FiUZzeb8wztrWHIbPfJJgXI3ldiOAwCAT3r0swIVl1fr3lP68LBeI1CS0WyMMZpyZqpaRgbr2umLVV5dazsSAAA+ZemGYk39ab0uHNxFfdrH2o7j1SjJaFZxkSF6ckx/FRaVafKHbFsNAIC71Lkc3T07V/FRobrhhCTbcbweJRnNbnCPeF1xTHfNyNqguTlbbMcBAMAnvLNgvXI2luiuEcmK4WG9RqMkw4obTkhS/44tdPv7Odq4u9x2HAAAvFrR3io9Oq9AR3RrpVPS2tuO4xMoybAiODBAz4xNl+NI109foto6tq0GAKChHvq4QBU1dbr/NB7WcxdKMqzp1CpCfzutr7LW7dYzX622HQcAAK+0oHCX3lu0UZce1U09EthZz12CbAeAfzstvYO+XbVDz321StEtOtmOAwCAV6mpc+nuD3LVoUW4rj6uh+04PoUrybDuvlP7qlNchEr7nqGS8hrbcQAA8Bpv/LhWK7aVatKoFEWEcO3TnSjJsC4qNEhPj02XKyRat7+fI8dh22oAAA5ka0mlnvx8pY7rnaATUtrYjuNzKMnwCGkdWyji5y/1Se5WTVuwwXYcAAA83v1zl6vW5ejeUTys1xQoyfAY4et+1JAe8bpvTp5WbSu1HQcAAI/17codmpuzRVcN7aFOrSJsx/FJlGR4DCNHT4xJU0RIkK6ZtliVNXW2IwEA4HGqaus06cM8dY2P1ISju9mO47MoyfAoCTFheuysVBVsLdUj8wpsxwEAwOO89M0aFRaVafIpfRQWHGg7js+iJMPjHNe7jcYP7qLXflirrwu2244DAIDH2LCrXM99vVoj+rXT0UmtbcfxaZRkeKTbT+qt3m2jdcu/lmpHaZXtOAAAWOc4jiZ9mKegAKO7R6bYjuPzKMnwSGHBgXpmXLpKK2t108ylcrkYCwcA8G+fL9+mrwq26/phSWobG2Y7js+jJMNjJbWJ1l0jU/Ttyh169YdC23EAALCmvLpWkz9arl5tojX+yC624/gFSjI82nmDOmlYchs9Mq9AuZtKbMcBAMCK575arU3FFbr/tL4KDqS+NQf+L8OjGWM05cxUtYwI0XXTF6u8utZ2JAAAmtXq7aV6+bs1Gp3RQQO7xtmO4zcoyfB4cZEhevLs/lpTVKb75+TbjgMAQLNxHEd3zspVREiQ7jw52XYcv0JJhlc4ske8JhzdTdMWrNe83C224wAA0CxmZm/UgsJduuOk3oqPCrUdx69QkuE1bjqhl1ITY3Xbe8u0paTCdhwAAJrUrrJqPfRxvjI7t9SYzI624/gdSjK8RkhQgJ4em66aOpeun75EdYyFAwD4sAc/zldpZa0eHN1PAQHGdhy/Q0mGV+kaH6nJp/TRT4W79OI3P9uOAwBAk5i/Zqf+lb1Rlx7dTUltom3H8UuUZHidMwckamRqOz3x+UotXr/bdhwAANyqqrZOE2ctU8e4cF17XE/bcfwWJRlexxijB07vp7YxYbpu+hKVVtbYjgQAgNu89M0a/byjTPed2lfhIYG24/gtSjK8Umx4sJ4a218bd5dr0uw823EAAHCLtUVlevbr1RrRr52G9kqwHcevUZLhtQ7rEqdrjuup9xdv0geLN9mOAwBAoziOo7tn5yo0MED3jEqxHcfvUZLh1a45rocyO7fUXR/kav3OcttxAABosA+XbtZ3q4p08196qU1MmO04fo+SDK8WFBigp8b2lzHSdTMWq6bOZTsSAACHrKSiRvfPyVdqYqzOO7yz7TgQJRk+ILFlhB44vZ8Wry/WM1+ush0HAIBDNmVegXaVVenB0/spkJnIHoGSDJ9wSlp7nTkgUc99vVrz1+y0HQcAgIO2aP1uvbNgvcYP7qq+HWJtx8E+lGT4jHtP6aPOcRG6YcYSlZQzFg4A4Plq6ly68/1lahsTphtPTLIdB79BSYbPiAoN0tNj07WjtEp3zMqR47BtNQDAs732Q6EKtpZq0qg+igoNsh0Hv0FJhk9J69hCN53YSx8v26p3szbYjgMAwB/auLtcT36+SsOSE/SXPm1sx8F+KMnwOZcd3U2Du7fSvR8u1ygaYJ4AACAASURBVM879tqOAwDA/3AcR/d++MtmWPee0kfG8LCep6Ekw+cEBBg9Maa/woIDdO20xaqqrbMdCQCA3/k0b6u+yN+uG07oqcSWEbbjoB6NKsnGmPuNMTnGmCXGmM+MMe3dFQxojLaxYXrkjFTlbd6jxz5dYTsOAAD/VVJRo3tm56lP+xhddGRX23HwBxp7JflRx3FSHcfpL2mOpHvckAlwixP7tNV5h3fSy98V6tuVO2zHAQBAkvTIvAIV7a3Sw6NTFRTIP+p7qkb9ZBzH2fObLyMlMU4AHmXiySnqmRClG99dqqK9VbbjAAD83ILCXXrnp/W66Miu6pfITGRPZho7JssY84CkCySVSBrqOE69l+yMMRMkTdj3+YCMjIxGnRe+Jz8/X8nJyW4/bm1UGxUfdqmCd61RzNJ3xKMR3qWp1gW8G+sC9fH0deGYQBUffoWcgCC1/M/zMi5m+jeH7OzsbMdxMg/1fQcsycaYLyS1redbEx3Hmf2b190hKcxxnEkHOmlkZKRTVlZ2qFnh4zIzM5WVldUkx37th0JN/mi5Jp/SRxcO7tIk50DTaMp1Ae/FukB9PH1dPPn5Sj395Sq9cdFAHZPU2nYcv2GMaVBJPuDUasdxhh3ksd6RNFfSAUsy0NzGD+6ib1fu0AMf52tQtzj1bhtjOxIAwI+s2laqF/69Wqf1b09B9hKNnW7R8zdfniKpoHFxgKZhjNGjZ6UpJixY105brMoaxsIBAJqHy+Xo9veXKTI0SHePTLEdBwepsY9UPmyMyTXG5Eg6UdJ1bsgENIn4qFA9PiZNK7ft1QNz823HAQD4iakL1it73W7dNSJFraJCbcfBQWrUJuGO45zhriBAczgmqbUuGdJVr3xfqKOTWuuEFLYBBQA0na0llXrkkwIN6RGvMzI62I6DQ8BwPvidW4b3Ukq7GN36r6XatqfSdhwAgA+b9GGual0uPXB6X7ae9jKUZPid0KBAPTMuXZU1Lt347hK5XIz3BgC437zcrfo0b5uuH5akzq0ibcfBIaIkwy/1SIjSpFEp+mH1Tr383RrbcQAAPmZPZY3umZ2rlHYxumQIW097I0oy/NbZh3XUSX3b6tFPVyhnY7HtOAAAH/LIJ/u2nj6jH1tPeyl+avBbxhg9NLqfWkeH6rrpS1RWVWs7EgDAByxcu0tTf1qvvx7ZVamJLWzHQQNRkuHXWkSE6Mmz+2vtzjLd+2Ge7TgAAC9XVVunO95fpg4twnXjCUm246ARKMnwe4d3a6Wrju2hmdkb9dHSzbbjAAC82N///bNWb9+rB07vq8jQRk3ahWWUZEDSdcN6qn/HFrpz1jJt3F1uOw4AwAut2Fqq57/+ZevpY3sl2I6DRqIkA5KCAwP0zNh0OY50/fQlqq1z2Y4EAPAitXUu3fqvpYoJC9Y9o/rYjgM3oCQD+3RqFaG/ndZXWet26/mvf7YdBwDgRV79oVBLN5Zo8ql9FBcZYjsO3ICSDPzGaekddHp6Bz395Uplrd1lOw4AwAsUFpXp8c9W6sSUNhrRr53tOHATSjKwn/tO7aMOLcN13fQlKqmosR0HAODBXC5Ht/0rR6FBAfrbaWw97UsoycB+osOC9fTYdG3dU6mJs5bJcdi2GgBQv6k/rdOCtbt018gUJcSE2Y4DN6IkA/XI6NRSNwzrqTk5W/Teok224wAAPNDG3eV6+JMCHdUzXmcNSLQdB25GSQb+wBXH9tCgrnG6Z3au1haV2Y4DAPAgjuPozlm5ciQ9eHo/brPwQZRk4A8EBhg9eXZ/BQcG6Nrpi1Vdy1g4AMAv3lu0Sd+u3KHbhvdWx7gI23HQBCjJwJ9o3yJcD4/up5yNJXryi5W24wAAPMD20krdP2e5DuvSUucf3tl2HDQRSjJwACf1a6dxAzvqxW9+1o+ri2zHAQBYds8HeaqoqdPDZ6QqIIDbLHwVJRk4CHePTFHX+Ejd8O4S7Sqrth0HAGDJx8u2aF7eVt0wLEndW0fZjoMmREkGDkJESJCeGZuuXWXVuu29HMbCAYAf2l1WrXtm56pfh1hdelRX23HQxCjJwEHq2yFWtw3vrc+Xb9PUn9bbjgMAaGb3z1mu4vIaTTkzVUGBVChfx08YOAQXHdlVR/WM1/1zlmvVtlLbcQAAzeTrgu16f/EmXTm0h5LbxdiOg2ZASQYOQUCA0eNj0hQVGqRrpi1WZU2d7UgAgCZWUl6j29/PUa820bpqaHfbcdBMKMnAIUqIDtNjZ6WpYGupHplXYDsOAKCJ3TdnuYr2Vuuxs9IUGhRoOw6aCSUZaIChvRM0fnAXvfbDWn1dsN12HABAE/li+Ta9t2ijrjq2u/olxtqOg2ZESQYa6PaTeqt322jdPHOptpdW2o4DAHCz4vJq3TFrmXq3jdbVx/W0HQfNjJIMNFBYcKCeGZeuvVW1uundpXK5GAsHAL5k8kfLtbusWo+PSVNIEJXJ3/ATBxohqU207hqZou9WFenVHwptxwEAuMmneVs1a/EmXX1cD/Vpz20W/oiSDDTSeYM66YSUNnpkXoFyN5XYjgMAaKTdZdWaOCtXKe1idNXQHrbjwBJKMtBIxhg9ckaq4iJDdN30xSqvrrUdCQDQCJM+zFNJxS+3WQSzaYjf4icPuEFcZIieGNNfa4rKdP+c5bbjAAAaaF7uFn24dLOuPa4nm4b4OUoy4CZH9ojXZUd317QFGzQnZ7PtOACAQ7Rzb5UmzspV3w4xuvxYNg3xd5RkwI1uOjFJ6Z1a6Pb3lmltUZntOACAQ3DPh3naU1mjx8/qz20WoCQD7hQcGKBnx6UrMMDo6mmLVFXLttUA4A3m5mzR3Jwtun5Yknq1jbYdBx6Akgy4WWLLCD12VppyN+3Rg3PzbccBABxA0d4q3T07V2mJsbrs6G6248BDUJKBJnBCShtdPKSr3vjPOn2ybIvtOACAP+A4jibOWqa9lbV67Kw0BXGbBfZhJQBN5LbhvZXWsYVufS9H63eW244DAKjHe4s26dO8bbr5L0nq2YbbLPD/KMlAEwkJCtBz49IlifuTAcADbdxdrns/zNPArnG6eAi3WeD3KMlAE+oYF6FHz0xTzsYSPfxJge04AIB9XC5HN89cKsdx9PhZaQoMMLYjwcNQkoEmNrxvW40f3EWv/bBWn+ZttR0HACDp1R8KNX/NLk0a1Ucd4yJsx4EHoiQDzeCOk3urX4dY3TJzqTbs4v5kALBp5bZSTfl0hYYlt9FZmYm248BDUZKBZhAaFKjnz8mQ40hXT1us6lqX7UgA4Jeqa126YcYSRYcG6eEz+skYbrNA/dxSko0xNxtjHGNMvDuOB/iiTq0iNOXMVC3dUKwp87g/GQBseParVcrbvEcPje6n+KhQ23HgwRpdko0xHSWdIGl94+MAvu2kfu10wRGd9cr3hfpi+TbbcQDAryxav1vPf71aZw5I1Il92tqOAw/njivJT0q6VZLjhmMBPu/Ok5PVp32Mbpq5VJuKK2zHAQC/UF5dqxtnLFG72HBNGpViOw68gHGchndbY8wpko53HOc6Y8xaSZmO4xT9wWsnSJqw7/MBGRkZDT4vfFN+fr6Sk5Ntx2gWdeFxKh50mQL37lBs9qsyDvco/xF/Whc4eKwL1OfP1sXeXiNUmZip2OzXFVy8rpmTwabs7Oxsx3EyD/V9ByzJxpgvJNX3bxITJd0p6UTHcUoOVJJ/KzIy0ikrKzvUrPBxmZmZysrKsh2j2czJ2ayr31msS4Z01V0juarxR/xtXeDgsC5Qnz9aF9+s3KELX12gS4/qqokj+PPW3xhjGlSSgw70Asdxhv3BCftJ6ipp6b4nQxMlLTLGDHQch2GwwAGMTG2vhYW79Mr3hRrQuaVO6tfOdiQA8DnF5dW6ZeZSJbWJ0k0n9rIdB16kwfckO46zzHGcBMdxujiO00XSRkkZFGTg4E0ckaL+HVvoln/lqLCIf10BAHdyHEcTP8jV7vJqPTGmv8KCA21HghdhTjJgUUhQgJ4/N0PBgUZXvJ2tiuo625EAwGe8t2iT5uZs0fXDktS3Q6ztOPAybivJ+64oH/B+ZAC/16FFuJ48u79WbCvVXR/kqjEP0wIAfrFuZ5kmzc7VoK5xuvyY7rbjwAtxJRnwAMf2StA1x/XUe4s2asbCDbbjAIBXq6lz6brpSxQYYPTk2f0VGMCuejh0lGTAQ1x3fE8d1TNe93yYp9xNJbbjAIDXevbLVVqyoVgPjU5V+xbhtuPAS1GSAQ8RGGD01Nn91SoyRFdOXaSSihrbkQDA6ywo3KXn9u2qNyKVqUFoOEoy4EFaRYXquXMytLm4Qje9u5T7kwHgELiCwnTDjCXqGBehe0/pYzsOvBwlGfAwAzq31J0nJ+uL/G36x7drbMcBAK9R1muEtu6p1FNn91dU6AG3ggD+FCUZ8EB/PbKLRvRrpynzCjR/zU7bcQDA481avFFV7VJ1/fE9ld6ppe048AGUZMADGWP08Bn91KVVpK5+Z7G276m0HQkAPNaGXeW6+4M8Be1epyuH9rAdBz6Ckgx4qOiwYP39vAHaW1Wjq6ctVm2dy3YkAPA4tXUuXTd9sYyRovPeY9wb3IaSDHiwXm2j9dDoflpQuEsPf1JgOw4AeJznvl6tReuL9cDp/RRYyfhMuA8lGfBwp6cnavzgLnrl+0LNXrLJdhwA8BhZa3fpmS9XaXR6B52S1t52HPgYSjLgBSaOSNZhXVrqtvdylL9lj+04AGBdcXm1rp22WIktIzT5VMa9wf0oyYAXCA4M0PPnZig2PFiXvZWtknI2GgHgvxzH0a3/ytGOvVV67px0RYcF244EH0RJBrxEQnSYXjh3gLaUVOi6GYtV52KjEQD+6a356/TZ8m26bXhvpSa2sB0HPoqSDHiRAZ1batKoPvr3ih16+ouVtuMAQLPL21yiv83J13G9E3TxkK6248CHUZIBL3PuoE4ak5moZ75arc/yttqOAwDNpqyqVte8s1gtI4P16JmpMoZxb2g6lGTAyxhjdN+pfZWaGKsb312qn3fstR0JAJrFPbPzVLizTE+dna5WUaG248DHUZIBLxQWHKi/nzdAIUEBuuytbO2tqrUdCQCa1PuLNuq9RRt1zXE9dUT3VrbjwA9QkgEv1aFFuJ47J11rduzVLTOXynF4kA+Ab1qzY6/u+iBXA7vG6drj2HYazYOSDHixwd3jdcdJyfokd6te/GaN7TgA4HZVtXW6+p3FCgkK0NNj+ysokOqC5sFKA7zcJUd11ai09nr00wJ9s3KH7TgA4FYPfVyg5Vv26LEz09QuNtx2HPgRSjLg5YwxeuSMfurVNkZXv7NIa3iQD4CPmJe7Ra//uFZ/PbKLhqW0sR0HfoaSDPiAiJAgvXzBAAUHBuiSN7O0p5Id+QB4t7VFZbplZo7SEmN1+0m9bceBH6IkAz4isWWE/n5uhtbvLNe109iRD4D3qqyp05VTFykgwOj5czMUGhRoOxL8ECUZ8CGDurXS5FN/2ZFvyqcFtuMAQINM/mi5lm/ZoyfGpCmxZYTtOPBTQbYDAHCvcwd1Vv6WPfrHN2uU3DZGp6V3sB0JAA7arMUbNW3Bel1+THcdn8x9yLCHK8mAD5o0qo8GdY3Tbe/laOmGYttxAOCgrNxWqjvf/2Ue8s0nJtmOAz9HSQZ8UHBggF44N0Oto0M14a0sbd9TaTsSAPypsqpaXTl1kSJDA/XcuHTmIcM6ViDgo1pFherlCzJVWlmry97OVmVNne1IAFAvx3E0cdYy/bxjr54Zm66EmDDbkQBKMuDLktvF6IkxaVq8vlh3fZDL1tUAPNI7C9brgyWbdcOwJA3uEW87DiCJkgz4vOF92+n6YT31r+yN+uf3hbbjAMDv5G4q0eQPl+vopNa6emgP23GA/2K6BeAHrj2up1ZsLdUDH+era3wkT4wD8Agl5TW6cuoixUWG6Kmz+ysgwNiOBPwXV5IBPxAQYPTEmP7q2z5W10xbrOWb99iOBMDPuVyOrpuxWFtKKvT8uRmKiwyxHQn4HUoy4CfCQwL1yoWZig0P1sVvLGTiBQCrnvpylf69YocmjeqjAZ1b2o4D/A9KMuBH2sSE6ZULM1VSUaNL3sxSRTUTLwA0vy+Wb9MzX67SWQMSde6gTrbjAPWiJAN+pk/7WD0zNl3LNpXoxneXyOVi4gWA5lNYVKYbZixRvw6xuv+0vjKG+5DhmSjJgB8altJGE09O1ie5W/XYZytsxwHgJ8qqanXZW1kKCjT6+3kZCgsOtB0J+ENMtwD81MVDuurnHWV64d8/q2t8pM7K7Gg7EgAf5jiObn0vR6u379WbFw1SYssI25GAP8WVZMBPGWN036l9NKRHvO6ctUzz1+y0HQmAD3vlu0LNzdmiW4f31pCebBgCz0dJBvxYcGCAnj83Q53iInT529kqLCqzHQmAD/pxdZEe+iRfJ/Vtq8uO7mY7DnBQKMmAn4sND9ar4w+TkXTx6wtVXF5tOxIAH7K5uEJXT1usbq2j9OhZaTyoB69BSQagzq0i9dIFmdq4u0KXvpmlyhpGwwFovMqaOl3+draqa136x/kDFBXKo1DwHpRkAJKkw7rE6fExaVq4drduencpo+EANIrjOLrtvRwt21SiJ8/ur+6to2xHAg4Jf6UD8F+j0tpra0mlHvg4X+1iw3TXyBTbkQB4qRe/WaPZSzbrlr/00gkpbWzHAQ5Zo64kG2PuNcZsMsYs2fdxsruCAbDjkqO6avzgLnrl+0K9+n2h7TgAvNCX+ds05dMCjUxtpyuP7W47DtAg7riS/KTjOI+54TgAPIAxRnePTNGWkgrdP3e52sWG6aR+7WzHAuAlVm8v1XXTl6hP+xg9eiYP6sF7cU8ygP8RGGD09Nh0pXdsoetnLFH2ul22IwHwAsXl1brkjSyFBQfqpfMzFR7CjnrwXsZxGv5wjjHmXknjJe2RlCXpJsdxdv/BaydImrDv8wEZGRkNPi98U35+vpKTk23HwG+4giNUnHmxnOAIxWa9oqDy5t9whHWB+rAuPI9jArSn/3mqadlZsdmvKbhkY7NnYF2gPtnZ2dmO42Qe6vsOWJKNMV9IalvPtyZKmi+pSJIj6X5J7RzHuehAJ42MjHTKyti0AL+XmZmprKws2zGwn3U7yzT6hR8VERqo9684Uq2jQ5v1/KwL1Id14Xkmf5Sn135YqylnpmqMpW3uWReojzGmQSX5gLdbOI4zzHGcvvV8zHYcZ5vjOHWO47gkvSxpYEPCA/BcnVtF6tXxh6motFoXv7FQZVW1tiMB8DDvLtyg135Yq4uHdLVWkAF3a+x0i98+zXO6pNzGxQHgidI6ttBz56Qrb/MeXf52tqpq2WwEwC/mr9mpiR8s01E943XHSb1txwHcprEP7k0xxiwzxuRIGirpBjdkAuCBjk9uo0fOSNV3q4p047tLVcdmI4DfW7Njry57K1ud4iL03DkZCgpkHgB8R6NGwDmOc767ggDwfGcOSNTusmo98HG+4iJCdN+pfRjvBPip3WXVuuj1hQoKMHpt/EDFhgfbjgS4FTvuATgklx7dTUVlVfrHN2sUFxmiG05Ish0JQDOrqq3TZW9la3NJpaZdOkidWkXYjgS4HSUZwCG7fXhv7S6r1tNfrlJcZIguHNzFdiQAzcRxHN3x3jItWLtLz4xL14DOcbYjAU2CkgzgkBlj9ODp/bS7vEb3fpSnFhHBOrV/B9uxADSD575arfcXb9JNJyTplLT2tuMATYY77AE0SFBggJ4dl67DusTppneX6puVO2xHAtDEZi/ZpMc/X6nR6R109XE9bMcBmhQlGUCDhQUH6pULM9WzTbQufytbi9bXu+EmAB+QvW6XbvlXjgZ2jdNDZ/TjoV34PEoygEaJCQvWGxcdpoSYUP31tYXK37LHdiQAbrZuZ5kmvJmt9rFh+sd5AxQaFGg7EtDkKMkAGi0hOkxvXzxI4cGBOv+fP2n19r22IwFwk6K9Vbrg1QVyOY5eHX+YWkaG2I4ENAtKMgC36BgXoamXDpIknffKT9qwq9xyIgCNVV5dq4tfX6itJZV65cLD1K11lO1IQLOhJANwm+6to/TWxYNUUVOncS/P15aSCtuRADRQbZ1LV7+zWMs2lejZceka0Lml7UhAs6IkA3Cr5HYxevOigSour9G5L/+kHaVVtiMBOESO4+iuD3L1VcF23XdqX53Yp63tSECzoyQDcLu0ji302l8P05aSSp3/z5+0u6zadiQAh+CZL1dr+sINunpoD513eGfbcQArKMkAmsRhXeL08gWZWlNUpgtfW6A9lTW2IwE4CDMWrteTX6zUGRmJuulEtp2H/6IkA2gyQ3rG6+/nZmj55j266LWFKquqtR0JwJ/4qmCb7pyVq6OTWuthZiHDz1GSATSp45Pb6Omx6Vq8oVh/fZ2iDHiqhWt36cqpi5TcLlovnJuh4EAqAvwbvwMANLkRqe301Nn9lbV2F0UZ8EDLN+/RRa8vVPvYcL3+14GKCg2yHQmwjpIMoFmMSmuvp8emU5QBD7O2qEwXvLpAUaFBevPigYqPCrUdCfAIlGQAzYaiDHiWrSWVOu+fP6nO5dJbFw9UYssI25EAj0FJBtCsKMqAZ9hdVv3fEY1vXDRQPRKibUcCPAolGUCzoygDdpVV1Wr86wu1ble5Xr4wU6mJLWxHAjwOJRmAFRRlwI6q2jpNeCtLuZtK9Ny4dA3uHm87EuCRKMkArPm1KGev263z//mTSirYcARoStW1Ll359iL9sHqnHh7dj+2mgT9BSQZg1ai09nr+nHQt21Sic16er517q2xHAnxSbZ1L105brC8Ltuv+U/vorMyOtiMBHo2SDMC64X3b6eULMrV6+16NfWm+tu2ptB0J8Cl1Lkc3vLtU8/K26u6RKTr/iC62IwEej5IMwCMc2ytBb1w0UJuLKzTmH//Rxt3ltiMBPsHlcnTrv3L00dLNum14b108pKvtSIBXoCQD8BiHd2ulty8ZpN1l1Rrz4n+0Zsde25EAr+Y4jiZ+sEzvLdqoG4Yl6Ypju9uOBHgNSjIAj5LeqaWmTzhCVbUujfnHfNVGJtiOBHglx3F074d5mrZgg64a2l3XHt/DdiTAq1CSAXiclPYxmnHZEQoMkEoG/FVLNxTbjgR4Fcdx9MDcfL3xn3W69KiuuvnEXjLG2I4FeBVKMgCP1CMhSjMvGyxTW6lxL8/Xd6t22I4EeAXHcTT5o+V65ftCjR/cRXeenExBBhqAkgzAY3VqFaEWWf9U51aRuuj1hZq9ZJPtSIBHc7kc3T07V6//uFYXD+mqSaNSKMhAA1GSAXi0gOq9mnHZ4cro1FLXTV+iV78vtB0J8Egul6M7Zy3T2/PX67JjuumuEVxBBhqDkgzA48WEBeuNiwZqeJ+2um/Ock2ZVyDHcWzHAjxGncvRre/laPrCDbp6aA/dPrw3BRloJEoyAK8QFhyo58/N0LiBnfTCv3/Wbe/lqLbOZTsWYF1tnUs3vbtE/8r+ZczbzX/hIT3AHYJsBwCAgxUYYPTg6X3VOjpUz3y5SrvKqvXMuHRFhPBHGfxTTZ1LN8xYojk5W3TLX3rpqqGMeQPchSvJALyKMUY3npCk+0/rq68KtmvsS/O1vZRtrOF/KmvqdMXb2ZqTs0V3nNSbggy4GSUZgFc6//DOeun8TK3atlenP/+jVm4rtR0JaDallTW68NUF+rJgu+4/tY8uO4ad9AB3oyQD8FrDUtpo5uVHqKbOpTNe+FHfryqyHQlocjv3Vumcl39S9rrdeurs/jr/iC62IwE+iZIMwKv17RCrWVcdqfYtwjX+tQV6N2uD7UhAk9lcXKEx//iPVm4r1UsXDNCp/TvYjgT4LEoyAK/XoUW4Zl5xhI7o3kq3/itHj326ghFx8DlrduzVmX//Udv3VOnNiwbquN5tbEcCfBolGYBPiAkL1qvjD9PYwzrqua9X6+ppi1VRXWc7FuAWORuLddaL/1FVrUvTJhyuQd1a2Y4E+DzmJgHwGcGBAXpodD91iY/UI/MKtG5nmV46P1PtW4TbjgY02NcF23Xl1EWKiwzRmxcPVPfWUbYjAX6BK8kAfIoxRpcf013/vDBT64rKdcpz3yt73S7bsYAGmbZgvS55M0vdEyI166rBFGSgGVGSAfik43q30ayrBisqNEhjX5qvdxfyQB+8h+M4evyzFbrj/WUa0iNeMyYcoYToMNuxAL/S6JJsjLnGGLPCGJNnjJnijlAA4A49EqI1+6ohOrxbK936Xo4mf5THVtbweNW1Lt00c6me/Wq1zs7sqFcuzFRkKHdHAs2tUb/rjDFDJZ0qKdVxnCpjTIJ7YgGAe8RGBOu18YfpwY8L9OoPhVq1ba+eGZeuuMgQ29GA/1FSUaOrpi7S96uLdMOwJF17fA8ZY2zHAvxSY68kXyHpYcdxqiTJcZztjY8EAO4VFBige0alaMqZqVqwdpdGPfu9lm4oth0L+J3CojKd/sIPmr9mp6acmarrhvWkIAMWmcbMEjXGLJE0W9JwSZWSbnYcZ+EfvHaCpAn7Ph+QkZHR4PPCN+Xn5ys5Odl2DHgYd6+Lmuj2Kk0dI1dotCJXfKKwTVmihngfX/vzorplV5WmjpEcRzE5M/6vvTuPrqo89zj+fTKTAIEMDCEBEgjzTCAqAjK0goi0qLRarmIdar2ita2tt9SrHbztbcvVCtZaKkqdikoRUFFUoAiCjMqcGJnCGDCEkMTM7/0jaS83jeTE5AxJfp+1WOucnH32flg87PPL3u95X0LzDvu7pCapufWFNI5t27Ztc86l1fd9dYZkM3sX6FTLS3OAR4DVwL3ACGAxkOLq2GlUVJQrLCysb63SzKWlpbF161Z/lyEBxht9kVdUyvcWf8TajNNMH9aFR742kFZhwY16DPGu5nS+eH7TYR5evofkuCievnkEXWMj/V1Sk9Wc+kIaj5l9qZBc55hk59zE8z5RkgAAECdJREFUixz0u8DfqkPxZjOrBOKA0/UtRETEV9pFhrHw5hHMW53FY+9lsvd4Pk/OHE5yXJS/S5MWpLyikl+8vpdFGw8zrnc8j98wlDYRof4uS0SqNXRM8mvAeAAz6wWEAWcaWpSIiLcFBRn3Tkzl2VtGcjK/mGvmreet3Sf8XZa0EGcLS5n1zBYWbTzM7aOT+fPNIxSQRQJMQ0PyQiDFzHYDfwVurmuohYhIIBnbK57XZ19OSnwUdz6/nTlLd1FcpuWsxXt2HT3H1fPWs/lgLr+5dhBzpvQjOEgj40UCTYOmgHPOlQIzG6kWERG/SGwfySt3XsbcVRk8te4AWw7lMv/GYfTq2MbfpUkzs3jLER5ctof41uG8cuelDE5q5++SROQLaMU9EREgLCSI/7iqL4u+PZLcwlKmzlvPCx8eRjfHpDEUl1XwwJKd/HjJLtKTY1gx+3IFZJEAp5AsInKBsb3iWXnvGEYmxzBn6W7uemE754rK/F2WNGFHzxZx/R838tct2dw9rifP3jJSi9mINAEKySIiNcS3CWfRLSP5yVV9eGfvKSb/fh0ffKrvJEv9rdpzkimPr+fQmUIW3JTGD6/srfHHIk2EQrKISC2Cgow7xvRgyXcvIzw0mBsXfMjDy/fweam+1Cd1Ky6r4KFlu7njuW10jYlkxezL+Uq/jv4uS0TqQSFZROQiBie14817RjPrsu48+8Ehpjz+PjuOnPV3WRLAsnIK+PofPmDRxsPcdnkyS757Gd01B7dIk6OQLCJSh1ZhwTx8TX9euC2d4rIKrn3yA373dgal5ZX+Lk0CiHOOl7dmM3Xeek7lF/PMrBH89Op+hIXoo1akKdL/XBERD43qGcdb941h+rBE5q/JYtoTG9h3It/fZUkAOFNQwnee28aPXt3JkKR2rLx3NOP6dPB3WSLSAArJIiL10DYilN9dP5gFN6Vx+nwxU+et53dvZ2gBkhbsrd0nufLRdazNOM1PrurD87el07FthL/LEpEGatBiIiIiLdVX+nUkrdtYfvnGPuavyeLNXSf41fSBpKfE+rs08ZFzn5fxsxV7+Nv2Ywzo0paXZgzRAjQizYiuJIuIfEnto8KYO2Mwz906ktKKSr7xp038ZOku8os1r3JztyYjh8mPrWPZR8e5Z0IqS+8apYAs0szoSrKISAONTo1n1X1jePSdTJ5ef5D39p3ioan9mTygE2aaE7c5OX2+hJ+/vpcVHx8ntUNr/jZzuFbOE2mmFJJFRBpBZFgIc6b0Y+rgBH68ZBd3vbCd0alxPHxNf3rEt/Z3edJAzjle2XaUR97Yx+elFdw3sRd3XpFCeEiwv0sTES9RSBYRaUSDEtux4u5RPL/pMHNXZTLpsXXcNjqF2eN7EhmmU25TdOB0AXOW7mbjgc8Y2T2G/5o+kJ4d9IuPSHOnM7aISCMLCQ5i1qhkpgxK4Ncr9/Pk2k95bccxHry6n4ZgNCHni8uYvzqLhRsOEhESzCNfH8ANI7oSpGWlRVoEhWQRES+JbxPO3BmDuTE9iQdf28NdL2znkpQYfjqlHwO6RPu7PPkClZWOpTuO8eu39nP6fAnXD0/k/km96dBG07qJtCQKySIiXja8WwwrZl/Oi5uP8Og7mUydv56vD+3C/Vf2pnN0K3+XJxf4ODuPn63Yw/YjeQxOaseCm9IYoi/mibRICskiIj4QHGT82yXdmDYkgT+s+ZSFGw7y5q4T3D46he+M7UHrcJ2O/enA6QLmrsrkjV0niGsdzm+vG8S1wxI1tEKkBdNZWUTEh9pGhPLA5D58K70rv3k7g3mrs3hpczb3TOjJN0d0JSxE09f7Uk5+MY+99wmLt2QTHhLEPRNSuX10Mm0iQv1dmoj4mUKyiIgfJMVEMu+GoXx7VHd+9eZ+/nPZHp76+wHunZjK9KFdCAlWWPam3MJS/vz+ARZuOEhFpWNmelfuHp9KfJtwf5cmIgFCIVlExI+Gdm3P4u9cwrpPzjB3VQY/enUnf/z7p9w3sRdTBnbW7f5GlpNfzIL3D/D8piMUl1cwdVACP/hqL7rFRvm7NBEJMArJIiJ+ZmaM7RXPmNQ4Vu09xf+symT2Szt4Yk0Ws8enMmlAJ4IVlhvk6Nkinvr7ARZvzaai0jFtcAJ3jetBzw5aSlpEaqeQLCISIMyMK/t3YmLfjry+8zi/f/cT/v3F7aTER3HXFT2ZNiSBUA3DqJftR87yzIZDrNx1AjO4bngid47toSvHIlInhWQRkQATHGRMG9KFqwcl8Nbuk8xfk8UPX/mYR9/J5M6xKVyflkREqJZD/iKl5ZWs3H2ChRsO8XF2Hm3CQ7j5su7cenkyCe005Z6IeEYhWUQkQAUHGVMGdeaqgZ1Yk5HD/NVZPLhsD3PfyWRGWhIz07vRNTbS32UGjE9PF/DylmyWbD/GmYISUuKi+Pm0/lw7LJEoTbEnIvWks4aISIAzM8b36ci43h3YfDCXv2w8zNPrD7Lg/QOM792BmZd2Y3TPuBY5I8b54jJW7j7Jy1uy2Xr4LMFBxvg+HbgxvStjU+P1xUcR+dIUkkVEmggzIz0llvSUWE6eK+bFDw/z4uZs3ntmC/Ftwpk2OIGvDe1C/4S2mDXfcFhYUk5JxwHc8ZetrM08TWl5JSlxUTwwuQ/Th3XR8tEi0igUkkVEmqBO0RF8/6u9uXt8Kqv3n2LpjmMs2niIP68/SK+OrZkyMIEJfTs0m8Cck1/M2szTrN6Xw9rMHIoHXs9H2XncOLIrUwd3ZljX9s3i7ykigUMhWUSkCQsLCWLSgM5MGtCZvKJS3th1gtd2HOOx9zJ59N1MEqIjGN+3AxP6dGRkckyTGZtbXFbBR9l5fJB1htUZOew+lg9Ax7bhXD88ieXzH2LTqlc1nEJEvKZpnC1FRKRO7SLD+FZ6N76V3o0zBSWs3p/Du3tPsWTbMZ7fdITgIGNAl2jSk2NIT45hWNf2tI8K83fZAOScL2bPsXy2HMpl88Fcdh49R2lFJUEGw7q25/4rezOudwf6dm6DmbHyl4cVkEXEqxSSRUSaobjW4cxIS2JGWhLFZRVsPlgVPj88+BnPbjjEn9YdACAhOoJ+CW3p17ktfTu3JTk+iqT2kV654uycI7ewlEOfFXEkt5DMUwXsPZ7PnuP5nCkoASAkyBiYGM0to7ozonsMI7rHEB0Z2ui1iIjURSFZRKSZiwgNZkyveMb0igeqhjLsOJLHzqN57D2Rz97j+azen0Ol+7/3xEaFkRgTSZd2EbSLDKN9ZCjtI8NoFxlGq9BgQoKN0GAjNDiIIDNKyisoKaukpLyS4rIKzn1exmeFpZwpKOGzglJOny/hSG4RBSXl/zxGSJCR2rENV/SOp1/ntvRPaMvAxGgiw/TRJCL+pzORiEgLExEazKU9Yrm0R+w/f1ZcVsEnpwo4nFvIkdwisnM/Jzu3iIyT58krKiPv8zIqLkzRHggLCSK+dTixrcPoFB3BiO7t6RYbRbfYSLrFRpEU04rwEC2KIiKBSSFZRESICA1mYGI0AxOja329stJxvqScs4WllJRXUlZR9ae80lFR6QgPCSIiNJjwkCDCQ4OJbhVKVFiwZpwQkSZLIVlEROoUFGREtwolupXGB4tIy9DylmcSEREREamDQrKIiIiISA0KySIiIiIiNSgki4iIiIjUoJAsIiIiIlKDQrKIiIiISA0KySIiIiIiNSgki4iIiIjUoJAsIiIiIlJDg1bcM7PFQO/qp+2APOfckAZXJSIiIiLiRw0Kyc65b/zjsZnNBc41uCIRERERET9rUEj+BzMzYAYwvjH2JyIiIiLiT40SkoHRwCnn3CdftIGZ3QHcUf2YtLS0Rjq0NBf79u1TX8i/UF9IbdQXUhv1hTQmc85dfAOzd4FOtbw0xzm3rHqbJ4Es59xcTw4aFRXlCgsL61urNHNpaWls3brV32VIgFFfSG3UF1Ib9YXUxsy2Oefq/dtTnVeSnXMT6zhwCDAdGF7fg4uIiIiIBKLGGG4xEdjvnDvq6RuKiooKzCyjEY4tzUucmZ3xdxEScNQXUhv1hdRGfSG16V33Jv+qMULyN4GX6vmejC9z2VuaNzPbqr6QmtQXUhv1hdRGfSG1MbMvNQanwSHZOTerofsQEREREQkkWnFPRERERKQGf4XkP/npuBLY1BdSG/WF1EZ9IbVRX0htvlRf1DkFnIiIiIhIS6PhFiIiIiIiNSgki4iIiIjU4NWQbGaTzCzDzLLM7IFaXg83s8XVr39oZt29WY8EBg/64vtmttfMdprZe2bWzR91im/V1RcXbHedmTkz0zRPzZwnPWFmM6rPF3vM7EVf1yi+58FnSFczW2NmO6o/R67yR53iW2a20MxyzGz3F7xuZvZ4dd/sNLNhde3TayHZzIKBJ4DJQD/gBjPrV2OzW4GzzrmewKPAf3urHgkMHvbFDiDNOTcIeBX4jW+rFF/zsC8wszbAPcCHvq1QfM2TnjCzVOA/gFHOuf7A93xeqPiUh+eKnwIvO+eGUrWWwx98W6X4ybPApIu8PhlIrf5zB/BkXTv05pXkkUCWc+6Ac64U+CswrcY204BF1Y9fBSaYmXmxJvG/OvvCObfGOVdU/XQTkOjjGsX3PDlfAPyCql+ain1ZnPiFJz1xO/CEc+4sgHMux8c1iu950hcOaFv9OBo47sP6xE+cc+uA3ItsMg34i6uyCWhnZp0vtk9vhuQuQPYFz49W/6zWbZxz5cA5INaLNYn/edIXF7oVWOnViiQQ1NkXZjYUSHLOve7LwsRvPDlX9AJ6mdkGM9tkZhe7iiTNgyd98TAw08yOAm8Cs31TmgS4+uaPRlmW+ovUdkW45nxznmwjzYvH/+ZmNhNIA8Z6tSIJBBftCzMLompI1ixfFSR+58m5IoSqW6dXUHXH6X0zG+Ccy/NybeI/nvTFDcCzzrm5ZnYp8Fx1X1R6vzwJYPXOnN68knwUSLrgeSL/esvjn9uYWQhVt0Uudqlcmj5P+gIzmwjMAa5xzpX4qDbxn7r6og0wAFhrZoeAS4Dl+vJes+bpZ8gy51yZc+4gkEFVaJbmy5O+uBV4GcA5txGIAOJ8Up0EMo/yx4W8GZK3AKlmlmxmYVQNnl9eY5vlwM3Vj68DVjutbtLc1dkX1bfVn6IqIGuMYctw0b5wzp1zzsU557o757pTNVb9GufcVv+UKz7gyWfIa8A4ADOLo2r4xQGfVim+5klfHAEmAJhZX6pC8mmfVimBaDlwU/UsF5cA55xzJy72Bq8Nt3DOlZvZ3cDbQDCw0Dm3x8x+Dmx1zi0HnqbqNkgWVVeQv+mteiQweNgXvwVaA69Uf4/ziHPuGr8VLV7nYV9IC+JhT7wNfNXM9gIVwP3Ouc/8V7V4m4d98QNggZndR9Xt9Fm6ANf8mdlLVA29iqsej/4QEArgnPsjVePTrwKygCLgljr3qb4REREREfn/tOKeiIiIiEgNCskiIiIiIjUoJIuIiIiI1KCQLCIiIiJSg0KyiIiIiEgNCskiIiIiIjUoJIuIiIiI1PC/4SIznYTd2NYAAAAASUVORK5CYII=\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",
    "\n",
    "fig.savefig('fem3_a.png',dpi=300)                              # save figure as png (dpi = number of pixels)"
   ]
  }
 ],
 "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
}