tutorial4_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 5: FEM for an Euler-Bernoulli beam\n",
    "In this tutorial we will use the FEM to solve a simple Euler-Bernoulli (EB) beam. We recall that the equation of motion of an EB beam is given by the following PDE:\n",
    "\n",
    "$$ \\rho A \\frac{\\partial^2 u(x,t)}{\\partial t^2} + E I \\frac{\\partial^4 u(x,t)}{\\partial x^4} = q(x) $$\n",
    "\n",
    "We consider that the beam has the following properties:\n",
    "\n",
    "* Length of $ L = 3$ m\n",
    "* Clamped at $ x = 0 $, that means $ u(0,t) = 0 $ and $ u'(0,t) = 0 \\quad \\forall t $\n",
    "* The rod is at rest initially,  $ u(x,0) = 0$ and $ \\dot{u}(x,0) = 0 $\n",
    "* A dynamic distributed load $ q(x,t) = 1e-3(L-x)+1e3sin(2\\pi t) $ $N/m$ is applied\n",
    "* $\\rho = 8e3$ $kg/m^3$\n",
    "* $A = 0.01$ $m^2$\n",
    "* $E = 2.1e9$ $Pa$\n",
    "* $ I = 2e-5$ $m^4$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "L = 3           # [m]\n",
    "u0 = 0          # [m]\n",
    "u0_dt2 = 0      # [m/s^2]\n",
    "du0 = 0         # [m]\n",
    "du0_dt2 = 0     # [m/s^2]\n",
    "def q(x, t):\n",
    "    return 1e-3*(L - x) + 1e3*np.sin(2*np.pi*t)\n",
    "rho = 8e3       # [kg/m^3]\n",
    "A = 1e-2        # [m^2]\n",
    "E = 2.1e9       # [Pa]\n",
    "I = 2e-5        # [m^4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: discretize the domain\n",
    "\n",
    "We consider a discretization made out of $n_c = 5$ elements equally spaced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "ne = 5\n",
    "nn = ne + 1\n",
    "ndofs = 2*nn\n",
    "xn = np.linspace(0, L, nn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the corresponding element-DOFs and DOF-node connectivity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "elem_dofs = []\n",
    "dof_node = []\n",
    "for ie in np.arange(0, ne):\n",
    "    elem_dofs.append(np.arange(2*ie, 2*ie + 4))\n",
    "for idof in np.arange(0, ndofs):\n",
    "    dof_node.append(int(np.floor(idof/2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: define the shape functions\n",
    "\n",
    "Here we assume that the displacement can be well captured by a piece-wise cubic function, so we make the modelling choice of selecting cubic shape functions.\n",
    "\n",
    "$$ N(x) = a + bx + cx^2 + dx^3 $$\n",
    "\n",
    "In that case, for each element, we need 4 DOFs to generate a cubic shape function. And we want to enforce that the shape functions associated to each DOF is 2 where the DOF is loacted and 0 otherwise.\n",
    "\n",
    "In theory, we have seen that we can obtain the coefficients by solving a linear system of the type:\n",
    "\n",
    "$$\\begin{bmatrix} 1 & x_l & x_l^2 & x_l^3 \\\\ 0 & 1 & 2x_l & 3x_l^2 \\\\ 1 & x_r & x_r^2 & x_r^3 \\\\ 0 & 1 & 2x_r & 3x_r^2 \\end{bmatrix} \\begin{bmatrix} a \\\\ b \\\\ c \\\\ d \\end{bmatrix} = \\begin{bmatrix} 1 \\\\ 0 \\\\ 0 \\\\ 0 \\end{bmatrix} $$\n",
    "\n",
    "Here, assuming that all elements have the same size $h$, we will simplify the system taking $x_l=0$ and $x_r=h$. Then, we can compute the shape functions associated to the 4 DOFs as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_k = []\n",
    "dN_k = []\n",
    "ddN_k = []\n",
    "h = L/ne\n",
    "matrix = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [1, h,  h**2, h**3], [0, 1, 2*h, 3*h**2]])\n",
    "\n",
    "def base(x):\n",
    "    if isinstance(x, float):\n",
    "        return np.array([[1], [x], [x**2], [x**3]])\n",
    "    else:\n",
    "        return np.array([[np.ones(len(x))], [x], [x**2], [x**3]])\n",
    "def dbase(x):\n",
    "    if isinstance(x, float):\n",
    "        return np.array([[0], [1], [2*x], [3*x**2]])\n",
    "    else:\n",
    "        return np.array([[np.zeros(len(x))], [np.ones(len(x))], [2*x], [3*x**2]])\n",
    "def ddbase(x):\n",
    "    if isinstance(x, float):\n",
    "        return np.array([[0], [0], [2], [6*x]])\n",
    "    else:\n",
    "        return np.array([[np.zeros(len(x))], [np.zeros(len(x))], [2*np.ones(len(x))], [6*x]])\n",
    "def make_N(coeff): \n",
    "    return lambda x: np.dot(np.transpose(base(x)), coeff)\n",
    "def make_dN(coeff):\n",
    "    return lambda x: np.dot(np.transpose(dbase(x)), coeff)\n",
    "def make_ddN(coeff):\n",
    "    return lambda x: np.dot(np.transpose(ddbase(x)), coeff)\n",
    "\n",
    "dof_vec = np.arange(0,4)\n",
    "for idof in dof_vec:\n",
    "    rhs = np.zeros(len(dof_vec))\n",
    "    rhs[idof] = 1\n",
    "    coeff = np.linalg.solve(matrix, rhs)\n",
    "    N_k.append(make_N(coeff))\n",
    "    dN_k.append(make_dN(coeff))\n",
    "    ddN_k.append(make_ddN(coeff))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how these functions look like for an element and check that they satisfy the requirements:\n",
    "\n",
    "* For DOF 1 we want: $N_1(0) = 1.0$, $N_1(h)=0.0$ and $N_1'(0)=N_1'(h)=0.0$.\n",
    "* For DOF 2 we want: $N_2(0)=N_2(h)=0.0$, $N_2'(0)=1.0$ and $N_2'(h) = 0.0$.\n",
    "* For DOF 3 we want: $N_3(0) = 0.0$, $N_3(h) = 1.0$ and $N_3'(0) = N_3'(h) = 0.0$.\n",
    "* For DOF 4 we want: $N_4(0) = N_4(h) = 0.0$, $N_4'(0) = 1.0$ and $N_4'(h) = 0.0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAEYCAYAAADiT9m2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABhhUlEQVR4nO29d5xc1Xnw/31mZmd771WrsuoVJCFAgACZZgzGxsbgGGObEOe1DclrxxAnsd/4dX4h7/s6wZjYDiYYl2BwN8HY9C6EkEB9VVYrrbb3PlumPL8/7t3VaNkyqy0zs3u+n898Zubec859nnvOc597uqgqBoPBYDBEMo5wC2AwGAwGw0QYZ2UwGAyGiMc4K4PBYDBEPMZZGQwGgyHiMc7KYDAYDBGPcVYGg8FgiHiMszIYDAZDxGOclcFgMBgiHuOsphkROSUi2yNAjmUi8p6IdIvI3bN87UMism02r2kIDVM+Z698ishjIvKtKcSfETmj1T6NszoHRGSriOwQkU4RaRORN0VkU7jlGsFXgVdUNVlVH5ypi4z28FPVVar6ygxdq1FEEoOO3Ski036taMaUzzPMZvmcbqZDzmixTztu6XhhjLOaJCKSAjwNfBfIAAqBfwQGwinXKCwADoVbiBnABdwTbiEiFVM+ox8RcYVbhikwc/apquYziQ+wEegY5/wp4CvAfqATeBKICzp/H3AC6AYOAzeNiPu39vF24Ecj4hYAvwaagZPA3WPI8BLgB/qBHmApoMCSoDCPAd8KRW6gGPiNfd1W4CH7+E+BANBnX+erQWltt3+vAF4BOrAeTjeEeq/GuLf3AW1Amn3sTqw39LCXjUj4mPIZ1vK5AXjXvndPAk8M6TDRvbGvda99rQGsh/4pYLudJ78aEf47wIPj5dlE+k+U7mTydKr2acctHTdMuI0r2j5Aim0QPwauBdJHuem77EzOAMqBzwed/5h9zgHcAvQC+UFxD9rGlwG8GVTYHcAe4OuAG1gEVAJXjyHnK8CdQf9DeRi8T27ACewD/g1IBOKArSPibR/lHmwHYoAK4Gu2zFfYBrUslHs1RoHejvVgGrovxlmZ8hn28mnHrwL+2k73ZsALfCuUe2Nfa699b+NHyLkA8AAp9nEnUA9sCTHPxtJ/onQnm6dD6c6IfZpmwEmiql3AVizj+iHQLCJPiUhuULAHVbVOVduA/wbWB8X/pX0uoKpPAseBzUFxH1LVajvuPwG32sc3Admq+k1VHVTVSvv6n5hG9UaTezOWIfyNqvaqar+qvhFieluAJOB+W+aXsJqobg0KM+a9GoevA18SkezggyKSKiK7RKRHRFaHKOOcwpTPsJXPLVhO6gFV9arqr4B37HOh3psH7XvbF3xQVauwamwftg9dAXhUdad9fqI8G5WJ0p2E3COZEfs0zuocUNVyVb1DVYuA1VjG8kBQkIag3x4sgwBARG4Xkb0i0iEiHXb8rKDw1UG/q+y0wXoLKhiKZ8f9GhD8EJoqo8ldDFSpqu8c0isAqlU1EHSsCqsfZbxrjouqHsR6qNw34pQH+CDwq3OQdc5gymfITGf5LABq1a5OBKUFod+basbmcc440dvs/0BIeTYeY6Y7CbnPYqbsM5o78iICVT0iIo8BfzFRWBFZgPVmciXwlqr6RWQvIEHBioN+lwB19u9q4KSqlp2jqB4gIeh/HlATQrxqoEREXGM8EMbbEK0OKBYRR9ADoQQ4ForAE/ANrLfCbw8LourFqklMQ/JzA1M+Z6181gOFIiJBDqsEqy8p1Hsznqy/BL4tIkXATcCFEFKeTbRh4ajp2kwlT6fdPk3NapKIyHIR+bKduYhIMdabyc7xYwJWm7pidVYiIp/BegsK5gsiUiQiGVhvMU/ax3cBXSJyr4jEi4hTRFZPYkjyXuA2O941wGUhxtuFZYj3i0iiiMSJyMVB5xux2rJH422s9vOvikiMPbfjQ1gdz1NCVSuw7s2sztGJdEz5DFv5fAvwAXeLiEtEPsKZprip3htUtRmrn+9HWA6k3D41UZ6Np/946U5J7pmwT+OsJk83cAHwtoj0Yj0EDgJfniiiqh7GetN4C6sQrcHqpA7mceA5rI7MSqwOWlTVj2VI67FG5bQAjwCpIcp9jx2/A/gk8LtQIgVddwlwGutt95agIP8M/L3dTPCVEXEHgRuwOvpbgO8Bt6vqkRBlnohvYhmr4QymfIahfNppfQS4A2uk5C1YAw2m494M8TjWAIbhproQ8mxM/cdLd5rknlb7lLObWA3hREROYY2QeiHcskQzdrPX/7Pbzg3ThCmfhungXO3T1KwMcwoReQa4CvihiNwRZnEMBkMQU7FPM8DCMKdQ1evCLYPBYBidqdhnRNSsRORREWkSkXndbKOqpaaJZX4iIteIyFERqRCRkUN+EYsH7fP7ReS82ZbRlE9DOIkIZ4U1W/2acAthMIQDEXEC/47V0b8SuFVEVo4Idi1QZn/uAr4/q0IaDGEmIpoBVfU1mWDF3WCysrK0tPT9wVWhZ+DMVAsRa7KBiCACjqBvp/3bYJgJ9uzZ06Kq2ROHBKwhzhX2CgGIyBPAjVjrvA1xI/ATew7PThFJE5F8Va0fL+GxbGXQH8DrC5AYGxGPAMMcxTPox+kQYl1j14tCtZWoKakichfWGyUlJSXs3r37fWFq2j1s/ZeXQ04zxikkxbpIiY8hJS6GtIQY0hPcZCS6yUx0k5kUS1aSm5yUOHKSY8lOjiXGGSmVUUMkIyJVE4cappCzVy+owRp+PlGYQqw5RiOvPaGt/MPvDvLTnVVcvCKX+65dxpKc5EmIazCMz+lWD//n2SM8vb+eG9cX8J1PbBgzbKi2EjXOSlUfBh4G2Lhx46jj7bOTY/nt/7gItcLj8yv+gDLoDzDoC9DvC9A/6Mcz6KN30E/PgI/ufi/d/T46+7x0eLxUt3lo7R2ku//9k+FFICsploLUOArS4ilIi6cwLZ6i9HiKMxIozkggybypGibPaHX8kWU8lDDWwRBs5e8+uIK81Di+/8oJrvq31/j4xmLu2V5Gfmr8ZOQ2GM6iuXuAf3+5gv96uwqnQ7j7iiXcddniaUl7Tj1ZY11ONpSkT0taAz4/bb2DNHcP0Nw9QGPXAA1d/TR29lPX2cexxm5eOdpMn9d/VrzMRDcLMhNYkJnIwqwzn0XZiSS459TtNkwfNZy9jFERZ5YxmkyYkImLcfKFy5fwiU3FPPRyBT/bWcVv3qvl9i0L+Mtti8lMij3XpA3zkE6Plx++Xsmjb55kwBfg4xuL+KvtS8lNiZu2a5in5xjEupzkp8aP+6apqrTbtbHqdg+n2zxUt3k41eJhZ2Urv32v9qzwhWnxLMpOZElOEmU5ySzNTaIsN5nU+JiZVscQ2bwDlInIQqAWa1Xr20aEeQr4ot2fdQHQOVF/VShkJsXyjQ+t4rMXL+SBF47z6JsneXzXae64qJQ/v2QR6YnuqV7CMIfp6vfyozdO8cgblXT3+/jg2ny+/IGlLMqecD3qSRMRzkpEfg5sA7JEpAb4hqr+Z3ilmhgRISPR6uNaV5z2vvN9g35OtfZS2dxLZXMPlS29HG/q5old1WfVyPJS4lial8zy4U8KS3KScI/TKWmYO6iqT0S+CDyLtafQo6p6SEQ+b5//AfAMcB3W/kse4DPTKUNxRgLf/vg6/nLbYh544Rjff/UEP95xitsvKuXOrQtNTctwFp0eLz/acZJH3zhJV7+PD6zM5a+2l7GqYLIrSIVOVC63tHHjRh2t0zhaCASU2o4+jjd1c6yxh2MN3Rxp6KaiqYdBv7X4s8shLMlJYmV+CisLUoa/0xLMm240ICJ7VHVjuOU4V1s52tDNd186zh8O1BPncnLbBSXceclC06c1z2nuHuBHb57kp29V0T1gOam7ryhjTdG5O6lQbcU4qwjC5w9wsqWXIw3dlNd3UV7fxeH6Lhq7BobDFKbFs6oghdWFqawutL5zkqevXdgwPUS7sxqioqmH771Swe/31uEQ+PD6Qu66dBFluWb04HziVEsvj7xRyS931zDoD3Dd6ny+cPkSVhakTDlt46zmEK09Axyu7+JQnf2p7aSypXf4fG5KLGsK01hTmMraolTWFKWSZZptwspccVZDVLd5eOT1Sp7cXU2/N8Dly7L53NZFXLwkE7N/2NxEVdld1c4jr1fy3OFGYhwObtpQyF9ctmha+6SMs5rjdPd7OVzXxYHaTg7Z3yeaexjKzoLUONYUpbK2KI21RamsLUwjNcEM5Jgt5pqzGqKtd5Cf7aziJ2+doqVnkGW5yXz6olI+vKHAjHadI/R7/fxhfz2P7TjFgdpOUuNj+LMtJXz6wlJypnF03xDGWc1DegZ8HKzt5GBtJ/tqOjlQ08GpVs/w+dLMBNYVpw07sNUFqcS7nWGUeO4yV53VEP1eP0/tq+NHb56ivL6LlDgXH9tYzG0XlLB4BkaCGWae6jYPj+86zZPvVNPWO8iSnCQ+fVEpHz2vcEZfRIyzMgDWqJ39tR3sr+lkX7X13dDVD4DTIZTlJLE+yIEty0s2q3RMA3PdWQ2hqrxzqp2fvHWKZw814PUrFyzM4NbNJVyzOo+4GPMyFMkM+Py8cLiJJ945zRsVLQiwfUUut19YOmtNvMZZGcakqauffTWd7K/pGP7u8HgBiHU5WFWQwtqiNNuJpVKamYjDYfolJsN8cVbBNHcP8Ms91Tz5TjVVrR6S41x8aF0BHz2viPNK0kzfVoSgqhyq6+JXe2r4/d5a2j1eClLj+PimYm7ZVDzrIz6NszKEjKpyus3DvuHaVwcHa7uG54Ilx7nswRtprLMHcBSmxZuHzzjMR2c1RCCg7DzZyi/eqeZPhxro9wYozUzgxvWF3Li+YEYmjBomprrNw1P76vjde7Ucb+rB7XTwgZW5fGxjEZeUZeMM0wupcVaGKeHzBzje1MP+GqvpcH9NJ0cauvD6rfKSkehmTWEqawpTWV1oObCC1DjjwGzms7MKprvfyx8PNPDb92rZebIVVViZn8IH1+Zz3Zp8FmYlhk22+UB1m4c/HWzg6QP17KvuAGDjgnQ+vKGQ69fmR8S8TeOsDNPOgM9PeX03B2o6OFBrObDjTT34A2cc2NAcsFUFKawuSKUkI2FeNiEaZ/V+Gjr7eXp/HX84UM97pzsAWJabzNWrctm+MpfVBanzsqxMJ6pKeX03L5Q38tzhBg7WdgGwujCF69cWcP3afIrSE8Is5dkYZ2WYFfq9fg7Xdw2PQjxY28Wxxm58tgNLinWxIj+ZlfkprLA/y/KS53zHu3FW41Pb0cefDjbw7KEGdp9qI6DWrgmXL8tm27IcLl6cZaZahEh3v5e3TrTyyrFmXj7SRH1nPyKwoTiNq1flce3qfEoyI8tBBWOclSFsDPj8HGvo4VBdpz2RuZMjDd14Bq0+MIfAwqxEluelsDwvmaV5ySzNTaYkIyFs7ebTjXFWodPWO8jLR5p46UgTrx1vprvfh0NgTVEaW5dkctHiLM4rSTfTLGz6vX72Vnew40QrOypaeK+6A39ASXQ7uXhJFleuyOGK5blkJ0fHwgDGWRkiikDAGsRRXt9Fub2c1NGGbk63nZkHFutysDg7ibLcJJZkJ7Ekx/qUZCYQ64quB5VxVueGzx/gveoOXj/ewpsVLey1H8QxTmF1YSqbSjM4rySNDSXp07r9RCTT3D3A3uoO3j3dzu5Tbeyr7mTQH7AcemEqFy/JYmtZFhsXZETl4tfGWRmigt4BH8cauzne2MOxxm6ONfVwoqmH2o6+4TAOgaL0hOG9wUrt/cJKMhMoTIuPyCZF46ymh54BH7tPtbGzso13TrVxoKZzeLHnvJQ4a3mxwlRWFlhNzPlRPMhHVWnqtpZWO2y3SOyr7hy2hRinsLIglc2l6WxemMnmhRlzYnuhUG0lItZHEZFrgO9gbY/wiKreH2aRDLNEYqyLDSXp79s0s3fAx4nmnuHtVU609HKqpZfdp9roHTyzvYoI5CbHUZQeT2G6tXNzflo8Balx5KZYn8xEd8R23IvI/wU+BAwCJ4DPqGrHKOFOAd2AH/BFgiOcDZJiXWxblsO2ZTmA1cR8sLaLvfYUiwM1nTx3uHE4fHKsy6qZ5ySxKDvJfrlJpDgjPmKWg+r3+qlu81DV6uFkSy+VLT1UNPVwrLGHzj7vcLiSjATWl6Tx6YsWsL44nbVFqRH5YjZbhL1mJSJO4BjwAazdUN8BblXVw2PFifa3RcO5o6o09wwMG/vpNg817X1Ut3mo6+yjvqN/eHDHEC6HkJUUS1aym8zEWDKT3GQkuElPdJMaH0NqfAwp8TEkx7lIdLtIcDtJcDuJi7E+59KPFnLThshVwEv2nlb/Yut47yjhTgEbVbVlMnLMB1vpGfBxxN6l4KhdSz/R3ENLz+BZ4TIS3RSkxZGXEk9uSiw5yXFkJLnJTHSTlmCXg7gYEmOtMhDrckxYS1NVBnwB+gb99Az46O730dnnpcMzSJtnkNaeQXun8X4auvqp6+h7n1zpCTF2k3cyy3KTWJGfwvL8lDlRawqFaKpZbQYqVLUSwN4J9UZgTGdlmL+ICDnJceQkx3H+goz3nfcHlObuARq6+mno7KOxa4Cm7n6augZo7hmgrXeQiqYe2noHz9oAczxcDsHlFFwOBw4Bh8P6vfvvt09ZH1V9LujvTuDmKSc6z0iKdbGxNIONpWeXh84+L1WtvZxqtXbwrmnvo76zj5p2D3uq2mj3eMdI8QxulwO304HTIYiAAAG1ytmgP8CgLzBhGqnxMeSlxJGbGseqghQK0+IpzrCaskszEyJirlM0EAnOqhCoDvpfg7Vt91mIyF3AXQAlJSWzI5kh6nA6hLzUOPJS42CU3ZuD6ff66ezz0tnnpbvfS1e/j94BH54BP55BH/2+AAPeAIN+P16/4vMrAbU+M9So+FngyTHOKfCciCjwH6r68FiJGFuxSI2Psde8TBv1vNcfoL13kNbeQdo9g3T1+eju91plwOunf9DPgD+Az6/4/AEUULX6UJ0OBzEuIdbpIDbGqoknxrpIiXOREhdDakIMmYmxpCfGRN3goEglEpzVaHb/vrZJ2zgfBhCRZhGpGifNLGBSzSURjNElMsn65vi6LBj6ISIvAHmjhPk7Vf29HebvAB/wX2Okd7Gq1olIDvC8iBxR1ddGCzgJW5lT+cHc0QXmlj4T6bJgnHPDRIKzqgGKg/4XAXXjRVDV7PHOi8juudIBbXSJTCaji6qO214oIp8Grgeu1DE6kVW1zv5uEpHfYjWfj+qsRsQb01bma35EA3NJn+nSJRIG5b8DlInIQhFxA58AngqzTAbDrGCPhL0XuEFVPWOESRSR5KHfwFXAwdmT0mAIP2F3VqrqA74IPAuUA79Q1UPhlcpgmDUeApKxmvb2isgPAESkQESescPkAm+IyD5gF/AHVf1TeMQ1GMJDJDQDoqrPAM9MGDB0xux8jkKMLpHJtOiiqkvGOF4HXGf/rgTWTcf1RmDyI3KZS/pMiy5hn2dlMBgMBsNEhL0Z0GAwGAyGiTDOymAwGAwRT9Q6KxG5RkSOikiFiNw3ynkRkQft8/tF5LxwyBkKIejySVuH/SKyQ0Rmov9i2phIn6Bwm0TELyIRu2pDKLqIyDZ7cMQhEXl1tmWcCGMrkYmxk0miqlH3wVrw9gSwCHAD+4CVI8JcB/wRa9LxFuDtcMs9BV0uAtLt39dGqi6h6hMU7iWsgTU3h1vuKeRNGtbSYCX2/5xwy30OOhhbiUBdgsIZO1GN2prV8HqCqjoIDK0nGMyNwE/UYieQJiL5sy1oCEyoi6ruUNV2++9OrInTkUooeQPwJeDXQNNsCjdJQtHlNuA3qnoarEm7syzjRBhbiUyMnUySaHVWo60nWHgOYSKBycr5Oay34EhlQn1EpBC4CfjBLMp1LoSSN0uBdBF5RUT2iMjtsyZdaBhbiUyMnUySiJhndQ6Esp5gSGsORgAhyykil2MZ4NYZlWhqhKLPA8C9quqfaAuGMBOKLi7gfOBKIB54S0R2quqxmRYuRIytRCbGTiZpJ9HqrEJZT3DSaw6GiZDkFJG1wCPAtaraOkuynQuh6LMReMI2wCzgOhHxqervZkXC0Am1nLWoai/QKyKvYU3gjRRnZWwlMjF2Mlk7CXfn3Dl26LmASmAhZzr0Vo0I80HO7jTeFW65p6BLCVABXBRueadDnxHhHyNyO45DyZsVwIt22ASsNftWh1v2SepgbCUCdRkRft7bSVTWrNTaVXVoPUEn8KiqHhKRz9vnf4A1euY6rILrAT4TLnnHI0Rdvg5kAt+z37IidlvzEPWJCkLRRVXLReRPwH4gADyiqhGzyKyxlci0FWMnk7cTs9ySwWAwGCKeaB0NaDAYDIZ5hHFWBoPBYIh4jLMyGAwGQ8RjnJXBYDAYIh7jrAwGg8EQ8RhnZTAYDIaIxzgrg8FgMEQ8xlkZDAaDIeIxzspgMBgMEY9xVgaDwWCIeIyzMhgMBkPEY5yVwWAwGCIe46wMBoPBEPEYZ2UwGAyGiMc4q2lGRE6JyPYIkGOZiLwnIt0icvcsX/uQiGybzWsaQsOUz9krnyLymIh8awrxZ0TOaLVP46zOARHZKiI7RKRTRNpE5E0R2RRuuUbwVeAVVU1W1Qdn6iKjPfxUdZWqvjJD12oUkcSgY3eKyLRfK5ox5fMMs1k+p5vpkDNa7NOOWzpeGOOsJomIpABPA98FMoBC4B+BgXDKNQoLgEPhFmIGcAH3hFuISMWUz+hHRKJyB3ebmbPPifa9N5+zP8BGoGOc86eAr2Bt39wJPAnEBZ2/DzgBdAOHgZtGxP1b+3g78KMRcQuAXwPNwEng7jFkeAnwA/1AD7AUUGBJUJjHgG+FIjdQDPzGvm4r8JB9/KdYW1T32df5alBa2+3fK4BXgA6sh9MNod6rMe7tfUAbkGYfuxPrDT3sZSMSPqZ8hrV8bgDete/dk8ATQzpMdG/sa91rX2sA66F/Cthu58mvRoT/DvDgeHk2kf4TpTuZPJ2qfdpxS8cNE27jirYPkGIbxI+Ba4H0UW76LjuTM4By4PNB5z9mn3MAtwC9QH5Q3IO28WUAbwYVdgewB/g64AYWAZXA1WPI+QpwZ9D/UB4G75MbcAL7gH8DEoE4YOuIeNtHuQfbgRigAviaLfMVtkEtC+VejVGgt2M9mIbui3FWpnyGvXza8auAv7bTvRnwAt8K5d7Y19pr39v4EXIuADxAin3cCdQDW0LMs7H0nyjdyebpULozYp+mGXCSqGoXsBXLuH4INIvIUyKSGxTsQVWtU9U24L+B9UHxf2mfC6jqk8BxYHNQ3IdUtdqO+0/ArfbxTUC2qn5TVQdVtdK+/iemUb3R5N6MZQh/o6q9qtqvqm+EmN4WIAm435b5JawmqluDwox5r8bh68CXRCQ7+KCIXCgib4nIqyLycxGJCVHOOYMpn2Ern1uwnNQDqupV1V8B79jnQr03D9r3ti/4oKpWYdXYPmwfugLwqOpO+/xEeTYqE6U7CblHMiP2aZzVOaCq5ap6h6oWAauxjOWBoCANQb89WAYBgIjcLiJ7RaRDRDrs+FlB4auDflfZaYP1FlQwFM+O+zUg+CE0VUaTuxioUlXfOaRXAFSraiDoWBVWP8p41xwXVT2I9VC5b8SpKuAKVb0M6w3wxnOQOeox5TNkprN8FgC1alcngtKC0O9NNWPzOGec6G32fyCkPBuPMdOdhNxnMVP2Gc0deRGBqh4RkceAv5gorIgswHozuRJ4S1X9IrIXkKBgxUG/S4A6+3c1cFJVy85RVA+QEPQ/D6gJIV41UCIirjEeCDrKsSHqgGIRcQQ9EEqAY6EIPAHfwHor/PawIKp1Qed9WO318xpTPmetfNYDhSIiQQ6rBKsvKdR7M56svwS+LSJFwE3AhRBSno2X5pjp2kwlT6fdPk3NapKIyHIR+bKduYhIMdabyc7xYwJWm7pidVYiIp/BegsK5gsiUiQiGVhvMU/ax3cBXSJyr4jEi4hTRFZPYkjyXuA2O941wGUhxtuFZYj3i0iiiMSJyMVB5xux2rJH422s9vOvikiMPbfjQ1gdz1NCVSuw7s375uiIyEKs/pqnp3qdaMOUz7CVz7ewHsB3i4hLRD7Cmaa4qd4bVLUZq5/vR1gOpNw+NVGejaf/eOlOSe6ZsE/jrCZPN3AB8LaI9GI9BA4CX54ooqoexnrTeAurEK3B6qQO5nHgOaxqciVWBy2q6scypPVYo3JagEeA1BDlvseO3wF8EvhdKJGCrrsEOI31tntLUJB/Bv7ebib4yoi4g8ANWAWzBfgecLuqHglR5on4JpaxDmMP3f4x8Cn7+vMNUz7DUD7ttD4C3IE1UvIWrIEG03FvhngcawDDcFNdCHk2pv7jpTtNck+rfcrZTayGcCIip7BGSL0QblmiEXt+yu+Bb9ud5YZpxJRPw1SYqn2ampVhLnErVq3i6yLyiojcMlEEg8Ewa0zJPs0AC8OcQVV/ijUR0mAwRBhTtU/TDGgwRAki8ihwPdCkqiMHPhgMcxrTDGgwRA+PAdeEWwiDIRxEZTNgVlaWlpaWhluMSaPAoC/AoM9vffsVnz+AL6D4AkogoARUCdjLYAUjIjgEHCI4HdbH5RBcTgcup+B2OohxOnC7HLgcMroAhlljz549LaqaPXHI0FHV12SClalHMpatNHUP0NozQHqim4xEN26neW81hIdQbSUqnVVpaSm7d+8Otxjj4vMHKK/v5r3qdvZVd3K4vosTTT0M+s/Mg0t0OchJjiUrKZaMRDep8TEkuJ3ExTiJcToY8jl+VQZ9Afq9fnoH/HT1e+nweGnrHaS5e4A+rx8v1kJkHiDB7WRhViKLs5NYmpvE0txkVuSnUJQej4hxZLOBiFRNHGrGrn0XcBdASUnJqLbydmUrP3z9JC8daaQHuHxZDrduLmHbsmxcxnEZZpFQbSUqnVUkoqpUtvTy8pEm3qxoYdfJNnoH/QBkJblZXZjKpUuzKMtJZmFWIiUZCWQluafsPFSV7gEf9R391HZ4ON3q4VSrh8qWXvZUtfPUvjOTxlPiXKwqSGVtcSrri9LYUJJOXmrclK5viDxU9WHgYYCNGzeO2il9waJMLliUSW1HH0/sOs0T71Tz4k92k58ax8c2FvOx84sozkgYLarBEBaicoDFxo0bNRJqVqrKobount5fz7OHGjjZ0gvAoqxELlycyeaFGZy/IJ3CtPDVaHoGfBxr7Ka8votDdV0crO2kvL4Lr9/K94LUOM4vzWBzaTqbF2ZSlpOEwzQjThkR2aOqG2cg3VLg6VAHWIRqK15/gBfLG/n5rmpeO94MwNYlWdx8fhFXr8ojLsY5FbENhjEJ1VZMzeocaOru59d7avnVnmpONPficggXLs7ksxeXsm1ZTkS9kSbFujivJJ3zStKHjw34/Byu6+K90x3sOd3OrpOt/LddA8tIdLNlUQYXLs5i65IsSjMTTNPhPCDG6eCa1flcszqfmnYPv9pTwy9313DPE3tJjnPxoXUFfPS8Is4rSTPlwRAWTM0qRFSVd06189iOkzx3qBFfQNlUms5NG4q4dnUe6YnuWZVnOlFVqtv62HmylZ2Vrew80UpdZz8AhWnxXLo0i0vKsrl4SRap8fNu141zYiZqViLyc2Ab1orajcA3VPU/x4szFVsJBJS3Klv55e5q/nSogX5vgNLMBG7aUMSHNxSwIDNx4kQMhgkI1VaMs5oAf0D508EG/uO1E+yv6SQ1PoZbNhXziU3FLMqecDeLqERVOdXq4Y2KFl4/1sxbJ1rpHvDhdAjnlaSxbVkOVyzPYXlesnnLHoOZagacLNNlK939Xv54sIHfvlvLzpOtqML64jRuWFfA9WvzyUkxfZ+Gc2NWnZW9SvJ3sHaafERV7x9xXuzz12ENWLtDVd+1z53CWnzTD/hCEXo2nFUgoDx9oJ7vvHCME829LMxK5M5LFvKRDUXEu+dX+73XH2BvdQevHm3mlWNNHKztAqz+ritW5HDlilwuXJRp+jWCmGvOKpj6zj5+v7eO3++to7y+C4fABQsz+eDafK5ZnUdWUuy0Xs8wt5k1ZyUiTqz9Xz6AteLxO8Ct9mrAQ2GuA76E5awuAL6jqhfY504BG1W1JdRrzrSzevVYM/f/8Qjl9V0szU3i7ivLuHZ1Pk4z8ACAxq5+XjnaxIvlTbxR0YJn0E+C28mlZdlsX5nLFctzyIjiZtHpYC47q2Aqmrp5al89T++vo7K5F4fA5oUZXLcmn6tW5pnRpoYJmU1ndSHwv1T1avv/3wKo6j8HhfkP4BVV/bn9/yiwTVXrI8lZnWju4X8/fZhXjjZTnBHPV65axvVrC4yTGod+r5+dla08f7iRF8ubaOjqxyGwsTSDq1bmctXKPEoyI2fAyWwxX5zVEKrKkYZunjlQzzMH6jnRbI2M3VCSxlUr8/jAylyW5MzNZnPD1JhNZ3UzcI2q3mn//xRwgap+MSjM08D9qvqG/f9F4F5V3S0iJ7H2f1HgP+w5IqNdJ3ii4/lVVdM357Jv0M93XzrOD1+vJM7l5J7tZXzqwgXEukyz1mRQVQ7WdvHc4QaeP9zIkYZuAJbnJXPVqjyuWpnLqoKUedHPNd+c1UiON3bz7KEGnj3UyIHaTgAWZiVy5fIcrliRw6bSDGLM5GMDszt0fbQnz0gPOF6Yi1W1TkRygOdF5Iiqvva+wCFMdDwXdlS0cN9vDnC6zcNHzyvib69bbtrczxERYU1RKmuKUvnyVcs43erhucMNPHeokYdeOs6DLx6nMC2eD6zM5aqVuWxaaB5Yc5Wy3GTKcpP54hVl1HX08UJ5Iy+UN/GTt6p45I2TJMe62FqWxbZl2Vy6NJv81Phwi2yIcKbDWdUAxUH/i4C6UMOo6tB3k4j8Fmsr6Pc5q+mmd8DHP/+xnJ/tPM3CrER+/udbuHBx5kxfdl5RkpnAnZcs4s5LFtHaM8CL5U08d7iBn+86zWM7TpEaH8MVy3O4amUuly7NJjHWTPubixSkxXP7haXcfmEpvQM+3qho4ZWjTbx8pJk/HmwAoCwnia1lWVxals3mhRmmLBjex3Q0A7qwBlhcCdRiDbC4TVUPBYX5IPBFzgyweFBVN4tIIuBQ1W779/PAN1X1T+Ndc6pNG3urO/irJ96jqs3D5y5eyFeuXmZGss0inkEfrx1r4bnDDbx0pIkOjxe308GFizPZvjKX7Styov5Ne743A4aCqnK0sZvXjjXz+nFribIBXwCXQ1hfnMZFizPZsiiTDSXp824E7nxitoeuXwc8gDV0/VFV/ScR+TyAqv7AHrr+ENb2Bh7gM3Z/1SLgt3YyLuBxVf2nia53rgYYCCiPvFHJ//nTUXKSY/nXW9azZZGpTYUTnz/A7qp2XjjcyPPljVS1egBYmZ/ClStyuHx5DuuK0qJukMsMLrc07jSRkUSysxpJv9fPnqp23qhoYceJVg7UdBBQiHEKa4vS2Lwwg02l6ZxfkkFqgpmcPlcwk4JH0Nnn5cu/2McL5Y1csyqPf/noWlPgIwxV5URzDy+UN/FieSN7qtoJqLUE1CV2/8YlZdlR0ac4QytYTDhNZCTR5KxG0t3v5Z1Tbbxd2cbbJ9s4WNuJL2A9r8pykjivJJ0NJWmsL0mjLCc56l5o5gOqitevuF1j902btQGDONrQzV0/3U1tex/f+NBK7riodF6MSIs2RIQlOcksyUnm85ctpsMzyGvHW3j5SBOvHmvm93utrtBVBSlsLcvikiXZbCxNn09NuJuBClWtBBCRJ4AbgTGdVTSTHBfDFctzuWJ5LmCN2t1b3cHuU23sOd3Os4cbeHJ3NWBti7O6wBrcs6YwldWFKSzMSjIObBZRVRq6+jlY28WB2k4O1nayv6aTj5xXyNeuWzHl9Oe8s3r+cCN/9cR7JMS6ePIvtnD+goxwi2QIkbQENzesK+CGdQUEAsrBuk5eO9bMa8dbePSNk/zHq5W4nQ42lKSxZVEmFyzKYEPxnO7fKASqg/7XYPUBn8XI/azmCvFuJxcuzhweCKWqnGzpZW91B/trOtlX08HPdlYx4LP2jIuLcbAsL4WV+cksz0thaW4yS3OTyIyCmnmk4xn0UdHUw9GGbo40dHOkoYvDdV20e7wAOAQWZydx2dJszl+QPkFqoTFnmwFVlR++Xsk///EIawtT+Y9PbTSz6ecQvQM+dp1sY8eJFnZWtnGornO4f2N1YerwSvPrS9IoSI2b9Zr0DDUDfgy4esScxs2q+qWx4kRzM+C54PMHONHcy4HaTg7XdXG4vpPy+m46+7zDYTIT3SzJSWJxThKL7E1KS7MSKUqPN1MpglBVGrsGONnSy8mWXk4093CiuYeKph5qO/oYch2xLgfL8pJZmZ/CyoIUVhWksCI/hQR3aHWhWW0GnOLagJPqMA4Ff0D5X08d4qc7q/jgmny+/fF186mpaF6QGOvi8uXWAAyw+iT3VLWx62Q7e6ra+OnOKv7zjZPAmc0vh4xoeV4yCzITo/HBFMo0kXmNy2k9OJflJcP51jFVpal7gCMN3Rxv7OZ4Yw/Hm7r5w/76s5yY0yEUpsVTkpFAcUY8RekJFKTFUZiWQH5qHLkpceP2vUQb/oDS2jNAfWc/dR191Hb0UdPeR027h9Nt1qffe2Zn81iXg4VZiWwoSefjG4spy0limW1Ls9HcOmVnZXf6/jtBnb4i8tSITt9rgTL7cwHwfeCCEONOin6vn3ueeI9nDzXyF5cu4t5rlpvNBOcB1pytM/0bg74A5fVd7Kvp4EBNJwdqO3n9eAt+u4M+ximUZiayKDuR0sxEijMSKEqPpyAtntzkOFLiXZHYr/kOUCYiC7GmiXwCuC28IkU+IkJuiuVsLluaPXxcVWnrHRyuOVS1eqiyH9LPH26kpWfwfWllJrrJSYkjJzmWrKRYspLdZCXGkp7oJiMxhtR4N6nxMaTGx5Ac55rVl2SfP0B3v4+ufi+dfV7aPV46PIO09w7S1jtIc88gzd0DNPcM0NzVT1P3wPCAlSES3U6KMxJYkJnIpWXZLMhMoDQrkYVZiRSkxof1WTodNatQOn1vBH6iVpvjThFJE5F8oDSEuCHT3e/lrp/s4a3KVr5+/Uo+u3XhOStliG7cLgfritNYV5w2fKzf66eiqYdjjd0cazzTpPHykWYG/YGz4zsdZCS6SUuIISUuhoRYJ/ExTtwuB06HEONw8C83r51VnVTVJyJfBJ7lzDSRQxNEM4yBiJCZFEtmUiwbS9/fl9036Ke2o4+6jj7qO/uo7+ynsWuApq5+mnsGONbYTWvP4PvKTjBup4PEWCcJbhcJbifxbidxLiexMQ7cTgcxTgcup+B0CE4RRIShd6SAKoGA4lfLEXn9yqA/wKDPT783QL/Xj2fQj2fQR++Anz6vf0w5HGKNqs1KiiU7OZbF2ZnkpcSRlxpHfmo8BWlxFKTGk5YQE4kvacD0OKtQOn1HC1MYYlwgtE7jobekB25Zz4c3FE5CBcN8IC7GyerCVFYXpp513B9QmrsHqGn3UNfZP/wwausZpKPPS3e/l7beQTyDfgZ9AfwBJVz2rKrPAM+E5+rzi3i3kyU5SeMuwKuqdA/4aOsZpN0zSGefd/jT3e+ju99H74CP3kEfngE//T4//V4/PQM+Bn0BvP4APr/iV8Uf0OF+IFVFxHZiDsHlEGKcDmKcQqzLSUp8DDnJscS7LUeYFOskOS6GpFgXKfExpMXHkJZgfTISY0mNj4n6kZHhXhswlLjWwaC1AUWkWUTGXMn2pr8jCwh5FfcIx+gSmWTJfePqsmDWJBmHPXv2tIxjK3MqP5g7usDc0mciXUKylXCvDegOIe77UNXs8c6LyO5IWOpmOjC6RCbRost4thItOoTCXNIF5pY+06XLdAxtGe70FRE3VqfvUyPCPAXcLhZbgE5VrQ8xrsFgMBjmOVOuWY3V6Ru8NiBWG/t1QAX22oDjxZ2qTAaDwWCYW0zLPKvROn1tJzX0W4EvhBp3Ghh1A8coxegSmcwFXeaCDkPMJV1gbukzLbpE5QoWBoPBYJhfzJ3p2AaDwWCYsxhnZTAYDIaIJ2qdlYhcIyJHRaRCRO4b5byIyIP2+f0icl445AyFEHT5pK3DfhHZISLrwiFnqEykT1C4TSLiF5GbZ1O+yRCKLiKyTUT2isghEXl1tmWcCGMrkYmxk0miqlH3wRo5eAJYhDVXax+wckSY64A/Yk083gK8HW65p6DLRUC6/fvaSNUlVH2Cwr2ENbjm5nDLPYW8ScNaHqzE/p8TbrnPQQdjKxGoS1A4YyeqUVuzGl6PUFUHgaE1BYMZXo9QVXcCQ+sRRhoT6qKqO1S13f67E2vydKQSSt4AfAn4NdA0m8JNklB0uQ34jaqeBlDVSNPH2EpkYuxkkkSrsxprrcHJhokEJivn57DegiOVCfURkULgJuAHRDah5M1SIF1EXhGRPSJy+6xJFxrGViITYyeTJFp3Cp7KeoSRRshyisjlWAa4dUYlmhqh6PMAcK+q+iN1hWebUHRxYe2cdCUQD7wlIjtV9dhMCxcixlYiE2Mnk7STaHVWU1mPMNIISU4RWQs8Alyrqq2zJNu5EIo+G4EnbAPMAq4TEZ+q/m5WJAydUMtZi6r2Ar0i8hqwDogUZ2VsJTIxdjJZOwl359w5dui5gEpgIWc69FaNCPNBzu403hVuuaegSwnWUlUXhVve6dBnRPjHiNyO41DyZgXwoh02ATgIrA637JPUwdhKBOoyIvy8t5OorFnpFNYjjDRC1OXrQCbwPfsty6cRuiJziPpEBaHooqrlIvInYD8QAB5R1YPhk/psjK1Epq0YO5m8nZjllgwGg8EQ8UTraECDwWAwzCOMszIYDAZDxGOclcFgMBgiHuOsDAaDwRDxGGdlMBgMhojHOCuDwWAwRDzGWRkMBoMh4jHOymAwGAwRj3FWBoPBYIh4jLMyGAwGQ8RjnJXBYDAYIh7jrAwGg8EQ8RhnZTAYDIaIxzgrg8FgMEQ8xllNIyJySkS2R4Acy0TkPRHpFpG7Z/nah0Rk22xe0xA6pozOXhkVkcdE5FtTiD8jckarjRpnNUlEZKuI7BCRThFpE5E3RWRTuOUawVeBV1Q1WVUfnKmLjPbgU9VVqvrKDF2rUUQSg47dKSLTfq1ox5TRM8xmGZ1upkPOaLFRO27peGGMs5oEIpICPA18F8gACoF/BAbCKdcoLAAOhVuIGcAF3BNuISIZU0ajHxGJyh3cbWbORifa9958znyAjUDHOOdPAV/B2rq5E3gSiAs6fx9wAugGDgM3jYj7t/bxduBHI+IWAL8GmoGTwN1jyPAS4Af6gR5gKaDAkqAwjwHfCkVuoBj4jX3dVuAh+/hPsban7rOv89WgtLbbv1cArwAdWA+mG0K9V2Pc2/uANiDNPnYn1tt52MtGpHxMGQ1rGd0AvGvfuyeBJ4Z0mOje2Ne6177WANZD/xSw3c6TX40I/x3gwfHybCL9J0p3Mnk6VRu145aOGybcxhVNHyDFNoYfA9cC6aPc8F12BmcA5cDng85/zD7nAG4BeoH8oLgHbcPLAN4MKugOYA/wdcANLAIqgavHkPMV4M6g/6E8CN4nN+AE9gH/BiQCccDWEfG2j3IPtgMxQAXwNVvmK2xjWhbKvRqjMG/HeigN3RfjrEwZjYgyasevAv7aTvdmwAt8K5R7Y19rr31v40fIuQDwACn2cSdQD2wJMc/G0n+idCebp0PpzoiNmmbASaCqXcBWLMP6IdAsIk+JSG5QsAdVtU5V24D/BtYHxf+lfS6gqk8Cx4HNQXEfUtVqO+4/AbfaxzcB2ar6TVUdVNVK+/qfmEb1RpN7M5YR/I2q9qpqv6q+EWJ6W4Ak4H5b5pewmqduDQoz5r0ah68DXxKR7OCDIpJr99O8KiIviUh+iHLOKUwZDVsZ3YLlpB5QVa+q/gp4xz4X6r150L63fcEHVbUKq8b2YfvQFYBHVXfa5yfKs1GZKN1JyD2SGbFR46wmiaqWq+odqloErMYylAeCgjQE/fZgGQMAInK7iOwVkQ4R6bDjZwWFrw76XWWnDdYbUMFQPDvu14DgB9BUGU3uYqBKVX3nkF4BUK2qgaBjVVh9KONdc1xU9SDWA+W+EadasN6oLwN+AnzuHGSeE5gyGjLTWUYLgFq1qxNBaUHo96aasXmcM070Nvs/EFKejceY6U5C7rOYKRuN5o68sKOqR0TkMeAvJgorIguw3kquBN5SVb+I7AUkKFhx0O8SoM7+XQ2cVNWycxTVAyQE/c8DakKIVw2UiIhrjIeBjnJsiDqgWEQcQQ+DEuBYKAJPwDew3gi/PSyIqj/ofDKm8x4wZZTZK6P1QKGISJDDKsHqSwr13own6y+Bb4tIEXATcCGElGfjpTlmujZTydNpt1FTs5oEIrJcRL5sZywiUoz1VrJz/JiA1Z6uWB2ViMhnsN6AgvmCiBSJSAbWG8yT9vFdQJeI3Csi8SLiFJHVkxiOvBe4zY53DXBZiPF2YRnh/SKSKCJxInJx0PlGrHbs0Xgbq+38qyISY8/r+BBWp/OUUNUKrHtz1vwcEVkvIm8DX8QylHmHKaNhK6NvAT7gbhFxichHONMUN9V7g6o2Y/Xz/QjLgZTbpybKs/H0Hy/dKck9EzZqnNXk6AYuAN4WkV6sB8BB4MsTRVTVw1hvGW9hFaA1WB3UwTwOPIfViVmJ1Tk79EbyIaz28pNY1elHgNQQ5b7Hjt8BfBL4XSiRgq67BDiN9aZ7S1CQfwb+3m4i+MqIuIPADVid/C3A94DbVfVIiDJPxDexDDX4mntV9QLgH7BGrc1HTBkNQxm10/oIcAfWSMlbsAYaTMe9GeJxrAEMw011IeTZmPqPl+40yT2tNipnN7EawoWInMIaHfVCuGWJRkQkVlUH7N9XY41Y+p9hFmtOYcqoYSpM1UZNn5VhrnCeiPwLZ+bvfDbM8hgMhrOZko0aZ2WYE6jqW8Cl4ZbDYDCMzlRtNCL6rETkURFpEpGD4ZYlXKhqqWlemZ+ISLGIvCwi5WItMnrPKGG2ibXW31778/XZltOUUUM4iYg+KxG5FGs5kJ+o6sjRRwbDnMaeHJmvqu+KSDLWqgEftjvPh8JsA76iqteHR0qDIbxERDOgqr4mE6y4G0xWVpaWloYc3GB4H35V/H7FH1Drd0AJ2L8DCgG1/iugQ8uS2b+HWJCZMFby7Nmzp0VVs8cMEISq1mMNv0ZVu0WkHGti6uFxI4aAsZXRUQVfIIDPr/gCI8qBnfeBMfIdQMSayCQiOAQcIjgcglPA6RCcDgdOh+ByCDFOByKjSWGA0G0lIpxVKIjIXcBdACUlJezevTvMEhkijUBAae0dpL6zj/rOfho6+2ns6qexa4DmngFaugdo7R2grXcQr3/sFgUHEOsU4lxOYmOcxLocxDjFfgiJ9WAS4Zl7LhkzDRGpGvPkONgvbRuw5gCN5EIR2Yc1mfUrqjrqpEpjKzDg83O61cPJll5Ot3mobvNQ095HbUcfjV39tHu8o8YbyvsEt4u4GAdulwO304FDrLwH60XGH1C8fmXA52fAF8Az4GfQH0CxJluNnJ2cGh9DXkocBWlxFKUnUJQez4LMREqzEijNTCQuxjmTtyOiCdVWosZZqerDwMMAGzduDH/bpSEsdPd7qWr1UNXq4XSb9alpP/MgGvQFzgrvcgjZybFkJ8eSnxrHmsJUMpLcZCS4SU90kxofQ2p8DCnxLpJiXSTHxhDvduJ2zX53rogkYa1w/VdqrfEXzLvAAlXtEZHrsOYhjbqywHyylUBAOdXay+H6Lo7Ud3O0sZvjjd2cbvMQCNI8OdZFYXo8hWnxnL8gndyUOLKTY8lKiiUj0U1aglUOkuNcxLrOzXEM+gJ093vp7PPS0eelvXeQlp4BmroGaOoeoL6zn7qOPvZUtdPVf8adiUBxegJlOUkszUtmeV4yqwpSWJiVNOwgDVHkrAzzB39AOd3moaKphxPNPVQ293CypZeTLb209AyeFTYj0U1xejwr81P4wMpcClLjKEiLpyAtntyUODIT3TiiwOBFJAbLUf2Xqv5m5Plg56Wqz4jI90QkS1VbZlPOcNPY1c+7Ve28V93BvuoODtV10TNgPfidDmFhViIrC1K4YV0Bi7KTKM1KpDQzgdT4GGSG2+LcLgeZSbFkJsVOGLazz2vV/Fp7OWGX8+ONPbx6rBmf7WXjY5ysKkhhbVEaG0rSOG9BOoVp8TOqQyRjnJUhbARsp3SkoZtjjd0cb+rheGM3lS29Z9WQspLcLMpK4srluZRmJbIwK4HijARKMhJIjosJowbTg1hP0f8EylX1X8cIkwc0qqqKyGasFqvWWRRz1lFVqtv62FnZys7KVnadaqOm3VqQ3O10sKIghZs2FLKmMJWVBSksyUmKmua01PgY1hSlsqbo7MUgBn0BTjT3cKiui4O1nRyo7eS/3q7i0TdPApCfGsem0gwuWJTBlkWZLMpKnHEnHClEhLMSkZ8D24AsEakBvqGq/xleqQzTSWefl/L6Lsrt5pojDV0ca+yhz3tmbcui9HjKcpK4bGk2i3OSWJydxJLsJFITot8hTcDFwKeAA/YipGCtu1cCoKo/wNof6S9FxIe1md4nNBKG8k4znR4vb1S08PrxZl4/3kJth+WcMhPdbCrN4DMXL+S8kjRWFqScc3NdJON2OViRn8KK/BRuPr8IAK8/wJH6bt493c6uU23srGzlqX3W+sF5KXFcvCSLS5dmcUlZNhmJ7nCKP6NExND1ybJx40adj53G0YCq0tDVz6HaLg7VdXGorpPD9V3Db8RgNd0tz0tmeV4Ky/KSWJaXQllOEomxEfHuNC2IyB5V3RhuOSLdVlSV4009vFDeyEvlTbx7up2AQnKci4sWZ3LxkiwuWpzJ4uykeVODmAhV5VSrhx0nWthxopU3K1ro8HgRgbVFaVy5PIcrV+SwMj8lKu5ZqLYyd54OhllHVant6BturjhYazVdtPZa/UoisDAzkfXFady6uYSVBSmszE8hJzk2KozIMDOoKnurO/jToQaePdjAqVYPAKsLU/jC5UvYtiybdUVpuJwRsWZBxCFi9c0tzErkkxcswB9QDtR28urRZl4+2sS/vXCMf33+GIVp8Vy9Ko9r1+Rxfkl6VPTdjodxVoaQGKox7a/p5EBNJ/trOzlY20mb7ZhcDqEsN5nLl+ewuiCF1YWprMhPmVO1JcO5o6ocquviqX11/GF/PbUdfbgcwkVLsrjzkkVsX5FLXmpcuMWMSpwOYX1xGuuL07hnexnN3QO8dKSR5w418jO7vysnOZbr1uRzw/oCNhSnReXLonmSGEalrXeQ/TUd7K/pZH9NB/tqOmnuHgAs41iam8z2FTmsKUxlTVEay/OSo6Zz2zB71LR7+N17tfz2vVpONPcS4xQuKcvmf35gKdtX5pIaP+f7I2ed7ORYbtlUwi2bSugZ8PFieSPPHKjn8V2neWzHKUoyEvjwhkI+sqGQ0qzEiROMEIyzMtA74ONgbSf7azrZW9PB/poOqtusPiYRWJSVyCVLslhbZDmmVQUpxjEZxqTf6+dPBxv4xe5qdpywBixesDCDOy9ZxLWr80hLmLuDACKNpFgXN64v5Mb1hXT1e3nuUCO/31vLd186zoMvHmfjgnQ+vrGYD67Nj/hWEDPAYp4x6AtwtKGbfTXWPJX9NZ0cb+oenkBZmBbPuuJU1halWc6pMHVODA+fbebjAIujDd08/nYVv3mvlu5+HyUZCXz0vCI+cl4hxRljL01lmH0aOvv57Xu1/HJPNZXNvSS6ndywvoBPXrCA1YWT3RNyapgBFgb8AaWyuYd9QU155fVdw3OYMhLdrC1K5erVeay3HVRWCBMaDYYhvP4Azx5q4Cc7qth1qg2308E1q/P4xOZitizMjPpO/blKXmocf7ltMZ+/bBHvnm7n57uq+e17tfx8VzXritO4fcsCrl+XH1HTA0zNao4QCChVbR4O1HZywHZMh2o76R205jElup2sLkxlXbFVY1pXlEZRenxUdrRGA3O9ZtXWO8jjb1fx051VNHYNUJKRwCcvKOFjG4vn9FyfuUxnn5ffvFvDz3ZWcaK5l8xEN5+8oIQ/27KAnJSZG/xialZzmKH10A4OzXKv6eRgXSfd9npjQxMLP3p+0XBz3uJss86YYepUNvfwyBsn+fWeGgZ8AS4py+L/u2kNly/LMbWoKCc1PobPXLyQOy4q5c2KVh7bcZLvvlzB9189wQ3rCvnzSxeyPC8lbPIZZxXhDPoCVDT1cKiuk0N1XRyu6+Jw/Zn10NxOByvyk/nQugLWFlrLtyzNTSbGzFExTCN7qzv4/isVPHe4kRing49sKORzWxdSlpscbtEM04yIsLUsi61lWZxq6eXRN0/yy901/PrdGrYty+YvL1vM5oUZs94qY5xVBNHSM8DRhm7K6y2HVF7fTUVT9/B2FvExTlbkJw+vh7aqMMU4JsOMoarsrGzjoZeP82ZFKylxLr6wbQmfvqiU7GTTtzkfKM1K5Js3ruavty/lZzureGzHKW55eCcbF6TzxSuWcNnS7FlzWsZZhYFOj5fjTdbCrUftRVyPNXaftaJ4TnIsy/NTuGxpNqsKrLXCFmYlmqY8w4yjquw40coDLxzjnVPtZCXF8rfXLueTWxaQFOHDmw0zQ3qimy9dWcadlyziyXdO8/Brldzxo3dYV5TKPdvLuHxZzow7LVPyZgh/QKlt76OypYfK5l4qW3o40dRLRXPP8ORasGpLS3OTuGJ5DsvyUliWm8zy/GQzKs8QFt6ubOXbzx9j18k28lLi+McbVnHLpmIzr84AQLzbyR0XL+S2Cxbwm3dreOjlCj772G7WFaXyP69axqVlWTPmtCLCWYnINcB3ACfwiKreH2aRQqKzz0ttex817R6q2/uobvNQ1dpLlb0zafButMlxLhZnWyuKl+UksSQniaW5yRSmxZuOacOENmBvI/Id4DrAA9yhqu9O1/X313Twf589yuvHW8hNieWbN1pOKpKGLhsiB7fLwSc2l/DR84v4zbs1PPhiBZ9+dBebSzP4m2uWsak0Y9qvGXZnJSJO4N+BDwA1wDsi8pSqHg6XTJ5BH6091i6fLT2DNHVbW6M3dfUPb5de19FH98DZm1cnxbooyUhgeV4yV63MY6G9ZfWi7CSyktxmmLhhVEK0gWuxdgYuAy4Avm9/T4mTLb3832eP8MyBBtITYvi761bwqQsXmJqUISRinA5u2VTCTRuKePKd03z3pQo+9oO3uGJ5Dl+9Ztm0jh4Mu7MCNgMVqloJICJPADcCk3ZWfYN+Xj3WhCoEFPyq+PwBvP4Ag74A/d4A/V4/Hq+f3gEfPQM+uvt9dPVZW1F39nlp9wzS7w28L20RyEy0tkYvyUxgy6IMe5vsBIrS4ynOSCA9YeZ3IzXMSUKxgRuBn9h7WO0UkTQRyVfV+nO5YGvPAN958TiPv30at8vB3VeW8eeXLDSrlRjOCbfLwacuLOXm84v50Y6TfP+VE1z3nde5+fwi/ucHlk3LIsWR4KwKgeqg/zWM8sYoIncBdwGUlJSMmlBr7wCf/9nELSMxTiHB7SIp1kVynIuUuBiK0hNYXRhDekIM6YluspJiyUqyvnOS48hMcptRd4aZIhQbGC1MIfA+ZxWKrXz3pQr+6+3T3Lq5mLuvLCMn2ax4bpg68W4n/2PbEm7dVMJDL1fwk7dO4fMr/3rL+imnHQnOarSqyPuW1VDVh4GHwZqVP1pCOclx/PGeSwBwiOB0gMvhIMblwO10EBvjIM7lxO0yTscQUYRiAyHZCYRmK1+8Ygl/tmUBS3KSJiOnwRAS6Ylu/uH6ldxxUem0jWCOBGdVAxQH/S8C6saLsGfPnhYRqRonSBbQMg2yRQJGl8hkIl0WTCKtUGxg0nYCE9rKfMqPaGMu6TMtthIJzuodoExEFgK1wCeA28aLoKrZ450Xkd2RsC7bdGB0iUymWZdQbOAp4It2f9YFQGco/VXj2YrJj8hlLukzXbqE3Vmpqk9Evgg8izVs91FVPRRmsQyGWWMsGxCRz9vnfwA8gzVsvQJr6PpnwiWvwRAOwu6sAFT1GSxjNBjmJaPZgO2khn4r8IXZlstgiBTm6kiDh8MtwDRidIlM5oIuc0GHIeaSLjC39JkWXaJyPyuDwWAwzC/mas3KYDAYDHOIqHVWInKNiBwVkQoRuW+U8yIiD9rn94vIeeGQMxRC0OWTtg77RWSHiKwLh5yhMpE+QeE2iYhfRG6eTfkmQyi6iMg2EdkrIodE5NXZlnEijK1EJsZOJomqRt0Ha8TUCWAR4Ab2AStHhLkO+CPWZMotwNvhlnsKulwEpNu/r41UXULVJyjcS1iDCm4Ot9xTyJs0rGWRSuz/OeGW+xx0MLYSgboEhTN2ohq1NavhtdRUdRAYWkstmOG11FR1J5AmIvmzLWgITKiLqu5Q1Xb7706sCaGRSih5A/Al4NdA02wKN0lC0eU24DeqehpAVSNNH2MrkYmxk0kSrc5qrHXSJhsmEpisnJ/DeguOVCbUR0QKgZuAHxDZhJI3S4F0EXlFRPaIyO2zJl1oGFuJTIydTJKImGd1DkzrWmphJmQ5ReRyLAPcOqMSTY1Q9HkAuFdV/RG+Sn0ouriA84ErgXjgLRHZqarHZlq4EDG2EpkYO5mknUSrs5qxtdTCQEhyisha4BHgWlVtnSXZzoVQ9NkIPGEbYBZwnYj4VPV3syJh6IRazlpUtRfoFZHXgHVApDgrYyuRibGTydpJuDvnzrFDzwVUAgs506G3akSYD3J2p/GucMs9BV1KsJbZuSjc8k6HPiPCP0bkdhyHkjcrgBftsAnAQWB1uGWfpA7GViJQlxHh572dRGXNSufQWmoh6vJ1IBP4nv2W5dMIXeQyRH2iglB0UdVyEfkTsB8IYG1JfzB8Up+NsZXItBVjJ5O3E7OChcFgMBginmgdDWgwGAyGeYRxVgaDwWCIeIyzMhgMBkPEY5yVwWAwGCIe46wMBoPBEPEYZ2UwGAyGiMc4K4PBYDBEPMZZGQwGgyHiMc7KYDAYDBGPcVYGg8FgiHiMszIYDAZDxGOclcFgMBgiHuOsDAaDwRDxGGdlMBgMhojHOKtpREROicj2CJBjmYi8JyLdInL3LF/7kIhsm81rGkLHlNHZK6Mi8piIfGsK8WdEzmi1UeOsJomIbBWRHSLSKSJtIvKmiGwKt1wj+Crwiqomq+qDM3WR0R58qrpKVV+ZoWs1ikhi0LE7RWTarxXtmDJ6htkso9PNdMgZLTZqxy0dL4xxVpNARFKAp4HvAhlAIfCPwEA45RqFBcChcAsxA7iAe8ItRCRjymj0IyJRuYO7zczZ6ET73pvPmQ+wEegY5/wp4CtYWzd3Ak8CcUHn7wNOAN3AYeCmEXH/1j7eDvxoRNwC4NdAM3ASuHsMGV4C/EA/0AMsBRRYEhTmMeBbocgNFAO/sa/bCjxkH/8p1vbUffZ1vhqU1nb79wrgFaAD68F0Q6j3aox7ex/QBqTZx+7EejsPe9mIlI8po2EtoxuAd+179yTwxJAOE90b+1r32tcawHronwK223nyqxHhvwM8OF6eTaT/ROlOJk+naqN23NJxw4TbuKLpA6TYxvBj4FogfZQbvsvO4AygHPh80PmP2eccwC1AL5AfFPegbXgZwJtBBd0B7AG+DriBRUAlcPUYcr4C3Bn0P5QHwfvkBpzAPuDfgEQgDtg6It72Ue7BdiAGqAC+Zst8hW1My0K5V2MU5u1YD6Wh+2KclSmjEVFG7fhVwF/b6d4MeIFvhXJv7Gvtte9t/Ag5FwAeIMU+7gTqgS0h5tlY+k+U7mTzdCjdGbFR0ww4CVS1C9iKZVg/BJpF5CkRyQ0K9qCq1qlqG/DfwPqg+L+0zwVU9UngOLA5KO5Dqlptx/0n4Fb7+CYgW1W/qaqDqlppX/8T06jeaHJvxjKCv1HVXlXtV9U3QkxvC5AE3G/L/BJW89StQWHGvFfj8HXgSyKSPdpJEblVRJpDlHHOYcpo2MroFiwn9YCqelX1V8A79rlQ782D9r3tCz6oqlVYNbYP24euADyqutM+P1GejcpE6U5C7pHMiI0aZzVJVLVcVe9Q1SJgNZahPBAUpCHotwfLGAAQkdtFZK+IdIhIhx0/Kyh8ddDvKjttsN6ACobi2XG/BgQ/gKbKaHIXA1Wq6juH9AqAalUNBB2rwupDGe+a46KqB7EeKPeNPCciDqw32uqR5+YTpoyGzHSW0QKgVu3qRFBaEPq9Ga/cPs4ZJ3qb/R8IKc/GY8x0JyH3WcyUjRpnNQVU9QhWc8XqicKKyAKst5IvApmqmobVpCJBwYqDfpcAdfbvauCkqqYFfZJV9boQRfUACUH/80KMVw2UjNPhq2McB0v2YrtwDlEC1IZ47fH4BvDnnP1QAcvYfoXVTm/AlFFmr4zWA4UiEnyvSoJkDOXejCfrL4FtIlIE3ITtVELIs/HSHDPdSco9GtNuo8ZZTQIRWS4iX7YzFhEpxnor2Tl+TMBqT1esjkpE5DO8/wHyBREpEpEMrDeYJ+3ju4AuEblXROJFxCkiqycxHHkvcJsd7xrgshDj7cIywvtFJFFE4kTk4qDzjVjt2KPxNlbb+VdFJMae1/EhrE7nKaGqFVj3Znh+jog4gY9z5p7NS0wZDVsZfQvwAXeLiEtEPsKZprip3htUtRmrn+9HWA6k3D41UZ6Np/946U5J7pmwUeOsJkc3cAHwtoj0Yj0ADgJfniiiqh4Gvo1VqBuBNVgd1ME8DjyH1YlZidU5i6r6sYxoPdaInBbgESA1RLnvseN3AJ8EfhdKpKDrLgFOAzVYHbhD/DPw93YTwVdGxB0EbsDq5G8Bvgfcbr/pTwffxDLUIf4M+MWIJp35iCmjYSijdlofAe7AGil5C9ZAg+m4N0M8jjWAYbj2E0Kejan/eOlOk9zTaqNydhOrIVyIyCms0VEvhFuWaERE/gVr6HAAuBD4sarO6soIcx1TRg1TYao2Gs2TzwyGYVT13qHfIrLbOCqDIbKYqo2aZkDDnENVN4ZbBoPBMDbnYqPT4qxE5BoROSoiFSIy2nBFEZEH7fP7ReS8UOPOF1S11DSvGMZDRB4VkSYRORiO65syaggnU3ZW9giPf8fqpFwJ3CoiK0cEuxYosz93Ad+fRFyDwWDxGHBNuIUwGMLBdPRZbQYq7NnNiMgTwI1Ya1QNcSPwE3vC3E4RSRORfKA0hLjvIysrS0tLS993PKBKp8cLYk0yEBEcAg4RnA7r43I4OGsmhMEwSbx+pWfAS3qCe8wwe/bsaVHVUWfwnyuq+ppMsDL1SMayFYNhNujq8+JyOkhwO8cME6qtTIezKuTs2cg1WENnJwpTGGJcAETkLqxaGSUlJezevft9YWraPWz9l5cnFDglzkVeahx5qfEUpcdTkpFAaWYCC7OSKM1KINY19o01zG9OtvTyZ4+8jb/Pyx+/so3s5NhRw4lI1agnZoFQbMVgmCmauvv5xTvV/HxXNYMdfVy/oZB/vWX9mOFDtZXpcFaj1VNGjocfK0woca2Dqg8DDwNs3Lhx1DB5KXHsuO8K/AHFF1AGfQH6vH56B3x09Xnp6PPS0j1Ac88AjV391HX0c6Cmg3aPdzgNp0NYkJnAstxkVuSnsKoghTVFqeQkx412ScM8Yl91B5997B0U+PmfbxnTUYWbUGzFYJhOAgHljYoWfr7rNM8fbsQXUC5eksk/XL+S7StypuUa0+Gsajh7CZYizizBMlEYdwhxQ8bldFCQFj/peN39XqpaPZxo7qGiqYdjjd0cru/ijwfPLAuWlxLHuuJUNpSkc/6CdNYUphIXY2pg84UXyxv54uPvkZnk5sef3czi7AmXMTQY5jx1HX38ak8Nv9hdTU17H+kJMXx260I+samYRdNsI9PhrN4BykRkIdaaWp/AWv8pmKeAL9p9UhcAnapaL9bKuxPFnXGS42JYXZjK6sKzJ2b3DPgor+9if00n+2s62FvdwbOHGgFwOx2sK05l88IMtizKZOOCDOLHaZc1RCeqyqNvnuKf/nCYVQWp/OcdG00t2zCv6ff6ef5wI7/cU8Prx5tRhYuXZPLVa5Zz9arcGetGmbKzUlWfiHwReBZrP5RHVfWQiHzePv8D4BngOqy9YzzAZ8aLO1WZpoukWBebSjPYVJoxfKy1Z4A9Ve3srmpn18k2fvBqJf/+8glinMKGknS2Lslia1kWawtTcTnNNLZoZsDn5+9/e5Bf7qnhmlV5/Ost60hwh28evYj8HNgGZIlIDfANVf3PsAlkmDeoKnuq2vn1u7U8vb+O7n4fBalxfOmKMj52fhHFGQkTJzJFonK5pY0bN2qkdBr3DPjYfaqNtypbebOihUN1Xahagzi2lmVx2dJsti3LITfFvI1HE3Udffzlf73LvuoO7r5iCX+1fSkOR+jDSEVkTyRMTo4kWzFEHxVN3fx+bx2/31vH6TYP8TFOrl2dx0fOK+KixZmTsomxCNVWzHJLUyQp1sW2ZTlsW2Z1Irb1DrLjRAuvHWvmtWMtPHPA6vdamZ/CFctzuHx5DuuL03BOQyYbZoaXjzbx5V/sY9AX4Ad/dh7XrM4Pt0gGw6xR3ebh6f31/Pe+Og7Xd+EQuHhJFvdcWcY1q/NIjA2P2zDOaprJSHRz/doCrl9bgKpytLGbl4808/LRJr7/6gkeermCzEQ325bl8IGVOVxSlh22zDeczYDPz7efO8bDr1WyPC+Zf//keWYghWFeUNPu4Y8HGnj6QD37qjsAWF+cxj9cv5IPrc0nJwJahsxTcgYREZbnpbA8L4W/3LaYTo+XV48382J5I88fbuDX79bgdjm4aHEmV63MY/uKnIgoFPORow3d/PWTezlc38UnLyjhH65faUZ7GuY0FU3dPHuokWcPNbC/phOANYWp3HvNcj64Jp+SzJnvh5oMU3JW9gZsT2KtRHEK+Liqto8S7hrgO1iDKB5R1fvt4/8LazfJZjvo11T1manIFMmkJsRww7oCblhXgNcf4J1TbbxwuIkXyhv52m8P8LXfwoaSNK5amcdVq3LNW/0s4PUHePi1Sh544RgpcTE8cvtGtq+czp3YDYbIwB9Q3j3dzgvljTx/uJHK5l7AqkHdd+1yrlsdeQ4qmCkNsBCR/wO0qer99iK06cHLwNthnMAx4ANY863eAW5V1cO2s+pR1f83mevOtU7joebC5w818tzhRg7UWm85i7MTuWpVHletzGVdUdq0dGYazvDu6Xb+7rcHKa/v4oNr8vnmjavITJqeib5mgIUhEhhqzXn5SBOvHG2i3ePF5RAuXJzJVStz2b4yl/zUyc9NnU5ma4DFjVhDaQF+jLU98r0jwoSyduC8Jri58EtXllHX0ccL5Vb1/IevVfL9V06QkxzL9pW5XLUylwsXZ5oloaZAS88A/+/ZozzxTjW5KbH8x6fO5+pVeeEWy2CYMv6Asr+mg9eOtfDa8WbeO91OQCE9IYZty3K4ckUOly7NJiUuJtyiTpqpOqtcVa0HsCf5jrauxkTr/31RRG4HdgNfHq0ZEd6/3tlcpiAtntsvLOX2C0vp9Hh5+WgTzx5q4Hfv1fL426dJinVx2dJsPrAyl8uX5ZCaEH0FLxz0e/386M1TfO/lCvq8fv78koXcs30pSWaAiyFKUVVOtvTy5olWdlS0sONEK519XkSs/qcvXL6EbcvmxgjkCa1URF4ARnvt/LsQrzHe+n/fB/63/f9/A98GPjtaIvN1vbPUhBg+vKGQD28opN/rZ8eJFp4/3MgL5U384UA9ToewqTSd7StyuXJFLguzEsMtcsQx6Avwi93VPPRSBQ1d/VyxPIe/++CKqOwTHKv/1zA/UFUqW3rZdbKNtytbeauylcauAQAKUuO4amUuW8uyuKQsm4zEsXcFiEYmdFaqun2scyLSKCL5dq0qH2gaJdiYaweqamNQWj8Eng5V8PlIXIyTK5bncsXyXP4poOyr6eCF8kZeONzEt/5Qzrf+UM7CrEQuX5bD5cuz2bwwY143F3oGffzinWr+47VK6jv7OX9BOv92y3ouXJwZbtHOiaD934b7f0XkKVU1TepzlAGfn8N1XdaqOafa2V3VRkvPIABZSbFsWZTBhYszuWhxFqWZCcgc3v9oqu0fTwGfBu63v38/Spgx1w4ccnR2uJuAsOyAGo04HNbyThtK0vmbq5dT3ebhpSNNvHSkiZ+9XcWjb54kPsbJRYszuXRpNpeUZbEwK3FOF+Yh6jr6+NnOKh7fdZoOj5eNC9K5/6NrubQsK9r1n7b+34O1nVS3eVhbnEZBaly035c5QSCgVLV5htch3VvdwaG6LgZ9AQBKMhK4tCybTQsz2Lwwg0XzxJ6HmKqzuh/4hYh8DjgNfAxARAqwmiium2D9v/8jIuuxmgFPAX8xRXnmLcUZCXz6olI+fVEpnkEfb51o5dVjzbx6rJkXj1gV3sK0eC5eksnFS7K4cFHmnJrTNegL8PLRJn7xTjUvH7X0/cDKXO66dBHnL8iYIHbUENL+b6H07/7m3VoeffMkAJmJbnsh5xRW5qeyIj+Z0sxEM/p0BvH6A1Q291Je38Whuk4O1nZxsK6T7n4fAHExDtYUpnLHRaWcV5LGeSXpc8pezwWzNuA8oKq1l9ePt/D68WbeOtFKl20Qi7IS2Wy/pW1ckEFxRnxUvakN+gK8fbKVZw7U88eDDXR4vOQkx/KxjUXcurmEovTwzRmZiaHrIvIx4GpVvdP+/ylgs6p+aaw4Y9lKv9fPkYZu9td0cKCmkwO1nVQ09eALWM+D+BgnS/OSWZabxNLcZJbkJLEkJ4mC1HjjxCaBP6BUt3msrYeaujne2MORhm5ONPUw6LdqTG6XY3jvvLWFqawtSmNpbtK8WQjbrA1oGGZBZiILMhP5sy0L8AeUw3Vd7DjRwtsn2/jDgXqeeMd6Wc9KimVDSRrri9NYW5TKmsJU0sbZuj0cVLd5hmuNrx1vprvfR6LbyfaVuXx4QyGXLMmay0Yeyt5xIREX42R9sZXXQ/R7/VQ09XC4vovy+i6ONXbz0pEmfrG7ZjhMfIyThVmJLMpOZGFWIiUZCSzItL5zkmPnpSPz+QPUd/ZT3e7hdKuHk629nGrppbK5l6pWz7BTAshPjWNpbjKXlmWxIj+FFfkpLM5OnMtldtowzmqe4XQIa4pSWVOUyl9cthh/QDnW2M3uqnbeq2rnveoOnj88PO6FwrR4luclsywvmbLcJJZkJ1OalUDyLMzT6OzzcqS+i4N1Xeyt7mDPqTbqOvsByEmO5drVeXxgZR6XlGXNl6WRQtk77pyJi3GOuq9ba88AFU09VNibk55s6eVAbSd/PNiAP3CmZcbtdFCYHk9hWjz5qXHk29+5KbHkJMeRkxxLRqI7qh7MXn+Alp4BmroGaOoeoKGrn8bOfuo6+qjt6KOus4/6jv7hGilY96E4I55F2UlcsTyHxdlJLMm1aqbROL8pUjDOap7jdMjwG96ntiwAoMMzyKG6Lg7UdnK4rosjDV28eqz5LIPMTHRTnJFAYXo8Balx5KbEkZVkPYxS42NIjnORFOsizu0kxuHA4QBBCKgy4Asw4PXTM+Cjq99HW+8Azd0DNHYNUNveR1VbLydbeoeH5II1LHdDSTp3laZz0ZIsynKSoqrJcjoI1/5vmUmxZCbFcsGis0dRev0Batr7qGrtpbq9j5o2DzXt1kP81WPNNPcMMLKXQQTSE9xkJFqf9IQY0uLdpCbEnFVuEmNdJLpdxLudxMc4iYtx4HbZH6cDl9OBUwQRcNjfqhBQJaCKP6D4/IrXH2DQH6DfG6Df62fA56d3wI9n0EfPgJ+efi/d/T66+r10eLy0e7x0eAZp8wzS2jNIZ5/3fffD6RDyUuLIT43jvJJ0itbFU5SeYNcyE8hPjY/6OU2RiHFWhveRluDm4iVZXLwka/jYoC9AVWsvFU09VLb0Ut3mobrdw6HaTl4sb6TfGxgnxdDJSoplQWYCW5dksyQnieX5yazKT5n3nctD2GtnRsT6mTFOBwuzEsec2+f1B2i2ayNNXQM09wzQ3NVPS+8gbT2WQzjZ0kuHp4OOPu/wqLdwEB/jJM12mOkJblbkp5CR4CYrKZasZDfZSbHkpcaRkxxHdnKscUZhwDgrQ0i4XQ7KcpMpy01+3zlVpavPR0vvAO29g7R7vHT3e/EM+un3+vH6rbddVcXpcBDjFOJinCTFukiJd5GWYD0MspNj50tz3rwgxumgIC2egrTQ1p7r9/rp7vfRM+Cjd8CHZ9CqAQ3VigZ9AQb8Aby+AL5AAH+A4XI1hIjgdAhOEVxOweV0EOt0EBvjINblJN7tJNFtfSfFukiOiyEp1oXbFT1Nk/OVqBwNKCLNQNU4QbKAllkSZ6YxukQmE+myQFWzZ0uYsZjAVuZTfkQbc0mfabGVqHRWEyEiuyNhxevpwOgSmcwFXeaCDkPMJV1gbukzXbqYuq/BYDAYIh7jrAwGg8EQ8cxVZ/VwuAWYRowukclc0GUu6DDEXNIF5pY+06LLnOyzMhgMBsPcYq7WrAwGg8EwhzDOymAwGAwRT9Q6KxG5RkSOikiFiNw3ynkRkQft8/tF5LxwyBkKIejySVuH/SKyQ0TWhUPOUJlIn6Bwm0TELyI3z6Z8kyEUXURkm4jsFZFDIvLqbMs4EcZWIhNjJ5NE7Rng0fTBWhftBLAIcAP7gJUjwlwH/BEQYAvwdrjlnoIuFwHp9u9rI1WXUPUJCvcS1tJBN4db7inkTRrW5ocl9v+ccMt9DjoYW4lAXYLCGTtRjdqa1fCOqao6CAztmBrMjcBP1GInkCYi+bMtaAhMqIuq7lDVdvvvTqytISKVUPIG4EvAr4Gm2RRukoSiy23Ab1T1NICqRpo+xlYiE2MnkyRandVoO6YWnkOYSGCycn4O6y04UplQHxEpBG4CfjCLcp0LoeTNUiBdRF4RkT0icvusSRcaxlYiE2MnkyRaF7IdbcnjkWPwQwkTCYQsp4hcjmWAW2dUoqkRij4PAPeqqj/Ct/kIRRcXcD5wJRAPvCUiO1X12EwLFyLGViITYyeTtJNodVah7Jg6bbuqzjAhySkia4FHgGtVtXWWZDsXQtFnI/CEbYBZwHUi4lPV382KhKETajlrUdVeoFdEXgPWAZHirIytRCbGTiZrJ+HunDvHDj0XUAks5EyH3qoRYT7I2Z3Gu8It9xR0KQEqgIvCLe906DMi/GNEbsdxKHmzAnjRDpsAHARWh1v2SepgbCUCdRkRft7bSVTWrHSMHVNF5PP2+R9gjZ65DqvgeoDPhEve8QhRl68DmcD37Lcsn0boiswh6hMVhKKLqpaLyJ+A/UAAeERVD4ZP6rMxthKZtmLsZPJ2YpZbMhgMBkPEE62jAQ0Gg8EwjzDOymAwGAwRj3FWBoPBYIh4jLMyGAwGQ8RjnJXBYDAYIh7jrAwGg8EQ8RhnZTAYDIaI5/8HCxMJl4REcksAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xplot = np.arange(0, h + h/100, h/100)\n",
    "\n",
    "fig, axs = plt.subplots(4, 2)\n",
    "axs[0, 0].plot(xplot, N_k[0](xplot))\n",
    "axs[0, 0].set_title(\"Shape function N$_1$\")\n",
    "axs[0, 1].plot(xplot, dN_k[0](xplot))\n",
    "axs[0, 1].set_title(\"Shape function derivative N$_1$'\")\n",
    "axs[1, 0].plot(xplot, N_k[1](xplot))\n",
    "axs[1, 0].set_title(\"Shape function N$_2$\")\n",
    "axs[1, 1].plot(xplot, dN_k[1](xplot))\n",
    "axs[1, 1].set_title(\"Shape function derivative N$_2$'\")\n",
    "axs[2, 0].plot(xplot, N_k[2](xplot))\n",
    "axs[2, 0].set_title(\"Shape function N$_3$\")\n",
    "axs[2, 1].plot(xplot, dN_k[2](xplot))\n",
    "axs[2, 1].set_title(\"Shape function derivative N$_3$'\")\n",
    "axs[3, 0].plot(xplot, N_k[3](xplot))\n",
    "axs[3, 0].set_title(\"Shape function N$_4$\")\n",
    "axs[3, 1].plot(xplot, dN_k[3](xplot))\n",
    "axs[3, 1].set_title(\"Shape function derivative N$_4$'\")\n",
    "\n",
    "# automatically fix subplot spacing\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: computation of the elemental matrices\n",
    "\n",
    "In the theory we have seen that the mass and stiffness matrices for the EB beam problem can be computed as:\n",
    "\n",
    "$$ M_{i,j}^k = \\int_{x_a}^{x_b} \\rho A N_i^k(x)N_j^k(x)dx, $$\n",
    "\n",
    "$$ K_{i,j}^k = \\int_{x_a}^{x_b} E I N_i''^k(x)N_j''^k(x)dx. $$\n",
    "\n",
    "Here we call the auxiliar function `elemental_matrices(P)` with `P` being a variable containting all the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.integrate as scp\n",
    "M_k = np.zeros((4, 4))\n",
    "K_k = np.zeros((4, 4))\n",
    "for idof in np.arange(0, 4):\n",
    "    for jdof in np.arange(0, 4):\n",
    "        def eqn_M(x):\n",
    "            return rho*A*N_k[idof](x)*N_k[jdof](x)\n",
    "        def eqn_K(x): \n",
    "            return E*I*ddN_k[idof](x)*ddN_k[jdof](x)\n",
    "        M_k[idof, jdof] = scp.quad(eqn_M, 0, h)[0]\n",
    "        K_k[idof, jdof] = scp.quad(eqn_K, 0, h)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----------------------------------------------------------------------------------------------------\n",
    "**Problem**: Implement directly the elemental matrices given in the theory and compare with the ones computed here.\n",
    "\n",
    "-----------------------------------------------------------------------------------------------------\n",
    "\n",
    "In this example we also have a distributed dynamic loading function. This will give us a contribution to the right-hand-side vector that reads:\n",
    "\n",
    "$$ Q_i^k = \\int_{x_a}^{x_b} q(x)N_i(x)dx $$\n",
    "\n",
    "Since we want to integrate a multi-variable function $q(x,t)$ over one variable only $x$, we will use a symbolic integral here. Note that in that case, $q(x)$ is not the same for all elements, so we'll have to store all the elemental contributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Implement directly here ...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def make_Q(eqn, xe):\n",
    "#     return lambda t: scp.quad(eqn, xe[0], xe[1], args=(1,))[0]\n",
    "\n",
    "def make_Q(eqn, xe):\n",
    "    return scp.quad(eqn, xe[0], xe[1])[0]\n",
    "\n",
    "Q2_k = []\n",
    "for ie in np.arange(0, ne):\n",
    "    dofs = elem_dofs[ie]\n",
    "    nodes = [dof_node[dofs[0]], dof_node[dofs[-1]]]\n",
    "    xe = xn[nodes]\n",
    "    Q_k = []\n",
    "    for idof in np.arange(0, 4):\n",
    "        def eqn1(x): \n",
    "            return (1e-3*(L-x))*N_k[idof](x)\n",
    "        def eqn2(x):\n",
    "            return 1e3*N_k[idof](x)\n",
    "        def make_Q1():\n",
    "            return lambda t: [make_Q(eqn1, xe), make_Q(eqn2, xe)]\n",
    "        #Q_k.append([make_Q(eqn1, xe), make_Q(eqn2, xe)])\n",
    "        Q_k.append(make_Q(eqn1, xe))\n",
    "    Q2_k.append(Q_k)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.0008460000000000003, 8.280000000000006e-05, 0.0007739999999999996, -7.919999999999996e-05], [0.001745999999999999, 0.00040679999999999964, -0.0004859999999999993, 0.0010007999999999992], [0.012365999999999993, 0.0031067999999999972, -0.011465999999999992, 0.005104799999999997], [0.026225999999999996, 0.006886799999999997, -0.025685999999999987, 0.009640799999999996], [0.019565999999999986, 0.005266799999999996, -0.01938599999999999, 0.006832799999999995]]\n"
     ]
    }
   ],
   "source": [
    "print(Q2_k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: global assembly\n",
    "\n",
    "To construct the global matrices we add all the elemental contributions to the corresponding nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = np.zeros((ndofs*ndofs))       # 1-D array for global stiffness matrix\n",
    "M = np.zeros((ndofs*ndofs))       # 1-D array for global mass matrix\n",
    "Q = np.zeros(ndofs)              # 1-D array for Q vector\n",
    "\n",
    "for ie in np.arange(0, ne):\n",
    "    # Get the nodes of the elements\n",
    "    dofs = elem_dofs[ie]\n",
    "    NodeLeft = dof_node[dofs[0]]\n",
    "    NodeRight = dof_node[dofs[-1]]\n",
    "    \n",
    "    # Get the degrees of freedom that correspond to each node\n",
    "    Dofs_Left = 2*(NodeLeft) + np.arange(0, 2)\n",
    "    Dofs_Right = 2*(NodeRight) + np.arange(0, 2)\n",
    "\n",
    "    # Assemble the matrices at the correct place\n",
    "    nodes = np.append(Dofs_Left, Dofs_Right)\n",
    "    for i in np.arange(0, 4):\n",
    "        Q[nodes[i]] = Q[nodes[i]] + Q2_k[ie][i]\n",
    "        for j in np.arange(0, 4):\n",
    "            ij = nodes[j] + nodes[i]*ndofs\n",
    "            M[ij] = M[ij] + M_k[i, j]\n",
    "            K[ij] = K[ij] + K_k[i, j]\n",
    "            \n",
    "# Reshape the global matrix from a 1-D array to a 2-D array\n",
    "M = M.reshape((ndofs, ndofs))\n",
    "K = K.reshape((ndofs, ndofs))\n",
    "\n",
    "# Transform the Q vector into functions dependent on time\n",
    "Q_vector = []\n",
    "for i in np.arange(0, ndofs):\n",
    "    Q_vector.append(lambda t: Q[i] + 0.0*t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the `spy` command from the MatPlotLib.PyPlot package and inspect the structure of the matrices. Since the nodes are only attached to two elements at each side, we will only have contributions from the neighbour nodes. Therefore, you see that the matrices have a tri-block-diagonal structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7faa60f3ba90>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD8CAYAAACxd9IeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAKpUlEQVR4nO3dX4ilB32H8efbXUU3Vow4iu6GJkJIK0KJO5RoQEpWIa3ietFChEgqwt74J4ogsTe59UJEL4qwxGjAEClrwCDFGqIihRI62QRMskokbpPV1Z0gVfEmBn+9mBO6mc5klvO+Z943+T0fCDNz5rDny2yefc85c+adVBWSXv7+bOoBkvaHsUtNGLvUhLFLTRi71ISxS03MJvYkNyb5aZKfJbltBnuuSPKDJGeSPJbk1qk3ASQ5kOThJN+ZwZbXJTmV5CeLr9M7J97z6cXf1aNJ7knyqgk23JnkQpJHL7rs9UnuT/LE4u3l+70LZhJ7kgPAvwB/B7wN+FCSt027iueAz1TVXwHXAR+bwSaAW4EzU49Y+DLw3ar6S+CvmXBXksPAJ4H1qno7cAC4aYIpXwdu3HbZbcADVXU18MDi4303i9iBvwF+VlVPVtWzwDeB41MOqqrzVXV68f7v2fof+fCUm5IcAd4H3DHljsWW1wLvBr4KUFXPVtX/TDoKDgKvTnIQOAT8cr8HVNWPgN9su/g4cNfi/buAD+7npufNJfbDwNMXfXyOicO6WJIrgWuBByee8iXgs8CfJt4B8FZgE/ja4mHFHUkum2pMVf0C+ALwFHAe+G1VfW+qPdu8qarOw9ZBBHjjFCPmEnt2uGwWr+NN8hrgW8Cnqup3E+54P3Chqh6aasM2B4F3AF+pqmuBPzDR3VOAxePg48BVwFuAy5LcPNWeOZpL7OeAKy76+AgT3AXbLskr2Ar97qq6d+I51wMfSHKWrYc5NyT5xoR7zgHnqur5ezun2Ip/Ku8Bfl5Vm1X1R+Be4F0T7rnYr5O8GWDx9sIUI+YS+38BVye5Kskr2Xpi5b4pByUJW49Hz1TVF6fcAlBVn6uqI1V1JVtfn+9X1WRHrqr6FfB0kmsWFx0DHp9qD1t3369Lcmjxd3eM+TyReR9wy+L9W4BvTzHi4BQ3ul1VPZfk48C/s/Us6p1V9djEs64HPgz8OMkji8v+uar+bbpJs/MJ4O7FP9BPAh+ZakhVPZjkFHCare+kPAyc3O8dSe4B/hZ4Q5JzwO3A54F/TfJRtv5R+sf93gUQf8RV6mEud+MlrZixS00Yu9SEsUtNGLvUxOxiT3Ji6g0Xc8/e5rbJPTubXezALL4wF3HP3ua2yT07mGPsklZgX19Uk8RX8Ozh6NGjL/r5zc1N1tbW9mnNpZnbps57zp49yzPPPLPTD5bN4+Wy+j8bGxtTT9BL2Pr6+q6f82681ISxS00Yu9TEoNjndkZYSbtbOvaZnhFW0i6GHNlnd0ZYSbsbEvuszwgr6YWGfJ/9ks4Iu3hd8CxeLih1NiT2SzojbFWdZHEuMF9BJ01nyN342Z0RVtLulj6yz/SMsJJ2Mei18YvTKntqZeklwFfQSU0Yu9SEsUtNGLvUhLFLTRi71ISxS00Yu9SEsUtNGLvUhLFLTRi71ISxS00Yu9SEsUtNGLvUhLFLTRi71ISxS00Yu9SEsUtNDDq7rMaX7PSLdqZT5e/1eLnwyC41YexSE8YuNWHsUhPGLjVh7FITS8ee5IokP0hyJsljSW4dc5ikcQ35PvtzwGeq6nSSPwceSnJ/VT0+0jZJI1r6yF5V56vq9OL93wNngMNjDZM0rlFeQZfkSuBa4MEdPncCODHG7Uha3uDYk7wG+Bbwqar63fbPV9VJ4OTiur72UprIoGfjk7yCrdDvrqp7x5kkaRWGPBsf4KvAmar64niTJK3CkCP79cCHgRuSPLL47+9H2iVpZEs/Zq+q/wDm9fOYknblK+ikJoxdasLYpSaMXWrC2KUmjF1qwtilJoxdasLYpSaMXWrC2KUmjF1qwtilJoxdasLYpSaMXWrC2KUmjF1qwtilJoxdasLYpSZG+fVPl+ro0aNsbGzs502+qK1T3+vFzO1rVOUvFVqWR3apCWOXmjB2qQljl5owdqmJwbEnOZDk4STfGWOQpNUY48h+K3BmhD9H0goNij3JEeB9wB3jzJG0KkOP7F8CPgv8abcrJDmRZCPJxubm5sCbk7SspWNP8n7gQlU99GLXq6qTVbVeVetra2vL3pykgYYc2a8HPpDkLPBN4IYk3xhllaTRLR17VX2uqo5U1ZXATcD3q+rm0ZZJGpXfZ5eaGOWn3qrqh8APx/izJK2GR3apCWOXmjB2qQljl5owdqkJY5eaMHapCWOXmjB2qQljl5owdqkJY5eaMHapCWOXmjB2qQljl5owdqkJY5eaMHapCWOXmjB2qYlRzi77UlVVU0/4f5JMPWHW5vj1meP/RzvxyC41YexSE8YuNWHsUhPGLjUxKPYkr0tyKslPkpxJ8s6xhkka19BvvX0Z+G5V/UOSVwKHRtgkaQWWjj3Ja4F3A/8EUFXPAs+OM0vS2IbcjX8rsAl8LcnDSe5IctlIuySNbEjsB4F3AF+pqmuBPwC3bb9SkhNJNpJsbG5uDrg5SUMMif0ccK6qHlx8fIqt+F+gqk5W1XpVra+trQ24OUlDLB17Vf0KeDrJNYuLjgGPj7JK0uiGPhv/CeDuxTPxTwIfGT5J0ioMir2qHgHWx5kiaZV8BZ3UhLFLTRi71ISxS00Yu9SEsUtNGLvUhLFLTRi71ISxS00Yu9SEsUtNGLvUhLFLTRi71ISxS00Yu9SEsUtNGLvUhLFLTRi71MTQU0lrZFU19YQXSDL1hNl7qXyNPLJLTRi71ISxS00Yu9SEsUtNDIo9yaeTPJbk0ST3JHnVWMMkjWvp2JMcBj4JrFfV24EDwE1jDZM0rqF34w8Cr05yEDgE/HL4JEmrsHTsVfUL4AvAU8B54LdV9b2xhkka15C78ZcDx4GrgLcAlyW5eYfrnUiykWRjc3Nz+aWSBhlyN/49wM+rarOq/gjcC7xr+5Wq6mRVrVfV+tra2oCbkzTEkNifAq5LcihbLw4+BpwZZ5aksQ15zP4gcAo4Dfx48WedHGmXpJEN+qm3qroduH2kLZJWyFfQSU0Yu9SEsUtNGLvUhLFLTRi71ISxS00Yu9SEsUtNGLvUhLFLTRi71ISxS00Yu9SEsUtNGLvUhLFLTRi71ISxS00Yu9SEsUtNGLvUxKBTSevlr6qmnvACW7+PRMvwyC41YexSE8YuNWHsUhN7xp7kziQXkjx60WWvT3J/kicWby9f7UxJQ13Kkf3rwI3bLrsNeKCqrgYeWHwsacb2jL2qfgT8ZtvFx4G7Fu/fBXxw3FmSxrbsY/Y3VdV5gMXbN443SdIqrPwJuiQnkmwk2djc3Fz1zUnaxbKx/zrJmwEWby/sdsWqOllV61W1vra2tuTNSRpq2djvA25ZvH8L8O1x5khalUv51ts9wH8C1yQ5l+SjwOeB9yZ5Anjv4mNJM7bnD8JU1Yd2+dSxkbdIWiFfQSc1YexSE8YuNWHsUhPGLjVh7FITxi41YexSE8YuNWHsUhPGLjVh7FITxi41YexSE8YuNWHsUhPGLjVh7FITxi41YexSE8YuNWHsUhOpqv27sWQT+O89rvYG4Jl9mHOp3LO3uW3qvOcvqmrHX720r7FfiiQbVbU+9Y7nuWdvc9vknp15N15qwtilJuYY+8mpB2zjnr3NbZN7djC7x+ySVmOOR3ZJK2DsUhPGLjVh7FITxi418b+hlXtMNni5QgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.spy(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To apply the boundary conditions, we will remove the rows associated to the fixed DOFs and add the contribution to the right-hand-side. First, we identify the free and fixed DOFs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "fixed_dofs = np.arange(0, 2)                    # fixed DOFs\n",
    "free_dofs = np.arange(0, ndofs)                 # free DOFs\n",
    "free_dofs = np.delete(free_dofs, fixed_dofs)    # remove the fixed DOFs from the free DOFs array\n",
    "\n",
    "# free & fixed array indices\n",
    "fx = free_dofs[:, np.newaxis]\n",
    "fy = free_dofs[np.newaxis, :]\n",
    "bx = fixed_dofs[:, np.newaxis]\n",
    "by = fixed_dofs[np.newaxis, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can re-order the matrices and vectors in blocks, such that it's easy to operate with the blocks corresponding with the fixed DOFs. We'll use the notation $_I$ to designate an interior DOF and $_B$ to designate a boundary node.\n",
    "\n",
    "$$ \\begin{bmatrix} \\bold{M}_{II} & \\bold{M}_{IB} \\\\ \\bold{M}_{BI} &\\bold{M}_{BB} \\end{bmatrix}, \\quad \\begin{bmatrix} \\bold{K}_{II} & \\bold{K}_{IB} \\\\ \\bold{K}_{BI} &\\bold{K}_{BB} \\end{bmatrix}, \\quad \\begin{bmatrix} \\bold{Q}_{I}  \\\\ \\bold{Q}_{B}  \\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mass\n",
    "Mii = M[fx, fy]\n",
    "Mib = M[fx, by]\n",
    "Mbi = M[bx, fy]\n",
    "Mbb = M[bx, by]\n",
    "\n",
    "# Stiffness\n",
    "Kii = K[fx, fy]\n",
    "Kib = K[fx, by]\n",
    "Kbi = K[bx, fy]\n",
    "Kbb = K[bx, by]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now operating with the different blocks, we can get a system of ODEs for the free DOFs, which is given by:\n",
    "\n",
    "$$ \\bold{M}_{II}\\bold{\\ddot{u}}_I + \\bold{K}_{II} \\bold{u}_I = \\bold{Q}_I - \\bold{M}_{IB}\\bold{\\ddot{u}}_B - \\bold{K}_{IB} \\bold{u}_B $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "ub = np.array([u0, du0])\n",
    "ub_dt2 = np.array([u0_dt2, du0_dt2])\n",
    "RHS = -np.dot(Mib, ub_dt2) - np.dot(Kib, ub)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5: solving the ODE\n",
    "\n",
    "At this point we have all the information to solve the ODE for each DOF. We can use what we learned in the previous modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Time [s]')"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAERCAYAAABrWly6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAi0klEQVR4nO3dd3yV5f3/8deHDGbY04QtS9mEAI46aitWW1zfKlWsolKstUP9tvZba1vtsLU/W0cd1IEoQnHvVat1IoS9BIKshBVGwgqZn98f59giBpKQc3Kf8X4+HnmcnHPf574/Rz1vr1z3dV+XuTsiIhL/GgVdgIiIRIYCXUQkQSjQRUQShAJdRCRBKNBFRBKEAl1EJEEEGuhm9oiZbTOzpRE63utmVmRmLx/y+vtmtjD8s8nMno/E+UREYknQLfSpwNgIHu8OYMKhL7r7ye4+1N2HAh8Dz0bwnCIiMSHQQHf394CdB79mZr3DLe154ZZ1/zoc721gz+G2m1kGcDrw/FGWLCISs1KDLqAaU4DJ7r7azEYB9xEK4Ug4D3jb3XdH6HgiIjEjpgLdzFoAJwBPmdnnLzcObzsfuLWatxW4+5m1PMV44KH61ikiEotiKtAJdQEVhfu6v8Ddn6Uefd9m1g7IIdRKFxFJOEFfFP2CcFfIWjP7HwALGRKhw/8P8LK7H4jQ8UREYkrQwxZnEBp10s/M8s3sSuAS4EozWwQsA8bV4XjvA08BXw0f7+CumIuBGZGrXkQktpimzxURSQwx1eUiIiJHL7CLou3bt/cePXoEdXoRkbg0b9687e7eobptgQV6jx49yM3NDer0IiJxyczWH26bulxERBKEAl1EJEEo0EVEEoQCXUQkQSjQRUQShAJdRCRBKNBFRBKEAl1EpAHd9c/VLNtUHJVjx9r0uSIiCevpefn85Z+rKKus5PhjWkX8+Gqhi4g0gOWbdvOL55Ywplc7fnJG36icQ4EuIhJlxSXlXDN9Hq2apnH3+GGkpkQnetXlIiISRVVVzo1PLaJgVwkzJ42mQ0bjqJ1LLXQRkSi6/99reGv5Vn7+jQFk92gb1XMp0EVEouTfqwr585sr+daQY5h4Yo+on0+BLiISBRt37ueHMxbQr1MGt18wCDOL+jkV6CIiEVZSVsn3Hp+Hu/PghBE0S2+Yy5W6KCoiEkHuzv89t4QVW3bzyHdH0r1d8wY7t1roIiIR9NhH63huQQE//mpfTuvfsUHPrUAXEYmQOWt38ttXVnDGgI5cd/qxDX5+BbqISARs3X2A70+fT9e2zbjzoqE0ahT9i6CHUh+6iEg9lVVUcc0T89hfVsH0q0bRsklaIHUo0EVE6um2l5czf0MR935nGP06ZwRWh7pcRETq4ancjTw+ez2TvtKLcwYfE2gtCnQRkaO0OL+IXzy/lBN6t+OnZ/YLuhwFuojI0di2+wCTps2jQ4vG3BPFGRTrQn3oIiJ1VFpRyeQn5lFcUs7T14yhXYvozaBYFwp0EZE6cHd++fxS5m8o4m/fGR6VlYeOVvB/I4iIxJHHPlrHrNx8rjv9WM4e3CXocr5AgS4iUksf5m3ntldWcMaATlFbRq4+FOgiIrWwYcd+rn1yPr3aN+cvFw0J5E7QmijQRURqsLe0gqun5eIOf78sm4yA7gStiQJdROQIqqqcG2YtZPW2Pdz7nWH0aN9w0+HWlQJdROQI7np7NW8s28ovzj6Ok/t0CLqcI1Kgi4gcxutLN3PX26u5YHhWg6wJWl8KdBGRaqzYvJvrZy1iaNfW/O68gQ2yJmh91RjoZtbVzN4xsxVmtszMflTNPmZmd5tZnpktNrPh0SlXRCT6tu05wFWP5ZLRJJUHJ4ygSVpK0CXVSm3uFK0AbnD3+WaWAcwzs7fcfflB+5wF9An/jALuDz+KiMSVA+WVTJo2j537ynhq8hg6tWwSdEm1VmML3d03u/v88O97gBVA5iG7jQOmechsoLWZxdYtVCIiNXB3fvr0YhZuLOIvFw1hYGbs3NZfG3XqQzezHsAw4JNDNmUCGw96ns+XQ19EJKbd/XYeLy7axE/H9mPswPhrk9Y60M2sBfAM8GN3333o5mre4tUcY5KZ5ZpZbmFhYd0qFRGJopcWbeIv/1zF+cMzueaU3kGXc1RqFehmlkYozKe7+7PV7JIPdD3oeRaw6dCd3H2Ku2e7e3aHDrE9nlNEkseCDbu48alFjOzRhj+cPyguRrRUpzajXAx4GFjh7nceZrcXgcvCo11GA8XuvjmCdYqIREVBUQlXT5tHx5aNeXBCNo1T42NES3VqM8rlRGACsMTMFoZf+z+gG4C7PwC8CnwDyAP2A1dEvFIRkQjbV1rBVY/lUlpeyYyrR9G2eXrQJdVLjYHu7h9QfR/5wfs4cG2kihIRibbKKudHMxeycstuHr0ihz6dMoIuqd60YpGIJKU/vv4p/1yxlVvHHc8pfRPjmp5u/ReRpDNzzgamvPcZl43pzmVjegRdTsQo0EUkqfx7VSG/eH4pX+nbgVvOOS7ociJKgS4iSWP5pt18/4l59O2UwX2XDCc1JbEiMLE+jYjIYWwuLmHi1Lm0bJrGo5ePpEXjxLuEmHifSETkEHsOlHPFo3PZW1rBU5PH0LlV/Ey4VRdqoYtIQiuvrOL70+eTt20v9186nAFdWgZdUtSohS4iCcvdufm5pby/ejt/umBwzC8hV19qoYtIwvrbO3n8I3cj151+LN8e2bXmN8Q5BbqIJKTnFxTw5zdXcd6wTK7/Wt+gy2kQCnQRSTizP9vBT59ezOhebfnjBYPjdvbEulKgi0hCydu2h0nTcunWrhkPXppNemryxFzyfFIRSXhbdx/gu4/MJT01hUcvH0mrZmlBl9SgFOgikhCKS8r57iNzKNpfxtQrRtK1bbOgS2pwGrYoInHvQHklk6blsqZwL49enhN3iztHigJdROJaZZVz/ayFfLJ2J3ddPJST+rQPuqTAqMtFROKWu/Obl5bx6pIt3Hz2AMYNzQy6pEAp0EUkbt337hqmfbyeSV/pxVUn9wq6nMAp0EUkLs3K3cgdb6zkvGGZ3DS2f9DlxAQFuojEnbdXbOXnzy7h5D7t+eMFg2nUKDluHKqJAl1E4sr8Dbu49sn5HNelJfdfOiKpbhyqif5JiEjcyNu2hyunzqVTyyY8kqCLVNSHAl1E4kL+rv1c+tAcUho1YtrEHDpkNA66pJijQBeRmFe4p5QJD89hf1kFj1+ZQ/d2zYMuKSbp7xURiWnFJeVc9sgcthQf4ImrchJ6xaH6UgtdRGJWSVklV06dS962PTwwYQQjurcNuqSYpha6iMSksooqrpk+j3kbdnHv+OGc0jexl4+LBLXQRSTmfD4/y7srC/n9eYM4e3CXoEuKCwp0EYkp7s4tLyzl5cWbuems/ozP6RZ0SXFDgS4iMeXPb65k+icbmHxKbyaf0jvocuKKAl1EYsaU99bwt3fWMD6nGz8b2y/ocuKOAl1EYsK0j9fx+1c/5ezBXfjtuQOTZmHnSFKgi0jgZs3dyC0vLOOMAZ3460VDSdFkW0elxkA3s0fMbJuZLT3M9lPNrNjMFoZ/bol8mSKSqF5YWMDPnl3MyX3a87dLhpGWonbm0arNOPSpwL3AtCPs8767nxORikQkaby+dDPXz1rEqJ5tmTIhm8apKUGXFNdq/F+hu78H7GyAWkQkibzz6Taum7GAIVmteOi7I2marjCvr0j9bTPGzBaZ2WtmdvzhdjKzSWaWa2a5hYWFETq1iMSbD/O2870n5tGvcwaPXpGjaXAjJBKBPh/o7u5DgHuA5w+3o7tPcfdsd8/u0EG38Yokozlrd3LVY7n0bNecxyeOolXTtKBLShj1DnR33+3ue8O/vwqkmVn7elcmIgln4cYiJk6dS5dWTXjiqlG0aZ4edEkJpd6BbmadLTxg1MxywsfcUd/jikhiWVpQzGUPf0Kb5mlMv3qUFqiIgho7rsxsBnAq0N7M8oFfAWkA7v4AcCFwjZlVACXAxe7uUatYROLO0oJiLnnoEzKapPHkVaPp0qpp0CUlpBoD3d3H17D9XkLDGkVEvmRpQTGXPvwJLRqnMuPq0XRt2yzokhKWRvCLSNQs2xQK82ZpKcy4ejTd2inMo0mBLiJRsWxTqJulWVoKMyeNUZg3AAW6iETc8k27FeYBUKCLSEQt37Sb7zw0O9TNMkndLA1JgS4iERNqmc+maTjMu7drHnRJSUWBLiIRsWJzKMybpKUwU2EeCAW6iNTb0oJixv9dYR40BbqI1Mv8DbsY//fZNE9P5R+TxijMA6QpzkTkqH3y2Q4mTp1Lh4zGTL96NJmtdQdokBToInJU3l9dyNXTcslq04zpV42iU8smQZeU9BToIlJnb6/YyjXT59OrfXOeuGoU7Vtooq1YoEAXkTp5felmrpuxgAFdWjJtYg6tm2kK3FihQBeRWnthYQHXz1rE0K6tefSKkbRsosUpYolGuYhIrcyau5Ef/2MhI3u0YdrEHIV5DFILXURqNPXDtfz6peWc3Kc9UyZka0HnGKVAF5HDcnfuens1f/3nar5+XCfuHj+MJmkK81ilQBeRalVVObe+vJypH63jwhFZ3H7+IFJT1EsbyxToIvIl5ZVV/OzpxTy7oICJJ/bk5rMH0KiRBV2W1ECBLiJfcKC8kh88uYB/rtjKDV/ryw9OP5bwOvAS4xToIvIfew6Uc/W0XGZ/tpNbxx3PZWN6BF2S1IECXUQA2LG3lMsfncvyzbv560VDOXdYZtAlSR0p0EWETUUlTHj4E/J3lTBlwgi+OqBT0CXJUVCgiyS51Vv38N1H5rDnQAXTJuYwqle7oEuSo6RAF0lic9ft5KrHcklLacSMSaMZmNkq6JKkHhToIknq9aVb+NHMBWS2bspjE3Po2laLOcc7BbpIEnr843Xc8uIyhnZtzcPfHUnb5poxMREo0EWSiLtzxxsrue/dNZwxoCP3jB+ueVkSiAJdJEmUV1Zx0zNLeGZ+PuNzunLbuIG6lT/BKNBFksC+0gqumT6f91YV8pMz+vLDr+ruz0SkQBdJcIV7Spk4NXTD0O3nD+LinG5BlyRRokAXSWCrtu7hikfnsmNfqW4YSgIKdJEE9f7qQr7/xHyapKfw1PdOYFCWxpgnOgW6SAKaMWcDNz+/lD4dW/Dw5SPJbN006JKkAdR4idvMHjGzbWa29DDbzczuNrM8M1tsZsMjX6aI1EZVlfOHV1fw82eXcNKx7Xlq8hiFeRKpzZilqcDYI2w/C+gT/pkE3F//skSkrkrKKrn2yfk8+N5nXDq6Gw9/N5sMLeScVGrscnH398ysxxF2GQdMc3cHZptZazPr4u6bI1WkiBzZtj0HuHraPBbnF3Hz2QO48qSeGpaYhCLRh54JbDzoeX74tS8FuplNItSKp1s3DZ0SiYTPR7Ls3FfGg5eO4OvHdw66JAlIJG4Tq64Z4NXt6O5T3D3b3bM7dOgQgVOLJLd/fbqV8+/7iPLKKmZ9b4zCPMlFooWeD3Q96HkWsCkCxxWRw3B3Hvj3Z/zpjU85/piWTJmQzTG6+Jn0IhHoLwI/MLOZwCigWP3nItFzoLySm55ZzPMLN3H24C78+cIhmmBLgFoEupnNAE4F2ptZPvArIA3A3R8AXgW+AeQB+4ErolWsSLLbUnyA7z2ey6L8Ym78el+uPU1zssh/1WaUy/gatjtwbcQqEpFqLdxYxKRpuewtreDBCSM4U/3lcgjdKSoSB55bkM/PnllCx4zGPHvlCfTv3DLokiQGKdBFYlhFZRV3vLGSB9/7jFE923L/pSO0upAclgJdJEbt2FvKD2cu4MO8HVwyqhu/+ubxpKdqQQo5PAW6SAxanF/E5MfnsX1fGX+6YDDfHtm15jdJ0lOgi8SYWXM3cvMLS+nQojFPTx7D4KzWQZckcUKBLhIjSisq+fWLy5kxZwMnHtuOe8YPV3+51IkCXSQGbCoq4Zrp81m0sYjJp/Tmxq/31QLOUmcKdJGAfbRmO9c9uYAD5ZXcf8lwzhrUJeiSJE4p0EUCUlXl3PduHne+tYoe7ZvzjwmjObZjRtBlSRxToIsEYMfeUn4yaxHvrSrkm0OO4Q/nD6JFY30dpX70X5BIA5u7bifXPbmAnfvL+N15A/lOTjfNxyIRoUAXaSBVVc6D733Gn99cSdc2TXn2mhMYmNkq6LIkgSjQRRrArn1lXD9rIe+sLOTsQV24/YJBWu9TIk6BLhJl89bv4ron57N9bxm3jTueS0d3VxeLRIUCXSRKKquc+97J469vryazdVOeueYEBmWpi0WiR4EuEgUFRSX8ZOZC5qzbybihx3DbuQNpqS4WiTIFukiEvbpkMzc9s5jKKucvFw3hvGFZQZckSUKBLhIh+8sq+M2Ly/lH7kaGdG3N3RcPpXu75kGXJUlEgS4SAUsLivnhjAWs3bGPa0/rzY/P6Eua5mKRBqZAF6mHyirnofdDY8vbNW/Mk1eNZkzvdkGXJUlKgS5ylDbs2M8NTy1k7rpdjD2+M384fxBtNN2tBEiBLlJH7s6MORv57SvLSTHjzm8P4bxhmRpbLoFToIvUwdbdB/jZM4t5d2UhJx7bjjsuHMIxrZsGXZYIoEAXqbWXFm3i5ueXUlpRyW++dTwTRnenUSO1yiV2KNBFalC0v4xfvrCMlxZtYmjX1tz57SH06tAi6LJEvkSBLnIEry3ZzC9fWEbR/jJu/HpfJp/SW0vDScxSoItUY9ueA/zqhWW8tnQLAzNbMm1iDscd0zLoskSOSIEuchB355n5Bdz28nJKyiv52dj+XH1yT7XKJS4o0EXCCopK+L9nl/DvVYVkd2/DHy8cTG/1lUscUaBL0quqcqZ/sp7bX/sUB41gkbilQJektnLLHn7x3BJy1+/i5D7t+f15g+jatlnQZYkcFQW6JKX9ZRXc9fZqHn5/LRlNUrnjwsFcOCJLd3tKXFOgS9J5e8VWbnlhGQVFJXw7O4ufnzVAc7BIQqhVoJvZWOAuIAV4yN1vP2T7qcALwNrwS8+6+62RK1Ok/jYXl/DrF5fxxrKt9OnYglnfG0NOz7ZBlyUSMTUGupmlAH8DvgbkA3PN7EV3X37Iru+7+zlRqFGkXioqq5j60Tr+8tYqKt356dh+XHVSL9JTNRRREkttWug5QJ67fwZgZjOBccChgS4Scz7K285vXlrOyq17OK1fB24dN1AXPSVh1SbQM4GNBz3PB0ZVs98YM1sEbAJudPdlh+5gZpOASQDdunWre7UitbRx535+/+oKXlu6haw2TXng0hGceXwnXfSUhFabQK/uG+CHPJ8PdHf3vWb2DeB5oM+X3uQ+BZgCkJ2dfegxROqtpKyS+/+9hgf/vYZGZtzwtb5c/ZVeNElLCbo0kairTaDnA10Pep5FqBX+H+6++6DfXzWz+8ysvbtvj0yZIkfm7ryyZDO/f2UFm4oP8M0hx/Dzs/prrnJJKrUJ9LlAHzPrCRQAFwPfOXgHM+sMbHV3N7McoBGwI9LFilRnaUExt728nE/W7mRAl5b85aKhjOqldT0l+dQY6O5eYWY/AN4gNGzxEXdfZmaTw9sfAC4ErjGzCqAEuNjd1aUiUVVQVMKf31jJcwsKaNMsjd+eO5DxOd1I0S37kqQsqNzNzs723NzcQM4t8a24pJz73s3j0Q/XATDxxJ5cc2pvWjVNC7YwkQZgZvPcPbu6bbpTVOJGWUUVT8xezz3/Wk1RSTnnDc3khjP7kal+chFAgS5xwN15dckW/vTGp6zfsZ8Tj23Hz88awMDMVkGXJhJTFOgSs9ydd1cVcuebq1hSUEy/ThlMvWIkp/TtoPHkItVQoEtM+njNDv7fmyvJXb+LrDZNuePCwZw/PEsXPEWOQIEuMWX+hl3c+eYqPsjbTqeWjfntuQP5dnZXzbsiUgsKdIkJyzYVc+ebq3j70220a57OL885jktGddMdniJ1oECXQC3aWMS97+Tx1vKttGySyv+e2Y/LT+hB88b6T1OkrvStkUDMXbeTe/6Vx3urCmnVNI0fn9GHK07sqbHkIvWgQJcG4+58mLeDe/61mk/W7qR9i3RuOqs/l47uTgu1yEXqTd8iibqqKuefK7Zy37trWLixiE4tG3PLOccxPqcbTdPVRy4SKQp0iZoD5ZU8Mz+fh95fy9rt+8hq05TfnTeQC0dk0ThVQS4SaQp0ibgde0t5fPZ6pn28np37yhic1Yp7vzOMscd3JjVFww9FokWBLhHzWeFeHv5gLU/Py6e0ooozBnTk6pN7kdOzre7sFGkACnSpl8oq592V23js4/W8t6qQ9NRGXDA8kytP6smxHTOCLk8kqSjQ5agU7S9jVu5GHp+9no07S+jUsjHXf60v43O60SGjcdDliSQlBbrUydKCYh7/eD3PLyygtKKKnJ5tuWnsAL5+fCfS1D8uEigFutRoz4FyXly0iX/M3cji/GKapqVwwYgsLhvTnf6dWwZdnoiEKdClWu7OvPW7mDl3I68s3kxJeSX9O2fwq28ex/nDs3RHp0gMUqDLFxTuKeX5BQXMnLuBNYX7aJ6ewrnDMrl4ZFcGZ7XSaBWRGKZAF/aXVfDmsq08t6CAD/K2U1nljOjehj9d2JuzB3XRRFkicULf1CRVUVnFh2t28PyCAt5YtoX9ZZVktm7K5FN6cd6wTA05FIlDCvQkUlkV6hd/dclmXl68me17S2nZJJVxQzM5b1gm2d3b0EgrAonELQV6gquorGLOup28tmQLry/bQuGeUtJTG3F6v46cOyyT0/p30LwqIglCgZ6AyiqqmP3ZDl5buoU3l21hx74ymqQ14vT+HTlrYBdO699R09WKJCB9qxNE4Z5S3lm5jXc+3cb7q7ezt7SC5ukpnD6gE98Y2JlT+nWgWbr+dYskMn3D45S7s2zTbv716Tbe/nQbizYWAdC5ZRO+NfQYTu/XkZP6tNeanCJJRIEeRwqKSvhw9XY+yNvOR2u2s31vGWYwJKs1N3ytL6cP6MhxXVpqrLhIklKgx7Ci/WXM/mwHH+Rt58O8Hazdvg+A9i0ac9Kx7Tnx2Pac2q+jJsMSEUCBHjPcnfxdJeSu38ncdbvIXbeTVVv3AtAsPYXRvdpx6ejunHRse/p2aqFWuIh8iQI9IAfKK1m+eTeLNxaRu34Xuet2sWX3AQAyGqcyvHsbvjXkGHJ6tmNo19akp2omQxE5MgV6AyitqGTllj0szi9mSX4xiwuKWbV1D5VVDkCXVk3I6dmWkT3aMKJ7W/p1ziBFN/iISB0p0COoqirUbbJy6x5Wbd3Dyi2hxzWFeymvDIV3m2ZpDMpqzVf7d2RQVisGZ7WiS6umAVcuIolAgX4UiveXs27HPtbt2Mf6HftZt2Mfa7btZdXWvZSUV/5nv8zWTenXOYNT+3VkcFYrBmW2IqtNU/V/i0hU1CrQzWwscBeQAjzk7rcfst3C278B7Acud/f5Ea61QVRWOdv3lrK5+ABbikvCjwfYVHyADTv3s37HPor2l3/hPV1aNaFn++ZcnNOVfp0y6Ns5gz4dW5DRRHOGi0jDqTHQzSwF+BvwNSAfmGtmL7r78oN2OwvoE/4ZBdwffgyEu1NaUcW+0gr2l1Wyr6yCfaWVlJRVsvtAOTv3lbFrXxk793/+WM6ufWVs31vKtj2l/+nb/lx6aiM6t2xCt7bNOHtQF3q0a073ds3o0b453do20807IhITatNCzwHy3P0zADObCYwDDg70ccA0d3dgtpm1NrMu7r450gW/u3Ibv31lBZVVTkVVFZWVTkWVU1nllFdWUVnllJRXckgmVyujSSptm6fTplk6HTIa07dTBl1aNaFzqyYHPTalTbM0dZOISMyrTaBnAhsPep7Pl1vf1e2TCXwh0M1sEjAJoFu3bnWtFYCMJmn06xQaBZLayEKPKeHHRo1IaWQ0TUuhWeMUmqen0iw9heaN//vYonEq7Zqn07pZuoYCikhCqU2gV9c0PbT9W5t9cPcpwBSA7OzsWrShv2xE9zaM6N7maN4qIpLQatNEzQe6HvQ8C9h0FPuIiEgU1SbQ5wJ9zKynmaUDFwMvHrLPi8BlFjIaKI5G/7mIiBxejV0u7l5hZj8A3iA0bPERd19mZpPD2x8AXiU0ZDGP0LDFK6JXsoiIVKdW49Dd/VVCoX3waw8c9LsD10a2NBERqQsN8xARSRAKdBGRBKFAFxFJEAp0EZEEYaHrmQGc2KwQWH+Ub28PbI9gOfFAnzk56DMnh/p85u7u3qG6DYEFen2YWa67ZwddR0PSZ04O+szJIVqfWV0uIiIJQoEuIpIg4jXQpwRdQAD0mZODPnNyiMpnjss+dBER+bJ4baGLiMghFOgiIgki7gLdzMaa2UozyzOzm4KuJ9rM7BEz22ZmS4OupaGYWVcze8fMVpjZMjP7UdA1RZuZNTGzOWa2KPyZfxN0TQ3FzFLMbIGZvRx0LQ3BzNaZ2RIzW2hmuRE9djz1oYcXrF7FQQtWA+MPWbA6oZjZV4C9hNZsHRh0PQ3BzLoAXdx9vpllAPOAcxP837MBzd19r5mlAR8AP3L32QGXFnVmdj2QDbR093OCrifazGwdkO3uEb+ZKt5a6P9ZsNrdy4DPF6xOWO7+HrAz6Doakrtvdvf54d/3ACsIrVGbsDxkb/hpWvgnflpbR8nMsoCzgYeCriURxFugH24xaklQZtYDGAZ8EnApURfuelgIbAPecveE/8zAX4GfAlUB19GQHHjTzOaZ2aRIHjjeAr1Wi1FLYjCzFsAzwI/dfXfQ9USbu1e6+1BCa/LmmFlCd7GZ2TnANnefF3QtDexEdx8OnAVcG+5WjYh4C3QtRp0kwv3IzwDT3f3ZoOtpSO5eBLwLjA22kqg7EfhWuE95JnC6mT0RbEnR5+6bwo/bgOcIdSVHRLwFem0WrJY4F75A+DCwwt3vDLqehmBmHcysdfj3psAZwKeBFhVl7v5zd89y9x6Evsv/cvdLAy4rqsysefhCP2bWHPg6ELERbHEV6O5eAXy+YPUKYJa7Lwu2qugysxnAx0A/M8s3syuDrqkBnAhMINRiWxj++UbQRUVZF+AdM1tMqOHylrsnxTC+JNMJ+MDMFgFzgFfc/fVIHTyuhi2KiMjhxVULXUREDk+BLiKSIBToIiIJQoEuIpIgFOgiIhEQ6Yn0zOxP4YnaVpjZ3eHhvEekQBcRiYypROhmMDM7gdDw3cHAQGAkcEpN71OgS9wzs3YHjVffYmYF4d/3mtl9UTjfVDNba2aTj7DPyWa2PJmmPU521U2kZ2a9zez18Lwt75tZ/9oeDmgCpAONCU3WtrWmN6XWsWaRmOPuO4ChAGb2a2Cvu/85yqf9X3d/+gg1vR++GUo3ByW3KcBkd19tZqOA+4DTa3qTu39sZu8AmwnNYXWvu6+o6X0KdElYZnYqcKO7nxMO+p6E7sjsC1wPjCY0QVIB8E13LzezEcCdQAtgO3C5u2+u4Tz/A/wKqASK3T1iky1J/ApPLncC8NRB3d+Nw9vOB26t5m0F7n6mmR0LDCA0XxXAW2b2lfBfAYelQJdk0hs4DTiO0HQKF7j7T83sOeBsM3sFuAcY5+6FZnYR8DtgYg3HvQU4090LPp+PRYRQl3ZReAbNLwhPOHekSefOA2Z/Pke+mb1GqAFyxEBXH7okk9fcvRxYAqQAn8+hsQToAfQjdAHqrfC85Dfz3xbSkXwITDWzq8PHFSE85fPa8F9wWMiQWr59A3CKmaWGZx49hdD8VUekFrokk1IAd68ys3L/70RGVYS+CwYsc/cxdTmou08O94+eDSw0s6Hhfn1JIuGJ9E4F2ptZPqFuuEuA+83sZkIXNmcCi2pxuKcJ9bUvIXSB9HV3f6mmNynQRf5rJdDBzMaEL0qlAX1rmtHTzHqHVxf6xMy+SWjOfgV6knH38YfZVOehjO5eCXyvru9ToIuEuXuZmV0I3G1mrQh9P/4K1DRF8x1m1odQC/9tatcCE4k4TZ8rUkdmNhV4+UjDFsP79Qjvl9BLyUns0EVRkborBm6r6cYi4CVCQx9FGoRa6CIiCUItdBGRBKFAFxFJEAp0EZEEoUAXEUkQ/x/o9yg8nrV0XwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Construct a matrix to reshape Q \n",
    "R = np.identity(ndofs)\n",
    "R = R[fx, 0:ndofs]\n",
    "\n",
    "# Set Dimensions and initial conditions of state vector\n",
    "nfdofs = len(free_dofs)\n",
    "udofs = np.arange(0, nfdofs)\n",
    "vdofs = np.arange(nfdofs, 2*nfdofs)\n",
    "q0 = np.zeros((2*nfdofs))\n",
    "\n",
    "# Time span (output purposes)\n",
    "tf = 5e-8\n",
    "tspan = np.arange(0, tf, 1e-10)\n",
    "\n",
    "# Solve\n",
    "def odefun(t, q):\n",
    "    return np.append(q[vdofs], np.linalg.solve(Mii, np.transpose(np.dot(R, Q).reshape(1,-1) + RHS - np.dot(Kii, q0[udofs])))).tolist()\n",
    "\n",
    "sol = scp.solve_ivp(fun=odefun, t_span=[tspan[0], tspan[-1]], y0=q0, t_eval=tspan)  \n",
    "\n",
    "# Plot results\n",
    "plt.figure()\n",
    "plt.plot(sol.t,sol.y[1])\n",
    "plt.xlabel('Time [s]')"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "8661b582be0f012a274b7385348a70dad23ce1ee0d922b32977c503b98c884f2"
  },
  "kernelspec": {
   "display_name": "Python 3.9.7 ('base')",
   "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.9.7"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}