diff --git a/docs/tutorials/2015-11-03-constitutive-modelling-using-juliafem.ipynb b/docs/tutorials/2015-11-03-constitutive-modelling-using-juliafem.ipynb new file mode 100644 index 0000000..45ebd40 --- /dev/null +++ b/docs/tutorials/2015-11-03-constitutive-modelling-using-juliafem.ipynb @@ -0,0 +1,342 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Constitutive modelling using JuliaFEM\n", + "\n", + "Author(s): Jukka Aho\n", + "\n", + "**Abstract**: This is reproduction of the results from notebook [*Ideal plastic Von Mises material*](https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-09-24-Ideal%20plastic%20Von%20Mises%20material.ipynb) made by Olli Väinölä. Small strain theory is used, see https://en.wikipedia.org/wiki/Flow_plasticity_theory. The purpose of this notebook is to show how one can easily design and simulate material models using JuliaFEM." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "using JuliaFEM: IntegrationPoint, Field, FieldSet, TimeStep, Increment\n", + "using ForwardDiff\n", + "using PyPlot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The atomic structure here is `IntegrationPoint`. It has identical `Field`-structure like elements and can store multidimensional variables which can be time-dependent also. That way one can easily store, for example, measured strain and easily run material simulation for measured data." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[])" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ip = IntegrationPoint([0.0, 0.0], 1.0)\n", + "\n", + "ip.fields[\"total strain\"] = Field()\n", + "ip.fields[\"plastic strain\"] = Field()\n", + "ip.fields[\"plastic potential\"] = Field()\n", + "ip.fields[\"elastic strain\"] = Field()\n", + "ip.fields[\"effective plastic strain\"] = Field()\n", + "ip.fields[\"stress\"] = Field()\n", + "ip.fields[\"young\"] = Field(200.0e9)\n", + "ip.fields[\"poisson\"] = Field(0.3)\n", + "ip.fields[\"yield stress\"] = Field(200.0e6)\n", + "ip.fields[\"plastic rate parameter\"] = Field() # material parameter" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Accessing parameters is done just like with `FieldSet`s in general:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: [2.0e11]\n", + "INFO: [2.0e8]\n" + ] + } + ], + "source": [ + "info(last(ip.fields[\"young\"]))\n", + "info(last(ip.fields[\"yield stress\"]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Our material model is" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "calculate_stress! (generic function with 1 method)" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\"\"\" Ideal plastic material model. \"\"\"\n", + "function calculate_stress!(ip::IntegrationPoint, strain::Matrix, time::Number)\n", + "\n", + " poisson = last(ip.fields[\"poisson\"])[1]\n", + " young = last(ip.fields[\"young\"])[1]\n", + " mu = young/(2*(1+poisson))\n", + " lambda = young*poisson/((1+poisson)*(1-2*poisson))\n", + "\n", + " # yield function\n", + " function f(stress)\n", + " stress_dev = stress - 1/3*trace(stress)*I # stress deviator tensor\n", + " stress_v = sqrt(3/2*trace(stress_dev*stress_dev')) # equivalent tensile stress\n", + " stress_y = last(ip.fields[\"yield stress\"])[1]\n", + " return stress_v - stress_y\n", + " end\n", + "\n", + " stress_trial = lambda*trace(strain)*I + 2*mu*strain\n", + "\n", + " if f(stress_trial) <= 0.0\n", + " info(\"time=$time: no yield\")\n", + " push!(ip.fields[\"stress\"], TimeStep(time, Increment(Matrix[stress_trial])))\n", + " push!(ip.fields[\"total strain\"], TimeStep(time, Increment(Matrix[strain])))\n", + " return true\n", + " else\n", + " info(\"time=$time: yield, f(stress_trial) = $(f(stress_trial))\")\n", + " dt = time - ip.fields[\"total strain\"][end].time\n", + "\n", + " # associated flow rule, plastic potential ψ(σ) = f\n", + " psi = f\n", + "\n", + " \"\"\" Calculate equations\n", + "\n", + " dσ - C (dϵₜ - dγ*dΨ/dσ) = 0\n", + " σₑ(σ) - σy = 0\n", + "\n", + " dϵₜ = total strain\n", + " dϵₑ = elastic strain\n", + " dϵₚ = plastic strain\n", + " \n", + " dϵₜ = dϵₑ + dϵₚ\n", + " => dϵₑ = dϵₜ - dϵₚ = dϵₜ - dγ*dψ/dσ\n", + " dσ = λ⋅tr(dϵₑ)I + 2μ⋅dϵₑ\n", + " => dσ - λ⋅tr(dϵₑ)I - 2μ⋅dϵₑ = 0\n", + " \"\"\"\n", + " function residual(params::Vector)\n", + " stress = reshape(params[1:9], 3, 3)\n", + " gamma = params[end]\n", + "\n", + " strain_prev = last(ip.fields[\"total strain\"])[1]\n", + " stress_prev = last(ip.fields[\"stress\"])[1]\n", + "\n", + " dstrain_total = 1/dt*(strain - strain_prev)\n", + " #dstress = 1/dt*(stress - stress_prev)\n", + " dstress = 1/dt*(stress_trial - stress)\n", + " # derivative of plastic potential ψ(σ) w.r.t 2nd order tensor σ(ϵ)\n", + " # https://en.wikipedia.org/wiki/Tensor_derivative_%28continuum_mechanics%29\n", + " # \"Derivatives of scalar valued functions of second-order tensors\"\n", + " dpsi_dstress = derivative(psi, stress)\n", + " dstrain_plastic = gamma*dpsi_dstress\n", + " dstrain_elastic = dstrain_total - dstrain_plastic\n", + " dstress_elastic = lambda*trace(dstrain_elastic)*I + 2*mu*dstrain_elastic\n", + "\n", + " stress_delta = dstress - dstress_elastic\n", + " return [vec(stress_delta); psi(stress)]\n", + " end\n", + "\n", + " \n", + " # solve equations using Newton iterations. Jacobian is calcualated\n", + " # using automatic differentiation as usual.\n", + " params = [vec(stress_trial); 0.0]\n", + " for i=1:5\n", + " A = ForwardDiff.jacobian(residual, params)\n", + " b = -residual(params)\n", + " dparams = A \\ b\n", + " params += dparams\n", + " norm(dparams) < 1.0e-7 && break\n", + " end\n", + " stress = reshape(params[1:9], 3, 3)\n", + " push!(ip.fields[\"stress\"], TimeStep(time, Increment(Matrix[stress])))\n", + " push!(ip.fields[\"total strain\"], TimeStep(time, Increment(Matrix[strain])))\n", + " end\n", + "\n", + "end ### END OF CONSTITUTIVE MODEL" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Small \"simulator\" to study the behavior of material model in `IntegrationPoint`" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: time=1.0: no yield\n", + "INFO: time=2.0: no yield\n", + "INFO: time=3.0: yield, f(stress_trial) = 2.0e8\n" + ] + } + ], + "source": [ + "function run(steps=3)\n", + "\n", + " # initialization\n", + " ip = IntegrationPoint([0.0, 0.0], 1.0)\n", + " ip.fields[\"total strain\"] = Field()\n", + " ip.fields[\"plastic strain\"] = Field()\n", + " ip.fields[\"plastic potential\"] = Field()\n", + " ip.fields[\"elastic strain\"] = Field()\n", + " ip.fields[\"effective plastic strain\"] = Field()\n", + " ip.fields[\"stress\"] = Field()\n", + " ip.fields[\"young\"] = Field(200.0e9)\n", + " ip.fields[\"poisson\"] = Field(0.3)\n", + " ip.fields[\"yield stress\"] = Field(200.0e6)\n", + " ip.fields[\"plastic rate parameter\"] = Field() # material parameter\n", + "\n", + " # run\n", + " for (time, strain_11) in enumerate(linspace(0, 200e-5, steps))\n", + " # set up strain etc.\n", + " time = Float64(time)\n", + " poisson = last(ip.fields[\"poisson\"])[1]\n", + " strain = zeros(3, 3)\n", + " strain[1, 1] = strain_11\n", + " strain[2, 2] = -poisson*strain_11\n", + " strain[3, 3] = -poisson*strain_11\n", + "\n", + " # calculate stress\n", + " calculate_stress!(ip, strain, time)\n", + " end\n", + "\n", + " return ip\n", + "end\n", + "\n", + "ip = run();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualize results:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: [0.0 200.00000000000003 200.00000000000006]\n", + "INFO: [0.0 1000.0 2000.0]\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoIAAAHqCAYAAACQvy+/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlcTfnjBvDn3rQJqYiyZewMk+wMBoPoRiUlsi8zQtZhMDOWyD62RmMPCcnSYsmWhGSLMZYZjGWQJZGtTff8/vDt/uZON9O9qnO793m/Xl6v6ZzPPee593u/eZzlcySCIAggIiIiIr0jFTsAEREREYmDRZCIiIhIT7EIEhEREekpFkEiIiIiPcUiSERERKSnWASJiIiI9BSLIBEREZGeYhEkIiIi0lMsgkRERER6ikWQiEjL2NnZoXr16mLHICI9wCJIpMOys7Oxdu1atG/fHpaWljAyMkKFChXwxRdfYPjw4YiMjFQaHxQUBKlUik2bNomUuPClp6dj8eLFaNGiBczNzWFsbAxbW1s0bdoUY8aMwYkTJ5TGz5w5E1KpNNfywiSRSCCRSIpsf0Skv0qIHYCICkd2djZkMhmio6NhYWEBmUyGypUrIzMzE7///jtCQkLwxx9/wNnZOddrdbWEvHnzBu3bt0diYiJsbGzQu3dvVKxYEW/evMGlS5ewZs0apKamol27dqLmPHbsmKj7JyL9wSJIpKO2bduG6Oho2NvbIzY2FqVLl1Zan5aWhrNnz6p8rSAIRRGxyC1btgyJiYno2rUrIiMjUaKE8q/Aly9f4saNGypfW5SfCU8LE1FR4alhIh11+vRpAMCgQYNylUAAMDU1Rfv27RU/f/XVVxgyZAgAYPDgwZBKpYo/9+/fB/D/p0ljY2MREhKCFi1aoFSpUkrF5d27d5g3bx7s7e1RqlQplC5dGq1bt8b27dtV5ty0aRNat26N8uXLw9TUFFWrVoWjoyNCQ0OVxv3222/w8vKCnZ0dTExMYG1tjSZNmmD8+PF4//69Wp/JyJEjc5VAAChbtixatmyp+NnOzg6zZ88GAHTo0EHpM8kxaNAgSKVS3LlzBytXrkSjRo1QsmRJdOjQAQCQlZWFgIAAdO/eHdWqVYOJiQmsrKzQuXNnHDx4UGVOVdcI/vO0fUxMDL766iuUKVMG5ubmkMlkeRbYjzl06BCcnZ1hbW0NExMTVK1aFS4uLjh69KjK/aoilUoV7zXHx74nCQkJkEqlcHNzyzNXvXr1YGJigpcvXyotj46ORvfu3VGuXDmYmJigZs2amDx5MlJTU9V+70T0AY8IEumocuXKAQD++OOPfI0fPHgwLCwsEB4eDhcXF9jb2yvWmZubK41dsmQJDh8+jB49eqBTp06Kv4hfvnyJjh074tKlS2jSpAmGDh0KuVyOgwcPom/fvrh69Sr8/PwU25k2bRrmz5+Pzz77DH369IG5uTkePXqEc+fOISwsDB4eHgA+lMAWLVrAwMAAPXr0QPXq1fHq1SvcvHkTgYGBmDt3rspi96mfyfjx47F3717ExsZi0KBBsLOzy3Ps2LFjERcXB5lMBplMBgMDAwDA8+fPMW7cOLRp0wZdu3ZF+fLl8ejRI0RGRqJ79+5Yu3Ythg4dmmt7eZ2ej4qKQnh4OLp3746RI0fi6tWr2L9/P86dO4dr167BysoqX+9txowZ8PPzQ+nSpeHi4oIqVarg4cOHOH36NLZu3YpOnTrlK8/H1qn6nrRo0QJ16tTB/v37kZKSAktLS6XXnD17Fn/88Qfc3d1RtmxZxfJZs2Zh1qxZsLKyUpTXy5cvY/Hixdi/fz/i4+NV/oOHiP6DQEQ6KTExUTAyMhKkUqnQv39/Yffu3cLdu3c/+pqNGzcKEolE2LRpk8r1M2bMECQSiVCqVCnh0qVLudYPHDhQkEgkwqJFi5SWp6enC46OjoJUKlV6naWlpVClShUhLS0t17aSk5MV/z1hwgRBIpEIERERuca9fPlSkMvlH31fOaKiogSJRCIYGxsLPj4+wr59+4RHjx599DU57zk2Nlbl+pz3XLlyZZWfb0ZGhvDw4cNcy1NTU4XPP/9csLS0zPX+q1WrJlSvXl1pWc7/NoaGhsKxY8eU1k2dOlWQSCTCwoULP/peckRHRwsSiUSoUaOGyvf/4MGDXPvN6zshkUiEDh06KC37r+/JvHnzBIlEIgQEBORa5+PjI0gkEiEqKkqx7NixY4JEIhHatGkjpKamKo0PCgoSJBKJMH78+I+/aSJSiaeGiXSUvb09goODUaFCBQQHB6NXr16oXr06rKys4ObmhqioKI23PWLECHzxxRdKy54/f47g4GA0a9YMkyZNUlpnbGyM+fPnQxAEhISEKJZLJBIYGhoqnWrNoerIlomJSa5l5ubm+b65xcnJCcuXL4epqSkCAwMhk8lQqVIl2NjYwNvbG3FxcfnajiqTJ09GtWrVci03MjKCra1truVlypTB4MGD8eLFC5w7dy7f++nTp0+uU7EjRowAgHxvZ+XKlQA+HLGzsbHJtb5SpUr5zvMxqr4nANC/f3+Vp5szMzOxfft2VKhQAd26dVMsX7FiBQBg7dq1KFOmjNJrBg4ciC+++AJbt24tkMxE+oanhol0WO/eveHq6oqYmBicOnUKiYmJOHnyJPbu3Yu9e/diwIABCAoKUnu7zZs3z7Xs3LlzkMvlAD5cI/ZvWVlZAIDr168rlvXr1w8rV65E/fr14eHhgfbt26Nly5a5TkX36dMHK1asgIuLC9zd3dGpUye0adMGNWrUUBp36dIl7N27V2mZhYUFxo4dq/h5zJgxGDZsGA4fPoz4+HgkJibi9OnTCAkJQUhICH788UfMmjVLvQ8Eqj+THFevXsWiRYtw4sQJPH78GOnp6UrrHz16lO/9NG3aNNeyypUrAwBevHiRr22cOXMGUqkUjo6O+d6vJvL6TCpVqoROnTrh8OHDuH79OurVqwcAiIyMxIsXLzBhwgSlfxzEx8fD0NAQoaGhKm/ayczMxLNnz/DixQtYWFgUzpsh0lEsgkQ6rkSJEujcuTM6d+4MAJDL5di1axeGDBmCzZs3w9XVFT179lRrmxUrVsy17Pnz5wA+FMK8jkxJJBK8fftW8fPSpUvx2WefYePGjZg/fz7mz5+PEiVKoHv37liyZImi6DVr1gxxcXGYO3cuwsLCsGXLFgBAnTp1MGPGDPTp0wcAcPnyZcXNHTns7OyUiiDw4UaZHj16oEePHgA+lNS1a9di7Nix8PPzg5ubm8ojWep+JsCH0tWxY0fI5XJ06tQJLi4uKFOmDKRSKRITExEeHo6MjIx87+ef183lyLk+Mjs7O1/bePnyJSwsLGBsbJzv/Woir88E+HCTzeHDh7Fp0ybMnz8fABRHCAcOHKg09vnz58jOzv5oQZdIJHjz5g2LIJGaeGqYSM9IpVL07t0b48ePBwDExMSovQ1Vp2JzjuJNmDABcrlc5Z/s7GylO1KlUinGjh2LS5cu4cmTJ9i1axdcXV0REREBR0dHZGZmKsa2bNkSkZGRePnyJU6dOoUff/wRT548Qd++fRXbHDhwYK59/vXXX//5fgwNDeHj4wMvLy8Ams3jl9fp6Tlz5iA9PR2HDh3Cvn378PPPP2PmzJn46aefPnoUsTCVLVsWL168yHVkUpWcI3Oq7sz+9129//axU/aurq4oU6YMgoODIQgCnj59igMHDsDe3h4NGzZUGmtubg5LS8s8v1c5360qVar85/shImUsgkR6qlSpUgCU58fLudM1v0eW/qlFixaf9ASO8uXLw9XVFTt27ECHDh1w+/ZtXL16Ndc4Q0NDtGrVCrNmzVJcOxYREaHRPv8t5zP5p0/5TADg1q1bsLKyUjlJdWxsrEbb/FStWrVS3M39X3KOsOVMIfRP58+f1ziDiYkJPDw88OjRIxw+fBghISHIzs7OdTQwJ29KSgquXbum8f6ISDUWQSIdtW3bNhw5ckTlNVWPHz/G2rVrAUCpoOTcoHHv3j2191e+fHn069cP58+fx5w5cxTXC/7T7du3cffuXQAfrus6depUrjFZWVlISUmBRCJByZIlAXyY/0/V0avHjx8DgGLcf/n111+RkJCgct2NGzewc+dOSCSSAvtMgA+TQz9//hxXrlxRWr5+/XocOnRIo21+qjFjxgAAJk6cqPL6xH8ua9asGaRSKUJCQpCWlqZYnpKSgsmTJ39SjkGDBgEANm/ejM2bN8PQ0BD9+vXLNS7n6PXw4cORlJSUa/3bt2/z/N+ViD6O1wgS6aizZ89i+fLlqFixIr788kvFHHh37tzBvn37kJ6eDhcXF/Tq1UvxmtatW6NkyZJYtmwZnj9/jgoVKgAAfH19c92tqUpAQABu3ryJn376CVu2bEGbNm1QoUIFPHr0CNevX8f58+exfft22NnZ4d27d2jbti1q1qwJBwcHVKtWDenp6Th8+DBu3LiBnj17ok6dOgCAhQsXIiYmBm3btoWdnR1KlSqFq1ev4uDBg7C0tFTcNftfoqOj4ePjAzs7O7Rp0waVK1dGRkYGbt68iejoaGRnZ8PX1xdNmjRRvKZjx46QSqWYOnUqrly5ojhC9sMPP+Rrn+PGjUN0dDS+/PJLeHh4oEyZMjh//jxOnToFd3d3hIWF5Ws7Balz58744YcfMGfOHNSrVw8uLi6oXLkynjx5gpMnT6JVq1bYuHEjgA/X+fXr1w9btmyBvb09unfvjlevXuHAgQNo3749Ll26pHGO1q1bo2bNmti5cyeysrLQo0cPxVyP/9SxY0fMnz8fU6dORa1atdC9e3fY2dnhzZs3uHfvHk6cOIG2bdti//79Gmch0luiTl7zL2fPnhVGjRol1K9fXzAzMxOqVq0qeHh4CH/++afSuJx5u/79p27duiq3u27dOqFu3bqCiYmJUKtWLWHlypVF8XaIRPX3338Lv/zyi+Dq6irUqVNHKFOmjGBkZCTY2toKTk5OwtatW1W+7uDBg0KrVq2EUqVKCRKJRJBKpcK9e/cEQRCEmTNnClKpNM859QRBEDIzM4WAgAChdevWgrm5uWBsbCxUq1ZN+Prrr4Xly5cLz58/FwRBELKysoSFCxcK3bp1E6pWrSqYmJgI1tbWQqtWrYTVq1cLWVlZim0eOnRIGDx4sFC/fn3B3NxcMDMzE+rWrSuMHTtWuH//fr4/kz///FNYsmSJ0K1bN6FmzZqCmZmZIl+vXr2Effv2qXxdcHCwYG9vL5iamio+kxyDBg1S+oxUiYqKElq2bCmULl1asLCwELp27SrExcUJQUFBglQqzTVHn52dXa55BPMam0PVfH7/Zf/+/YKjo6NgaWkpGBsbC1WrVhXc3NyEmJgYpXEZGRnCd999J1SuXFkwMjISatWqJcyfP194//69yv3m53uSY86cOYrPdPfu3R8de/LkScHDw0OwtbUVjIyMBGtra6Fx48bCxIkThQsXLqj13onoA4kgaM9DRd3d3REfH4/evXujUaNGSEpKQkBAAN68eYMzZ86gQYMGAD6cTtixYwfWr1+v9Hpzc3M4OTkpLVu9ejVGjhwJd3d3dO3aFSdOnMCWLVswf/78Tz6tQURERFScaVURjI+PR7NmzZQeFXXr1i00bNgQ7u7uiikjBg0ahN27d+PVq1cf3V5aWhqqVKmC1q1bK11M3r9/f+zduxd///23yqkYiIiIiPSBVt0s0qpVq1zPC61Zsybq16+f64HqgiBALpd/tAzGxMQgJSUFPj4+SstHjRqFt2/fYt++fQUXnoiIiKiY0aoiqIogCHjy5EmuC4jfvXuHMmXKoGzZsrCyssLo0aOVJqoFgMTERAC5Z+J3cHCAVCr9pIuciYiIiIo7rb9reOvWrXj06BHmzJmjWGZra4spU6bAwcEBcrkcBw4cwKpVq3D58mUcP35cMe9XUlISDAwMcpVIIyMjWFlZqfVYJyIiIiJdo9VF8MaNGxg1ahRat26tNMmov7+/0jgPDw/Url0b06dPR1hYGDw9PQF8uEbQyMhI5baNjY2V5sQiIiIi0jdaWwQfP34MJycnWFhYICws7KOPKgI+TDj6448/4ujRo4oiaGpqqvSIqn9KT0+HqampynXJycmIjo6GnZ1dnmOIiIiICkJaWhru3r2Lrl27qpxLszBpZRFMTU1Ft27d8OrVK8TFxX30weU5TExMYGlpiZSUFMUyGxsbZGdnIzk5WemDzczMREpKCmxtbVVuKzo6Gt7e3p/+RoiIiIjyKTg4WOXTdQqT1hXB9PR0ODs749atWzhy5Ajq1q2br9e9fv0aycnJKF++vGJZ48aNAQDnzp1Dt27dFMvPnz8PuVwOe3t7ldvKeQJDcHAw6tWrp+E7IV02fvx4LF26VOwYpIX43aCP4feDVLl+/Tq8vb0V/aMoaVURzM7OhqenJxISEhAeHo4WLVrkGpORkYHMzEyULl1aabmfnx8AwNHRUbGsY8eOsLS0RGBgoFIRDAwMhJmZWa7Jp3PknA6uV68eHBwcPvl9ke4xNzfnd4NU4neDPobfD/oYMS5H06oiOHHiRERGRsLZ2RnJyckIDg5WWu/t7Y2kpCQ0btwYffv2VTyHNDo6GgcOHEC3bt3Qs2dPxXgTExP4+flh1KhR8PDwQJcuXRAXF4etW7fC39+fk0kTERGRXtOqInj58mVIJBJERkYiMjJSaZ1EIoG3tzcsLCzg7OyMw4cPY9OmTcjOzkatWrUwb948TJo0Kdc2R44cCUNDQyxZsgQRERGoWrUqli1bBl9f36J6W0RERERaSauKYExMzH+OMTc3x+bNm9Xa7rBhwzBs2DBNYxERERHpJK1/sgiRNvLy8hI7AmkpfjfoY/j9IG3DIkikAf4yp7zwu0Efw+8HaRsWQSIiIiI9xSJIREREpKdYBImIiIj0FIsgERERkZ5iESQiIiLSUyyCRERERHqKRZCIiIhIT7EIEhEREekpFkEiIiIiPcUiSERERKSnWASJiIiI9BSLIBEREZGeYhEkIiIi0lMsgkRERER6ikWQiIiISE+xCBIRERHpKRZBIiIiIj3FIkhERESkp1gEiYiIiPQUiyARERGRnmIRJCIiItJTLIJEREREeopFkIiIiEhPsQgSERER6SkWQSIiIiI9xSJIREREpKdYBImIiIj0FIsgERERkZ5iESQiIiLSUyyCRERERHqKRZCIiIhIT7EIEhEREekpFkEiIiIiPcUiSERERKSnWASJiIiI9BSLIBEREZGeYhEkIiIi0lMsgkRERER6ikWQiIiISE+xCBIRERHpKRZBIiIiIj3FIkhERESkp1gEiYiIiPQUiyARERGRnmIRJCIiItJTLIJEREREeopFkIiIiEhPsQgSERER6SkWQSIiIiI9xSJIREREpKdYBImIiIj0FIsgERERkZ5iESQiIiLSUyyCRERERHqKRZCIiIhIT7EIEhEREekpFkEiIiIiPcUiSERERKSnWASJiIiI9BSLIBEREZGeYhEkIiIi0lMsgkRERER6ikWQiIiISE+xCBIRERHpKa0qgufOncPo0aPRoEEDlCpVCtWqVYOnpydu3ryZa+z169fh6OiI0qVLw8rKCgMGDEBycrLK7a5fvx716tWDqakpateujYCAgMJ+K0RERERar4TYAf5pwYIFiI+PR+/evdGoUSMkJSUhICAADg4OOHPmDBo0aAAAePDgAdq1awcLCwvMmzcPr1+/xuLFi3HlyhWcPXsWhoaGim2uXr0aI0eOhLu7OyZNmoQTJ07A19cX7969w+TJk8V6q0RERESi06oiOHHiRDRr1gwlSvx/LE9PTzRs2BDz58/Hli1bAAD+/v5IS0tDYmIiKleuDABo3rw5OnfujKCgIAwfPhwAkJaWhunTp0MmkyE0NBQAMHToUMjlcvj5+WHEiBEoW7ZsEb9LIiIiIu2gVaeGW7VqpVQCAaBmzZqoX78+bty4oVi2a9cuyGQyRQkEgE6dOqF27dqKwgcAMTExSElJgY+Pj9I2R40ahbdv32Lfvn2F9E6ISJ88e/YMgwcPRoMGDVCnTh00aNAAgwcPxrNnz8SORkRaLOd3h7u7u2gZtOqIoCqCIODJkydo2LAhAODhw4d49uwZmjZtmmtss2bNcODAAcXPiYmJAJBrrIODA6RSKS5duoR+/foVYnoi0nVPnz5F69atcfv2baXl165dQ1xcHOLj41G+fHmR0hGRtsrrd0dR06ojgqps3boVjx49gqenJwAgKSkJAGBjY5NrrI2NDVJSUpCVlaUYa2BggHLlyimNMzIygpWVFR49elTI6YlI102ZMiXPX+S3b9/mtchEpNLHfncUJa0+Injjxg2MGjUKrVu3xsCBAwF8uO4PAIyNjXONNzExUYwxNDREWloajIyMVG7b2NhYsS0iIk2dPXv2o+tjY2Nx8eLFIkpDRMVFbGys2BEAaHERfPz4MZycnGBhYYGwsDBIJBIAgKmpKQAgIyMj12vS09OVxpiamiIzM1Pl9tPT0xXjiIg09f79+4+uv3PnDpo0aVJEaYiI1KOVRTA1NRXdunXDq1evEBcXh4oVKyrW5ZwSzjlF/E9JSUmwsrJSTB9jY2OD7OxsJCcnK50ezszMREpKCmxtbT+aY/z48TA3N1da5uXlBS8vL43fGxHpln/f4PZv1atXR1hYWBGlISJtd/DgQURHR+P8+fN49+6d2HG0rwimp6fD2dkZt27dwpEjR1C3bl2l9ZUqVUL58uVx7ty5XK89e/Ys7O3tFT83btwYwIeJqrt166ZYfv78ecjlcqWxqixduhQODg6f8naISMc1b94c165dy3N9+/bt+XuEiBQcHBwwbdo0DB48GEFBQWLH0a6bRbKzs+Hp6YmEhATs3LkTLVq0UDmuV69eiIqKwoMHDxTLjh49ips3b6J3796KZR07doSlpSUCAwOVXh8YGAgzMzM4OTkVzhshIr3xz39k/luNGjWwcOHCIkxDRMXFwoULUaNGDbFjaNcRwYkTJyIyMhLOzs5ITk5GcHCw0npvb28AwLRp07Bz50506NABY8eOxevXr7Fo0SI0atQIgwcPVow3MTGBn58fRo0aBQ8PD3Tp0gVxcXHYunUr/P39OZk0EX2Sa9euYfjw4ejYsSMqV66M8+fP4/379yhRogSaN2+OhQsXcuoYIlKpfPnyiI+Px+TJkxEbG4s7d+6IkkMiCIIgyp5V6NChA06cOAFVkSQSCbKzsxU/X7t2DRMmTMDJkydhbGwMJycnLFmyROUv3XXr1mHJkiW4c+cOqlatitGjR8PX1zfPHBcvXkSTJk1w4cIFntIhIpWePXuGFi1awMzMDKdOnUKZMmXEjkRExZSYvUOrjgjGxMTke2z9+vVx8ODBfI0dNmwYhg0bpmksIiIl6enpcHFxwbt37xATE8MSSETFllYVQSIibScIAoYOHYqLFy/i+PHjqFatmtiRiIg0xiJIRKQGPz8/hISEYMeOHXne0EZEVFxo1V3DRETabPv27ZgxYwb8/Pzg4eEhdhwiok/GIkhElA/x8fEYNGgQ+vfvj+nTp4sdh4ioQLAIEhH9h7t376Jnz55o1qwZ1q5dq3jkJRFRccciSET0EampqZDJZChdujT27NkDY2NjsSMRERUY3ixCRJSH9+/fw9PTEw8ePMCZM2eUnllORKQLWASJiPIwbtw4HDlyBAcPHsz13HMiIl3AIkhEpMLKlSvxyy+/YPXq1fj666/FjkNEVCh4jSAR0b/s378f48aNw4QJEzBixAix4xARFRoWQSKif7hy5Qr69OkDmUyGhQsXih2HiKhQsQgSEf3P48ePIZPJUKNGDWzduhUGBgZiRyIiKlQsgkREANLS0uDi4oKsrCxERkaiVKlSYkciIip0vFmEiPSeXC7HoEGD8NtvvyEuLg6VK1cWOxIRUZFgESQivTdjxgyEhoZi165daNKkidhxiIiKDE8NE5Fe27JlC+bMmYP58+fDzc1N7DhEREWKRZCI9NbJkycxbNgwDBkyBJMnTxY7DhFRkWMRJCK9dPv2bbi4uKB169YIDAyERCIROxIRUZFjESQivfPy5UvIZDJYWlpi165dMDIyEjsSEZEoeLMIEemVrKwsuLu748mTJ0hISIClpaXYkYiIRMMiSER6QxAEjB49GidOnMChQ4dQq1YtsSMREYmKRZCI9MbSpUuxZs0abNiwAV999ZXYcYiIRMdrBIlIL0RERGDSpEmYMmUKBg8eLHYcIiKtwCJIRDrv0qVL6Nu3L1xdXeHv7y92HCIircEiSEQ67dGjR5DJZKhbty62bNkCqZS/9oiIcvA3IhHprLdv36JHjx6QSCSIiIhAyZIlxY5ERKRVeLMIEekkuVyOAQMG4MaNGzh58iRsbW3FjkREpHVYBIlIJ02bNg179uxBeHg47O3txY5DRKSVWASJSOds3LgRCxYswJIlS+Ds7Cx2HCIircVrBIlIpxw/fhwjRozAiBEjMH78eLHjEBFpNRZBItIZf/75J9zc3PDVV18hICAAEolE7EhERFqNRZCIdEJKSgpkMhkqVKiAnTt3wtDQUOxIRERaj9cIElGxl5mZCTc3N6SkpCAhIQFly5YVOxIRUbHAIkhExZogCPj2228RHx+Po0ePokaNGmJHIiIqNlgEiahYW7hwITZu3IgtW7bgyy+/FDsOEVGxwmsEiajY2r17N77//nv88MMP8Pb2FjsOEVGxwyJIRMXS+fPn4e3tDQ8PD8yaNUvsOERExRKLIBEVOw8ePECPHj3QqFEjBAUFQSrlrzIiIk3wtycRFStv3ryBs7MzDA0NER4eDlNTU7EjEREVW7xZhIiKjezsbPTt2xe3b9/GqVOnUKFCBbEjEREVayyCRFRsTJ48Gfv27UNUVBQaNmwodhwiomKPRZCIioU1a9bg559/xooVK9CtWzex4xAR6QReI0hEWu/IkSPw8fHBqFGjMGbMGLHjEBHpDBZBItJq169fh7u7Ozp37oxly5aJHYeISKewCBKR1kpOToZMJkPlypWxfft2lCjBq1mIiAoSf6sSkVbKyMiAq6sr3rx5g6NHj8Lc3FzsSEREOodFkIi0jiAIGDZsGM6dO4eYmBjY2dmJHYmk+slqAAAgAElEQVSISCexCBKR1pk7dy6Cg4Oxbds2tGrVSuw4REQ6i9cIEpFWCQ0NxY8//ohZs2ahT58+YschItJpLIJEpDUSEhIwcOBA9OvXDz/++KPYcYiIdB6LIBFphXv37qFHjx5wcHDAunXrIJFIxI5ERKTzWASJSHSvXr2CTCaDmZkZ9u7dCxMTE7EjERHpBd4sQkSiev/+Pfr06YP79+8jPj4e5cuXFzsSEZHeYBEkIlFNmDABhw4dwoEDB1C/fn2x4xAR6RUWQSISzS+//IKVK1ciMDAQnTt3FjsOEZHe4TWCRCSKgwcPwtfXF+PGjcO3334rdhwiIr3EIkhERe7333+Hh4cHunfvjsWLF4sdh4hIb7EIElGRevLkCWQyGapXr46QkBAYGBiIHYmISG+xCBJRkUlLS4OLiwsyMjIQFRWF0qVLix2JiEiv8WYRIioSgiBgyJAhuHz5MmJjY1GlShWxIxER6T0WQSIqEjNnzsT27duxc+dONGvWTOw4REQEnhomoiKwdetWzJ49G/7+/nB3dxc7DhER/Q+LIBEVqlOnTmHIkCEYOHAgvv/+e7HjEBHRP7AIElGh+euvv+Di4oKWLVtizZo1kEgkYkciIqJ/yPc1gqVLl4ZEIoEgCPneeM54iUSCqKgotGvXTqOQRFT8vHz5EjKZDGXLlsXu3bthZGQkdiQiIvqXfBfBt2/fom3btqhevbpaO8jIyMCOHTvw/v17tcMRUfGUlZUFDw8PJCUl4cyZM7CyshI7EhERqaDWXcMjRoxAv3791NpBcnIyduzYke/xb9++xcKFC5GQkICzZ8/i5cuX2LhxIwYOHKg0btCgQdi8eXOu19epUwfXr1/PtXz9+vVYvHgx7t69iypVqsDX1xejR49W670Q0X8TBAG+vr6IiYnBoUOHUKdOHbEjERFRHvJdBH19fVGvXj21d1CyZEn4+vrme86wZ8+ewc/PD9WqVYO9vT2OHz+e53VFxsbGWL9+vdIyc3PzXONWr16NkSNHwt3dHZMmTcKJEyfg6+uLd+/eYfLkyWq/JyLK2/Lly/Hrr79i3bp16NChg9hxiIjoI/JdBJctW6bRDkqWLKnWa21tbfH48WNYW1vjwoULH51vzNDQEH379v3o9tLS0jB9+nTIZDKEhoYCAIYOHQq5XA4/Pz+MGDECZcuWzXc+IspbVFQUJkyYgO+++w5Dhw4VOw4REf0Hrbtr2MjICNbW1gDwnzemCIIAuVyOV69e5TkmJiYGKSkp8PHxUVo+atQovH37Fvv27fv00ESEy5cvo0+fPujZsyfmz58vdhwiIsqHT36yyN9//41Lly4hNTUVcrk81/oBAwZ86i7y9O7dO5QpUwbv3r2DhYUFvLy8sGDBApiZmSnGJCYmAgCaNm2q9FoHBwdIpVJcunRJ7eseiUhZUlISZDIZ6tSpg+DgYEilWvdvTCIiUkHjIpieno4BAwZg165dHz1yV1hF0NbWFlOmTIGDgwPkcjkOHDiAVatW4fLlyzh+/DgMDAwAfPgLysDAAOXKlVN6vZGREaysrPDo0aNCyUekL969e4cePXpALpcjIiJC6R9iRESk3TQugtOmTcOePXvg7++PVq1a4auvvkJQUBBsbGywfPlyPHz4UOVdvQXF399f6WcPDw/Url0b06dPR1hYGDw9PQF8uEYwr/nLjI2NkZaWVmgZiXSdXC7HgAEDcO3aNcTFxaFSpUpiRyIiIjVofP4mLCwMgwYNwpQpU1C/fn0AQOXKldG5c2dERkaibNmy+OWXXwosaH6MHz8eUqkUR48eVSwzNTVFZmamyvHp6ekwNTUtqnhEOueHH37A7t27ERISAgcHB7HjEBGRmjQ+Ivj06VO0aNECABRl6u3btwA+PFGkV69emD17Nn799dcCiJk/JiYmsLS0REpKimKZjY0NsrOzkZycrHR6ODMzEykpKbC1tc1ze+PHj881HY2Xlxe8vLwKPjxRMRMUFIR58+Zh0aJF6Nmzp9hxiIiKhW3btmHbtm1Ky1JTU0VK8wlFsEKFCnj+/DkAwMzMDGXLlsWNGzfg7OwMAHj9+jXS09MLJmU+vX79GsnJyShfvrxiWePGjQEA586dQ7du3RTLz58/D7lcDnt7+zy3t3TpUh7lIFIhNjYWI0aMwLBhwzBx4kSx4xARFRuqDihdvHgRTZo0ESWPxqeGmzdvjpMnTyp+7tGjBxYvXozg4GBs3rwZS5cuRcuWLQsk5L9lZGTg9evXuZb7+fkBABwdHRXLOnbsCEtLSwQGBiqNDQwMhJmZGZycnAolI5GuunnzJtzc3NC2bVusWrUqzwnfiYhI+2l8RNDX1xdhYWFIT0+HiYkJZs+ejfj4eMVdwjVq1MCKFSs02nZAQABevnypuKM3IiIC9+/fV+w3JSUFjRs3Rt++fRWPr4qOjsaBAwfQrVs3pdNUJiYm8PPzw6hRo+Dh4YEuXbogLi4OW7duhb+/PyeTJlJDSkoKZDIZypcvj7CwMBgaGoodiYiIPoFE+K9Zm9WQnZ2NK1euwMDAAPXq1UOJEpr1zOrVq+PevXsfAv7vaIMgCJBIJLhz5w7Mzc0xZswYnDlzBo8ePUJ2djZq1aqFfv36YdKkSYqpY/5p3bp1WLJkCe7cuYOqVati9OjR8PX1Vbn/nEO0Fy5c4Klhov/JzMyEo6MjLl++jISEBNSsWVPsSEREOkHM3qF2Uzt48CCWLVuGu3fvwsrKCr1798a4ceMAAAYGBh+95i6/7ty5859j1J2aZtiwYRg2bJimkYj0miAI8PHxwcmTJ3HkyBGWQCIiHaFWEYyNjYWTkxMEQYCVlRVu3bqF+Ph4PHz4EIsWLSqsjEQkssWLF2P9+vXYtGkT2rVrJ3YcIiIqIGrdLOLv7w9ra2v89ttvePbsGZ49e4aOHTti1apVnJiZSEft3bsXU6ZMwbRp0wr1kZFERFT01CqCv//+O3x8fPD5558DACwsLODv74+0tDRcvXq1UAISkXguXryIfv36oVevXoq78omISHeoVQQfP36Mzz77TGlZ9erVAUDldC5EVHw9fPgQzs7OaNCgATZt2gSpVOPZpoiISEup9Zs9585dpQ387y+HArz5mIhE9ubNGzg7O8PAwAAREREoWbKk2JGIiKgQqH3X8ObNm3HmzBnFzznXBgYEBGDv3r25xms6lyARiSM7Oxve3t64efMmTp06hYoVK4odiYiIConaRfDQoUM4dOhQruWqSiDAIkhU3Hz//feIjIxEREQEGjVqJHYcIiIqRGoVQblcXlg5iEgLrFu3DosXL8ayZcv4+EUiIj3Aq7+JCABw9OhRjBw5EiNHjszzqTtERKRbWASJCDdu3IC7uzs6duyIFStW5LopjIiIdJNap4adnZ3V/gsiIiJCrfFEVLSSk5Mhk8lga2uL0NBQjZ8RTkRExY9av/H37dsHY2NjxV2E/zVlDI8qEGm3jIwMuLm54dWrV0hISIC5ubnYkYiIqAipVQQrVaqEhw8foly5cujXrx88PT1hY2NTWNmIqBAJgoARI0bg7NmzOHbsmGJyeCIi0h9qXSN4//59HDt2DI0bN4afnx+qVKmCr7/+Ghs2bMCrV68KKyMRFYJ58+Zh8+bN2LBhA1q3bi12HCIiEoFaRVAikeCrr77CmjVrkJSUhF27dsHS0hKjR4+GtbU13NzcEBYWhvT09MLKS0QFYOfOnZg+fTpmzJiBvn37ih2HiIhEovFdw0ZGRujZsydCQ0Px5MkTrF69Go8fP4anpycWLVpUkBmJqACdPXsWAwYMgJeXF2bMmCF2HCIiEtEnTx+TkZGBQ4cOITw8HImJiTAxMUG1atUKIhsRFbD79++jR48eaNy4MTZs2MAbuoiI9JxGRTA7OxsHDx7EwIEDYW1tDS8vL6Snp2Pt2rV4+vQpBgwYUNA5iegTvX79GjKZDKampti7dy9MTEzEjkRERCJT667hU6dOISQkBDt37sTz58/RqlUrzJs3Dx4eHihXrlxhZSSiT5SdnQ0vLy/cu3cPp0+fhrW1tdiRiIhIC6hVBNu2bQtTU1N069YNXl5esLOzg0Qiwf3793H//n2Vr3FwcCiQoESkuYkTJ+LgwYPYt28fGjRoIHYcIiLSEmo/QiAtLQ27d+/G7t27/3OsRCJBdna2RsGIqGAEBgZi+fLl+OWXX9C1a1ex4xARkRZRqwhu2LChsHIQUSGIjo7GmDFj4OvrCx8fH7HjEBGRllGrCA4aNKiQYhBRQbt69So8PDzQtWtX/Pzzz2LHISIiLfTJ08cQkfZ5+vQpZDIZqlWrhu3bt8PAwEDsSEREpIXyXQR3796Nv//+W+0dZGVlYffu3Xj27JnaryUi9aWnp8PFxQVpaWmIiopC6dKlxY5ERERaKt9F0N3dHXFxcWrvIDU1Fe7u7rhy5YraryUi9QiCgCFDhiAxMRERERGoWrWq2JGIiEiLqXWNYFxcHN6/f6/WDl6/fq3WeCLS3OzZs7Ft2zaEhoaiefPmYschIiItp1YRXL16NVavXq3RjgRB0Oh1RJQ/27Ztw8yZMzFnzhz07t1b7DhERFQM5LsI/vXXX5+0Ixsbm096PRHl7fTp0xg8eDAGDBiAadOmiR2HiIiKiXwXQTs7u0KMQUSaunPnDlxcXNC8eXOsWbMGEolE7EhERFRMcPoYomIsNTUVMpkMZcqUwe7du2FsbCx2JCIiKkbUfsQcEWmH9+/fw8PDA48ePUJ8fDzKlSsndiQiIipmWASJiiFBEDB27FgcO3YMBw8eRN26dcWORERExRCLIFExtHLlSqxatQpr1qxBp06dxI5DRETFFK8RJCpm9u3bh/Hjx2PixIkYPny42HGIiKgYK9AjgnK5HDExMcjMzMSXX37JR1sRFbDffvsNffr0gUwmw4IFC8SOQ0RExZzGRwSnT5+ODh06KH4WBAFdunRB586d4eTkhM8//xy3b98ukJBEBDx+/BgymQw1a9bE1q1bYWBgIHYkIiIq5jQugrt27UKzZs0UP4eFheHYsWOYO3cuoqKikJ2djRkzZhRISCJ9l5aWhp49eyI7OxuRkZEoVaqU2JGIiEgHaHxq+OHDh6hVq5bi5927d6NevXqYOnUqAMDHxweBgYGfnpBIz8nlcgwcOBC///47Tpw4gcqVK4sdiYiIdITGRwRLlCiBjIwMAB/+ojp69CgcHR0V662trfHs2bNPT0ik53766SeEhYUhODgYTZo0ETsOERHpEI2LYIMGDbBlyxakpKQgKCgIycnJcHJyUqy/f/8+J7gl+kSbN2/G3LlzMX/+fLi6uoodh4iIdIzGp4ZnzJgBmUymKHtt2rRRunlk3759StcQEpF64uLiMGzYMAwZMgTfffed2HGIiEgHaVwEO3fujIsXL+Lw4cOwsLCAp6en4mH3L168QNu2bdGzZ88CC0qkT27dugVXV1d8+eWXCAwMVPx/i4iIqCB90jyCDRo0QIMGDXItt7CwwLJlyz5l00R668WLF5DJZLCyskJYWBiMjIzEjkRERDpK4yL46tUrpKamokqVKoplDx8+xK+//orMzEz06tULzZs3L5CQRPoiKysL7u7uePbsGc6cOQNLS0uxIxERkQ7TuAh+8803uHPnDs6cOQPgQzFs2bIlHj58CKlUiuXLl+PAgQNK1w0SUd4EQcCoUaMQFxeHw4cPK03PREREVBg0vmv45MmTSncJBwcHIykpCadOncKLFy/QqFEjzJ07t0BCEumDn3/+GWvXrsWaNWvQvn17seMQEZEe0LgIJicnK01sGxERgTZt2qBVq1YoXbo0+vfvj0uXLhVISCJdFx4eju+++w7ff/89Bg0aJHYcIiLSExoXwbJly+Lx48cAPjz+Ki4uDl26dFGsL1GiBN69e/fpCYl0XGJiIvr27Qs3NzceRScioiKl8TWCrVu3xqpVq1CnTh0cPHhQ8SzUHDdv3kSlSpUKJCSRrnr48CGcnZ1Rv359bN68GVKpxv82IyIiUpvGRXD+/Pno2rUr3N3dAQATJkzA559/DgB4//49QkNDlR45R0TK3r59ix49ekAikSAiIgIlS5YUOxIREekZjYtgrVq1cOPGDVy7dg3m5uaoXr26Yl1aWhoCAgJgb29fICGJdI1cLoe3tzf++OMPnDp1CjY2NmJHIiIiPfRJE0obGRmpLHulS5eGi4vLp2yaSKdNnToV4eHhCA8PxxdffCF2HCIi0lOfdEFSamoq5s2bhy5duqBx48Y4e/YsACAlJQU///wzbt26VSAhiXTJ+vXrsXDhQixZsgTOzs5ixyEiIj2m8RHBBw8eoF27dnjw4AFq1qyJGzdu4M2bNwA+PGJu9erVuHv3LlasWFFgYYmKu5iYGHz77bf45ptvMG7cOLHjEBGRntO4CH733Xd48+YNEhMTUaFCBVhbWyvWSSQSuLi4ICoqqkBCEumCP//8E7169UKHDh2wcuVKSCQSsSMREZGe0/jU8KFDhzBmzBg0aNBA5frq1avj/v37Ggcj0iXPnz+Hk5MTKlasiNDQUBgaGoodiYiISPMjgmlpaUpHAf/t9evXmm6aSKdkZmaiV69eePnyJRISElC2bFmxIxEREQH4hCOC9erVQ2xsbJ7rw8PD0bhxY003T6QTBEHAN998g/j4eOzduxefffaZ2JGIiIgUNC6C48ePx44dO7BgwQKkpqYCALKzs3Hz5k14e3vj9OnTGD9+fIEFJSqOFixYgKCgIGzYsAFt2rQROw4REZESjU8Ne3t74969e5g+fTqmTZsGAHB0dIQgCDAwMMC8efPg6upaYEGJiptdu3Zh6tSp+PHHH9GvXz+x4xAREeXySRNKT58+Hd7e3ti9ezdu3rwJuVyOmjVrws3NjafASK+dP38e/fv3h6enJ2bNmiV2HCIiIpU0KoJv375F27ZtMWLECHz77bc8BUz0D3///TecnZ3xxRdfYOPGjZwmhoiItJZG1wiamZnh7t27/AuO6F/evHkDZ2dnGBsbY+/evTA1NRU7EhERUZ40vlnE0dER0dHRBZmFqFjLzs6Gl5cX/vrrL0RFRaFChQpiRyIiIvoojYvgjz/+iD///BPe3t44efIkHj58iJSUlFx/iPTFd999h/379yM0NBSff/652HGIiIj+k8ZFsEGDBrh27RpCQkLQrl07VKlSBeXKlVP6U758ebW3+/btW8yYMQOOjo6wtLSEVCrFpk2bVI69fv06HB0dUbp0aVhZWWHAgAFITk5WOXb9+vWoV68eTE1NUbt2bQQEBKidjSgvq1evxtKlS7F8+XI4OjqKHYeIiChfNL5r+KeffirIHArPnj2Dn58fqlWrBnt7exw/flzltYgPHjxAu3btYGFhgXnz5uH169dYvHgxrly5grNnzyo9wmv16tUYOXIk3N3dMWnSJJw4cQK+vr549+4dJk+eXCjvg/TH4cOHMWrUKIwePRqjR48WOw4REVG+aVwEZ86cWYAx/p+trS0eP34Ma2trXLhwAc2aNVM5zt/fH2lpaUhMTETlypUBAM2bN0fnzp0RFBSE4cOHA/jwKLzp06dDJpMhNDQUADB06FDI5XL4+flhxIgRfOQXaez69evo3bs3unTpgqVLl4odh4iISC0anxoeMmQIEhIS8lx/9uxZDBkyRO3tGhkZKZ5hLAhCnuN27doFmUymKIEA0KlTJ9SuXVtR+AAgJiYGKSkp8PHxUXr9qFGj8PbtW+zbt0/tjETAh6PXTk5OqFKlCrZv344SJT5pWk4iIqIip3ERDAoKwu3bt/Nc/9dffyEoKEjTzX/Uw4cP8ezZMzRt2jTXumbNmiExMVHxc85//3usg4MDpFIpLl26VCgZSbdlZGTA1dUVb9++RVRUFMqUKSN2JCIiIrUV2iGMpKSkQptDLSkpCQBgY2OTa52NjQ1SUlKQlZUFQ0NDJCUlwcDAAOXKlVMaZ2RkBCsrKzx69KhQMpLuEgQBw4YNw/nz53H8+HFUq1ZN7EhEREQaUasIhoeHIzw8XHHKds2aNThy5EiucS9evMCRI0fyvL7vU6WlpQEAjI2Nc60zMTFRjDE0NERaWhqMjIxUbsfY2FixLaL8mjNnDoKDg7F9+3a0bNlS7DhEREQaU6sIXr16FaGhoYq7eBMSEnDhwgWlMRKJBGZmZmjfvj1+/vnngkv6DzlHGjMyMnKtS09PVxpjamqKzMxMldtJT0/nkx9ILTt27MBPP/2E2bNnw9PTU+w4REREn0StIjht2jRMmzYNACCVSrFu3Tr069evUIJ9TM4p4ZxTxP+UlJQEKysrxfQxNjY2yM7ORnJystLp4czMTKSkpMDW1jbP/YwfPx7m5uZKy7y8vODl5VUQb4OKmTNnzmDgwIHw9vbGDz/8IHYcIiIqhrZt24Zt27YpLUtNTRUpzSdcIyiXywsyh1oqVaqE8uXL49y5c7nWnT17Fvb29oqfGzduDAA4d+4cunXrplh+/vx5yOVypbH/tnTpUjg4OBRgciqu7t69i549e6Jp06ZYt24dn7NNREQaUXVA6eLFi2jSpIkoeTS+a/jfrl+/jtmzZ8PHxwfLly/Hq1evCmrTKvXq1QtRUVF48OCBYtnRo0dx8+ZN9O7dW7GsY8eOsLS0RGBgoNLrAwMDYWZmBicnp0LNScXfq1ev4OzsjFKlSmHPnj0qr00lIiIqjtQ6IhgQEIAVK1bg9OnTSqdZIyMj4e7ujqysLMWyFStWICEhIdfduvndz8uXLxV39EZEROD+/fsAAF9fX5QpUwbTpk3Dzp070aFDB4wdOxavX7/GokWL0KhRIwwePFixLRMTE/j5+WHUqFHw8PBAly5dEBcXh61bt8Lf35+TSdNHvX//Hp6envj7778RHx+v0WMTiYiItJaghs6dOwtdu3ZVWpaVlSVYW1sLpqamQlBQkHDlyhVhwYIFQokSJYSxY8eqs3kFOzs7QSKRCBKJRJBKpYJUKlX897179xTjrl69KnTt2lUwMzMTLC0thf79+wtPnz5Vuc21a9cKdevWFYyNjYVatWoJy5cvz3P/Fy5cEAAIFy5c0Cg/6Y7Ro0cLBgYGwuHDh8WOQkREOkrM3qHWEcFr164pHt2WIyYmBs+ePcPUqVMxcOBAAMDnn3+Oy5cv48CBA1i2bJna5fTOnTv5Gle/fn0cPHgwX2OHDRuGYcOGqZ2F9FdAQAACAgLw66+/4uuvvxY7DhERUYFT6xrB58+fo0qVKkrLcuYRdHV1VVreunVr3Lt37xPjEYnjwIEDGDt2LMaPH49vvvlG7DhERESFQq0iWKFCBTx+/FhpWVxcHEqWLJnr7lsjI6M8J3Im0ma///47PD094eTkhEWLFokdh4iIqNCoVQSbNm2KTZs2Ke4Ivnr1Ks6dO4euXbuiRAnls8x//PEHKleuXHBJiYrAkydPIJPJ8NlnnyEkJAQGBgZiRyIiIio0al0jOGPGDDRt2hS1a9dGgwYNcP78eQDA1KlTc43ds2cPOnToUDApiYpAWloaevbsiczMTERGRqJUqVJiRyIiIipUah0RbNiwIWJiYtCkSRM8fPgQrVq1woEDB9C0aVOlcTExMTA1NVWaz49Im8nlcgwePBi//fYbIiIicl0LS0REpIvUfrJI69atsW/fvo+O6dChA37//XeNQxEVtZkzZ2LHjh3YtWtXrn/YEBER6aoCe7IIUXEVHBwMPz8/zJs3D25ubmLHISIiKjIsgqTXTp48iaFDh2Lw4MGYMmWK2HGIiIiKFIsg6a2//voLrq6uaNWqFX799VdIJBKxIxERERUpFkHSSy9fvoSTkxMsLCywa9cuznlJRER6Se2bRYiKu6ysLPTu3RtPnjzBmTNnYGVlJXYkIiIiUbAIkl4RBAFjxozB8ePHcejQIdSuXVvsSERERKJhESS9smzZMqxevRrr16/nhOdERKT3eI0g6Y3IyEhMnDgRkydPxpAhQ8SOQ0REJDoWQdILly5dgpeXF1xcXDBv3jyx4xAREWkFFkHSeUlJSXB2dkadOnWwZcsWSKX82hMREQEsgqTj3r17hx49ekAQBERGRsLMzEzsSERERFqDN4uQzpLL5ejfvz+uXbuGkydPwtbWVuxIREREWoVFkHTW9OnTsWfPHuzduxeNGzcWOw4REZHWYREknbRx40bMnz8fixcvRo8ePcSOQ0REpJV4jSDpnNjYWHzzzTcYPnw4JkyYIHYcIiIircUiSDrl5s2bcHNzQ7t27fDLL79AIpGIHYmIiEhrsQiSzkhJSYGTkxPKly+PnTt3wtDQUOxIREREWo3XCJJOyMzMRK9evZCSkoKEhARYWFiIHYmIiEjrsQhSsScIAkaOHIlTp07h6NGjqFGjhtiRiIiIigUWQSr2Fi1ahA0bNmDTpk1o27at2HGIiIiKDV4jSMXanj178P3332P69OkYMGCA2HGIiIiKFRZBKrYuXLiAfv36wd3dHbNnzxY7DhERUbHDIkjF0oMHD+Ds7IyGDRti06ZNkEr5VSYiIlIX//akYufNmzdwdnaGoaEhwsPDYWpqKnYkIiKiYok3i1Cxkp2djX79+uHWrVs4deoUKlasKHYkIiKiYotFkIqVKVOmICoqCpGRkWjUqJHYcYiIiIo1FkEqNtauXYslS5Zg+fLl6N69u9hxiIiIij1eI0jFwtGjR+Hj4wMfHx+MGTNG7DhEREQ6gUWQtN6NGzfQq1cvdOrUCcuXL4dEIhE7EhERkU5gESStlpycDCcnJ1SqVAk7duxAiRK8moGIiKig8G9V0loZGRlwdXXF69evceTIEZibm4sdiYiISKewCJJWEgQBw4cPx7lz5xATE4Pq1auLHYmIiEjnsAiSVvL398eWLVsQEhKCVq1aiR2HiIhIJ/EaQdI6oaGh+OGHHzBz5kx4eXmJHYeIiEhnsQiSVklISMDAgQPRt29f/PTTT2LHISIi0mksgqQ17iwAa4YAACAASURBVN27h549e8LBwQHr16/nNDFERESFjEWQtMKrV6/g7OwMU1NT7NmzByYmJmJHIiIi0nm8WYRE9/79e3h5eeHevXuIj4+HtbW12JGIiIj0AosgiW7ixImIjo7G/v37Ub9+fbHjEBER6Q0WQRLVqlWrsGLFCqxatQpdunQROw4REZFe4TWCJJro6Gj4+vpi7NixGDlypNhxiIiI9A6LIIni6tWr8PDwgKOjI5YsWSJ2HCIiIr3EIkhF7unTp5DJZLCzs8O2bdtgYGAgdiQiIiK9xGsEqUilp6fDxcUFaWlpiI2NRenSpcWOREREpLdYBKnICIKAIUOGIDExEbGxsahatarYkYiIiPQaiyAVmVmzZmHbtm3YuXMnmjdvLnYcIiIivcdrBKlIhISEYNasWZg7dy7c3d3FjkNERERgEaQicPr0aQwe/H/t3XlYlFX/BvD7GURWEWRRUIsg0bIMRFRMfYVcayAXQMUFFDJz39fMFiOX1F4y91JeMc0gV3IJldJAwMQ0tVxyCbWEQCRk5/z+KOcXDijiDM8s9+e6uGTOnDlzP3h8+HqeZUYgLCwMs2fPljsOERER/YOFIGnV5cuX0bdvX3To0AFr1qyBJElyRyIiIqJ/sBAkrcnLy4NSqYSNjQ2++uormJmZyR2JiIiI/oUXi5BWlJWVISQkBDdu3MCxY8fg4OAgdyQiIiK6DwtB0jghBCZMmIBDhw5h//79aNmypdyRiIiIqAosBEnjoqOjsWrVKqxbtw7+/v5yxyEiIqJq8BxB0qiEhARMmTIF06ZNQ2RkpNxxiIiI6AFYCJLGnDp1CoMGDUJAQAAWLlwodxwiIiJ6CBaCpBE3b96EUqlEixYtsHnzZpiYmMgdiYiIiB6ChSA9trt37+LVV19FeXk5du/eDSsrK7kjERERUQ3wYhF6LBUVFQgLC8OZM2dw5MgRNG3aVO5IREREVEMsBOmxzJs3D/Hx8fjqq6/Qtm1bueMQERHRI2AhSLUWExODqKgoLF68GH379pU7DhERET0ivT1HMCkpCQqFosqvtLS0Sn3PnTuH3r17o0GDBrC3t8fw4cORnZ0tU3LD8N133+G1115DREQEpk2bJnccIiIiqgW9XxGcOHEifHx8KrW5u7urvs/MzETXrl1hZ2eHDz74APn5+fjwww9x+vRppKWlwdTUtK4j672LFy+iX79+6Ny5M1auXAlJkuSORERERLWg94Vgly5d0L9//2qfj4qKQmFhITIyMtCsWTMAQPv27dGjRw9s3LgRr732Wl1FNQi5ubl45ZVX4ODggPj4eNSvX1/uSERERFRLento+B4hBPLz81FWVlbl8/Hx8VAqlaoiEABeeukleHh4YNu2bXUV0yCUlpYiKCgI2dnZSEhIgJ2dndyRiIiI6DHofSE4YsQINGzYEBYWFvD398cPP/ygeu769evIyspCu3bt1F7n4+ODjIyMuoyq14QQGDNmDI4cOYLt27fj6aefljsSERERPSa9PTRsZmaGoKAgvPzyy3BwcMCZM2fw4YcfokuXLkhOToanpydu3rwJAHB2dlZ7vbOzM3JyclBaWsrzBGtg6dKlWL9+PTZu3IiuXbvKHYeIiIg0QG8LQV9fX/j6+qoeK5VKBAUFoU2bNpg9ezb27t2LwsJCAH8XjfczNzcHABQWFrIQfIgdO3ZgxowZmD17NsLCwuSOQ0RERBqi94eG/83d3R2BgYE4fPgwhBCwsLAAABQXF6v1LSoqAgBVH6raiRMnMGTIEAwYMAALFiyQOw4RERFpkN6uCFanefPmKCkpQUFBgeqQ8L1DxP928+ZN2NvbP3A1cPLkyWjYsGGltsGDB2Pw4MGaDa2jrl+/joCAALRu3RoxMTFQKAzq/w1ERER1bsuWLdiyZUultry8PJnSGGAh+Ouvv8LCwgLW1tawtraGo6Mj0tPT1fqlpaXB09PzgWMtX77caD82raCgAAEBATAxMcHOnTthaWkpdyQiIiK9V9WC0okTJ+Dt7S1LHr1d4snKylJr+/HHH7Fr1y707NlT1TZgwADs2bMHmZmZqraDBw/iwoULCA4OrpOs+qaiogJDhw7FhQsXsHv37iovtiEiIiL9p7crggMHDoSlpSV8fX3h5OSEs2fPYu3atbC2tsbChQtV/ebMmYMvv/wSfn5+mDhxIvLz87FkyRK0adMGI0aMkHELdNesWbOwa9cu7Ny5Ey+88ILccYiIiEhL9LYQ7NevHzZv3ozly5fjzp07cHJyQlBQEObPnw83NzdVv2bNmuHbb7/FlClTMGvWLJiZmUGpVGLp0qW8WrgK69evx5IlS7B8+XIolUq54xAREZEW6W0hOH78eIwfP75GfZ999lns27dPy4n036FDh/DGG29g9OjRmDhxotxxiIiISMv09hxB0qxffvkFAwYMgJ+fH6KjoyFJktyRiIiISMtYCBL+/PNPvPLKK3B2dsa2bdt4yJyIiMhI6O2hYdKM4uJi9O/fH3l5eUhNTYWtra3ckYiIiKiOsBA0YkIIvP766zh27BgOHTpU6SIbIiIiMnwsBI3YwoULERMTg9jYWLz44otyxyEiIqI6xnMEjVRcXBzmzJmDt956C0OGDJE7DhEREcmAhaARSk9Px7BhwzBo0CC8/fbbcschIiIimbAQNDLXrl1DYGAgPD09sWHDBt4mhoiIyIixEDQi+fn5CAgIgJmZGXbs2AFzc3O5IxEREZGMeLGIkSgvL8fgwYNx+fJlJCcno3HjxnJHIiIiIpmxEDQS06ZNw969e5GQkIDnnntO7jhERESkA1gIGoHVq1fjo48+wooVK9C7d2+54xAREZGO4DmCBu7AgQMYN24cxo8fj7Fjx8odh4iIiHQIC0EDdvbsWQQHB6Nnz55YtmyZ3HGIiIhIx7AQNFBZWVlQKpV44oknsHXrVtSrx7MAiIiIqDJWBwaoqKgIffv2xd27d3H48GHY2NjIHYmIiIh0EAtBAyOEQEREBE6cOIGkpCQ8+eSTckciIiIiHcVC0MC89957+Pzzz/HFF1+gQ4cOcschIiIiHcZzBA3I1q1bMX/+fLz33nsICQmROw4RERHpOBaCBiIlJQXh4eEYNmwY5s6dK3ccIiIi0gMsBA3AlStX8Oqrr8LHxwfr1q2DJElyRyIiIiI9wEJQz+Xl5UGpVKJBgwbYvn07zMzM5I5EREREeoIXi+ixsrIyDBw4EJmZmTh27BgcHBzkjkRERER6hIWgHps0aRISExOxb98+tGrVSu44REREpGdYCOqpjz/+GJ988gnWrFmD7t27yx2HiIiI9BDPEdRDX3/9NSZNmoQpU6Zg1KhRcschIiIiPcVCUM+cPn0agwYNglKpxOLFi+WOQ0RERHqMhaAe+f3336FUKuHu7o7NmzfDxMRE7khERESkx1gI6onCwkL07dsXpaWl2L17N6ytreWORERERHqOF4vogYqKCoSHh+PUqVM4cuQImjVrJnckIiIiMgAsBPXA/PnzsW3bNsTHx8Pb21vuOERERGQgeGhYx23atAkLFizAwoUL0b9/f7njEBERkQFhIajDjh49isjISIwcORIzZsyQOw4REREZGBaCOurSpUvo27cvOnXqhFWrVkGSJLkjERERkYFhIaiDbt++DaVSiUaNGiE+Ph7169eXOxIREREZIF4somNKS0sRFBSEP/74A6mpqWjUqJHckYiIiMhAsRDUIUIIjBs3Dt999x0OHDiAFi1ayB2JiIiIDBgLQR2yfPlyrF27Fp999hm6desmdxwiIiIycDxHUEfs2rUL06ZNw8yZMzFixAi54xAREZERYCGoA06ePInQ0FD069cPUVFRcschIiIiI8FCUGY3btyAUqlEq1atsGnTJigU/CshIiKiusGqQ0YFBQUIDAyEJEnYtWsXLC0t5Y5ERERERoQXi8ikoqICw4cPx88//4yjR4/CxcVF7khERERkZFgIymTOnDnYvn07du7cCU9PT7njEBERkRFiISiDDRs2YNGiRVi6dCkCAgLkjkNERERGiucI1rGkpCSMGjUKo0aNwuTJk+WOQ0REREaMhWAdOn/+PPr3749u3bphxYoVkCRJ7khERERkxFgI1pGcnBwolUo0btwYX375JUxNTeWOREREREaO5wjWgZKSEvTv3x85OTlITU2Fra2t3JGIiIiIWAhqmxACo0ePRkpKCg4ePAh3d3e5IxEREREBYCGodYsXL8aGDRuwadMmdO7cWe44RERERCo8R1CLvvrqK8yaNQtvvvkmhg4dKnccIiIiokpYCGrJ8ePHMXToUISEhOCdd96ROw4RERGRGhaCWpCZmYnAwEC0adMGGzduhELBHzMRERHpHlYoGvbXX38hICAApqam2LlzJywsLOSORERERFQlXiyiQeXl5QgNDcWlS5fw/fffo3HjxnJHIiIiIqoWC0ENmjFjBhISErBnzx48//zzcschIiIieiAWghqydu1aLFu2DNHR0ejTp4/ccYiIiIgeiucIakBiYiLGjBmDsWPHYvz48XLHISIiIqoRFoKP6dy5cwgKCkKPHj3w0UcfyR2HiIiIqMZYCD6G7OxsKJVKNGvWDFu3bkW9ejzSTkRERPqDlUstFRcXo1+/fvjrr79w8OBBNGzYUO5IRERERI+EhWAtCCEQGRmJ9PR0HD58GK6urnJHIiIiInpkLARr4f3330dsbCy2bNkCX19fueMQERER1QrPEXxE27Ztw7x58/DOO+9g0KBBcschIiIiqjUWgo8gNTUVYWFhGDJkCObNmyd3HCIiIqLHYhSFYHFxMWbOnAkXFxdYWlqiY8eOSExMfKQxrl69isDAQLRt2xbr16+HJElaSktERERUN4yiEAwPD8fy5csxbNgwREdHw8TEBC+//DK+//77Gr3+zp07UCqVsLKywo4dO2Bubq7lxERERETaZ/CFYFpaGr744gssXLgQixYtQmRkJA4dOoQnn3wSM2bMeOjry8rKMGjQIFy7dg179uyBo6NjHaQmXbdlyxa5I5CO4tygB+H8IF1j8IVgXFwc6tWrh1GjRqnazMzMEBERgZSUFFy/fv2Br58yZQoOHDiAuLg4PPvss9qOS3qCO3OqDucGPQjnB+kag799TEZGBjw8PGBtbV2p3cfHBwBw8uRJNG3atMrX9uzZE3/++ScWL16MHj16aD0rERERUV0y+BXBmzdvwtnZWa39XtuNGzeqfe2ff/4JAFizZg2ysrK0E5CIiIhIJgZfCBYWFsLMzEyt/d4FH4WFhQ8d49KlSzU6n5CIiIhInxj8oWELCwsUFxertRcVFamev19VxeG3336LEydOaD4g6aW8vDzOB6oS5wY9COcHVeXcuXMAarY4pWkGXwg6OztXefj35s2bAAAXFxe1565cuaLWdvnyZXh7e2s8H+kvzgeqDucGPQjnB1XnypUrePHFF+v0PQ2+EPTy8kJSUhLy8/PRoEEDVXtqaioAwNPTU+01vXr1QmxsLFxdXatcMSQiIiLSlMLCQly5cgW9evWq8/eWhBCizt+1DqWlpaFjx45YsmQJpk6dCuDvTxp57rnn4OjoiOTkZJkTEhEREcnD4FcE27dvj+DgYMyePRu3bt2Cu7s7YmJicO3aNWzYsEHueERERESyMfgVQeDvFcB58+YhNjYWubm5eOGFF/Dee+/x3oBERERk1IyiECQiIiIidQZ/H0EiIiIiqhoLwX8pLi7GzJkz4eLiAktLS3Ts2BGJiYlyxyItSUpKgkKhqPIrLS2tUt9z586hd+/eaNCgAezt7TF8+HBkZ2dXOe6nn36KZ555BhYWFvDw8MCKFSvqYnOolgoKCjB//nz07t0bjRo1gkKhQExMTJV9tTEPbt++jVGjRsHR0RHW1tbw9/dHRkaGxraPHk9N50d4eHiV+5JnnnmmynE5P/Rfeno6xo0bh9atW8Pa2hpPPvkkBg4ciAsXLqj11el9hyCVQYMGCVNTUzFjxgyxbt060alTJ2FqaiqOHj0qdzTSgsOHDwtJksSkSZPE5s2bK31lZ2er+v3222/CwcFBtGjRQnz88cciKipKNGrUSHh6eoqSkpJKY65evVpIkiSCg4PF+vXrxfDhw4UkSWLRokV1vXlUQ5cvXxaSJAlXV1fh5+cnJEkSMTExav20MQ/Ky8tFp06dhLW1tXj33XfFJ598Ilq3bi1sbGzEhQsXtLrdVDM1nR9hYWHC3NxcbV+yZ88etb6cH4ZhwIABwsXFRUycOFF8+umnYsGCBaJJkybC2tpa/PTTT6p+ur7vYCH4j9TUVCFJkli6dKmqraioSDz99NOiU6dOMiYjbblXCMbHxz+w3xtvvCGsrKzEb7/9pmpLTEwUkiSJtWvXqtru3r0r7O3tRUBAQKXXDx06VFhbW4vc3FzNbgBpRHFxsfjjjz+EEEIcP3682l/02pgHX3zxhdoczMrKEnZ2diI0NFRj20i1V9P5ERYWJho0aPDQ8Tg/DEdycrIoLS2t1HbhwgVhbm4uhg4dqmrT9X0HC8F/TJ8+XZiamor8/PxK7R988IGQJElkZmbKlIy05V4hGBcXJ+7cuaP2D/oeJycnMXDgQLX2li1biu7du6seJyQkCEmSxN69eyv1S0lJEZIkidjYWM1uAGlcenp6tb/otTEPgoODhbOzs9qYr7/+urCyslJbLSB5PWh+hIWFCWtra1FeXi7y8vKqHYPzw/C1bdtWtGvXTvVY1/cdPEfwHxkZGfDw8IC1tXWldh8fHwDAyZMn5YhFdWDEiBFo2LAhLCws4O/vjx9++EH13PXr15GVlYV27dqpvc7Hx6fS+Rj3vr+/b9u2baFQKDiH9Ji25kFGRgbatm1b5Zh3797F+fPnNbUJVAfu3r0LGxsb2Nrawt7eHuPGjUNBQUGlPpwfhk0IgT/++AMODg4A9GPfwULwHzdv3oSzs7Na+722qj6vmPSbmZkZgoKCEB0djV27dmHBggU4ffo0unTpovoHd+8zqaubGzk5OSgtLVX1NTExUe0A7qlfvz7s7e05h/SYtuYB9zuGw8XFBTNnzsTGjRuxdetWBAYGYuXKlejduzfKy8tV/Tg/DNvmzZtx48YNDBw4EIB+7DsM/pNFaqqwsBBmZmZq7ebm5qrnybD4+vrC19dX9VipVCIoKAht2rTB7NmzsXfvXtXf+8PmhqmpKQoLC1G/fv0q38vMzIxzSI9pax4UFRVxv2MgoqKiKj0OCQmBh4cH5s6di7i4OFVhwPlhuH7++WeMHTsWnTp1QlhYGAD92HdwRfAfFhYWKC4uVmsvKipSPU+Gz93dHYGBgTh8+DCEEKq/95rMDQsLC5SUlFQ5blFREeeQHtPWPOB+x7BNnjwZCoUCBw8eVLVxfhim33//Ha+88grs7OwQFxcHSZIA6Me+g4XgP5ydnatcSr23rOvi4lLXkUgmzZs3R0lJCQoKClTL7Pfmwb/dvHkT9vb2MDU1BfD3HCovL1e7N1RJSQlycnI4h/SYtuYB9zuGzdzcHI0aNUJOTo6qjfPD8OTl5aFPnz64c+cO9u3bhyZNmqie04d9BwvBf3h5eeH8+fPIz8+v1J6amgoA8PT0lCMWyeDXX3+FhYUFrK2t0bRpUzg6OiI9PV2tX1paWqV54eXlBQBqfY8fP46KigrOIT2mrXng6emJEydOQNz3SZ+pqamwsrKCh4eHJjeD6lh+fj6ys7Ph6OioauP8MCxFRUUICAjAxYsXsWfPHrRq1arS83qx76jx9cUG7t59BD/88ENV2737CPr6+sqYjLTl1q1bam0nT54Upqamom/fvqq2N954Q1haWlZ5D6g1a9ao2goLC3kfQT33oNuDaGMe3LsXWFxcnKotKytL2NraisGDB2ty00gDqpsfRUVF4s6dO2r9p0+fLiRJEjt27FC1cX4YjrKyMhEYGCjq16+vdsuXf9P1fYckxH3lpBEbOHAgtm/fjsmTJ8Pd3R0xMTE4fvw4Dh48iM6dO8sdjzTM398flpaW8PX1hZOTE86ePYu1a9fCzMwMKSkpaNmyJQAgMzMTXl5esLW1xcSJE5Gfn48lS5bgiSeeQHp6umpZHwBWrVqFsWPHIigoCD179sSRI0ewadMmREVFYdasWXJtKj3EihUrcPv2bdy4cQOrV69G//79Vf/7njBhAmxsbLQyDyoqKtC5c2f89NNPmD59Ouzt7bFy5UpkZmYiPT0dLVq0qPOfBal72PzIycmBl5cXQkNDVfuN/fv3Y+/evejTpw8SEhIqjcf5YRgmTZqE6OhoBAQEIDg4WO35oUOHAtDO7xCNzo1HKhsNXFFRkZg+fbpwdnYW5ubmokOHDuLAgQNyxyItiY6OFh06dBD29vbC1NRUNG3aVAwfPlxcunRJre+ZM2dEr169hJWVlWjUqJEYNmxYlSuKQgixbt060apVK2FmZiZatGgh/vvf/2p7U+gxubq6CkmShCRJQqFQCIVCofr+6tWrqn7amAe5ubkiMjJSODg4CCsrK+Hn5yd++OEHrWwn1c7D5sft27fFsGHDRIsWLYSVlZUwNzcXzz//vFi4cKEoKyurckzOD/3XrVs31Vy4/0uhUFTqq8v7Dq4IEhERERkpXixCREREZKRYCBIREREZKRaCREREREaKhSARERGRkWIhSERERGSkWAgSERERGSkWgkRERERGioUgERERkZFiIUhERERkpFgIEhERERkpFoJERFXYuHEjFAoFrl27ppHxkpKSoFAoVF8nTpzQyLiaYmtrq8o2fvx4ueMQUR1hIUhEOuf06dMICgqCq6srLCws0KxZM/Ts2RMrVqyo1C8qKgo7d+7UWg5JkjQ+5ty5cxEbG4unnnpK42MDwNmzZ/H222/j6tWrj/S69evXY9OmTQC0s91EpJskIYSQOwQR0T3Jycnw8/ODq6srwsLC0KRJE/z2229ISUnBr7/+ivPnz6v6WltbIyQkBJ999pnGc1RUVKCsrAz169fXyHhJSUnw9/dHUlISunbtqpExqxIXF4eQkJBav49CocC4ceMQHR2thXREpGvqyR2AiOjf3n//fdjZ2SE9PR02NjaVnsvOzq70WJIk1PT/sgUFBbCysqpxDoVCobEiUA41+bkIIVBcXAxzc/M6SEREuoiHholIp1y6dAmtW7dWKwIBwMHBQfW9QqFAQUEBYmJiVOe2jRw5EgDw9ttvQ6FQ4Ny5cwgNDUWjRo3QpUsXAMCpU6cQHh4ONzc3WFhYwNnZGREREcjJyan0XlWdI+jq6oqAgAAcPXoU7du3h4WFBdzd3VWHVGurW7du8PPzU2sPDw9XO4S8detWeHt7w8bGBg0bNkSbNm1Uq3cbN25ESEgIAMDPz0/1c/nuu+8q5d+/fz/atWsHS0tLrFmz5rGyE5F+44ogEekUV1dXpKSk4MyZM2jdunW1/TZt2oTIyEh06NABo0aNAgC4u7tX6hMcHAwPDw988MEHqhWyxMREXL58GREREWjSpAl++uknrF27FmfOnMGxY8cemE2SJFy8eBHBwcGIjIzEiBEj8OmnnyI8PBze3t549tlna7XNkiRVe17ev9u/+eYbhIaGonv37njttdcA/H1OYHJyMiZMmID//Oc/mDBhAqKjozF37lw888wzAKD6U5Ik/PLLLwgNDcXo0aPx+uuvo2XLlrXKTESGgYUgEemUadOmoU+fPvD09ISPjw+6dOmC7t27w8/PD/Xq/f8ua8iQIRg9ejTc3NwQGhpa5Vienp6IjY2t1DZmzBhMmTKlUlvHjh0xePBgHD16FJ07d642mxACv/zyC44cOYIXX3wRwN/FZvPmzbFhwwYsWbKkVtsshKi2EPz3Id6EhAQ0bNgQ+/fvr7L/U089hc6dOyM6Oho9evRQO0dQCIGLFy9i//796NGjR62yEpFh4aFhItIp3bt3R0pKCgIDA3Hq1CksWbIEvXr1QtOmTbF79+5HGmv06NFqbf8+H66oqAjZ2dno0KEDACAjI+OhY7Zu3VpVBAJ/H65u2bIlLl++/EjZasPOzg5//fUXDhw4UOsx3NzcWAQSkQoLQSLSOe3atUN8fDxu376NtLQ0zJ49G/n5+QgKCsK5c+dqPE5Vt2jJycnBxIkT0bhxY1haWsLJyQlubm4AgLy8vIeO+cQTT6i12draIjc3t8a5amvMmDHw8PBAnz590Lx5c0RERGD//v2PNIa2bltDRPqJhSAR6ax69eqhXbt2eP/997Fq1SqUlpbiyy+/rPHrLSws1NpCQkKwfv16jBkzBtu3b8c333yDffv2Afj7ljEPY2JiUmX749yJq7rDwuXl5ZUeOzo64uTJk9i1axcCAwNx+PBh9OnTB+Hh4TV+r6p+JkRkvHiOIBHpBW9vbwDA77//rmp71Bsf5+bm4tChQ3j33Xfx5ptvqtovXLigmZC1ZGdnV+Wh5atXr6pto6mpKZRKJZRKJYQQGDNmDNasWYO33noLbm5uvBk0ET0SrggSkU45fPhwle1ff/01AFS6ytXKyuqRDsneW827f+Xvo48+etSYGvX000/j559/rnSfxB9//BHff/99pX733+JGkiQ8//zzAIDi4mIAUN0rsS4OVROR/uOKIBHplPHjx6OwsBD9+vVDy5YtUVJSguTkZGzbtg1PPfUURowYoerr7e2NxMRELF++HM7OznBzc0P79u2rHdvGxgZdu3bF4sWLUVpaChcXFxw4cABXrlx57NyPc2h45MiRWLZsGXr16oWRI0fi1q1bWLNmDZ577jncuXNH1S8iIgK5ubnw9/dH06ZNcfXqVXz88cfw8vJS3SLGy8sLJiYmWLRoEW7fvg0zMzO89NJLcHR0fOxtJCLDwxVBItIpS5cuhZ+fH77++mtMnToVU6dOxfHjxzF27FikpqZWutH0smXL4O3tjTfffBOhoaFYvXq16rnqDpF+/vnn6NWrFz755BPMmTMHZmZm2Lt3tR4HhwAAAR9JREFUb5V97x/jQff6e5xDsq1atcL//vc/5OXlYerUqdizZw9iY2PRtm3bSuMOGzYM5ubmWLlyJcaOHYtNmzZh8ODBlfI3btwYq1evxq1btxAZGYkhQ4aoLrDhYWMiuh8/a5iIqA7c+6zhHTt2oFOnTrCzs6v2whM55OTkoKKiAk5OTvysYSIjwhVBIqI61LdvXzg5OeHHH3+UO0olbm5ucHJy4qohkZHhOYJERHXA09MTiYmJqsceHh4yplG3e/dulJaWAgCaN28ucxoiqis8NExERERkpHhomIiIiMhIsRAkIiIiMlIsBImIiIiMFAtBIiIiIiPFQpCIiIjISLEQJCIiIjJSLASJiIiIjBQLQSIiIiIjxUKQiIiIyEj9H949EwYoqt39AAAAAElFTkSuQmCC", + "text/plain": [ + "PyPlot.Figure(PyObject )" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "function plot_stress_strain_curve(ip)\n", + " steps = length(ip.fields[\"total strain\"])\n", + " eps11 = zeros(steps)\n", + " sig11 = zeros(steps)\n", + " for i=1:steps\n", + " # extract from integration points -- too complicated\n", + " # field -> timestep -> increment -> (vector of tensors, take first) -> (first component)\n", + " eps11[i] = ip.fields[\"total strain\"][i][end][1][1,1]\n", + " sig11[i] = ip.fields[\"stress\"][i][end][1][1,1]\n", + " end\n", + " eps11 *= 1.0e6\n", + " sig11 *= 1.0e-6\n", + "\n", + " PyPlot.figure(figsize=(7, 5))\n", + " PyPlot.plot(eps11, sig11, \"-ko\")\n", + " PyPlot.title(\"Stress-Strain curve\")\n", + " PyPlot.xlabel(\"Strain [ustr]\")\n", + " PyPlot.ylabel(\"Stress [MPa]\")\n", + " PyPlot.ylim([0, 250])\n", + " info(sig11')\n", + " info(eps11')\n", + "end\n", + "\n", + "plot_stress_strain_curve(ip)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 0.4.0", + "language": "julia", + "name": "julia-0.4" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "0.4.0" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/src/integrate.jl b/src/integrate.jl index 9571a20..5b8497f 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -1,6 +1,7 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + function get_default_integration_points(element::Quad4) [ IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0), diff --git a/src/types.jl b/src/types.jl index 77482a3..7e3dea4 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,28 +1,33 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -### INTEGRATIONPOINT ### +function ForwardDiff.derivative{T}(f::Function, S::Matrix{T}, args...) + shape = size(S) + wrapper(S::Vector) = f(reshape(S, shape)) + deriv = ForwardDiff.gradient(wrapper, vec(S), args...) + return reshape(deriv, shape) +end """ Integration point -xi :: Array{Float64, 1} +xi (dimensionless) coordinates of integration point -weight :: Float64 - Integration weight -attributes :: Dict{Any, Any} - This is used to save internal variables of IP needed e.g. for incremental - material models. +weight + integration weight +fields + FieldSet what can be used to store internal variables, stress, strain, ... """ type IntegrationPoint xi :: Vector weight :: Float64 - fields :: Dict{ASCIIString, FieldSet} + fields :: FieldSet end + function IntegrationPoint(xi, weight) IntegrationPoint(xi, weight, Dict()) end -call(N::Basis, ip::IntegrationPoint) = N(ip.xi) - - +function Base.convert(::Type{Number}, ip::IntegrationPoint) + return ip.xi +end