beam_buckling-checkpoint.ipynb 50.4 KB
Newer Older
Jeremy BLEYER's avatar
Jeremy BLEYER committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Eulerian buckling of a beam\n",
    "\n",
    "In this numerical tour, we will compute the critical buckling load of a straight beam under normal compression, the classical Euler buckling problem. Usually, buckling is an important mode of failure for slender beams so that a standard Euler-Bernoulli beam model is sufficient. However, since FEniCS does not support Hermite elements ensuring $C^1$-formulation for the transverse deflection, implementing such models is not straightforward and requires using advanced DG formulations for instance, see the `fenics-shell` [implemntation of the Love-Kirchhoff plate model](http://fenics-shells.readthedocs.io/en/latest/demo/kirchhoff-love-clamped/demo_kirchhoff-love-clamped.py.html) or the [FEniCS documented demo on the biharmonic equation](http://fenics.readthedocs.io/projects/dolfin/en/2017.2.0/demos/biharmonic/python/demo_biharmonic.py.html).\n",
    "\n",
    "As a result, we will simply formulate the buckling problem using a Timoshenko beam model.\n",
    "\n",
    "## Timoshenko beam model formulation\n",
    "\n",
    "We first formulate the stiffness bilinear form of the Timoshenko model given by:\n",
    "\\begin{equation}\n",
    "k((w,\\theta),(\\widehat{w},\\widehat{\\theta}))= \\int_0^L EI \\dfrac{d\\theta}{dx}\\dfrac{d\\widehat{\\theta}}{dx} dx +  \\int_0^L \\kappa \\mu S \\left(\\dfrac{dw}{dx}-\\theta\\right)\\left(\\dfrac{d\\widehat{w}}{dx}-\\widehat{\\theta}\\right) dx\n",
    "\\end{equation}\n",
    "where $I=bh^3/12$ is the bending inertia for a rectangular beam of width $b$ and height $h$, $S=bh$ the cross-section area, $E$ the material Young modulus and $\\mu$ the shear modulus and $\\kappa=5/6$ the shear correction factor. We will use a $P^2/P^1$ interpolation for the mixed field $(w,\\theta)$. "
   ]
  },
  {
   "cell_type": "raw",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
24 25 26
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
Jeremy BLEYER's avatar
Jeremy BLEYER committed
27 28 29 30 31 32
   "source": [
    "For issues related to shear-locking and reduced integration formulation, we refer to the :ref:`ReissnerMindlinQuads` tour."
   ]
  },
  {
   "cell_type": "code",
33
   "execution_count": 1,
Jeremy BLEYER's avatar
Jeremy BLEYER committed
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook\n",
    "\n",
    "L = 10.\n",
    "thick = Constant(0.03)\n",
    "width = Constant(0.01)\n",
    "E = Constant(70e3)\n",
    "nu = Constant(0.)\n",
    "\n",
    "EI = E*width*thick**3/12\n",
    "GS = E/2/(1+nu)*thick*width\n",
    "kappa = Constant(5./6.)\n",
    "\n",
    "\n",
    "N = 100\n",
    "mesh = IntervalMesh(N, 0, L) \n",
    "\n",
    "U = FiniteElement(\"CG\", mesh.ufl_cell(), 2)\n",
    "T = FiniteElement(\"CG\", mesh.ufl_cell(), 1)\n",
    "V = FunctionSpace(mesh, U*T)\n",
    "\n",
    "u_ = TestFunction(V)\n",
    "du = TrialFunction(V)\n",
    "(w_, theta_) = split(u_)\n",
    "(dw, dtheta) = split(du)\n",
    "\n",
    "\n",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
66 67
    "k_form = EI*inner(grad(theta_), grad(dtheta))*dx + \\\n",
    "         kappa*GS*dot(grad(w_)[0]-theta_, grad(dw)[0]-dtheta)*dx\n",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
68 69 70 71 72
    "l_form = Constant(1.)*u_[0]*dx"
   ]
  },
  {
   "cell_type": "raw",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
73 74 75
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
Jeremy BLEYER's avatar
Jeremy BLEYER committed
76
   "source": [
Jeremy BLEYER's avatar
Jeremy BLEYER committed
77
    "As in the :ref:`ModalAnalysis` tour, a dummy linear form :code:`l_form` is used to call the :code:`assemble_system` function which retains the symmetric structure of the associated matrix when imposing boundary conditions. Here, we will consider clamped conditions on the left side :math:`x=0` and simple supports on the right side :math:`x=L`."
Jeremy BLEYER's avatar
Jeremy BLEYER committed
78 79 80 81
   ]
  },
  {
   "cell_type": "code",
82
   "execution_count": 2,
Jeremy BLEYER's avatar
Jeremy BLEYER committed
83
   "metadata": {},
Jeremy BLEYER's avatar
Jeremy BLEYER committed
84
   "outputs": [
85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/ffc/uflacs/analysis/dependencies.py:61: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  active[targets] = 1\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/ffc/uflacs/analysis/dependencies.py:61: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  active[targets] = 1\n"
     ]
    },
Jeremy BLEYER's avatar
Jeremy BLEYER committed
116 117 118
    {
     "data": {
      "text/plain": [
119 120
       "(<dolfin.cpp.la.PETScMatrix at 0x7fdbf38b0990>,\n",
       " <dolfin.cpp.la.Vector at 0x7fdbf38b00f8>)"
Jeremy BLEYER's avatar
Jeremy BLEYER committed
121 122
      ]
     },
123
     "execution_count": 2,
Jeremy BLEYER's avatar
Jeremy BLEYER committed
124 125 126 127
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
Jeremy BLEYER's avatar
Jeremy BLEYER committed
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
   "source": [
    "def both_ends(x, on_boundary):\n",
    "    return on_boundary\n",
    "def left_end(x, on_boundary):\n",
    "    return near(x[0], 0) and on_boundary\n",
    "\n",
    "bc = [DirichletBC(V.sub(0), Constant(0.), both_ends),\n",
    "      DirichletBC(V.sub(1), Constant(0.), left_end)]\n",
    "\n",
    "K = PETScMatrix()\n",
    "assemble_system(k_form, l_form, bc, A_tensor=K)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Construction of the geometric stiffness matrix\n",
    "\n",
    "The buckling analysis amounts to solving an eigenvalue problem of the form:\n",
    "\n",
    "\\begin{equation}\n",
    "(\\mathbf{K}+\\lambda\\mathbf{K_G})\\mathbf{U} = 0\n",
    "\\end{equation}\n",
    "\n",
    "in which the geometric stiffness matrix $\\mathbf{K_G}$ depends (linearly) on a prestressed state, the amplitude of which is represented by $\\lambda$. The eigenvalue/eigenvector $(\\lambda,\\mathbf{U})$ solving the previous generalized eigenproblem respectively correspond to the critical buckling load and its associated buckling mode. For a beam in which the prestressed state correspond to a purely compression state of intensity $N_0>0$, the geometric stiffness bilinear form is given by:\n",
    "\n",
    "\\begin{equation}\n",
    "k_G((w,\\theta),(\\widehat{w},\\widehat{\\theta}))= -\\int_0^L N_0 \\dfrac{dw}{dx}\\dfrac{d\\widehat{w}}{dx} dx\n",
    "\\end{equation}\n",
    "\n",
    "which is assembled below into the `KG` `PETScMatrix` (up to the negative sign)."
   ]
  },
  {
   "cell_type": "code",
164
   "execution_count": 3,
Jeremy BLEYER's avatar
Jeremy BLEYER committed
165
   "metadata": {},
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/ffc/uflacs/analysis/dependencies.py:61: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  active[targets] = 1\n"
     ]
    }
   ],
Jeremy BLEYER's avatar
Jeremy BLEYER committed
183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204
   "source": [
    "N0 = Constant(1e-3)\n",
    "kg_form = N0*dot(grad(w_), grad(dw))*dx\n",
    "KG = PETScMatrix()\n",
    "assemble(kg_form, tensor=KG)\n",
    "for bci in bc:\n",
    "    bci.zero(KG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we made use of the `zero` method of `DirichletBC` making the rows of the matrix associated with the boundary condition zero. If we used instead the `apply` method, the rows would have been replaced with a row of zeros with a 1 on the diagonal (as for the stiffness matrix `K`). As a result, we would have obtained an eigenvalue equal to 1 for each row with a boundary condition which can make more troublesome the computation of eigenvalues if they happen to be close to 1. Replacing with a full row of zeros in `KG` results in infinite eigenvalues for each boundary condition which is more suitable when looking for the lowest eigenvalues of the buckling problem.\n",
    "\n",
    "##  Setting and solving the eigenvalue problem\n",
    "\n",
    "Up to the negative sign cancelling from the previous definition of `KG`, we now formulate the generalized eigenvalue problem $\\mathbf{KU}=-\\lambda\\mathbf{K_G U}$ using the `SLEPcEigenSolver`. The only difference from what has already been discussed in the dynamic modal analysis numerical tour is that buckling eigenvalue problem may be more difficult to solve than modal analysis in certain cases, it is therefore beneficial to prescribe a value of the spectral shift close to the critical buckling load."
   ]
  },
  {
   "cell_type": "code",
205
   "execution_count": 5,
Jeremy BLEYER's avatar
Jeremy BLEYER committed
206 207 208 209 210 211
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
212 213 214 215 216
      "Computing 3 first eigenvalues...\n",
      "Critical buckling loads:\n",
      "Exact:    0.31800  FE:    0.31805  Rel. gap 0.01%%\n",
      "Exact:    0.93995  FE:    0.94033  Rel. gap 0.04%%\n",
      "Exact:    1.87267  FE:    1.87415  Rel. gap 0.08%%\n"
Jeremy BLEYER's avatar
Jeremy BLEYER committed
217 218 219
     ]
    },
    {
220 221 222 223 224 225
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/bleyerj/.local/lib/python3.6/site-packages/matplotlib/font_manager.py:1241: UserWarning: findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.\n",
      "  (prop.get_family(), self.defaultFamily[fontext]))\n"
     ]
Jeremy BLEYER's avatar
Jeremy BLEYER committed
226 227 228
    },
    {
     "data": {
229
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEACAYAAABYq7oeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXdYFOfWwH/D0gVciliw4NpL1CyWqMQSIcU0TVBjyr1JvggxxdQrMcm90ZRrIL0awUTTo2CKid4oWDAYYoSNvQZEsRdYupTd+f4YWXrdOji/5/Fxl3nnnTNnZ+bM+57zniOIooiCgoKCgkJNnOwtgIKCgoKC46EYBwUFBQWFeijGQUFBQUGhHopxUFBQUFCoh2IcFBQUFBTqoRgHBQUFBYV6KMZBQUFBQaEeinFQUFBQUKiHYhwUFBQUFOrhbG8BWoqfn5+o0WjsLYZDUF5ejqurq73FcAgUXVSj6KIaRRfVZGRkXBBFsVNr95ONcejUqRPp6en2FsMh2LJlC5MmTbK3GA6BootqFF1Uo+iiGkEQjrVlP9lMK3l6etpbBIdh+PDh9hbBYVB0UY2ii2oUXZiPbIyDwWCwtwgOQ2Fhob1FcBgUXVSj6KIaRRfmIxvjUFZWZm8RHIasrCx7i+AwKLqoRtFFNYouzEc2xkFBQUFBwXbIxji4ubnZWwSHITg42N4iOAyKLqpRdFGNogvzkU20kkqlsrcIDoOfn5+9RXAY5KYLo9HIiRMnKC4utkrfBw4csHi/cuRK0YWLiwuBgYH4+PhYvG/ZGIeSkhJ7i+Aw6HQ6JUzvMnLTxYULFxAEgQEDBuDkZNmBe2FhId7e3hbtU65cCboQRZHS0lJOnjwJYHEDIZtpJQWF9oBer6dz584WNwwKVx6CIODp6UlQUBDnzp2zeP+yuUKdnWUzyLE6vr6+9hbBYZCbLgwGAy4uLlbpW5l6reZK0oWHhwcVFRUW71c2xsHDw8PeIjgMygKfauSoC0EQrNKvslC0mitJF9a6nmRjHJRFLdWkpKTYWwSHQdFFNco9Uo2iC/ORjXFQqEYURXuL4DAoulBQsA4WNw6CICS1oI1GEIT5giCEXf5fbWk52jPWGkbKEUUX8iU5OZnY2NgGVzPrdDpiY2PR6XR2kMwyJCYmEhsbi16vt7cobcJixuHygz4SCGtB8wRRFGNFUUwG4oD45nZo72FprWHixIn2FsFhUHRRjTXukaysLKKjoxEEgdjYWOLi4oiLiyMqKorExESz+g4LC2PHjh0NGgCtVktmZmatTMwhISEkJye3qG9HeF5ERESwdOlSs/oIDw+3kDStx2IhQJcf9AiC0KQ2BEHQArk19tMLgtCsQSktLTVbxvbCrl27ZOmItQaKLqopKSmxuCNWo9GwYMECYmNjmT9/fq1tgiCQkZGBVqs1q//G6NOnT63vCQkJTbaviTV00RbaWoMmOTmZrKysFhtDa2CP+FANUHeclSsIglYUxUbHkJWVldaVSkbk5eXZWwTLYjRAbhboj0HxRSi5AJfyQVCByhmcXKBDAHTsLv3z6Q7OUiGXdqcLM7B15mKtVktycrJZxqE1tOZBK/cszmFh0vtydHS03WSwh3Focb6Dy9NUkQCBgYFs2bIFkC4Sb29vdu3aBYC/vz9Dhgxh69atgLQmIjQ0FJ1OR0FBAQAjR47k7Nmz5OTkANCvXz/c3NzYu3cvVf3379+f1NRUQMrlNHbsWNLT0ykqKgJgzJgxnDhxwrQiccCAAahUKvbv3w9Aly5d6N27N2lpaYAUfjtmzBi2b99uGvmMHTuWo0ePcubMGQAGDx6MwWDg0KFDAAQFBdG9e3e2b98OgJeXFyNHjiQtLc2UmVYURfbv329a+DJ06FDKyso4cuQIAD169KBz586mIbmPjw9arZbU1FSTkZ0wYQL79u3j4sWLgBQSWlhYaJr/DQ4Oxs/PzzTk9/X1Zfjw4aSkpCCKIoIgMHHiRHbt2mV6QGu1WnJzc8nOzm7yd0pL/gXfvJ345e+lC+cwntmHk7G8pZcFosqVgg69KfTuj49bMJQOJ31fpsP9TqGhoRw+fLjW71RZWWmKpHFxcWHgwo0tPm9LcmhRGOXlks5dXV1RqVSmc1epVHh4eJj0WSVvcXExRqMRgIMHD5KXl8fs2bMpLCwkNTWVxx57jCeeeIIHHniAn3/+mUcffZQ1a9YwYsQIBEHg8OHDrFu3jmHDhpGfn8/NN99MZWUlpaWlFBYW8tRTT6HT6Vi4cCGTJ0+moqICo9FIYWEhO3fuZN68ecydO5e7776bTZs28dJLLxEVFUWPHj3Izc1l8+bNLF26FFEUMRgMrFmzhpycHFQqFYcOHWLy5MkcO3aMxx9/nKKiIlMwwx9//MH8+fN54IEHCA4OxtXVlcTERObNm0d2djY5OTkUFhby6KOPArB792727dtHly5dyM7OZtKkSVx11VWUlpaSkZHBli1bGD16NLm5uVy4cIHCwkI8PDxwdnZm7dq1ZGdno9Fo+O2333jxxRdNv4m3tzclJSUmw1Y18qnSf1O/06VLl0hNTW3wuddWBEtHewiCIIqi2KiXUBCECCBKFMXwGn/LA6Y0NXK4+uqrxb/++suissqVgoICq+RSsSr5J2D3Ktj/E5zeBdS57jr2AD+NNELwDAD3jlIbYyUYKqDorNRH/gnIz6m9r5MzBIfCwFtg8DTwanVFRJtx4MABBg0aZPoe/Nxau8iR/frNLW6r1+vx9fUlIyMDkPwQK1euZMGCBbVGDdHR0fTp04fIyEhAmi+PiYlBq9Wi1+uZMmWKqY/Y2FgALl68yKhRowgLCyM5OZmIiAhTf7GxsajValN/db9HR0ej1+tN8/ohISGmqSeDwUBAQIDpxaVPnz5kZGSgVjcc+9JQXwsWLDDJ4+vrS15eHllZWURFRZGUVB13ExISQkZGRr1zrHvcuvsmJiaSlZVVb7quJlXHbY6611VNBEHIEEWx1VbCHiOHLBoYPTRlGED+w0RLkpubKw/jYKiAfT/CX1/C0a2YDILKDXqNBc1k6DEaAgeDRysC1kr1cEoHJzIoPbgBjzMZkLVF+rf+eclAjHpI6tvBo5la85BujrKyMqtmL64yBFqtloiICMLDw4mKijI9QP39/Rvdd9WqVbXeYqseiNHR0WRlZTFlyhQ2bmzdKMrf37+WX8LPz88UGXThwoVaSRnVajW5ubmNGoeG+mpoGisxMbHeNJpGozH5COq+qdc8XmJiImq1upYDvjF5HAGbGAdBEDRAriiKelEUdTXDDy9va9brohT7qSY7O9uxUxKXFULG5/DHEig4If1N5QYDp8JVM0EzCVzNcBZ6qKHPddDnOraLo5h03zA4vB72/QB/J8GeVdK/rsNh4nMw4CaHNxKWoLy83Kap7WfMmEF0dHStt/22oFariYmJYcaMGbXeyM3B09Oz3oO3rc7h5mhNqKqfn18tI+vIWDKUVSsIwvzLn2PqRCDFADNrfJ9Ttc4BiADmWEoOBTtSXgK/vQXvDIENL0iGIaA/3Pw2PHsYZqyQDIQ5hqEhPP1gxGy4ZxU8sQtCn5ampk7vgu9mQ/x1cCQZlAVzFiUpKaneA67mg7JmGOrMmTPrRd5UvUFrNBrCwsLQaDTExcVZTL5Zs2YRFxdHYmJiq0cljREREVEv9DY3N5ewsDBmzpxZ65xB0kdubq5p37rb7RmN1ByWDGXVATogtoFtMxppCy0YNYBS7Kcm1noDajOGStj1DWxeDIWnpL/1HAfj50G/G8CKGUjr6ULdE8JegonRkLEcfntbmoL6+k7oGwY3xYJ/n4Y7kzmurq4W7zMrK8s0D1/zwZ2RkYGfn1+tOP7IyEiio6NND7ywsDAWL15MfHw8arWahIQEYmNjTT6IgIAA03SMRqMhJCSEqKgo9Ho9YWFhrFy5Ej8/P8LCwtDr9U1+z8rKIj09naVLlxITE4OHh4dpnYSfnx8rV64kKirKFAVUE51O12Rfq1atQq/Xm8J5o6OjTeeh0+lISEgApBFQfHx8rXMEaeosPj4ejUZDTExMre2Njbp0Oh3Jycno9Xqio6MJDw9vUHZrYnGHtLVQHNLV5OXlOU420pMZ8POTcGa39L3LMAh/GfpMtsnhm9VFeTHsWAZb34KyfGl6K/QpCH0SXGyfzLEpx6G5VFZWKtmLL7N+/XqOHTvGzJkzTc7g6Ojoek709oI1HNKyya2kFPuppio01K5cKoB18yF+imQYOvaAO5ZBZIrNDAO0QBeuHWD8E/B4OgyfDYYySHkdPrkWTso3NUNDKAtFq/n1118JCwsz+R00Gg2zZs1qMFWHQsMorxkKrScrBX6cCwUnpYVqYx+FSc9JD2JHxSsQpn8C2n9II50Lh+DTcEnu0KfB6crJ/38l8PLLL7NkyRLUajV+fn6mef+qMFiF5pGNcVCGy9U0FTJoVSouwaZXIO1D6Xs3Ldz2PnS5yj7y0AZd9BoHUSmQvBC2fwKbXpWc1TOWg083q8hoK66kAjfNoVKpmlw/oNA8splWUor9VDNkyBDbH/T8ISnqJ+1DabQw+QX4vyS7GgZooy5cPOCmGLjvB/DuCjl/wNIJkL3N8gLaEOUeqUbRhfnIxjgoxTuqqUoTYjP2roa4yXBuH/j1kYzCxPlS3iM7Y5Yu+lwHD6dC7wlQfB6+uA3++ES2Ia9V6S4UFF1YAtkYBwU7UFkO/3sOEh+EimK4agZEbYXuIfaWzHJ0CIB7f4Bxj0upOn6NhjWPS6u7FRSuYGRjHJSiLtXYxP9SfAG+uB22L5Gyok59E+6IBzcv6x+7FVhEFypnuP5ViFgOzh5Suo9vZkkrvRUUrlBkYxy8vBzroWRPQkNDrXuAcwck/8Lx36U5+Qf+B6PnOGQKCovqYugd8M+fwdMfMjfC8pug4LTl+rcyjlDgxlFQdGE+sjEOyjqHaqxaOvFIMnx6vVRbodvVMGcz9BhlveOZicV10WMUPJQs+VbO7JF0kSuP2Pji4mJ7i+AwKLowH9kYByUrazVVudotju4L+GYmlBVImU3vXwc+Xa1zLAthFV34aSSne9BIyD8Oy6dK0VoOTlWtBQVFF5ZANsZBwYqIIqS8ITliRQNc+6w0/27pBHlyooM//ONH6BUKhaclA3Fmj72lalckJycTGxvb4KplnU5HbGysdUfJViYxMZHY2NhWZW11JGRjHDp0cODVtzbGnOpO9TAaYN2zsPlVQICb34Ip/7ZqsjxLYlFd1MXNG+5JgD5TpNKlK26GUzutdzwzsUbN5KqcRIIgEBsbS1xcHHFxcURFRZGYmGhW32FhYezYsaNBA6DVak2J86oICQlpcRZTR6gfHRERUSsxYWvQ6XTExcURGxvLjBkz7JL2w/6B6i2kokIJLazi7NmzlnHQGyrg+0jY972UkO7OZTD4NvP7tSEW00VjuHrC7G8h4QE4tBa+nA73r4XOg613zDZSUVFh8VXSGo2GBQsWmDKS1kQQBDIyMsxKZNdUhuGaxXcAU5W3lmANXbSFtmRQ1uv1pKenm1J9JCcnEx4eTmZmpqXFaxJ5vB6Cqd6tAqY62GZRcQlW3icZBldvuO972RkGsJAumsPZTapF0e96KM2VFstdOGL947YSW79AabVam9YjaM2DVs4vk1lZWcTExJi+jxw5kqysLJtPT8lm5KBgQcqL4dvZcDQFPHzh3tUQ1I4WtlkDZ1eY+SV8O0sqR/r5rVKIr19v8/pd2NEi4gG0KnhzYb5Zx9LpdOj1+lpvt1FRUURHRxMZGUliYiJz5sxh48aNppFFVY2CqloGdesTREVFkZ6eTkxMTL1tOp2OOXPmEBUVRWRkJMnJyURHRxMVFYVGo0Gv15OUlFRrGqeqVgRAZmYm4eHhZGVl1Uu+V7cvkEYpVSVMqx7MVSMnnU5Heno6Go2GrKwsU6Gihs6xKuFfXZk0Gg1JSUm1jEAVWq22Vh3q9PR01Gq1zUuKysY4uLu721sEh6Ffv35t37msCL6eIa1h6BAoOV072yFXk4UwSxetxcUd7vpG0t+xbdIU0/9tkDK+XgFU+QaysrJYuXIlCQkJpgdWWFhYrcI1defb9Xo9c+bMMT30qvwXNbeHh4c3Okev1WqZNWuW6XtYWBhhYWFkZGSYHvaLFy82PXjd3NyYMWMGeXl5gDRFtWDBggYL5jTUV3R0NDqdznROvr6+zJ8/3+SDqVnONCQkhIyMjHrnWNVPFVUjgqp9axYQqktNQ7B06VLi4+Mb1Is1kY1xUFZIV9PmqnhlRfB1BBxPA+9u0oKvgL6WFc7G2LxCoGsHuHul5Jw+vUvS5/1rJed1WzDzDb4mFRUVuLi4WKy/utSsfRwREUF4eDhRUVGmB2hTGXJXrVpVK3ig6oFY9XY+ZcqUVpfy9Pf3r+WX8PPzM0295Ofn4+fnZ9qmVqvJzc1t9O27ob4amsZKTEys52PRaDSmEUHdAImax0tMTEStVtdywDc3GoiLi2PWrFlm1+luC7LxOSiFTKrZu3dv63eqaxju/0X2hgHaqAtzcfOGexLBt7dkIFbeK+WhsjOXLl2y6fFmzJhR6824rajVamJiYpgxY0bzjVuIu7t7vQevtcrrtsYX4Ofnh1arNRnYpupLJCcno9Fo7GIYQEbGQcEMGjIM7bSOss3wCpSc+B06ST6IH+fCFbbwKikpqd5bdM0HZc0w1JkzZ9ZzXle9QWs0GtO8fc2pJnOZNWsWcXFxJCYmtnpU0hgRERH1Qm9zc3MJCwtj5syZtc4ZqOV3iIiIqLe9MYe+Tqcz1bQGzA4bbguymVay5nBZbgQGtmKOu6IUvpvdbg1Dq3Rhafw0kjN/+VTYmyh9v+4Fu4ljjYSMWVlZJj9AzQd3RkYGfn5+tXwEkZGRREdHmx54YWFhLF68mPj4eNRqNQkJCcTGxpqctQEBAbUctCEhIURFRZmc1StXrjQ9IPV6fZPfs7KySE9PZ+nSpcTExODm5mZaJ+Hn58fKlSuJiopq0Oeg0+ma7GvVqlW1/APR0dGm89DpdCQkJADSCCg+Pr7WOYI0dRYfH49GoyEmJqbW9oZGBVXTbDUNrT1GEIIok9z1ISEhYk1Hz5VMiwvJV5bBd/fA30ng1VmKrmlHhgFaoQtrciRJSjsiGmHaJzBidqNNmyoEby6iKCq+ucskJSVx9OhRZs6ciVqtNjmSFyxYYNa6DEelqetKEIQMURRbvVpUNtNKSvGOalJTU5tvZKiU6jD8nSRlGf3HT+3OMEALdWFt+oXDTbHS5zWPQ7Z9ZFLukWrWrVtHWFiYye+g0WiYNWuWXVYayxXZTCsptAKjEX56FA7+Au4dpXKYgdZ5W1W4zOg5cDFTqn/x3T0wZ1O7NMZy4eWXX2bJkiWo1Wr8/PxM8/5NOYAVaiMb4+Akk1w/tqDJ8E1RlKqZ7f4OXDrAPauh63DbCWdjbB7K2hQ3vAZ5R+Hwr/Dd3VLq77aGuLYBZUqpGkEQGlw/oNByZPPEVRLvVTN27NjGN27+L/wZBypXmP2NQ9disARN6sLWOKmkankBA+D8Qfg+yqYRTEpBrGoUXZiPbIyDUuynmrrhcCbSPoatsSA4QcRnoJlkS7HsQqO6sBfuPlKiPveOUqK+lNdtdmilwE01ii7MRzbGQSn2U02DjsddK2H9Aunz7R/BoFttK5SdcEgnrH8fyTgLTpASA/vX2OSwSoGbahRdmI9sjINCExxJgp8ekT5f/xqMuNu+8ihA3zAIf1n6/OMjDpnFVUGhKWRjHBSfQzVjxoyp/nIiHVb9A4yVMP4JGPeY/QSzA7V04WiMfQyGTIfyQinFRpl1RznKPVKNogvzkY1xkHN+dktz4sQJ6cP5w1JajIoSGHEPhC2yr2B2wKQLR0QQ4LYPqh3UP8+z6uGUmifVKLowH9kYB+XHrubkyZNQcBq+ugNK86D/jXDr+9LD6Arj5MmT9hahady8YdZX4OoFe1dDWaHVDqW8QFWj6MJ8ZGMcFKpRVRZLI4b8HOg+CiKWg0o2S1auPDr1l4IEAEr1UrElBZKTk4mNjW1w1bJOpyM2NrbB+tJyITExkdjYWJtXcLMUsjEOSrGfy1SWMTrrXTi7F/z7wuyVUp3jK5QBAwbYW4SWMWQajJkLiJCXLaU3sTDWWBBYlZNIEARTgZ64uDiioqLMzhQaFhbGjh07GjQAWq3WlDivipCQkBaXJXWExZF1Cx61hqqKcomJiURFRdkl7YdsXjeV1Z9IC6p+eBi3U39KifTu/R46NF5g5UrAEYrIt5jwl0H3OxjKQX9cKjFqwevaGveIRqNhwYIFDVYsEwSBjIwMsxLZNVVjoWbxHZBKd7a0JoOjPC/aWkNiypQpHD161FSkKCoqqlb1OVsgm5GDUuwHSPo37PueSpXH5WIzvewtkd3Zv3+/vUVoOc6ukjEXVFCWD8XnLdq9rYv9aLXaFr/JW4LWPGhtrQtLU2UYgFq5oWyJbEYOVzxpH0Pah+DkzL4hzzG86zB7S6TQFpycQd1VysFUcIqrVk+xixh7/rnHrP11Oh16vd6UyC45OZmoqCiio6OJjIwkMTGROXPmsHHjRtPIomqqpKqWQd3aClFRUaSnpxMTE1Nvm06nY86cOURFRREZGUlycjLR0dFERUWh0WjQ6/UkJSXVmsapqhUBkJmZSXh4OFlZWfWS79XtC6RRSlUJ06ysLPR6vWnkpNPpSE9PR6PRkJWVZSpU1NA51n2o16xfkZSURExMTIP6rVtDurF21sSixkEQBA0QAegALRAnimKD3hhBEGKA+YAeSAeiRFFsdGLtii72s+8HWP+89Pn2j3Fza7+J9FpLly5d7C1C6/FQQ3kni48crE2VbyArK4uVK1eSkJBgeoiFhYXVKkZTd75dr9czZ84cqmqyVPkvam4PDw9vdI5eq9Uya9Ys0/ewsDDCwsLIyMgwPewXL15sevA6OzszY8YM8vLyAGmKasGCBQ0W+2mor+joaHQ6nemcfH19mT9/vskHU3OKJyQkhIyMjHrnWNVPFVlZWcTExJj2rVlAqCGysrJITEwkPDy8QbmtjaVHDgmiKIYACIKQDsQDjRWGzRRFscUTg47gYLIL2dvg+0hAhLCFMHwWvcvK7CyU49C7d297i9A2fLpBWRF7wr8AD19Q9zLb/2A0Gq2avbhqBFBV/zg8PJyoqCjTA9Tfv3H/16pVqxg5srreTNUDsertfMqUKa0u5env71/LL+Hn52eKDCotLcXPz8+0rWruvm5d6ab6amgaKzExsZ6PRaPRmEYENc+x6rg191Wr1bUc8I3JU9Xv/PnziYuLIzw8XL4+B0EQtIBpDHV5xGAxc+eQOXSszbkDUolPQzmMmgPjnwQgLS3NzoI5DrLVheAEvsHS/6V5UGr+nLKtk83NmDGj1ptxW1Gr1cTExDBjRmPvka3HxcWl3oO3rc7h5mhNqKqfnx9ardZkYBuqL5GVlUVsbKzpe1X9bVtHLFly5KBBmiKqSa4gCFpRFBsKVlYLghBxeZ9wYHHdKShBECKBSJBqBW/ZskU6kEaDt7c3u3btAiSrP2TIELZu3QpItXRDQ0PR6XQUFBQAMHLkSM6ePUtOTg4A/fr1w83Njb1791LVf//+/U2Vxdzc3Bg7dizp6ekmwzRmzBhOnDhhWng1YMAAVCqVySnapUsXevfubXpgeXh4MGbMGLZv325yqI8dO5ajR49y5swZAAYPHozBYODQoUMABAUF0b17d/7aug6tbj7uZfkw8BbSfKdRlpICSOUg9+/fz7lz5wAYOnQoZWVlHDki5e/p0aMHnTt3NoUB+vj4oNVqSU1NpbJSCqGcMGEC+/bt4+LFiwAMHz6cwsJC0wUYHByMn5+f6S3H19eX4cOHk5KSYipHOXHiRHbt2mUaumu1WnJzc8nOzrbZ71T1QLTX77R9+3ZAShE9cuRI0tLSKLs8sgsNDeXw4cO1fqfKykoKC6WFcC4uLrh6B+FUkIOoz+GSwQkPb1/T9qp+S0tLTYknPTw8MBgMpkWhrq6uODs7U1JSgsFgoKSkBE9Pz1p9eHt7m7YDeHp6UllZWasPlUplOneVSoWHh4dJn1V9FRcXmxLaeXp68uuvvzJs2DAKCwtxc3PDYDBw9uxZCgsLcXZ2Jj09neLiYgoLC5k6dSoxMTG1+jh06BCVlZV06dKFMWPG0KtXL5YsWcK9994LSAvZjEaj6fhVei0qKkIURcrKyjAajZSWllJZWYnBYDCdl8FgYNq0aXz00Ud06tSJH3/8kcLCQpycnOjQoYOpjyodV1RUcOnSJQoLC3F3d8doNJpkr5rSLiws5MYbb+TJJ5+spZcLFy4wZcoUBg8ezJIlSygsLDT9Trm5uRw/fhwPDw+mT5/OzJkzKSwsRKVS4enpyZo1a5g8eXKt32nv3r2cPn3adD579uxBrVYTFBREZWVlg7/TpUuXSE1NbfB+ajOiKFrkH9JDPKHO3zIBbSPt1TU+a4GMpvofOHCgeMVQmi+KH48XxZd8RDE+TBTLS2pt/uOPP+wkmOMhN13s37+/9h+MRlHMPSqKJ3WieO6AKBoNbe67sLDQPOEaIDMzU5w/f74IiEuXLjX9i4yMFCMjI2u1zcvLEyMjI8WkpCQxKSlJjIiIECMiIsS8vDxRFEUxIyNDjImJEZOSksSEhARx8+bNolarFSMiIsTMzExx6dKlIiDGxMSIGRkZolarFcPCwsTMzMxmvyclJYlqtVqMjIwU8/LyxMLCQjEyMtLUJiIiQkxKSmrwHJvrq6ZcoiiKSUlJpvOIiYkRMzMza/VV8xw1Gk0tHdTcNyEhoVG9JyQkmHQdEREhZmRkNPk71buuagCki214pgviZetpLpdHAVGiKIbX+FseMEVseORQd38R8BUbcWCPHDlSdLjc/dagslxa/Xw0RVrk9uCGK34tQ3uiwULwRoOUe8lQDl6B4BNkH+HaEVXTMDNnzkStVpscyQsWLDBrXYaj0uB1dRlBEDJEUWz1EMKS3qsswK/uHxsyDIIgaAVByGigbaOTd1dE8Y6q2s9HU6BDINy7ukHDUDWVodBOdOGkkhzSAEXn2px/6Yr0yzXC2rVrCQsLM/kdNBoNs2bNsstKY7liMZ+DKIqsZGPGAAAgAElEQVS6mqsSL4e1Jtf5nnvZAGQBS2tsCwOaXIt/RRTv2LgQ9qySkrTdkyA5LBtAWRBYTbvRhZsXeHWBojOQdww6DWx1vixLzQK0B15++WWWLFmCWq2utYisIQewQsNYOpR1jiAI86le5zCnxrYYIInLax8EQci67HAG6FOn7ZXHH5/AtvekRVIzv4BuI+wtkYKt8e4CZQVSCvb8nMvRTI6RBkKONLZ+QKFlWNQ4XJ5CqppGSq6zbUad761ad9+uC4bv/R5+fU76fNsH0LfpVbNjx461gVDyoF3pQhAkg3D+IFzSS+Gtni33NykFbqpRdGE+ssmtVNZeF34d3Qo/RGFa5NaCEp9Hjx61tlSyod3pwtkNOnaXPuefgMqWX/ft9h5pA4ouzEc2xqFdFu84swe+u0eKUhnzsGmRW7O7XY69V5CnLpr1DXj4gXtHEI2gPwYt9CVUrWFRuLJ0YS1fk2yMQ7sj9yh8dac0xzxkOtywWJlfvgJQqVTNv+gIAnTsKfmfyoulCCYFhUYoLS21Su452RgHDw8Pe4tgOQrPwpfToegsBF8L05dCK3LiDB482IrCyQu56UKtVnP27Nnmo+9UztXhrYWnobyk2b6VgljVXAm6EEWRkpISTp48SWBgoMX7l03K7nYTpncpXxox5B2FrsPhrm+keeZWUJUGQUF+uggICODEiROmNBzNUloqrXs4kSuFujYxujQYDPIqfmRFrhRduLi40LlzZ3x8fCzet2yMg9yLdwBQUQrfzoaze8CvD9yzGtxb/6MeOnSIrl27WkFA+SE3XTg5OdGzZ8+W71BeDJ9cC7mZMO5xuP7VRptu2bKFSZMmmS9kO0DRhfnIZlpJ9lSWw6p/wrFt4N0V7vsBvDrZWyoFR8e1A8bpn1DmpKIo7SPyj6zHKF4BC0IV7I5sRg6urq72FqHtGA1SuOqR9VLu/nu/N6vEZ1CQknunivaoi8N5h9lxZgd7Luxh74W9HCs4Br0un+fvz+Kx3Z3eag2ajhqGBgwlvFc4gZ6B7VIXbUXRhflYLPGetdFqtWLNIhmyQRTh5ydA9zm4esM/f4KgELO6LC0tbV8OejNoL7oorihm3dF1rD68mn0X99Xb7uLkgouhAkE0UlwneEFAYFSXUYQFhTF94HTcndu/M7Y52st1YQnamnhPNsZhwIABYoudeI6CKMKGF6Xaz87u0ogheLzZ3SrzqdXIXRdlhjK+3P8ly/Yso7hCSi7p7erNdT2uY1inYQwNGEo/3364OLnAuYOwdAL5YgVZU/9Lpk8nfjvxG7+d/I0KoxQe28mjE5HDIrmz3524qK7c0rpyvy4sSVuNg2ymlWSHKMLGlyXD4OQCM7+0iGFQaB+Iosim45t4I/0NThZJRYm0gVoi+kcQ3iu84bf/wIEQ9hId1z/P1Zvf5upH/iCifwSF5YVsPL6RpTuWcqL0BK9tf40V+1bwdMjTXB98vY3PTKG9IBvjILuwtJQYSH0bBBXMWA79LXeTtus8U61EjrrIL8vn39v+zeaczQD0VfclenQ013S9pvmdx8yFg+vgWCr88iTM/AJvV2+m9Z1Gt7xu5Afm8+FfH5KVn8UzKc8w9fhUnh/zPB3dOlr5rBwLOV4XjoZsppVkVeznt7ekUYPgBHd+CkPvsLdECg7CnvN7eDblWU4Vn8Lb1Zt5V88jon8Ezk6teE/LOwZLxkF5EdwRD8Nm1tpsMBpIOJzA2xlvU1pZSqBnIK+Me4VxQeMsfDYKcsARiv1YFdkU+9n6hmQYEGDaJ1YxDFW1jxXkpYtvD37LP379B6eKTzHUfygJtyZw18C7WmcYQIp0u3Gx9Hnds5AvTUtV6ULlpOKugXeReGsiwzsN51zJOR5OfpgVe1e0n8WkzSCn68JRkY1xkEWxny2vw6ZXkQzDxzB8llUOo2ScrEYOuhBFkXcy3uG/2/9LpbGSuwfezec3fU6QlxnhllffB/1vlFbcr3kcRLGeLnr69OTzGz/nkRGPICLyVsZbLExbSIWhHSaxrIMcrgtHRzbGwaERRdj0GmxZLE0lTV/aotTbCu0fg9HAorRFfLb3M5wFZ/4b+l8WjFmAq8rMdTuCALe+L2VwzdwI6Z812EzlpGLu8Lm8Pelt3FXufH/keyKTIskvyzfv+ArtHtn4HEJCQsSMjHplp+1PzXBVwUmaA74qwqqHrKysxNlZNrEEVsWRdVFhqOC5355jw7ENuKnceHvS20zoPsGyB9n3AyTcDy6eVM5JwTmwf6NN917Yy+ObHudC6QUG+Q0i/vr4duuoduTrwta0e5+DQw4TjQZpSJ/2oZRe+c5PrW4YAA4fPmz1Y8gFR9WFwWhgQeoCNhzbgJeLF0vDl1reMICU7n1oBFSUUL7q/6RrshGGBgzlm6nf0MO7BwdyD/DQhofQX9JbXiYHwFGvCzkhG+PgcMV+Kssh8UH460tpgdtd39osKuncOSW/fxWOqAtRFFn852LWZ6/Hy8WLT2/4lJDO5q2Kb5Kpb4BXFzwv7IbfP2iyaVevriy/YTm9fHpxMPcgD214iLxLedaTzU444nUhN2RjHByKSwXwzUzY/yO4+UhJ9Cy4jkFB3ny862NWHlqJq5Mr71/3PoP9rVxzwtMPbv9I+rz5NThbP/1GTTp36MxnN3xGsE8wh/IOEZUURVF5kXVlVJAdsjEODpMnpeAULJ8KWZvBMwD++TP0sm38+NChQ216PEfG0XTx3cHv+GTXJzgJTrwx8Q1GdRllmwP3C6N0yGyp5Oz3UdLItgkCPQP57IbP6OndkwO5B3hqy1PtKorJ0a4LOSIb4+AQjvNzB2BZuFSPwb8vPJQM3UbYXAyH9L/YCUfSxfbT23n9z9cBWDh2Idf1vM6mx88b9TT4BkvXZ8rrzbbv5NmJT8I+wc/djz9O/8F/fv9Pu0kH7kjXhVyRjXGwe7Gfwxvg0+uh4AT0GAP/lwR+ve0iypEjR+xyXEfEUXRxsugkz6Y8i0E08NBVDzG933Sby3A4+5QURo0Aqe9Azp/N7tPDpwcfT/kYD2cPfsn6hXd171pfUBvgKNeFnJGNcbAbogip70o+hrICGDwN/vGTNM+roACUVpbyxKYn0JfpCQ0K5bERj9lPmJ7XwPh5IBqlGiLlzWcWGBIwhLcnvY2z4Mzyvcv54cgPNhBUwdGRjXGwS7Gf8hL4PhKSXwJEmPwCzFgBLvb1f/To0cOux3ck7K0LURR5adtLHMo7RC+fXsRMiEHlZJ8kkSZdTH4BAodAbpa0BqcFhAaF8sI1LwDwyh+vsPv8bmuJaRPsfV20B2RjHFxcbJyb/vxhWDYF9qwClw4w6yuYOL/JAu+2onPnzvYWwWGwty4SDifwv+z/4ensyXuT38PH1fKF3luKSRfObnBHHKhcpZXTR5JatH9E/whmDZhFhbGCpzY/xfmS81aU1rrY+7poD8jGONg08d6ulRA3Cc7tv+x4ToJBt9ru+M0gm+y0NsCeuvg7729id8QCsHDcQvqo+9hNFqijiy5DpREEwE+PQklui/qIHhWNNlDLudJzPLXlKcoNTUc9OSrKPWI+sjEONqFUDz88DD9EQkUxXDUDIrdA5yH2lkzBwbhUeYl/bf0XZYYypvWdxk29b7K3SPUZ9zj0HAdFZ6XaDy2I+HNRufDWpLfo7NmZXed38caON2wgqIIjIhvjYPViP38nSznyd30rrXi+9X0pT5Kbt3WP2wZ8fOw3deFo2EsXb6a/yd/6vwn2CWbB6AV2kaEu9XThpILpS8DVC/b/BLtXtqifAI8A3pv8Hs5Oznx36Ds2ZG+wgrTWRblHzEc2ifesVuynJBeSF4Luc+l70EiY/gkE9LP8sRTaBVtytvD4psdxcXLh66lfM8h/kL1Fapq/vpKmlly9Ye42qR5EC/j6wNe8/ufreLl4seqWVfTwUZy8cqTdJ94rKrLw8n6jETI+hw9CJMOgcoUpL8GD6x3eMKSmptpbBIfB1rrIL8tnUdoiAJ7QPuFQhqFRXYy4R/KZlRdK4a1NJOeryd0D72ZKzykUVRTx7NZnZeV/UO4R85GNcbDoCCc7FT4Ng5/nQWkuBF8LUb/BtU+DyvHT/FZWVtpbBIfB1rqI+TOGC6UX0AZquW/wfTY9dnM0qgtBgFveA68ucDwNtrVsoZsgCLw8/mWCvILYf3E/b2e8bUFprYtyj5iPbIyDRcjZAV/cDituhpMZ4N0VIj6T8iMFDrS3dAoOzpacLfyc9TPuKndeGf8KToKMbp8O/jCtKjnff+HUXy3azcfVhzcmvIGz4MzXB74m9aTyRn6l0P59DkYDHNkA2z+BrC3S39x8pEiOa+Y6pMO5OYxGI05OMnowWRFb6SK/LJ9pP03jQukFokdFc+/ge61+zNbSIl2smw9/LgW/PhC1Fdy8WtT3sj3LeE/3HgEeAay+bTV+7o6dIUC5R6pp9z6H0tLS1u2gz4Ft78H7V8O3d0mGwaUDXPsMPLFLWtAmQ8MAsG9f0ymZryRspYvYHbGm6aS7BzlmCdgW6SJ8EQQOhtxM+PW5Fvf9wJAH0AZquVB6gUW/L3KMRJhNoNwj5mPRCXZBEDRABKADtECcKIoNlppqTVtowRyi0QjnD0qptPf9CCdqJB1T94LRc+Dqe8HDt3Un5YBcvHjR3iI4DLbQxfbT21mTuQY3lRsvj3/ZYaeTWqQLFw+pYmHcJKlQVd8wGDKt2d1UTioWX7uYO9fcyaacTXx/5Hvu7H+n+UJbCeUeMR9Le18TRFEMARAEIR2IB2ZYoG01leVQXgT643Dxb7iYCad3wfHfobRGRStnD+h/AwybJf1vp3w3CvKm3FDOq3+8CkDksEh6+bQsDNSh6TwYrn8V/vcvKSij+0jo2L3Z3bp5deOFa15gwW8LiNkRw+guo5XwVgenqKywzftazDgIgqAFTGv0RVHUC4IQZm7bKk6WHef2+IEEGAz4Gwz0L69gXGkpA8srTHNjxW6B5HYaRUXfm/AZdjP+vr4IDpALydIMHz7c3iI4DNbWxad7PiW7IBtNRw0PDHnAqscyl1bpYvQcaeHnkfWweo4UlNGCSL2be99MSk4Kv2b/yovbXmT5jcsdciSl3CMSi766q837WnLkoAHqTgvlCoKgFUVR15a2giBEApEA7sHuZLm6kIWUgO9/wHuo8ah0xre4MwW5EzmdfxXkC/A38Gsa7s4CnT0FgrycuGZgDwJcK/Aqz8XHVaBfv364ubmxd+9eAAIDA+nfv78pPtrNzY2xY8eSnp5uWmMxZswYTpw4wcmTJwEYMGAAKpWK/fv3A9ClSxd69+5NWloaIFWvGzNmDNu3bzf5TMaOHcvRo0c5c+YMAIMHD8ZgMHDo0CEAgoKC6N69O9u3bwfAy8uLkSNHkpaWZipg0qtXL06fPm2qkzt06FDKyspMOex79OhB586dTfllfHx80Gq1pKammqbnJkyYwL59+0zD7+HDh1NYWEhWVhYAwcHB+Pn5odNJP4evry/Dhw8nJSUFURQRBIGJEyeya9cu8vKkEZtWqyU3N5fs7GzpR9Zo8Pb2ZteuXQD4+/szZMgQtm7dCoCzszOhoaHodDoKCgoAGDlyJGfPniUnJweg2d/JaDRy3XXXWeV32rZ/G3G74wB4fODjXDh3oVW/U2hoKIcPH7bZ71ReXk7nzp1b/Dud7HQ3I4/9idvx3ynd8DLbPa5r0e80yTCJP93+RHdOx6KfFzHZZ7LD3U979uzB1dXVIX8nsM39dGxPAlu8pH3bgsWilS4/yMNFUZxR42+ZwIy6xqE1bavo008jfp2UwPmyAs6WnGPPxV3svLCd3LLqQuKdnIfQyXg9xfp+5OSWkl/acNnDILUHw3t0JKSXH2N6+zGoqw8qJ/mMMLZs2cKkSZPsLYZDYC1diKLInKQ5bD+9nWl9p/HK+FcsfgxL0yZdZG2BL6ZJayH+sQZ6X9ui3TYf38y8zfNwU7mReGsiwR2DWyuuVbnS75HUn1eQnL2I1T5e7L1/b5uilSw5csgF1HX+1li8W2vaAuDs5MI1wSE1/jILURQ5mn+U7498T+KRRM5X7OM8+xjZbyTvj/03fi7dyTxfxKEzRRw+W8j+UwXsPZXPSX0pJ/WlrNsjvW14uzkzurcf1/YLYEL/TvQO6NAup6MUWs767PVsP70dtZuap0Oetrc41kMzSYrg++1N+H4OPLxNWhPRDJN7Tua2PrexJnMNL257kc9v/NxudSwUarM9ZR0Bu6L5sXsnnGj7c8ySIwctEF/lZL78tzxRFOuFB7WmbRXDhg0Td+9uvABJYXkhqw+vZvm+5eReysXZyZn/G/p/zBk2BzeVm6mdwSiSeb6Incf1/Jmdy59HczmeW1Krrx5+HoQN6sz1g7swKtgXZ5VjzalmZ2cTHBxsbzEcAmvooqSihNt+vI2zJWd5aexLRPSPsGj/1qLNujBUSgtDc/6AfjfA7O+gBWsE8svyueOnOzhXeo6nQ57mgaGO45O5Uu8Rne5Pev80nVc6u5PUwZM7+93JovGL2jRysOgiuMuLLaoikDRATNXU0eXvuVXhqk21bYirr75a/Ouv5ld15pfl807GO6w+shqAvuq+vDPpnSaHvaf0paT+fYGth8+T+vcF9CXV01FqTxeuH9yZW4d3Y6zG3yEMRUFBgZJ18jLW0MX7uveJ3xPPYP/BfDP1G9m8EZulC30OfBIKl/QQ/jKMf6JFu/124jce2fiIw00vXYn3yP7Dh/D5+mby3PK5O6gLbio31k5fSxevLg6xCG6OIAjzL0ceRQBzamyLAWa2sG09SkpKmtpsoqNbRxaOW8iKG1cQ7BPM3/q/mb12NpuOb2p0n25qD2aO7MGHd2vJeDGc1XPHEjVRg6ZTB/QlFaxKP8F9n/7JmP9u5D8/7WX3Cb1dFwFVObUULK+LnIIcVuxbAcCC0QtkYxjATF2oe8C0JdLn5EVwLK1Fu13b/Vpu63MbZYYyXvr9JYyise0yWJAr7R7JPHEap29mEiSc5/WAbgDcM+geOndoe0U8ixoHURR1oijGiqKYfPl/fY1tM0RRjGtJW0sQ0jmE7275jvBe4RRVFPHE5id4X/d+sxevykkgpJcfC24axKZnJpH89ASemNIPTUAHLhaX80XaMW77cBvXv7OVpSmZXCgqs6TYCnYmdkcsFcYKbutzGyMCR9hbHNsycKqUVkY0QOIDUHyhRbvNHzUff3d/dOd0rDzUspoRCpbj1MV8Ln46k4Fk84tXELvdjHi7evPg0AfN6tf+cyQtxNm59b7zDi4deGviWzwd8jROghPxe+J5IfUFKowNRzE1RN9Ab54K78/GZybyy+OhPDi+N34dXDlyrojF/zvI2MUbefzbv0jLvGiz0YSvr/xXeVsKS+oi9WQqW05soYNLB57UPmmxfm2FRXQx5SXocQ0UnpYc1C1I793RrSMvXvMiAO9kvMPJopPmy2EmV8o9crGwlP1L7mW0uJtcJzXLNVKp2sirIuno1tGsvtt/4r3L/H7qd57c/CSllaVM6jGJNye+WctR3RoqDEY2HzzHqvQcNh08h/GyCvt39uL+cb2ZfnUQHq7ymY5QgApjBXeuuZOj+Ucdzrlqc/JPwtJroeQiTIyGyc+3aLdntjzDhmMbGNt1LEvDlyoRf1amsLScTe89xO2XfqIED365YQGvHF5Glw5d+GX6L6bnW7tPvFdY2PZl4ADjuo1j2fXL8HH1YUvOFh5JfoTiiuI29eWicuL6IV1Y9s9RpEZfx7wp/Qj0duPw2SKe/2EP1yzeyOv/O8i5gktmydwYKSkpVulXjlhKF4mHEzmaf5Qe3j24Z9A9FunT1ljsuugYBHcuAwRIiYFDv7ZotwVjFtDRrSNpp9NYk7nGMrK0kfZ+j1yqMPDzx89x+6WfqMCZgunLWHZyLQCPjXiszS++NZGNcbAEwzoNY8WNK+jk0Yk/z/zJoxsf5VKleQ/wbmoPng7vT2r0dbx31wiG91CTX1rBJymZhMZsJjpxN3+fs2wVO7mM9myBJXRRUF7Axzs/BuDpkKdxVbma3ac9sOh10ec6mPJv6fP3kVIOs2YI8Ahg/qj5QHUWW3vRnu+RSoORb+Ne5+7CTwHIv+ED1jtf5HTxafr59uMWzS0WOc4VZRwA+vn24/MbPyfQI5CMsxk8k/JMq3wQjeHq7MTtI4L46dHxfP/IOG4a2oUKo5GV6TmEv5PC3K8y2Hcq3wJngDJcr4EldBG/Ox59mZ6RnUcypecUC0hlHyx+XYQ+DQNvgbJ8WHkvlDc/0r5Vcyvjuo2joLyA1/983bLytIL2eo8YjSKfr1jCfefeBODs+EW4htxC/J54AJ7UPmmxCLsrxudQl0x9Jv/89Z/kl+UztfdUFl+72OIJxLLOFxH/21FWZ5yg3CBFSYUNCmTelH4M6153gbiCPcgpyOG2n27DYDTw3S3fMdh/sL1FciwuFUD8dXDxCAyZDhHLpVQbTXCy6CTTf5pOaWUp709+n8k9J9tI2PaNKIqs+OZr7j78JG5CBaeHP0bX6a/xdsbbLN+7nJGdR/LZDZ/VM4zt3ufQ6mI/zdBH3YdPwj7B09mTdUfX8fqfr1t8KKrp5MXiO65i6/zJPDi+N+4uTiQfOMdtH24j8ot0DpwuaFO/VYm3FMzXxTu6d6g0VnJrn1tlbxiscl24+8BdX4OrF+z7QUqz0QxBXkHMu3oeAK/+8SqF5eb5C9tCe7xHvvnpFyIOPysZhn5303Xaq5wqOsXX+78G4JmRz1h0xCQb42CNguFDA4bywXUf4OrkyrcHv+Wbg99Y/BgAXTq6859bB5MafR1REzS4uzixYf9ZbnrvNx77RsfRC61zjFdlbVQwTxc7z+0k6VgS7ip308NMzljtuug0oNpBvelVOPBLs7vMHjibYZ2Gca70HO9mvGsduZqgvd0jq9dv4oa/HsFbKOV096l0nf0hCALv//U+5cZypvaeytCAoRY9pmyMg7UY3XW0KeNm7I5YqxZQD/ByY8HUQWydP5n7xwXjqnLil92nCX87hRd/3MO5QutENynURxRF3kp/C4B/DPmHWStJrwgG3ART/iN9/j4SzuxtsrnKScWisYtwdnJm1eFVpJ+x3JTwlcbarWmM+/3/CBAKON1pPF3v/xycVOy7sI+1WWtxcXJhntbyLzeyMQ6enp5W63uqZipRw6Iwikb+lfIvMvXNR2aYQ6C3OwtvG8KWf01i1sgeGEWRr/44zsTYLbyTdJiS8qZHSVqt1qryyYm26mLT8U3sPL8TP3c/s1eSOgpWvy5Cn4KrZkJFMXw7G4rONdm8r29f5lwlZcVZlLaIMoPtsgm0l3tk846dXLXxProKuZxRa+k6JxGcXaWXmwzp5ebeQfcS5BVk8WPLxjgYDM2v1DSHR0Y8Ykq18djGx9Bfsmg2jwbppvYgJmIY65+cQPjgzpRWGHhv4xEmvbGFVTtyMBgb9oHk5uY2+PcrkbboosJYwTu6dwCYO3wuHVw6WFosu2D160IQ4LYPICgE8o/Dt3dBedM5zx666iE0HTVkF2SzdNdS68pXg/Zwj6TtPkjPX2bTUzjPGa9BdHn4J3CVXpJTTqSw48wOOrp15KFhD1nl+LIxDlVVm6yFk+DEa6GvMdh/MCeKTrAgdYHNkoj16+xN/D9GsipqLMO7d+RcYRnzV+/m1g9S2Z5Vv1B6VWUohbbpYvXh1RwrOEYvn17c2f9OywtlJ2xyXbi4Sym91T3hZEazKTZcVa4sGrcIAYHle5dzKPeQ9WVE/vfIXwf/xm/1DPoIpzjrrqHz3LVScADSy03VlOjc4XPxcbVO9lnZGAdb4OHswbuT3qWjW0dST6YSvzvepscf3duPHx4Zz3t3jaBbR3f2ny5gVtwfPPqNjhN5LctKq9A0xRXFLNklZR99UvskLk4udpZIhngFwj2J4N4RDv4CG15ssvmIwBHcNfAuKsVK/vP7f6g0Wj64pD2x/+9jeHx7BwOE45x160ngo/9DqFGA6buD35FdkE2wTzAz+89soifzkI1xcHMzfzl4S+jq1ZXXr30dAYGPdn5E2qmWpS62FE5OArePCGLjM5N4Kqw/7i5OrN19milvpfDBxiNcqjCg0WhsKpMj01pdrNi3gtxLuQzvNFzWC94awqbXRacBcNc34OQCf3wMv3/YZPMntE/QtUNX9l/cz5f7v7S6eHK9R/4+ngNf3c5A4RhnXXoQ8MgGBO8upu15l/JMLzfPjnwWF5X1Xm5kYxxUKtslsgsNCiVqeBQiItFbozlTfMZmx67Cw1XFE2H92PjMJG4Z1pWySiNvJR3mhne3svOc8uZVhbe3d4vbXii9wOf7PgekNBntbRVta3RhEYJDq2tAbHgBdjYeCt7BpQP/GStFO3208yOy87OtKprNdWEBjufkULH8dgZzlLPO3fCd+yuqjl1rtflo50cUlhcytutYJnSfYFV5ZGMcWlrsx1I8POxhxnYdS15ZHtFbozG0IHWxNQhSe/Dh3Vq+mTOGfoFeHLtYwrzEA0R9mc4pvWUXBsqR1ix2WrprqSkrr7Zz+4hmqYldFn4NmwE3Xk6T8dNjcHBto01Dg0JNhYEWpi20qk9PbovgTp86TtlntzBIzOSMqisd567H1a97rTZH8o6QcDgBlaBi/qj5Vn+5kY1xsDUqJxWvT3idAI8AdOd0LNuzzK7yjOsTwLonruWFqYNwV8H6fWcJezuF+K1ZVBgco/qWI3O84DiJhxNxEpx44uqWlcBUaCHXzIUJ86UiQQkPwNHfGm36r5H/ws/dj4yzGSQeTrShkI7LhVPHKIu/iX5iNidUPfCem4S7f89abURR5I0db2AUjczoP4O+vn2tLpdsjENbiv2Yi5+7H6+FvgbAkl1L2Hlup81lqImLyok5EzQsua0bNw3tQkm5gdfWHeDWD1LZmWP90FtHxN/fv/lGwAd/fUClWMntfZRBKbIAACAASURBVG63yY1lD1qqC6sw+XkY9RAYyuCbWY2WGVW7q3lhzAsAvJX+FqeKTllFHLvqohXoT/1N2bIbCBZPkK3qhc/c9XQI6FGv3abjm0g7nYa3qzePjHjEJrLJxjh4eHjY5bjjuo3j/iH3YxANPPfbc3bJE1OXiaOGs+TeEJbfP4ruvh4cPFPI9I+3sXDNPorKrix/xJAhQ5pts+/CPn7N/hVXJ1eb3Vj2oCW6sBqCADe9AcPukhbJfR0Bx7c32DS8VzjhvcIpqSzhpd9fskp6bbvqooUUHttNZXw4QcbTHHHS0PHh9fgE1F/MVlpZSsyOGAAev/pxfN1tU+VONsbB3GI/5jDv6nkM8hvEyaKTvLb9NbvJUcXWrVsBmDwwkKSnJhI1UYOTILDi92zC305h08GzdpbQdlTpoimqFrzdPehuunTo0kxr+dISXVgVJyeY9jFcNQPKi+CrOyFnR71mgiDw/JjnUbup+eP0H6w+striothdF81Q/Pc2WHETAWIuO52GoJ67Ht9OXRtsG787ntPFpxnoN9Cqoat1kY1xsCcuKhdiJ8Ti4ezB2qy1rM9eb2+RTHi4qlhw0yDWPDaeYd07cjr/Eg+uSGfet39xsch26QoclbRTaWw/vR1vV28euso6K0kVauCkgmmfwNA7obwQvroDjv1er1mARwDPj5HKj76Z/iani07bWlK7cWnfWpy/mo63WMRWp9EEPrKWTp0CG2x7rOAYK/atAOCFMS9YrFZDS5CNcbB32GFwx2CeCXkGkNIQ27PKVUP+lyHdOvLDI+N58eZBuLs4sWbXKcLeTuGnnSfbdVWspnxRRtHIuzopI+iDQx80u+C6o2MPv1yDqJxhehwMuQPKCuDL6XB4Q71mNwbfSFjPMIorii0+veQwuqhD+bYluCbcixtl/OJ0Hb0fWU23gIaniURRZPH2xVQYK7i9z+2MCBxhU1llYxy8vLzsLQIzB8xkfLfx6Mv0LPx9od0euqGhoQ3+XeUk8NC1GjY8OZFxffzJK6ngie92MueLDKvVs7Y3jekCYMOxDey/uJ9Aj0DZ1oVuDU3pwuaonKU039p/QuUl+G427KkdnSQIAi9c84Kp7nTC4QSLHd6hdAFgNFC59l+4Jj2HE0aWOs1i8MNf0COg8dQXG49vZNupbXi7ePNUyFM2FFZCNsbB1uscGkIQBBaNW4S3qzcpJ1L48e8f7SKHTqdrcntPf0++fmgMr99xFd5uziQfkMJeEzNOtLtRRGO6qDBW8IHuAwAeHvEwHs72CWiwJc1dFzbHSQW3vgfjnwBjJax+CP5YAjWuwQCPAF4cI6XfeDP9TY4VHLPIoR1KF6V6DN/chfOOOMpFFf92mseUh99GE9j4Qr3C8kL+u/2/ADx29WP4e9g++ko2xsHaWVlbSucOnU1zpTE7YqwWitcUBQXNV5ATBIG7Rvdk/VMTmDSgEwWXKnk2YRcPfZ7O2XY0imhMFz8c+YHjhccJ9glmet/pNpbKPrTkurA5ggDhL0PYQkCEX5+Ddc+CoTqq7sbeN3JT75sorSzl+dTnLZJ7yWF0ce4gYvx1qP7egF7swCOq/3Bv5Hz6BjY9E/JuxrucLz3PsIBhzBowy0bC1kY2xsGRuLn3zYT3CrfKXKml6ab2YPn9o3hzxnC83Z3ZePAc4W+n8L2u/Y0iqiipKDHln3n86sdxdnLM+ecritCn4I5loHKFHcvgm5lwKd+0+YUxLxDoGcju87v5bO9ndhTUguxfg7hsCkJuJgeMPbnHKYanH3qQAV2aTu2hO6tj1eFVOAvOvDTuJZs6oWsiG+PQoYPj5NwXBIEXxrxg1VC8phg5snW1wgVBICKkO0lPTWTy5VHE06t2EfllBucL5R3R1JAuvj7wNRdKLzDUfyjhvcLtIJV9aO11YXOGzYB//gKeAZC5EZaFwbmDAHR062iqyLhk5xL2X9xv1qHsqouKS7BuPqy6D6G8iJ8M43jA6TViHrqNwd2aTq9dbihnYdpCAB4Y+gD9ffvbQOCGkY1xqKiosLcItfD38GfB6AWANFdqy+R8Z8+2bR1Dl47ufHb/KN6IGIa3mzNJ+89y/TsprNsj3zDCurrQX9Kb3jyfCnnK7lFutqSt14VN6TkG5myEToPgwmGInwy7vgOkBaezB86mUqwkems0JRVt9zPaTRfnDsKyKfDnUipx5pWKe/mP6kmWzZnI0KDmo+WW7VnG0fyjBPsEEzU8ygYCN45sjEN5ebm9RajHTb1vYnKPyRRXFLMwzXbRSzk5OW3eVxAEZozswfqnJhDaN4C8kgoe+VrHvG//Ql/ieDpujrq6iN8TT1FFEeO7jWd019F2kso+mHNd2BTfYMlADLsLKkrghyhY8ziUF/N0yNP0VfcluyCbxX8ubvMhbK4LoxG2L4W4SXB2L2edg5hetpBE19v5es41LTIM+y7uM9WQ+c/Y/+Cmsk2ZgsaQjXFwRARB4N/X/BsfVx+2ndxmt+ilttBN7cGX/zeaV6YNxcNFxZpdp7jh3a1sOdR0XeD/b+/M46Kq9///PMOw7yCbiCCoiGAgooiSS2luuZVaN22xwrLFa96blS16rW+ldjPbXW519VaaZYv9ShMT0MQFSVNQQUBWAdn3ZWbO748ZFjfWgZmB83w8eCjMOWfe8zrLez7vz/vzfuszORU5fH3hawCWj1iuY2skWsTEEuZ+qs5mMjKF+O3wyVjMsk6xYdwGTI1M+eHSD/yS+ouuLW2dq0nw+TT4dSUoqjlqOYk7KtaSYebLl4+Htskx1CprefnwyyhEBQ8MeYCRriO7wfCWMRjnYGZmpmsTboqThRMvjHoBgA1xG7hadbXL33PQoEFaOY4gCDw42pNf/347IzztySur5ZHPT7Lq+7NUGkiNpuZafHT6I+pV9czwnsEQhyE6tEo3aOu66DYEAUY8oh5FOPtDcRp8MYOBx7ayMljt3NceW0tmeftHAd2iRX01RK+HT8Mh8xiipQvv91nNA4WPYmJhw1cRbXMMAB/Ef0BKaQpeNl5688XGYJyDPseOZ3rPJNw9nPK68k4NhduKtrviefWx5Jsnwnhx2hBMjGR8dTyD6e8f5lR6sVbfpyto0CKpOIm9KXuRy+Q8E/SMjq3SDd3VLVHruA6DJVHqst+CDI5/yvxf/4/Jdn5U1leyMnoldcr2hTy7VAtRhHPfwYcj4dD/gbIWReBCnrD9mHezfHGwNOGriNH4922bY4jLjWN74vbGPvb6sibHYJxDdbX+NrZpCC9ZyC04kH6AA+kHuvT9zp07p/VjGskEnhzvw4/PjGWIqzXphVXM//QoG/ZfoE6hv/0iGrTYeGojIiL3+d5HP+t+rezVM+mK66LbkJvAHS9DxO/gFohQnsPqM7/RV5RxrvAcb594u12H6xItRBFSDsFnU+DbR6E0E1wCqH7gBx7If5DfUmvpY2XC1xGj8XNrOSupgfK6cl754xVERB4f9ji3Od2mfbs7iME4B32nr1XfxuHgm8ffpLS2tJU99BM/Nxt+fGYsT4z3RgQ+OpTC3I//IDlP96XKb8XxK8c5kn0ES2NLlty2RNfmSHSGvkEQEQWzPsTW3Il3c3IwUYnsTtrN98fWX7O6utsQRUjar0693TEHMo+DpRPMfJ/Shw7yQKQJJ9KKcLExZeeSsFbXMTQdVmT10dVkV2Tj5+DHk7c92cUfpH0YjHMwNu66Rtra4j7f+xjuPJyC6gLeiXuny97H2fnmFRy1halcXel115Iw+tmbk5BTxt0fHOHzP9JQqfRr4Vwfpz68e+pdAB4LeAwHMwcdW6Q7uvq66DZkMgh+EJ49hf+oZ3i5tBKAN85vJ/GzO9Spr3Utp7lqRYvKAvjjffhghHrRXnYcWDjCna/Bs/EUDfkbCz87yZ8ZJbjbmfPNE2Gtrnxuzs6LOzmQfgBLY0s2jN+AsZF+PeMEQ1klO2LECPHUqVO6NqNV0krTmPfTPOpUdWyZvIWwvmFafw+FQtFtVSfLa+p5/edEvonLAmDsQEfemR+Im61+xEX3XtrLqj9W4WzuzM/3/Kw38Vpd0J3XRbdSXcyaXx7lu4pLuNcr+DonF3u5JQTcA0Nng+dYML42YaXDWlQWqEcJF39R/6vSrK+y7gthT0HIo2BiSV5ZDYu2HSc5v4L+DhZ8FRFKP3uLNr9NQmECD/7yIPWqejaM38BUr6ntt7WNCIJwShTFdq8K1JpzEATBG5gHxAPBwBZRFG/au1IQhHXASqAEiAOeEEUxtaXj+/r6ihcvXtSKrV3N1r+28v6f79PPqh97Zu/R+gMrKiqKCRMmaPWYrbE/IZeX9pylqLIOGzM5b8wdxqzAvt1qw/XUKeu4a+ddFCoK+deYf3HPoHt0ao+u0cV10V3UKmt5+JcHSSg6T7DKmK3pKZg0vCg3B69wGHA7uASA6zCi4hJb10KlgtIMyI6H7FOQeUI9OhA1c2yCDAbdpc6oGjhZXWkWyCyqYuG242QUVTHI2Yr/PR6Ki03bsynL6spYsHcB2RXZ3Od7H6+MfqW9crSLjjoHbX7N2C2K4giNMXHAVmD+LbZNEUVRf9OPOskjAY+w7/I+koqT+Pj0x/wj5B+6NqnTTPF3ZXh/O1787iy/X8hn2dd/EpmYx+uzA7C10M1weNfFXRQqChloN5BZPrN0YoNE92BqZMqmOz7ggV8eIL4qn9fGLuQtwQXhUiTk/gWXDqh/NIwxtoHEfuowkIWjukKsSqEu+FdTop5MLrvSNDJoQGYM3hNhyHTwnQ42134BupRfzqJtJ8gtq2GYuy3/fXQUDpYmtBWFSsHKmJWN8wwrR67slC5diVZGDoIgBAPrRFGc3OxvxaIo3rSLhSAIS0RR3NKG4y4BlgA4OzuP2LVrFwDe3t5YW1tz5swZQN1M3N/fv7E1oFwuJzw8nPj4+MbqjCEhIeTl5TWunBw0aBCmpqaNWQ3Ozs4MHjyYI0eOAOpUuLCwMOLi4qioqAAgNDSUrKwssrOzAfD19cXIyIjERHUdGFdXVwYMGEBsbCzpten8O/ff6paInqtwEV0ACAsLIy0tjdxcdbmNoUOHolQqaRgVubu7069fP44fV/fftbKyIiQkhNjYWGpr1XWQjI2Nsbe3Jz9fvWAtICCA2tpakpOTAfDw8MDFxYW4uDgAbGxsCA4O5siRIygU6vUL48aNIyEhgcLCQgACAwMpLy8nNVU9gPPy8sLBwaGx9LG9vT2BgYFERUURlVnPVxfqqFOCo4URj/kbM9TRiODgYIqKirh8+XKXn6f90ftZm7OWKlUVH97xIZZ5lh06T6DuTx4aGsrx48cbs+K0cZ7Cw8NJSkrqtvNUWVlJv379CAwMJDo6GlEUEQSB8ePHc+bMGYqL1anJ3XmetHU/NZwn20G2LPp5EbViLdNsp/H6tNfJSDyJ4uJv2JQl4aTKR154AVm9ep6iNeqM7ZB7jqTA1JNcWV9Kbf3wCxp10/P0zcHjvBtXQ0U9jPSyJ8JXiYmgbNd5WhW5iujyaKyNrPlm9jekxKd0+XmytrbWXVhJEIR5wH2iKM5v9rcUYL4oijcUVhcEYSWQijqsNBl461YhqAZCQkLEhpvIUNhwcgPbE7czxGEIX834CmOZfk04dYa0gkqe23Wa05nq0/ZY+ACen+KLmXH3VJBcf3I9OxJ3EOoayta7tur1OhgJ7RKdGc2yQ8tQiSreGPsGswfOvnYDlQoq8qCqQD2HUKV+YCMzApkcTK3B1kM9KjBuW8j36KUCIrbHUVmnZKKvEx8vHIG5Sfuu9a8vfM2bx9/EWGbMtru2EewS3K79O0pHw0raylZqb4rIFlEUvxVFMRLYBRxsbQd9aPbTXp4Oehp3K3cuFF3gvwn/1dpx9cFJDuhjybdPhrFi8mCMZAL/OZLGzA+OkJDT9Sm8GWUZfH3hawQEpltMlxyDBn24LrqD8R7jG8Mxrx197cae7jIZcUnZ6sV1PhNh2Dz1j/9c8JsJ3hPA0afNjmHfuVwe+fwklXVKZgf1ZctDIe12DEezj7LuxDoA/jXmX93mGDpDi85BEIQlgiCsa+FnkmbTIsDuut1v6TCajxI0I4tgQRCu3/8a9KXZT3uwMLZgddhqQF2GOK00TSvHbRiW6xq5kYxldw5iz9IxeDtZkpxfwZyP/uCTqBSUXZjyuvHURhQqBbMHzsZB2XtTV69HX66L7mCh30KeDHwSlajixZgXic6MvuZ1bWnxv2PpPPXlKeqUKh4O82TjgiCMjdr3nfrP/D9ZHrUcpagkYlgEM31masW2rqbFTymK4hZRFF9o4SdSs2kqN3EGtwgpBQuCcENOamthJUMlrG8YcwbOoU5Vx5qja1CJ+rvauKMEetjx/569nYfCPKlXiqzbd4H7t8SSUaj90d6pvFNEZkRiLjfn2eHPav34EobDU4FP8Yj/IyhEBc9FPcfRnKNaO7Yoirx7IIlXfjiHSoR/TB7Mmln+yGTtG6WevXqWpZFLqVZUM9N7Js8MN5zSLloJK13vBDRprZHNf282MkgFNjd7bRJwbefxm6BPzX7ayz9D/omjmSPx+fHsvtj5JuqhoaFasEq7mJsYsXZ2AF8sHomztSknLxczbVMMO09kaK2UuUpU8c5J9eLCxf6LcbZw1kstdEVv00IQBFaMWMH9vvdTr6pn2e/LOJihjlB3RguFUsWq78/y/sFkZAK8fc8wnr1zULvDl+cLz/NE5BNU1lcy1Wsqa8euRSYYzLpjra6QjhAEYaXmYT8PiGj22jpgATSOEFI1IaslqCekI2442nXoW7Of9mBratvYd3pj/MZONwbKysrShlldwgRfZ/YvH8eMYW5U1il5cc9ZIrbHkV/e+b7VP6X8xLnCczibO/Ow/8OAfmvR3fRGLQRB4KXQl1gweAG1ylpWRK1gd9LuDmtRWasgYnscX5/IxFQu49NFI7h/VP92H+d0/mmWHFhCeV05d/a/kzdvf9Pg2tVqzTmIohgviuJ6URQjNf82n1eY3zx1VbPNlmZhq1ZDSvrY7Kc9TPaczJ3976SyvpK1sWs79W26IfVPX7G3NOHDB4az6f4gbMzkRJ7PZ8rGGH7tRMe5iroK3jv1HgDPhTyHhbF6Naq+a9Gd9FYtZIKMV0a/wlNBT6ESVayNXcvWhK3tvsfyy2u4b0sshy5exd7CmK8iQrnL37Xd9uy/vJ/H9j9GSW0JE/pNYMO4DQaZqWg4YxwDp6HvtLWJNYezD/Nz6s+6NqlLEQSB2UHu7Fve1HFu6ZfxPLfrNKXV7R8FbvlrC4U1hQQ5BTFjwIwusFjCkBEEgaWBS3kt7DVkgoxfS39lRdQKyurK2rR/Ul45cz86yrnsMjwdLdjz1FhGeLYv2UEURT4/9zn/jP4ndao65g+ez8aJG/WuZlJbMRjnoK/NftqDk4VTYwre2yfepqC6oEPH8fX11aZZXUpfO3O2PzqKtbP9MTOW8f2f2UzZGEN0UtubIl0uvcyO8zsQEHhx1IvXxH4NSYuuRtIC5g+ez6aJm7CQWxCZEcmCvQs4e/Vsi/vEJF3l3o+Pkl1STZCHHd8tHcOAPu2b4yytLWVlzMrGIpArRqzg1dGvGlwoqTkG4xx6Si77bJ/ZjHUfS1ldGW8ce6ND4SUjo+5ZaKYtZDKBh8K8+GXZ7Qzvb0duWQ0Pf3aCVd+fpaINHec2xG1AoVIwZ+Ac/Pv4X/OaoWnRlUhaqJngMYHNYzfj7+hPdkU2D+17iE/PfEqtsvaGbb88ns7iL05SXqtgxjA3di4ZTR+r9jUKOpJ9hLk/zmXf5X2Yy815Z/w7LA5YbPDPLINxDvrc7Kc9CILA6tGrsTS25GDGQfan7299p+toKC9gaHg7WbH7iTBWTvXF2Ejgq+MZTH0vhqMptx5BxWTFEJMVg5WxFcuCl93wuqFq0RVIWjRRcrmE7dO2s9BvIQqVgo9Of8TsH2ZzMP0goiiiUKpYuzeRl78/h1Il8vREHz742/B2rfDPqcjh1T9eZWnkUq5WX2W483C+m/kdU7ymdOEn6z4Md8xjwLhZubFixApeP/Y6bx57k5EuI3E0d9S1Wd2C3EjGUxMGcscQZ/7xzRkScsp4YOtxHg7z5IVpQ7AwabokqxXVvHn8TQCeDHySPuZ9dGW2hAFiYmTCi6Ne5A6PO3jrxFtcKrnE8qjlBPUJpiRvFGeTPDA2kvPm3GHMD/Fo83HzKvPYenYr3yV/h0KlQC6T83TQ0yz2X4yRrOeM3gxm5GAIzX7aw/zB8wl1C6W4tpjXj73ervCSq2v7Myj0jSGuNvzw9FiemzQYuUzgv7HpTH3vMLEphY3bbP1rK9kV2Qy2H8xCv4U3PU5P0EJbSFo00VyLUW6j2D1zN6tCV2FlbM3pgnguG32K9aB1zJt0Dm+PfGoULadaF1QX8M3Fb4j4LYIp301h18VdKFVKZnjP4PtZ3/P4sMd7lGMAqdmPTsmpyOGen+6hsr6SN8PfbPOy+traWsNtJn8TEnJK+efuvzh/RZ1Z8uBoT+4fa8qifQtQqBTsmLaDIOegm+7b07ToDJIWTdxMi8jEPJ7bHUuN2QmsnE6gMMprfE0ukzPUcShulm5YyC2wMLagRlFDRnkGmeWZ5FXmIaJ+VhoJRkzynMTSwKX42Pl06+fqCDpv9tPVGFKzn/bwffL3vHb0NayNrdkzew+ulq1/++uJTV3qFCo+jrrEh79fQqFSYef9GUrTZO4ddC9rxqy55X49UYuOImnRRHMtVCqRTQeT2XRQXYJ7+jBXNsy7jcTi0+y/vJ8/8/8kuTi58eF/M4xlxozpO4ZJnpOY6DERW1Pb7vgYWkEfmv1IdIA5A+dwMOMg0VnRrDm6hk8mfWLwWQ4dwUQuY/mkwUzxd2XpD1spNE1GpbCgKGsSJVV12Fm0vaGKhEQDpdX1rNh1moMX8pEJ8PyUITw53htBEBjpOpKRriMBKK8r51zBOUpqS6iqr6JKUYVcJqe/dX/6W/fH1crVIBeydQaDcQ4ymcFMj7QLQRBYHbaauT/N5Y+cP9h1cRf3D7m/xX3MzXtun2QXewVK+x+hFpSFM/gpuZSjSTGsne3PtADXGxxnT9aivUhaNGFubs7ZrFKe+uoUmUXV2Job88HfhjNusNNNt7c2se6Sfu+GjME8cQ258F5rOFk4NfaRfSfuHVJLWmyn3aMLrL19/G1KaosJdQvl50eeY5SXAwUVtTz1ZTwR20+RU3JtSnNP1qK9SFqoEUWRZNGVez85SmZRNcPcbfn52fBbOgaJm2MwzqGysm1t/wyVqV5TmeUzi1plLStjVlKnvHUtqYbWlD2N3zN+59fLv2IuN2dN2BoGOluzc8loXp8TgLWpnMjzeUx+N5ov/khr7BfRU7XoCJIWUF5Tz993nuaVH85Rp1SxaHR/dj8ZhoeDha5NMzgMxjmoVD2vD8L1rApdhYe1BxeLL/J+/Pu33K6nLAhsTsOKcYC/B/+dftb9APXq6gdHe3JgxXim+rtSWadkzd5E5nz0B2ezSnukFh2lt2txJrOEGe8f4aczOZgawab7g3hjzrBua13b0zAY59AbsDS25O3b38ZIMOK/if/VavMSfWfDyQ1crb5KkFMQfxvytxted7U149MHR7D5wRG42ZpxNruU2R8d4cvztZTVGG45d4nOo1KJbIlJ4d5PjpJRVMVQNxv+Ncac2UHuujbNoDGYVNaeuM7hVmw+s5kPT3+Io5kju2fuxsni2lhpT8tnP5hxkOWHlmMiM2H3rN1423q3uH1lrYKNB5L4/OhllCoRJ2tTXp7ux+ygvr0y06uBnnZdtIWckmr+8c0ZYlPViycfGePFS9OHgFLR67S4FR1NZTWYkUNt7Y1Fs3oqjw97nFGuoyisKWRlzEoUqmuL06WlaacXtT5wteoqa46uAWD5iOWtOgYAS1M5r9w9lJ+eGYufkylXy2tZvus09205xoXctpVo7on0pOuiLfx4Opsp78UQm1qIo6UJ/3k4hDWz/DGVG/U6LboCg3EOhtwJrr0YyYxYN24dfcz7EJcXx0enP7rm9dzcznWS0xdUoopX/niFktoSxvQdc8sSGbfCv68tzwcbsX7ebThYmnAirYjpmw7z2o/nKKky7OZQHaGnXBetoc5eO8Xfd56mvEbBJD9n9j83jjv9XBq36S1adCUG4xx6G33M+7B+3HpkgoxtZ7cRkxWja5O0zlfnv+JozlHsTO14Y+wbHeqvKxMEFoR4cOgfE3gozBOA7bHpTHgnih2xl1Eoe34iQ2/i579yuGtjDL+czcXSxIi37hnG1odC2l1mW6J1DGbOISgoSDx9+rSuzeh2tp3dxqb4TdiY2LBzxk48bDzIz8/H2dlZ16Z1igtFF1j4/xZSp6rjvYnvcWf/Ozt0nOu1OH+ljH/tTeBYahEAg5ytWDXDj4m+hq1XW+gJ18WtyC2tYc1PCexLUI8Ixvg4su7e226ZotqTtWgvPX7OwVCcmLZ5NOBRJvSbQFldGc/+/iwVdRUolUpdm9UpSmtLee7Qc9Sp6pg3eF6HHQNwgxZ+bjZ8HTGajxcG4+FgTnJ+BYs/P8mD/zneWNivp2Lo18XNUKlEdsReZvK70exLUI8W3pgTwJePh7a4dqEnatHdGIxzqKlpuaRuT0UmyHjr9rfwsfUhpTSFFw6/wPkL53VtVodRiSpeOvwSWRVZ+Dn48cLIFzp1vJsVYxQEgenD3IhcMZ5V04dgbSbncHIB098/zIpdp8kqrurUe+orPa0w5bnsUuZ9epRXf0ygvFbBJD8XDqwYz6LRnq1mpfU0LXSBwTiH3oyViRUf3PEBtqa2xGTFsLdkr65N6jCbz2zmcPZhbE1t2ThxI2byrusNbio3Ysk4H6Kfn8gjY7yQywT2/JnNHe9Es3ZvIgUVvScDzpAoqarj1R/OMevDI8RnlOBsbcqni4LZ+tAI+tpJ9aO6C4NxDiYmvbsqp4eNmPRcWwAAEKFJREFUB++Ofxe5ICeyLJI9yXt0bVK7icmK4ZMznyAgsP729bhbdX6Rkrt768dwsDRhzSx/Dq6YwOygvtQpVXz2Rxrj1h9iw/4LlFb1jEy4tmihzyiUKnYcS+eOf0ez41g6giDwWPgADv5jPFMD3Nq1hsXQtdAHDGZCOjg4WIyPj9e1GTrnm4vf8Pqx15EJMt6d8G6n4vXdSUJhAov3LaZaUc2y4cuIuC1CK8etrq5udzXShJxS3v0tiYMX8gGwNpWzeKwXj4YPMOjS4B3RQh8QRZGoi1f5v1/Ocym/AoDQAQ6snR2Ar6t1h45pqFp0BT1+QrqnF95rKwt8FzDVdioqUcXK6JWczD2pa5NaJas8i6cjn6ZaUc3d3nfz+LDHtXbsjhSb8+9ry38eGcmep8YwdqAj5bUK3v/9EmPf/p11+y4YbLjJEAvvxWcUs3DbcRZ/cZJL+RV4OlrwycJgdi4Z3WHHAIaphb5hMP0cJJqYbjsdW1dbdl3cxbLfl/HZlM/wc/TTtVk3pbimmKWRSymsKWS022jWjlmrNyUugvvb8+Xjozl5uYj3DyZzOLmAT6JS+OxIGgtCPIi43Zv+jlI1z67gQm4Z//4tiQOJ6ladNmZylt05iAfDPDGVS4Xy9AGDcQ5GRtIF04C1tTUvBb9ESW0J+y/vJ+JABJsnb8bf0V/Xpl1DRV0Fz/z+DJfLLjPYfjAbJ2zE2Ei73bSsrKw6fYyRXg7seCyUPzOK+ejQJSLP57PjWDpfHk9n2jA3HgsfQHB/ey1Y27VoQ4uu5lx2KR/8nsz+BLVTMDc24tFwL5bc7oOthfauDUPQQt8xmDmHkJAQMS4uTtdm6BV1yjpWRK0gOisaK2MrPp70McOdh+vaLEC9lmFp5FLOFpzFzdKNHdN24GLp0vqOekBSXjmfRqfw0+kcFJq+EUEediwe68W0ADdM5AYTjdULRFHkRFoRm2NS+V0zz2Mil/HAqP48NdEHZ+uuy1iT6Picg8E4Bz8/P/H8ecPN79cmsbGxhIWpWxrWq+p56fBL7L+8H3O5OZsmbtJ5u8PimmKWHFjChaILuFu5s+2ubY39GbRNcy20TU5JNTuOpfPV8QxKq9UZTX2sTLhvpAd/G9Wffvb6FXLqSi06Qr1SxS9nr/CfI2n8lVUKqEcKi0b3J+J2b5xtus4p6JsWuqSjzsFgwkq9odlPW2leodZYZsy629dhZmTGjyk/8vTBp3l19KvMHTRXJ7blVuayNHIpl0ou4Wnjyba7tuFq6dpl79eV1Xr72pnzwtQhPHvHQPbEZ7MjNp2LeeV8dCiFj6NSCB/Yh/khHtw11EUvGsroS+Xi3NIavj6Rwc6TGeSVqW1ysDRh0WhPHg7zxLEb6iDpixaGjME4B4lbYyQzYu3YtVibWPO/8//jtaOvcb7oPM+PfB5jmXZj/C1xKu8UK6JWUFRThI+tD1vv2npDLwpDxMJEzqLRniwM7U9cejFfHkvnl7O5HE4u4HByATZmcu4O7MuswL6M8nJAJtOPCffupE6hIupiPt+eyuLghfzGNq4+TpY8Fu7NPcHueuFAJdqOwYSVelOzn9ZQKBTI5Tf363uS9/DGsTeoV9UT4hLChvEb6GPep0vtEUWRnRd3sv7EehSiglC3UN4Z9w52ZnZd+r7QshZdSXFlHT+dyWH3qUzOZTfVbHK1MePu29yYGuDK8P72GHWjo+huLVQqkT8zS9h7JoefzuRQVKkuky6XCUzxd2Xh6P6EeTvqJDtNV9eFPtLj5xwCAgLEc+fO6doMvSAxMZGhQ4fe8vXT+adZEbWCq9VXsTax5vmQ55kzcE6X3KRXq67y9om3+S39NwAeHvowy0csRy7rnhuzNS26g/NXyvjxdA57z+SQXdLUx7mPlSmThzpzxxAXxvg4YmnatZp0hxb1ShWn0ovZn5DLvnO5XCltqnk22MWKe4P7MTfYXeeTzPpwXegLPd45+Pr6ilIxLTVRUVFMmDChxW3yq/JZfXQ1R7KPADDabTSvjn6V/jb9tWKDSlSx++JuNsVvory+HHO5OavDVjPDe4ZWjt9W2qJFdyGKIvEZxfxyNpf9CblkFTc5CmMjgZFeDoQP6sNob0eGudtibKTdrKeu0iK7pJrYlEIOXcwnJukq5TVNnQn72poxbZgbc4e749/XRm/WsOjTdaFrevyEtET7cLZw5uM7P+bn1J9Zf3I9x64cY9YPs5g2YBqPBTzGQPuBHTquQqUgMj2SzxM+J7EwEYDb3W/n5dEva6VWkiEjCAIjPB0Y4enAKzP8SLxSxoHEPKKTrnI6s4SjKYUcTVH3OrYwMWKEpz3DPewI1PzoQ8MahVJFcn4FZzJLiM8oJja1kMyi6mu28Xay5M4hzkwb5kZQP7teOcfSG9D6yEEQhAOiKE5uZRtvYB4QDwQDW0RRLGlpn97a7OdmFBQU0KdP2+cRCqsLeS/+Pfam7EUpquvcj+83nkmekxjXbxwOZg4t7i+KImmlaRzKPMTOizvJrVQ3XHEyd+LFUS8y2XOyzr4xtlcLXVFcWceRSwXEphZyLLWQ1Ks3loNxsTHF19WGIa7WDHS2wsvREi9HC5ysTdukb3u0qFequFJSQ0ZRFcn55STllXMxt5zzV8qprr+2F4K1mZzQAQ6ED+zDxCHOeDpatu1D6xBDuS66A52HlQRBmAR4A5tFUWzxStYYO0LzfztgqyiK81vaJzAwUDxz5oxWbDV0srOzO1R1Mrsimy/OfcGe5D3UqdSThwICw5yGMdBuIK6WrrhZuiETZJTVllFaV0p6WTonc09SUF3QeBwvGy8W+S1ips9MLIx1m+vfUS10TX5ZDXHpxZzJLOF0Zglns0upqrt5gxozYxmuNmY425jhYmOGvYUxtubqHwsTOSZyGcZGAuWlJdja2aMSRRRKkap6JZW1CiprFZRU1XO1vJarFbXkldVwpbSmMaPoevo7WKhHM/1sCR3gyNC+Nt06sa4NDPW66Ap07hyaGSK25BwEQQgG1jUfXQiCUCyKYov1CaQ5hyY6G08tqC7gQPoBojOjOZF7gnpV6yWrHc0cGeU6irt97ibcPbxD/Z67gp4SW1apRDKKqriQq/4Gn1pQweXCKtILKynpgpLigqDOrPKwt8DH2QpfFysGu1ozxNUGB0vDrUzbQE+5LrSBIc05eAPXh5CKBEEIFkXxmprcgiAsAZYAODs7ExUVpT6AtzfW1tY0jCQcHR3x9/cnJiYGALlcTnh4OPHx8ZSVqdMMQ0JCyMvLIzMzE4BBgwZhampKQwaUs7MzgwcP5sgR9QSuqakpYWFhxMXFUVGhKSMcGkpWVhbZ2dkA+Pr6YmRkRGKiOvbu6urKgAEDiI2NBcDc3JzQ0FCOHz9OdbU6bhsWFkZaWhq5uerQzNChQ1EqlY2dq9zd3enXr19jVUkrKytCQkKIjY1tXNgjiiKJiYnk56tLEQQEBFBbW0tycjIAHh4euLi40FBuxMbGhuDgYI4cOYJCoZ5MvG/cfQTUBzDHaA6ptalYuFqQWZpJWkEaIiKudq44WTtRlV+Fj6kPfs5+BAUFER0dTcylGARBYPz48Zw5c4bi4mIAgoODKSoq4vLly912nhqq9erjeQoPDycpKald52lqcDBWJZcIdFWAK4wbN4kTf54l9UohJbUi1s79yC+pJDOvgMp6ERNzKwS5CfkFhdTWK7CysMDBwZ6Cq/mYysDMWMDPZwDlRXmYKGuwNRUIG+6PqbKKK1kZQC3e3u6N5+mvzJ5xP1VUVBAVFdVl5+n6+2ncuHEkJCRQWKieUwoMDKS8vJzU1FQAvLy8cHBwoKHtgL29PYGBgURHRyOKYpfeTx1FFyOHJcDk5mEkQRBSgPnXO4fmDBs2TDx79qxWbTVUUlJS8PHx0bUZeoGkRROSFk1IWjTRJSMHzYO8JYUPiKIY2c73LAKuXx3V8owoYGzcfSt99R0XF8MoYNcdSFo0IWnRhKRF52nROYiiuKUL3jOVmziDlkYNIDX7aU5cXJwUT9UgadGEpEUTkhadp1tmFQVB8NZkJd3gBDRpre0dfUhISEhIdCFacw6CIAQLgrBS8/91mtTWBtYBC5r9HiEIwkrNNvOAVhsKS81+mrCxsdG1CXqDpEUTkhZNSFp0HoMpnyE1+5GQkJBoPx2dkNaPZPU20JD+JkFjeqCEpEVzJC2akLToPAbjHAxlhNMdNORWS0haNEfSoglJi85jMM5BQkJCQqL7kOYcDBCVSoVMJvl1kLRojqRFE5IWTfT4OYeG5fISkJCQoGsT9AZJiyYkLZqQtOg8BuMcpBhiEw31WyQkLZojadGEpEXnMRjnICEhISHRfRiMc7Cw0G3fAH0iMDBQ1yboDZIWTUhaNCFp0XkMxjkolTdvhNIbKS8v17UJeoOkRROSFk1IWnQeg3EODbXXJWisES8hadEcSYsmJC06j8E4BwkJCQmJ7sNg1jkIglAOSH1C1fQBClrdqncgadGEpEUTkhZN+IqiaN3enXTRJrSjXOzIQo6eiCAIcZIWaiQtmpC0aELSoglBEDq0elgKK0lISEhI3IDkHCQkJCQkbsCQnENXtCw1VCQtmpC0aELSoglJiyY6pIXBTEhLSEhISHQfhjRykJCQkJDoJiTnICEhISFxA5JzkDB4BEHYrGsbJCT0AUEQDtzkb96CIKwUBGGS5l+7Nh1L3+ccBEHwBuYB8UAwsEUUxRLdWqUbBEEIBkIAO2Ak8IIoir26ToAgCJOAA6IoCrq2RZcIgjCv+e+iKH6rK1t0SbN7BNT3ybe94R7R3AfewObr7wVNs58Rmv/bAVtFUZzf2jENYRHc7mYfLA7YCrT6wXoampMaIoriFs3vk4ADgI9ODdMhGk2KgF75ZaEBQRBWAqmiKH6r0eQg0CudAzBJFMX1Db9oRpVP6NCebkEUxUi4cRStcZZFzbYr0Tw7WkWvw0o3+2BAmz5YD8QbeKHZ73GAd1uHiD2USaIoxuvaCF2iOf8vNYwURFEsafgy1Ut5opffE9fjzY1fnoo0z9YW0feRwy0/WG97KIiiGC8IwuRmfwoBSnpxiG0SEKlrO/SAECBVE1YqQR167RWhlFuwDkgTBKHhi9QLLW3cC3Do6I767hw6/MF6Itfd8E8AEbqyRZdo5qGKeqtjvA5v1A4hUhMyiANO0UvDjaIobtGMHBpCSZH07rBjEeq5l+a06bmq12ElOvHBejKCICwBdvXWSUfUD8MQQRCWaLSw0/zfW9eG6YBU1PMNJdAYevXupVogCMJKURTXa0Jrm1HPy/VmUrnJM7MtkRd9Hzl0+IP1VDThlNSGCajeyPVOURCEzQ0T9b2Qm4WPeuU3Zc290fhs0IwifHpjGLoBTTi68XfNl4Y2PTv0euRw/QltzwfriTRM0DfLTJjXyi49GkEQ7DSZOmjyt3vdt2VNqLGkYRJW829qL51zKEI9qryG3uAYBEEIbnYvrLsuIymiYZ0D6mUBbQpHG8I6h2DUGUq9ep2D5sGXct2fU0VR7JWxZYkmNNfGE6ivDx/Uue690Tk0fGFqiDbYoZ6L6fHOoSvQe+cgISEhIdH96HVYSUJCQkJCN0jOQUJCQkLiBiTnICEhISFxA5JzkJCQkJC4Ack5SEhISEjcgOQcJCQkJCRuQHIOEhISEhI3IDkHCQkJCYkbkJyDhISEhMQN/H9vduwosyz+hAAAAABJRU5ErkJggg==\n",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
230
      "text/plain": [
231
       "<Figure size 432x288 with 1 Axes>"
Jeremy BLEYER's avatar
Jeremy BLEYER committed
232 233
      ]
     },
234 235 236
     "metadata": {
      "needs_background": "light"
     },
Jeremy BLEYER's avatar
Jeremy BLEYER committed
237 238 239 240 241 242 243 244 245 246 247
     "output_type": "display_data"
    }
   ],
   "source": [
    "eigensolver = SLEPcEigenSolver(K, KG)\n",
    "eigensolver.parameters['problem_type'] = 'gen_hermitian'\n",
    "eigensolver.parameters['spectral_transform'] = 'shift-and-invert'\n",
    "eigensolver.parameters['spectral_shift'] = 1e-3\n",
    "eigensolver.parameters['tolerance'] = 1e-12\n",
    "\n",
    "N_eig = 3   # number of eigenvalues\n",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
248
    "print(\"Computing %i first eigenvalues...\" % N_eig)\n",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
249 250 251 252 253 254 255 256 257 258
    "eigensolver.solve(N_eig)\n",
    "\n",
    "# Exact solution computation\n",
    "from scipy.optimize import root\n",
    "from math import tan\n",
    "falpha = lambda x: tan(x)-x\n",
    "alpha = lambda n: root(falpha, 0.99*(2*n+1)*pi/2.)['x'][0]\n",
    "\n",
    "plt.figure()\n",
    "# Extraction\n",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
259
    "print(\"Critical buckling loads:\")\n",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
260 261 262 263 264
    "for i in range(N_eig):\n",
    "    # Extract eigenpair\n",
    "    r, c, rx, cx = eigensolver.get_eigenpair(i)\n",
    "    \n",
    "    critical_load_an = alpha(i+1)**2*float(EI/N0)/L**2\n",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
265 266
    "    print(\"Exact: {0:>10.5f}  FE: {1:>10.5f}  Rel. gap {2:1.2f}%%\".format(\n",
    "           critical_load_an, r, 100*(r/critical_load_an-1)))\n",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
    "    \n",
    "    # Initialize function and assign eigenvector (renormalize by stiffness matrix)\n",
    "    eigenmode = Function(V,name=\"Eigenvector \"+str(i))\n",
    "    eigenmode.vector()[:] = rx/np.max(np.abs(rx.get_local()))\n",
    "\n",
    "    plot(eigenmode.sub(0), label=\"Buckling mode \"+str(i+1))\n",
    "\n",
    "plt.ylim((-1.2, 1.2))\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above, we compared the computed FE critical loads with the known analytical value for the Euler-Bernoulli beam model and the considered boundary conditions given by:\n",
    "\n",
    "\\begin{equation}\n",
    "F_n = (\\alpha_n)^2 \\dfrac{EI}{L^2} \\quad \\text{with }\\alpha_n \\text{ solutions to } \\tan(\\alpha) = \\alpha\n",
    "\\end{equation}\n",
    "\n",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
289
    "In particular, it can be observed that the displacement-based FE solution overestimates the exact buckling load and that the error increases with the order of the buckling load."
Jeremy BLEYER's avatar
Jeremy BLEYER committed
290 291 292 293
   ]
  }
 ],
 "metadata": {
Jeremy BLEYER's avatar
Jeremy BLEYER committed
294
  "celltoolbar": "Raw Cell Format",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
295 296 297 298 299 300 301 302
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
303
    "version": 3
Jeremy BLEYER's avatar
Jeremy BLEYER committed
304 305 306 307 308
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
309 310
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
Jeremy BLEYER's avatar
Jeremy BLEYER committed
311 312 313 314 315
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}