{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Axial, bending and shear stiffnesses of generic beam cross-sections\n",
    "$\\newcommand{\\dS}{\\,\\text{dS}}\n",
    "\\newcommand{\\deg}{^\\circ\\text{C}}\n",
    "\\newcommand{\\bsig}{\\boldsymbol{\\sigma}}\n",
    "\\newcommand{\\be}{\\boldsymbol{e}}\n",
    "\\newcommand{\\bn}{\\boldsymbol{n}}$\n",
    "\n",
    "In this tour, we investigate the computation of stiffness moduli (axial, bending and shear) of generic beam cross-sections. In particular, we discuss the [computation of shear correction factors](#Computing-shear-correction-coefficients).\n",
    "\n",
    "The illustrative application considers a temperature distribution on the cross-section reducing the material mechanical properties and affecting the computation of the corresponding stiffness moduli. The animation illustrates the evolution of the transverse shear stress distribution as a function of the temperature increase.\n",
    "\n",
    "![](shear_stress_evolution.gif)\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "Complete sources including mesh files can be obtained from :download:`cross_section_analysis.zip`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Mechanical properties of beam cross-sections\n",
    "\n",
    "### Symmetric homogeneous cross-section\n",
    "\n",
    "We consider a Timoshenko beam oriented along the $x$ axis and consisting of a symmetric cross-section in the $(y,z)$ plane and made of a isotropic homogeneous elastic material (Young modulus $E$ and shear modulus $\\mu$). We denote by $(y_0,z_0)$ the location of the beam axis. The beam linear elastic beahviour for planar bending in the $(x,z)$ plane is described by the following constitutive laws:\n",
    "\\begin{align}\n",
    "N &= EA \\epsilon \\\\\n",
    "M_y &= EI_y \\chi_y\\\\\n",
    "V_z &= \\kappa_z \\mu A \\gamma_z\n",
    "\\end{align}\n",
    "where $\\epsilon$ denotes the axial strain, $\\chi_y$ the bending curvature strain and $\\gamma_z$ transverse shear strain. $A$ is the cross-section area, $I_y$ the bending inertia and $\\kappa_z$ is the shear correction factor described in the [next section](#Computing-shear-correction-coefficients).\n",
    "\n",
    "### Generic heterogeneous cross-section\n",
    "When the cross-section is made of a heterogeneous material, the above relation must be replaced with equivalent stiffness moduli. Moreover, when the cross-section is no longer symmetric with respect to the $z=z_0$ axis, an axial/bending coupling appears. The generic constitutive relations are therefore:\n",
    "\\begin{align}\n",
    "N &= H_{N} \\epsilon + H_{NM}\\chi_y \\\\\n",
    "M_y &= H_{NM}\\epsilon + H_{M} \\chi_y\\\\\n",
    "V_z &= H_V \\gamma_z\n",
    "\\end{align}\n",
    "where $H_N$ denotes the axial stiffness modulus, $H_M$ the bending stiffness modulus, $H_{NM}$ the axial/bending coupling stiffness modulus and $H_V$ the shear stiffness modulus. These moduli are defined as follows:\n",
    "\\begin{align}\n",
    "H_N &= \\int_S E(y,z) \\dS \\\\\n",
    "H_{NM} &= \\int_S E(y,z)(z-z_0) \\dS \\\\\n",
    "H_M &= \\int_S E(y,z)(z-z_0)^2 \\dS \\\\\n",
    "H_V &= \\kappa_z \\overline{H}_V = \\kappa_z \\int_S \\mu(y,z) \\dS\n",
    "\\end{align}\n",
    "\n",
    "In the homogeneous case, we indeed retrieve $H_N = E A$, $H_{NM}=EA(z_G-z_0)$, $H_M=EI_y$ and $H_V=\\kappa_z\\mu A$ with $z_G$ being the position of the cross-section centroid. As a result, we have no coupling ($H_{NM}=0$) when the beam axis passess through the cross-section centroid. \n",
    "\n",
    "> Note that we always have $\\Delta_H = H_NH_M - H_{NM}^2 \\geq 0$ due to the Cauchy-Schwartz inequality. This ensures that the above behaviour is well-posed.\n",
    "\n",
    "For complex geometrical shapes, the computation of the above moduli must be done numerically, for instance using a 2D mesh of the cross-section.\n",
    "\n",
    "### FEniCS implementation\n",
    "\n",
    "We illustrate the computation of these various moduli for a typical box girder bridge cross-section subject to fire conditions. We consider here only the concrete cross-section and neglect the presence of steel rebars and pres-stress cables in the computation. The cross-section geometry is defined and meshed using `gmsh`. We readily compute the cross-section area and the centroid position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAADbCAYAAACcGb3IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOy9WXBcZ3ag+d3c9z2RCxI7QBAkuIG7RElUlaSoKrtccrTb81oRnnHH+N3t1/LTRNvz5qeefvGLHTGebke4vFWVSyqSoiiJIkGABAhiR2YiF+S+I/c7D1ReayWJEsBMXecXwVCAAhLn8Pz3nv8/2y+IokifPn369OnzLBTdFqBPnz59+vQ+fWfRp0+fPn2eS99Z9OnTp0+f59J3Fn369OnT57n0nUWfPn369HkufWfRp0+fPn2ei6rbAhwWDodDHB8f77YYR0a9Xkej0XRbjCNDzvrt7e3h8Xi6LcaRIWfbgfz1u3//fkoURffzvk82zsLtdnPv3r1ui3Fk3Lhxg+vXr3dbjCNDzvr99Kc/5a//+q+7LcaRIWfbgfz1EwQh+CLfJ5swlNFo7LYIR8rly5e7LcKRImf9BgcHuy3CkSJn24H89XtRZOMsGo1Gt0U4UnZ3d7stwpEiZ/0KhUK3RThS5Gw7kL9+L4psnEW9Xu+2CEdKJBLptghHipz1KxaL3RbhSJGz7UD++r0osnEWffr06dPn6JCNs9DpdN0W4UiZnp7utghHipz1czqd3RbhSJGz7UD++r0osnEWgiB0W4QjRalUdluEI0XO+ikUsnnMvhY52w7kr9+LIptVvL+/320RjpTHjx93W4QjRc76JZPJbotwpMjZdiB//V4U2TiLPn369OlzdMjGWajV6m6LcKR4vd5ui3CkyFk/k8nUbRGOFDnbDuSv34siG2eh1Wq7LcKRMjY21m0RjhQ562ez2botwpEiZ9uB/PV7UWTjLEqlUrdFOFI++uijbotwpMhZP7k3dcnZdiB//V4U2TiLPn369OlzdAiiKHZbhkNhZmZGXFlZeSm/KxqNEgwG2dvbQxRFRkdHSSQS7O3t4fV6cTqdrK2t0W63OXHiBPl8nmKxSKFQ4Pjx4ySTSRQKBc1mE6PRSC6Xw2q1IggCzWaTYrHI8PAwqVSKdrtNq9VCrVbTaDTIZDLMzMyQTqfZ39/HZrPhcDjY2NhApVIxMDBAq9UilUqh0Wjw+Xzs7e0Rj8cJBAIIgsDu7i4ul4vBwUGi0Sj1eh2n04koimxtbeHxePB4PMTjcZLJpKTT5uYmrVaLyclJKpUK2WyWRqPB8PAw8Xgcg8GA2Wwml8tRq9WwWq3odDrW1tZwOBw4nU7q9Tq5XI5KpcLU1BQrKysMDAxQrVYRBAFRFLFYLGi1Wvb39ymXy3i9XkqlErVaDaVSycDAAMFgEFEU8Xq9KBQKdnd3qVar0r9NvV4nm80yNjZGPp8nmUxy8uRJgsEgKpUKURSx2+2k02kEQWBiYoJIJEKtVsPhcJDJZNDpdIiiiCAI6HQ68vk85XKZiYkJgsGns9dGRkbY2dnBZDKRzWaZmJggmUySzWaZnJwkk8nwN3/zN/zRH/0RFouFdDpNMplkfHycSqVCrVbDYDBI+lYqFZxOJ7Vajd3dXYaGhlCr1ZRKJRKJBJOTk9TrdSqVCplMhvPnz/PkyROUSiXNZpPBwUFisRitVouJiQlWVlYwm82YzWZSqRRWq5Vms8nAwAAbGxuSzex2O5ubmwiCwNjYGLlcjkKhgFqtxufzEYvFUCgUUv4llUpRKpXw+XySLOPj4yQSCXQ6HUqlEofDwd7eHul0mvHxcZRKJclkEqVSicfjoVQqkc/nUalU6PV6yuUyRqMRk8nE+vo6JpMJp9NJu90mGAzidDoZHBwkmUwSiUQYHBzE4XCwu7uLIAiMjo4SjUZJJBI4nU5JJ7vdTiAQIJFI0Gg0UKvVuN1uYrEY9XqdoaEhSYdms8n09DQPHjzAbrfj8XhIJBKIoohCoWBwcJBQKIRWq0WpVGIwGEgkElgsFqrVKlqtlnK5jMPhIB6PYzQaUSqVOJ1OHj58iNVqZWBggGazSTAYZGBgALfbTSQSQRAE1Go1169ff6l9Y4Ig3BdF8cLzvk82J4uXOUiwXq/TbrfZ29tDpVJx69Yt9vf3EUWRRqPBnTt3MJlMWCwWPv74YxKJBCqVCqVSya1bt3C73WxublIul7l79y7T09PMz8+zsbHBgwcPmJ2d5b333mNvb4/l5WVmZmaYn5+nVqvRbreJRCKSw1laWqJarZJKpRBFkdu3b+NwOAgGgzQaDX71q1/RbDbRarVkMhlyuZz0svzVr35Fo9HA5/Px8OFDMpkMFouFer3Ov/3bv1EqlWi1WlSrVT744AMCgQBGo5FgMMijR48kx7a6ukqr1WJzc5PNzU2KxSKDg4MsLCyQSqUQBIFGo8HS0hIrKyu4XC5UKhUbGxvSw14ul6nX6+j1ehYXF8lkMiwsLDA7O8uNGzfIZrOsrq4yNDTE+++/T71ex2q1sr6+zsLCAjabTbKFUqkkk8kA8OGHH9JqtXA4HKyvr1OtVkkkEoyMjFAsFqUH++///u9pNBp4vV7m5+dpNpuEw2HppbCwsIDb7UatVnPnzh2cTiejo6NSiCIajaLX60kmkzQaDfb399ne3iaXy9FoNGg0GiwuLlKr1XA6nQSDQTKZDAaDgWAwyMbGBqVSCZvNxubmJvl8HrvdzsLCAvF4nGKxiMVi4f79+yQSCcrlMjqdjpWVFfb398nn87hcLrRaLblcjmq1SjgcxmQyIQgCd+/epV6vs7q6it/v5+bNm7RaLQAymQwffvghBoMBk8nEBx98gMFgYG9vj3K5zPvvv8/IyAjr6+tks1k+/PBDVCoV7XZb+vd0uVwEg0Ha7bbkzB48eEC1WkWtVpPNZrlz5w4WiwWNRsPW1hYPHz6kVCqxu7tLJpNBEAS2t7fZ3t7GZDLhcDhYWVkhk8nQbrepVCq8//77VCoVrFYrlUqFmzdv4nK5MBqNvPfee5JOtVqNaDSKRqNBEAR++ctfUq1WWV1dZXJyUlr7oVAInU7H3bt3aTabxGIxarUa9Xqder3O+++/z/T0NGtra6TTaf71X/+V2dlZHj58yM7ODrdu3eL06dN8+umnlEolHj58yNTUFB999BHNZpO9vT00Gg23bt1CpVJhMpnY3NykUChI/1Y3btyg2WySSCRot9u02+2X9i47CP2TxW9BMBiUXt4AKpWKTCaD3+/nyZMn/OQnP+HGjRsA0o5qa2tL2rVlMhlKpRIajQaXy0UqlaLRaFAsFjl79iyRSIR2u008HufChQuEQiEsFgvtdpvJyUkePnyITqejWCzi9/uJRCJotVqKxSJOp5NcLse1a9eYn59nbm6Of/qnf+L111+XHNnJkydZXV3ld3/3d/nggw9wOp0YjUbW19c5deoUi4uLvPbaa+zs7ODz+Zifn2dgYIBAIMDa2hrxeJzh4WEqlQperxeTyUQikZBOAna7HYfDQTKZJJfLMTU1JclYrVY5d+4cT548IZPJMDExwdbWFkNDQwQCAR4/fsyZM2eYn5/n4sWLLC4uMjU1xcbGBna7nUgkwsTEBJubm7hcLiqVCkNDQ8zPz9NqtTh16hSJRIJsNovH40Gv1/Ppp5/i9Xqp1+vYbDZp11ur1YhEIrz22mu0Wi0++eQTXnnlFe7du4fL5SKdTjMyMsLjx4+xWCycPn2ax48fIwgCOzs7/OhHP+LDDz9kcHCQVCqFz+dDrVYTj8cZHBxkb2+PRqPB3/3d3/Huu+/idDpptVo0Gg1EUaRSqWAwGDAYDJRKJelk0dkJWywWWq0WGxsbjI6OSqeunZ0dVCoVZrMZk8lErVZjYWEBvV7P0NAQyWSS/f19isUiKpWKRCLB9PQ0VquVUCiESqUiEAiwubkJwPDwMO12m/39fWq1GhMTE2xsbOB0Otnf3ycQCBAOh5menmZhYYELFy5w//59xsbGMJvNPHz4EIfDwdbWFm63G5fLJe3GAVwuF5ubm9hsNlqtFj6fj5WVFaxWK0ajEZVKhVarZXV1lenpaR4+fIjZbObkyZPk83kikQizs7OkUim8Xi93797l1KlTPHr0CLvdLp2WJicn+Zd/+RdGR0cJhUK88sorqFQqbt68yY9+9CPpebh16xZOp5N0Oo3NZiOfzxMIBEgmk9hsNrLZLJcvX+bRo0fYbDbS6TRut5vHjx8zMDBApVKRnkOr1Uq9XsflcklOsnNqrdVqFItFTCaT9G5oNBq88sorvPfee9I69nq95PN5bDYbGo2Gt99+G4PB8FLeZfDiJwvZOIvp6WlxdXX1pfyuzjH4xo0b+Hw+6vU6hUKBmZkZ6QXkcDgYHBxkdXUVQRCw2WwYDAYikQhqtRqn00mz2UQURdrttlT6m81msVqtmEwmWq0WhUIBr9dLIpEAwGAwIIoiiUSCsbEx9Ho9zWaTxcVF1Go1oihy5coVPvnkE2ZnZ9nY2GBiYgJBEFhbW0OhUHDy5Em0Wi23bt2SHg6n08nAwACLi4u4XC6azSaBQICFhQWuXbvG0tKSFAYJBAJsb29jt9uJRqMYjUaq1SoOhwOTyUQ6nSaVSjE7O8vS0hJjY2M4HA7+4R/+gTfeeAOn08mvfvUrxsfHKZVK0ktxcnKS9fV1SqUSFy5cYH5+nrGxMRKJhBS2cDqdKJVK9vf3SafTDA0NEQwGUSgUmM1mdnZ20Gq1jIyMEI1GKRaL2O12LBYLsViMUCjEO++8w8LCAjqdjpMnT7K8vEy1WsXpdGKz2VCr1XzyySdcu3aNUChENpvF7/czMDDAp59+ikajwe/3o1arEQSBx48fc+XKFba3tymVSmi1WsbHx0kmkywsLLC6usqPf/xjKpUKgUAArVYrnTreeust1tbWCIVCBAIBAoEAGxsbNJtNnE4njUaDbDaLwWBgdHSUe/fuYbVacbvdKJVKlpaWcLvdHD9+nGg0SjqdZnh4mK2tLfR6PWazGavVSiwWo91uMz09zc7ODuVyGZfLhdlsplwuS+EWm82GXq+n0WiQy+VQKpVSaG55eZnLly+zuLgofXbn1JTP56WT3uuvv44oivziF7/g+vXrmEwmHj58yPT0NOvr61itVvR6PaIoEo1GcTgciKKIy+UiEokwPT1NOBxma2sLrVbLW2+9xS9/+Us8Hg9ut1sK+XZCurVajeHhYZLJJOfOnaNUKhEMBqlWq8zOziIIAgsLC5w4cYLNzU1OnjzJ/Pw8oihSrVa5cOGCFO5dWVlhYmJCciQA5XKZZrOJyWRCo9FQKpXI5XLo9XrsdjvFYlE6EXTCT53Tj91uJx6PUy6XGRkZQRRFwuEwFosFv9/Po0eP0Gq1mM1mtre3uXbtGufOnXup1Z19Z3GE7Ozs8ODBA1wuF4uLi1JMP5VKsb+/L+2ki8UiAwMD6PV60uk08DTG3Wq1SKfT0lHabDZTLBYRBEHKXVSrVZrNJoAUfqpUKqTTaZxOJ0NDQ2xvb1OtVtHpdIyOjkpxz3A4zKuvvsri4iKDg4NYLBa2trZoNBrSCUWhUHD8+HE+/vhjjEYjo6Oj3L59m6mpKdrtNuVyGUEQcLlcmEwmrFYr8Xic/f19MpkMp06dYmtrC6VSyaNHj/jRj35EuVxGpVLx4MEDLl++TCQSIRAIEAqF2N3dxWq1SqcppVJJMBikUCjwxhtvEA6H0el0FAoF9Ho9y8vLXLx4EYVCwfb2NiqVCqfTidPpZGVlBaPRiM/nk+LTSqWSXC6HwWCgWq2STCYxGo14vV6i0SiBQIC9vT0qlQp6vR6j0Yher2dlZQW9Xi+dcDpORxRFKZ9ht9tRKBTs7++j0Wiw2WxEIhGUSiU2m00KxZjNZoaGhkilUtIuOZvN8otf/IJ3332X0dFRwuEw4XCYiYkJ8vk8ly9fplgs8vDhQ9xuN0+ePOHixYvk83mUSiWlUgmr1UqhUMBgMNBqtRBFkWKxiM1mo1QqYTKZ2Nvbw+fzodfrKRQKVKtVKea+srKCz+eTHKzT6ZTyMZ0wWygUkk4l+Xye4eFhBEGQXnROp5ORkRF+9atfMTU1RS6X48KFC4TDYQApHzA0NCSFcZRKJdvb20xMTDA0NEQ0GmVgYICPPvqIixcv0mq1MBqN/OIXv2BmZkY6Od+7dw+Px4PZbMZisZBKpaQwascpKZVKnjx5wvXr19nZ2aFQKDA7O0uhUGB1dZW5uTlWV1dRKBS43W5pIzQ0NMTOzg7j4+Nks1kCgQDLy8t4PB7K5TLj4+PSpkqpVDI0NEQmk8FoNCIIAiqVSlpjjUZD2iR1QripVAq/3y+Fi/f396X1HovFUCqV0olCoVDg8XjI5/Osr69z+vRpyuUyP/jBD3ryZCGbnMXLbHxqNBoYjUYeP37M4OAgrVZLSgiOj49Lidjx8XEGBwcxm82o1WqsVqu0yPb391lfX0en05HJZKhUKlIcMxKJEI1GpdBNpVIhmUyi0+k4fvw4lUqF+/fvMzo6itPppFKpcPfuXYaHhykUCpjNZh49esSlS5coFos8efJEckQAQ0NDiKLIRx99xNzcHI1Gg+3tbc6fP08sFmNwcJCpqSlsNhuhUEhKdJtMJoxGIx6Phzt37jA+Po5Go+EnP/kJNpuNx48fUygU+P73v08ymSQQCHD37l1yuRzj4+O43W5GR0dJp9NsbW0BcO7cOfL5vBQ7HhgYYHNzk6mpKSkH4vF4GB0dJRaLsbq6ypkzZ6hWq6ysrCAIAlqtlkgkIp1SfD4fJ06cwOFwEIvFsNvtrKysSLvHzsOfy+VwOp34fD6Wl5cRRZHJyUmSySSJRILR0VEpMZpKpRgbG6NQKLC+vo7T6cTj8bC6ukq5XObcuXMIgsDy8jLNZpN3331X2p22220GBwdZXl5GpVLxgx/8QJIhk8mwvLyM3W4nkUjw7rvvSiE9h8NBIBAgFouh1+ulcGGpVJJerh0n4HY/vRVTq9WSz+fZ39/HYrFIuZZ8Po9er+fUqVP4fD4p2Z5IJCiVSgwMDPDkyRNSqRRzc3OkUik2NjbQarWcPHmSUqnEvXv3mJycJBgM4vF4pMRzu91mbGwMm81GOBymWCwyNjZGOp3m6tWrlMtlPvzwQzwej7Q+isUijx49wmq18vu///u0Wi1mZma4c+eOJJ8gCFLeofMCP3PmDG63m729PS5fvszW1hbFYpFLly7x5MkTYrEY169fJxKJYDKZpBDpysoKV69eJZVKoVarpTW2uLiIQqGg0WhIoaXOpsrv90unxc4pPx6PS+G7dDpNq9WSNkILCwuoVCqy2SxKpRKFQoHD4cDlcuF2uxkcHESn00mn9nK5TD6fZ2dnh9dee42trS0p99mLKH/2s591W4ZD4a/+6q9+9id/8icv5XdVq1Xm5+dRqVQoFAopnnrs2DFisRiZTAa1Wi0lsBYXF7FYLBgMBulFabPZ8Pl87OzskM1mcTgcaDQaVldXaTab+Hw+Go2GVOXk8XioVCpSgtLn85HL5dje3mZmZkY6nRQKBebm5tDr9UQiEex2OydOnJAqUcbGxvD5fBQKBWKxGCaTSQonCIKAQqGgXC5Lu+VLly5JL5VOVczOzg4/+MEP+PWvf40gCBiNRhYXFzEYDAQCAekB393dlcIn2WwWhULBnTt38Pv91Ot1jEYj+/v7bGxscP78eR48eEAwGGR0dJRkMsmFCxeYmJiQKsvGxsaoVCosLS3h9Xrxer3s7OxQKpUYHBykWq3y6NEjBEGQQikul0t6AIeHh6UksiiKOBwOVCoVOp2OM2fOoNfrqdfrBAIBpqenpeTs7OysdILweDzMzs5KG4Lz589jMpmoVqu4XC5Onz6NUqkkGo1y/vx59Ho99+7d4+rVq5w5c0barZ85cwaz2Uy9XpecodVqlXIDHo9H+h1nzpyRXiCdHX7nRDE2NoZSqUSlUkmhq2q1ikqlkmL5arVaCuO1Wi0qlQr37t2TKuB0Oh3Ly8t4vV78fr/0wjx37hzVapXd3V0uXrwoJeZdLhfb29tEIhFef/11wuEwtVqNRqNBu93GYrHw4YcfcuzYMcLhMOPj42i1WqLRKFqtllqths1mo9FosLa2RqlUwuFwcPPmTX7nd36HhYUFGo0Gg4ODBINB8vk8V69elZLffr9fsqlCocBgMFAul0mn08zMzEhTfjOZDAMDA8zMzKBSqaTnYXx8XMoTCoLA7Owsy8vLqNVqWq0WJ06c4MGDB+h0OhQKBRaLhZWVFanKy2w2s7S0hNFoxGKxUCqV2N7eZmRkBJ1Ox+7uLolEQqp66hRWmEwmKYTVbDal0+WlS5d48OCBtIE5f/48KtXLu/H6z//8z2M/+9nP/p/nfV8/DPVb0Gq1aDabVCoVaYBhvV4nnU6jVqux2+1kMhnpQiaPx0MkEpEeXLfbzc7ODq1WS8pP7Ozs0G638Xq9iKLI7u4u9Xpdiuvn83lSqRTnzp0Dnl6ok8lkOHbsGDqdjmQySTwex+fzSYnMSqVCvV7n9OnTUsVSZxf8m9/8RvpsnU4nleIGAgHpdDIyMsJHH33EtWvXWF1dJRaLodVqeeONN1heXpZKZvf390kmk4yMjGC329nY2MDn80mnlcXFRTweD5lMBrfbLe3uZmZm+PWvf834+DitVguVSkWxWJRKLTt5G41GQy6Xw+Vy4fF42N/fZ3l5mXq9zrlz58jlcuzv72MwGKTEskajYWBgQCrHdDqdGAyGrkwn7tYd3KIoSi/RQqGA0+kkGo1SrVYlR5XP5wGw2+3Mz8/jdrsZHx9HrVazu7tLo9HAZrMRDAaZmJggGo0SDoelnEOnumljY4NXX32VRCJBPp+XSnUnJyfZ2dnB5XLRaDQwGAxsbW0xMzODQqHg0aNHGAwGhoaGcDgc7OzscOzYMe7cuUOtVsPlcnH27Fnef/99Tp06RbvdZnt7G5vNJoUnjx07BkCz2WRzc5M33niDjY0NFAoFoihy4sQJaf0bjUbphBgKhRgdHZX+LTobsU64aXFxEbfbjUKhwOv1sr29LZUaOxwOKQTndrvR6/VSebHP56PZbFIoFKRy22w2S7PZlMqGO++Hzmc1Gg10Oh16vR6tVvtS12k/ZyEjRFHk5s2bXLp0CVEUpcqZarUqhRNarRaCIEgJtk5CTq1WS1NP7XY71WqVTCaDXq9HrVZLcf6lpSVGR0exWCxSpUg8HmdkZER6UfzjP/4jp0+flnaUH3zwAQAzMzMMDg7y6aefSiGB3/zmNxw/flyK2apUKum4XiqVGBsb48MPP8TtdtNoNKQQmSAIUuK6U9ff6RXRaDQ4HA50Op3Ui/FdoFvO4rehWq2STqdRKBQkEgnphdZoNKSQVS6XQ61Wo1arcTgc5PN5arUaiUSCmZkZ6VTZqQJKJBLkcjlOnDjBzZs3mZiYAJ5eN/vKK6/w8ccfo1AoKBQKXLp0iU8++UTKTV27dk3qc4nH4zidTlwul5R/CoVCTExMSPmcTk+KXq/H5XJJL2mDwYDVapXyUBaLBZVKRalUQqlUSkUWSqUSnU6H2WwGnob1PvjgA65fv95Fqxwt32lnIQjCHPDWZ19eBP4PURRzz/qZs2fPigsLC0cuW7foHGsPk4cPH5LNZikWixSLRenlu7u7y8jICM1mU4qPRyIRhoeH+ed//mf+8A//kA8++AClUkmr1eLs2bN8+umnUkI3n8+jUCikcFOxWMRqtZJIJIjH47z99tvcu3cPg8GASqVCrVYTCARIpVJ4PB5OnDghq+F7f/qnf8pf/uVfdluMQyMWixGPx4nH43g8HtbW1qQXdSaT4eTJkzx48ACDwYBGo8FoNGK320mlUlKCXqvVYrfbefToEdeuXZP6dgRB4OLFi/z85z/n8uXL6PV6qVS5EzIMBoMEAgGpAMBms9FsNrl48SJ+v//Q9T2KZ6+X+M4muAVBsAEXRFH8C1EU/wL4f4H3nvdzvej0DpNOs9Fhcvr0aV555RWpXj+dTmM2m5mYmGB9fV16yNfX1/H5fASDQebm5tjc3JTKHc+fP080GpVyA+12m2PHjlEsFsnlcrTbbaLRKI1GA5VKxfj4uJRn6exQdTodqVSKt99+m0uXLsnKUQA922T12+Lz+Th37pyUrzAYDOj1etbX1zl79izBYJCRkRH0ej0mk4lwOEwul0OhUBCLxRgfH5eKPHw+H6FQiNnZWfR6vdSwd+3aNSnEGw6H2d/fZ3Jyko2NDUZGRjCbzSSTSaxWKxaLhR//+MdH4ijgaJ697yI95yyAC8Cffe7rXwNznzmRb6RarR6pUN3mqEJsarWaN954g6tXr0plgltbW5w/f55IJCKNiuiUG9ZqNdxuN4VCgYGBAW7fvs3Q0BBjY2Pcv3+f6elpIpEIZ8+eleQ+c+YM4XCY4eFh6vU64XAYhULBlStX8Pv9vPXWW5jNZtlejdspm5YbRqORn/zkJ9JYjjfffJN4PC6NpfH5fJTLZaanp6XS7cuXL7O6usrs7Cz37t1jamoKo9HI6uoqJpOJRqOB3+8nk8kwOjrKzs4OoigyNDTE4uIir7/+OslkklQqxdDQEK+++iqvvfbakd5GKNfw9kF5eSn3F0QUxV8LgvCfP/dX45/9/TPDUH2+HS6Xi7feeosPPviA48ePS7vDThlgLpeTSik1Gg0zMzNSb0E+nyeTyfC9731PGr3RyZecPn0al8vF1atX2djY4OTJk2xubjIwMIDf75diw32+uygUCt544w3i8Thut1vqul9ZWeHYsWPSCJJOFZPD4eDBgwf83u/9Ho8ePcLpdEqJ3eHhYSkRf+fOHTwej1T4oNVqpTVVq9V49dVXZX9lbS/Rc84CQBTF+c99+b8Bf/G8n9FoNEcnUA8wODh45L9DpVLx5ptvsrKygs1mY3BwkPfffx+XyyU1W3Xm6dhsNmlsydLSErOzs5RKJQRBIBgMks1mpfEclUqFYrGIx+MhmUzy7rvvfuWyqpehX7eQu0Ps2K5Tzjw0NMRvfvMbbDYby8vLCIKA2WwmGAxiNpulHojHjx9z6tQpPv30U6rVKgaDgbGxMWnqgUajQa1WYzQapYkAoijidDqlJPnL1GQS69AAACAASURBVO8/Oj3pLDp8FnqaE0Xx7W/4/38M/DGA3+//wjwms9nM4uIi8LQ2/eTJk9y6dQt4+lLszE4qFAoAXLhwgb29PakcbmpqCq1Wy9LSEgADAwMcO3aM27dvA0+rJK5evcq9e/ekuzQuX77M7u4ukUgEgOnpaZRKpXSHr9frZWxsTBo+p9fruXz5Mp988okUn+3Uk8fjcQBOnDghNf50Jm0GAgE++eQT4Gkz4oULF/joo4+kWVXXrl1jbW1NGhHS6QtYX18HnjbleTwe7t27B4DFYmFubo7bt29LXeMXL15kYWEBr9crlSsODw8DT/NDjx8/xuPx4HQ6pUGBnVElWq0WURQ5fvw4y8vLwNOTy6VLlygUCnz44YdfsVNn9MJ33U6dkMXn7VQoFLh3796R2On1119neXlZCnWdOXOGYrEoNT2Ojo7icDiYn3+6/7Lb7Zw5c4abN29KE3XfeOMNFhcXyWazAMzNzZHJZNjZ2fmKneCrz5NCoWBqauoLdnrttdekcvF6vU6r1ZLGz7jdbilPFYvF8Hq9qFQqGo0G0WgUi8VCo9HgzJkzLC4uSp37nXLo1dVVwuHwodvpm56n8+fP8/jx4++8nb7peXpRerIaqoMgCP8d+LMXCUHJuXQW4MaNG10p3+tMC52amuK9997j1KlTrK2toVKpqNfrTE1Nkc1mabVa0njqVquF2WymUCjw9ttvMzAw8NzxBd3S72XwXSqd/W14nu2KxSKPHz+WRuW3Wi0sFovUdNipbAqHwzQaDWq1mlSBd+HCBTY3N7l27RoWi+XlKfU55Lw24TtcDdVBEIT/ymeO4nnJ7T5Hh91u54c//CHRaJS5uTlpltLw8DBWq1Uap12pVKQEpV6vR6VSMT09jclkeqlzbvr0HmazWWpc6wzY7ExdVigURKNRdnZ2UKvVjI+Po9PpiMVijI2NkUql+OEPf9g1R9Hn3+lJZyEIwh8A//NzJ4q3nvX9AEql8miF6jLdLCdVKBRcv34dr9fLwMAALpeLcrlMu93G4/FgMBioVCo0Gg2mpqakJGZnxPWLILdy2c8j93zai9iuMyfJ4/GQSqWYmJiQBg6KoihdetQZdWMymRgZGeHVV1/tStf955Hz2jwIPecsBEEYB/4/YFMQBFEQBBH4b8/7ObnvXg8SWzwqRkZGuHjxIplMRupricfj0o5xYGCAXC6H3W4nGAxSLpeloWvPoxf0OyqOqv6/V3gR2z158oR2u00oFJJClJ1R+CdPnqRQKJDNZgkGg1itVi5evNgziWU5r82D0HPOQhTFLVEUhS/9eW7pQ7lcfhnidY1euTQ+l8tJF/c0Gg1GRkYol8u0Wi1isZg09jmTyRAKhV64tLFX9DsKOsl4ufIitlOpVASDQXZ2djh+/DjZbFYa3/H48WOMRiNut1sq5niZg/Seh5zX5kHoOWfx2yK3Ltkv06nM6DZ+vx+LxYJSqcRutxOLxaSZVJ25VSqVinfeeYd2u/3C1932in5Hgdw7gF/EdkajkUajwe/8zu8ATzd3oihKd4t0ZkMNDg5Kd5f0CnJemwdBNs6iz8vBarViMBiIRqOoVCrpus90Oo3X68Vut7O0tMT29jaBQEAqT+3zH5toNCpNY37w4AF2ux2/30+pVKJer+P3+6V70j0eT9fzFH2+Su+c9b4lck9CXbt2rdsiAEg9FaOjo6ytrUn3SIyPjyMIAl6vl3A4TKFQIJFIYLfbpTrxZ9Er+h0Fnf4UufI823VmhsXjcak6rpOP6FxDnEgk2NnZeanNdi+KnNfmQZDNyULuR8W1tbVui/AFFAoFb775Jul0munpaTKZDMVikfv372O1WnE6nSgUCpLJpNRE9Cx6Tb/DRK6zoTo8z3adprHOdb0zMzPcu3ePWCwmzY7KZDJ8//vfR6/X99ycNzmvzYMgG2fRaDS6LcKR0uke7QXUarVU6dQZOGgymYhGo0xMTGC329nZ2cHhcGC1Wl/IkfeSfoeN3Isvnmc7i8VCKBRieHiYdrstdfenUikMBgMrKytcunRJOpH2WqmxnNfmQZCNs+jz8ggEAl+4frVWq0m3p3U6uQuFAm63W7oHu89/XEKhEKdOnZLyE6lUina7jdVqZWNjA6VSSaVSkarqhoaGui1yn69BNs5Cr9d3W4QjZXZ2ttsiSOTzeWZnZ0mn0wSDQek+5LGxMen+8R/+8IesrKzQaDSoVCrP/cxe0u+wkfPFOfB825XLZWl22KVLl0gkEqjVanw+H1qtlnK5TDgcJhwOc+nSpRcKW75M5Lw2D4JsnEUvz7g6DHopJ9MZWT43N0e73Uan09FoNNjd3WVycpJAIMDf//3fYzAY2N/ff6EwTC/pd9h0hsnJlefZrlarUa1WUavV/PznP2d6ehqn08ne3h71eh2dTiddpBWPx3uuwVbOa/MgyMZZ9FpS7LDpTLjsBSwWC9FolGQyiV6vl6pbyuUyDx48oNls8uabb0p3G7tcrud+Zi/pd9hkMplui3CkPM92DoeDVquF2+3m7bffJpFIEAwGKRQK0miPXC4nXYzVax3vcl6bB0E2zqLPy6MzI6rVapFOp7Hb7USjUYxGIxMTE2xublKpVCgUCmQyGYLBYLdF7tMlWq0W0WiUZrNJKBQimUxKjZvDw8Ps7u5is9lQKpVks1m8Xu8LN3L2ebnIxln0WgXFYdNLST+9Xk+5XCYWizExMcHu7i5arZbx8XHp+ky73Y5KpSIcDmM0Gp87ULCX9Dts5D4x9Vm2C4VCeL1eEokECoUCl8uFSqUiEokgiiJ+v59QKMTIyIh0uui1hjw5r82DIBtn8eWb1+SGx+PptggSnca8q1evUiqVSKfTmEwm8vk8w8PDZLNZNjc3cbvdTE5Okslknlt+2Ev6HTZybxh9lu3q9TqxWIzBwUGpYTMWi2Gz2RBFEYPBQC6Xo1KpcPny5Z7M78h5bR4E2TgLudeyd27h6hU8Hg+hUAhRFJmbm6NUKrGxsUGr1eLEiRNUq1U2Njaw2+2Sc3kWvabfYRKNRrstwpHyPNslk0kmJyfZ3t6mVCrh8/kYHBwkHA4TCoW4ePEizWaTYDDYc/kKkPfaPAiycRZ9Xi6CILCzs4PT6WRnZwev14vNZmNvb4/t7W0GBwcZHx8nHo8zNTVFMpnstsh9ukAsFuPVV19lfn4el8vF1NQUrVaLpaUlTCYTDoeDUCiE3+8nEonIfujidxnZOAu5X37Ua3HvSqXCiRMnyGQyUhmkRqNhZGQEURTZ3d2lXq+j0WikPMaz6DX9DhOtVtttEY6UZ9lOFEVWVlaw2+0olUq2t7fJZDK4XC6USiUqlQqTyUQoFOLs2bM9WaYq57V5EGTjLHqtNvuwmZub67YIX8BgMOB2uxkeHiYSiWAwGMhkMtItaHa7na2tLYaGhtje3n5u3LfX9DtMfD5ft0U4Up5lO5fLxebmJqOjo1IhxNjYGGq1mmKxyP7+vnSFqsVi6clKKDmvzYMgG2dRKpW6LcKRcvv27W6L8AX8fj8ajYb3338fm81GLpdjamqKyclJQqEQJpOJCxcuEI/HOXv2LE+ePHnm5/WafodJKBTqtghHyjfZThRFFhYWeOWVV9ja2uLChQtSZZTP52NycpJ0Oo3D4eD27duoVKqe7HaX89o8CLJxFnLv4O61KpFisUg2m2VsbIxKpYLT6SSTyXD37l3UajWlUol8Ps/e3h6CINBqtZ4Zj+41/Q4TuV/M9U22SyQSOBwOyuUyqVSKYrFIo9FAEAQWFhZIpVJ4vV7K5TJOp5NYLNaTJfByXpsHQTbOos/Lxe/3E4vFyGQyDA4OIooiSqUSh8MhzYrKZDIIgsD29jbDw8Oy32H3+SKda1OLxSJarVbq+hcEgXq9jsfjod1uo9frEQSBXC6HzWbrtth9vgHZOAuz2dxtEY6U119/vdsifAGTyYRSqUQURXK5HO12m1QqxalTp8jlcuTzeXw+nxSu2tzcfGYvTK/pd5iMjIx0W4Qj5Zts1xlVXyqV0Gq1nDp1Srpj+/Tp00SjUQRBQKvV0mg0UCgUPdeQB/JemwdBNs7ieR3C33WWl5e7LcIXEAQBl8uF2+0mHA7TaDRwuVx89NFH+Hw+jEYjmUyGer2OwWDAZDIRj8e/8fN6Tb/DRO73IXyT7TpjPqxWK263m/v372OxWJiYmGBnZweNRkO9Xmd9fR2n04nX633Jkr8Ycl6bB0E2zkLuccVevG2tVquRz+f53ve+Rz6fJ51O43K5yOfzaDQaqtUq5XKZSCTC9PT0M+/j7kX9Dgu5b2S+yXbhcJhXXnmFaDRKPB4nm83i8/nY29vDaDRKJ9Dvfe97lEqlFxpl3w3kvDYPgmycRZ+Xj1Kp5Ny5cywsLGCz2fD7/bTbbc6ePYsoing8HsrlMu+88w4PHz7E7XZ3W+Q+LxGn08nS0hJvvvkmtVqNubk5tre3uXTpEu12m4GBAex2O8vLy5w4caInk9t9/h3ZOAu591mcOXOm2yJ8hcnJSZaWlrDZbGi1WkZGRqjX66hUKux2O61WC6vVyr/+67+iVCqfefrrRf0OC7nPFvo624miKDXY3bx5E6PRSDabZXR0FJVKRaPRIBAISNesFotFxsbGXrboL4Sc1+ZBkI2zkPuYgGKx2G0RvkIikZBeAkNDQ3zyyScAUk9Fs9lkYGCAyclJUqkUW1tb31ji3Iv6HRb1er3bIhwpX2e7QqFAqVQikUjg8XiwWCw4HA5qtRpLS0solUru3LnD5OQkWq2WbDbbs7kdOa/NgyAbZ9GLYwIOk62trW6L8BXUajVutxudTsf9+/dxOBwMDg6yv79PMpmkXC6zvb1No9FgeHgYq9XK3t7e135WL+p3WPTaNaGHzdfZbnd3F6VSycTEBCaTiWQyyd7eHplMhmKxiE6n4/Tp0zx48IBKpcLAwEDPVjTKeW0eBFW3Bejz3cVut0uxZ6VSSTQaJZPJMDk5STKZRKvVotfrKZVKVKtVjEaj7HfZfZ5iMBgolUo4HA7a7TZ2ux2DwUAqleLMmTNsb2/z6NEjLBYLXq+XSqXC4OBgt8Xu8wxkc7KQ+7C20dHRbovwFQYGBsjn85RKJcbHx7Hb7ej1eh4/fky9Xsfv9zM8PIxer2d5eZlAIEChUPjaz+pF/Q4LuTeafZ3tcrkcgUCA9fV1tFotw8PDuFwu9Ho9i4uL1Go1PB4PExMT5HI51Go1Op3u5Qv/Ash5bR4E2TgLuU+ddTgc3RbhK+h0OjY3N9nc3MRisVAqlSgWi7zyyitUq1VyuRw3b94kn8/j9XopFArfGJLpRf0OC71e320RjpSvs102m0Wn06HVasnn8/ziF7+QuvovX76MUqkkmUxitVoJBoMEg8GebMgDea/NgyAbZ9GrNdqHxfz8fLdF+AqCIOB0OnG73Tx58gSdTsfJkye5f/8+IyMjtFotfD4f5XKZ6elpNjY2vnFOUi/qd1jEYrFui3CkfJPt1tfXOXbsGPl8ntOnT7O+vs7MzAwrKyuMjY1htVr55JNPCAQCPTlAsIOc1+ZBkI2z6NMdfD4farWabDbL8PAw29vbnDp1img0SqFQwGazYTKZmJ+fx2Qy9eQI6j6Hj1arxWazsbi4KIWIbTYboVCIsbExdnZ2GBoaQqlUsr+/L/vyYjkgG2ehUsk7V2+327stwtciCAKiKDIxMUGr1cLtdrO+vs7o6ChKpVK6WnV2dpZ0Os3KysrXfk6v6ncY9Gos/rD4su3q9TrBYJBkMonH42F4eJhisSjd0b6xsYHNZkOlUknrpJcn88p5bR4E2bxh5R4X7tXGoEQiwczMDE+ePKFcLmMymaSJtHa7nVqtRiQSoVarMTo6SrVapVAofOX2sV7V7zDo1ZlHh8WXbRcMBtFqtQwNDbGzs0MqlUKn0zE8PMzOzg4jIyNEo1Gpr+L8+fPPnBvWbeS8Ng+CbE4Wcm+cuXnzZrdF+FoCgQBra2tEo1GuXbtGoVBgY2MDg8Eg7RYdDgdarZZUKoXP52N3d/crn9Or+h0GOzs73RbhSPmy7crlMmq1mlQqRbPZxOl0YrFYpObNtbU1CoUCc3NzlMtlVlZWetqhynltHgTZOAu506uXO3XmPV25coW7d+8yODiIXq8nFAqRzWaZmppCEAQpr6HX67/2TuNe1a/P8/my7ex2OyqVilqthk6nQ6lU4vF4aDabPH78GKVSidVqJRwOc/r0aQRB6OkEd39tPqXvLL4j9GpZYT6fp9lssr6+TiAQoFKp8OTJE2l43Pvvv4/P58NmszEyMkKpVCIcDn/lc3pVvz7P58u2297elprwrFYrVquVu3fvksvlePPNN9na2iKXyzE0NMSTJ0/Y29vr6SGC/bX5FNk4i14dFXBYvPHGG90W4WsZHR2l3W6TTCbx+/1ks1m+//3vc+PGDTweD2fOnCGVSpFOpxEEgXg8/rVlzr2q32Eg96auL9uuWq0SjUZRKpUkEgmy2Sx+v5+TJ09y584d3nzzTZRKJYIgUC6X0ev1Pd24KOe1eRBk4yzkfmfA4uJit0X4WoxGIwMDA0xMTHD79m2p+c7r9RIOh0kmk+j1eum6VYVC8bXVQb2q32HQy8nbw+DLttPpdAiCQLVaRaPR4HK5AHj06BE2m41isYjL5eLu3buMjo4yOTnZDbFfGDmvzYMgm2oouV9+1MvD6FQqFaIoMj09LTmCzjWZuVyOtbU1nE4n4+PjLC0toVB8dY/Sy/p9W6rVardFOFI+b7tWq0UulyOTyTA2NoYoimxsbCAIAoODg9hsNgwGA61WC7vdTjwe7/mydzmvzYPQ21bq850gHo9Ls30WFxexWCz4/X4MBgM6nY5kMkk2myUSiTA8PEy73aZer/d0nLrPb0csFkMQBEZGRtjY2ECpVKLVavF6vWg0GhQKBRsbGxQKBUZGRrBYLLI/eckF2YSh5H750dzcXLdF+EYGBwdRq9Xcvn2bK1eukM1m2draotFokEgkaDQaDAwM4HQ62dvbY3R0lGAw+IXP6GX9vi0+n6/bIhwpn7ddMpnE6/WyubmJ0WhkZGQEtVpNNBpFFEW2trak8R87OzuUy+We//eR89o8CLJxFnK//CiTyXRbhG9kZGSEZDLJ6OgohUKBRqOB2WwmHo8zNTWFx+OhVquRSqWkvMWXmyh7Wb9vi9zzaZ+3ncvlQqlUUq/XyWaz5HI5HA4HJ0+eZHt7G4vFQrlcJpPJMDQ0RLlcZmhoqIvSPx85r82DIBtnIffLj3q5sSuZTHLt2jV2d3d59OgRQ0NDbGxscPz4cSKRCIlEgna7zYkTJxgaGmJ3d/croYde1u/bksvlui3CkfJ52+3u7hKJRBgaGmJmZgZBEEgmk6yurjI3N8fm5iYjIyNEIhGWl5e5evUqqVSqe8K/AHJemwdBNs6iT/fQaDTUajXGxsaw2+3kcjkGBwdZXFyk2Wxit9sxmUyUSiX29vYQBIFSqdRtsfscAYVCAY1GQzQapVgsSkMCPR4Pn376KU6nk1QqhdVqZW5ujng8jtVq7bbYfV4A2TgLuV9+ND4+3m0RvpFOWWyj0SAajQJPS2qPHz9OsViUxj44HA6uXLlCs9n8SoVJL+v3bZH7ILrP267TpHn58mUp0Z1IJNjZ2WFgYIChoSEMBgOhUIhUKkW1WsXv93dR+ucj57V5EHrSWQiCMC4Iwn8VBOGtz/773I4duV9+1MtNhy6Xi3w+j1qt5p133mFychKVSkWj0eD06dOoVCoSiQR7e3tS7uLLV2j2sn7fFrlXfXVsJ4oiPp+PZDLJ/v4+giCwublJvV7n3Llz2O12CoUC4+PjvPXWW6hUKgRB6PmNnpzX5kE4kLMQBGH0aMT4Cv9dFMW/EEXx18D/BP7b835A7pcf9XJjkE6nIxKJkM/nSSaTlMtldnd3pT8DAwMMDAxQLpf5+OOP8fv9hMPhL4yl7mX9vi17e3vdFuFI6dguk8lQKpUIBAIsLCxI97B7PB7K5TJLS0sUi0X29vYol8tUq9WvHSrZa8h5bR6Eg54s7IIg/O+dLwRB+E+CIHzvMAUSBOELZz5RFLeAPzzM39Hn8PF6vZjNZj7++GNKpRKjo6PE43F0Oh2hUIh8Po9Go5FGgvj9fiKRSLfF7nOIRKNR/H4/oVAIu93OwMAASqWS3d1disUiiUQCk8mE2+3m9u3bGI3Gnp422+eLHMhZiKL4AJjsOAhRFP8X8I4gCH90iDLNAV8pH/myE/kyvd4F+m1xOp3dFuGZDA8PE4/HmZ6eplqt0mw2uX79Ou12m/Pnz7O/v49arSYej0t3XXw+yd3r+n0b5H7XSsd2giCgUCioVqvs7e1Rr9dJp9NcuXIFgB//+MfS3x87dox4PM7IyEg3RX8h5Lw2D8JBw1C/BETgviAIZz/76/8L+ItDlOnrbkfPAM/MW8j9gTx58mS3RXgmlUpFSlSWy2Xq9TofffQR5XIZURQ5d+4coiji9XoZHh5maWnpC7H8Xtfv29DL47cPg47tNBoNCwsLDA8PMzw8DMC5c+ekNXDnzh28Xi9ra2vs7e0xNzf3nbiHRs5r8yAcdDs+AfwXURTzwlOswP/47M9LRxCEPwb+GJ4moX76058CT6tPNBqNFCvW6/UMDAxIXcMKhYLh4WFisZjUn+H3+ymVShQKBeBphU8nMQtPq3ucTiehUAh4mlAfGhoiGo1Sr9eBp53MhUJBegCcTicKhYJkMgmAyWTCZrNJcVqVSkUgEGB3d1eabRUIBMjlctKu2+120263icViaDQazGYzFotFCuF0QjvhcFhqTBweHiadTlMul4GnL6tmsyk1F1ksFkwmk1S5pNVq8fl8hEIhKY/QqWLpNJR5PB6p0Qqe3qes1+uJxWLAv8/majabiKKIQqGg1WrRarWkU59Op0MURfb391Eqlfz85z/H5/NJu02r1fqdt1M6naazHjt2un//PrFYrCfspNPp8Hq9X+gd6IQMOzOsfD4f+/v7Un/I856nzklhbW2NdrtNq9VCr9dL9oenfVBKpVJaCx3nUavVUKvVPWEn+PrnSRRFTCbTd95O3/Q8vSjCQS72+CwU9KfAn4miWBAE4T8BY6Io/t8v/CHP/x1/wFOH9Pbn/i4LnP8sf/G1TE9Pi6urq4clRs9x48YNrl+/3m0xvpH5+XmSySSiKFIul8lms1gsFqampnj48KFU9eL3+9FqtYTDYURR5A/+4A+A3tfv2/DTn/6Uv/7rv+62GEdGx3Z/+7d/i0qlYmRkhEKhQK1Wo1AoUKlUOHbsGPv7+8RiMSqVCqdOnaLVauFyuZidne22Cs9EzmsTQBCE+6IoXnje9x00Z7EliuL/KYpi4bOv/xeQP+Qk9zxfE4p6lqMA+V9Q0us5Ga/XS6vVwmQy0W63OXbsGEqlkvv373P27Fm0Wi2lUolkMintMD8/5qHX9fs2fN2UXTnRmTo8ODiIVqtFrVZTKpXIZDLU63WuXLnC1tYWmUyGyclJ/H4/zWaTer3+nUhwy3ltHoRvvYpFUfwfwPYhyNL5vC84hc9OM3/3vJ8zmUyHJUJPcu3atW6L8EwKhQJarZZIJMLMzAyxWIz9/X3MZjOffvopFouFoaEhisUit27dQhAEcrmcdGVlr+v3bejE7+XKtWvX2N/fl8I2t27dIpfLYTAYOHbsGHfv3kWv19Nut4lGo0xMTEgXYvV6jwXIe20ehEPZ8oiieGjO4jP+c6cpD/gDURT/y/N+QO59FvPz890W4ZmMjY1RLBZRKpXcvHkTrVaL2WxGqVRiNpuJRqNSrHR8fFwqSOjEfXtdv29DJw4tV+bn5wmFQpjNZrRaLWNjY1SrVRqNBsvLy+h0OgqFAsPDw+h0On7zm99IwyS/Cw1vcl6bB6Enz8efhbv+QhTFX4ui+EKVVnKfOttJ6PYqWq0Wl8tFvV7H7Xaj0WjQarWMj4/jcDg4e/asNCuo0Wjw8OFDpqampMRir+v3bZD7kMtCoUCpVMLhcPDkyRN2d3exWq3s7e1x7tw5vF4vFy5coFgsSgnkz1fP9TpyXpsHoSedRZ/vJhqNhnQ6zejoKIlEAr1ez82bN1Eqlej1eq5cuUK5XEatVnP69GkePHiAw/F1ldJ9vmu43W4++OADJicn8Xg8FItF3nzzTemK1Rs3bmA2mwmHw/j9fur1ej8X8B1DNtYyGo3dFuFIuXDhucUKXSebzRIIBIhEIthsNur1OuP/f3tv1hzXdd57/zZ6Qs8juhuNqdEACAIgCRKERFKiZsmK5Ypt2crJTS5ykTe5z0VS5xOciitfIH5vcps3ubBdlUSWRJuSSJEUARAAAWKeG2ig53nu3u8FiV20j20OJtLc2/tXxapGswk+q/da+7/Xs54hFGJ/f5/t7W3Jh+10Otnc3MRqtRIOh+nu7pbF+J4XuTxBPy+Tk5PMzs5it9vZ3Nzk6tWrNBoNHjx4QL1ep1KpMDg4SDgc5syZMxweHuJwOGRTul3Jc/NZUMzOolartdqEE0UO9YUMBoOU4zA4OMju7i6xWAyv18vQ0BDlchmn00kkEqG/v1960gR5jO95UXo59qOjI7a3t2k2m7z++uusr6/jdrupVquYzWbOnj1LOp2W8hSSySRut1s2ibRKnpvPgmLE4jiRR6kc31RfZvr6+igUCly8eJHZ2VksFgsXL15kdXWVRCLBwMAAiUQCj8eD1Wolm81Kbig5jO95UbrPe29vD5fLhSAICIKAy+UiHo/jcDjw+/3cuXOHgYEBXC4X4XCYUChEpVKhu7u71aY/FUqem8+CYsRCpfXs7e1htVoplUpMTEyg0WhYX1/njTfeoNlsEo1G0Wq1JBIJ7t69K2Vtq8gbURTxeDxUKhXu3LlDIpGg2Wyi1+uJx+O89tprHB0dUa/XOXPmDKVSCZvNJmU8q8gDxYhFe3t7q004UYaGhlptwhNxOp0kEgnMZjO/2//VjgAAIABJREFU+tWvCAQCCILAzZs3GRwcZHNzk0wmQ7PZxOVyMTQ0xOHhIcViURbje16UfojvcDhoNBr09/fj9/vJ5/OIosjy8jI2m421tTXq9TqBQIDPPvsMvV5PLBajs7Oz1aY/FUqem8+CYsRC6Rncckhe8ng8BAIBIpEIXq+X3d1dvF4vb775Jvfv3+fP/uzPKBaLZDIZ+vr6uH37NqFQiJ2dHVmM73lRetRPqVSio6ODO3fu4HA4aG9vJxaL8f7773N0dER/f7/Ue72zs5NUKoXZbJaNWCh5bj4LihGL4wJdSmVhYaHVJjyR4yfora0tzp07R6FQIBaLEY/HOXv2LPPz83R1daHT6dja2uL06dPs7u6i0WhkMb7n5bjIoVLJZrPcvXuXc+fOkc/nKZfLBINB7t+/z9jYGNvb2zQaDbLZLGNjY+zt7WG322Ujokqem8+CPK6WiiwwGAwcHR0RCoX49ttv6e3tpa2tjaOjI6nqqsFgIJlM4nQ6KRaLUqavinwRRRGdTsfq6iqXLl0iGo3i9XqpVqssLCxQr9c5PDzk3LlzLC4uEgwGicViivcGKA3F7Cx0Ol2rTThR5NITQRRFHA4HoihycHCAVqulvb1der25uYnf72dnZwedTkdfXx8PHjyQzfieB6XnAOXzeYxGIxcuXODu3bv4fD52d3dpa2tje3sbi8WCXq+XSoU7HA5ZFVdU8tx8FuRzxZ6A0v2Kp06darUJT8Vxjwav14sgCGxubtLf3893v/tddnZ2pHr7o6OjxGIx9vf3cbvdshnf86D0TmsWi4Vqtcr8/DzBYJBkMinV/PrOd77D4OAg2WyWg4MDXC4XxWLxNyoOv+woeW4+C4oRC6UnPt24caPVJjwVR0dHWK1WQqEQZ8+eRaPRMD8/z97eHnq9nkajwcDAAKVSicHBQVKpFBaLRTbjex6OGzEpkUajQbPZJJlMMjIygiiKeL1edDod8XicWq3GwsIClUqFV199lYGBgd9o0CMHlDw3nwXFiIXKy4HFYkEURRYXF0kkEnR0dGCxWNje3qavrw+NRsP9+/fJZDLEYjGpUY6KPDk+e+ju7iYWi5FMJllaWsLj8eD3+1lbW5N6l9y7d49MJkO9XpdFtVmV30QxYiEnH+jzIBc3m9FoRK/X4/F4yOfzmM1motEor7/+Ovfv38flctFsNslmsxiNRvL5PFarVdFVgzUaTatNODGOjo7QarVkMhlMJhO1Wo1sNkswGGRlZYXh4WHy+Tx6vR6j0Ugmk8Hr9cqm1AfIZ+2dNIq5wyr9EPHKlSutNuGp8Pl8BAIB7t69S1dXF7u7u4RCIUwmE8FgkP39fXp6ehBFkVqtJvXqdjqdrTb9xJCTf/5ZqVarVCoVOjo6aDQa1Go1Lly4wOzsLMFgELPZTF9fH5ubmwwPD7OysoLRaJRNjgXIZ+2dNIoRC6U3P5qammq1CU9FrVajUCjQ29vLwsICbrebZDJJJBKhra2NZrNJPp+n0WhQrVZxOp3s7u6Sy+VabfqJoeSyFg6Hg1qtRiwWk3aI6XRautbLy8vUajV6enqYnZ2lq6uLg4MDWe0s5LL2ThrFiIWS3RggnwP8QCDA7u4u7e3t6PV6tre3sdls3Llzh2w2i8Viobu7G4fDQTweJxqNcurUKdmUq34elFzkcnZ2Fp1OR0dHBw8ePMBmsxEIBNDpdOzu7lKpVKhUKuzt7SGKIhaLhVwuJyuxkMvaO2kUIxYqLwc6nQ6z2czh4SEjIyPY7Xb29/f54IMPSKfTkn+7VCoxOTnJ0dGRFJOvIj+sVivNZpP5+XnGxsYQBIGDgwMcDgd2u52hoSESiQRarVZ6rfRaWUpFMStU6WcWly5darUJT41Go6Gnp4d6vc7Y2BgGg4GlpSUuXLhALpcjn8/T39/P4uIiIyMjkgtDqXR1dbXahBOh2Wyi0+lIJBKcPXuWzc1NdDodFouFer1OX18f29vbiKLI6OgoJpOJjo6OVpv9zMhp7Z0kihELpZeMCIfDrTbhqWk2m9LTpsVi4dSpUxiNRm7duoXD4aCnp4disUhvb+9vuGiazWYLrT45lBoanM/n8fl86PV6arUanZ2dtLW14fF4MBgMzM/Po9FoOHfuHHq9nvv378vyGstp7Z0kihELJfuFAfb391ttwlNzfF7R399POBxmZmaGcrnMhQsXiEajbG5uks1mmZmZQRRFzp07Ry6Xk1Wi1rOg1MP7cDhMuVzG4XBQKBS4d++eVEhwb28Pp9OJx+Ph1q1b7O7u0t/fj9lsll07ATmtvZNEMWKh8vIQCoVwu90sLCwQiUQIhUKYzWZ8Ph/nzp0DHt5AOzo6EASBvb09NBqNYsVCqRQKBZLJJM1mE6PRSCAQoNlssru7y4ULFxgeHiYWizE4OIhWq+Xbb7/F4/HQ29vbatNVngPFiIXcnlaeleHh4Vab8NQcHh6STqc5c+YMDoeDu3fvks/nOTw8lOLs9Xo95XJZcmXU63VMJlOrTT8RlFobymw2UygUaGtrQ6PRUKvVqNfrOBwOZmZmiMfjGAwGpqenqdfrdHd3I4oi8Xi81aY/E3JaeyeJYsRC6eWO5ZQF3N7eTrlclqKf3nnnHaxWK1999RUWi0XqjGez2aQcDIfDodjtvlIjvXZ3d3G73YiiyM7OjhTxVCwW8Xq93Lp1C1EU+eCDDygUCpJQyE085bT2ThLFzGKlNz968OBBq014ajweDxqNRipJfdw57cqVK8TjcU6fPs3GxgZarZaJiQkikQipVEqxQQqxWKzVJpwIjUaDTCYjib/T6SSdThMMBkkkEly+fJlEIkEikaCrq4t8Pk8sFsPn87Xa9GdCTmvvJFGMWKi8PLhcLkwmEy6Xi0ajgSiK6PV6KpUKIyMjbGxsSK6Lvb09BgYGqNfr6hOcjBBFEa1WSzwex2g08uDBAzweD9lslnA4zKlTp6jVajSbTbq7uwmHw7hcLrq6umTTIU/lN1GMWCi9+ZHf72+1CU+NXq8nkUhIGbt2ux2dTke1WmVtbQ2r1UogEKBSqVAsFnG73eRyOcmloTQsFkurTXjhlMtlOjs7pYKQHo+Hra0turq6MJvNbGxsUCqVsFqtmEwm7HY7zWZTlpn6clp7J4lixELplSH7+/tbbcIzodVq6e7uplqtsrOzg8/nkxLzKpUK09PT6PV6gsEgR0dHdHV1EY/HFZmT4HA4Wm3CC2d/f59KpcLp06cpl8v4/X6MRiNzc3Ok02nGx8fZ2dnB4XCwu7tLo9EgEAig1+tbbfozI7e1d1IoRiyUXr/l1q1brTbhmRgeHiYSidBsNonFYqRSKUKhEBsbGxwdHfH2229TqVS4efMmtVqNUqmEXq9XZAKUEseUyWRIp9MYDAbq9Tp3794lmUzy1ltvUa1WWV9fx+12Uy6XSSQShMNh2traGBoaarXpz4zc1t5JoRixUHm5iMViZLNZtFoto6Oj3Lt3D7/fT7PZ5M0332R6ehqdTsepU6dwOBxoNBrq9TqZTKbVpqs8Bblcjkajgd1ul8q7dHR0sLCwwPnz56nX6wwODjI9PY3T6ZTKgUSj0VabrvKcKEYslBqeeIycqnTCwwNQo9EoCcaPfvQjjEYj7e3tfPHFFzgcDqxWK+VymZ2dHarVKkajUXZhlU+DEg90PR4PVquVnZ0daffg9/vR6XTcuHEDnU6HyWTiBz/4AXa7nVgshkajkWUQg9zW3kmhmDusWkjw5UIQBEKhEOVymbt375LNZjk8POTevXv09/djsVhIpVLk83kuXbrE7u4uR0dHinTZdHd3t9qEF044HCYWi7G2tsbIyAhWq5W1tTVMJhNdXV3s7OywvLxMtVpldXVVaqErR7GQ29o7KRQjFoVCodUmnCh37txptQnPRHd3N4VCgb6+PiqVCpFIhHw+z1/91V8xODhILBbDZrPR1dXF3NwcwWCQcrlMvV5vtekvHCUKoCiKlEolRkZGWFpaor29nYGBARKJBKOjo3z88cdUq1X29/dJpVL4fD4EQaCvr6/Vpj8zclt7J4VixEKO1SyfBbklHQqCgMFgIBwOc/XqVba2tvB6vTx48IDZ2VkcDgdGo5GDgwMEQUCn05HL5TCZTIoLn1WaAJbLZaxWq5SprdVqSSaTHB0dMTAwwNTUFDMzMwSDQTY3N3nnnXdIpVIkk0lZegDktvZOCsWIhcrLhd/vZ29vj7Nnz3Lr1i1Onz5NJpNheXmZvr4+BEFgZWUFrVZLMBikWCzidDppa2tTfAVhuZNMJrFarTidTmZnZ9Hr9RiNRqn214ULFzg4OCAWizE0NMTNmzc5d+4cxWJRlqGzKg9RjFgoMfHpceTWNF6r1dLX18fKygomkwmbzUYkEuF73/see3t7VCoVafexs7NDKpXC4/EQDofZ3d1ttfkvFKWdWcTjcSlju729nVQqRb1ep1gsYrPZiMfjvPXWWxweHlKr1RgYGGBpaYlQKNRq058Lua29k0IxYlGpVFptwomytbXVahOemUqlQi6X49SpU1Lpj1/+8pdSj2afz0cikaBcLjM4OMj29jahUIiDg4NWm/5CkWPW8h/i8PCQYDBIJBLh9OnT6PV60uk0HR0d6PV6rFYrX3/9NdVqlZGREak6rVx3jHJceyeBYsRCqUXojjk8PGy1Cc/MsU97fX2dw8NDLl26RH9/P/fu3QPA6/UyOTlJPB4nmUxKdaOU9iSutITR7u5uqZnR7OwsgiDw+uuv43A48Hq9XL9+nc7OTl5//XX29/fZ29ujo6NDtkEoclx7J4FixELl5aNWq9Hf308kEiGTyTAzM0N/fz/j4+McHR2xtLREPp+nUqlgMBhYXl6mVCoptkqrUkgmkxQKBebn5+ns7EQURSlbf2FhgbGxMcbGxpieniadTkvC0mg0Wm26yh+BYsRC6Ykzo6OjrTbhmRkZGWFra4t3332Xnp4eBEGgUCgQCAQQRRGTycTc3Bw9PT1Uq1Wq1apUvlxJdHR0tNqEF8pxKY/jhmNtbW2srKxQLBZpNpt0dXWRy+UQRZG+vj7efvtt9vb2OHXqVIstfz7kuPZOAsWIhdLCLX8bOT6VJRIJ4GHo4e3bt3G5XGxsbBCJRLh8+TIajYbe3l6y2SyZTIbx8XF2d3cxm82KCoVW0lhKpRJ2u51cLseZM2ekPIrBwUF0Oh2vvPIKqVRKEpTPP/8crVZLrVYjlUq12vznQo5r7yRQjFiUy+VWm3CirKystNqEZ6a9vR2DwcD6+jr9/f34fD4KhQL379/n2rVrlEolzGYzly9flqqTbm9v09nZyfXr1xUT334smnInmUxy+/Zturu72djYYHt7m66uLgKBgBTVtry8zOLiIplMBq/Xy5kzZ1haWqJSqeD1els9hOdCjmvvJFCMWKi8fHR2dpJIJBgbG5MS8jweDxaLhc7OTprNJi6XC41GgyAIlMtl6ZA7Eolw7do12T6NKo2DgwM+++wzwuEwyWSS0dFR6vW6FFgSCAQwm81YLBb0ej0dHR0sLS3R2dlJMBikUCgozh33p4ZixELpyT5dXV2tNuGZsdvt9Pf3c+3aNXZ2drh48SLxeJwLFy5QKpXQarU0Gg3W1tYolUrodDq0Wi1nzpyR/OBTU1Oyj0axWq2tNuG5EUWR1dVVNjY2pJBmp9OJzWZDo9FQLBYpFApSSY+joyNeeeUVYrEY/f39pFIpvvzyS0ZHR2VZFwrkufZOAsWIhdI75ckxnFSr1ZJIJPjhD38IwLfffss777zDr371K8bGxqRqpNFo9DeqBs/OzlIoFLhy5Qo7OztMT08zPz8vW9+/zWZrtQnPRaFQYHp6mnA4zOrqKq+++ir1ep2pqSmKxSImkwmDwSBFR6XTad5++22uXbvGhQsXSCaTpNNpPv74Y5LJZKuH89zIce2dBIoRC7nGcD8tci1mZjab+fnPf47JZGJ8fFzqwz0/P4/P56NarTI8PIzVaqW/v590Oo3NZkOn07GwsEBbWxvJZJL19XWuXbtGsVhs9ZCemf39/Vab8MxEIhF+/etfk8lkiEQiADx48ACNRkMgECCbzeL3++nu7qarq4tYLMbly5f59a9/zdDQEC6Xi46ODsxmM//1X/8ly5pQx8h17b1oXjqxEARhQhCEf3j0598FQVBeT8o/Iex2OwMDA6RSKXZ2dqhUKoiiyMjICEajkUQiIbXdvH37NhcvXmRubg673U5XVxc2m41qtUqhUMBisXD79m3i8Xirh6VoVldX2drawuFwEI/HqVQqOBwOqb/23Nwcly5d4tatW+Tzeen6aDQaqRNesVgkHo9zcHBAMBhUZJ+SPzVeKrF4JAyToij+RBTFnwD/Blx7mn8rV3/o0yLX2lf1ep3d3V0uX75MNBplaWkJn8/H119/TaPR4LXXXsNms5FKpbh8+TKrq6t897vfpdFosLKyQjAYpFarkcvlmJ+fJ5/P8+2338qq7LecztPu37/P4uIihUKBjY0NUqkUtVqNvr4+lpaWAPj+97/PrVu3ePvtt7l3756UwV2tVrl79y46nY79/X329/d5++232d/fl3W0olzX3ovmpRILYBL4x8d+/gKYeJrdhclkOjGjXgYmJydbbcJz0Wg0eOWVV2g2m7z//vsYjUYWFhYYHh7miy++wGAw4Ha78Xq93L9/n66uLjY3N0mn07zzzjtMT08zMTGBRqOROrPt7++ztbXF8vKyLMp/BwKBVpvwRMrlMlNTU1I+xP3799HpdNTrdS5dusS9e/d4++23yefzPHjwgIsXL3L9+nXOnDlDIBBAEAS++OILBgYGODw8pFQq8Z3vfIdms8n58+cRBKHVQ3xu5Lr2XjTCy5bMJgjChCiKM8evgWlRFH/nTBME4W+BvwXwer0X/+3f/g2AUCiE1Wplbm4OALfbzdjYGF999RXw8OD16tWrzMzMkM1mgYcT4ujoiL29PQCGhoYwGAwsLCzw6Pdz6tQpbty4AYDBYODKlStMTU1JtX8uXbpEOByWfNTDw8NoNBoePHgAPCzb3d/fLzWANxqNXLp0iTt37kg5BVeuXGFra0uKABodHaXRaDAzM4PZbKarq4vu7m7Jj2qxWJicnOTWrVtSMcWrV6+yuroq9Ts+c+YMlUqFtbU1AHp6evD5fExNTQEPD2AnJia4ceOGdPN98803WVxclHIExsfHyeVybG5uAhAMBnG5XMzMzAAPI2TGx8f58ssvEUURQRB46623uHbtmtSvOZPJ0N3dTTqdRq/X02w2sdlsZDIZNBoNJpOJVCqFRqPBaDRSrVaxWq0UCgW0Wq1UiE6r1ZLP57Farej1ek6dOsX29vZLcZ2OY/Ifv04//elP+fu///uX9joNDAxw7949tFot6XQai8VCs9mkXq/T3t4uXYd8Po9OpyOfz0vXrlarIYoiOp2OWq2GTqejWCzS1tYmBS0IgkBXVxehUOilWU+/6zr9vvWk0+lwOp0tv05vvfUWc3NzUjj5xMQEyWRSmvvPe9+zWq3Toig+URFfOrF4HEEQ/glAFMV/fNJnh4eHRSUnz1y/fp2333671WY8M/v7+8zMzFCpVKTueB0dHdRqNaLRKM1mk0AgQLVaRRAEjo6OsFqtaLVafD4f0WiUhYUFfvzjH7O2tobFYmF3dxeXy0UkEqFWqzExMUF/f/9LG7Xy13/91/zrv/5rq834nSwtLZHJZKTif06nU2piFI/HOXXqFD//+c8ZHx+n0Wjgcrk4Ojoil8sRCAQ4PDyks7OTSCSCKIpYrVb6+vpYXl6mVqtRrVa5cOGCFP0mR+S69p4WQRCeSixeNjeUxCPX08TTCIXKy0tbWxvNZhOr1YrFYkGn00ld8c6cOcPo6CjNZhO9Xi91UtNoNJRKJRKJBPV6nffee4///u//pq+vj729PZrNJh6Ph4GBAbxeL6urq+zs7EhPeipPx9TUFIlEgu3tbYLBIJ2dnZRKJQKBAEtLS4yMjHDt2jUuXbpEKpXC5XKxvLws5VkcHh7y2muvUa1WGR0d5cKFC9TrdZLJJBaLRdoJx2IxWUdDqTxE+z/xnzxyFw38gY98LoriF7/13j8Bf/G0/4fSD6GuXr3aahOeC5/PRzKZJJvN8t577/HZZ59hMBjQ6XTcvXuXV199lXK5TCwWo729na6uLqLRKGazmVqtRiKRoFqt0tfXx9bWFhqNhtOnT9PW1katVmNoaAir1cr169c5f/48169fJ5/Po9frcbvdrK2tYbPZcLlcpFIpyZ0VCARYX1/HaDRitVrxer1sbW0hiiI+n49wOIxWq5UELpfLcXR0RFdXF4IgUCwW0Wq19PT0cHBwQKVSwel0YjAYSKfT1Ot1aQdVrVb5xS9+QXd3N5ubm1Kms8FgIJFI0Gw2pd8Ti8Xo7e3FZrOxv79Pe3s7Wq0WvV5PPB7HZDJJv6fRaODxeKhWq+zt7eH1evH7/USjUanz4LHw5vN5+vv7yefzpFIpGo0GoVCI2dlZrly5Qq1WY2tri8HBQQqFAg6Hg8XFRXp6ekgmk5RKJaLRKB0dHUSjUU6fPs3q6iqffvopwWAQURSZn5/H4XBIv+ujjz7iyy+/xGQyodX+j9xqTgS5rr0XzUvphhIE4R+An4qimBYEwSGK4hO7x5w5c0Y89ocqkQcPHsi2+uUvf/lLdnd3GR8fx2az8cUXX3Dp0iUWFxelIm1er1c629Dr9cRiMTweD6FQiEqlwvb2NtFolIGBAer1Ovl8nnw+z+TkJDdv3qSzsxOHw8Hq6ipvvPEGy8vLmM1mqb/3/fv36enpIZVKMTExQTgcxmg0IggC6XSaRqNBV1cXbW1tZLNZNBoN8XicWCxGZ2cn8NB3LQgC2WxWKlMyOzvLxYsXsdlszM3NEQwGSSaTnD59mqOjI9bX17lz5w4/+tGPyOfz+Hw+nE4ni4uLeL1estksY2Nj3LlzR6rkGo1G0el0XLp0if/8z/8kGAyyt7fHe++9RzgcZmtri7GxMbRaLeFwGJ1OJ50rPHjwgLGxMTo7O/nmm2/o7e0lGo3yxhtv8Ktf/QqTycTg4CArKyu43W7K5bJUEbZcLlOtVqlUKpKYms1mSZy3t7eJxWL4/X7K5TIajQa3200qlZLynAYGBtje3ua1115jZWWFSqWCz+fjjTfeaNn8+2OR89p7GmTrhhIE4RPgPx4TiPef5t8pvfnR8eGaHKnX64yPj9Pe3s5XX31FX18fCwsLGI1GjEYjoigiiiIGg4FQKITFYsFms3FwcCD9GRoaYmBgAFEUSafTOBwOnE6ntNsYHR1lZmaG8fFxdnZ20Gq1uFwunE4ny8vLjI6OYjab6e7uZnp6moGBAXw+H6urq1itVoLBIPl8no2NDfr7+1lbW0MQBEZGRsjlckxMTOB2u1laWuLq1auEw2HK5TLvvPMO0WiUg4MD3nrrLaLRKHa7HZ1Ox/b2NlevXqXZbJLL5fD7/WSzWWZnZ3nvvfeYm5ujv78fgFgshslkQqPRoNFo6Ozs5PPPP+fcuXM0m00GBweJRqNoNBoGBgYwGo3cu3dP2vlotVpWVlZ4//33WVtbo9FoSE//xwUcJyYmsFgsUmmOtbU1enp6sNvtJBIJ0uk0VqsVl8tFsViku7ubsbExNjY2SCQSHB4eYrPZMJlMBINBzGaz1B7XYrHg8XiIRCK43W6++uorvF4vg4ODsi8IKee19yJ5qcRCEIQQ8O/AhiAIoiAIIg/dUSoyRqfTMT8/z8LCAhMTE2SzWbRaLcVikXw+j9vtlrrqCYLA/v4+Pp+Pnp4eYrEYoVCIlZUVcrkcyWQSj8eDw+FAq9VSr9cJhUJEIhF8Pp/UkKfZbLK7u0s8HueNN96gWCySy+Xw+XxcvnyZnZ0dFhcXGR4eRhAE7t69S1tbmxSl4vP58Pv9ZDIZent7mZqaIpVK8cknnzA1NUV3dze1Wo2joyN6e3sxm83Mzs4yMjJCb28vCwsLXLhwgb29PTQaDaFQiEKhQLPZpKOjg5mZGT7++GP29/dZW1vj448/JpVKcXBwgM/nY3d3l7feeotqtUosFsPn80lP+vl8nvX1dSYnJ8nlchweHlIoFOjq6mJpaYlgMMjR0RHRaJRXXnmFtrY2yuUy29vb0ve0u7tLvV5Hp9PR0dFBe3s7zWYTh8OBxWKhUChIYbLDw8PSjsJsNpNKpTAYDLhcLpaWljAYDGSzWRKJBJVKhWQyyfDwMEdHR9y7d09WeSYqv5+XSixEUdwURVH4rT9/6KxDQunNj86cOdNqE56b46fOSqXCwsICbrcbu92O3+/HbreTTCYxmUwsLy+zv7+P2Wxmf38fg8HA6OgoN2/exGq14nA4sNls5HI5FhcXCQaDVCoVarUaRqORV155Bb/fz97eHrVajWKxiNVq5fbt25RKJQwGA/Pz82xsbBCNRvH7/cDDQ/jOzk7pfOHYtZPL5RgYGCCXy2Gz2fB6vXzzzTf09fURCoVIJBIIgkC1WkWj0aDX66UkNa/XK7m6Go2G9NlYLMbw8DAdHR1MTU3hcDioVCrMzs7icrmw2+3s7OwwPj7OnTt3MBqNBAIBwuEwk5OTRCIRCoUC3d3daDQahoeHKRaLdHR0UC6XJbdPtVrFbDZz//59uru7pZv/cc8Qu93ORx99RDqdJpvNksvlePfdd1lYWCAcDmMymbDb7bhcLr755htcLhednZ3k83lEUWRnZ4fV1VVGR0fJZDKYzWbMZjPBYJBAIMDOzg7RaBSXyyX7HCg5r70XyUslFn8ML+PZy4vkOOZbjnR2dnJwcMDIyAjj4+NYLBYcDgelUgmHw0FPTw8WiwWLxUKpVKJer2MwGPB6vezv73PhwgW2trbI5/OUSiUpfnxhYQGz2UypVOLg4ICpqSkCgQD7+/sMDAzQ0dHB6uoqgUAAv98vHXB7vV4++OADaSeTTCZJJBJYLBYODw9pb2+np6eHXC7H3t5TWLjGAAATKklEQVQeJpMJs9nMnTt3cDqdpFIprFYrZ86cYW9vj3Q6zcbGBuVyWTo0zmQydHR0sLm5Sb1eRxRFDg4OOH/+PBaLhUwmg8ViYX5+Hr1eT3t7O1tbW4yPjzM4OMjh4SFOp5P9/X06Ozulzw4MDGAymaTcgY2NDelw/viMJhqNkkqlqNfr+P1+IpEIm5ubuN1utre3OTg4oLe3l2+//ZajoyMpJPbmzZvSeUc+n5d6j4yMjNBsNjGZTFIgSaVSob29nY6ODoxGIz09PVLE0/E5x+joKEdHR9KZj1yR89p7kShGLORcTuBpkHNYaCKRwO/3s7q6KhUFPG6cc3BwgNVqJZfLSaU9TCYTZ8+elcI6k8kk77//vhRtVKvVyGQyVCoV3G430WiUUqnE8PAwiUSCYDBINBrF6XTS2dlJMpkkk8kwMTGB3W5nb2+Pg4MDvvrqK6md69WrV8lkMqTTaXp7e3nw4IFUOTUYDCIIAqFQiJ6eHuLxOKVSCYvFgt1uJ5/Pc/bsWXZ2dsjn83i9XjweD8lkEr1ej8lkwu1243K5sFgsUt2kzs5O6T2z2Ux7ezs6nQ6TyUQ0GiUUCmEwGICHZ3JarRadTkdvby8XLlyQvrNisYjBYMDv9zMwMEA6ncZut9PR0UGj0WBvbw+Px4PdbsdkMtHT00MikcBut2Oz2Ugmk1ICXj6fp16vS9FkH374Idlslrm5OQqFAuPj4+j1eiqVCqdPn2ZlZQWbzcba2hput5twOCxlgs/Ozkr2yBk5r70XiWLEQuXl5Tjc0+12s76+zqlTp7BarXz66ad4PB4WFxfRarUUCgUpFPY4k3dwcJCtrS2mp6cJhUL09/dTLpcpl8tSe8+LFy/idruJx+Mkk0kajQbd3d3s7e3R1dWF3W6nVquxsrJCvV5neHiY5eVlenp6yGQy+P1+5ufnqVarnD9/nqmpKekG2Nvby/T0NMlkEqfTye3bt5mcnCQajRIOh+nt7WVkZASHw8GPf/xj9vb26Onpwe/3Ew6HGRwclDKGe3p6pEP0S5cuMTU1xblz5xgcHESr1TI5Ocn6+jpra2tcuXJFqsw7MzNDJBKRoqJisRhffPEFqVSKtrY2KbM5Fouxvr7O4OAglUqFra0tyuWydIC/ubkp9UKPx+PodDoajQaXLl1Cr9dTr9cpl8tkMhmGhobo6elhZWWFra0tzp8/z/LyMtPT03R0dDA0NMT29jYTExNSY6PPPvsMvV7P4OAg6+vrjI6OSqHIKvJHMWKh9EO0np6eVpvw3AQCAbRaLdvb2/j9fmq1GqVSidHRUckvnslkcDgcUmRRJpNBr9fj8Xh49dVXsVqtfPnll8TjcQ4PD3G5XLS3tyOKIrOzs1Ikj8lkwuPxUK/XuXjxIl9//TVut5szZ86QSqUYGRlBq9Xi9/vJ5XLUajXcbje5XI7BwUGsViudnZ0Ui0UymQwul4tGo0G1WqVYLEolOvr7+6Xon8PDQ2ZmZtjZ2ZGyzsvlMhMTE9y9exdBEKT8h1gsJoW4BgIBjo6OmJqaQqPRMDc3h1ar5dy5czx48IBLly6Ry+Wkg/FCoSBlwg8ODhIMBmk0Guzv72Oz2cjn81QqFTo6OhgZGaFQKHD69Gl8Ph/Xrl1jdHQUo9GIw+HAYDBIB9m3bt2ira0Ni8UiicXe3h5ra2uYzWYmJyclUTo6OmJjY4PDw0N6e3u5f/8+oiii0WiYmJigXC7jdrux2WwsLy/T3t6Ox+Np9RT8o5Dz2nuRKEYslN78SM5PZ1arlWQyicFgoL29nevXr2O326UYfbvdTnd3N1NTU+h0OprNJm+99RblcplwOCwdFPv9fkRRxOl0YrFY6OnpYXFxkaGhIUZHR4nFYhiNRjY2NiiVSgiCwODgILFYDIPBwPj4ONPT09JT+nHE0+3bt+nv7+fw8JBUKoXH48Hj8dDd3c3NmzdxOBx4PB4CgQDXr19HEAQSiQQLCwt4vV6pqu5xtNfu7i42m42ZmRnefPNNGo2GVLW1v79fyvcYGxujvb1dKvvt8Xiw2WzSzTsajZJIJKRzjp2dHVwuF729vXR3d9Pd3c3q6irZbBaPx0Oz2aS/v59vv/1W6kjYbDaJx+OMjY1htVoplUpSuHA4HOb06dOMj48zNzdHIBBgZGQEk8lEuVymra1NSkA8zqi/fPmydJh/nFdyfEh/7Ar89NNPpTDger3+G42t5Iic196LRN5X8TGU3vzouEiZHDluv2k2m2lra+O73/0ulUoFvV5Pb28vW1tbJJNJHA6HlEMxOzvL0NCQFOk0PT2N3+/H5/NhsVh48OABdrudH/zgB2SzWWKxmNRJ79ittLOzQzAY5ODggLm5Ofb396VkvoODAyk34bi4Ya1W48aNGzQaDdxuN6urq/T19REIBFhdXZUOzo1GI4uLi3R0dEiuq6+//ppXX32VX/ziF0xOTjI9PU0wGCQWi6HRaNje3mZ0dJTbt2+TTCY5f/48CwsLlEolCoUCo6Oj6PV6qf+D2Wxmd3cXk8lEOBzGbrdjtVrZ2trC5/Nx7949rl+/js/nw2azSaHHx5nWiUSCSCTC3t4ee3t7UvjsrVu3cLvd6HQ6jEYjMzMztLe38+d//uc4nU5mZ2ex2+10dnYyMDDA3NwczWaTTCZDIBBge3sbq9VKNpvFZDJxdHREOp2mp6eH3t5ednZ2+OEPfyhlkB8Xf5Qzcl57LxLFiIXKy00ul0MURVKpFHq9Hq1WKx2EHlegHRgYoFqtYrFYaGtrY21tTbrhvvvuuxwcHLC8vExfXx+1Wo3p6Wk2NjbQarVUKhW8Xi+1Wg2DwcDp06c5ODggm83y0UcfcXh4iN1ul85Ejl07H330EZFIhP7+fvb39/n+979PJBJhZ2eHwcFBurq62N/fp1QqEQqFiMVilMtldDod/f39xGIx2traJEEbGhpiYWGBwcFBqejeMSaTifb2dlZXV6lWqxgMBlZWVqTfdxzZFQ6HKRaLGI1GhoaGpJLfvb29vxFye+7cOfL5PLVajWw2K5UaeeONN6RWtSaTiQsXLlAoFNjd3WVoaEiKdnK5XPh8Pg4PD9nc3GRubk5KHkwmk2xsbPDqq69yeHjI+Pg4h4eH5PN5bDab1PEwk8kgCAK1Wo179+5JJUwMBoOUQKmiDBQjFkpvfiTXPs7HnDt3TnpK1mg0Ujb12toaHR0d6PV6qbZQs9lkcnKSYrFIOBzG6XRKN8VkMokgCLzyyitUq1Xy+Tw9PT2Uy2Xi8bgUGXQcCdTR0cH09DQfffQR1WoVr9fL4eGh5Fq6ceMG7733HrFYjO9973t8/fXXDA8Pk8lksNvtUhLchx9+iN1u58yZM7S3tzM0NITb7ebUqVNsbm6i1+txOp3AQ2E8PkQ+d+6cVMLE5XLx2muv4Xa72d3dxel08sknn3D69GlKpRLxeJzu7m4mJycxGAwEg0EsFguffPIJmUyGzs5OqT2tXq/HbDZjNBpxu91Swb6zZ89SKBT48MMPpfOYcDjMuXPnODw8xO12o9frqdVqxGIxKZnv+Kzi3XfflVxXe3t7mM1mOjs72dzcJJlMcuXKFXK5HF1dXdy5c0c6/4hEIvj9fpxOJ06nk5WVFQRB4Pz58y2eeX88cl97LwrFiIXcE3+exMTERKtN+KM4DhcNhUKkUimSySTpdJqOjg5sNpuU59DZ2Uk0GmV9fZ333nuPYrGIyWSiUqmg0WiYnJxkc3NTyhw2Go3Mz89z/vx58vk87e3tLC0tsbq6ytjYGNvb25w6dYpPP/2UZrPJr3/9a6kGVCwWY2RkhM8++4xms8nPfvYzgsEg6XSavr4+MpmMdEO8d+8eBwcHHB0dSTf2aDQqVcHV6XSMjY2RSCR4/fXXWV1dBR52nqtUKpw6dYrp6WnW1tbw+/0YjUZ2d3fJ5/Nsbm5SKpWwWq3EYjG2t7el/h+pVIqf/exneL1eqQf23bt36e7uJpfLEYlEMBgM5PN5TCYTd+7coaOjg2vXrqHX6/nmm28IBAKkUilpZzQ7Oyu5BsfHx6UoMYfDwcLCApFIhJ6eHrq7u6V8kGKxyIcffsjNmzep1+tYrVYCgQDxeBytVktnZyepVIrV1VVqtRrd3d1UKhVyuVyLZ94fj9zX3ovipSwk+DycPn1aXF5ebrUZJ8aNGzdkXf1ycXGRtrY2UqmUJADHbppGoyE9BR/vHI6T244zqI/Lg2g0GjweD+vr69hsNqnsx/LystSzO5vNUiqVaGtrw+l0SofrJpNJ2o0AeDwe4vE47e3t0t8VCgUajQZ+v5+joyMMBgMOh0MqTWIymfD5fOzs7CCKIl1dXWSzWfL5PJlMRnJ/CYKAz+ejUqnwz//8z/zN3/wNwWCQcDhMrVbD4/EgiiKJRELaHRw3+XG73dLORBCE/6uyrdFoJBKJ0NbWRjAYZH9/X8rYdjgcUvb74+MVBAG3200sFkOn0+HxeMjlclL4r91uJxKJkEwmGRoaolAoUCwWqVarUmmT9fV1zp8/L/URaW9vx+fzsbKygslkks5w0uk0giDQ29sLwODgYMvm3YtA7mvvSTxtIUHFiIXa/EjeKHl8L3PzoxeBkq8dKH98sq06q6KioqLy8qEYsbBara024UR58803W23CiaLk8fX19bXahBNFydcOlD++p0UxYiH3mvlPYnFxsdUmnChKHp/S+yEo+dqB8sf3tChGLOr1eqtNOFESiUSrTThRlDw+pT/IKPnagfLH97QoRixUVFRUVE4OxYiF0vMsxsfHW23CiaLk8Sm9tpCSrx0of3xPi2LEotFotNqEE0UJyU1/CCWPTwn1kf4QSr52oPzxPS2KEQuld7Pa3NxstQknipLHl0qlWm3CiaLkawfKH9/TohixUFFRUVE5ORSTwS0IQg5Qbgo3eIB4q404QZQ8PiWPDdTxyZ1hURSfmKim/Z+w5H+IladJWZcrgiBMqeOTJ0oeG6jjkzuCIDxVww7VDaWioqKi8kRUsVBRUVFReSJKEoufttqAE0Ydn3xR8thAHZ/cearxKeaAW0VFRUXl5FDSzkJFRUVF5YRQUjSUispLhSAIIeATYAaYAH4qimK6tVa9OARBmADef/TjK8D/o6TxPY4gCP8iiuLftdqOF4kgCJ8ALmATQBTFL/7g55XohnrWL0GuKHQCK+YGJAjC56IofvDodQj4R6VcL0EQHMD/EkXxp49+/gT436IoXmytZS+eR3NyWhRFodW2vCgeXa+QKIo/eTQ3//1J105xYvE8X4IcUegEVswN6NHc+5djsXj0XkoURWcLzXphCILwPg/HN/DoZweQApxyFfffx6N5+P8q5doBCIKwcXztnhYlnln8kyiKPwEQRXFTjjeapyQEKGpRApPAPz728xfAxKMbkdyY4Hdcn0ciInse7db/4rG3Qo/eV9ScFAThE1EU/6PVdrxIHj1oph97/VQoSiye90uQG0qcwKC4G5Drd7yXBOQofL8TURRnHvvxL4GftMqWk+CRsCuximAISD7aMW0KgvAPj3aKfxBFiQXP+SXICQVPYED5NyAl8mjnNyGK4j8+8cPyYuK35qNScAHvi6L4H48exH4K/PuT/pHSoqGOv4TjQ8WfAluAYnyNPJzAittV/DaP3YA+eOKHX05+1y7ChfJchwD/xG/uCGXPo4dMRQbG8HBuSiIoimJaEASHIAghURR/74OoLMRCEIS/Bf7QYcznj1wYz/UltJqnHZ9cJ/AzXL/HkfsNaIbf4Yp6mefh8yAIwj/wMMorLQiCQ6Yuw9/H/xIEKX7E8Wgef6GAa/i7dktpHt4/fy+yEIvj6Jin4Lm+hFbzDOMDGU7gZxyfIm5AoihuPnadjt2H/1/rLHrxPHL3/sdj1+d9QBG73t9+eHkUpq6Ish+P5mbyeG0du7aftM6UGDr7OfAXj30JigydBRAEQVRS6CxIN6CZY/GT82H+byflHUfpKYFHY9v4rbc3nzUc82XnkTv0b3m40/0JD8OFX9oHs6fl0bj+Nw+v4QDwf/4UxeKZvwS5oeAJ/CdxA1JRkSOKEwsVFRUVlReP0kJnVVRUVFROAFUsVFRUVFSeiCoWKioqKipPRBULFRUVFZUnooqFioqKisoTUcVCRUVFReWJqGKhoqKiovJEVLFQUVFRUXkiqlioqKioqDwRWRQSVFGRI4/1gk8e17d6vO1vS41TUXlG1J2FisoJ8KhT4wwPKx7/5WN/9ZcouHmVinJRxUJF5WRIPyru+JfAvz32vix7kqioqIUEVVROEEEQUqIoOh+9DvGw0ZNaRVdFdqg7CxWVE+JRZ8PHXU7qrkJFtqhioaJysjwuFn8HfN4qQ1RU/hhUN5SKygkiCMK/8FAgXMC/AE6lNeNS+dNADZ1VUTkhBEEIiaL4d49eT/CwV7oqFCqyRHVDqaicAI8Os6cfe+ufeOiGUlGRJaobSkXlhBAE4W95mGcRAv5DCX3SVf50UcVCRUVFReWJqG4oFRUVFZUnooqFioqKisoTUcVCRUVFReWJqGKhoqKiovJEVLFQUVFRUXkiqlioqKioqDyR/x9rWbwnMpVZ2AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Area: 9.41 m^2\n",
      "Centroid axis location: 0.09 m\n",
      "Bending inertia: 35.70 m^4\n"
     ]
    }
   ],
   "source": [
    "from dolfin import *\n",
    "import ufl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "mesh = Mesh()\n",
    "with XDMFFile(\"box_girder_bridge.xdmf\") as in_file:\n",
    "    in_file.read(mesh)\n",
    "htop = max(mesh.coordinates()[:, 1])\n",
    "hbot = min(mesh.coordinates()[:, 1])\n",
    "wleft = min(mesh.coordinates()[:, 0])\n",
    "wright = max(mesh.coordinates()[:, 0])\n",
    "    \n",
    "    \n",
    "dx = Measure(\"dx\", domain=mesh)\n",
    "x = SpatialCoordinate(mesh)\n",
    "\n",
    "A = assemble(Constant(1.)*dx)\n",
    "y_G = assemble(x[0]*dx)/A\n",
    "z_G = assemble(x[1]*dx)/A\n",
    "\n",
    "z_0 = z_G\n",
    "I_y = assemble((x[1]-z_0)**2*dx)\n",
    "\n",
    "plt.figure()\n",
    "plot(mesh, linewidth=0.2)\n",
    "plt.xlabel(\"$y$\", fontsize=16)\n",
    "plt.ylabel(\"$z$\", fontsize=16)\n",
    "ax = plt.gca()\n",
    "ax.axhline(y=0, color='k', linewidth=0.5)\n",
    "ax.axvline(x=0, color='k', linewidth=0.5)\n",
    "plt.show()\n",
    "\n",
    "print(\"Area: {:.2f} m^2\".format(A))\n",
    "print(\"Centroid axis location: {:.2f} m\".format(z_G))\n",
    "print(\"Bending inertia: {:.2f} m^4\".format(I_y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considering that the cross-section is made of a concrete material with $E=50$ GPa at ambient temperature, we assume that the bottom surface is subject to an imposed temperature $T_{fire}$ while the top surface remains at the ambient temperature $T_{amb} = 20\\deg$. The remaining surfaces are assumed to have zero flux condition. We solve here a stationary heat transfer problem in order to determine the temperature field on the cross-section. This temperature field will decrease the value of the Young and shear moduli through a common stiffness reduction factor, the expression of which is assumed to be:\n",
    "\\begin{equation}\n",
    "k(T[\\deg]) \\approx \n",
    "\\exp\\left(-\\frac{T-20}{211}\\right), \\quad \\text{for } T \\geq 20^\\circ\\text{C}\n",
    "\\end{equation}\n",
    "that is:\n",
    "$$E(y,z)=k(T(y,z))E_{amb}, \\quad \\mu(y,z)=k(T(y,z))\\mu_{amb}$$ \n",
    "where $E_{amb}$, $\\mu_{amb}$ are the stiffness moduli at ambiant temperature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAADiCAYAAACftsmhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2de5AdV33nv2ded16SRyNbFnrg0QiNsKQgPL5GjKOVTJDBAUIqRBaELIEKWIY8IJWwUly7xVao3aWkIlULywISbCqEDRtbDgRYKINFShIyQvFojECSV5Y1kpFkJOEZyZr3487ZP7r7Tt++/b6nu885/ftU3Zq5/Tjn9+vvOb/+3dOnuxnnHARBEARBEETy1GVtAEEQBEEQRF6gxIsgCIIgCCIlKPEiCIIgCIJICUq8CIIgCIIgUoISL4IgCIIgiJSgxIsgCIIgCCIlGrI2oLHQwZtblyZaR0uBY2KKJVpHWpAv8qGLH0B6vozePPsK5/yOxCtKmHrWzhtZZ+L1tC2ow9jIXOL1JI0ufgDki6yk4csUv1RT/Mo88WpuXYret/xdonW8/5038I3vLUq0jrQgX+RDFz+A9Hw58q3ffCnxSlKgkXViZeMnE6/nE/95KT732NXE60kaXfwAyBdZScOXF6c/UVP8yjzxSoPD/W1ZmyAM8kU+dPED0MsXnfjuP9zIpN423hhp+zE247s+Kz+SIE++iG4HSaKCLrlIvNpa9BhCBcgXGdHFD0AvX3RiQUd94nVEPbmGLcN+Ek7Dj7TQ1Zc02kGSqKBLLibXF9dPZG2CMMgX+dDFD0AvX3TigXcvTKTcNt5Y/nhv0xDqE1QHkJwfWaCbL2m1AxGJnR8q6JKLES+CIAiikqCTbPTy5vcZY7Ou9dVBj5tQdMJPl1rbAeDdFrK8HJk1mSdea1e34/A37wcAbH3PTxKp4/nzhUTKDUuhrSXUdlNjwaMNSfliaZAm58+fx6MfXp1qnUm0sbCaiGwHSZFW+2KanH839C5Bf//HAQBrCp9PrJ4TR8aElueWdHmdZNtDjlCM2k6kVlnOk+7JwxOJnHTPTX1caHlhOH/+PP7sL9ONX6LbmNUOTh6ecCyvbgtx2oG9LGdbSCr5Et1XLOxtjLFP1FSW0MSLMdYLYJv59T4Aj3DOb4bdP6mT/+joKNrb213XJZXshT3Juu3jd+K9cKUpsJwskqg43HnnnanXGfXYhGkffpok1Q5qwe8Y+PUV3ak1fiV58vfTJerJOGzSFfZE69zeLwE799xkqLKySKSikkX8inpcwrYNuy7OthC3HQDVbUFE8hV0DFSIYcISL8ZYB4Ai53yv+X07gB8BuFdUHXHp7+/HAw884LouTpISdDKOc7J129/txPuOLSPl2/1VSbC88NNFFryOsb0N2DWxk2Q7sCOyHaigSRLIHL8Af12inozf2PSliu9BJ9pWHm6y8jgrVezvdtLd/heLsP+vXyknfz+b/lgk22VChb7i1zbs7WBel/m2UGs7sMoISr5EJ9kq6CJyxKsIYDeA/eb3gwAOMMY6ovxqVIHD37zfNfmq9UTrVp7bSVf1hEsHgi6Pi2wL1A5SITfxy0p23tj0Jd+ky3mibfE48U6YJ1pre3sC5rzs5LSByA57OwDgmXSJbgdW8pXnNiDsrkbO+UEAD9sWdZvLMw9aCxeKv8vh8DfvT/3Ed/ib9+O+e16Tap1JkoQuaeOmiegE3K3MpNqeDprEQeb4BSSji/PE53WybeH1nidb+3prG/u+9jJ/Nv0xvPmB1dqccHXpK5YuFrW2A7d9naNnSbYBFXRhnPNkCmZsDwBwznf7bVcsFnl/f38iNqTF1vf8JPBk29LeGljOxOh41bKpsQka2VAIv7YQ1Abc9LejU1tgjJ3gnBeztsOLPMUvAHh7w9+X/7dOmJUn0eDf6OPMeAacNfJhjXj8YPZDgqwk0sBqC852EKUNANXtANCnLdQavxK5q9GcL9HLOX/QY/1OADsBYNmyZTh06BAAoLu7GwsWLMDJkycBAIsXL8b69etx5MgRw9iGBmzevBkDAwO4desWAKBYLOLatWu4dOkSAGDNmjUoFAo4deoUAGDJkiUYGhpCqWSIXygU0NfXh/7+foyOjgIANm3ahMuXL+PKlSsAgLVr16K+vh5nzpwBACxduhSrVq3CsWPHAAAtLS3YtGkTjh8/jomJCfzNx4G+vnvx6b89hq7XGJeEfnr6NtTXcWzaYNRx4eoCDP5qId56j1HHrfFGHDq5HG+79xKamwzbvv9vr8X6lVexcokx0fH3fuc+TE1NlY/PypUrcfnyZVjJ8sKFC9Hb24ujR49idta4br5lyxacPn0aQ0NDAICNGzdiZGQEg4ODAICuri50dnZiYGAAALBo0SJs3LgRhw8fBuccjDFs3boVJ0+exI0bxhOAe3t7MTw8jIsXLwrV6YUXXgAzb29bsmQJenp6cPTo0cR0AoC+vj5cuHABV68ar5RYt24dSqUSzp49CwBYvnw5VqxYgePHjwMA2tvbUSwWcezYMUxNTQEANm/ejBdeeAHXr18HAHDO0dPTg3PnzuFvPg6sXHkn/vqzQ3jbmwwNXp1oxo9/sQy//aaX0FhvaPfdn96F+9Zex9JFhl0/OXMnOpbNYN1dNzBXmsPpC+24OtSEtxaHAQDXbzThP/7lWxPXiTGGFStW+PYnETrJjGzxq6enB0899RTa29sT6xffGfsDfOZD/4KVRWNk4tT/nkF9E7B+h/GogZefKeHGj0v4jb82biSZ+BXH6b+dwcZPNaJxIcMcgJ/8pyl0PVyPpRubMQfgdz+4heKXAvFrw4YNOHnyJOrrjQTri2f/He688048vu+HAICRX87h1Odncf9/aUJjs+HviV3TWP2hBnSsM5Kxs/tm0LCCYdU7jbuiz/1wBkNnGYp/bnx/5cUSdn/p4VR0OnbsGBobG8s6yRi/EhnxYoztA7A7zDB9Gr8YDx06lNpku7d94Lny/2FGuZxMjI7jh1+/x3N9mr4kjS6++Pnxtg88F7sdWPi1B9GkpYnMI16yxS8g3b7y/rpvlP9v5XVosz3jqcXjMSAT5mlkDBxf5e/zLFuXPg/kx5ePsH8CALSBBeoPGG0AmB/9+sbc+4XZGYY0dKk1fgl/cj1jbBfMoGX+cswVP/z6Pb4nykJrwfeT5kmWSJ4ffv0efPtLa301d8NK1qg9pEve4xdgnChbeV1F0tXC5pOulnqXj7neL+ki1OSr/H24nRlJl5v2QGX7sNpMK69LPelSBaEjXuYt2AOc80HrO+f8Sb990vjFODc3h7q69N+O9LsfM4Z/vU6ubjzxt12+67PyJQl08SWMHzv+6mKosqbGjcsB3/7S2lrNikVamsg44iVr/AKy6yufYI+XT7h2mhvc3+m5e+IPfMvTpc8D+fJlT8v/cV0+OVu5z0TJGP36HH+vUPuikIYu0ox4Mca6ARwAcJ4xxhljHMAeUeXXwunTpzOp99tfWiv8BJqVL0mgiy9BfoRNugAjSc8q6QL00SQqMscvIDtdPsffi5Z6I9Gyf1oLJddPEDq1rzz54qW3s1201GebdAFq6CJscr35K1HKF4FYEzWzpqWtueYyZPFFBLr4EuRHGN0nxsI90TtpdNEkKjLHLyBbXTrb3J/F1Vyofgff11/7FXzgl494lqVT+8qTL25aT041VCXbYZLvpFFBl8zf1ZgHRCRchN5QGyFkxe2kayzP70uO88SB130RzYVorw0i/MlF4rVx48asTfClpS38HDDZfYmCLr4E+WHXd2JsKmlzakIXTXQjS138EqzmpmjJl07tK0++OHWenG6UNvFWQRc9ZgYGMDIykmn9LW0F34+dj+654VtW1r6IRBdfovgR1Ba+9ulsn3Gliya6kaUuzU0znh8AaG6ervj4oVP7yosvbvoGtYksUUGXXCRe1sP3dIB8kQ9d/AD08kUnstTFeeL1SrSaCtNoKvgnXjq1r7z4YulqaRu2PWSFCrrk4lKjrDS30nXzPOCm8+R49r8MCSIMv3PqL/CDe/eG2rYg6eUnIj6FwgympowYFpRYE+HIReLV1dWVaf0iE6ysfRGJLr74+fEX/2PUdbmsSbcumuiGLLrUmljJ4ocI8uRLkO5WYvb2E7tEmRQbFXTJReLV2dmZtQnCIF/kQxc/AL180YmsdXE78cYZ/cjaD5HkxRc/naenjHd3yjTSqYIuuZjjZb1QNSv++5+3V3xvbW30/fiRtS8i0cUXPz9q0ToLdNFEN7LWxT7PxzmXq1CYrviceOgxz3Ky9kMkefHFqa8drzaRJSrokosRLxmQ8SRLZA+1C0IFnCdcPxolOQETtfPz3/0rNBaAGXNkC4jWFgh3cpF4LVq0KGsTQtHcbLwQ7dP/OItP/aG7NKr4EgZdfPHzw9J0cjL7JzqHQRdNdEMmXWpJrGTyo1by5EuQ5vbELGtU0CUXiZcsD1SzTsK1IIsvItDFlzB+hNFehuRMF010I2td/E68jRHm92Tth0jy4IubtjNT1aP0Mo1yqqBLLuZ4HT58OGsThCRdgBy+iEIXX0T58d8+HP4NBkmhiya6kbUujYUZz4/btl5k7YdI8uCLl75h20IWqKBLLka8OOdZm+BLa0v4/Fd2X6Kgiy9efnz2m3NV2o5PzKVhUmx00UQ3ZNKllhOsTH7USl58CdLbbQQsS1TQJReJF2MsaxMiJVd+yOCLKHTxJYofotpBUuiiiW5krcvdT3weL37gY77bhEnIsvZDJHnwxa6pV4Ily0iXhQq6sKyzw2KxyPv7+zO1IQ0++81oIx2ffI/cJ2giGNLcG8bYCc55MWs7aiUv8QsAXvrIR0Jv29g8jWVf+IcErSHSIIrmd331qwlaIhe1xq9cRPqTJ09mbYInzc3M9eOFzL5ERRdfvPzw0tZP36zRRRPdkFmXxubpqo8XMvsRFd19efnP/iiStrKggi65uNR448aNrE3AJ99Thy98v/bRRRl8EYUuvsTxQ9bkSxdNdEMGXUScdGXwQxR59UX25EsFXXKReKnIF77P8WfvkPPkTBAEAQANEj1GgBCLn7azEj23S0VykXj19vZmbUIVrTFHPGT0JS66+OLlR5DG45Py3X2jiya6IYMuIpIsGfwQRZ59kTnhVkGXXMzxGh4eztoEAMaJ2PrERRZfRKCLL25+/N2/Bu9nbw+1tAmR6KKJbsisS0PzjOvHDZn9iIruvkTRVSZU0CUXidfFixezNiESLQWOloL7aIhqvvihiy9uflgaeuno5I9/S7BRMdFFE92QQRcRJ2IZ/BBFXn2RPSFTQZdcXGqUjbAnY0IPgvSemJJjtIsg4tLg8iynW59+NxZ+6jsZWEOIwE3TWckelqoquRjx6u7uztoEAP4jWWGRxRcR6OJLrX7IlIjrooluyKBL53/9ZzQUZlw/YZHBD1Ho7MutT7/bdTsv/aO0gaRRQZdcjHgtWLAgaxNC0dzofOBm9fsdVfElDLr44uaHU8vJGTV+4+iiiW6opEtDq/fEa5X8CEJnX+wazo4H38Eo08imCrqocTaoEZkeqNbcOOf5cfKtfytVLZPJl1rRxZcwfvjp7qZ9VuiiiW7IpEtD67Tvxw+Z/KiVvPhSi95ZoIIuuUi8ZOEPfpPm8hAEoS4ynmiJbKE2EZ1cXGpcvHhx1iYE0tJUPbplUHm5UQVfwqKLL25+uOk5MV196Vg2dNFEN1TQpS7EnW0q+BEWnX3x0nJuUv7J9SrokovEa/369VmbUMY7wQqHTL7Uii6+OP146mfuvwD9tZcjKdNFE92QRZcwyZUfsvghgjz6Uqv+aaCCLrm41HjkyJGsTRAG+SIfuvgB6OWLTuiiiy5+AOSLrKjgSy5GvFSguaG2kTBCHvy0nJytHtl66I303jNCbeo9RkJmvroFjR+R/0RIzDPz1S2ob3ZfV1LgUqMK5CLxamiQx81aEyyZfKkVXXyJ4ofsCbYumuiGLLp4JVhhkcUPEeTRl1r1TwMVdGGcZ/vwxmKxyPv7+zO1IU0OnZqIvM8DG1oSsIRIiqga51FfxtgJznkxaztqJW/xCzBGRKJAI15qEVVfIH8a1xq/cjHHa2BgIGsTfGlumPb92JHdlyjo4ovTjyh6yoYumuiGKrrUNc9Ufeyo4kcYdPXFTUO/SfWyJV0q6CL/mJwAbt26lbUJZWo98crkS63o4ktUP6rbgDwjXrpoohsy6VLLnW0y+VErefNFhTsaATV0ycWIl0y8+fW3Rdq+UK9GYycMfvr/Xs3aBIIgCEJicjHiVSzKO5UkamIlsy9R0cUXux9h9JwqyXtnkC6a6IaqurC2ytFdVf1wQ1dfWNs0+Ji6d1qroEsuRryuXbuWtQkVFOpnyp+oyOZLLejiS1Q/7PrLNqKpiya6IbMurG3a8wMAc9+6t7ytzH5ERWdf/DR1JtOyoYIuuUi8Ll26lLUJsWmum6r4rrIvTnTxxe5Hc91UlWYqoYsmuiGTLrWchGXyo1Z09MWeJPshcxKmgi65uNQoM2FO0s+fv4K7Vy9PwRpCFEG6Ts4VAAD3rLk9DXMIgiAISchF4rVmzZqsTaiglhER2XypBV18ieOHrKNiumiiG7rooosfgJ6+sJYZ8Al556CGQQVdhCZejLFuANsBDADoBbCfc35TZB1xKBQKWZsQiwImq5cp6osbuvhi98PSbAoe79yQHF00iYOs8QuQS5e63ztRdUmKtYSbqyiTH7Wiqy9+WqqQlKmgi+g5Xvs453s55wcBPAlgj+DyY3Hq1KmsTXClgEnfjxuy+hIHXXyx/Bg8f768LI62MqCLJjGRMn4B8unCWmYqPmGRzY9ayKMvTt3rfu9EwpZFRwVdhCVe5q/FMpzzQQA7RJWvE7KffInkoTYgFxS/CIJIC5EjXr0AqoblnQEtC5YsWZK1CRV0r14darum0kT5YyGbL7Wgiy+WH256qYYumsRA2vgFqK8Lf+oNANT3ww75Iicq+CJyjleny7JhAB0C64hFT09P1ia4EucELasvcdDFFzc/grSdrpfnNUF2dNEkBtLGL0BhXVpmK74q64cLuvnCD6+bf3vZhLr33amgSyZHlzG2E8BOAFi2bBkOHToEAOju7saCBQtw8uRJAMDixYuxfv16HDlivISzoaEBmzdvxsDAQPl9TMViEdeuXSs/u2PNmjUoFArl67xLlizBtWvXwBgDYEy86+vrQ39/P0ZHRwEAmzZtwuXLl3HlyhUAwNq1a1FfX48zZ84AAJYuXYpVq1bh2LFjAICWlhZs2rQJx48fx8SEcYLt6+vDhQsXcPXqVQDAunXrUCqVcPbsWQDA8uXLsWLFChw/fhwAUMdL6F5+J8798mXMluYAAD13LcPVV27g1phR5oo7F2NmtoRrQzfx4suvYOXKlbh48SLq6+sBAAsXLkRvby+OHj2K2VkjwG3ZsgWnT5/G0NAQAGDjxo0YGRnB4OAgAKCrqwudnZ3lF4kuWrQIGzduxOHDh8E5B2MMW7duxcmTJ3Hjxg0AQG9vL4aHh3Hx4kWhOj3//PNoaGgo69TT04OjR49KpVN7ezuKxSKOHTuGqSnjTsTNmzfjhRdewPXr1wEAs7OzuPvuu/H8hcvG8bitHQvbW3HhirG+pdCErmVLcPalK5ib4wCA13ctx+XrQ/ileTxk0alUKqGrq8u3P4nQSWXSjl89PT146qmn0N7eLk2/aJvZgWLjEzhW+gCmeSsA4Dfbv4IXJt+CX8++DgCwoeV7mJxegBentgAAVp4/T/FLwvi1YcMGDAwMoGn6Tw2dmp7DkqazODH+PkOnumu4p+1JPDPyCGZhPNH+39V/BWdKb8eQ2fZl0umZZ54pT7CXNX4xznnNhQAAY2w7gEc55w/alt0AcK85X8KVYrHI+/v7hdjgxaFDh/DAAw8kWkdULr8QPAGwMDte8f2OdW+S0pe46OKL5cevz/xbedlUQ2uofVf0bEjKrFikpQlj7ATnXJp3e8gcvwD5+go/vC7yPmzrGen8qAXdfNnK/iTyfmzrmQSsqY00dKk1fokc8RqAy3C9X9BKC1lvL3UmVqH2kdSXOOjii5sfQdpONbRKl3QB+mgSA2njF6CPLrr4AWjoi3wPoY+FCroIm1zvDFDmpNQnRJVfC319fVmbUEWcpAuQ05e46OJLX19fxWhXGOLqnzS6aBIVmeMXoI8uuvgBkC+yooIvop/j9TBjbBdjbBuA7ZzzRwWXH4s0LgUkQcPMRMUHUNcXN3Txpb+/31UrFdFFk5hIGb8A+XSJe4lJNj9qgXyRExV8ETq53vzVuNf8elBk2bVgTZqTjTgnaFl9iYMuvoyOjgKLWMWyIG1nG+W8q1EXTeIga/wC1NSFt5Yqvz+7FqNj+zKyRjwqauLF6OgoeFulXmy83ncfGed3AWroou49o4qj8qgIUTukP6E6zsSKUJjx00Bb5aIgfZnvWsIP0ZcapWTTpk1Zm1DFoo1bQ23XMD1W/gBy+hIXXXzZtGlThU6WViqiiya6oYsuuvgB6OXLfZ1fy9oEYaigSy4Sr8uXL2dtgifOE7bfCXzk2e9L7UtUdPHFzY8ousqELprohi666OIHoJcvV8bfmLUJwlBBl1wkXtaD0GQjzglYVl/ioIsvcfyQNfnSRRPd0EUXXfwA9PLl5cmNWZsgDBV0oTleklM/7ZwLFO7BnES62HUqNck5cZ4gCILInlwkXmvXrs3aBE+qEyt/ZPYlKrr40jU3XPE9SNNSUwsW3PeOJE2KjS6a6IYuuujiB6CXLz0LpLqJtyZU0CUXiZf1bjDZiJp0AfL6EgddfKmPeHtPHN3TQhdNdENVXeba5iq+q+qHGzr5wpqnK7SqG1N3FpIKuqh7dCNgvdxSVeomx8of1X2xo4svL47ysj6qo4smuiGjLnX3ncVc25zvx8mZAWleBlAzMmoSh9kza3B26LcrlgXp6qatLKigSy5GvGRHhxM2YRBGy7nmtsBtCIIgZKVh3bmsTVCaXCReS5cuzdoEV+IkXLdPXEvAkmyQVZeo3NEwG2l7mRNtXTTRDV10WdIm/2hEWHTRBCBd0iYXlxpXrVqVtQmuNL/1j8JvPDkOTI7jtaVXkzMoZWTVJSqvLb1a1kd1dNFEN3TR5a7bfpK1CcLQRROAdEmbXCRex44dy9qEYKwTt9fH5ERpUYZGikUJXUJQoUlIHWVFF010Qxddnn35I1mbIAxdNAFIl7TJxaVGqVHgZEx4M/m9LwNYHGEH0psgCCLP5GLEq6VFkwdajo+heW46ayuEoYUulibj8s7bioIWmmiILro0N+gzVUIXTQDSJW1yMeKlwkszAYQ6ed8zezUFQ9JBGV0CKGsSpF+rcTdj8zs/mrBF8dFFE93QRZfisr/H7Jm/1+KuOF00AQxddEEFXXIx4nX8+PGsTfBmfGz+E4LnGuS/YyMsUusSgdCaRNA5K3TRRDd00eXZax9CqTXiE4clRQdNpi72oNTK0P/yh7I2RRgq6JKLxGtiQt4nhUdlYq4O41/7dNZmCEEHXebGxzAxp0830kETHZFVF+fIVamV+X4mZ27LyFLxyKpJHCbmbvPVTSVU0CUXlxpVYi5oRKQ9HTuIaPjpVtdKD0wl9EW1EzMRHbvGha4XMrRED3KRePX19WVtgieBiZaD3rHBhCxJH5l1iUKQJlE1zhJdNNENXXQpdv+vrE0Qhi6aAKRL2uhzjcSHCxcuZG2CJ60f/FSk7S813Y45BYZSwyCzLmGwLvlearo9Y0vEobomuqKLLr985f6sTRCGLpoApEva5CLxunpVrTsB5yYmPD+/rtfnWqNqujixa+Kll2qoromu6KLL9Vt3Z22CMHTRBCBd0iYXlxplR8UTNBEO0pYgqpm62ENzhYjckosRr3Xr1mVtgjBWj1zK2gRh6KILaUIkjS669LzmqaxNEIbqmkxd7Cn/T7qkSy5GvEqlUtYmxKI0Wv16mdm2VpTG9HjtjKq6WFj6eGlS395ataz9o59J3K5aUF0TXdFFl7k5fU45umgCkC5pk4sRr7Nnz2Ztgi+l0XHXjxsXO7sAAK9+9hMpWpgMsusSFksTJ2E1lQldNNENXXR58dq2rE0Qhi6aAKRL2uQi8SIIgiCSgeZqEUQ09Blf9GH58uVZmxCb2fHKydm3D1+pWqYqKusCzGvj1KShVf6XtHqhuia6oosuSzt+nrUJwtBFEyC8Liok2SrokovEa8WKFVmbEEjYZOqOmy8nbEl6qKCLF0Of3ln+36mJyomxyprojC66vKbjuaxNEIYumgCkS9rk4lKj7C/NvO2Tnwu97fNd9yZoSbrIrktYSBMiaXTR5bmLHyz/P359bYaW1I4umgCVuqiOCrrkYsRLB2ZGjInZvDRX/p/IDrsGliaNC6rvYiSIvGN/zx+vU++lyzoyfn0t4KFL/TjPyqzckIvEq71djae9h0moCpOjKViSDqroEoSlSZB+KiRmumiiG7LrEjaZam1+JWFL0kN2TaJg10X1xFgFXXJxqbFYLGZtQiBhR7HWnJ0fRr36H/59Uuakggq6hMGuiR+LP7U/YUtqRxdNdEMXXd7Q842sTRCGLpoApEva5CLxOnbsWNYmCGFmbALPv/5+zIxNlD8qo7Iudg2ef70+L5hVWROd0UWXgTMfztoEYeiiCUC6pE0uLjVOTU1lbUIk/BKq2cZCipYki2q6eDHbWPDUrLFNrUdL6KKJbuiiy/RsW9YmCEMXTQDSJW1ykXipgOqjV3ni0p/8fuhtSVciD7QuOav8XYoEkRa5SLw2b96ctQnCeO2Pvo3J2dmszRCCqrpMvlqZTK159ocZWSIeVTXRHV10Ka7/ctYmCEMXTYBwurQukf9VPIAauuRijtcLL8j/tF0vJl+dqPj8+u7erE0Shsq62Lm8/O4KjVRGF010QxddLlz+rYrvKo+SqaqJ2zF36qIyKuiSi8Tr+vXrWZsQyMov/nNVkuV2Eh+7c/6pvJOvTuAX73oLfvGut6RpqjBU0MXJL971lipd7JoA1cmySgmZiprkAV10GXq1p2qZqsmXipp4HWs3XVRFBV1ykXjphtuJXNXkSyXsx1ilZIogZEfV5Esl6BjLQy4Srw0bNmRtghAmX7U2afAAAB1BSURBVJ3AomcOua6bGplE/9a+dA2qEVV06d/ah/6tfZgamaxaZ2miSxKmiiZ5Qxdd1nZ913OdaomBSpoEHVs/XVRDBV1ykXipcHupG26Xqkqt1bf92hMClZIvFXRxHk+35MvSxO+y4pp//H4yBgpGBU3yiC66TE0v9F0/fn2tMgmYCpqEPZ5BuqiECroIS7wYY72MsV3m5wBjrENU2bVy7ty5rE0IRZg5QTffMP9C5qmRSddEwBqhkR3ZdfE6hs5jbtfEQsX5XYD8miQJxbDkufjy1lDbqZB8ya5JlGMYVhcVkF0XQFDiZQaoIud8L+d8L4DHAfxIRNmEO24Jl8X40CTGhyZxZMM9KVqkF0c23IPxIe9j7JX0EmpCMUw+VEi+ZIWOndyIGvEqAtht+34QQK8svxhXrlyZtQlCmBqZRMupnwcmXXZkTr5k1OXIhnsqjplf8gXMa6ILMmqSEhTDaiTMc56W3XEiUpkyJxCyahLnmEXVRWZk1cWOkMSLc34QwMO2Rd3m8psiyq+VO++8M2sTYmONrFjJVvOFQc9tvZIEWZMv2XTxOk5ByVfzhcEKnVQeCZNNk7SgGJYOizuiP4RT1nlfsmlSy3EK0kWVh6cC8unihrAn13POB2xf3wtgr9e2jLGdAHYCwLJly3Do0CEAQHd3NxYsWICTJ08CABYvXoz169fjyJEjhrENDdi8eTMGBgZw69YtAMabyK9du4ZLly4BANasWYNCoYBTp04BAJYsWYJr166BMQYAKBQK6OvrQ39/P0ZHRwEAmzZtwuXLl3HlyhUAwNq1a1FfX48zZ84AAJYuXYpVq1aVX77Z0tKCTZs24fjx45iYMObv9PX14cKFC7h69SoAYN26dSiVSjh71miwy5cvx4oVK3D8+HEAQHt7O4rFIo4dO1aeDLiooRGvrO/F1Gu7AAC3Hf5XlF7ThtHiJgBA6/OnMLZhI1ipBABoHHoFi37wPfx6x/tRYoaUC778VUw89CBmu+4y9vnO91C643b83y9+CQCw4R2/jc7OTgwMGHItWrQIGzduxOHDh8E5B2MMW7duxcmTJ3Hjxg0AQG9vL4aHh3Hx4kWhOj3//PNoaGgo69TT04OjR49motOJb/0L8Cc70fSL02j6+S8w+ofvAwDUDQ2j/fEnce1d28EWGpPobz/wDYxsur+s01xjExYe+3GFTrMXBjH8jneXdfoNAEePHsWs+daBLVu24PTp0xgaGgIAbNy4ESMjIxgcNBLrrq6uTHQqlUro6ury7U8idJKRsDEs7fjV09ODQ4cOob29Xer4tXnzZpx76aHyM6HWdn0XU9MLy/OHlt1xApevbUJdndEH2luvYsPrnsCzpz6K0lwTAOBNv/EFnHvpnbhxaxUA4O5V38LYxBL88upvAvif6H7dOyl+uej0/OlvAvg4li4+iaW3/ww/O/tBAEBr8yt4Q883MHDmw+X3MRbXfxkXLv9WhU6nXtyBxoaJsk6LO87iF+feX9bpTUvUiF/FYhHPPPMMCoVCWScZ4xfjnNdcSEWBxtD8Ac75g2G2LxaLvL+/X6gNTg4dOoQHHngg0TpEETQp/vr7P4gl3/ha+XvQaIydiV9NAwDePnQ6nnGCkUWXHyxeDwBoeU1TqO1bFzdXfHdq4kbx8LF4xqVMWpowxk5wzouJVxSDKDEsjfgFyNNXgggacfnpzz+ON7/h8zXVIcvoiyyaiBgNDNJFlmMehjR0qTV++Y54mb/sVvts8rQ5RG9nDyqH7DNn4UJ9bpVtHHoFQLSEC5hPuoD5RCNrpv/49/GD3//TrM0oM/Gr6VDJl3XsrQTM0kQHdOorAMUw2WhvvVpzGaIuO85WP5knEi2tO3Br7NHY+zeM1Va/SEToIgsq9BWhI16MsV0A9nPObzLGOsLMj0jrF6Mq+I14RU22LOxJFxFM2JEvO85RMCeqjHilhawjXlFjGMWvSmSci+VFrYlXrciUeAWh0ohXGtQav0Q+x2s7gCdtgWqbqLJrxbqeqxLWIyHsHwAY+ciHIpUjc9I1tfuRrE1wJeoxG/nIhzz1AtRKulTsK6KgGJY8z576aNYmCGOgXx9fdNJFhb4iZHI9Y6wbwAHzf2vxIIAnRZRfK9aEQBUIGtXiTeFGY2ROuMoUGrO2wJOwlx0Bb03ijlBmiUp9RSQUw9LBmkSfNSJGu+ZKtfky2ybPqJcsuohAhb4iJPHinA8CYIEbEomjRMKlCNaxjHPpkVALimEEQaSF8Lsao5LGHIm5uTnU1anxWsqgZ27xujqwubmq5SomXF6+yIxbEhbkx5ZTzyVpklDS6iuyzvGKSlpzvFSKYX7zvOZ4HepYtn1e1Nyuubk61NXV7osMo15Buqg0xyuNviLNHC+ZOX1ajscniGDiIeMO94lfTVd8VGR2x0NZmxAZt+NuaaIDOvUVndBFl3MvvTNrE4Rx/pwYX7Ke5A/466JS0gWo0VeEPUBVZqyHvKmGW0I1vXwlmKKJlpO5nq6sTaiZiV9Nu2qi6uVJVfuK7uiii/Vg1KwQmeS8ejNbX0SStS4iUaGv5CLxUoktp56T5jlbRHxUHYUkCB2RYVTJD8s+GS47EsmTi0uNGzduzNoEYTR+/TtZmyAMXXzRxQ9Ar76iE7rocveqb6VeZ1JJV8/rxfuSVYKYhS5JoUJfyUXiNTIykrUJwuCvuSNrE4Shiy+6+AHo1Vd0QhddxiaWpFLPbNv8JynGxpLxJQ3bnaSlSxqo0FdykXhZL+7Ugdlt/u9yVAldfPHzQ5b3YoZFp76iE7roYrzsOhnSTliuXErOF4u0fEpSl7RRoa/QHC+CIAhCGWSfr5UUXn7TvDD1yPw5XoyxXwN4KeFqbgegy5uMyRf50MUPID1f7uKcK3+NNqX4BejTxnTxAyBfZCUNX2qKX5knXmnAGOvX4WGNAPkiI7r4Aejli07ooosufgDki6yo4Esu5ngRBEEQBEHIACVeBEEQBEEQKZGXxGt/1gYIhHyRD138APTyRSd00UUXPwDyRVak9yUXc7wIgiAIgiBkIC8jXgRBEARBEJlDz/EiUoUx1g1gO4ABAL0A9nPOb2ZrVTwYY70Atplf7wPwiKq+2GGM7eOcP5q1HQQhI7rEMIpf2ZHbS40qiOOHqp2GMfY05/xB8/9uALtV1IEx1gFgB+d8v/l9O4DHOOf3ZmtZbZjt6gTnnGVtC+GPyjFM1fgF6BHDKH5lSy4vNZri7MzajriYnabIOd/LOd8L4HEAP8rYrEDMIFWGcz4IYEdG5tRKEcBu2/eDAHpNbVSmG4ASJ8A8o3IMUzV+AVrFMIpfGZLLxAuKiOODqp2mFy7H3RnMVIBzfhDAw7ZF3eZyZdsVY2w75/zJrO0gQqFyDFM1fgGaxDCKX9mSu8RLJXG8ULjTdLosGwagQsCtgnM+YPv6XgB7s7KlVswTh/xvlyWUj2EKxy9AoxhG8Ss7cjW5XjVx/NCp06iO+Uu915r3oSi9Kp/M84IuMYzilzxQ/EqfXCVeUEycMCjWadx+GXZC3UsmFntQ+QteKRhj22Bc7iHkR6sYplj8AvSMYRS/Ukb5xIsxthPAap9NnuacH1RBnLC+OJap1GkG4DJUb05QVRLG2C4YdzXdZIx1KHK5xI0djJVvBOow2+JBlbVRBV1iWA7iF6BZDKP4lQ25eZyEGbTsEyD3AXgUEosThNlp9qvUaRhjJ6xbllW9FdvCvAV7wGo/qs+9sWCMcdlvx84jusUwFeMXoE8Mo/iVHblJvJyoII4fqnYa58MHzdvJlcP047xj8SDn3O8Xv9SYl312whiF2Atgn4on9LygcgxTNX4BesQwil/ZkrvESyVxvNCx0xAEEQ7VYxjFLyLv5C7xIgiCIAiCyIrcPceLIAiCIAgiKyjxIgiCIAiCSAlKvAiCIAiCIFKCEi+CIAiCIIiUoMSLIAiCIAgiJSjxIgiCIAiCSAlKvAiCIAiCIFKCEi+CIAiCIIiUoMSLIAiCIAgiJSjxIgiCIAiCSAlKvAiCIAiCIFKCEi+CIAiCIIiUoMSLIAiCIAgiJSjxIgiCIAiCSImGrA3ovPPNfGb6ZtVyVueeE9bVe+eKdV771DHPfTzr8dnHa51HUb77MO9qPMvzMc1zXR3jnvswj3Wx9oHPPphzX87dlxvlRd+HzXms89vHa51XWT7ruOB9vNbxkth95vzKK7nrOjfrU96Md1s4V5r8Aef8Ic8NFGHbg218aKhUtZx5HBav5R7NHKy6aNs+Hp2deyz32t5jOffa3m8fj7q9yvKsw8sHv7K4e8D03j6irT42zUWtI2Ldcx6+xal7zmu5Vwj008LLJs/lHnFE0PYAUPLax+Oc5VVHyaeWWX6lpviVeeI1M30TvW/5u6rlhbYW1+1b2ls9yyq0Ftz3aWv23KelzX2f5tZGz31aPdY1N9d779Pi3nGam70bdavHupaC9wmtudG9sbQ0eUfw5gb3dc0N0577FOpn3Pepm/LeB5Ouy5tKE977zI67Lm+Y8d6nYXrMdXn9tPc+dZPu+2DSvX4AwLj7PnMeywFgbsLdhtKodz2z4+77zIx47zMz5r7P5Kvex8Bv3dSIu3bjQ+7LAWDiV97t56HhM7d7rlSIoaESDv/4rqrlDR5NoH7cve96La8b8z7hsnGPeDPhHtb5hHvc4mNNrsvnJr1jYMljndc+s+PudcxOeWzvU/fslHtZM5Meyz3qiL7cvXwAmPJYN+2xfMqjDq/tJz18A4DJafeyJj3qmJxybx/jU+7taXLWuw1OeJxWJjxOUWMeSc64xy+SCY9fHuM+v0hGmfu5aYzNeix33/4W8z6XXZ/aVVP8okuNBEEQBEEQKUGJF0EQBEEQREpQ4kUQBEEQBJESlHgRBEEQBEGkBCVeBEEQBEEQKUGJF0EQBEEQREpQ4kUQBEEQBJESlHgRBEEQBEGkBCVeBEEQBEEQKUGJF0EQBEEQREowzr1fP5OKAYydAjzeJSOW2wG8QvVQPVSPFPW8osO7GlOMX26kpT3VTXVT3ZXUFL8yf1cjgEnOeTHpShhj/VQP1UP1yF2PgqQSv9zIUhOqm+qmuuNDlxoJgiAIgiBSghIvgiAIgiCIlJAh8dpP9VA9VA/VoyhZHheqm+qmuhWsO/PJ9QRBEARBEHlBhhEvgiAIgiCIXECJF0EQBEEQREqk9jgJxlgvgG3m1/sAPMI5v2mu6wawHcAAgF4A+611Ecr+Cuf8Xsdyz3Lj1BlQzzYAwwC6BdXjeqzM9dsBdAIYBADO+cE4dTHGrDo6zHoe55wPuGy3j3P+qMPfWHq5lCWsXTDG9gA4D+AJADsADHPOnwwqK4F6RLYFV6391sXw5wAcbcxjO2HtQDai+BK3LXmti6GXVxzy60uubTai35HL9/IbRrsNW69nnPJru4J8FhpTYmjt1ceT1jrycRXld0AMjaO3a3/xqFdYv/arC5zzxD8wOsxO2/ftAE7Yvj9t+78bwL4IZW8zneUu6zzLjVpnQD27HN/31FBP0LHabtVnlhf7OAK4AaDDrR7bNlU+x9XLWZbodgFgDwBu+uXURGRb8KtHZFvw01pkO+AuH6cfwtqBjJ8ovsRtS17rItbtGodC9CXXNhu27rjle9UR0WfPOOXXdmv1OY5fgrV27eNJax33uAr02y+GRrILPudtl3qF9mvfuoI2EPExnT9v+95hHqAO09CnHdvfiFEHd3z3LLeWOt0EdClrX9x6/I6V+f28x35x6uq2/b/Tub+5fLu9nBqPnbMsoe0CwPaox0ZkPQm0BVetRbYDc/ttjmU7XbYT1g5k+0TxJW5b8loX9zg641CIuFHVZiP6Hbl8nzpejdpG7W3T2jeo7dbqcwy/hGrt08eT1jrycRXst18MjWyXuY57rYurZ5z2ZH1SmePFjeHRh22Lus3lN2Fko27Dht01VutXrvA6GWNPM8Y6zGHxAyFscMXvWJlDptbwZq9j1zh1Ddq+PgzjF4N93+3cHEKupR6vspJqFxGPjch6rOU1twU/rQW3g2FeffnyCce+wtqBpETxJW5bcl0H4KEIdXsS0Jfs5drbS2i/Y5bvVUc9gJkw9Zp1eMWpwLbrYlPc+FVzTEEErf36eNJaI95x9awDMdq4RwyNa1cYhPbroPaUSuIFALxy7tB7Aew1/+902XwYRhZfC37liq7zYbPMCwB6bY0jVj0+x6obwLDZ4AYZY7ts8x9i1cUY62aM7QJwwNGou2HOK3AQuR6fskS3i27zeAwyxvbYOqDotuBVDyCuLfhpLawd2IM1Y6wDQKdjmbB2IDFRfInblrzWLYtQty8+fQlwb7NR20rU8uFRxziAxrD1Au5xKqjtetgUuR8KjClRtPbr44lqHfO4wqeOqG3cNYbWYFcYRPdr3/6b+rsazQPWyzl/MO26E6QIYDfM67uMMXDO9wbsE4jLseqEMdT6oLl+P4zGuShuHZzzQbOcPY6RjV6XUY64BJYlol3YjzljbB+ApwGsjltezHpEtQU/rYW3A5PHAHzGsUxkOyBSwK0vebTZPS67iyxfWN/ziVMWVW1XhM9pxRQXQvXxpLVG+OMq8piEiaFZ2CWM1Ea8bOxB5TCpW3bYCfchyyj4lSusTmu4kXN+kHO+H4bQj5kdotZ63I5V+ZeOme13mDbErsss5wCAA7bh3YMem0eqJ6AsOzW3C/OYAyhfnrCGe4W2Ba96BLeFIK2FtwMYgd7+q1JYO5CcKL7EbUte616OUHdYnH3Jq83G1TBs+fCooxVAS9R6nXHKsbqi7frYFDV+iYwpUbT26+N2ktY67HG1bK7J74AYGteuMIju177HNdXEyxwq3m3OV7KMHYDLcJ3jun4c/MoVWWfFycksw7pdOnY9PsfKiSV+pLoYY9sYYydsi/rNv1YZOxhjOxljO2F0+J1mp4jjk1dZQb5G8gfAjzzqF9YWAuoR2RaCtPZbF7nNmX4Nu6wS2Q5kJYovcduS6zoAT0WoOxC3vuTTZiNrGLF8rzpKqE68XOsNEadc264In0XHFETT2q+PW/YlrXWU4+pZB6L57RdD49oVBqH9Oqj/ppZ4mdeqn7RlqduAagPNoF41YS4qfuUKrvMg5p+nUlF/3HoCjtWwrZN1AxjknN+MUdcwgMdt34tmWYPWrw3rY9a9P45PfmWF8DWKP/2wDT1b5QaVJbIeCGwLIbQW1Q4sqiaJimwHMhPkCzPmF3UEbRtnXZS6g/DqS/Bos1Hrjlq+j9//FKFezzhlW+Y2wVmEz0JjSox25trHHf4korVJ6OMqym/4xNC4dnmRZL/2qxdAOu9qNI0571g8yDlfbVtffgAZjzAnxsx0e2EMue6FcWun64MkeeU14Eh1BtSzHcaw5iCM7PegLbGIWk/QseqAcX37PIxh2M/YOmPUuqyH8wHAvTCelzJoW98B4/Zty+d93JhrEVkvt7LMVcLaBZt/qOBNAKs557tt60S2Bb96RLYFP62FtQNzn50wbkev2lZkO5CVgPZxAEZ/3x9i28jrItbtGodCxA3XNhu27rjle9UR0eegOOXadmv1OY5fQesi1u3ax5PWOu5xFei3ZwyNalfAeTvRfu0HvSSbIAiCIAgiJbKYXE8QBEEQBJFLKPEiCIIgCIJICUq8CIIgCIIgUoISL4IgCIIgiJSgxIvIPYyxXhb9FRNu5Ww3ywl1+z9BEISOBMXUvMfK1F8ZRETDvHX2Qc75w4yxp2Hcyv+kuW4bjNtkH4dxC+2jMJ57cwDG7cfbAdzLHU/4JaooovJxDNbrNcq3Hlsw46GFwPzDDAcBFM3nWz3JGAPUfYI7kQMct90DxjOQurnPq6GY8RDT3bzyfa47YfSBmwDsr/p6nFe+S1BqTD/2AFiVRKx0O3YCy94D49ESro+CyZCKmOok77GSEi/5eQLzrySw3l9l0QHgYdtzoh4E8KztuSRPm9srEwQB4/k1WSWL1jNyzGcX7YHtKcoeJ59eGO8Te4ISXEIRHuOcl18zYzt5w7bM2QcfsSdT5rOWVttizT7M//BTqh+Yff1RUeUFHTuB9Ww3/30CLg8djVBOZvE2r9ClRskxO8SA+f8Aqp8u7vlqAjNBcHuFg7SYiU/sIFIr5vFcbZ6MnrbZtQdAv/NXq6kJvUSaUAKPSz+fcWxT1QddEof7YOsfMOJSP+d8wC8m6U7IYyeKbhg/tG/6jVb6kXW8zSuUeEkOm3/BMhhj2xwn/jBD1/3Bm0jF7uBNkoVzvtv82I/vLsw/bd/J4x7LCUI2BgFsM6cpACj/uLO/5iTzPqgwqh071ezVArrUKDmOES/naEuY4eGi+atmGMB95msUegF8BUYiYb3FfbX5vQPm/A9zCH6bufxJGL9wO6xy7JWY2znrsfa15kxZc9Wsl5wWYZubYCvjQcZYJ4yTQbdlq2nPdvP7WznnAwF1VNgT4li5wuZf6O36Sz7ur02CSBvzlTO7ATzNGBuE8eNtnzUqE7IPWq9h6bD1jQ4AOxljAzD6nDO+3Ms5L1/O84gXHQB2mPt0mOv3O5d5vNYqUhywzXMbMMu1v5OyHB/dYo7L/oDxA7cYdOxs+xZtx+ag+Qout7hccdwcNj4I4CZjrNOmS1VcjWKv2T687HM9xk7biBBwzumjyQfGpPpdtu/dMCaIW993Athp/r8HRjCw1p0AsN32/Ybt/33Wfub37QAOhKxnn1UPjPdYAfOTQa312xx1bXf4tcdR/9NWWW51+Nnjcdysd395re8GwP22cRyb7qzbAn3o4/eBmSiZMYO79CdnH9zl6IN7HPHihL1/eMSXbvN/1/5p1mG3Y7vbMh+fQsUB0/cTjn3LMclmv2vMce5v3zbo2Dltso5NmOPm4q9TE9e4GsXeEPZVxXMP23xjqk3fXMZKGvHSm+0wfhHZ53Xca/4dQuWEWusOPTesO/gAlO9IORCynptWPXx+rsO9fP5lr52w/dr0YChgfUUd5p2HXvZEhhu/9gCPGxUYY908x/NaCHWwJlJzY7R8P4D95l19jwEQOXrhFl+sfu4VL/ZgfiTuaT7/Mu2KZT51RokDzikYzqsHfjFnh31/Hm00fTsq58YBwLBtGonfcQvCK65GsTfIPrd47ot5rrBsGeQuI3h5gxIvvVkMo6FbHUTIJE+XZ68E1TPs+P4YY2wIxuVLz4QlYkJjryMJv/fCuGvLLWj0wscPgpCIbsaY86T5BDySroR+VLj2TzMpXG0mSu81T9iPOJdx/8tbgXHATDQTx+PYLU6wylBx1QszYQtjnzOe+xKgVy6hyfV68ziMpKCMfVJtROx3R26D8Ws5cj3W/BDO+V4zKHUA6LRtb//l6Xx0hkUxwFaRfgMo/0osOssxk1C6FZtQiT2O792ovkPRvs6NWh586dU/HwOMkRTbqIzbslrreQLVMcTNH6+Yc9C5f4j4ZbHPaROATl7jM74C4moUexOxj6iERrw0wDYZshfGL9qbMCZEDphPB96F+UmkB61fjzCGkPthdLgigEfNibc7YEyc3WUb2l9tdlZrcn155Cegnm2WTWbn7Ycx/L8dRqc/AOPXtjUUvg/AbvPSnvVLdT+APbZg0Q/j190jpu0VdXjZU+tx5pzfyxjbZfplXX69SUGJUIx9bP7hp52ofvhmRR90xIuDMPpTEcaJfQDz86keM5/n1QGf+OLTP3ttz6a6adrhtqwKt1jjVY95OW63eQz6Mf+jco9pn3UZ1jXmmFMPdnvEF99jZ+67x6x70Dx2D9t88Dtu5STJse0gfOJqFHtD2OeM50QMmDnJjSA8YcYzrJ7lmt69ZwYZIQ9ANQNfrp9lRBBEvgkTU/McK+lSIxGGXL5PiyAIgiBEQ5caCV8cw8u5/HVCEARBEKKgES/CF3NS62rO+YMaJ13W/I2aJuCbw+v3IeJdPwRBEJrhG1PzHitpjhdBEARBEERK0IgXQRAEQRBESlDiRRAEQRAEkRKUeBEEQRAEQaQEJV4EQRAEQRApQYkXQRAEQRBESlDiRRAEQRAEkRL/HzyRbwuHubPhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "V = FunctionSpace(mesh, \"CG\", 1)\n",
    "v = TestFunction(V)\n",
    "u = TrialFunction(V)\n",
    "T = Function(V, name=\"Temperature\")\n",
    "T_amb = Constant(20.)\n",
    "T_fire = Constant(800.)\n",
    "\n",
    "a = dot(grad(u), grad(v))*dx\n",
    "L = Constant(0.)*v*dx\n",
    "bc = [DirichletBC(V, T_amb, \"near(x[1],{}) && on_boundary\".format(htop)), \n",
    "      DirichletBC(V, T_fire, \"near(x[1],{}) && on_boundary\".format(hbot))]\n",
    "solve(a == L, T, bc)\n",
    "\n",
    "k = exp(-(T-20)/211)\n",
    "\n",
    "plt.figure(figsize=(10,4))\n",
    "plt.subplot(1, 2, 1)\n",
    "pl = plot(T, cmap=\"coolwarm\")\n",
    "cbar = plt.colorbar(pl, orientation=\"horizontal\")\n",
    "cbar.set_label('\"Temperature [°C]')\n",
    "plt.subplot(1, 2, 2)\n",
    "pl = plot(k, cmap=\"plasma_r\")\n",
    "cbar = plt.colorbar(pl, orientation=\"horizontal\")\n",
    "cbar.set_label(\"Stiffness reduction factor [-]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the bottom of the cross-section is strongly affected by the stiffness reduction. We now compute the various moduli at ambient temperature and the corresponding values for this temperature field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Axial stiffness:\n",
      " H_N_amb = 470.42 GN\n",
      " H_N = 195.80 GN\n",
      " relative = -58.38%\n",
      "\n",
      "Bending stiffness:\n",
      " H_M_amb = 1784.97 GN.m^2\n",
      " H_M = 836.46 GN.m^2\n",
      " relative = -53.14%\n",
      "\n",
      "Coupling stiffness:\n",
      " H_NM_amb = 0.00 GN.m\n",
      " H_NM = 362.28 GN.m\n"
     ]
    }
   ],
   "source": [
    "E = Constant(50.)\n",
    "\n",
    "print(\"Axial stiffness:\")\n",
    "H_N_amb = float(E)*A\n",
    "H_N = assemble(k*E*dx)\n",
    "print(\" H_N_amb = {:.2f} GN\".format(H_N_amb))\n",
    "print(\" H_N = {:.2f} GN\".format(H_N))\n",
    "print(\" relative = {:+.2f}%\".format(100*(H_N/H_N_amb-1)))\n",
    "\n",
    "print(\"\\nBending stiffness:\")\n",
    "H_M_amb = float(E)*I_y\n",
    "H_M = assemble(k*E*(x[1]-z_0)**2*dx)\n",
    "print(\" H_M_amb = {:.2f} GN.m^2\".format(H_M_amb))\n",
    "print(\" H_M = {:.2f} GN.m^2\".format(H_M))\n",
    "print(\" relative = {:+.2f}%\".format(100*(H_M/H_M_amb-1)))\n",
    "\n",
    "print(\"\\nCoupling stiffness:\")\n",
    "H_NM_amb = float(E)*(z_G-z_0)\n",
    "H_NM = assemble(k*E*(x[1]-z_0)*dx)\n",
    "print(\" H_NM_amb = {:.2f} GN.m\".format(H_NM_amb))\n",
    "print(\" H_NM = {:.2f} GN.m\".format(H_NM))\n",
    "#print(\" relative = {:+.2f}%\".format(100*(H_NM/H_NM_amb-1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing shear correction coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Timoshenko cross-section kinematics assumes a constant distribution of transverse shear strains. The associated shear stress are also uniform for a homogeneous cross-section and therefore do not satisfy traction-free boundary conditions on the border. For a rectangular cross-section, the shear stress field induced by a uniform shear force is parabolic across the height of the section. The corresponding shear stiffness is therefore not equal to $\\mu A$ as would be the case for a uniform stress field but $\\kappa_z \\mu A$ where $\\kappa_z$ accounts for the non-constant distribution of the exact shear stress field.\n",
    "\n",
    "### Mindlin's energetic method\n",
    "\n",
    "Mindlin proposed to compute such a shear correction coefficient by equating, for a given value of the applied shear force, the exact elastic shear energy $U_1$ with that obtained using the simplified constant shear strain kinematics, noted $U_2$, corrected by $\\kappa_z$. More precisely:\n",
    "\\begin{align}\n",
    "U_1 &= \\dfrac{1}{2}\\int_S \\dfrac{\\sigma_{xz}^2 + \\sigma_{yz}^2}{\\mu(y,z)} \\dS \\\\\n",
    "U_2 &= \\dfrac{1}{2} V_z \\gamma_z = \\dfrac{1}{2} \\dfrac{V_z^2}{\\kappa_z \\overline{H}_V}\n",
    "\\end{align}\n",
    "where $\\sigma_{xz}(y,z)$ and $\\sigma_{yz}(y,z)$ are the shear stress fields on the cross-section corresponding to an imposed vertical shear force $V_z$, that is:\n",
    "\\begin{equation}\n",
    "V_z = \\int_S \\sigma_{xz}(y,z)\\dS\n",
    "\\end{equation}\n",
    "\n",
    "Equating $U_1=U_2$ as proposed by Mindlin's method therefore yields:\n",
    "\\begin{equation}\n",
    "\\kappa_z = \\dfrac{1}{2} \\dfrac{V_z^2}{U_1 \\overline{H}_V}\n",
    "\\end{equation}\n",
    "\n",
    "### Corresponding variational formulation\n",
    "\n",
    "It can be shown (see [[BAT90]](#References)) that the shear stress field $\\bsig(y,z) = \\sigma_{xy}(y,z)\\be_y + \\sigma_{xz}(y,z)\\be_z$ is the solution to the following linear system of equations defined on the cross-section:\n",
    "\\begin{equation}\n",
    "\\begin{cases}\n",
    "\\operatorname{div} \\bsig + f(y,z) = 0 & \\text{on } S \\\\\n",
    "\\bsig(y,z) = \\mu(y,z) \\boldsymbol{\\nabla}u(y,z) & \\\\\n",
    "\\bsig\\cdot\\bn = 0 & \\text{on } \\partial S \n",
    "\\end{cases}\n",
    "\\end{equation}\n",
    "where $f(y,z)$ is a given source term and $u(y,z)$ is an unknown function which can be related to the warping deformation of the cross-section induced by the imposed shear force.\n",
    "\n",
    "> Note that since the above problem is subjected to pure Neumann boundary conditions, the solution $u$ is defined up to an arbitrary constant. We must therefore enforce an additional equation (fixed value at some point, zero average, etc.) for removing this indeterminacy.\n",
    "\n",
    "This antiplane shear elasticity problem is akin to the stationnary heat equation so that the corresponding variational formulation is straightforward: Find $u\\in V$ such that:\n",
    "\\begin{equation}\n",
    "\\int_S \\mu \\nabla u \\cdot \\nabla v \\dS = \\int_S f u \\dS \\quad \\forall v\\in V\n",
    "\\end{equation}\n",
    "where the source term associated with a given vertical shear force $V_z$ is given by:\n",
    "\\begin{equation}\n",
    "f(y,z) = \\dfrac{E(y,z)}{\\Delta_H}((z-z_0)H_N - H_{MN})V_z\n",
    "\\end{equation}\n",
    "where $\\Delta_H = H_NH_M-H_{NM}^2$.\n",
    "\n",
    "> Note that for the above problem with homogeneous Neumann conditions to be well posed, the source term must be of zero average on $S$. This property is ensured by construction owing to the definitions of $H_N$ and $H_{NM}$.\n",
    "\n",
    "### Rectangular homogeneous cross section\n",
    "\n",
    "For a rectangular homogeneous cross section $S=[-b/2;b/2\\times [-h/2;h/2]$ of width $b$ and height $h$, we have:\n",
    "\\begin{equation}\n",
    "f(y,z) = \\dfrac{z}{I_y}V_z\n",
    "\\end{equation}\n",
    "The solution of the above problem is therefore trivial and given by:\n",
    "\\begin{equation}\n",
    "\\bsig = \\dfrac{1}{2 I_y}V_z\\left(\\frac{h^2}{4}-z^2\\right)\\be_z\n",
    "\\end{equation}\n",
    "and we have:\n",
    "\\begin{align}\n",
    "U_1 &= \\dfrac{1}{2}\\dfrac{1}{4 I_y^2 \\mu}V_z^2 \\int_S \\left(\\frac{h^2}{4}-z^2\\right)^2 dz = \\dfrac{6}{5}\\dfrac{1}{\\mu bh}V_z^2 \\\\\n",
    "\\overline{H}_V &= \\mu b h\n",
    "\\end{align}\n",
    "which yield the widely known shear correction factor for a rectangular cross-section:\n",
    "\\begin{equation}\n",
    "\\kappa_z = \\dfrac{5}{6}\n",
    "\\end{equation}\n",
    "\n",
    "### FEniCS implementation\n",
    "We now define and solve the corresponding variational problem for the reduced shear modulus $\\mu(y,z)=k(y,z)\\mu_{amb}$. We use a mixed function space involving a scalar Lagrange multiplier to enforce a zero-average condition on $u$ as [discussed here](https://fenicsproject.org/olddocs/dolfin/latest/python/demos/neumann-poisson/demo_neumann-poisson.py.html). We also perform a sanity check that the source term is indeed of zero average."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average f: -9.338840589517491e-11\n"
     ]
    }
   ],
   "source": [
    "nu = Constant(0.2)\n",
    "mu_amb = E/2/(1+nu)\n",
    "mu = k*mu_amb\n",
    "\n",
    "Ve = FiniteElement(\"CG\", mesh.ufl_cell(), 1)\n",
    "Re = FiniteElement(\"R\", mesh.ufl_cell(), 0)\n",
    "W = FunctionSpace(mesh, MixedElement([Ve, Re]))\n",
    "\n",
    "v, lamb_ = TestFunctions(W)\n",
    "u, lamb = TrialFunctions(W)\n",
    "\n",
    "V_z = Constant(1.)\n",
    "Delta_H = H_N*H_M-H_NM**2\n",
    "a_shear = dot(mu*grad(u), grad(v))*dx \n",
    "a_tot = a_shear + (u*lamb_+v*lamb)*dx\n",
    "f = (H_N*(x[1]-z_0)-H_NM)/Delta_H*k*E*V_z\n",
    "L_shear = f*v*dx\n",
    "\n",
    "print(\"Average f:\", assemble(f*dx))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now solve the problem and extract the displacement function. The stress components are then obtained and fields are saved to a `.xdmf` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = Function(W)\n",
    "solve(a_tot == L_shear, w)\n",
    "\n",
    "u = w.sub(0, True)\n",
    "u.rename(\"Displacement\", \"\")\n",
    "\n",
    "sig_xy = Function(V, name=\"Shear stress xy\")\n",
    "sig_xy.assign(project(mu*u.dx(0), V))\n",
    "sig_xz = Function(V, name=\"Shear stress xz\")\n",
    "sig_xz.assign(project(mu*u.dx(1), V))\n",
    "\n",
    "with XDMFFile(\"results.xdmf\") as out_file:\n",
    "    out_file.parameters[\"functions_share_mesh\"]=True\n",
    "    out_file.write(u, 0)\n",
    "    out_file.write(sig_xy, 0)\n",
    "    out_file.write(sig_xz, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the shear reduction coefficient and the corresponding shear stiffness are computed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kappa_z = 0.198\n",
      " shear stiffness relative reduction: -91.74%\n"
     ]
    }
   ],
   "source": [
    "energy = 0.5*assemble(ufl.energy_norm(a_shear, w))\n",
    "H_V = assemble(mu*dx)\n",
    "kappa = float(V_z)**2/2/H_V/energy\n",
    "print(\"kappa_z = {:.3f}\".format(kappa))\n",
    "print(\" shear stiffness relative reduction: {:+.2f}%\".format((kappa*H_V/float(mu_amb)/A-1)*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The corresponding shear stress distribution looks as follows:\n",
    "\n",
    "![](shear_stress.png)\n",
    "\n",
    "corresponding to the following warping function:\n",
    "\n",
    "![](warping.png)\n",
    "\n",
    "## References\n",
    "\n",
    "[BAT90] Batoz, J. L., & Dhatt, G. (1990). Modélisation des structures par éléments finis - vol 2 - poutres et plaques. HERMES."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Format de la Cellule Texte Brut",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}