JuliaFEM.jl/notebooks/2015-09-24-Ideal plastic Von Mises material.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ideal plastic Von Mises material\n",
    "\n",
    "Author(s): Olli Väinölä <olli.vainola@student.oulu.fi>\n",
    "\n",
    "In this notebook is an small tutorial, how to create a Von Mises material without any hardening. Equations are formulated into rate depended form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# imports\n",
    "using PyPlot\n",
    "using ForwardDiff\n",
    "using NLsolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a isotropic Hooke material."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6x6 Array{Float64,2}:\n",
       " 2.69231e5  1.15385e5  1.15385e5  0.0        0.0        0.0      \n",
       " 1.15385e5  2.69231e5  1.15385e5  0.0        0.0        0.0      \n",
       " 1.15385e5  1.15385e5  2.69231e5  0.0        0.0        0.0      \n",
       " 0.0        0.0        0.0        1.53846e5  0.0        0.0      \n",
       " 0.0        0.0        0.0        0.0        1.53846e5  0.0      \n",
       " 0.0        0.0        0.0        0.0        0.0        1.53846e5"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "\n",
    "\"\"\"\n",
    "Create a isotropic Hooke material matrix C \n",
    "\n",
    "More information: # http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_isotropic.cfm\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "    E: Float\n",
    "        Elastic modulus\n",
    "    ν: Float\n",
    "        Poisson constant\n",
    "\n",
    "Returns\n",
    "-------\n",
    "    Array{Float64, (6,6)}\n",
    "\"\"\"\n",
    "function hookeStiffnessTensor(E, ν)\n",
    "    a = 1 - ν\n",
    "    b = 1 - 2*ν\n",
    "    c = 1 + ν\n",
    "    multiplier = E / (b * c)\n",
    "    return Float64[a ν ν 0 0 0;\n",
    "                   ν a ν 0 0 0;\n",
    "                   ν ν a 0 0 0;\n",
    "                   0 0 0 b 0 0;\n",
    "                   0 0 0 0 b 0;\n",
    "                   0 0 0 0 0 b].*multiplier\n",
    "end\n",
    "\n",
    "# Pick material values\n",
    "E = 200.0e3\n",
    "ν =  0.3\n",
    "C = hookeStiffnessTensor(E, ν)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# Defining equations for the calculation\n",
    "\n",
    "Functions are defined for strain controller simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "calculate_stress (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "\"\"\"\n",
    "Equivalent tensile stress. \n",
    "\n",
    "More info can be found from: https://en.wikipedia.org/wiki/Von_Mises_yield_criterion\n",
    "    Section: Reduced von Mises equation for different stress conditions\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "    σ: Array{Float64, 6}\n",
    "        Stress in Voigt notation\n",
    "\n",
    "Returns\n",
    "-------\n",
    "    Float\n",
    "\"\"\"\n",
    "function σₑ(σ)\n",
    "    e1 = (σ[1] - σ[2])^2\n",
    "    e2 = (σ[2] - σ[3])^2\n",
    "    e3 = (σ[3] - σ[1])^2\n",
    "    e4 = σ[4]^2 \n",
    "    e5 = σ[5]^2\n",
    "    e6 = σ[6]^2\n",
    "    return sqrt((e1 + e2 + e3 + 6 * (e4 + e5 + e6)) / 2.)\n",
    "end\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Von Mises Yield criterion\n",
    "\n",
    "More info can be found from: http://csm.mech.utah.edu/content/wp-content/uploads/2011/10/9tutorialOnJ2Plasticity.pdf\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "    σ: Array{Float64, 6}\n",
    "        Stress in Voigt notation\n",
    "    k: Float64\n",
    "        Material constant, Yield limit\n",
    "\n",
    "Returns\n",
    "-------\n",
    "    Float\n",
    "\"\"\"\n",
    "function vonMisesYield(σ, k)\n",
    "    σₑ(σ) - k\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "Function for NLsolve. Inside this function are the equations which we want to find root.\n",
    "Ψ is the yield function below. Functions defined here:\n",
    "\n",
    "    dσ - C (dϵ - dλ*dΨ/dσ) = 0\n",
    "                 σₑ(σ) - k = 0\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "    params: Array{Float64, 7}\n",
    "        Array containing values from solver\n",
    "    dϵ: Array{Float64, 6}\n",
    "        Strain rate vector in Voigt notation\n",
    "    C: Array{Float64, (6, 6)}\n",
    "        Material tensor\n",
    "    k: Float\n",
    "        Material constant, yield limit\n",
    "    Δt: Float\n",
    "        time increment\n",
    "    σ_begin:Array{Float64, 6}\n",
    "        Stress vector in Voigt notation\n",
    "\n",
    "Returns\n",
    "-------\n",
    "    Array{Float64, 7}, return values for solver\n",
    "\"\"\"\n",
    "function G(params, dϵ, C, k, Δt, σ_begin)\n",
    "    # Creating wrapper for gradient\n",
    "    yield(pars) = vonMisesYield(pars, k)\n",
    "    dfdσ = ForwardDiff.gradient(yield)\n",
    "    \n",
    "    # Stress rate\n",
    "    dσ = (σ_begin - params[1:6]) / Δt\n",
    "    \n",
    "    # Calculating plastic strain rate\n",
    "    dϵp = params[end] * dfdσ(params)\n",
    "    \n",
    "    # Calculating equations\n",
    "    function_1 = dσ - C * (dϵ - dϵp[1:6])\n",
    "    function_2 = yield(params)\n",
    "    [function_1[:]; function_2]\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "Function which calculates the stress. Also handles if any yielding happens\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "    dϵ: Array{Float64, 6}\n",
    "        Strain rate vector in Voigt notation\n",
    "    Δt: Float\n",
    "        time increment\n",
    "    σ: Array{Float64, 6}\n",
    "        Last stress vector in Voigt notation\n",
    "    C: Array{Float64, (6, 6)}\n",
    "        Material tensor\n",
    "    k: Float\n",
    "        Material constant, yield limit\n",
    "\n",
    "Returns\n",
    "-------\n",
    "    Tuple\n",
    "    Plastic strain rate dϵᵖ and new stress vector σ\n",
    "\"\"\"\n",
    "function calculate_stress(dϵ, Δt, σ, C, k)\n",
    "\n",
    "    # Test stress\n",
    "    σ_tria = σ + Δt * C * dϵ\n",
    "    \n",
    "    # Calculating yield\n",
    "    yield = vonMisesYield(σ_tria, k)\n",
    "    \n",
    "    if yield > 1\n",
    "        \n",
    "        # Yielding happened\n",
    "        # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values\n",
    "        initial_guess = [σ[:]; 0.0]\n",
    "        f(σ_)  = G(σ_, dϵ, C, k, Δt, σ_tria)\n",
    "        df     = ForwardDiff.jacobian(f)\n",
    "        \n",
    "        # Calculating root \n",
    "        result = nlsolve(not_in_place(f, df), initial_guess).zero\n",
    "        \n",
    "        σ      = result[1:6]        \n",
    "        dλ     = result[end]\n",
    "        \n",
    "        # these are not necessary to calculate for our calculation but why not?\n",
    "        yield_f(σ_) = vonMisesYield(σ_, k)\n",
    "        dfdσ_ = ForwardDiff.gradient(yield_f)\n",
    "        dϵᵖ    = dλ * dfdσ_(σ)\n",
    "    else\n",
    "        dϵᵖ = zeros(Float64, (6, 1))\n",
    "        σ = σ_tria\n",
    "    end\n",
    "    return (dϵᵖ, σ)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Defining strain history"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "linspace(-0.0,-0.0006,200)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "steps = 200\n",
    "\n",
    "ϵ_tot = zeros(Float64, (steps, 6))\n",
    "\n",
    "# Adding only strain in x-axis and counting for the poisson effect\n",
    "ϵ_tot[:, 1] = linspace(0, 0.002, steps)\n",
    "ϵ_tot[:, 2] = linspace(0, 0.002, steps).*-ν\n",
    "ϵ_tot[:, 3] = linspace(0, 0.002, steps).*-ν"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation\n",
    "\n",
    "Ok, we're good to go! Now we just need to define yield limit and the main loop.\n",
    "\n",
    "This simulation is not time dependent, but since it's already defined in the equations we'll give it value 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtcAAAI6CAYAAADsVGQBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xt8jvXjx/H3fbNTDrPN+ayvQyjNKSH54os55TRzzlkh852kkKSJpRS1kqKcySnHWA6LkcMUKodSiS+mzJxtxnb9/rhz/9xtmLm36972ej4e9yO77s997X3d3a23a5/rc1kMwzAEAAAA4IFZzQ4AAAAAZBeUawAAAMBJKNcAAACAk1CuAQAAACehXAMAAABOQrkGAAAAnIRyDQAAADgJ5RoAAABwEso1AAAA4CSUawBAmpUtW1blypUzOwYAuCzKNQBTJCUl6dNPP1XDhg3l6+srd3d3FSlSRI8//rgGDBigNWvWOIyfPXu2rFar5syZY1LijJeQkKB33nlHderUkbe3tzw8PFS8eHHVqlVLQ4cO1bZt2xzGv/7667JarSm2ZySLxSKLxZJp3w8AsprcZgcAkPMkJSWpdevWioiIkI+Pj1q3bq2SJUsqMTFRP/30kxYuXKiff/5Zbdq0SfHa7Frsrly5ooYNG2rfvn0qVqyYOnXqpKJFi+rKlSvav3+/PvnkE128eFFPP/20qTm3bNli6vcHAFdHuQaQ6RYtWqSIiAj5+/tr69atypcvn8Pz8fHx2rNnT6qvNQwjMyJmuqlTp2rfvn1q3ry51qxZo9y5HX88X7hwQUeOHEn1tZn5njAlBADujmkhADLdt99+K0nq3bt3imItSV5eXmrYsKH963//+9/q27evJKlPnz6yWq32x4kTJyT9/xSJrVu3auHChapTp47y5s3rUAavXbumSZMmyd/fX3nz5lW+fPlUr149LV68ONWcc+bMUb169VSoUCF5eXmpdOnSCggI0JIlSxzG/fDDD+ratavKli0rT09PFS5cWDVr1lRISIhu3rx5X+/JoEGDUhRrSSpQoICefPJJ+9dly5bVG2+8IUlq1KiRw3tyS+/evWW1WnXs2DF98MEHqlatmh566CE1atRIknTjxg2Fh4erZcuWKlOmjDw9PeXn56emTZtqw4YNqeZMbc717VN2IiMj9e9//1v58+eXt7e3Wrdufce/FNzN119/rTZt2qhw4cLy9PRU6dKl1a5dO23evDnV75saq9VqP9Zb7vY52b17t6xWqzp06HDHXJUrV5anp6cuXLjgsD0iIkItW7ZUwYIF5enpqfLly2vkyJG6ePHifR87gKyNM9cAMl3BggUlST///HOaxvfp00c+Pj5atWqV2rVrJ39/f/tz3t7eDmOnTJmijRs36plnnlGTJk3s5ebChQtq3Lix9u/fr5o1a6pfv35KTk7Whg0b1K1bNx08eFChoaH2/YwePVphYWF6+OGH1aVLF3l7e+v06dOKjo7WsmXLFBQUJMlWrOvUqaNcuXLpmWeeUbly5XTp0iUdPXpU06dP15tvvplqWX7Q9yQkJEQrV67U1q1b1bt3b5UtW/aOY4cNG6aoqCi1bt1arVu3Vq5cuSRJ586d03//+1/Vr19fzZs3V6FChXT69GmtWbNGLVu21Keffqp+/fql2N+dpuasXbtWq1atUsuWLTVo0CAdPHhQX331laKjo3Xo0CH5+fml6djGjRun0NBQ5cuXT+3atVOpUqV06tQpffvtt1qwYIGaNGmSpjx3ey61z0mdOnVUqVIlffXVV4qLi5Ovr6/Da/bs2aOff/5ZgYGBKlCggH37+PHjNX78ePn5+dn/QnDgwAG98847+uqrr7Rz585U/xIJIJsyACCT7du3z3B3dzesVqvRs2dPY8WKFcYff/xx19d8/vnnhsViMebMmZPq8+PGjTMsFouRN29eY//+/Sme79Wrl2GxWIy3337bYXtCQoIREBBgWK1Wh9f5+voapUqVMuLj41PsKzY21v7n4cOHGxaLxVi9enWKcRcuXDCSk5Pvely3rF271rBYLIaHh4cxePBgY926dcbp06fv+ppbx7x169ZUn791zCVLlkz1/b1+/bpx6tSpFNsvXrxoPProo4avr2+K4y9TpoxRrlw5h223/t24ubkZW7ZscXhu1KhRhsViMSZPnnzXY7klIiLCsFgsxr/+9a9Uj//kyZMpvu+dPhMWi8Vo1KiRw7Z7fU4mTZpkWCwWIzw8PMVzgwcPNiwWi7F27Vr7ti1bthgWi8WoX7++cfHiRYfxs2fPNiwWixESEnL3gwaQrTAtBECm8/f31/z581WkSBHNnz9fHTt2VLly5eTn56cOHTpo7dq16d73wIED9fjjjztsO3funObPn6/atWtrxIgRDs95eHgoLCxMhmFo4cKF9u0Wi0Vubm4O0yxuSe0MrKenZ4pt3t7eab4As1WrVpo2bZq8vLw0ffp0tW7dWiVKlFCxYsXUo0cPRUVFpWk/qRk5cqTKlCmTYru7u7uKFy+eYnv+/PnVp08fnT9/XtHR0Wn+Pl26dEkxDWPgwIGSlOb9fPDBB5JsZ5aLFSuW4vkSJUqkOc/dpPY5kaSePXumOtUkMTFRixcvVpEiRdSiRQv79vfff1+S9Omnnyp//vwOr+nVq5cef/xxLViwwCmZAWQNTAsBYIpOnTqpffv2ioyM1I4dO7Rv3z5t375dK1eu1MqVK/Xss89q9uzZ973fJ554IsW26OhoJScnS7LNuf2nGzduSJIOHz5s39a9e3d98MEHqlKlioKCgtSwYUM9+eSTKaahdOnSRe+//77atWunwMBANWnSRPXr19e//vUvh3H79+/XypUrHbb5+Pho2LBh9q+HDh2q/v37a+PGjdq5c6f27dunb7/9VgsXLtTChQs1duxYjR8//v7eEKX+ntxy8OBBvf3229q2bZvOnDmjhIQEh+dPnz6d5u9Tq1atFNtKliwpSTp//nya9rFr1y5ZrVYFBASk+fumx53ekxIlSqhJkybauHGjDh8+rMqVK0uS1qxZo/Pnz2v48OEOf+HauXOn3NzctGTJklQvLE1MTNTZs2d1/vx5+fj4ZMzBAHAplGsApsmdO7eaNm2qpk2bSpKSk5O1fPly9e3bV3PnzlX79u3Vtm3b+9pn0aJFU2w7d+6cJFvJvtMZVIvFoqtXr9q/fu+99/Twww/r888/V1hYmMLCwpQ7d261bNlSU6ZMsZfn2rVrKyoqSm+++aaWLVumefPmSZIqVaqkcePGqUuXLpKkAwcO2C9AvKVs2bIO5VqyXcz5zDPP6JlnnpFkK/6ffvqphg0bptDQUHXo0CHVM673+55ItiLbuHFjJScnq0mTJmrXrp3y588vq9Wqffv2adWqVbp+/Xqav8/t85BvuTXfPCkpKU37uHDhgnx8fOTh4ZHm75sed3pPJNuFoBs3btScOXMUFhYmSfYz2b169XIYe+7cOSUlJd31Lz0Wi0VXrlyhXAM5BNNCALgMq9WqTp06KSQkRJIUGRl53/tIbRrGrbPNw4cPV3JycqqPpKQkh5UorFarhg0bpv379+vPP//U8uXL1b59e61evVoBAQFKTEy0j33yySe1Zs0aXbhwQTt27NDYsWP1559/qlu3bvZ99urVK8X3/P333+95PG5ubho8eLC6du0qKX3rTN9pasqECROUkJCgr7/+WuvWrdO7776r119/Xa+99tpdz3ZnpAIFCuj8+fMpzqCn5tYZ5NRWZPnnah7/dLfpOu3bt1f+/Pk1f/58GYahv/76S+vXr5e/v78ee+wxh7He3t7y9fW94+fq1merVKlS9zweANkD5RqAy8mbN68kx/Wbb61wkdYzoLerU6fOA93JsFChQmrfvr2++OILNWrUSL/99psOHjyYYpybm5vq1q2r8ePH2+firl69Ol3f859uvSe3e5D3RJJ+/fVX+fn5pXpjmq1bt6Zrnw+qbt269lVc7uXWmeBbyzHebu/evenO4OnpqaCgIJ0+fVobN27UwoULlZSUlOKs9a28cXFxOnToULq/H4DshXININMtWrRImzZtSnWO6pkzZ/Tpp59KkkPpu3UR4fHjx+/7+xUqVEjdu3fX3r17NWHCBPv869v99ttv+uOPPyTZ5snu2LEjxZgbN24oLi5OFotFDz30kCTb+tSpnWU9c+aMJNnH3cvHH3+s3bt3p/rckSNHtHTpUlksFqe9J5LthjDnzp3Tjz/+6LB91qxZ+vrrr9O1zwc1dOhQSdKLL76Y6nzv27fVrl1bVqtVCxcuVHx8vH17XFycRo4c+UA5evfuLUmaO3eu5s6dKzc3N3Xv3j3FuFu/ZRkwYIBiYmJSPH/16tU7/nsFkD0x5xpAptuzZ4+mTZumokWL6qmnnrKv0Xzs2DGtW7dOCQkJateunTp27Gh/Tb169fTQQw9p6tSpOnfunIoUKSJJCg4OTrFKQ2rCw8N19OhRvfbaa5o3b57q16+vIkWK6PTp0zp8+LD27t2rxYsXq2zZsrp27ZoaNGig8uXLq0aNGipTpowSEhK0ceNGHTlyRG3btlWlSpUkSZMnT1ZkZKQaNGigsmXLKm/evDp48KA2bNggX19f+2oZ9xIREaHBgwerbNmyql+/vkqWLKnr16/r6NGjioiIUFJSkoKDg1WzZk37axo3biyr1apRo0bpxx9/tJ/JffXVV9P0Pf/73/8qIiJCTz31lIKCgpQ/f37t3btXO3bsUGBgoJYtW5am/ThT06ZN9eqrr2rChAmqXLmy2rVrp5IlS+rPP//U9u3bVbduXX3++eeSbPOmu3fvrnnz5snf318tW7bUpUuXtH79ejVs2FD79+9Pd4569eqpfPnyWrp0qW7cuKFnnnnGvhb57Ro3bqywsDCNGjVKFSpUUMuWLVW2bFlduXJFx48f17Zt29SgQQN99dVX6c4CIIsxdSHAf9izZ48xZMgQo0qVKkaePHmM0qVLG0FBQcYvv/ziMO7W2q3/fDzyyCOp7nfmzJnGI488Ynh6ehoVKlQwPvjgg8w4HAB38L///c/48MMPjfbt2xuVKlUy8ufPb7i7uxvFixc3WrVqZSxYsCDV123YsMGoW7eukTdvXsNisRhWq9U4fvy4YRiG8frrrxtWq/WOaz4bhmEkJiYa4eHhRr169Qxvb2/Dw8PDKFOmjPGf//zHmDZtmnHu3DnDMAzjxo0bxuTJk40WLVoYpUuXNjw9PY3ChQsbdevWNWbMmGHcuHHDvs+vv/7a6NOnj1GlShXD29vbyJMnj/HII48Yw4YNM06cOJHm9+SXX34xpkyZYrRo0cIoX768kSdPHnu+jh07GuvWrUv1dfPnzzf8/f0NLy8v+3tyS+/evR3eo9SsXbvWePLJJ418+fIZPj4+RvPmzY2oqChj9uzZhtVqTbGGdNmyZVOsc32nsbektt70vXz11VdGQECA4evra3h4eBilS5c2OnToYERGRjqMu379uvHSSy8ZJUuWNNzd3Y0KFSoYYWFhxs2bN1P9vmn5nNwyYcIE+3u6YsWKu47dvn27ERQUZBQvXtxwd3c3ChcubFSvXt148cUXje++++6+jh1A1mYxjFR+L2uSwMBA7dy5U506dVK1atUUExOj8PBwXblyRbt27VLVqlUl2X5d98UXX2jWrFkOr/f29larVq0cts2YMUODBg1SYGCgmjdvrm3btmnevHkKCwt74F8bAgAAALdzqXK9c+dO1a5d2+FWwb/++qsee+wxBQYG2pe46t27t1asWKFLly7ddX/x8fEqVaqU6tWr53BRUc+ePbVy5Ur973//S3XpKAAAACA9XOqCxrp16zoUa0kqX768qlSpoiNHjjhsNwxDycnJdy3YkZGRiouL0+DBgx22DxkyRFevXtW6deucFx4AAAA5nkuV69QYhqE///wzxYUk165dU/78+VWgQAH5+fnphRdecLgBhCTt27dPUsq7htWoUUNWq/WBLnYBAAAA/snlVwtZsGCBTp8+rQkTJti3FS9eXC+//LJq1Kih5ORkrV+/Xh999JEOHDigb775xr72a0xMjHLlypWimLu7u8vPz+++busLAAAA3ItLl+sjR45oyJAhqlevnsPi/RMnTnQYFxQUpIoVK2rMmDFatmyZOnfuLMk259rd3T3VfXt4eDisiwoAAAA8KJct12fOnFGrVq3k4+OjZcuW3fVWtZJtIf+xY8dq8+bN9nLt5eXlcIvi2yUkJMjLyyvV52JjYxUREaGyZcvecQwAAADMEx8frz/++EPNmzdPdR16s7hkub548aJatGihS5cuKSoqSkWLFr3nazw9PeXr66u4uDj7tmLFiikpKUmxsbEOb3piYqLi4uJUvHjxVPcVERGhHj16PPiBAAAAIEPNnz8/1TuomsXlynVCQoLatGmjX3/9VZs2bdIjjzySptddvnxZsbGxKlSokH1b9erVJUnR0dFq0aKFffvevXuVnJwsf3//VPd1625x8+fPV+XKldN5JEDahISE6L333jM7BnIAPmvILHzWkBkOHz6sHj162Hubq3Cpcp2UlKTOnTtr9+7dWrVqlerUqZNizPXr15WYmKh8+fI5bA8NDZUkBQQE2Lc1btxYvr6+mj59ukO5nj59uvLkyZPihjO33JoKUrlyZdWoUeOBjwu4G29vbz5nyBR81pBZ+KwhM7naFF6XKtcvvvii1qxZozZt2ig2Nlbz5893eL5Hjx6KiYlR9erV1a1bN1WqVEmSbRrH+vXr1aJFC7Vt29Y+3tPTU6GhoRoyZIiCgoLUrFkzRUVFacGCBZo4cSI3kAEAAIBTuVS5PnDggCwWi9asWaM1a9Y4PGexWNSjRw/5+PioTZs22rhxo+bMmaOkpCRVqFBBkyZN0ogRI1Lsc9CgQXJzc9OUKVO0evVqlS5dWlOnTlVwcHBmHRYAAAByCJcq15GRkfcc4+3trblz597Xfvv376/+/funNxYAAACQJi5/h0Ygu+vatavZEZBD8FlDZuGzhpyMcg2YjP8JIbPwWUNm4bOGnIxyDQAAADgJ5RoAAABwEso1AAAA4CSUawAAAMBJKNcAAACAk1CuAQAAACehXAMAAABOQrkGAAAAnIRyDQAAADgJ5RoAAABwEso1AAAA4CSUawAAAMBJKNcAAACAk1CuAQAAACehXAMAAABOQrkGAAAAnIRyDQAAADgJ5RoAAABwEso1AAAA4CSUawAAAMBJKNcAAACAk1CuAQAAACehXAMAAABOQrkGAAAAnIRyDQAAADgJ5RoAAABwEso1AAAA4CSUawAAAMBJKNcAAACAk1CuAQAAACehXAMAAABOQrkGAAAAnIRyDQAAADgJ5RoAAABwEso1AAAA4CSUawAAAMBJKNcAAACAk1CuAQAAACehXAMAAABOQrkGAAAAnIRyDQAAADgJ5RoAAABwEso1AAAA4CSUawAAAMBJKNcAAACAk1CuAQAAACehXAMAAABOQrkGAAAAnIRyDQAAADgJ5RoAAABwEso1AAAA4CSUawAAAMBJKNcAAACAk1CuAQAAACehXAMAAABOQrkGAAAAnIRyDQAAADgJ5RoAAABwEso1AAAA4CSUawAAAMBJKNcAAACAk1CuAQAAACehXAMAAABOQrkGAAAAnIRyDQAAADgJ5RoAAABwEso1AAAA4CSUawAAAMBJKNcAAACAk1CuAQAAACehXAMAAABOQrkGAAAAnIRyDQAAADgJ5RoAAABwEso1AAAA4CSUawAAAMBJKNcAAACAk1CuAQAAACehXAMAAABOQrkGAAAAnMSlynV0dLReeOEFVa1aVXnz5lWZMmXUuXNnHT16NMXYw4cPKyAgQPny5ZOfn5+effZZxcbGprrfWbNmqXLlyvLy8lLFihUVHh6e0YcCAACAHCi32QFu99Zbb2nnzp3q1KmTqlWrppiYGIWHh6tGjRratWuXqlatKkk6efKknn76afn4+GjSpEm6fPmy3nnnHf3444/as2eP3Nzc7PucMWOGBg0apMDAQI0YMULbtm1TcHCwrl27ppEjR5p1qAAAAMiGXKpcv/jii6pdu7Zy5/7/WJ07d9Zjjz2msLAwzZs3T5I0ceJExcfHa9++fSpZsqQk6YknnlDTpk01e/ZsDRgwQJIUHx+vMWPGqHXr1lqyZIkkqV+/fkpOTlZoaKgGDhyoAgUKZPJRAgAAILtyqWkhdevWdSjWklS+fHlVqVJFR44csW9bvny5WrdubS/WktSkSRNVrFjRXqIlKTIyUnFxcRo8eLDDPocMGaKrV69q3bp1GXQkAAAAyIlcqlynxjAM/fnnnypYsKAk6dSpUzp79qxq1aqVYmzt2rW1b98++9e3/vzPsTVq1JDVatX+/fszMDkAAAByGpcv1wsWLNDp06fVuXNnSVJMTIwkqVixYinGFitWTHFxcbpx44Z9bK5cuezF/BZ3d3f5+fnp9OnTGZweAAAAOYlLl+sjR45oyJAhqlevnnr16iXJNo9akjw8PFKM9/T0dBgTHx8vd3f3VPft4eFhHwcAAAA4g0td0Hi7M2fOqFWrVvLx8dGyZctksVgkSV5eXpKk69evp3hNQkKCwxgvLy8lJiamuv+EhAT7OADICb7/Xlq+XPr7RyUAZEmNG0utWpmd4s5cslxfvHhRLVq00KVLlxQVFaWiRYvan7s1HeTW9JDbxcTEyM/Pz74UX7FixZSUlKTY2FiHqSGJiYmKi4tT8eLF75ojJCRE3t7eDtu6du2qrl27pvvYACAzJSVJCxdKH34o7d4tFS4s+fmZnQoA7s/Fi4t06dIiSdKGDdKMGba+6IpcrlwnJCSoTZs2+vXXX7Vp0yY98sgjDs+XKFFChQoVUnR0dIrX7tmzR/7+/vavq1evLsl2c5oWLVrYt+/du1fJyckOY1Pz3nvvqUaNGg9yOABgqhEjpKlTpWbNpJUrbWd7crvcT34AuJeufz/+3/fff6+aNWuaE+cuXGrOdVJSkjp37qzdu3dr6dKlqlOnTqrjOnbsqLVr1+rkyZP2bZs3b9bRo0fVqVMn+7bGjRvL19dX06dPd3j99OnTlSdPHrVy5d8pAMADmjvXVqzff1+KiJDatqVYA0BGc6kfsy+++KLWrFmjNm3aKDY2VvPnz3d4vkePHpKk0aNHa+nSpWrUqJGGDRumy5cv6+2331a1atXUp08f+3hPT0+FhoZqyJAhCgoKUrNmzRQVFaUFCxZo4sSJ3EAGQLYVHS0NHCj16SO98ILZaQAg53Cpcn3gwAFZLBatWbNGa9ascXjOYrHYy3XJkiW1detWDR8+XK+88oo8PDzUunVrTZkyxeHW55I0aNAgubm5acqUKVq9erVKly6tqVOnKjg4ONOOCwAy05kzUvv2kr+/NH269Pf14ACATOBS5ToyMjLNY6tUqaINGzakaWz//v3Vv3//9MYCgCwjMVEKDLRdyLhihZTKqqUAgAzkUuUaAPBggoNtU0K++Ua6x4JIAIAMQLkGgGxixgzbY+ZMqW5ds9MAQM7kUquFAADSZ/t2aehQacgQqV8/s9MAQM5FuQaALO7kSds867p1pffeMzsNAORslGsAyMLi420rg7i7S0uXSv9YMAkAkMmYcw0AWZRhSM8/L/30k7Rjh+3W5gAAc1GuASCLmjbNdhfGBQukGjXMTgMAkJgWAgBZ0ubN0ogRtke3bmanAQDcQrkGgCzm2DEpKEhq0kQKCzM7DQDgdpRrAMhCrl6V2rWTChSQFi2ScuUyOxEA4HbMuQaALMIwpD59pN9+k3btknx9zU4EAPgnyjUAZBFhYbbl9pYvlx591Ow0AIDUMC0EALKAdeukMWOksWOlDh3MTgMAuBPKNQC4uJ9/tq0I0qaN9PrrZqcBANwN5RoAXNjFi1LbtlLx4tK8eZKVn9oA4NKYcw0ALio5WerZUzpzRtqzR8qf3+xEAIB7oVwDgIsaN05au9b2qFjR7DQAgLSgXAOAC1q+XJowQZo0SWrZ0uw0AIC0YvYeALiYH3+UevWy3YXx5ZfNTgMAuB+UawBwIXFxtjswli8vffaZZLGYnQgAcD+YFgIALuLmTalLF9sKIZs3S3nymJ0IAHC/KNcA4CJeeUXaskX6+mupbFmz0wAA0oNyDQAuYMECacoUaepUqXFjs9MAANKLOdcAYLLvvpP697ddxBgcbHYaAMCDoFwDgIn++ktq31567DHp44+5gBEAsjrKNQCYJDFRCgy0/XPFCsnT0+xEAIAHxZxrADBJSIi0a5cUGSmVLGl2GgCAM1CuAcAEM2dKH30kzZgh1a9vdhoAgLMwLQQAMtm330qDB0vPPy8NHGh2GgCAM1GuASATnToldewo1akjTZtmdhoAgLNRrgEgkyQkSB06SLlzS8uWSe7uZicCADgbc64BIBMYhjRokPTDD1JUlFSkiNmJAAAZgXINAJkgPFyaPVuaO1eqVcvsNACAjMK0EADIYJGRtmX3QkKknj3NTgMAyEiUawDIQH/8IXXqJDVqJE2ebHYaAEBGo1wDQAa5dk1q107Kn19avNh2ISMAIHvjRz0AZADDkPr2lY4elXbulPz8zE4EAMgMlGsAyACTJ0tffCEtXSpVq2Z2GgBAZmFaCAA42YYN0qhR0pgxUmCg2WkAAJmJcg0ATnT0qNSli9SypfTGG2anAQBkNso1ADjJpUtS27ZS0aLSggWSlZ+wAJDjMOcaAJwgOVl69lnp1Clpzx7J29vsRAAAM1CuAcAJ3nhDWr3a9qhUyew0AACzUK4B4AF9+aU0frw0YYLUurXZaQAAZmJGIAA8gIMHbdNBAgOl0aPNTgMAMBvlGgDS6fx52x0Yy5WTPv9csljMTgQAMBvTQgAgHZKSbEvuxcVJ0dFS3rxmJwIAuALKNQCkw+jR0qZNUkSE9PDDZqcBALgKyjUA3KdFi2y3N58yRfrPf8xOAwBwJcy5BoD7sG+f1K+f1KOHFBJidhoAgKuhXANAGp09a7uAsUoV6ZNPuIARAJAS5RoA0uDGDalTJykhwbautZeX2YkAAK6IOdcAkAbDh0s7dkhbtkilSpmdBgDgqijXAHAPn30mhYdL06dLDRqYnQYA4MqYFgIAd7FrlzRokDRggPSJKSn9AAAgAElEQVTcc2anAQC4Oso1ANzB6dNShw5SrVq2M9dcwAgAuBfKNQCk4vp1qWNHW6Fevlxydzc7EQAgK2DONQD8g2FIQ4bY1rTetk0qWtTsRACArIJyDQD/8NFH0qxZ0uzZ0hNPmJ0GAJCVMC0EAG6zdav03/9Kw4ZJvXqZnQYAkNVQrgHgbydO2G4U06CB9PbbZqcBAGRFlGsAkHTtmu3W5nnySEuWSG5uZicCAGRFzLkGkOMZhm0d6yNHpJ07pYIFzU4EAMiqKNcAcrwpU6SFC6XFi6XHHzc7DQAgK2NaCIAcLSJCevll6ZVXpM6dzU4DAMjqKNcAcqxff5W6dJGaN5cmTDA7DQAgO6BcA8iRLl+2XcBYqJBtSkiuXGYnAgBkB8y5BpDjJCfb1rA+cULavVsqUMDsRACA7IJyDSDHefNN6csvpVWrpMqVzU4DAMhOmBYCIEdZvVp67TVp/HjpmWfMTgMAyG4o1wByjMOHpR49pPbtpVdfNTsNACA7olwDyBEuXJDatpVKl5bmzJGs/PQDAGQA5lwDyPaSkqRu3aSzZ6XoaClfPrMTAQCyK8o1gGzv1VdtN4tZv14qX97sNACA7IxyDSBbW7JECguT3n5batbM7DQAgOyOWYcAsq0DB6Q+faSuXaUXXzQ7DQAgJ6BcA8iWYmNtFzBWqiTNnClZLGYnAgDkBEwLAZDt3LwpBQVJ165J27ZJDz1kdiIAQE5BuQaQ7YwYIUVFSZs22ZbeAwAgs1CuAWQrc+ZI06ZJ4eFSw4ZmpwEA5DTMuQaQbezZIz33nNSvnzR4sNlpAAA5EeUaQLZw5ozUoYNUvbr04YdcwAgAMIfLleurV69q3LhxCggIkK+vr6xWq+bMmZNiXO/evWW1WlM8KleunOp+Z82apcqVK8vLy0sVK1ZUeHh4Rh8KgExy/brUsaOUnCwtXy55eJidCACQU7ncnOuzZ88qNDRUZcqUkb+/v7755htZ7nAKysPDQ7NmzXLY5u3tnWLcjBkzNGjQIAUGBmrEiBHatm2bgoODde3aNY0cOTJDjgNA5gkOlvbulbZulYoXNzsNACAnc7lyXbx4cZ05c0aFCxfWd999p9q1a99xrJubm7p163bX/cXHx2vMmDFq3bq1lixZIknq16+fkpOTFRoaqoEDB6pAgQJOPQYAmefjj6VPPpFmzZKefNLsNACAnM7lpoW4u7urcOHCkiTDMO461jAMJScn69KlS3ccExkZqbi4OA3+x9VNQ4YM0dWrV7Vu3boHDw3AFFFR0tCh0gsvSH37mp0GAAAXLNf349q1a8qfP78KFCggPz8/vfDCC7p69arDmH379kmSatWq5bC9Ro0aslqt2r9/f6blBeA8//ufFBgo1a8vvfuu2WkAALBxuWkhaVW8eHG9/PLLqlGjhpKTk7V+/Xp99NFHOnDggL755hvlypVLkhQTE6NcuXKpYMGCDq93d3eXn5+fTp8+bUZ8AA8gPl5q317y9JSWLpXc3MxOBACATZYt1xMnTnT4OigoSBUrVtSYMWO0bNkyde7cWZJtzrW7u3uq+/Dw8FB8fHyGZwXgPIYhDRwoHTok7dghFSpkdiIAAP5flp4W8k8hISGyWq3avHmzfZuXl5cSExNTHZ+QkCAvL6/MigfACaZOlebPt13AWL262WkAAHCUZc9cp8bT01O+vr6Ki4uzbytWrJiSkpIUGxvrMDUkMTFRcXFxKn6XdbtCQkJSLO3XtWtXde3a1fnhAdzTpk3SiBHSSy9J/GcIADnHokWLtGjRIodtFy9eNCnN3WWrcn358mXFxsaq0G2/J67+96mt6OhotWjRwr597969Sk5Olr+//x33995776lGjRoZFxhAmv3+u9S5s9S0qTRpktlpAACZKbWTm99//71q1qxpUqI7y5LTQq5fv67Lly+n2B4aGipJCggIsG9r3LixfH19NX36dIex06dPV548edSqVauMDQvggV25IrVrJ/n4SIsWSX9frwwAgMtxyTPX4eHhunDhgn0lj9WrV+vEiROSpODgYMXFxal69erq1q2bKlWqJEmKiIjQ+vXr1aJFC7Vt29a+L09PT4WGhmrIkCEKCgpSs2bNFBUVpQULFmjixIncQAZwcYYh9ekjHTsm7dplK9gAALgqlyzXU6ZM0fHjxyVJFotFX375pVasWCGLxaJnn31WPj4+atOmjTZu3Kg5c+YoKSlJFSpU0KRJkzRixIgU+xs0aJDc3Nw0ZcoUrV69WqVLl9bUqVMVHByc2YcG4D5NnCgtWyZ9+aVUtarZaQAAuDuXLNfHjh2755i5c+fe1z779++v/v37pzcSABOsXSuNHSuNG2ebFgIAgKvLknOuAWR/R45I3btLzzwjvfaa2WkAAEgbyjUAl3Pxou1MdYkS0ty5kpWfVACALMIlp4UAyLmSk21nrM+ckaKjpfz5zU4EAEDaUa4BuJTXXpO++sr2qFDB7DQAANwfyjUAl7FsmfTmm1JYmHTbcvUAAGQZzGQE4BJ++EHq1ct2F8aRI81OAwBA+lCuAZju3DnbBYwVKkizZkkWi9mJAABIH6aFADDVzZu2s9WXL0tbtkh58pidCACA9KNcAzDVyy9L33wjbdwolS1rdhoAAB4M5RqAaebNk959V5o2TWrUyOw0AAA8OOZcAzDF3r3SgAFS797S0KFmpwEAwDko1wAy3Z9/Su3bS48/Lk2fzgWMAIDsg3INIFMlJkqBgbYLGVeskDw9zU4EAIDzMOcaQKYaNkzavdt2EWOJEmanAQDAuSjXADLNJ59IH38sffqpVK+e2WkAAHA+poUAyBQ7dkgvvCANGiT17292GgAAMgblGkCGO3lS6thRevJJaepUs9MAAJBxKNcAMlRCgtShg+TmJi1bJrm7m50IAICMw5xrABnGMKTnn5d+/FHavl0qXNjsRAAAZCynlevk5GRFRkYqMTFRTz31lPLly+esXQPIot5/X5ozR5o/X6pZ0+w0AABkvHRNCxkzZowa3XavYsMw1KxZMzVt2lStWrXSo48+qt9++81pIQFkPVu2SC++aHt07252GgAAMke6yvXy5ctVu3Zt+9fLli3Tli1b9Oabb2rt2rVKSkrSuHHjnBYSQNZy7JgUFCQ1biyFhZmdBgCAzJOuaSGnTp1ShQoV7F+vWLFClStX1qhRoyRJgwcP1vTp052TEECWcvWq1K6d5O0tLV4s5ebKDgBADpKuM9e5c+fW9evXJdnmWm/evFkBAQH25wsXLqyzZ886JyGALMMwpL59pd9+k1aulHx9zU4EAEDmSle5rlq1qubNm6e4uDjNnj1bsbGxatWqlf35EydOqGDBgk4LCSBreOstackS20WMjz1mdhoAADJfun5hO27cOLVu3dpeoOvXr+9wgeO6desc5mQDyP6++koaPVp69VXbDWMAAMiJ0lWumzZtqu+//14bN26Uj4+POnfuLIvFIkk6f/68GjRooLZt2zo1KADX9csvUrduUqtW0vjxZqcBAMA86b7UqGrVqqpatWqK7T4+PprK/Y2BHOPSJdsFjEWL2taztnLfVwBADpaucn3p0iVdvHhRpUqVsm87deqUPv74YyUmJqpjx4564oknnBYSgGtKTpZ69JBOnZL27LGtEAIAQE6WrnL93HPP6dixY9q1a5ckW9l+8sknderUKVmtVk2bNk3r1693mIcNIPsZP15au1Zas0aqVMnsNAAAmC9dv8Ddvn27w+og8+fPV0xMjHbs2KHz58+rWrVqevPNN50WEoDrWbFCeuMNacIE21xrAACQznIdGxurkiVL2r9evXq16tevr7p16ypfvnzq2bOn9u/f77SQAFzLTz9Jzz4rdeok/X3vKAAAoHSW6wIFCujMmTOSpPj4eEVFRalZs2b253Pnzq1r1645JyEAlxIXJ7VtKz38sPT559LfCwUBAAClc851vXr19NFHH6lSpUrasGGD4uPjHZbeO3r0qEqUKOG0kABcw82bUteu0oUL0qZNUp48ZicCAMC1pKtch4WFqXnz5goMDJQkDR8+XI8++qgk6ebNm1qyZInD7dABZA+jRkmbN0sREVK5cmanAQDA9aSrXFeoUEFHjhzRoUOH5O3trXK3/V82Pj5e4eHh8vf3d1pIAOZbuFB65x3pvfekJk3MTgMAgGtK901k3N3dUy3Q+fLlU7t27R4oFADX8v33Ur9+Us+e0rBhZqcBAMB1pfteahcvXtSkSZPUrFkzVa9eXXv27JEkxcXF6d1339Wvv/7qtJAAzPPXX7Y7MD76qDRjBhcwAgBwN+k6c33y5Ek9/fTTOnnypMqXL68jR47oypUrkmy3P58xY4b++OMPvf/++04NCyBz3bhhW27v+nXpyy8lLy+zEwEA4NrSVa5feuklXblyRfv27VORIkVUuHBh+3MWi0Xt2rXT2rVrnRYSgDlCQqSdO6UtW6TblrYHAAB3kK5pIV9//bWGDh2qqlWrpvp8uXLldOLEiQcKBsBcs2ZJH34offCB9NRTZqcBACBrSFe5jo+Pdzhb/U+XL19OdyAA5tu5Uxo8WHruOdsDAACkTbrKdeXKlbV169Y7Pr9q1SpVr1493aEAmOf0aaljR6l2bYnLJgAAuD/pKtchISH64osv9NZbb+nixYuSpKSkJB09elQ9evTQt99+q5CQEKcGBZDxEhKkDh0kq1Vatkxydzc7EQAAWUu6Lmjs0aOHjh8/rjFjxmj06NGSpICAABmGoVy5cmnSpElq3769U4MCyFiGIQ0ZIu3fL0VFSUWLmp0IAICsJ903kRkzZox69OihFStW6OjRo0pOTlb58uXVoUMHPfzww87MCCATfPih9Nln0pw5tikhAADg/t13ub569aoaNGiggQMH6vnnn2f6B5ANfPON9N//2h7PPmt2GgAAsq77nnOdJ08e/fHHH7JwmzYgWzh+3HajmIYNpbffNjsNAABZW7ouaAwICFBERISzswDIZNeuSe3bS3nzSl98IeVO90QxAAAgpbNcjx07Vr/88ot69Oih7du369SpU4qLi0vxAOC6DEPq10/6+Wdp1SqpYEGzEwEAkPWl6zzVrTszHjp0SAsXLkx1jMViUVJSUvqTAchQ77wjLV4sLVkiVatmdhoAALKHdJXr1157zdk5AGSiiAjplVekUaNs860BAIBzpKtcv/76606OASCz/Pqr1KWLFBAghYaanQYAgOwlXXOu+/btq927d9/x+T179qhv377pDgUgY1y+LLVtKxUuLC1YIOXKZXYiAACyl3SV69mzZ+u333674/O///67Zs+end5MADJAcrJtDev//U9auVIqUMDsRAAAZD8ZsvBWTEyMvLy8MmLXANIpNNS2KsiqVVLlymanAQAge0pzuV61apVWrVolwzAkSZ988ok2bdqUYtz58+e1adMm1eb+yYDLWLVKev11W8Fu08bsNAAAZF9pLtcHDx7UkiVL7Hdm3L17t7777jv784ZhyGKxKE+ePGrYsKHeffdd56cFcN8OHZJ69JA6dJBGjzY7DQAA2Vuay/Xo0aM1+u//M1utVs2cOVPdu3fPsGAAHtz587YLGMuWlebMkazpusoCAACkVbrmXCcnJzt8ffjwYS1btkwxMTGqVKmS+vTpo/z58zslIID0SUqSunWTzp2ToqNttzgHAAAZK83lOjw8XO+//76+/fZbFbztPslr1qxRYGCgbty4Yd/2/vvva/fu3Q7jAGSuMWOkr7+WNmyQ/vUvs9MAAJAzpPmXxKtXr9bDDz/sUJhv3ryp/v37K1euXPr888/1ww8/KCwsTCdOnNCECRMyJDCAe1u8WHrrLWnyZKlpU7PTAACQc6S5XB86dEh169Z12BYZGamzZ88qJCREvXr10qOPPqqRI0cqKChI69evd3pYAPe2f7/Ut6/Uvbs0fLjZaQAAyFnSXK7PnTunUqVKOWy7tRRf+/btHbbXq1dPx48fd0I8APfj7FmpXTvbOtaffir9vbgPAADIJGmec12kSBGdOXPGYVtUVJQeeugh+fv7O2x3d3eXu7u7cxICSJMbN6SgIOnaNWnbNon7OAEAkPnSfOa6Vq1amjNnji5duiTJtu51dHS0mjdvrty5HTv6zz//rJIlSzo3KYC7GjFC2r5dWrZMKl3a7DQAAORMaT5zPW7cONWqVUsVK1ZU1apVtXfvXknSqFGjUoz98ssv1ahRI+elBHBXn38uvf++9OGH0tNPm50GAICcK81nrh977DFFRkaqZs2aOnXqlOrWrav169erVq1aDuMiIyPl5eWlTp06OT0sgJR275aef17q318aNMjsNAAA5Gz3dROZevXqad26dXcd06hRI/30008PFApA2sTE2G5rXrOmFB7OBYwAAJiNmyEDWdT161LHjrY/L18ueXiYmwcAAKTz9ucAzGUY0gsvSN99Z1sZpFgxsxMBAACJcg1kSR9/LM2cKX32mVSnjtlpAADALUwLAbKYbduk4GBp6FCpTx+z0wAAgNtRroEs5MQJKTBQatBAmjLF7DQAAOCfKNdAFhEfL7VvLz30kPTFF5Kbm9mJAADAPzHnGsgCDEMaMEA6fFj69lupUCGzEwEAgNRQroEs4N13pQULpMWLJX9/s9MAAIA7YVoI4OI2bpRGjpReflnq3NnsNAAA4G4o14AL++03W6Fu1kx6802z0wAAgHuhXAMu6soVqV07yc9PWrhQypXL7EQAAOBemHMNuCDDkHr1kv74Q9q9W/LxMTsRAABIC8o14ILefFNasUJauVKqUsXsNAAAIK1cblrI1atXNW7cOAUEBMjX11dWq1Vz5sxJdezhw4cVEBCgfPnyyc/PT88++6xiY2NTHTtr1ixVrlxZXl5eqlixosLDwzPyMIB0W7NGGjtWev11qW1bs9MAAID74XLl+uzZswoNDdXPP/8s/7/XHLNYLCnGnTx5Uk8//bR+//13TZo0SSNGjNC6devUtGlT3bhxw2HsjBkzNGDAAD322GMKDw9X3bp1FRwcrMmTJ2fKMQFpdeSI1L27ba712LFmpwEAAPfL5aaFFC9eXGfOnFHhwoX13XffqXbt2qmOmzhxouLj47Vv3z6VLFlSkvTEE0+oadOmmj17tgYMGCBJio+P15gxY9S6dWstWbJEktSvXz8lJycrNDRUAwcOVIECBTLn4IC7uHDBdqa6VClp7lzJ6nJ/9QUAAPficv/7dnd3V+HChSVJhmHccdzy5cvVunVre7GWpCZNmqhixYr2Ei1JkZGRiouL0+DBgx1eP2TIEF29elXr1q1z8hEA9y8pyXbG+q+/pFWrpHz5zE4EAADSw+XKdVqcOnVKZ8+eVa1atVI8V7t2be3bt8/+9a0//3NsjRo1ZLVatX///owNC6TBa69JGzbY7sBYvrzZaQAAQHplyXIdExMjSSpWrFiK54oVK6a4uDj7vOuYmBjlypVLBQsWdBjn7u4uPz8/nT59OuMDA3exdKk0caI0aZLUvLnZaQAAwIPIkuU6Pj5ekuTh4ZHiOU9PT4cx8fHxcnd3T3U/Hh4e9nGAGQ4ckHr3lrp2lV56yew0AADgQWXJcu3l5SVJun79eornEhISHMZ4eXkpMTEx1f0kJCTYxwGZ7dw526ogFStKM2dKqSyKAwAAshiXWy0kLW5NB7k1PeR2MTEx8vPzk5ubm31sUlKSYmNjHaaGJCYmKi4uTsWLF7/j9wkJCZG3t7fDtq5du6pr167OOAzkYDdvSkFBtlucf/ON9NBDZicCAMB1LVq0SIsWLXLYdvHiRZPS3F2WLNclSpRQoUKFFB0dneK5PXv22NfHlqTq1atLkqKjo9WiRQv79r179yo5Odlh7D+99957qlGjhhOTAzYvvSRt3Spt2iSVKWN2GgAAXFtqJze///571axZ06REd5Ylp4VIUseOHbV27VqdPHnSvm3z5s06evSoOnXqZN/WuHFj+fr6avr06Q6vnz59uvLkyaNWrVplWmZAsq1hPXWq7fHvf5udBgAAOJNLnrkODw/XhQsX7Ct5rF69WidOnJAkBQcHK3/+/Bo9erSWLl2qRo0aadiwYbp8+bLefvttVatWTX369LHvy9PTU6GhoRoyZIiCgoLUrFkzRUVFacGCBZo4cSI3kEGmio6WBg6U+vSRhgwxOw0AAHA2i3G3O7WYpFy5cjp+/Lik/7/1uWEYslgsOnbsmEqXLi1JOnTokIYPH67t27fLw8NDrVq10pQpU1SoUKEU+5w5c6amTJlif/0LL7yg4ODgVL//rV8zfPfdd0wLgdOcOSPVqiWVLGmbEpLKYjcAACCNXLWvueSZ62PHjqVpXJUqVbRhw4Y0je3fv7/69+//ILGAdEtMlAIDbXdiXLGCYg0AQHblkuUayG6Cg21TQr75RrrLAjUAACCLo1wDGWzGDNtj5kypbl2z0wAAgIyUZVcLAbKC7duloUNtFy/262d2GgAAkNEo10AGOXnSNs+6bl3pvffMTgMAADID5RrIAPHxUvv2kru7tHSp9PcNQwEAQDbHnGvAyQxDev556aefpB07pMKFzU4EAAAyC+UacLJp02x3YVywQHKhZTcBAEAmYFoI4ESbN0sjRtge3bqZnQYAAGQ2yjXgJMeOSUFBUpMmUliY2WkAAIAZKNeAE1y9KrVrJxUoIC1aJOXKZXYiAABgBuZcAw/IMKQ+faTffpN27ZJ8fc1OBAAAzEK5Bh5QWJhtub3ly6VHHzU7DQAAMBPTQoAHsG6dNGaMNHas1KGD2WkAAIDZKNdAOv38s21FkDZtpNdfNzsNAABwBZRrIB0uXpTatpWKF5fmzZOs/JcEAADEnGvgviUnSz17SmfOSHv2SPnzm50IAAC4Cso1cJ/GjZPWrrU9KlY0Ow0AAHAllGvgPixfLk2YIE2aJLVsaXYaAADgapgpCqTRjz9KvXrZ7sL48stmpwEAAK6Icg2kQVyc7Q6M5ctLn30mWSxmJwIAAK6IaSHAPdy8KXXpYlshZPNmKU8esxMBAABXRbkG7uGVV6QtW6Svv5bKljU7DQAAcGWUa+AuFiyQpkyRpk6VGjc2Ow0AAHB1zLkG7uC776T+/W0XMQYHm50GAABkBZRrIBV//SW1by899pj08cdcwAgAANKGcg38Q2KiFBho++eKFZKnp9mJAABAVsGca+AfQkKkXbukyEipZEmz0wAAgKyEcg3cZuZM6aOPpBkzpPr1zU4DAACyGqaFAH/79ltp8GDp+eelgQPNTgMAALIiyjUg6dQpqWNHqU4dado0s9MAAICsinKNHC8hQerQQcqdW1q2THJ3NzsRAADIqphzjRzNMKRBg6QffpCioqQiRcxOBAAAsjLKNXK08HBp9mxp7lypVi2z0wAAgKyOaSHIsSIjbcvuhYRIPXuanQYAAGQHlGvkSH/8IXXqJDVqJE2ebHYaAACQXVCukeNcuya1ayflzy8tXmy7kBEAAMAZqBXIUQxD6ttXOnpU2rlT8vMzOxEAAMhOKNfIUSZPlr74Qlq6VKpWzew0AAAgu2FaCHKMDRukUaOkMWOkwECz0wAAgOyIco0c4ehRqUsXqWVL6Y03zE4DAACyK8o1sr1Ll6S2baWiRaUFCyQrn3oAAJBBmHONbC05WXr2WenUKWnPHsnb2+xEAAAgO6NcI1t74w1p9Wrbo1Ils9MAAIDsjnKNbOvLL6Xx46UJE6TWrc1OAwAAcgJmnyJbOnjQNh0kMFAaPdrsNAAAIKegXCPbOX/edgfGcuWkzz+XLBazEwEAgJyCaSHIVpKSbEvuxcVJ0dFS3rxmJwIAADkJ5RrZyujR0qZNUkSE9PDDZqcBAAA5DeUa2caiRbbbm0+ZIv3nP2anAQAAORFzrpEt7Nsn9esn9eghhYSYnQYAAORUlGtkeWfP2i5grFJF+uQTLmAEAADmoVwjS7txQ+rUSUpIsK1r7eVldiIAAJCTMecaWdrw4dKOHdKWLVKpUmanAQAAOR3lGlnWZ59J4eHS9OlSgwZmpwEAAGBaCLKoXbukQYOkAQOk554zOw0AAIAN5RpZzunTUocOUq1atjPXXMAIAABcBeUaWcr161LHjrZCvXy55O5udiIAAID/x5xrZBmGIQ0ZYlvTets2qWhRsxMBAAD8X3t3Hh1Vff9//DUDIRkS9jURNASEFlzCEvbQQlXAUmQJYECFSEBrjbQuwVhQkS8CBsR6pBUQFUE2Qb9Uj1sRF5CYIISCGooiFMH4JcjSGBIC5PP7Y36ZMmYhwGTu3Mnzcc4cwmfuvfOa4XM+eXt9z73eKK5hG3/9q7R0qfTyy1L37lanAQAAKIu2ENjCxx9Lf/yjNGWKNH681WkAAADKR3GNgHfwoPtGMfHxUnq61WkAAAAqRnGNgHbqlPvW5uHh0tq1UkiI1YkAAAAqRs81ApYx7utY79kjZWRITZtanQgAAKByFNcIWPPnSytXSqtXS9dfb3UaAACAC6MtBAHpvfekqVOlhx+WxoyxOg0AAEDVUFwj4HzzjXTrrdLAgdL//I/VaQAAAKqO4hoBJT/f/QXGZs3cLSG1almdCAAAoOrouUbAKClxX8P64EEpM1Nq2NDqRAAAABeH4hoBY9Ys6Y03pA0bpF/+0uo0AAAAF4+2EASEv/9devRRacYMaehQq9MAAABcGoprWC4nR7rtNmn4cGnaNKvTAAAAXDqKa1jqxAnpllukK6+Uli2TnMxIAABgY/RcwzLnzkljx0p5edK2bVK9elYnAgAAuDwU17DMtGnum8W8847Urp3VaQAAAC4fxTUssXatNGeOlJ4u3XST1WkAAAB8gw5X+N0//yklJUmJidIDD1idBgAAwHcoruFXR4+6v8DYoYP0wguSw2F1IgAAAN+hLQR+c/asNHq0dOqU9MknUt26VicCAADwLYpr+M2DD0qbN0sbN7ovvQcAABBsKK7hF8uWSX/5i/Tcc9KvfmV1GgAAgHPZYGwAABv8SURBVOpBzzWqXVaWdNdd0sSJ0j33WJ0GAACg+lBco1r98IM0YoTUubO0cCFfYAQAAMHNtsX1Rx99JKfTWe4jKyvLa9ucnBwNGjRI9erVU5MmTXTHHXfo6NGjFiWvOU6flkaOlEpKpPXrpdBQqxMBAABUL9v3XE+ZMkVxcXFeY23btvX8fOjQIfXr10+NGjXS7NmzlZ+fr3nz5mn37t3KyspSSEiIvyPXGPfdJ33+ufTxx1JUlNVpAAAAqp/ti+v4+HiNGDGiwueffPJJFRYWKjs7W61atZIkde/eXTfeeKNefvllTZo0yV9Ra5Tnn5cWL5aWLpV69rQ6DQAAgH/Yti2klDFG+fn5Onv2bLnPr1+/XkOGDPEU1pL0m9/8Ru3bt9fatWv9FbNG2bxZSkmR7r1XuvNOq9MAAAD4j+2L66SkJDVo0EAul0sDBgzQ9u3bPc8dPnxYeXl56tatW5n94uLilJ2d7c+oNcJ330kJCVKfPtLTT1udBgAAwL9s2xYSGhqqhIQE3XzzzWratKm+/PJLzZs3T/Hx8dq6datiY2OVm5srSYqMjCyzf2RkpI4dO6YzZ87Qd+0jhYXS8OFSWJj02msSHysAAKhpbFtc9+rVS7169fL8fciQIUpISNB1112ntLQ0vfPOOyosLJTkLsR/LiwsTJJUWFhIce0DxkiTJ0tffSV9+qnUrJnViQAAAPzP9m0h52vbtq2GDh2qDz/8UMYYuVwuSdLp06fLbFtUVCRJnm1weZ55Rlqxwv0Fxs6drU4DAABgDdueua5I69atVVxcrIKCAk87SGl7yPlyc3PVpEmTSs9a/+lPf1KDBg28xhITE5WYmOjb0Da3caP04IPSQw9JfDQAAMDXVq1apVWrVnmNnTx50qI0lQu64vrbb7+Vy+VSRESEIiIi1KxZM23btq3MdllZWYqNja30WAsWLFCXLl2qK2pQ+PZbacwY6cYbpdmzrU4DAACCUXknN3fs2KGuXbtalKhitm0LycvLKzP2z3/+U3//+9910003ecZGjhypt956S4cOHfKMffDBB/r66681atQov2QNVj/9JA0bJjVqJK1aJdWqZXUiAAAAa9n2zPWYMWNUt25d9erVS82bN9dXX32lxYsXKyIiQnPmzPFs98gjj+i1115T//79NWXKFOXn5ys9PV3XXXedkpKSLHwH9maMlJQk7d8vffaZu8AGAACo6WxbXA8fPlyvvvqqFixYoP/85z9q3ry5EhIS9NhjjykmJsazXatWrfTxxx/r/vvv18MPP6zQ0FANGTJE8+fP5yohl+HJJ6V166Q33pA6dbI6DQAAQGCwbXGdkpKilJSUKm3bsWNHvfvuu9WcqOZ46y1p+nTpscfcbSEAAABws23PNayxZ480bpw0dKj06KNWpwEAAAgsFNeospMn3Weqr7hCeuUVycnsAQAA8GLbthD4V0mJ+4z1Dz9I27ZJ9etbnQgAACDwUFyjSh59VHr7bffj6qutTgMAABCYKK5xQevWSbNmSXPmSIMGWZ0GAAAgcNE1i0rt2iWNH+++C2NqqtVpAAAAAhvFNSr044/uLzBefbW0dKnkcFidCAAAILDRFoJynT3rPludny9t2iSFh1udCAAAIPBRXKNcU6dKH30k/eMfUnS01WkAAADsgeIaZSxfLj39tPSXv0j9+1udBgAAwD7ouYaXzz+XJk2SJkyQqnh3eQAAAPx/FNfw+L//k4YPl66/Xvrb3/gCIwAAwMWiuIYkqbhYSkhwf5Hx9delsDCrEwEAANgPPdeQJE2ZImVmur/EeMUVVqcBAACwJ4praPFi6fnnpSVLpN69rU4DAABgX7SF1HCffirde6/0+99LyclWpwEAALA3iusa7NAhaeRIqWdP6ZlnrE4DAABgfxTXNVRRkTRihBQSIq1bJ9WpY3UiAAAA+6PnugYyRrr7bmn3bmnLFql5c6sTAQAABAeK6xro2WelZcukFSukrl2tTgMAABA8aAupYTZtkh54wP0YN87qNAAAAMGF4roG2b9fGj1aGjBAmjPH6jQAAADBh+K6higokIYNkxo0kFavlmrTEAQAAOBzlFg1gDHSnXdK+/ZJGRlS48ZWJwIAAAhOFNc1wNy50tq17kvuXXut1WkAAACCF20hQe7tt6VHHpGmTXPfMAYAAADVh+I6iO3dK40dK/32t9KMGVanAQAACH4U10HqP/9xf4GxZUv39ayd/EsDAABUO3qug1BJiXTbbdLhw1JWlvsKIQAAAKh+FNdBaMYM6a23pDfflDp0sDoNAABAzUFxHWRef1164glp1ix3rzUAAAD8h07cIPLFF9Idd0ijRklpaVanAQAAqHkoroPEsWPSLbdIMTHSSy9JDofViQAAAGoe2kKCwNmzUmKidOKEtHGjFB5udSIAAICaieI6CKSlSR98IL33ntSmjdVpAAAAai6Ka5tbuVKaN09asED6zW+sTgMAAFCz0XNtYzt2SBMnSrffLk2ZYnUaAAAAUFzb1JEj7jswXnONtGgRX2AEAAAIBBTXNnTmjPtye6dPS2+8IblcVicCAACARM+1Lf3pT1JGhrRpk9SqldVpAAAAUIri2maWLpUWLpSef17q29fqNAAAADgfbSE2kpEh3XOPdNdd7gcAAAACC8W1TXz/vTRypBQXJz37rNVpAAAAUB6KaxsoKpJGjJCcTmndOqlOHasTAQAAoDz0XAc4Y6Q//EHauVPavFlq2dLqRAAAAKgIxXWAW7hQevFFadkyd0sIAAAAAhdtIQHso4+kP/7R/bjjDqvTAAAA4EIorgPUv//tvlHMr34lpadbnQYAAABVQXEdgE6dkoYPlyIipDVrpNo07wAAANgCZVuAMUaaOFH617/c17Vu2tTqRAAAAKgqiusAM2+etHq1tHatdN11VqcBAADAxaAtJIC895708MNSWpq73xoAAAD2QnEdIL75Rrr1VmnQIGnmTKvTAAAA4FJQXAeA/Hzplluk5s2lV1+VatWyOhEAAAAuBT3XFispcV/D+rvvpMxMqWFDqxMBAADgUlFcW2zmTGnDBvfjl7+0Og0AAAAuB8W1hTZskB5/3F1g/+53VqcBAADA5aLn2iJffSXddps0YoT0yCNWpwEAAIAvUFxb4Phx9xcYo6OlZcskJ/8KAAAAQYG2ED87d04aO1b68Udp2zb3Lc4BAAAQHCiu/ezPf5bef196912pbVur0wAAAMCXKK79aPVqae5c9y3Ob7zR6jQAAADwNbp9/WTnTunOO6Vx46T777c6DQAAAKoDxbUf5OVJw4a5r2O9ZInkcFidCAAAANWBtpBqduaMNHq0dOqU9MknkstldSIAAABUF4rravbgg9KWLdIHH0hXXml1GgAAAFQniutq9NJL0rPPSgsXSv36WZ0GAAAA1Y2e62qSmSndfbeUnCz9/vdWpwEAAIA/UFxXg9xc923Nu3aVnnuOLzACAADUFBTXPnb6tDRypPvn9eul0FBr8wAAAMB/6Ln2IWOke++Vtm93XxkkMtLqRAAAAPAnimsfev556YUXpBdflHr0sDoNAAAA/I22EB/55BPpvvuklBQpKcnqNAAAALACxbUPHDwoJSRI8fHS/PlWpwEAAIBVKK4vU2GhNHy4VLeutGaNFBJidSIAAABYhZ7ry2CMNGmSlJMjbd0qNWtmdSIAAABYieL6Mjz9tPTqq9Lq1VJsrNVpAAAAYDXaQi7RP/4hpaZKU6dKY8ZYnQYAAACBgOL6Euzb5y6ob7pJmjXL6jQAAAAIFBTXF+mnn6Rhw6QmTaSVK6VataxOBAAAgEBBz/VFMEYaP146cEDKzJQaNbI6EQAAAAIJxfVFmDVLev116X//V+rY0eo0AAAACDQ1oi3k9OnTmjp1qqKiolS3bl317NlTGzduvKhjvPmmNH269Pjj0i23VE9OAAAA2FuNKK4nTJigBQsW6Pbbb9ezzz6rWrVq6eabb9ann35apf337JHGjXP3Wk+fXs1hAQAAYFtBX1xnZWVpzZo1mjNnjubOnavk5GRt2rRJV111lVJTUy+4/4kT7jPVrVtLr7wiOYP+E4O/rVq1yuoIqCGYa/AX5hpqsqAvFdetW6fatWtr8uTJnrHQ0FBNnDhRGRkZOnz4cIX7njvnPmN95Ii0YYNUr54/EqOm4ZcQ/IW5Bn9hrqEmC/riOjs7W+3bt1dERITXeFxcnCRp586dFe77/PPSu++678DYrl21xgQAAEAQCPqrheTm5ioyMrLMeOnY999/X+G+L74ozZ0rDRxYbfEAAAAQRIL+zHVhYaFCQ0PLjIeFhXmer8jAgdJDD1VbNAAAAASZoD9z7XK5dPr06TLjRUVFnud/rrTgHj06R9nZ1ZsPOHnypHbs2GF1DNQAzDX4C3MN/pCTkyOp8hOlVgj64joyMrLc1o/c3FxJUlRUVJnnDhw4IEmaOPG2as0GlOratavVEVBDMNfgL8w1+MuBAwfUp08fq2N4BH1x3blzZ3300UfKz89XvfMu95GZmSlJio2NLbPPwIEDtWLFCkVHR5d7ZhsAAADWKiws1IEDBzQwwL4c5zDGGKtDVKesrCz17NlT6enpeuCBByS579h4zTXXqFmzZtq6davFCQEAABAsgv7Mdffu3TVq1CilpaXpyJEjatu2rZYtW6aDBw/qpZdesjoeAAAAgkjQn7mW3Geqp0+frhUrVuj48eO6/vrrNXPmTN14441WRwMAAEAQqRHFNQAAAOAPQX+dawAAAMBfbFVcnz59WlOnTlVUVJTq1q2rnj17auPGjVXa98SJE5o8ebKaNWumiIgIDRgwQNkVXMR669at6tu3r8LDwxUZGakpU6aooKCgzHbGGD311FNq06aNXC6Xrr/+eq1evbrcYx4+fFijR49Wo0aN1KBBAw0bNkz79++v+puHX9l1rj3++ONyOp1lHlz1JnAF2lybNWuWhg4dqhYtWsjpdGrGjBkVvj7rmr3Yda6xrtlPIM21PXv2KDU1VbGxsapfv76ioqI0ZMgQbd++vdxj+mRdMzZy6623mpCQEJOammqWLFlievfubUJCQsyWLVsq3e/cuXOmd+/eJiIiwjzxxBNm4cKFplOnTqZ+/frm66+/9to2OzvbhIWFma5du5pFixaZadOmmbCwMDN48OAyx3344YeNw+Ewd911l3nhhRfMkCFDjMPhMKtXr/baLj8/31x99dWmZcuWJj093SxYsMBceeWVpnXr1ubHH3+8/A8GPmfXufbYY48Zh8NhFi1aZF599VXP4+fbIXAE2lxzOBwmKirKDBo0yDgcDjNjxoxyX591zX7sOtdY1+wnkObaAw88YBo1amQmTZpklixZYtLT0027du1M7dq1zcaNG7229dW6ZpviOjMz0zgcDjN//nzPWFFRkWnXrp3p3bt3pfuuWbPGOBwOs379es9YXl6eadSokRk7dqzXtoMHDzZXXHGFyc/P94y98MILxuFwmPfff98zdujQIRMSEmJSUlK89u/Xr59p3bq1OXfunGds7ty5xuFwmM8//9wztmfPHlO7dm3zyCOPVPETgL/Yea6V/hKiuLGHQJtrxhjz73//2xhjzNGjRysteFjX7MXOc411zV4Cba5t377dFBQUeO37448/mubNm5u+fft6jftqXbNNcf3QQw+ZkJAQrw/RGGNmz55tHA6HOXToUIX7jho1ykRGRpYZv+uuu0x4eLgpLi42xhhz8uRJExISYqZOneq1XXFxsalXr55JTk72jC1cuNA4HA6Tk5Pjte2qVauMw+Hw+q+zuLg406NHjzKvP3DgQNOuXbtK3jWsYOe5VvpLKC8vz5w8edKUlJRU/Y3D7wJtrp0vLy+v0oKHdc1e7DzXWNfsJZDn2vlGjBhhmjRp4jXmq3XNNj3X2dnZat++vSIiIrzG4+LiJEk7d+6sdN8uXbqUGY+Li9OpU6e0d+9eSdLu3bt19uxZdevWzWu7kJAQxcbGevX8ZGdnKyIiQr/4xS8qzVNSUqJdu3aVOWbptvv27Su3Fw3WsetcO19MTIwaNmyo+vXr6/bbb9eRI0cqe8uwSKDNtapiXbMfu86187Gu2YNd5toPP/ygZs2aef7uy3XNNsV1bm6uIiMjy4yXjn3//feXvW9ubq7X+Platmzp9Rq5ublq0aLFBY957NgxFRcXX3J2+J9d55okNW7cWCkpKVq8eLHWr1+v5ORkrVmzRvHx8crPz68wN6wRaHOtqljX7Meuc01iXbMbO8y1zZs367PPPtOYMWM8Y75c12xzh8bCwkKFhoaWGQ8LC/M8X5GioqIq7Vv6Z0Xbnv8aVc1zoWNeKDv8z65zTZLuu+8+r22GDx+u7t27a9y4cfrrX/+qqVOnVpgd/hdoc+1icld2zPO3QWCw61yTWNfsJtDn2pEjRzR27FjFxMQoNTXVK3dlx7xQ9vPZ5sy1y+XS6dOny4wXFRV5nr/cfUv/rGjbunXreh2zdP/LOeaFssP/7DrXKpKYmKiWLVvqgw8+qHQ7+F+gzbWLyV3ZMc/fBoHBrnOtIqxrgSuQ51pBQYGGDBmigoICbdiwoczv2sqOeaHs57NNcR0ZGVnu6fjS/zUQFRV12fuWnvYvHf/5tue/RmRkpH744YcLHrNx48YKDQ2t8JgXyg7/s+tcq0yrVq107NixC24H/wq0uVZVrGv2Y9e5VhnWtcAUqHOtuLhYI0aM0BdffKENGzaoY8eOXs/7cl2zTXHduXNn7d27t0x/VWZmpiQpNja2wn1jY2O1Y8cOmZ/d6T0zM1Ph4eFq3769JOmaa65R7dq1tW3bNq/tiouLtXPnTq/X6Ny5s06dOqWcnJxK8zidTl177bVljlm6bdu2bRUeHl7pe4d/2XWuVcQYowMHDnh9cQOBIdDmWlWxrtmPXedaRVjXAlcgzrWSkhLdcccd+vDDD7Vy5UrFx8eXeW2frmtVvq6IxUqvmzhv3jzPWOl1E3v16uUZy83NNTk5OebMmTOesdLrJq5bt84zlpeXZxo2bGgSExO9Xmfw4MEmKiqq3Osmvvfee56xQ4cOmTp16ph7773XM1ZSUmLi4+NN69atvS4VVNl1E9PS0i71I0E1sfNcO3LkSJn3U3opv2eeeeZiPwpUs0Cba+e70OXRWNfsxc5zjXXNXgJxrt1zzz3G4XCYJUuWVJrdV+uabYprY4wZPXq0544/ixYtMr179zZ16tQxmzdv9mwzfvx443A4PBenN8Z9x59evXqZevXqed3xp0GDBmbv3r1er7Fjxw4TFhZmunTpYv72t7+ZP//5z8blcplBgwaVyZOamuq5a96SJUvMb3/7W+NwOMyqVau8tsvPzzft2rUzLVq08Nzxp3Xr1qZVq1bm6NGjPv6U4At2nWsul8skJSWZ+fPnm4ULF5rExETjdDpNly5dTGFhoY8/JfhCoM21V155xcycOdOkpaUZh8NhBgwYYGbOnGlmzpzp9fqsa/Zj17nGumY/gTTXFixYYBwOh+ndu7dZsWKFWb58udfj/BvM+Gpds1VxXVRUZB566CETGRlpwsLCTI8ePcrc8WnChAnG6XR6/WMZY8zx48dNcnKyadq0qQkPDzf9+/c327dvL/d1tmzZYvr06WNcLpdp0aKFSUlJMT/99FOZ7UpKSszs2bNNdHS0CQ0NNddee61ZuXJlucc8dOiQGTVqlGnQoIGpV6+eGTp0qNm3b98lfhKobnada5MmTfLcKrZOnTqmffv2Ji0trdxjIjAE2lz79a9/bRwOh3E4HMbpdBqn0+n5+eOPP/balnXNXuw611jX7CeQ5lrp65TOr/Mf5b2+L9Y1hzE/a2wBAAAAcEls84VGAAAAINBRXAMAAAA+QnENAAAA+AjFNQAAAOAjFNcAAACAj1BcAwAAAD5CcQ0AAAD4CMU1AAAA4CMU1wAAAICPUFwDALy8/PLLcjqdOnjwoNVRAMB2KK4BIMDs3r1bCQkJio6OlsvlUqtWrXTTTTfpueee82zz5JNPasOGDdWWweFwVNuxASCYOYwxxuoQAAC3rVu3qn///oqOjtb48ePVsmVLfffdd8rIyNC3336rvXv3SpIiIiI0evRovfjiiz7PUFJSorNnz6pOnTo+PzYABLvaVgcAAPzXrFmz1KhRI23btk3169f3eu7o0aOenx0Oh6p6bqSgoEDh4eFVzuB0OimsAeAS0RYCAAFk37596tSpU5nCWpKaNm0qyV38FhQUaNmyZXI6nXI6nbrzzjslSY8//ricTqdycnI0duxYNW7cWPHx8ZKkXbt2acKECYqJiZHL5VJkZKQmTpyoY8eOeb1OeT3X0dHR+t3vfqctW7aoe/fucrlcatu2rZYvX15dHwUA2BJnrgEggERHRysjI0NffvmlOnXqVO42y5cvV3Jysnr06KHJkydLktq2beu1zahRo9S+fXvNnj3bc4Z748aN2r9/vyZOnKiWLVvqiy++0OLFi/Xll1/qs88+qzSXw+HQN998o1GjRik5OVlJSUlaunSpJkyYoK5du6pjx44+ePcAYH8U1wAQQB588EENHjxYsbGxiouLU3x8vG644Qb1799ftWu7l+xx48bp7rvvVkxMjMaOHVvucWJjY7VixQqvsXvuuUf333+/11jPnj2VmJioLVu2qG/fvhXmMsboX//6lzZv3qw+ffpIchfwrVu31ksvvaT09PTLedsAEDRoCwGAAHLDDTcoIyNDQ4cO1a5du5Senq6BAwfqiiuu0Jtvvlnl49x9991lxsLCwjw/FxUV6ejRo+rRo4ckKTs7+4LH7NSpk6ewltxtKh06dND+/furnAsAgh3FNQAEmG7dumn9+vU6ceKEsrKylJaWpvz8fCUkJCgnJ6dKx2jTpk2ZsWPHjmnKlClq0aKF6tatq+bNmysmJkaSdPLkyQse88orrywz1rBhQx0/frxKmQCgJqAtBAACVO3atdWtWzd169ZN7du3V1JSktatW6fp06df8EohLperzNjo0aOVkZGh1NRUxcbGKiIiQufOndOgQYNUUlJywTy1atUqd5wrugLAf1FcA4ANdO3aVZKUm5sr6eJv8nL8+HFt2rRJTzzxhKZNm+YZ//rrr30XEgBAWwgABJIPP/yw3PG3335bktShQwdJUnh4+EW1Y5Sedf75GepnnnnmUmICACrAmWsACCApKSkqLCzU8OHD1aFDBxUXF2vr1q1au3at2rRpo6SkJEnuvuyNGzdqwYIFioyMVExMjLp3717hcevXr69+/frpqaee0pkzZxQVFaX3339fBw4cuOzMtIUAwH9RXANAAJk/f75ee+01vf3221q8eLGKi4t11VVX6Q9/+IOmTZvmubnM008/rcmTJ2vatGkqLCzUhAkTPMV1RS0jK1euVEpKihYuXChjjAYOHKh33nlHUVFRZbb9+TEqOqbD4bjoFhUACGYOwykHAAAAwCfouQYAAAB8hOIaAAAA8BGKawAAAMBHKK4BAAAAH6G4BgAAAHyE4hoAAADwEYprAAAAwEcorgEAAAAfobgGAAAAfITiGgAAAPARimsAAADARyiuAQAAAB+huAYAAAB85P8BW8LXIVlBDbEAAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f31a0ba76d0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(0,250)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ϵ_last = zeros(Float64, (6))\n",
    "ϵᵖ = zeros(Float64, (6))\n",
    "σ = zeros(Float64, (6, 1))\n",
    "Δt = 1.0\n",
    "σy =  200.0\n",
    "ss = zeros(Float64, steps)\n",
    "ee = zeros(Float64, steps)\n",
    "for i=1:steps\n",
    "    ss[i] = σ[1]\n",
    "    ee[i] = ϵ_last[1]\n",
    "    dϵ = (reshape(ϵ_tot[i, :, :], (6, 1)) - ϵ_last) / Δt \n",
    "    dϵᵖ, σ = calculate_stress(dϵ, Δt, σ, C, σy)\n",
    "    ϵ_last += dϵ * Δt\n",
    "end\n",
    "PyPlot.plot(ee, ss)\n",
    "PyPlot.title(\"Stress-Strain curve\")\n",
    "PyPlot.xlabel(\"Strain\")\n",
    "PyPlot.ylabel(\"Stress\")\n",
    "PyPlot.ylim([0, 250])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.0-dev",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}