{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6d4ffa5a",
   "metadata": {},
   "source": [
    "Tutorial: prescribed displacement at the tip of a cantilever\n",
    "===\n",
    "\n",
    "A cantilever of 10 `BeamC` elements, clamped at $x = 0$, whose free end is\n",
    "given the transverse deflection $\\delta$ with no applied load. The prescribed\n",
    "deflection is imposed by partitioning the degrees of freedom (dofs) into the\n",
    "free ones, `bu`, and the known ones, `bk`, with the prescribed values `uk`,\n",
    "see the user guide page *Linear static analyses*. The deflected shape and\n",
    "the reactions, at the tip where the deflection is prescribed and at the\n",
    "clamp, are compared with the closed form of the Timoshenko beam.\n",
    "\n",
    "The problem is solved in two ways: by the partition written by hand with\n",
    "SciPy, and with `Analysis(..., bu=bu, uk=uk)` in a single linear step.\n",
    "\n",
    "Reference:\n",
    "\n",
    "- Bathe, K. J., *Finite Element Procedures*, 2nd ed., K. J. Bathe,\n",
    "  Watertown, MA, 2014: Sec. 4.2.2, Eqs. (4.42), (4.43) and (4.45), the\n",
    "  partition of the prescribed displacements; Sec. 1.2, Eq. (1.2), the tip\n",
    "  deflection of a cantilever with shear deformation.\n",
    "\n",
    "The model and the checks are those of\n",
    "`tests/test_prescribed_displacement_linear.py`, imported below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5a185c32",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-30T23:31:28.132454Z",
     "iopub.status.busy": "2026-10-30T23:31:28.132302Z",
     "iopub.status.idle": "2026-10-30T23:31:28.691424Z",
     "shell.execute_reply": "2026-10-30T23:31:28.690731Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "sys.path.insert(0, os.path.join('..', '..', 'tests'))\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.sparse import csr_matrix\n",
    "\n",
    "from pyfe3d import DOF\n",
    "from pyfe3d.solver import solve, Analysis\n",
    "from test_prescribed_displacement_linear import (cantilever_prop, cantilever,\n",
    "    boundary_conditions, timoshenko_cantilever, tip_force, internal_forces,\n",
    "    check_reactions)\n",
    "import test_prescribed_displacement_linear as tpd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa1d9be7",
   "metadata": {},
   "source": [
    "Model\n",
    "---\n",
    "\n",
    "Rectangular aluminium cross section, width $b = 0.1$ m along $y$, height\n",
    "$h = 0.2$ m along $z$, length $L = 1$ m ($L/h = 5$), $E = 70$ GPa,\n",
    "$\\nu = 0.3$, prescribed tip deflection $\\delta = 10$ mm along $Z$.\n",
    "`BeamProp` applies no shear correction factor, so the factor $k = 5/6$ of the\n",
    "rectangle is put in the transverse shear modulus of the bending plane $xz$,\n",
    "`prop.G13` $= k G$. `BeamC` uses it as $k G A$ = `G13*A`, in its shear\n",
    "stiffness and in the shear parameter $\\alpha_z = -12 E I_{yy}/(G_{13} A\n",
    "L_e^2)$ of its shape functions, with $L_e$ the element length. The bending\n",
    "stiffness of the deflection $w$ is $E I_{yy}$, with $I_{yy} = b h^3/12$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6702d831",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-30T23:31:28.693573Z",
     "iopub.status.busy": "2026-10-30T23:31:28.693358Z",
     "iopub.status.idle": "2026-10-30T23:31:28.699875Z",
     "shell.execute_reply": "2026-10-30T23:31:28.699197Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = 0.02 m2, Iyy = 6.66667e-05 m4, E = 7e+10 Pa, G13 = k G = 2.24359e+10 Pa\n",
      "N = 66 dofs, 59 free, 7 known; uk[tip_w] = 0.01 m\n"
     ]
    }
   ],
   "source": [
    "prop, L, delta = cantilever_prop()\n",
    "print('A = %g m2, Iyy = %g m4, E = %g Pa, G13 = k G = %g Pa'\n",
    "      % (prop.A, prop.Iyy, prop.E, prop.G13))\n",
    "nodes, ncoords, elements, KC0 = cantilever(prop, L, num_elements=10)\n",
    "N = KC0.shape[0]\n",
    "bk, bu, uk, tip_w = boundary_conditions(N, delta)\n",
    "f = np.zeros(N) # no applied load\n",
    "print('N = %d dofs, %d free, %d known; uk[tip_w] = %g m'\n",
    "      % (N, bu.sum(), bk.sum(), uk[tip_w]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9939cee7",
   "metadata": {},
   "source": [
    "Closed form\n",
    "---\n",
    "\n",
    "For a Timoshenko cantilever loaded by a tip force $P$, the bending moment is\n",
    "$M(x) = P (L - x)$ and the shear force $Q = P$. With the rotation of the\n",
    "cross section $\\theta$, $E I \\theta_{,x} = M$, and the shear strain is\n",
    "$\\gamma = w_{,x} - \\theta = P/(k G A)$. Integrating with the clamp\n",
    "conditions $w(0) = \\theta(0) = 0$:\n",
    "\n",
    "$$\n",
    "\\theta(x) = \\frac{P}{E I}\\left(L x - \\frac{x^2}{2}\\right), \\qquad\n",
    "w(x) = \\frac{P}{E I}\\left(\\frac{L x^2}{2} - \\frac{x^3}{6}\\right)\n",
    "     + \\frac{P x}{k G A}\n",
    "$$\n",
    "\n",
    "At $x = L$, $w(L) = P \\left(\\frac{L^3}{3EI} + \\frac{L}{kGA}\\right)$, which is\n",
    "Eq. (1.2) of Bathe (2014). Prescribing $w(L) = \\delta$ gives the tip\n",
    "reaction\n",
    "\n",
    "$$\n",
    "P = \\frac{\\delta}{\\dfrac{L^3}{3 E I} + \\dfrac{L}{k G A}}\n",
    "\\quad \\xrightarrow{k G A \\to \\infty} \\quad\n",
    "P_{EB} = \\frac{3 E I \\delta}{L^3}\n",
    "$$\n",
    "\n",
    "the Euler-Bernoulli value. The clamp reactions balance it: the force\n",
    "$-P$ along $Z$ and the moment $M_Y = P L$, from the moment equilibrium about\n",
    "the clamp, $M_Y \\boldsymbol{e}_Y + (L, 0, 0) \\times (0, 0, P) = \\boldsymbol{0}$. In\n",
    "pyfe3d the rotation about $Y$ is $r_y = -\\theta$, since the axial\n",
    "displacement of a point at the height $z$ is $u + z r_y$ (*Conventions*,\n",
    "sign conventions). The derivation is verified with SymPy by\n",
    "`test_closed_form_timoshenko_cantilever` of the test file, run here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d5342a87",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-30T23:31:28.701869Z",
     "iopub.status.busy": "2026-10-30T23:31:28.701723Z",
     "iopub.status.idle": "2026-10-30T23:31:29.173052Z",
     "shell.execute_reply": "2026-10-30T23:31:29.172247Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Timoshenko      P    = 1.357642e+05 N\n",
      "Euler-Bernoulli P_EB = 1.400000e+05 N, 3.12 % larger\n",
      "shear/bending flexibility 3EI/(L^2 kGA) = 0.0312\n"
     ]
    }
   ],
   "source": [
    "tpd.test_closed_form_timoshenko_cantilever()\n",
    "P, P_EB = tip_force(delta, prop, L)\n",
    "print('Timoshenko      P    = %.6e N' % P)\n",
    "print('Euler-Bernoulli P_EB = %.6e N, %.2f %% larger'\n",
    "      % (P_EB, 100*(P_EB/P - 1)))\n",
    "print('shear/bending flexibility 3EI/(L^2 kGA) = %.4f'\n",
    "      % (3*prop.E*prop.Iyy/(L**2*prop.G13*prop.A)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7394d3a",
   "metadata": {},
   "source": [
    "The shear flexibility of this beam is 3.1 % of its bending flexibility,\n",
    "$3 E I/(L^2 k G A) = \\frac{E}{4 k G}(h/L)^2$ for the rectangle: the\n",
    "Euler-Bernoulli tip force is 3.1 % larger than the Timoshenko one, a\n",
    "difference that is visible below but small.\n",
    "\n",
    "1. Partition by hand with SciPy\n",
    "---\n",
    "\n",
    "The dofs need not be reordered: the Boolean arrays select the blocks. The\n",
    "free dofs solve\n",
    "\n",
    "$$\n",
    "\\boldsymbol{K}_{{C0}_{uu}} \\boldsymbol{u}_u = \\boldsymbol{f}_u - \\boldsymbol{K}_{{C0}_{uk}} \\boldsymbol{u}_k\n",
    "$$\n",
    "\n",
    "and the reactions follow from the second row of the partition,\n",
    "$\\boldsymbol{K}_{{C0}_{ku}} \\boldsymbol{u}_u + \\boldsymbol{K}_{{C0}_{kk}} \\boldsymbol{u}_k - \\boldsymbol{f}_k$,\n",
    "equal to the internal forces of `update_fint` at the known dofs minus the\n",
    "loads applied there, $\\boldsymbol{F}_{{int}_k} - \\boldsymbol{f}_k$. These are Eqs. (4.43)\n",
    "and (4.45) of Bathe (2014), with the known set written $k$, as in pyfe3d."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1e889111",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-30T23:31:29.175894Z",
     "iopub.status.busy": "2026-10-30T23:31:29.175574Z",
     "iopub.status.idle": "2026-10-30T23:31:29.182539Z",
     "shell.execute_reply": "2026-10-30T23:31:29.181886Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tip reaction  FE 1.3576415826e+05 N, closed form 1.3576415826e+05 N, rel. diff. 1.4e-14\n",
      "clamp forces  [N]    [      0.               0.         -135764.15826221]\n",
      "clamp moments [N m]  [     0.         135764.15826222      0.        ]  P L = 135764.1582622188\n"
     ]
    }
   ],
   "source": [
    "KC0uu = KC0[bu, :][:, bu]\n",
    "KC0uk = KC0[bu, :][:, bk]\n",
    "KC0ku = KC0[bk, :][:, bu]\n",
    "KC0kk = KC0[bk, :][:, bk]\n",
    "\n",
    "u = np.zeros(N)\n",
    "u[bk] = uk[bk]\n",
    "u[bu] = solve(KC0uu, f[bu] - KC0uk @ uk[bk])\n",
    "\n",
    "R_partition = KC0ku @ u[bu] + KC0kk @ uk[bk] - f[bk]\n",
    "R = internal_forces(elements, prop, u, ncoords) - f\n",
    "assert np.allclose(R[bk], R_partition, rtol=1e-12, atol=1e-12*P*L)\n",
    "\n",
    "print('tip reaction  FE %.10e N, closed form %.10e N, rel. diff. %.1e'\n",
    "      % (R[tip_w], P, abs(R[tip_w]/P - 1)))\n",
    "print('clamp forces  [N]   ', R[0:3])\n",
    "print('clamp moments [N m] ', R[3:6], ' P L =', P*L)\n",
    "check_reactions(R, nodes, L, P, rtol=1e-9)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c32518ee",
   "metadata": {},
   "source": [
    "`check_reactions` asserts, with a relative tolerance of $10^{-9}$, the tip\n",
    "reaction $P$, the clamp force $(0, 0, -P)$, the clamp moment $M_Y = P L$ with\n",
    "the moment equilibrium of all the reactions about the origin, and zero\n",
    "forces at the free dofs. `BeamC` solves the homogeneous equations of the\n",
    "Timoshenko beam exactly, so the agreement is to round-off.\n",
    "\n",
    "The deflected shape: the nodal deflections and rotations against the closed\n",
    "form, with the tolerance $10^{-10} \\delta$, and $w(x)$ plotted with the\n",
    "Euler-Bernoulli shape of the same tip deflection,\n",
    "$w_{EB}(x) = \\delta \\left(\\frac{3 x^2}{2 L^2} - \\frac{x^3}{2 L^3}\\right)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "aa422581",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-30T23:31:29.184368Z",
     "iopub.status.busy": "2026-10-30T23:31:29.184221Z",
     "iopub.status.idle": "2026-10-30T23:31:29.334505Z",
     "shell.execute_reply": "2026-10-30T23:31:29.333635Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max |w_FE - w| / delta          = 2.2e-15\n",
      "max |r_y,FE + theta| / theta(L) = 7.4e-15\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAGGCAYAAACNCg6xAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbBRJREFUeJzt3XdcU9f7B/BPEvaUPQRRFPfeew9UXFXrKqJ1z6ptHdU6Wq2j31q1VazbWldt1Q4VpdZZFRS1DnAhDioKArJn7vn9YclPBG2Yl8Dn/XrlpTn35uTJQ0Iezj33XIUQQoCIiIiI/pNS7gCIiIiIdAULJyIiIiItsXAiIiIi0hILJyIiIiItsXAiIiIi0hILJyIiIiItsXAiIiIi0hILJyIiIiIt6ckdQHGRJAlPnjyBubk5FAqF3OEQERFRCSGEQEJCApydnaFUvn1MqcwUTk+ePIGrq6vcYRAREVEJ9fjxY7i4uLx1nzJTOJmbmwN4mRQLC4tC71+SJERFRcHOzu4/q1UqXMy9PJh3+TD38mHu5VHUeY+Pj4erq6umVnibMlM4ZR2es7CwKLLCKTU1FRYWFvwwFTPmXh7Mu3yYe/kw9/IorrxrM5WHP3UiIiIiLbFwIiIiItISCyciIiIiLZWZOU7aUqvVyMjIyPPjJElCRkYGUlNTedy7mDH38iiMvOvr60OlUhVyZERERYeF07+EEHj69ClevHiR78dLkoSEhASuE1XMmHt5FFbey5UrB0dHR/7siEgnsHD6V1bRZG9vDxMTkzz/EhdCIDMzE3p6evwCKGbMvTwKmnchBJKTkxEZGQkAcHJyKuwQiYgKHQsnvDw8l1U02djY5KsPfnnLh7mXR2Hk3djYGAAQGRkJe3t7HrYjohKvREwIOX36NHr16gVnZ2coFAocPHgw23YhBBYuXAhnZ2cYGxujffv2uHnzZqE9f9acJhMTk0Lrk4i0k/W5y8/cQiIq/dSSwIX70Th2KwYX7kdDLQlZ4ykRI05JSUmoV68eRo4cif79++fYvmLFCqxcuRLbtm1D1apVsXjxYnTp0gW3b9/WapVPbXG0gqj48XNHRG/idyMCiw+FIDw25d+WMLhYGWNezxrwrC3P4f0SUTh1794d3bt3z3WbEAKrVq3C3Llz8c477wAAtm/fDgcHB+zatQvjxo0rzlCJiIioGPjdiMCEnZfRqbo9Vg+qBytlKmIlI/ieuo8JOy/Dd1hDWYqnEnGo7m3CwsLw9OlTdO3aVdNmaGiIdu3a4dy5czJGphsWLlyI+vXrl4kYSsJrze1Qc0mLQQiBsWPHwtraGgqFAlevXi222IiItKGWBBYfCkGn6vbY4N0YFUwlbNv0Heq7lsMG78boVN0eSw6HyHLYrkSMOL3N06dPAQAODg7Z2h0cHPDw4cM3Pi4tLQ1paWma+/Hx8QBerj0jSVK2fSVJghBCc8uvrMcWpI+8+K+1c3x8fPDNN99g8uTJxRZTboojL6/2LedrzXr+khzDkSNHsG3bNpw4cQLu7u6wtbXNd7yF8bPNijW3zyblLut3FvNV/Jj74hFwPxrhsSlYPagewsMfo1mzZnj69Cn09PQwb948jG/njgHrLyDg/nM0d8/fSV2vysvPs8QXTllenwchhHjr3IilS5di0aJFOdqjoqKQmpqarS0jIwOSJCEzMxOZmZn5ik8IAbVanWusReXRo0ea/+/btw+LFi3CjRs3NG3GxsYwMjKCkZFRvl9XYcj6RVNUMWTlvqifR1tqtbpEx3D37l04OTmhadOmmra8xpuV56z3ekHe85mZmZAkCdHR0dDX1893P2WJJEmIi4uDEIKLvhYz5r543AuPAQDEPghGxz5eiI+Ph5mZGdq2bYvIyEhYKdX/7hcFdzN1gZ8vISFB+51FCQNAHDhwQHM/NDRUABCXL1/Otl/v3r3F8OHD39hPamqqiIuL09weP34sAIjY2FihVquz3ZKSksTNmzdFcnKykCQp37e0tLQCPb4gty1btghLS8sc7fPnzxf16tXT3Pfx8RF9+vQRixcvFvb29sLS0lIsWLBApKeniw8//FBYWVmJ8uXLi02bNmXr5++//xYdOnQQRkZGwtraWowePVrEx8drtv/555+iSZMmwsTERFhaWoqWLVuKsLCwbDFs375duLm5CQsLCzFo0CARFxenebxarRbLli0TlSpVEkZGRqJu3brixx9/zNY/AOHv7y8aNWokjI2NRYsWLURISIgm96+/1tDQUFG5cmUxbtw4kZmZKaKjo4W3t7coV66cMDY2Fp6enuL27dt5zvWmTZtEzZo1hYGBgXB0dBQTJ07UbAMg9u/fXyh5kyRJ/PLLL6Jhw4bC0NBQVKpUSfOzytp++/Zt0aZNG2FoaChq1Kghjh49miOGV28+Pj4CgObm5uYmJEkSKSkpYvLkycLOzk4YGhqKVq1aiYCAgBz5P3LkiGjUqJHQ19cXx48fF23bthWTJk0SU6dOFeXKlRP29vZi/fr1IiEhQYwYMUKYmZkJd3d3cejQoTfmMzk5Wdy8eVMkJSXl+GzylvstIyNDPHnyRGRkZMgeS1m7MffFc/vrbqRwm/W7sK3WWAAQNjY2IjAwUJP3i2HPhdus38VfdyML5fliY2MFABEXF/efdUqJH3GqVKkSHB0d4e/vjwYNGgAA0tPTcerUKSxfvvyNjzM0NIShoWGOdqVSmeOvBKVSCYVCobkB/784n7ZEIa4llJ8FON/0l39u7X/++SdcXFxw+vRp/PXXXxg1ahQuXLiAtm3bIiAgAHv37sWECRPQtWtXuLq6Ijk5Gd27d0fz5s1x8eJFREZGYvTo0ZgyZQq2bduGzMxM9OvXD2PGjMHu3buRnp6OwMDAbHkNDQ3FL7/8gt9//x2xsbF49913sXz5cixZsgQAMG/ePOzfvx++vr7w8PDA6dOn4e3tDXt7e7Rr104T/7x58/DVV1/Bzs4O48ePx6hRo3D27Nlsr0+hUODGjRvo2rUrfHx8sHTpUgDAyJEjcffuXfz666+wsLDArFmz0LNnTwQHB2s90uHr64sZM2Zg2bJl6N69O+Li4vDXX3/leH6FQlHgvB09ehTe3t5Ys2YN2rRpg9DQUIwdOxYKhQILFiyAJEno378/bG1tceHCBcTHx2PatGnZYnjd6tWrUblyZWzYsAEXL16ESqWCQqHArFmzsH//fmzfvh1ubm5YsWIFPD09ce/ePc1cKACYNWsW/ve//8Hd3R2WlpYAgO+//x4zZ85EYGAg9u7di4kTJ+KXX35Bv3798Mknn+Drr7/G8OHD8ejRo1yX/MiKNbfPJr0ZcyYf5r7oJT+8BnV8JBS1uqOGKgV/Hv8DQNYUFQXWn7oPV2tjNHO3hVJZ8KM8efpZ5nNgqFAlJCSIK1euiCtXrggAYuXKleLKlSvi4cOHQgghli1bJiwtLcX+/fvF9evXxZAhQ4STk5OIj4/X+jni4uLeWE2mpKSI4OBgkZKSomlLTEzM9pd5cd4SExPznMOtW7cKS0vLHO0LFiwQ9erV09z38fERbm5uQq1Wa9qqVasm2rRpo7mfmZkpTE1Nxe7du4UQQmzYsEFYWVlli+vQoUNCqVSKp0+fiujoaAFAnDx5MtfYFixYIExMTLL9vD7++GPRrFkzTa6NjIzEuXPnsj1u1KhRYsiQIUIIIU6cOCEAiD/++CNbDABEcnKySE9P14w4nTt3TlhbW4svv/xSs++dO3cEAPHXX39p2p4/fy6MjY3Fjz/+mGvcuXF2dhZz585943a8MmJa0Ly1adNGfPHFF9naduzYIZycnIQQQhw9elSoVCrx+PFjzfYjR47kGLV93ddffy3c3Nw09xMTE4W+vr7YuXOnpi09PV04OzuLFStWCCH+P/8HDx7U7CNJkmjbtq1o3bq1pi3rvePt7a1pi4iIEADE+fPnc40nt88fvZ1arRYRERHZPsdUPJj7ovfTTz8JpVIpjKu2EG4zfxPDN54TF8Oei3sPw8XFsOdi1LZAUXH27+LI9SeF9pxvqxFeVyJGnC5duoQOHTpo7s+YMQPAy8nN27Ztw8yZM5GSkoKJEyciNjYWzZo1w7Fjxwp1DaeypFatWtmqawcHB9SuXVtzX6VSwcbGRnMpjJCQENSrVw+mpqaafVq1agVJknD79m20bdsWI0aMQLdu3dClSxd07twZ7777brZLaFSsWDHbz8vJyUnTf3BwMFJTU9GlS5dscaanp2tGGbPUrVs3Wx/Ay1WnnZ2dAbyc99W5c2csXrwY06dP1+wbEhICPT09NGvWTNNmY2ODatWqISQkRKu8RUZG4smTJ+jUqZNW+xc0b0FBQbh48aJmVA54OXcpNTUVycnJCAkJQYUKFeDi4qLZ3qJFC61ie1VoaCgyMjLQqlUrTZu+vj6aNm2aIzeNGzfO8fg6depo/p/13nm1LevEjqyfNxHRm2zYsEGzzFCdcmp8MqgOvvQPxYD1FzT7uFoby7YUAVBCJoe3b9/+rWflKBQKLFy4EAsXLiy2mExMTJCYmKj1/qKQD9UVpdcPSykUilzbss4yEG+ZiJ/VvnXrVkydOhV+fn7Yu3cv5s2bB39/fzRv3vyNz5nVf9a/hw4dQvny5bPt9/rh1lf7yXruV8+GsLOzg7OzM/bs2YNRo0bBwsJC8xpy87bX9rqsy4Noq6B5kyQJixYt0qxf9iojI6NcX1N+rxmX22Nzi//VIjDLf72fcvs5ERG97osvvsDcuXMBAF27dsXvv/8OfX19eNWvgID7z3EvPApVXOzQzN0WqkI4PJdfPED7BgqFAqamprLcStpKyjVr1sTVq1eRlJSkafvrr7+gVCpRtWpVTVuDBg0wZ84cnDt3DrVr18auXbu07t/Q0BCPHj1ClSpVst1cXV3zFKuxsTF+//13GBkZoVu3bpozJWrWrInMzEwEBARo9o2OjsadO3dQo0YNrfo2NzdHxYoVcfz4ca1fV0Hy1rBhQ9y+fTtHTqpUqQKlUomaNWvi0aNHePLkiaav8+fPaxXbq6pUqQIDAwOcPXtW05aRkYFLly5pnRsiovwSQmD+/Pmaomno0KE4cuSI5g8wlVKB5u426FrdGs3dbWQtmgAWTqSFYcOGwcjICD4+Prhx4wZOnDiBKVOmwNvbGw4ODggLC8OcOXNw/vx5PHz4EMeOHctzQfLRRx9h+vTp2L59O0JDQ3HlyhWsXbsW27dvz3O8pqamOHToEPT09NC9e3ckJibCw8MDffr0wZgxY3D27Fn8/fffeO+991C+fHn06dNH674XLlyIr776CmvWrMHdu3dx+fJlfPPNN7nuW9C8zZ8/H99//z0WLlyImzdvIiQkRDMqBQCdO3dGtWrVMHz4cPz99984c+aM5hdPXvM1YcIEfPzxx/Dz80NwcDDGjBmD5ORkjBo1Ks/9ERFpS5IkTJo0CZ9//jmAl1N1fvjhhxI98b7kRkYlhomJCY4ePYqYmBg0adIEAwYMQKdOnfDtt99qtt+6dQv9+/dH1apVMXbsWEyePDlPl8P5/PPPMX/+fCxduhQ1atRAt27d8Ntvv6FSpUr5itnMzAxHjhyBEAI9evRAUlIStm7dikaNGsHLywstWrSAEAKHDx/OcVhp27Ztb+zXx8cHq1atwrp161CrVi14eXnh7t27ue5b0Lx169YNv//+O/z9/dGkSRM0b94cK1euhJubG4CXZ4EcOHAAaWlpaNq0KUaPHp1tPlReLFu2DP3794e3tzcaNmyIe/fu4ejRo7CysspXf0RE/yU9PR1eXl7w9fUFAKxduxZfffVViTvq8jqFeNvkolIkPj4elpaWiIuL08x7yZKamoqwsDBUqlQJRkZG+eq/MOc4Ud4UVu4fPHgADw8PBAcHw8PDoxAjLJ0KK++F8fkrayRJQmRkJOzt7Uv0X+alEXNfOJKTk9G5c2fN9AJfX1+MHz/+jfsXdd7fViO8rkRMDicqCfz8/DB27FgWTURERSg2NhatW7dGcHAwlEol9u7diwEDBsgdltZYOBH9621/7RARUcFFRESga9euCA4Ohp6eHo4cOYLOnTvLHVaesHAiIiKiInf//n107twZYWFhsLW1hZ+fHxo1aiR3WHnGwomIiIiK1LVr19C6dWskJCSgYsWK+OOPP1C5cmW5w8oXzmwjIiKiInPmzBk0bdoUCQkJsLW1xV9//aWzRRPAwomIiIiKyG+//YYOHTogLS0Nbm5uuH37tuYSWbqKhRMREREVuh9++AF9+/aFWq1GrVq1EBwcDGtra7nDKjAWTkRERFSoVq1aBW9vb0iShJ49e+LKlStFfh3W4sLCiYiIiAqFEAKzZ8/G9OnTAQBTp07Fr7/+muNi4LqMZ9URERFRganVaowePVpz2ar3338fq1atKnVX0+CIEwEA2rdvj2nTpskdRon2eo7+6z4RUVmRnp6Ofv36aYqm2bNnY/PmzaWuaAJYOJUKI0aMgEKhyHHz9PSUOzQAOeOzsbGBp6cnrl27JndohWr//v2aK3yXRk+fPsXQoUPh6OgIAwMDlC9fHitXrpQ7LCKSWWJiInr16oXffvsNALBy5UosXbpU5qiKDg/VlRKenp7YunVrtjZDQ8Nie/709HQYGBi8cfur8T19+hTz5s2Dl5cXHj16VGTPWdxKw9kibzNu3DikpaXhjz/+gJWVFZ4+fYro6Gi5wyIiGUVHR6NHjx4IDAyEsbExvvvuO3h7e8sdVpHiiFMpYWhoCEdHx2w3KysrAEDFihWxatWqbPvXr18fCxcufGN/QgisWLEC7u7uMDY2Rr169fDTTz9ptrdv3x6TJ0/GjBkzYGtriy5dumgdX/369TFr1iw8fvwYUVFRBX7ODh06YPr06Zg5cyasra3h6OiY47WlpaVh6tSpsLe3h5GREVq3bo2LFy9qtucnR6/L66G6n376CXXq1IGxsTFsbGzQuXNnJCUlabb7+fmhdevWKFeuHGxsbODl5YXQ0NBszzdlyhRMmzYNVlZWcHBwwIYNG5CUlISRI0fC3NwclStXxpEjR7I973/l+U3S0tLw4MEDnD9/Hunp6WjYsCE6dOig9eslotLln3/+QePGjREYGIhy5crhzz//LPVFE8DCid5g3rx52Lp1K3x9fXHz5k1Mnz4d7733Hk6dOqXZZ/v27dDT08Nff/2F7777Tuu+ExMTsXPnTlSpUgU2NjaF8pw7duyAqakpAgICsGLFCnz22Wfw9/fXbJ85cyZ+/vlnbN++HZcvX0aVKlXQrVs3xMTE5DdFBRIREYEhQ4bg/fffR0hICE6ePIl33nkHQgjNPklJSZgxYwYuXryI48ePQ6lUol+/fpAkSbPP9u3bYWtri8DAQEyZMgUTJkzAwIED0bJlS1y+fBndunWDt7c3kpOTNY/RJs+vy8zMhKenJ/bu3QtPT0+sXbsWvXr1QkJCQtEkiIhKtDt37qBhw4Z48OABDAwMcOTIETRv3lzusIqHKCPi4uIEABEXF5djW0pKiggODhYpKSn57l+SJJGeni4kSSpImPni4+MjVCqVMDU1zXb77LPPhBBCuLm5ia+//jrbY+rVqycWLFigud+uXTvxwQcfCCGESExMFEZGRuLcuXPZHjNq1CgxZMgQzf7169fPV3wAhJOTkwgKCtLsU5DnbNeunWjVqlW23Ddp0kTMmjVL07e+vr7YuXOnZnt6erpwdnYWK1asEELkPUfa3H+boKAgAUA8ePBAq/2FECIyMlIAENevX9c8X+vWrTXbMzMzhampqfD29ta0RURECADi/PnzQgjt8pybiRMnin379mVrq1ixoliyZEmB3/OF8fkra9RqtYiIiBBqtVruUMoc5v7l7y8LCwsBQJiYmIjg4OAif86izvvbaoTXcY7Tf4iIiEBERES2NisrK1SqVAmpqakIDg4G8PLwh1qthkql0lzt+fbt29kOvQAvDwlZW1sjKioKjx8/zrbN3NwcHh4e+YqzQ4cO8PX1zdaW3zk3wcHBSE1NzXH4LT09HQ0aNNDcb9y4cbbtO3fuxLhx4zT3jxw5gjZt2uSILyYmBuvWrUP37t0RGBgINze3fD9nljp16mS77+TkhMjISABAaGgoMjIy0KpVK812fX19NG3aFCEhIW9PRhGpV68eOnXqhDp16qBbt27o2rUrBgwYoDm8CryM+9NPP8WFCxfw/PlzzUjTo0ePULt2bQBA3bp1NfurVCrY2Nhky4WDgwMAaHKhbZ5fdeXKFfzwww85DmVaWlrm+GwQUel28uRJeHl5ISkpCeXKlcP169fh4uIid1jFioXTf/juu++waNGibG3Dhg3DDz/8gPDwcE2R9Crx7+GWESNG4MKFC9m27dixA++99x5+/PFHTJ48Odu2rl274ujRo/mK09TUFFWqVMl1m1KpzHYICAAyMjLe2FfWF/ShQ4dQvnz5bNtenXBuamqabVvv3r3RrFkzzf1XH/t6fI0aNYKlpSU2btyIxYsX5/s5s7y+uJpCodD0mfXaXz8tVgihactrjgpKpVLB398f586dw7Fjx/DNN99g7ty5CAgIQKVKlQAAvXr1gqurKzZu3AhnZ2dIkoTatWsjPT1d009ur/vVtqzXl5ULbfP8qv3796Nq1arZ+k1OTsbt27cxduzY/KaAiHTMwYMHMWjQIM0cx6wTRcoaFk7/Ydy4cejdu3e2tqw3iouLC4KCggBkH3HKsm3btlxHnADg3XffRYsWLbJtMzc3L+zwAQB2dnbZRgbi4+MRFhb2xv1r1qwJQ0NDPHr0CO3atdP6eczNzbV+DQqFAkqlEikpKQV6Tm1UqVIFBgYGOHv2LIYOHQrgZVF06dIlzWTuvOaoMCgUCrRq1QqtWrXC/Pnz4ebmhgMHDmDGjBmIjo5GSEgIvvvuO82o3dmzZwv8nPnJc2xsbI738caNGyGEwDvvvFPgmIio5NuyZQtGjRoFAOjTpw/27NkDIyMjmaOSBwun/+Dk5AQnJ6dctxkZGaFhw4YAXhZOmZmZ0NP7/5RWq1btjf3a2dnBzs6u0OJMS0vD06dPs7Xp6enB1tYWHTt2xLZt29CrVy9YWVnh008/zVbgvc7c3BwfffQRpk+fDkmS0Lp1a8THx+PcuXMwMzODj49PgeKLjY3Ft99+q1n7o6ieM4upqSkmTJiAjz/+GNbW1qhQoQJWrFiB5ORkzS+CvOaooAICAnD8+HF07doV9vb2CAgIQFRUFGrUqAHgZXFuY2ODDRs2wMnJCY8ePcLs2bML/Lz5yXOzZs2wdu1afP311/Dy8sLRo0cxe/ZsrFmzJtvkfiIqnZYvX675/dOoUSPs27evVF1CJa9YOJUSfn5+OQq8atWq4datW5gzZw7u378PLy8vWFpa4vPPP//P0ZTPP/8c9vb2WLp0Ke7fv49y5cqhYcOG+OSTTwocn7m5OapXr459+/ahffv2Rfacr1q2bBkkSYK3tzcSEhLQuHFjHD16VDN6mJ8c/Zdt27Zh5MiROQ4BAoCFhQVOnz6NVatWIT4+Hm5ubvjqq6/QvXt3AC8PHe7ZswdTp05F7dq1Ua1aNaxZsyZbvvIrr3l+77338OjRI6xZswYLFixA7dq1sW/fPvTs2ROZmZkFjoeISiYhBD7++GN89dVXAF5OJ/Hz8yuVq4HnhULk9lu9FIqPj4elpSXi4uJgYWGRbVtqairCwsJQqVKlfA89vjriVNbfVMWtpOZ+4cKFOHnyJE6ePCl3KEWisPJeGJ+/skaSJERGRsLe3h5KJVeVKU5lJfeZmZkYP348Nm/eDODlH1A7duyQLZ6izvvbaoTXld6fOpHMjh49ihUrVsgdBhFRnqSmpmLgwIGaa83NnTtX1qKppOGhOqIicv78eblDICLKk/j4ePTo0QN//fUX9PX18eOPP6Jv375yh1WisHAiIiIiREZGomPHjrh58yYUCgXWr1/PoikXLJyIiIjKuIcPH6Jdu3Z4+PAhlEolDh48qDnrmbJj4URERFSGBQcHo0uXLnjy5An09fVx4sSJbFdaoOw4OZyIiKiMCggIQOvWrfHkyRO4u7vjypUrLJr+A0eciIiISjm1JBAYFoPIhFTYmxuhaSVrHP/DH71790ZaWhqaNGmCI0eOcFFbLbBwIiIiKsX8bkRg8aEQhMemaNqsDCTc3bccaWlpcHR0zLYgML0dD9URERGVUn43IjBh52VUdzTH/oktcXNRN7zn+BT/3LwIm96z4N72HYSFhbFoygOOOBEREZVCaklg8aEQdKpujw3ejaFQAEuWLMGSTz8FoIDH+1/Cocck6BsYyh2qTuGIE9F/2LZtG8qVKyd3GEREeRIYFoPw2BRM7FAFgMC0adPw6aefAgB8fIZj96cj8Tg2BYFhMfIGqmNYOOm4ESNGQKFQaG42Njbw9PTEtWvX5A5NIz09HStWrEC9evVgYmICW1tbtGrVClu3bkVGRobW/WzYsAHt27eHhYUFFAoFXrx4kWOf2NhYeHt7w9LSEpaWlvD29s51v9Lq5MmTb8wNEZUtkQmpAIBKVoYYOnQo1qxZAwBYvXo1tm3bhmpOFtn2I+2wcCpkakngfGg0frn6D86HRkMtFf01lD09PREREYGIiAgcP34cenp68PLyKvLn1UZ6ejq6deuGZcuWYezYsTh37hwCAwMxadIkfPPNN7h586bWfSUnJ8PT0xOffPLJG/cZOnQorl69Cj8/P/j5+eHq1avw9vYujJdCRKRT7M1fXjS713tjsXfvXgDAjBkzMHXqVADA7acJ2fYj7bBwKkRHbz5D+/+dxJCNF/DBnqsYsvEC2n15An43Ior0eQ0NDeHo6AhHR0fUr18fs2bNwuPHjxEVFaXZ559//sGgQYNgZWUFGxsb9OnTBw8ePNBsv3jxIrp06QJbW1tYWlqiXbt2uHz5crbnUSgU+O677+Dl5QUTExPUqFED58+fx71799C+fXuYmpqiRYsWCA0N1Txm1apVOH36NI4fP45Jkyahfv36cHd3x9ChQxEQEAAPDw+tX+e0adMwe/ZsNG/ePNftISEh8PPzw6ZNm9CiRQu0aNECGzduxO+//47bt2+/sd/09HTMnDkT5cuXh6mpKZo1a4aTJ0++NZbffvsNjRo1gpGREdzd3bFo0SJkZmYWKFfa9rtp0yb069cPJiYm8PDwwK+//goAePDgATp06AAAsLKygkKhwIgRIwAAP/30E+rUqQNjY2PY2Nigc+fOSEpKeutrJCLd5m4uQZkSizt67gAUWL58Ob766isAgCQJ+J68B1drYzStZC1voLpGlBFxcXECgIiLi8uxLSUlRQQHB4uUlJR893/42hNRcfbv4v1tgSLoYYxITM0QQQ9jxKhtgaLi7N/FketPChL+G/n4+Ig+ffpo7ickJIhx48aJKlWqCLVaLYQQIikpSXh4eIj3339fXLt2TQQHB4uhQ4eKatWqibS0NCGEEMePHxc7duwQwcHBIjg4WIwaNUo4ODiI+Ph4Td8ARPny5cXevXvF7du3Rd++fUXFihVFx44dhZ+fnwgODhbNmzcXnp6emsfUrVtXdO3atVBf84kTJwQAERsbK4QQQpIkkZ6eLjZt2iQsLS1z7G9paSm2bNnyxv6GDh0qWrZsKU6fPi3u3bsnvvzyS2FoaCju3LkjhBBi69at2fr18/MTFhYWYtu2bSI0NFQcO3ZMVKxYUSxcuFCzT35ypW2/Li4uYteuXeLu3bti6tSpwszMTERHR4vMzEzx888/CwDi9u3bIiIiQrx48UI8efJE6OnpiZUrV4qwsDBx7do1sXbtWpGQkJDPn8BLWXmXJKlA/RTG56+sUavVIiIiQvMZp+KjK7l/+PChqFatmjCu2kJUmPmr6Pr5fnHpQYxISM0Qlx4U/XdTYSvqvL+tRngdCydR8F/cmWpJtFp2XIzcEiAyM7P/UNVqSYzaFihaLz8uMtUF+4LJjY+Pj1CpVMLU1FSYmpoKAMLJyUkEBQVp9tm8ebOoVq1ati+4tLQ0YWxsLI4ePZr7a8rMFObm5uK3337TtAEQ8+bN09w/f/68ACA2b96sadu9e7cwMjLS3Dc2NhZTp04tlNea5U2F0+LFi4WHh0eO/T08PMQXX3yRa1/37t0TCoVC/PPPP9naO3XqJObMmSOEyFk4tWnTJkd/O3bsEE5OTpr7+clVfvpNTEwUCoVCHDlyJNfcCCFEUFCQACAePHiQaw7yi4WTfHTly7s00oXc37hxQzg7OwsAwtnZWSzZfki0WnZcuM36XXNrvfy4zhRNQpSswonLERSCrDMXVg6sA6VSkW2bUqnAhPZV0N/3HALDYtCicuGvytqhQwf4+voCAGJiYrBu3Tp0794dgYGBcHNzQ1BQEO7duwdzc/Nsj0tNTdUcKoqMjMT8+fPx559/4tmzZ1Cr1UhOTsajR4+yPaZu3bqa/zs4OAAA6tSpk60tNTUV8fHxsLCwgBACCkX2nBSl3J7rbTFcvnwZQghUrVo1W3taWtobV9ANCgrCxYsXsWTJEk2bWq1GamoqkpOTYWJiAiDvucpPv6ampjA3N0dkZGTuCQFQr149dOrUCXXq1EG3bt3QtWtXDBgwgOu2EJVC586dg6enJxISEuDh4YHjx4/D1dUVs3JZOVylLL7fzaUJC6dCkHVGQlV7s1y3V3M0z7ZfYTM1NUWVKlU09xs1agRLS0ts3LgRixcvhiRJaNSoEXbu3JnjsXZ2dgBenp0XFRWFVatWwc3NDYaGhmjRogXS09Oz7a+vr6/5f1YxklubJEkAgKpVqyIkJKSQXunbOTo64tmzZznao6KiNIXL6yRJgkqlQlBQEFQqVbZtZma5/zwlScKiRYvwzjvv5NhmZPT/kyzzmqv89JvVT1YfuVGpVPD398e5c+dw7NgxfPPNN5g7dy4CAgJQqVKlNz6OiHTLoUOH8M477yA9PR2mpqY4dOgQXF1dAQAqpaJI/nAvi1g4FYKsMxLuRCaiSaWcZycU95kLCoUCSqUSKSkvl9dv2LAh9u7dC3t7e1hYWOT6mDNnzmDdunXo0aMHAODx48d4/vx5gWMZOnQoPvnkE1y5cgUNGjTIti0zMxNpaWkwNTUt8PMAQIsWLRAXF4fAwEA0bdoUwMsLWMbFxaFly5a5PqZBgwZQq9WIjIxEmzZttHqehg0b4vbt29mK1cJQGP0aGBgAeDlS9SqFQoFWrVqhVatWmD9/Ptzc3HDgwAHMmDGjQDETUcnw/fffY+TIkZAkCTY2Nrh58+Yb/2CkguFZdYWgaSVruFgZY/2pMEivLT9QHGcupKWl4enTp3j69ClCQkIwZcoUJCYmolevXgCAYcOGwdbWFn369MGZM2cQFhaGU6dO4YMPPkB4eDgAoEqVKtixYwdCQkIQEBCAYcOGwdjYuMCxTZs2Da1atUKnTp2wdu1a/P3337h//z5+/PFHNGvWDHfv3tW6r6dPn+Lq1au4d+8eAOD69eu4evUqYmJeLt5Wo0YNeHp6YsyYMbhw4QIuXLiAMWPGwMvLC9WqVcu1z6pVq2LYsGEYPnw49u/fj7CwMFy8eBHLly/H4cOHc33M/Pnz8f3332PhwoW4efMmQkJCsHfvXsybNy+P2Sn8ft3c3KBQKPD7778jKioKiYmJCAgIwBdffIFLly7h0aNH2L9/P6KiolCjRo0CxUtEJcP//vc/+Pj4QJIkVKhQAWFhYSyailKRzLIqgYrzrLriPHPBx8dHANDczM3NRZMmTcRPP/2Ubb+IiAgxfPhwYWtrKwwNDYW7u7sYM2aMJh+XL18WjRs3FoaGhsLDw0Ps27dPuLm5ia+//lrTBwBx4MABzf2wsDABQFy5ckXTltvk5NTUVLF06VJRp04dYWRkJKytrUWrVq3Etm3bREZGhhDi5QTs/3o7LliwINtrzbpt2bJFM0k5OjpaDBs2TJibmwtzc3MxbNiwbLHkJj09XcyfP19UrFhR6OvrC0dHR9GvXz9x7do1TWyvn63n5+cnWrZsKYyNjYWFhYVo2rSp2LBhQ4Fzldd+hXh51uDWrVs19z/77DPh6OgoFAqF8PHxEcHBwaJbt27Czs5OGBoaiqpVq4pvvvnmrTnRBieHy0cXJiiXViUp95IkiY8//ljzu3DIkCGaM6VLm5I0OVwhhCj6FRpLgPj4eFhaWiIuLi7H4arU1FSEhYWhUqVK2eaS5IUQAof+/gfLjt7JdgVqV2tjzO1RA561nQoUf2m3cOFCnDx58j/XT8qNEAKZmZnQ09Mr1onoZV1h5b0wPn9ljSRJiIyMhL29PZRKHjgoTiUl9xkZGRg9ejS+//57AMDy5csxc+ZM2eIpakWd97fVCK/TmTlOmZmZWLhwIXbu3ImnT5/CyckJI0aMwLx580rML45utRzgWccZFx/E8syFPDp69ChWr14tdxhERCVecnIyBg4cqJlOMHDgwFJdNJU0OlM4LV++HOvXr8f27dtRq1YtXLp0CSNHjoSlpSU++OADucPT4JkL+XP+/Hm5QyAiKvFiYmLQs2dPXLhwAQAwadIkfPvttzJHVbboTOF0/vx59OnTBz179gQAVKxYEbt378alS5dkjoyIiKjohYeHw9PTU3ONzyVLlrz12p1UNErGMS4ttG7dGsePH8edO3cAAH///TfOnj2rOX2eiIiotLp16xZatmyJmzdvwsnJCXv27GHRJBOdGXGaNWsW4uLiUL16dahUKqjVaixZsgRDhgzJdf+0tDSkpaVp7sfHxwN4OcHs9cUCJUmCeHn5GRRkrnzWY8vIfPsShbmXR2HkPetzl9tnk3KX9TuL+Sp+cuQ+MDAQ3bt3x4sXL+Dm5oYTJ07Azc2tTP38izrveelXZwqnvXv34ocffsCuXbtQq1YtXL16FdOmTYOzszN8fHxy7L906VIsWrQoR3tUVBRSU7Ov4J2RkQFJkpCRkQE9vfylRAihWXSQZ3YVL+ZeHoWV96zPX3R0dI5V0Sl3kiQhLi4OQogSc3JMWVHcuT9x4gRGjhyJtLQ0GBgYwNfXF8bGxm+9zFJpVNR5T0hI0HpfnVmOwNXVFbNnz8akSZM0bYsXL8YPP/yAW7du5dg/txEnV1dXxMbG5jjVUK1W4+7du7C3t3/j9cm0kZGRwV/8MmHu5VEYeY+OjkZkZCQ8PDxyXPaGcidJEqKiomBnZ8fCqZgVZ+537dqFESNGQK1Ww8zMDEFBQYV+xQJdUdR5j4+Ph5WVVelajiA5OTlHslQq1RuH1wwNDWFoaJijXalU5uhHqVTCysoKUVFRUCgUMDExyfNf0Flr2qjVao56FDPmXh4FzbsQAsnJyYiKioKVlRUL3zzKurQSC6fiVxy5X716NaZNmwbg5TVFb9y4AXt7+yJ7Pl1QlHnPS586Uzj16tULS5YsQYUKFVCrVi1cuXIFK1euxPvvv18o/Ts6OgJAvoc/s469KpVKfnkXM+ZeHoWV93Llymk+f0RlnRACc+fOxdKlSwG8vGTWd999V2jX9KSC05nC6ZtvvsGnn36KiRMnIjIyEs7Ozhg3bhzmz59fKP0rFAo4OTnB3t4eGRkZeX581hwNGxsb/gVYzJh7eRRG3vX19Xl4juhfmZmZGD9+PDZv3gwA+OyzzzBv3jz+QVjC6EzhZG5ujlWrVmHVqlVF+jwqlSpfv8glSYK+vj6MjIz45V3MmHt5MO9EhSclJQWDBw/Gr7/+CgDw9PRk0VRC8bcdERGRjGJiYtCpUydN0eTj44MjR46waCqhWDgRERHJ5NGjR2jdurXmslOffPIJtm3bJm9Q9FY6c6iOiIioNLlx4wa6deuGJ0+ewNnZGWvXrkXfvn3lDov+AwsnIiKiYnb69Gl4eXkhISEB5cuXx/nz5+Hq6ip3WKQFHqojIiIqRj///DO6dOmChIQEGBgYYM+ePSyadAgLJyIiomKybt06DBgwAOnp6TA1NcXff/+N1q1byx0W5QELJyIioiImhMC8efM0lw2ztbXF3bt3Ub16dZkjo7xi4URERFSEMjMzMXr0aCxZsgQAMHPmTNy7dw9OTk4yR0b5wcnhRERERSQpKQmDBg3CoUOHAABff/215hp0pJtYOBERERWB58+fo2fPnggMDATwcjXwKVOmyBwVFRQP1RERERWyBw8eoGXLlpqiadSoUTh8+DCvzVgKsHAiIiIqRH///TdatGiBu3fvAgAWLlyITZs28RIqpQQP1RERERWSEydOoE+fPkhISECdOnWwefNmNGnSRO6wqBCxcCIiIioEP/74I9577z1kZGTAzc0Np0+fRrly5eQOiwoZD9UREREV0Jo1azBo0CBkZGTA2NgY+/fvZ9FUSnHEiYiIKJ+EEJgzZw6WL18OALCyssLly5dRsWJFeQOjIsPCiYiISAtqSSDgfjTuhcegSqIKDV0tMW7sGHz//fcAADc3NwQFBcHGxkbmSKkosXAiIiL6D343IrD4UAjCY1P+bQmDfloc/rlwFyqVCr6+vnjvvfdgbGwsa5xU9Fg4ERERvYXfjQhM2HkZnarbY/WgelDHhGPMRwsQYVUXdn3nYJBLAsaMGSp3mFRMODmciIjoDdSSwOJDIehU3R4bvBvDNDUKA/r2wrUTvyBq/2LYpT/D+WRHqCUhd6hUTFg4ERERvUFgWAzCY1MwsUMVXLhwHq1atcLjx48BAIsWLcR30/rjcWwKAsNiZI6UigsP1REREb1BZEIqAOBO4EmM9B6K1NRU6OvrY/369Xj//feRmJaZbT8q/TjiRERE9Ab25kYAAJ8P5iA1NRU9e/bE1atXMWLECADA7acJ2faj0o+FExERUS4kScLP3y1H5ounsGz+LqpWq479+/fD2tr63+0CvifvwdXaGE0rWcscLRUXHqojIiJ6TVpaGkaOHIndu3fDuGoL2PWdg0adO+PakwRYKdW4/CgW60/dx/FbkfAd1hAqJS/gW1awcCIiInrFixcv0K9fP5w8eRIAUF6Kwhf9qmP1qccYsP6CZj9Xa2P4DmsIz9pOMkVKcmDhRERE9K/Hjx+je/fuuHnzJgwMDNCwYUMcP34cJiYmeKdpFQTcf4574VGo4mKHZu62HGkqg1g4ERERAbh27Ro8PT0REREBZ2dnHDp0CHXq1IFKpQIAqJQKNHe3gbuZGvb2NlCyaCqTODmciIjKvOPHj6NVq1aIiIgAAOzYsQP169fXFE1EWVg4ERFRmbZjxw5069YNiYmJUCgUWLduHTp27Ch3WFRC8VAdERGVSUIILFu2DJ988gkAQE9PDwcPHkTPnj1ljoxKMhZORERU5mRmZmLKlClYv349AMDNzQ379+9Hw4YNZY6MSjoeqiMiojIlKSkJ77zzjqZoWrNmDR48eMCiibTCESciIiozIiMj0atXLwQGBgIAOnbsiClTpsgcFekSFk5ERFQm3L17F927d0doaCiAl0XToUOHZI6KdA0P1RERUal37tw5tGzZUlM0+fj4wN/fH0ZGvDgv5Q0LJyIiKtX27duHjh074vnz56hcuTIWLFiAbdu2QankVyDlHQ/VERFRqSSEwP/+9z/MnDkTAODl5YU9e/bA1NRU5shIl7FwIiKiUuf15QYMDQ2xcuVKFk1UYCyciIioVElMTMSgQYNw+PBhAEC5cuVw9uxZeHh4yBwZlQYsnIiIqNR48uQJvLy8cOXKFQAvF7Y8d+4cnJ2dZY6MSgvOjCMiolLh+vXraNasGa5cuQJbW1uMHTsW165dY9FEhYojTkREpPP8/f3xzjvvIDExEa6urjh58iTc3d3lDotKIY44ERGRTtu8eTO6d++OxMREAIC3tzeLJioyLJyIiEgnCSEwb948jB49Gmq1GgqFAmvXrsWSJUvkDo1KMR6qIyIinZOWlob3338fu3btAgAYGBjgl19+gaenp8yRUWnHwomIiHRKTEwM+vbtizNnzkBPTw/z589H7969Ua9ePblDozKAhRMREemM+/fvo0ePHrh9+zZUKhUOHjyInj17yh0WlSE6Ncfpn3/+wXvvvQcbGxuYmJigfv36CAoKkjssIiIqBhcuXEDz5s1x+/ZtAEDz5s3Rtm1bmaOiskZnRpxiY2PRqlUrdOjQAUeOHIG9vT1CQ0NRrlw5uUMjIqIi9vPPP2PYsGFIS0sDAAwbNgzbtm2Dnp7OfI1RKaEz77jly5fD1dUVW7du1bRVrFhRvoCIiKjICSGwfPlyzJkzR9O2aNEifPrpp1AoFDJGRmWVzhyq+/XXX9G4cWMMHDgQ9vb2aNCgATZu3Ch3WEREVETS09MxatQoTdE0depUXL58GfPnz2fRRLLRmRGn+/fvw9fXFzNmzMAnn3yCwMBATJ06FYaGhhg+fHiO/dPS0jRDugAQHx8PAJAkCZIkFXp8kiRBCFEkfdPbMffyYN7lUxZyHxMTg4EDB+LkyZMAgMGDB+Prr78GAFlfd1nIfUlU1HnPS786UzhJkoTGjRvjiy++AAA0aNAAN2/ehK+vb66F09KlS7Fo0aIc7VFRUUhNTS2S+OLi4iCEgFKpMwN5pQJzLw/mXT6lPff379/H8OHDERoaCoVCATs7O0yZMgWRkZFyh1bqc19SFXXeExIStN5XZwonJycn1KxZM1tbjRo18PPPP+e6/5w5czBjxgzN/fj4eLi6usLOzg4WFhaFHp8kSZoPOD9MxYu5lwfzLp/SnPvTp0+jf//+iImJgUKhQJ06dXDs2DHY2dnJHRqA0p37kqyo825kZKT1vjpTOLVq1UpzCmqWO3fuwM3NLdf9DQ0NYWhomKNdqVQW2ZtdoVAUaf/0Zsy9PJh3+ZTG3G/fvh1jxoxBRkYGzM3N0b59e+zduxfGxsZyh5ZNacy9LijKvOelT535qU+fPh0XLlzAF198gXv37mHXrl3YsGEDJk2aJHdoRERUAJIkYd68eRgxYgQyMjIwcOBA3L17FwcPHixxRROR1iNOv/76a54779KlS6G96Zs0aYIDBw5gzpw5+Oyzz1CpUiWsWrUKw4YNK5T+iYioaKklgcCwGEQmpMLe3AhNK1kjPS0VPj4+2LdvHwDAwsIC27Ztg4mJiczREuVO68Kpb9++eepYoVDg7t27cHd3z2tMb+Tl5QUvL69C64+IiIqH340ILD4UgvDYFE2bk4UBUs7txN+HXhZNRkZG2LNnD4smKtHyNMfp6dOnsLe312pfc3PzfAVERESli9+NCEzYeRmdqttjzZAGqOZgjiPn/sZHW45B1OoPk7t3YR4XCn9/f9SpU0fucIneSuvCycfHJ0+H3d57770iOXuNiIh0h1oSWHwoBJ2q22ODd2MolQr4+fnh/XffRUJCIioMWwwHzwk4M6sTyjs7yR0u0X/SenL41q1b8zSK5OvrC1tb23wFRUREpUNgWAzCY1MwsUMVKJUKrFu3Dj169EBCQgLatm2DbbOGQTK2wqMUA7lDJdKKzpxVR0REuicy4eWCw1VsTTBt2jRMmjQJQggAwLRp09Ckavls+xGVdPlexyk1NRXXrl1DZGRkjqXKe/fuXeDAiIhI99mbv1xYsK/PRJz4aQsAQKVSYfPmzejXrx+CHsZm24+opMtX4eTn54fhw4fj+fPnObYpFAqo1eoCB0ZERLrPDnFQJMfghlQeCoUSpqYmOHToENq2bQtJEvA9eQ+u1sZoWsla7lCJtJKvQ3WTJ0/GwIEDERERoblobtaNRRMREQHAmTNn0KJ5Mzw76gvjKk3QcNp32Ol3Dg2btUTQw1iM3XEJx29FYm6PGlApFXKHS6SVfI04RUZGYsaMGXBwcCjseIiIqBTYsmULxo8fj4yMDLhbWOCzHu7wvWCKqb89An57BABwtTaG77CG8KzNs+lId+SrcBowYABOnjyJypUrF3Y8RESkw9RqNWbOnImVK1dq2rp06YJh7WphcJuaOVYO50gT6Zp8FU7ffvstBg4ciDNnzqBOnTrQ19fPtn3q1KmFEhwREemO+Ph4DBkyBIcPH9a0LV++HB9//DEAQKVUoEVlG7nCIyoU+Sqcdu3ahaNHj8LY2BgnT56EQvH/fzEoFAoWTkREZcz9+/fRq1cvBAcHQ6VSQaVSYc+ePejXr5/coREVqnwVTvPmzcNnn32G2bNnQ6nkUlBERGXZqVOn0L9/f0RHR8PZ2Rnr16+Hs7MzGjVqJHdoRIUuX4VTeno6Bg0axKKJiKiM27x5M8aPH4/MzEyYmJjg2LFjqFWrltxhERWZfFU+Pj4+2Lt3b2HHQkREOkKtVmPGjBkYPXo0MjMzAQCdO3eGm5ubzJERFa18jTip1WqsWLECR48eRd26dXNMDn/1bAoiIipd4uLiMGTIEBw5ckTT9sknn+Dzzz/nkQgq9fJVOF2/fh0NGjQAANy4cSPbtlcnihMRUekSGhqKXr16ISQkBAYGBlAoFNiyZQuGDh0qd2hExSJfhdOJEycKOw4iIirhTpw4gYEDByI6OhpOTk747bff4OrqCnt7e7lDIyo2hX6RX4VCgV69ehVKcEREJD8hBNauXYtp06ZBrVZDqVRi6tSpPGuOyqR8X+TX29sb0dHRObbxIr9ERKVHWloaJk2ahM2bN2va+vbty/X6qMzK90V+3333XV7kl4ioFHv69Ck6duyYrWhauHAhfvrpJ5iYmMgYGZF8eJFfIiLK4dKlS+jXrx/Cw8NhYWGBhg0bYuLEiRg4cKDcoRHJihf5JSKibHbu3InRo0cjNTUVrq6u+OOPP1C1alW5wyIqEXiRXyIiAvByjb45c+bgyy+/BAAolUpUq1aNRRPRK3iRXyIiwosXLzB06NBsi1oOHjw42/wmIuJFfomIyrxbt26hT58+uHPnDlQqFSRJwpIlSzB79mwuakz0Gl7kl4ioDDt8+DCGDBmC+Ph4uLq64quvvoKRkRHX4yN6A17kl4ioDBJCYPny5fDy8kJ8fDzKlSuH06dPY+DAgSyaiN6CF/klIipjkpOTMXr0aOzevRvAy7mpTZs2hZWVlcyREZV8vMgvEVEZEhYWhnfeeQdXr16FQqGAEAIzZ87EkiVLoFKp5A6PqMTjRX6JiMoIf39/DB48GDExMbCwsEBaWhq2bduGwYMHyx0akc7I90V+iYhINwghsGLFCnzyySeQJAmNGzfG/v37oaenBycnJ7nDI9IpWk8Ov3btGiRJ0rrjmzdvIjMzM19BERFR4UhMTMS7776L2bNnQ5IkKJVKzJgxA66uriyaiPJB68KpQYMGiI6O1rrjFi1a4NGjR/kKioiICu7u3bto3rw5fvrpJ83yMePHj8eAAQNkjoxId2l9qE4IgU8//VTrK2Knp6fnOygiIiqY33//He+99x7i4uJgYGAASZLw3XffYfTo0XKHRqTTtC6c2rZti9u3b2vdcYsWLWBsbJyvoIiIKH8kScLixYuxYMECAEDz5s1hZ2eHWbNmoVWrVjJHR6T7tC6cTp48WYRhEBFRQcXFxWH48OH49ddfAQD9+/fHrl27YGBgIHNkRKUHz6ojIioFgoOD0a9fP8315tRqNSpXrsyiiaiQsXAiItIhakkg4H407oXHoEqiCs3cbfHLwQPw8fFBYmKiZj7TunXrMH78eLnDJSp1WDgREekIvxsRWHwoBOGxKf+2hMEUqQg78BVSEhNhZmYGU1NT7N+/Hy1btpQ1VqLSKl+F0+3bt1GtWrXCjoWIiN7A70YEJuy8jE7V7bF6UD2IF08wfeGXuK2qCLu+c9BCfRMz3u2IChUqwNHRUe5wiUotrddxelXdunXRo0cPHDt2rLDjISKi16glgcWHQtCpuj02eDdGesQd9OnpibO/7EDML0uh9+wW/rFrikaNm7BoIipi+SqcwsLC0LRpU4wcORI1atTA+vXrkZycXNixERERgMCwGITHpmBi+8pYv94X7dq1w5MnT2BoaAh9fT2MaumCx7EpCAyLkTtUolIvX4WTs7MzFi5ciIcPH2LRokXYs2cPXFxcMHPmTDx8+LCwYyQiKtMiE1IBAF/Om4FJkyYhIyMD+vr6cHR0xPnz5zHVZ2C2/Yio6ORrjlNKSgpiY2MRGxsLZ2dnzJgxA2fPnoWvry9WrVrFVcOJiApR2otIAMD+Py9ApVLBy8sLcXFx2LdvH2xtbRH0MBYAYG9uJGeYRGVCvgonU1NTWFhYwM7ODubm5rCwsICFhQX69OkDCwuLwo6RiKjM+vnnnzH6/fdhMWQlbNsMxc4J7dCmdWs8e/YM1tbWkCQB35P34GptjKaVrOUOl6jUy1fhNHDgQBw7dgyenp744IMPUKVKlcKOi4ioTMvIyMCcOXPw1VdfAQBUf+2ARY8Psfm2CmZuL2ClFLj8KBbrT93H8VuR8B3WECqlQuaoiUq/fM1x2rt3L65fvw5TU1M0b94cvXv3xokTJwo7NiKiMikiIgKdOnXSFE36+vpwlqKwoFN5hMWmY8D6C+i07ioGrL+A288S4DusITxrO8kcNVHZkO8FMF1cXLBs2TLMnz8f27dvx4QJE2BoaIhp06Zh5MiRhRkjEVGZcerUKQwaNAjPnj2DkZERUlNT0bdvX2zZsgVmZmbw6SwQcP857oVHoYqLHZq523KkiagY5atwWr16NRISEpCYmKj5t3r16vjzzz8xevToIi+cli5dik8++QQffPABVq1aVaTPRURUHIQQ+N///oc5c+ZArVajdu3a2LFjB65evQofHx8oFC+LI5VSgebuNnA3U8Pe3gZKFk1ExSpfhdOePXtQrlw5WFlZoVy5cnBxcUHt2rXh5eUFKyurwo4xm4sXL2LDhg2oW7dukT4PEVFxiYuLw4gRI3Dw4EEAgJmZGb7//nvUr18f9evXlzU2IsouX4XT+fPnCzsOrSQmJmLYsGHYuHEjFi9eLEsMRESF6cqVK3j33Xdx7949qFQqSJKEFi1awMXFRe7QiCgX+ZocLpdJkyahZ8+e6Ny5s9yhEBEViBACvr6+aN68Oe7duwcjIyNIkoT58+fjyJEjsLOzkztEIspFvieHF7c9e/bg8uXLuHjxolb7p6WlIS0tTXM/Pj4eACBJEiRJKvT4JEmCEKJI+qa3Y+7lwbznX3x8PMaNG4cff/wRANChQweEhYVh7dq18PT0BIC35pW5lw9zL4+iznte+tWJwunx48f44IMPcOzYMRgZabcy7tKlS7Fo0aIc7VFRUUhNLfzLEkiShLi4OAghoFTq1ECezmPu5cG858+NGzcwbtw43L9/H0qlEjNnzsTUqVOhVquhp6eHyMjI/+yDuZcPcy+Pos57QkKC1vsqhBCi0CMoZAcPHkS/fv2gUqk0bWq1GgqFAkqlEmlpadm2AbmPOLm6uiI2NrZIVjeXJAlRUVGws7Pjh6mYMffyYN7zRgiBDRs2YPr06UhLS4OJiQmSk5OxZcsW+Pj45Kkv5l4+zL08ijrv8fHxsLKyQlxc3H/WCDox4tSpUydcv349W9vIkSNRvXp1zJo1K0fRBACGhoYwNDTM0a5UKovszZ5VyPHDVPyYe3kw79pJSEjAuHHjsHv3bgAvz5oTQmD37t0YPHhwvvpk7uXD3MujKPOelz51onAyNzdH7dq1s7WZmprCxsYmRzsRUUly7do1DBw4EHfu3IFSqYSenh4qVKiAn376CTVq1JA7PCLKI50onIiIdI0QAps2bcLUqVORmpoKFxcX7N27FxEREfD09ISpqancIRJRPuhs4XTy5Em5QyAiylViYiLGjx+PnTt3Anh5aG7cuHFo2bKlzJERUUHpbOFERFQSXblyBYMHD8adO3egUCigp6eHihUrYsCAAXKHRkSFgDPbiIgKgRACq1evRvPmzXHnzh2YmJhACAEfHx8EBASgevXqcodIRIWAI05ERAX0/PlzjBw5Er///jsAoE+fPpAkCYMGDcKwYcNkjo6IChMLJyKiAjh58iSGDRuGJ0+eQE9PDxMmTMDq1auhUCjkDo2IigAP1RER5UNmZibmz5+Pjh074smTJ7CwsEBmZibMzc1ZNBGVYhxxIiLKo0ePHmHYsGE4e/YsAMDS0hJqtbpAC1oSkW5g4URElAcHDhzAqFGjEBsbCzMzM1hYWMDBwQF79+6Fh4eH3OERURFj4UREpIWUlBR89NFHWLduHQCgXr162L9/P4QQcHFxyfUST0RU+rBwIiL6D8HBwRg8eLDmmpmmpqaoUKEC3N3dZY6MiIobJ4cTEb2BEAJr165Fo0aNcP36dZiYmAAAWrZsiQ0bNsgcHRHJgYUTEVEunj17Bi8vL0yePBmpqamwsrJCeno6vvzyS/j5+cHR0VHuEIlIBjxUR0T0mkOHDmHkyJGIioqCgYEBvvzySzg4OMDd3R1NmjSROzwikhELJyKif6WkpODjjz/G2rVrAQDlypVDly5dMHXqVJkjI6KSgoUTEZU5akkgMCwGkQmpsDc3QtNK1rh+7W8MHToUISEhAF6uzZSZmYnevXvLHC0RlSQsnIioTPG7EYHFh0IQHpuiaTNXpuPhwa+REBICMzMzpKSkwMPDA7t370aVKlVkjJaIShoWTkRUZvjdiMCEnZfRqbo91gxpAHN1AkZO+wR39Nxh1etjNGzYEFWMEmFvb4+FCxfCwMBA7pCJqIRh4UREZYJaElh8KASdqttjg3djHDx4AGPGjEFMTAwMDAzhOuEbqOv0xvqP2kNPxROOiSh3/O1ARGVCYFgMwmNT4NPEEe+/PxL9+/dHTEzMv8sMpKFKeigex6bg4oNYuUMlohKMI05EVCZEJqQCAIZ5dcTjsHtQKBQwMzODEAJ79uxBz779UXvBUc1+RES5YeFERKVeSkoKfvx+E2DQGM/SVHB3d4eFhQWsra2xfft2uLi4IOjhy5Eme3MjmaMlopKMhRMRlWoXL17E8OHDcev2HZQfuwHV+03FmaXeUGdmwtLSEkqlEpIk4HvyHlytjdG0krXcIRNRCcY5TkRUKmVkZGDhwoVo0aIFbt26BXMzU8Sf3o440wqYsf8W7scDyRkSgh7GYuyOSzh+KxJze9SASqmQO3QiKsE44kREpU5wcDCGDx+OoKAgAICDgwOePXuG6Z6N0X5QXazwv4f+vuc0+7taG8N3WEN41naSK2Qi0hEsnIio1JAkCatXr8acOXOQlpaGcuXKwczMDADg7++Pzp07AwB61nfNsXI4R5qISBssnIioVHjw4AFGjhyJkydPAgA6dOiAH374AVFRUahYsSIsLS01+6qUCrSobCNTpESky1g4EZFOE0Jg06ZN+PDDD5GQkABDQ0Po6enB0dERzs7OcHZ2ljtEIipFWDgRkc56+PAhRo8ejT/++AMAYGNjg+joaPTu3Rtr1qyROToiKo1YOBGRzpEkCRs2bMDHH3+MxMREGBkZwcTERLOY5aBBg+QOkYhKKRZORKRTwsLCMHr0aPz5558AgBYtWmD79u0IDg5G06ZN4eTEM+OIqOiwcCIinSBJEnx9fTFr1iwkJSXB0NAQxsbGaNu2LTw8PODh4SF3iERUBnABTCIq8UJDQ9GxY0dMnjwZSUlJcHJyQlpaGho1aoSJEyfKHR4RlSEsnIioxJIkCd988w3q1q2LU6dOwdjYGFZWVoiPj8e3336LY8eOoUKFCnKHSURlCA/VEVGJdO/ePbz//vs4c+YMAKBdu3bYsmULDhw4gP79+6NixYryBkhEZRILJyIqUTIzM7Fy5UosWLAAqampMDQ0hKGhId5//324u7vjww8/lDtEIirDWDgRUYkRFBSE0aNH4+rVqwAAe3t7REZGolu3bprLpRARyYmFExHJLikpCQsWLMDXX38NSZJgYWEBSZIgSRJ2796NQYMGQaHgteSISH6cHE5Esjp27Bjq1KmDr776CpIkYfDgwbh58yamTJmC4OBgDB48mEUTEZUYHHEiIlk8f/4cM2bMwI4dOwAAVlZWUKvV+Oijj+Di4oIvvvhC5giJiHLiiBMRFSshBHbu3IkaNWpoiiYXFxfExsZiwIABcHd3lzlCIqI344gTERWbBw8eYMKECfDz8wMAVKhQAU+fPoWBgQGOHTuGLl26yBwhEdHbccSJiIpc1hIDtWrVgp+fH/T19fH555/j3LlzmD17Nq5fv86iiYh0AkeciKhInTt3DhMmTMC1a9cAAI6OjkhLS8OECRNgY2ODRYsWyRwhEZH2OOJEREUiOjoaY8aMQatWrXDt2jWYmprC3NwcycnJ+OKLL2BlZSV3iEREecbCiYgKlSRJ2LJlC6pVq4ZNmzYBAGrWrImkpCR4enri1q1bGD9+PJRK/vohIt3DQ3VElGdqSSDgfjTuhcegSqIKzdxtoVIqcP36dUyYMAF//fUXAKBatWrYuHEjzMzM8OTJE/Ts2VPmyImICoaFExHlid+NCCw+FILw2JR/W8JQ3tIIFWIuYd/Xc6FWq2FkZIRy5crBzMwMrVu3hkKhQIMGDWSNm4ioMHCsnIi05ncjAhN2XkZ1R3P8PL45jk+ohynVUxF+MwDnlDVhULkpKlSogNTUVLi4uGDDhg1c9ZuIShUWTkSkFbUksPhQCDpVt8cG78awzIzFuFEj8NHIAXi0cx4UT27ArvNYJCQmYv369bhw4QIaNmwod9hERIVKZwqnpUuXokmTJjA3N4e9vT369u2L27dvyx0WUZkRGBaD8NgUvN+8PBYsmI/atWvj+PHjUKlU+OCDqdg1bwQU5nbY5R+IcePGQaVSyR0yEVGh05nC6dSpU5g0aRIuXLgAf39/ZGZmomvXrkhKSpI7NKIyITI+FQAwuHs7LF68GGlpaXB0dIQkSTAzM0PdinYAgDSlkZxhEhEVKZ2ZHJ51iYYsW7duhb29PYKCgtC2bVuZoiIqG65fv46lC74Eag7CszQVrK2tkZGRgeTkZKxevRoTJ07E1fB4AIC9OQsnIiq9dGbE6XVxcXEAAGtra5kjISq9YmNjMXXqVDRo0AABv++EOu4Z6g+Zifj4BPTt2xdnz57FpEmToFAo4XvyHlytjdG0Ej+TRFR66cyI06uEEJgxYwZat26N2rVr57pPWloa0tLSNPfj41/+NSxJEiRJKvSYJEmCEKJI+qa3Y+4Ln1qtxpYtWzBv3jw8f/4cAFC1alXM6lcfn5+IwIDV/pjsWQvGqlRcehCD706H4c/bkVg7pAEUEJAkIfMrKN34npcPcy+Pos57XvrVycJp8uTJuHbtGs6ePfvGfZYuXZrrNbCioqKQmppa6DFJkoS4uDgIIbgicjFj7gvXhQsXsGDBAs215WxtbZGYmIjnz5/D3SAeX/R0x5rT4Xh3Q4DmMc4WBviipzsa2isRGRkpV+hlBt/z8mHu5VHUeU9ISNB6X4UQQqf+NJwyZQoOHjyI06dPo1KlSm/cL7cRJ1dXV8TGxsLCwqLQ45IkCVFRUbCzs+OHqZgx94Xj/v37mDVrFvbv3w8AMDExgbGxMeLi4jBp0iR8+umnmuvLZa0cHvokCpWd7dDM3QYqJddrKi58z8uHuZdHUec9Pj4eVlZWiIuL+88aQWdGnIQQmDJlCg4cOICTJ0++tWgCAENDQxgaGuZoVyqVRfZmVygURdo/vRlzn39xcXFYvHgx1qxZg/T0dCgUCowdOxb16tXDwYMHsWrVKtSoUSPbY5RKoGUVW1SxkGBvb8u8y4Dvefkw9/IoyrznpU+dKZwmTZqEXbt24ZdffoG5uTmePn0KALC0tISxsbHM0RHpnszMTGzcuBHz58/XzGOqUKECateujfXr10MIgfHjx3PlbyKiV+hM4eTr6wsAaN++fbb2rVu3YsSIEcUfEJEO8/Pzw4cffojg4GAAgL29PZKTkxEdHY02bdpACMGCiYgoFzpTOOnYVCyiEunmzZv46KOPNOuiWVtbw9DQEM+ePcPo0aOxaNEiODo6yhwlEVHJpTOFExHlX3h4OObPn4/t27dDkiSoVCpMmjQJCxcuxIEDB9C0adM3Lu1BRET/j4UTUSn24sULLFu2DKtXr9Ysw+Hq6orHjx+jZs2asLKywvvvvy9zlEREuoOFE1EplJaWhrVr12LJkiWIiYkBALi4uCAyMhIZGRnYtGkT5wYSEeUDz6UkKkUkScKOHTtQrVo1fPjhh4iJiUGtWrXw6aefIjo6GrNnz8bdu3cxatQoqFQqucMlItI5HHEiKgWEEDh27BhmzZqFv//+GwBQrlw51KlTBydOnIBSqcSkSZPg4OAgc6RERLqNI05EOu7MmTNo164dPD098ffff8PExAT29vZ48eIFHB0dkZGRAYVCwaKJiKgQsHAi0lEXL15Et27d0LZtW5w5cwYGBgaoUKECkpOTUadOHVy8eBE//vgjjIyM5A6ViKjU4KE6Ih1z/fp1zJ8/HwcPHgQAqFQq+Pj4YNGiRfDz84Obmxu6dOkib5BERKUUCyciHXH37l0sWLAAe/bs0azs7ebmhgcPHqBGjRpwcXHB6NGj5Q6TiKhU46E6ohLuwYMHGDVqFGrUqIHdu3dDCIEKFSoAeHkW3YYNG/DBBx/IHCURUdnAwolIJmpJ4HxoNH65+g/Oh0ZDLWW/rFBoaChGjRoFDw8PbNmyBWq1Gl5eXti4cSPUajXWrVuHu3fvYsyYMdDX15fpVRARlS08VEckA78bEVh8KAThsSmaNhcrY8zrWQOVDROxZMkS/PDDD1Cr1QAAZ2dnODo64tdff4VCocDw4cNhYGAgV/hERGUWCyeiYuZ3IwITdl5Gp+r2WDOkAao5mOP2swQs/+0qxu8IwvNfliHp9l8AXl4eJSIiAhkZGRg2bJjmOnMsmoiI5MHCiagYqSWBxYdC0Km6PTZ4N4ZSqUBwcDC+/Pxz7Nv7I2z7zYVl+5FoW7kcAgMuIDk5GUuWLMGkSZNgamoqd/hERGUeCyeiYhQYFoPw2BSsGdIAV65cxrJly/Dzzz9DiJfzm4zCTkHpMROffrwd4tlt1KlTB+bm5jJHTUREWVg4ERWjyPhUAMDMse/h+NHDmnZnZ2dERERAxP7zcr+EVPRp2VKWGImI6M14Vh1RMZAkCfv378eCWdMAAGeuh0KlUqFevXoAAGNjY2zYsAE/+Z8DANibc7VvIqKSiCNOREUoLS0NO3fuxIoVK3D79m1AoUT5yt1Rvd9UHPiwO17ExiIkJASDBg2CUqnC2B2X4GptjKaVrOUOnYiIcsERJ6IiEB8fj5UrV6Jy5coYNWoUbt++/fLiu3a2iPlzE+JMK+CL088hbCqhz4BB+PufBIzdcQnHb0Vibo8aUCkVcr8EIiLKBUeciArRgwcPsGbNGmzevBnx8fEAADs7O6SnpyMuLg7t27fHxx9/jBSbalhyOAT9fc9pHutqbQzfYQ3hWdtJrvCJiOg/sHAiKiAhBM6fP4+vv/4a+/fvhyRJAF5O+P7888/xzjvvYO7cuZg4cSJq1aqleVzXWo4IDItBZEIq7M2N0LSSNUeaiIhKOBZORPmUkZGBn3/+GV9//TUCAwM17U5OToiKikJmZiaGDh0KIyMjrF27NsfjVUoFWlS2Kc6QiYiogDjHiSiPXrx4gS+//BKVK1fGkCFDEBgYCAMDAzg7OwMAVCoVli1bhrt378LIiGfHERGVJhxxItLS1atXsW7dOuzcuRPJyckAADMzM0yfPh2TJ0/Gt99+izp16qBv37686C4RUSnFwonoLdLS0vDTTz9h7dq1OH/+vKbd2toaiYmJSE9PR69evWBvb4/PPvtMxkiJiKg4sHAiysXDhw+xfv16bN68GVFRUQAAPT09uLm5ITQ0FIaGhpg+fTrGjh0Le3t7maMlIqLiwsKJ6F+SJMHf3x9r167FoUOHNGfHmZubY/jw4Zg3bx4CAwORnJyM/v3783AcEVEZxMKJyrxHjx5h69at2Lp1Kx4+fKhpt7W1RUxMDIQQaNmyJRwdHdG7d28ZIyUiIrmxcKIyKS0tDb/++is2bdoEf39/CCEAAOXKlUPVqlURGBiIChUqYMmSJRgyZAjMzc1ljpiIiEoCFk5Upty4cQObN2/Gjh07EB0drWm3tbVF165dsWnTJjx//hyRkZFo1KiRjJESEVFJxMKJdJpaEgi4H4174TGokqhCM3fbHKtvx8TEYN++fdi6dSsCAgI07aamphBCIDk5GS4uLujTpw+MjY3h6uoKV1fX4n4pRESkA1g4kc7yuxGBxYdCEB6b8m9LGFysjDGvZw20r2KFw4cPY8eOHTh06BAyMjIAvFycsnfv3mjQoAG++eYbvPfee/Dx8UG9evXkeyFERKQzWDiRTvK7EYEJOy+jU3V7rB5UD1bKVESrDfHFgUsYvyMISf6r8fzKH5r9LSwskJiYiLZt22quJzd79myeGUdERHnCS66QzlFLAosPhaBTdXts8G4Mg4QnWPW/5RjQvhEOftgDyfcCYdB0MGzt7DSXPKlatSpWr16Nffv2AQCUSiWLJiIiyjOOOJHOCQyLRnhsCuomBaFOnREIDg7WbDMwMIBt1GUkezTDNv9A/LnnOwwbNgy1a9eWMWIiIiotWDiRThBC4MaNG9i3bx92n7sLNH4P65fNh8hIhVKphIGBAVJTU2FtbY0+7RpjdwaQLPSxdOlSuUMnIqJShIUTlVhCCFy7dg379u3Dvn37cOfOHQCAoWsdODYGOvQbhh5Na2DZsmXo0aMHRowYgbZt2+JqeDx2+56DvbmRzK+AiIhKGxZOVKKkp6fj9OnT+O233/Dbb78hLCxMs00zshR+E8qUWFTqMQ7TvRtj6tSpiI6O/veacQr4nrwHV2tjNK1kLd8LISKiUomFE8kuOjoahw8fxm+//YajR48iPj5es83Q0BDNmzfHqVOn4OLign79+qFv375IsvLA5D1XMXbHJYxv5w4rpRqXH8Vi/an7OH4rEr7DGuZYz4mIiKigWDhRsRNC4NatW5pRpXPnzmkuqAu8nOCtVqtRuXJlBAUFwcTEBNeuXUO9evWgUPx/MaRSKbH4UAgGrL+gaXO1NobvsIbwrO1UrK+JiIjKBhZOVCyioqJw/PhxHDt2DP7+/ggPD8+23d3dHffv34dKpULLli3h5eWFfv36wczMDABQv379HH161nZCl5qOCLj/HPfCo1DFxS7XlcOJiIgKCwsnKhKpqak4e/Ys/P394e/vjytXrmTbrlAooFAoULt2bfz666+ws7ODv78/OnbsmKcL6qqUCjR3t4G7mRr29jZQsmgiIqIixMKJCkVaWhouXbqEM2fO4MSJEzh9+jRSU1Oz7VOhQgU8evQI+vr6aNOmDby8vNC9e3e4ubkBAPr06SNH6ERERFpj4VTGqSWBwLAYRCakwt7cCE0rWWt1qCshIQHnz5/HmTNncObMGQQEBOQolJRKJSRJQocOHbB7926oVCoEBgaiXbt2MDU1LaqXREREVGRYOJVhOS+SC81Fcl+fXP3PP/8gICAAZ8+exZkzZ3DlyhWo1eps+9jZ2cHBwQE3b95EjRo14OnpiY4dO6J169awtLQEAPTo0aPoXxgREVERYeFURr16kdw1QxqgmoM5bj9LwLoT9zDhh8uYWM8A6odBCAwMREBAAJ48eZKjD4VCASEElEolPv30UyxYsABxcXEAgHLlyhXzKyIiIip6LJzKoFcvkvt5N1dc+zsQf+z5G1evXsXlK1eRVGsAVr1ww5MNcwHx/8sEVK1aFR07dkRMTAzi4uLQsWNHtGjRAo0bN4axsTEAFkxERFS6sXAqIxITE3Hr1i2EhITgz5vhCEdd3P/hU2wecS7Hvgbx++Dk/T8YutSC8vk91KhRA23btsWUKVPg7u4uQ/REREQlg84VTuvWrcOXX36JiIgI1KpVC6tWrUKbNm1kjUktCQTcj8a98BhUSVTJtpaQJEn4559/EBoairt37yIkJAQhISEIDg7Go0ePNPuZ1GgLu9518fT25WyPt7CwwNatW9GgaQt0+PYyvtm4HaO71s+26CQREVFZplOF0969ezFt2jSsW7cOrVq1wnfffYfu3bsjODgYFSpUkCWmnBOsw944wbqgJEnC8+fP8fjxY4SHh+Phw4cIDQ1FaGgo7t27h7CwMKSnp7+1DwsLC3Tr2AaBALoM9IFn42po1KgRatWqBRsbGwBA0MNYAEDtKhVYNBEREb1CpwqnlStXYtSoURg9ejQAYNWqVTh69Ch8fX2xdOnSYo/n1QnWqwfVg5UyFbGSEXxP3ceEnZe1vvRHSkoKoqKiEBkZmeP29OlTPHr0CA8fPsSTJ0+QmZn5n/1ZW1ujcePGsLS0RGhoKGrWrImmTZuiUaNGqFq1KqysbdDuyxNwrTEKH3g3zrZopCQJXiSXiIjoDXSmcEpPT0dQUBBmz56drb1r1644dy7nPJ2i9uoE6zUDa+OXXw7i2bNn0NPTQ9XkFNxWOOPDHX/hDxGIzPR0REVF4fnz53jx4gXi4+MRHx+PtLQ0JCUl5Vj/SBu1a9eGh4cHEhISkJmZCQ8PD9StWxeNGjVCjRo1/nOS9ryeNTBh52WM3XEJE9pXQTVHc9x+mgDfk/d4kVwiIqI30JnC6fnz51Cr1XBwcMjW7uDggKdPn+bYPy0tDWlpaZr78fHxAF4e7nr1grL5FXA/GuGxKVg9qB6SkhIxdOjQbNsNnKvDyft/+Ha3H9IeX9eqTz09PRgaGsLY2Bj29vbw9PSEg4MDwsLCULFiRVSvXh21a9eGi4sL9PX139rXf73GrjUdsHZIA3xx5Bb6+/5/4elqZYy1Qxqga02HQslTcZAkCUIInYm3tGDe5cPcy4e5l0dR5z0v/epM4ZTl9Tk3Qohc5+EsXboUixYtytEeFRWVrxGe190LjwEAWClTkZiYjubNmyM6OhpGRkYwMDCAnrEZngLo5NUPlfXaIykpCXp6eihXrhysra1hbW0NZ2dn2Nvbw8LCAhYWFlrPJ4qNjS1w/ADQ0F6Jvd41cPWfREQnZcDGVB/1y5tBpVQgMjKyUJ6jOEiShLi4OM2aUlQ8mHf5MPfyYe7lUdR5T0hI0HpfnSmcbG1toVKpcowuRUZG5hiFAoA5c+ZgxowZmvvx8fFwdXWFnZ0dLCwsChxPlUQVgDDESkZo4OaEM2fOICoqCnZ2dlAqlbj8KBYD1l/A3BmT0dzdpsDPV5ScHHPmT5dIkgSFQqHJPRUP5l0+zL18mHt5FHXejYyMtN5XZwonAwMDNGrUCP7+/ujXr5+m3d/fP9eLwxoaGsLQ0DBHu1KpLJSkN3O3hYuVMXxP3ccG78ZQKpVQKBT/9q3A+lP34WptjGbuttkmX1PRyMo9f5EVL+ZdPsy9fJh7eRRl3vPSp0791GfMmIFNmzZhy5YtCAkJwfTp0/Ho0SOMHz++2GNRKRWY17MGjt+KxNgdl3D5USyS0tW4/CgWY3dcwvFbkZjbowYnWBMREZUiOjPiBACDBg1CdHQ0PvvsM0RERKB27do4fPgw3NzcZInHs7YTfIc1xOJDIRiw/oKm3dXaWOulCIiIiEh36FThBAATJ07ExIkT5Q5Dw7O2E7rUdETA/ee4Fx6FKi52sq0cTkREREVL5wqnkkilVKC5uw3czdSwt7fhnCYiIqJSSqfmOBERERHJiYUTERERkZZYOBERERFpiYUTERERkZZYOBERERFpiYUTERERkZbKzHIEQggAL69ZVxQkSUJCQgKMjIy4DH8xY+7lwbzLh7mXD3Mvj6LOe1ZtkFUrvE2ZKZyyrnzs6uoqcyRERERUEiUkJMDS0vKt+yiENuVVKSBJEp48eQJzc3MoFIW/QGV8fDxcXV3x+PFjWFhYFHr/9GbMvTyYd/kw9/Jh7uVR1HkXQiAhIQHOzs7/OaJVZkaclEolXFxcivx5LCws+GGSCXMvD+ZdPsy9fJh7eRRl3v9rpCkLD9ASERERaYmFExEREZGWWDgVEkNDQyxYsACGhoZyh1LmMPfyYN7lw9zLh7mXR0nKe5mZHE5ERERUUBxxIiIiItISCyciIiIiLbFwIiIiItISCyciIiIiLbFwyoN169ahUqVKMDIyQqNGjXDmzJm37n/q1Ck0atQIRkZGcHd3x/r164sp0tIlL3nfv38/unTpAjs7O1hYWKBFixY4evRoMUZbuuT1PZ/lr7/+gp6eHurXr1+0AZZiec19Wloa5s6dCzc3NxgaGqJy5crYsmVLMUVbeuQ17zt37kS9evVgYmICJycnjBw5EtHR0cUUbelx+vRp9OrVC87OzlAoFDh48OB/Pka271hBWtmzZ4/Q19cXGzduFMHBweKDDz4Qpqam4uHDh7nuf//+fWFiYiI++OADERwcLDZu3Cj09fXFTz/9VMyR67a85v2DDz4Qy5cvF4GBgeLOnTtizpw5Ql9fX1y+fLmYI9d9ec19lhcvXgh3d3fRtWtXUa9eveIJtpTJT+579+4tmjVrJvz9/UVYWJgICAgQf/31VzFGrfvymvczZ84IpVIpVq9eLe7fvy/OnDkjatWqJfr27VvMkeu+w4cPi7lz54qff/5ZABAHDhx46/5yfseycNJS06ZNxfjx47O1Va9eXcyePTvX/WfOnCmqV6+erW3cuHGiefPmRRZjaZTXvOemZs2aYtGiRYUdWqmX39wPGjRIzJs3TyxYsICFUz7lNfdHjhwRlpaWIjo6ujjCK7Xymvcvv/xSuLu7Z2tbs2aNcHFxKbIYywJtCic5v2N5qE4L6enpCAoKQteuXbO1d+3aFefOncv1MefPn8+xf7du3XDp0iVkZGQUWaylSX7y/jpJkpCQkABra+uiCLHUym/ut27ditDQUCxYsKCoQyy18pP7X3/9FY0bN8aKFStQvnx5VK1aFR999BFSUlKKI+RSIT95b9myJcLDw3H48GEIIfDs2TP89NNP6NmzZ3GEXKbJ+R1bZi7yWxDPnz+HWq2Gg4NDtnYHBwc8ffo018c8ffo01/0zMzPx/PlzODk5FVm8pUV+8v66r776CklJSXj33XeLIsRSKz+5v3v3LmbPno0zZ85AT4+/WvIrP7m/f/8+zp49CyMjIxw4cADPnz/HxIkTERMTw3lOWspP3lu2bImdO3di0KBBSE1NRWZmJnr37o1vvvmmOEIu0+T8juWIUx4oFIps94UQOdr+a//c2unt8pr3LLt378bChQuxd+9e2NvbF1V4pZq2uVer1Rg6dCgWLVqEqlWrFld4pVpe3veSJEGhUGDnzp1o2rQpevTogZUrV2Lbtm0cdcqjvOQ9ODgYU6dOxfz58xEUFAQ/Pz+EhYVh/PjxxRFqmSfXdyz/LNSCra0tVCpVjr86IiMjc1S8WRwdHXPdX09PDzY2NkUWa2mSn7xn2bt3L0aNGoV9+/ahc+fORRlmqZTX3CckJODSpUu4cuUKJk+eDODll7kQAnp6ejh27Bg6duxYLLHruvy8752cnFC+fHlYWlpq2mrUqAEhBMLDw+Hh4VGkMZcG+cn70qVL0apVK3z88ccAgLp168LU1BRt2rTB4sWLeWShCMn5HcsRJy0YGBigUaNG8Pf3z9bu7++Pli1b5vqYFi1a5Nj/2LFjaNy4MfT19Yss1tIkP3kHXo40jRgxArt27eJcg3zKa+4tLCxw/fp1XL16VXMbP348qlWrhqtXr6JZs2bFFbrOy8/7vlWrVnjy5AkSExM1bXfu3IFSqYSLi0uRxlta5CfvycnJUCqzf42qVCoA/z/6QUVD1u/YIp9+Xkpknaa6efNmERwcLKZNmyZMTU3FgwcPhBBCzJ49W3h7e2v2zzpVcvr06SI4OFhs3ryZyxHkQ17zvmvXLqGnpyfWrl0rIiIiNLcXL17I9RJ0Vl5z/zqeVZd/ec19QkKCcHFxEQMGDBA3b94Up06dEh4eHmL06NFyvQSdlNe8b926Vejp6Yl169aJ0NBQcfbsWdG4cWPRtGlTuV6CzkpISBBXrlwRV65cEQDEypUrxZUrVzRLQZSk71gWTnmwdu1a4ebmJgwMDETDhg3FqVOnNNt8fHxEu3btsu1/8uRJ0aBBA2FgYCAqVqwofH19izni0iEveW/Xrp0AkOPm4+NT/IGXAnl9z7+KhVPB5DX3ISEhonPnzsLY2Fi4uLiIGTNmiOTk5GKOWvflNe9r1qwRNWvWFMbGxsLJyUkMGzZMhIeHF3PUuu/EiRNv/d1dkr5jFUJwPJGIiIhIG5zjRERERKQlFk5EREREWmLhRERERKQlFk5EREREWmLhRERERKQlFk5EREREWmLhRERERKQlFk5EREREWmLhRERERKQlFk5EVCa0b98eCoUCCoUCV69eLVBfI0aM0PR18ODBQomPiHQDCyciKjPGjBmDiIgI1K5du0D9rF69GhEREYUUFRHpEj25AyAiKi4mJiZwdHQscD+WlpawtLQshIiISNdwxImIdNLu3bthZGSEf/75R9M2evRo1K1bF3FxcVr30759e0yZMgXTpk2DlZUVHBwcsGHDBiQlJWHkyJEwNzdH5cqVceTIkaJ4GUSkY1g4EZFOGjx4MKpVq4alS5cCABYtWoSjR4/iyJEjeR4N2r59O2xtbREYGIgpU6ZgwoQJGDhwIFq2bInLly+jW7du8Pb2RnJyclG8FCLSISyciEgnKRQKLFmyBJs2bcIXX3yB1atXw8/PD+XLl89zX/Xq1cO8efPg4eGBOXPmwNjYGLa2thgzZgw8PDwwf/58REdH49q1a0XwSohIl3COExHpLC8vL9SsWROLFi3CsWPHUKtWrXz1U7duXc3/VSoVbGxsUKdOHU2bg4MDACAyMrJgARORzuOIExHprKNHj+LWrVtQq9Wa4iY/9PX1s91XKBTZ2hQKBQBAkqR8PwcRlQ4snIhIJ12+fBkDBw7Ed999h27duuHTTz+VOyQiKgN4qI6IdM6DBw/Qs2dPzJ49G97e3qhZsyaaNGmCoKAgNGrUSO7wiKgU44gTEemUmJgYdO/eHb1798Ynn3wCAGjUqBF69eqFuXPnyhwdEZV2HHEiIp1ibW2NkJCQHO2//PJLvvo7efJkjrYHDx7kaBNC5Kt/IipdOOJERGXGunXrYGZmhuvXrxeon/Hjx8PMzKyQoiIiXaIQ/DOKiMqAf/75BykpKQCAChUqwMDAIN99RUZGIj4+HgDg5OQEU1PTQomRiEo+Fk5EREREWuKhOiIiIiItsXAiIiIi0hILJyIiIiItsXAiIiIi0hILJyIiIiItsXAiIiIi0hILJyIiIiItsXAiIiIi0hILJyIiIiItsXAiIiIi0tL/ASr0SzgMUbOrAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "largest difference Timoshenko - Euler-Bernoulli: 0.058 mm at x = 0.42 m\n"
     ]
    }
   ],
   "source": [
    "w_ref, theta_ref = timoshenko_cantilever(nodes, P, prop, L)\n",
    "w_fe = u[2::DOF]\n",
    "ry_fe = u[4::DOF]\n",
    "print('max |w_FE - w| / delta          = %.1e' % (abs(w_fe - w_ref).max()/delta))\n",
    "print('max |r_y,FE + theta| / theta(L) = %.1e'\n",
    "      % (abs(ry_fe + theta_ref).max()/theta_ref[-1]))\n",
    "assert np.allclose(w_fe, w_ref, rtol=0, atol=1e-10*delta)\n",
    "assert np.allclose(ry_fe, -theta_ref, rtol=0, atol=1e-10*theta_ref[-1])\n",
    "\n",
    "xs = np.linspace(0, L, 201)\n",
    "w_cf, _ = timoshenko_cantilever(xs, P, prop, L)\n",
    "w_eb = delta*(3*xs**2/(2*L**2) - xs**3/(2*L**3))\n",
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "ax.plot(xs, 1e3*w_cf, 'k-', label='Timoshenko, closed form')\n",
    "ax.plot(xs, 1e3*w_eb, 'k--', lw=1, label='Euler-Bernoulli, same $\\\\delta$')\n",
    "ax.plot(nodes, 1e3*w_fe, 'o', mfc='none', label='BeamC, 10 elements')\n",
    "ax.set_xlabel('$x$ [m]')\n",
    "ax.set_ylabel('$w$ [mm]')\n",
    "ax.grid(True, alpha=0.3)\n",
    "ax.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "i = np.argmax(abs(w_cf - w_eb))\n",
    "print('largest difference Timoshenko - Euler-Bernoulli: %.3f mm at x = %.2f m'\n",
    "      % (1e3*(w_cf - w_eb)[i], xs[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9384d2c1",
   "metadata": {},
   "source": [
    "2. `Analysis` with `bu` and `uk`\n",
    "---\n",
    "\n",
    "`Analysis` receives the callables of the full model and the partition: the\n",
    "known dofs are those of `~bu`, with the prescribed values `uk[~bu]`, imposed\n",
    "exactly. The linear analysis, `static(NLgeom=False)`, solves the same system\n",
    "$\\boldsymbol{K}_{{C0}_{uu}} \\boldsymbol{u}_u = \\boldsymbol{f}_u - \\boldsymbol{K}_{{C0}_{uk}} \\boldsymbol{u}_k$ in a\n",
    "single step and returns `cs[0]` of size `N`, with the prescribed values at\n",
    "the known dofs. `calc_fext` returns zeros: the analysis is displacement\n",
    "controlled. The reactions are `calc_fint(u) - calc_fext()` at the known dofs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0558a060",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-30T23:31:29.337352Z",
     "iopub.status.busy": "2026-10-30T23:31:29.337215Z",
     "iopub.status.idle": "2026-10-30T23:31:29.343308Z",
     "shell.execute_reply": "2026-10-30T23:31:29.342771Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "increments [1.0] , u[tip_w] = 0.01\n",
      "max |u_Analysis - u_partition| / delta = 8.7e-15\n",
      "tip reaction 1.3576415826e+05 N\n"
     ]
    }
   ],
   "source": [
    "def calc_fext(inc=1., silent=True):\n",
    "    return np.zeros(N)\n",
    "\n",
    "def calc_kC(c=None, NLgeom=False, silent=True):\n",
    "    return KC0\n",
    "\n",
    "def calc_fint(c, silent=True):\n",
    "    return internal_forces(elements, prop, c, ncoords)\n",
    "\n",
    "def calc_kG(c=None, NLgeom=False, silent=True):\n",
    "    return csr_matrix((N, N)) # not used by the linear analysis\n",
    "\n",
    "an = Analysis(calc_fext, calc_fint, calc_kC, calc_kG, bu=bu, uk=uk)\n",
    "increments, cs = an.static(NLgeom=False, silent=True)\n",
    "u_an = cs[0]\n",
    "print('increments', increments, ', u[tip_w] =', u_an[tip_w])\n",
    "print('max |u_Analysis - u_partition| / delta = %.1e'\n",
    "      % (abs(u_an - u).max()/delta))\n",
    "assert np.allclose(u_an, u, rtol=1e-12, atol=1e-14*delta)\n",
    "R_an = calc_fint(u_an) - calc_fext()\n",
    "print('tip reaction %.10e N' % R_an[tip_w])\n",
    "check_reactions(R_an, nodes, L, P, rtol=1e-9)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f8f7a98",
   "metadata": {},
   "source": [
    "Both ways give the closed form to round-off. The same `bu` and `uk` serve the\n",
    "geometrically non-linear analyses, `static(NLgeom=True)`, where the\n",
    "prescribed displacements grow with the load factor, $\\lambda \\boldsymbol{u}_k$, see\n",
    "*Non-linear static analyses*.\n",
    "\n",
    "The checks are also run as tests:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a572f9ee",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-30T23:31:29.344901Z",
     "iopub.status.busy": "2026-10-30T23:31:29.344696Z",
     "iopub.status.idle": "2026-10-30T23:31:29.351671Z",
     "shell.execute_reply": "2026-10-30T23:31:29.351035Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P = 1.3576415826e+05, FE 1.3576415826e+05, Euler-Bernoulli 1.4000000000e+05\n",
      "tests passed\n"
     ]
    }
   ],
   "source": [
    "tpd.test_partition_by_hand()\n",
    "tpd.test_analysis_bu_uk()\n",
    "print('tests passed')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}