mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
Jukka Aho
0f0c49da62
A lot of old files from old documentation systems etc. is in package. These are now removed or moved. Old notebooks are in docs/tutorials. This PR closes issue #124.
921 lines
67 KiB
Plaintext
921 lines
67 KiB
Plaintext
{
|
||
"cells": [
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"#One element solution with and without nonlinear geometrics.\n",
|
||
"\n",
|
||
"Author: Jukka Aho\n",
|
||
"\n",
|
||
"Email: <jukka.aho@kapsi.fi>"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 1,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"0x00023c5cd89a9c57"
|
||
]
|
||
},
|
||
"execution_count": 1,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"tic()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"We consider here linear quadrangle element and evaluate solution for it. It has been verified using Elmer that the correct solution for minimum displacement in tip is $u_\\mathrm{min,nonlinear}$ = -2.22224475 for geometrically nonlinear case and $u_\\mathrm{min,linear} = -2.18056991$ for linear case.\n",
|
||
"\n",
|
||
"Solution domain is $\\Omega = [0, 10] \\times [0, 1] \\in \\mathbb{R}^2$ and we are having $F=-2 \\mathrm{N}$ nodal force in $(10, 1)$. $E=90$, $\\nu=0.25$. Left side is fixed, $u|_{x=0}=0$."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Geometry, material definitions and shape functions:"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 2,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"(36.0,24.0)"
|
||
]
|
||
},
|
||
"execution_count": 2,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"X = [[0.0 10.0 10.0 0.0],\n",
|
||
" [0.0 0.0 1.0 1.0]]\n",
|
||
"\n",
|
||
"# Partial derivatives of bilinear Lagrange polynomials\n",
|
||
"dNdξ(ξ) = [[-(1-ξ[2])/4.0 -(1-ξ[1])/4.0],\n",
|
||
" [ (1-ξ[2])/4.0 -(1+ξ[1])/4.0],\n",
|
||
" [ (1+ξ[2])/4.0 (1+ξ[1])/4.0],\n",
|
||
" [-(1+ξ[2])/4.0 (1-ξ[1])/4.0]] \n",
|
||
"\n",
|
||
"E = 90\n",
|
||
"ν = 0.25\n",
|
||
"μ = E/(2*(1+ν))\n",
|
||
"λ = E*ν/((1+ν)*(1-2*ν))\n",
|
||
"λ = 2*λ*μ/(λ + 2*μ)\n",
|
||
"μ, λ"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Internal force vector:"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 3,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"f_int (generic function with 2 methods)"
|
||
]
|
||
},
|
||
"execution_count": 3,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"function f_int(X, u, dNdξ, λ, μ, dim=2)\n",
|
||
" T = zeros(size(X))\n",
|
||
" I = eye(dim)\n",
|
||
" \n",
|
||
" function J(ξ)\n",
|
||
" Jᵀ = X*dNdξ(ξ)\n",
|
||
" ∇N = inv(Jᵀ)*dNdξ(ξ)'\n",
|
||
" ∇u = u*∇N'\n",
|
||
" F = I + ∇u # Deformation gradient\n",
|
||
" E = 1/2*(∇u' + ∇u + ∇u'*∇u) # Green-Lagrange strain tensor\n",
|
||
" P = λ*trace(E)*I + 2*μ*E # PK1 stress tensor\n",
|
||
" S = F*P # PK2 stress tensor\n",
|
||
" return S*∇N*det(Jᵀ)\n",
|
||
" end\n",
|
||
"\n",
|
||
" a = 1/sqrt(3)\n",
|
||
" ipoints = [[-a -a], [a -a], [a a], [-a a]]\n",
|
||
" iweights = [1 1 1 1]\n",
|
||
"\n",
|
||
" for m = 1:length(iweights)\n",
|
||
" w = iweights[m]\n",
|
||
" ξ = ipoints[m, :]\n",
|
||
" T += w*J(ξ)\n",
|
||
" end\n",
|
||
"\n",
|
||
" return T\n",
|
||
"\n",
|
||
"end"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Solution"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 4,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"Iter f(x) inf-norm Step 2-norm \n",
|
||
"------ -------------- --------------\n",
|
||
" 0 2.000000e+00 NaN\n",
|
||
" 1 1.000142e+00 2.257867e-03\n",
|
||
" 2 9.933192e-01 4.422932e-03\n",
|
||
" 3 9.961075e-01 8.847632e-03\n",
|
||
" 4 9.855335e-01 1.769324e-02\n",
|
||
" 5 9.738205e-01 3.538698e-02\n",
|
||
" 6 9.612786e-01 7.077382e-02\n",
|
||
" 7 9.404424e-01 1.415511e-01\n",
|
||
" 8 9.493268e-01 2.831322e-01\n",
|
||
" 9 1.157757e+00 5.665408e-01\n",
|
||
" 10 6.914588e-01 2.840319e-01\n",
|
||
" 11 6.914588e-01 0.000000e+00\n",
|
||
" 12 5.805109e-01 2.845511e-01\n",
|
||
" 13 5.805109e-01 0.000000e+00\n",
|
||
" 14 4.993552e-01 2.851202e-01\n",
|
||
" 15 4.993552e-01 0.000000e+00\n",
|
||
" 16 4.100475e-01 2.858343e-01\n",
|
||
" 17 4.100475e-01 0.000000e+00\n",
|
||
" 18 3.195632e-01 2.866934e-01\n",
|
||
" 19 3.195632e-01 0.000000e+00\n",
|
||
" 20 2.284246e-01 2.877431e-01\n",
|
||
" 21 1.444984e-01 3.012556e-01\n",
|
||
" 22 8.015667e-05 4.216354e-03\n",
|
||
" 23 5.099410e-10 1.611875e-05\n",
|
||
"Results of Nonlinear Solver Algorithm\n",
|
||
" * Algorithm: Trust-region with dogleg and autoscaling\n",
|
||
" * Starting Point: [0.0,0.0,0.0,0.0]\n",
|
||
" * Zero: [-0.39914506093725427,-2.1779892316773526,-0.0722858269436904,-2.2222447543690937]\n",
|
||
" * Inf-norm of residuals: 0.000000\n",
|
||
" * Iterations: 23\n",
|
||
" * Convergence: true\n",
|
||
" * |x - x'| < 0.0e+00: false\n",
|
||
" * |f(x)| < 1.0e-08: true\n",
|
||
" * Function Calls (f): 24\n",
|
||
" * Jacobian Calls (df/dx): 19\n"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"2x4 Array{Float64,2}:\n",
|
||
" 0.0 -0.399145 -0.0722858 0.0\n",
|
||
" 0.0 -2.17799 -2.22224 0.0"
|
||
]
|
||
},
|
||
"execution_count": 4,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"using NLsolve\n",
|
||
"\n",
|
||
"function solve()\n",
|
||
"\n",
|
||
" free_dofs = [3, 4, 5, 6]\n",
|
||
" u = zeros(2, 4)\n",
|
||
" F = zeros(2, 4)\n",
|
||
" F[2, 3] = -2.0\n",
|
||
"\n",
|
||
" function f(u_)\n",
|
||
" u[free_dofs] = u_\n",
|
||
" T = f_int(X, u, dNdξ, λ, μ)\n",
|
||
" R = T-F\n",
|
||
" return R[free_dofs]\n",
|
||
" end\n",
|
||
"\n",
|
||
" # Go!\n",
|
||
" sol = nlsolve(not_in_place(f), zeros(4), store_trace=true, show_trace=true)\n",
|
||
" println(sol)\n",
|
||
" u[free_dofs] = sol.zero\n",
|
||
" return u\n",
|
||
"end\n",
|
||
"\n",
|
||
"u = solve()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Difference to Elmer solution is"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 5,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"4.3690935314089074e-9"
|
||
]
|
||
},
|
||
"execution_count": 5,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"abs(u[2, 3] - -2.22224475)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Linear solution is achieved by assuming small deformation, i.e. :\n",
|
||
" \n",
|
||
" F = I + ∇u ≈ I\n",
|
||
" E = 1/2*(∇u' + ∇u + ∇u'*∇u) ≈ 1/2*(∇u' + ∇u)\n"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"**OFFTOPIC:** Here we should put something like\n",
|
||
"\n",
|
||
" using JuliaFEM\n",
|
||
" test_allclose(u[2, 3], -2.22224457, tol = ...)\n",
|
||
"\n",
|
||
"and get a nice green \"**PASSED**\" if everything is ok at this point, otherwise assertionerror"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Notes about linearization"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"For more robust solution, we shoud get analytical Jacobian for nonlinear system of equations. If not given it's probably some sort of finite difference approximation, i.e."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 6,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"8x8 Array{Float64,2}:\n",
|
||
" 119.461 13.4859 61.8587 … -26.0246 -113.693 2.80636\n",
|
||
" 13.486 309.575 15.6443 -148.518 -3.10971 -316.139 \n",
|
||
" 61.8587 15.6443 130.663 -53.2775 -60.0524 1.5741 \n",
|
||
" 9.73232 155.082 36.0592 -303.876 1.51268 -163.427 \n",
|
||
" -67.6267 -26.0207 -132.47 63.0313 56.5391 10.2937 \n",
|
||
" -26.0246 -148.518 -53.2777 … 299.44 16.2711 152.956 \n",
|
||
" -113.693 -3.10985 -60.0524 16.2711 117.207 -14.6739 \n",
|
||
" 2.80624 -316.14 1.57406 152.956 -14.674 326.611 "
|
||
]
|
||
},
|
||
"execution_count": 6,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"function Lin(f, h=1.0e-6)\n",
|
||
"\n",
|
||
" function D(x)\n",
|
||
" J_approx = zeros((8, 8))\n",
|
||
" for i=1:8\n",
|
||
" Δx = zeros(2, 4)\n",
|
||
" Δx[i] += h\n",
|
||
" Δf = f(x+Δx) - f(x)\n",
|
||
" J_approx[i, :] = Δf[:]/h\n",
|
||
" end\n",
|
||
" return J_approx\n",
|
||
" end\n",
|
||
"\n",
|
||
" return D\n",
|
||
"\n",
|
||
"end\n",
|
||
"\n",
|
||
"J_approx = Lin(u -> f_int(X, u, dNdξ, λ, μ))\n",
|
||
"J_approx(u)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"A more elegant solution is to use automatic differentiation to get analytical tangent stiffness matrix"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 7,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"8x8 Array{Float64,2}:\n",
|
||
" 119.461 13.486 61.8587 … -26.0246 -113.693 2.80625\n",
|
||
" 13.486 309.576 15.6443 -148.518 -3.10973 -316.14 \n",
|
||
" 61.8587 15.6443 130.663 -53.2776 -60.0524 1.57406\n",
|
||
" 9.73233 155.082 36.0593 -303.877 1.51267 -163.427 \n",
|
||
" -67.6267 -26.0206 -132.47 63.0312 56.539 10.2937 \n",
|
||
" -26.0246 -148.518 -53.2776 … 299.439 16.271 152.956 \n",
|
||
" -113.693 -3.10973 -60.0524 16.271 117.207 -14.674 \n",
|
||
" 2.80625 -316.14 1.57406 152.956 -14.674 326.611 "
|
||
]
|
||
},
|
||
"execution_count": 7,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"using ForwardDiff\n",
|
||
"\n",
|
||
"function T!(u, T)\n",
|
||
" T[:] = f_int(X, reshape(u, 2, 4), dNdξ, λ, μ)\n",
|
||
"end\n",
|
||
"\n",
|
||
"J_analytical = forwarddiff_jacobian(T!, Float64, fadtype=:dual, n=8, m=8)\n",
|
||
"J_analytical(reshape(u, 8))"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Difference is in smooth situations quite small, but analytical is of course faster."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 8,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"0.000998738079715485"
|
||
]
|
||
},
|
||
"execution_count": 8,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"norm(J_approx(u) - J_analytical(reshape(u, 8)))"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"If we use this analytical Jacobian we should get quadratic convergence when using Newton-Rhapson method:"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 9,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"Iter f(x) inf-norm Step 2-norm \n",
|
||
"------ -------------- --------------\n",
|
||
" 0 2.000000e+00 NaN\n",
|
||
" 1 1.985767e+00 3.756746e-01\n",
|
||
" 2 1.933311e+00 3.652956e-01\n",
|
||
" 3 1.855451e+00 3.896874e-01\n",
|
||
" 4 1.749547e+00 4.460338e-01\n",
|
||
" 5 1.581001e+00 5.270850e-01\n",
|
||
" 6 1.331440e+00 8.104083e-01\n",
|
||
" 7 8.751161e-02 2.235404e-01\n",
|
||
" 8 2.570554e-04 1.112294e-02\n",
|
||
" 9 3.202487e-09 3.961907e-05\n",
|
||
"Results of Nonlinear Solver Algorithm\n",
|
||
" * Algorithm: Newton with line-search\n",
|
||
" * Starting Point: [0.0,0.0,0.0,0.0]\n",
|
||
" * Zero: [-0.3991450608223195,-2.1779892314495743,-0.072285826868806,-2.22224475412425]\n",
|
||
" * Inf-norm of residuals: 0.000000\n",
|
||
" * Iterations: 9\n",
|
||
" * Convergence: true\n",
|
||
" * |x - x'| < 0.0e+00: false\n",
|
||
" * |f(x)| < 1.0e-08: true\n",
|
||
" * Function Calls (f): 110\n",
|
||
" * Jacobian Calls (df/dx): 18\n"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"2x4 Array{Float64,2}:\n",
|
||
" 0.0 -0.399145 -0.0722858 0.0\n",
|
||
" 0.0 -2.17799 -2.22224 0.0"
|
||
]
|
||
},
|
||
"execution_count": 9,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"function solve2()\n",
|
||
"\n",
|
||
" free_dofs = [3, 4, 5, 6]\n",
|
||
" u = zeros(2, 4)\n",
|
||
" F = zeros(2, 4)\n",
|
||
" F[2, 3] = -2.0\n",
|
||
"\n",
|
||
" function f(u_)\n",
|
||
" u[free_dofs] = u_\n",
|
||
" T = f_int(X, u, dNdξ, λ, μ)\n",
|
||
" R = T-F\n",
|
||
" return R[free_dofs]\n",
|
||
" end\n",
|
||
"\n",
|
||
" function g(u_)\n",
|
||
" u[free_dofs] = u_\n",
|
||
" J = J_analytical(reshape(u, 8))\n",
|
||
" return J[free_dofs, free_dofs]\n",
|
||
" end\n",
|
||
"\n",
|
||
" # Go!\n",
|
||
" sol = nlsolve(not_in_place(f, g), zeros(4), method = :newton, store_trace=true, show_trace=true)\n",
|
||
" println(sol)\n",
|
||
" u[free_dofs] = sol.zero\n",
|
||
" return u\n",
|
||
"end\n",
|
||
"\n",
|
||
"u = solve2()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Clearly the last iterations are converging quadratic speeds. Own naive implementation of solver would be:"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 10,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"norm = 3.0900221367289094\n",
|
||
"norm = 0.3212131602153428\n",
|
||
"norm = 0.04043178194002758\n",
|
||
"norm = 0.0009291101052123099\n",
|
||
"norm = 1.5638899074064188e-7\n",
|
||
"norm = 1.036219989655546e-14\n",
|
||
"norm = 3.130358566071589e-16\n"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"2x4 Array{Float64,2}:\n",
|
||
" 0.0 -0.399145 -0.0722858 0.0\n",
|
||
" 0.0 -2.17799 -2.22224 0.0"
|
||
]
|
||
},
|
||
"execution_count": 10,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"function solve3()\n",
|
||
"\n",
|
||
" free_dofs = [3, 4, 5, 6]\n",
|
||
" u = zeros(2, 4)\n",
|
||
" F = zeros(2, 4)\n",
|
||
" F[2, 3] = -2.0\n",
|
||
"\n",
|
||
" function f(u_)\n",
|
||
" u[free_dofs] = u_\n",
|
||
" T = f_int(X, u, dNdξ, λ, μ)\n",
|
||
" R = T-F\n",
|
||
" return R[free_dofs]\n",
|
||
" end\n",
|
||
"\n",
|
||
" function g(u_)\n",
|
||
" u[free_dofs] = u_\n",
|
||
" J = J_analytical(reshape(u, 8))\n",
|
||
" return J[free_dofs, free_dofs]\n",
|
||
" end\n",
|
||
"\n",
|
||
" u2 = zeros(4)\n",
|
||
" for i=1:7\n",
|
||
" du = g(u2) \\ -f(u2)\n",
|
||
" u2 += du\n",
|
||
" println(\"norm = \",norm(du))\n",
|
||
" end\n",
|
||
" u[free_dofs] = u2\n",
|
||
" return u\n",
|
||
"end\n",
|
||
"\n",
|
||
"u = solve3()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"I don't know what NLsolve does because it takes a bit longer to converge?"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Comparison of results"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"At end we do compare results between linear and nonlinear systems and put all together."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 11,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"2x4 Array{Float64,2}:\n",
|
||
" 0.0 -10.0735 -8.83486 0.0\n",
|
||
" 0.0 -10.4167 -11.8171 0.0"
|
||
]
|
||
},
|
||
"execution_count": 11,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"function f_int(X, u, dNdξ, λ, μ, dim=2; nlgeom=true)\n",
|
||
" T = zeros(size(X))\n",
|
||
" I = eye(dim)\n",
|
||
"\n",
|
||
" function J(ξ)\n",
|
||
" Jᵀ = X*dNdξ(ξ)\n",
|
||
" ∇N = inv(Jᵀ)*dNdξ(ξ)'\n",
|
||
" ∇u = u*∇N'\n",
|
||
" F = I # Deformation gradient for small deformations\n",
|
||
" E = 1/2*(∇u' + ∇u) # Small strain tensor\n",
|
||
" if nlgeom # .. if nonlinear geometry\n",
|
||
" F += ∇u # .. add nonlinear part to deformation gradient, and\n",
|
||
" E += 1/2*(∇u'*∇u) # add nonlinear part to strain tensor => Green-Lagrange strain tensor\n",
|
||
" end\n",
|
||
" P = λ*trace(E)*I + 2*μ*E # PK1 stress tensor\n",
|
||
" S = F*P # PK2 stress tensor\n",
|
||
" return S*∇N*det(Jᵀ)\n",
|
||
" end\n",
|
||
"\n",
|
||
" a = 1/sqrt(3)\n",
|
||
" ipoints = [[-a -a], [a -a], [a a], [-a a]]\n",
|
||
" iweights = [1 1 1 1]\n",
|
||
"\n",
|
||
" for m = 1:length(iweights)\n",
|
||
" w = iweights[m]\n",
|
||
" ξ = ipoints[m, :]\n",
|
||
" T += w*J(ξ)\n",
|
||
" end\n",
|
||
"\n",
|
||
" return T\n",
|
||
"\n",
|
||
"end\n",
|
||
"\n",
|
||
"function Lin(f, h=1.0e-6)\n",
|
||
"\n",
|
||
" function D(x)\n",
|
||
" J_approx = zeros((8, 8))\n",
|
||
" for i=1:8\n",
|
||
" Δx = zeros(2, 4)\n",
|
||
" Δx[i] += h\n",
|
||
" Δf = f(x+Δx) - f(x)\n",
|
||
" J_approx[i, :] = Δf[:]/h\n",
|
||
" end\n",
|
||
" return J_approx\n",
|
||
" end\n",
|
||
"\n",
|
||
" return D\n",
|
||
"\n",
|
||
"end\n",
|
||
"\n",
|
||
"function solve4!(F, u; nlgeom=true, analytical_jacobian=false)\n",
|
||
"\n",
|
||
" free_dofs = [3, 4, 5, 6]\n",
|
||
"\n",
|
||
" # Tangent stiffness matrix\n",
|
||
" function T!(u, T)\n",
|
||
" T[:] = f_int(X, reshape(u, 2, 4), dNdξ, λ, μ, nlgeom=nlgeom)\n",
|
||
" end\n",
|
||
" if analytical_jacobian\n",
|
||
" Jacobian = ForwardDiff.forwarddiff_jacobian(T!, Float64, fadtype=:dual, n=8, m=8)\n",
|
||
" else\n",
|
||
" Jacobian = Lin(u -> f_int(X, u, dNdξ, λ, μ))\n",
|
||
" end\n",
|
||
"\n",
|
||
" # Assembly of residual vector R(u) = T(u) - F(u)\n",
|
||
" function f(u_)\n",
|
||
" u[free_dofs] = u_\n",
|
||
" T = f_int(X, u, dNdξ, λ, μ, nlgeom=nlgeom)\n",
|
||
" return (T-F)[free_dofs]\n",
|
||
" end\n",
|
||
"\n",
|
||
" # Derivative of residual vector w.r.t displacement field u, \"tangent stiffness\"\n",
|
||
" function g(u_)\n",
|
||
" u[free_dofs] = u_\n",
|
||
" if analytical_jacobian\n",
|
||
" J = Jacobian(reshape(u, 8))\n",
|
||
" else\n",
|
||
" J = Jacobian(u)\n",
|
||
" end\n",
|
||
" return J[free_dofs, free_dofs]\n",
|
||
" end\n",
|
||
"\n",
|
||
" # Go!\n",
|
||
" #sol = nlsolve(not_in_place(f, g), zeros(4), method=:newton) # not working as expected\n",
|
||
" #sol = nlsolve(not_in_place(f), zeros(4)) # working, but not quadratic convergence\n",
|
||
"\n",
|
||
" # working, very well :)\n",
|
||
" u2 = zeros(4)\n",
|
||
" for i=1:100\n",
|
||
" du = g(u2) \\ -f(u2)\n",
|
||
" u2 += du\n",
|
||
" if norm(du) < 1.0e-6\n",
|
||
" #println(\"norm = \",norm(du), \" converged in \", i, \" iterations\")\n",
|
||
" break\n",
|
||
" end\n",
|
||
" end\n",
|
||
"\n",
|
||
" #u0[free_dofs] = sol.zero\n",
|
||
" u0[free_dofs] = u2\n",
|
||
" return u0\n",
|
||
"end\n",
|
||
"\n",
|
||
"u0 = zeros(2, 4)\n",
|
||
"F = zeros(2, 4)\n",
|
||
"F[2, 3] = -20.0\n",
|
||
"solve4!(F, u0; nlgeom=true, analytical_jacobian=true)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 12,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"Solving 1 F = -0.000 0.000 0.000\n",
|
||
"Solving 2 F = -0.500 -0.545 -0.549\n",
|
||
"Solving 3 F = -1.000 -1.090 -1.103\n",
|
||
"Solving 4 F = -1.500 -1.635 -1.662\n",
|
||
"Solving 5 F = -2.000 -2.181 -2.222\n",
|
||
"Solving 6 F = -2.500 -2.727 -2.783\n",
|
||
"Solving 7 F = -3.000 -3.291 -3.343\n",
|
||
"Solving 8 F = -3.500 -3.822 -3.899\n",
|
||
"Solving 9 F = -4.000 -4.361 -4.449\n",
|
||
"Solving 10 F = -4.500 -4.906 -4.992\n",
|
||
"Solving 11 F = -5.000 -5.451 -5.526\n",
|
||
"Solving 12 F = -5.500 -5.997 -6.047\n",
|
||
"Solving 13 F = -6.000 -6.542 -6.554\n",
|
||
"Solving 14 F = -6.500 -7.087 -7.045\n",
|
||
"Solving 15 F = -7.000 -7.632 -7.516\n",
|
||
"Solving 16 F = -7.500 -8.177 -7.965\n",
|
||
"Solving 17 F = -8.000 -8.722 -8.388\n",
|
||
"Solving 18 F = -8.500 -9.267 -8.783\n",
|
||
"Solving 19 F = -9.000 -9.813 -9.145\n",
|
||
"Solving 20 F = -9.500 -10.358 -9.474\n",
|
||
"Solving 21 F = -10.000 -10.903 -9.767\n",
|
||
"Solving 22 F = -10.500 -11.448 -10.025\n",
|
||
"Solving 23 F = -11.000 -11.993 -10.250\n",
|
||
"Solving 24 F = -11.500 -12.538 -10.443\n",
|
||
"Solving 25 F = -12.000 -13.083 -10.610\n",
|
||
"Solving 26 F = -12.500 -13.629 -10.753\n",
|
||
"Solving 27 F = -13.000 -14.174 -10.878\n",
|
||
"Solving 28 F = -13.500 -14.719 -10.987\n",
|
||
"Solving 29 F = -14.000 -15.264 -11.084\n",
|
||
"Solving 30 F = -14.500 -15.809 -11.171\n",
|
||
"Solving 31 F = -15.000 -16.354 -11.250\n",
|
||
"Solving 32 F = -15.500 -16.899 -11.323\n",
|
||
"Solving 33 F = -16.000 -17.445 -11.390\n",
|
||
"Solving 34 F = -16.500 -17.990 -11.452\n",
|
||
"Solving 35 F = -17.000 -18.535 -11.512\n",
|
||
"Solving 36 F = -17.500 -19.080 -11.568\n",
|
||
"Solving 37 F = -18.000 -19.625 -11.621\n",
|
||
"Solving 38 F = -18.500 -20.170 -11.673\n",
|
||
"Solving 39 F = -19.000 -20.715 -11.723\n",
|
||
"Solving 40 F = -19.500 -21.261 -11.771\n",
|
||
"Solving 41 F = -20.000 -21.806 -11.817\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"F = zeros(2, 4)\n",
|
||
"N = 41\n",
|
||
"F[2, 3] = -20.0\n",
|
||
"uall = zeros(2, 4, 2, N)\n",
|
||
"u0 = zeros(2, 4)\n",
|
||
"t = linspace(0, 1, N)\n",
|
||
"for i = 1:N\n",
|
||
" @printf \"Solving %2i F = %8.3f \" i t[i]*F[2, 3]\n",
|
||
" u0 = solve4!(t[i]*F, u0; nlgeom=false, analytical_jacobian=false)\n",
|
||
" uall[:, :, 1, i] = u0\n",
|
||
" @printf \"%8.3f\" u0[2, 3]\n",
|
||
" u0 = solve4!(t[i]*F, u0; nlgeom=true, analytical_jacobian=true)\n",
|
||
" uall[:, :, 2, i] = u0\n",
|
||
" @printf \"%8.3f\" u0[2, 3]\n",
|
||
" println()\n",
|
||
"end"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 13,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stderr",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"INFO: Loading help data...\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"import PyPlot"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 14,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": [
|
||
"iVBORw0KGgoAAAANSUhEUgAAAroAAAItCAYAAADWntniAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XmcjXX/x/HXmTtr9qXSjahuUXclUlEpsoyloVIiN0ZptSWUiigR0t0y2hUSI2QP2RnUrbG0MLLMkO3O2MJYZ87vj+s2P2NmmLnmnOt7neu8n4/HPKbOueY679On6/i45nt9Lp/f7/cjIiIiIuIxEaYDiIiIiIgEgxpdEREREfEkNboiIiIi4klqdEVERETEk9ToioiIiIgnqdEVEREREU9SoysiIiIinqRGV0REREQ8SY2uiIiIiHiSGl0RERER8SRXNbqrV6+mS5cu3HDDDRQpUoSrrrqK1q1bs3nz5kzbbty4kcjISIoWLUrp0qVp3749ycnJBlKLiIiIiBv5/H6/33SIs1q1asWqVat4+OGHuemmm9izZw8xMTEcPXqUH374gRtuuAGAnTt3csstt1CyZEm6devGkSNHePvtt6lYsSL/+c9/yJcvn+F3IiIiIiKmuarRXbVqFbVq1eKSSy5Jf2zLli3ceOONtGrViq+++gqAZ599lrFjx5KQkED58uUBWLhwIQ0bNuSTTz6hc+fORvKLiIiIiHu4qtHNTs2aNYmIiGD16tUAXH755dSrV4/Y2NgM21WtWpUKFSowf/58EzFFRERExEVctUY3K36/n//+97+UKVMGgF27drFv3z5uvfXWTNvWqlWLtWvXOh1RRERERFzI9Y3u119/ze7du2ndujUAe/bsAaBcuXKZti1XrhwHDhzg9OnTjmYUEREREfdxdaObkJDAc889R506dejQoQMAx48fB6BAgQKZti9YsGCGbUREREQkfF1y8U3M2Lt3L82aNaNkyZJMnjwZn88HQKFChQA4efJkpp85ceJEhm3OlZyczLx586hUqVKWz4uIiIiIWcePHycpKYnGjRunL1vNC1c2uocPH6ZJkyb89ddfLF++nCuuuCL9ubNLFs4uYTjXnj17KF26dJbjxebNm0e7du2CF1pEREREAmLcuHE89thjed6P6xrdEydOcP/997NlyxYWLFhA1apVMzz/97//nbJly6ZPYDjXf/7zH6pXr57lfitVqgRY/+GqVasW8NziPq1atWLy5MmmY4hDQq3eH38MU6fCvHmmk4SmUKu35E1I1Hv1anjhBShfHt5/HwJwNjIcbdy4kXbt2qX3bXnlqkY3NTWV1q1b8+OPPzJ9+nRuv/32LLd76KGHGDNmDDt37swwR3fz5s288MILWf7M2eUK1apVo0aNGsF5A+Iq+fLlU63DSKjV2+eDihUhhCK7SqjVW/LG9fWeMAG6doV69WDyZCha1HSikBeoZaauanRfeOEFZs6cyf33309ycjLjxo3L8PzZpQcvv/wykyZNol69enTv3p0jR44wfPhwbrrpJqKjo01EFxe67rrrTEcQB4VavffsgXNWZUkuhVq9JW9cW2+/H0aMgN69oX17+Pxz0N1ZXcVVje769evx+XzMnDmTmTNnZnjO5/OlN7rly5dn6dKl9OzZk5deeokCBQrQvHlzRowYodv/ikhI2LsXbrrJdAoRsS01FXr2tJYpvPIKvPGG9asacRVXNbqLFy/O8bbXX389c+fODWIaEZHg2bMHGjc2nUJEbDlxAtq1sxbaf/QRPP206USSDVc1uiKB1Lx5c9MRxEGhVO+0NOuMrpYu2BdK9Za8c1W9DxyAli3hp5+sRjcqynQiuQBX3zBCJC9mzZplOoI4KJTqfeAAnDkDWdzgUXIolOoteeeaem/fDnfdBRs2wMKFanJDgM7oimcNGDDAdARxUCjV++wYcDW69oVSvSXvXFHv9euhSRMoWBBWroQqVUwnkhxQo3uelJQUEhISTMeQAFmzZo3pCAFTtWpVChcubDqGa7l69NB59u61vmvpgn2hVG/JO+P1XrAAHnzQam5nz4bLLzebR3JMje55EhISqFmzpukYIpnEx8eb/7CXgNAZXZEQMm4cREdDgwYwaRIUKWI6keSCGt1s6A5q4hZn7xIj3rFnDxQvbv0GVERcyu+HoUOhb1+r0f3kE83IDUFqdLOhO6iJhJZRo0bx+OOPm46RI3v2wJVXmk4R2kKp3pJ3jtc7NRW6d4eRI6FfPxg4UDNyQ5SmLoiIJ4TSeuw9e7RsIa9Cqd6Sd47W+/hxePhhaz7uJ5/A66+ryQ1hOqMrIp4wcuRI0xFybPduqFTJdIrQFkr1lrxzrN7791sjw9auhenTwU3ze8UWndEVEXGYzuiKuFBSEtx5J/z+OyxZoibXI9Toiog4yO+3zuhqja6Ii6xdC7VrW3dyWbUKbrvNdCIJEDW6IiIO+usvawmgzuiKuMT330PdulC+vHUjiGuvNZ1IAkiNrmSSlJREREQE0dHR6Y917NiRiIgIduzYYTCZSPaiQuRWnJqhGxihUm8JjKDVe+xYaNbManSXLIHLLgvO64gxanQlW75zrjL1+XwZ/l3Ebbp06WI6Qo7s3m1919KFvAmVektgBLzefj8MHgwdOkDHjtaFZ5deGtjXEFdQoys5MmTIEDZu3MiV+tNZXKpRo0amI+SIzugGRqjUWwIjoPVOTYVnn4VXXoEBA+DTT+ESDaHyKlVWcuSKK67giiuuMB1DJOTt3g3FiunkkYgRKSnQpg3Mng2ffw666Yjn6Yyu5EhWa3TPXcublJTEo48+SpkyZShUqBC1atVi9uzZ2e5vwoQJ1KtXjxIlSlCoUCGuv/563nzzTU6dOpVp22nTptGuXTuqVKlCkSJFKFKkCLfeeisffPABfr8/26yJiYl88MEH3HTTTRQuXJh69eoF5j+GSB5otJiIIcnJcN99sGABzJihJjdMqNGVHMtuje727du5/fbb2bFjBx06dKB169b8+uuvtGjRgiVLlmTavlOnTjz22GNs27aNhx9+mC5dulCqVCn69etHZGQkqampGbbv27cv69ato3bt2nTr1o0OHTpw9OhRunfvTocOHbLN2717d/r378/NN99Mjx49uOuuu/L0/sXdpk2bZjpCjuzerUY3EEKl3hIYea73tm1Qpw5s3WpddNa0aUByiftp6YLkWFZnTwGWLFnCwIED6devX/pjbdu2JTIykuHDh3PvvfemPz569GhGjx7Ngw8+yNdff02BAgXSnxs4cCADBw5k5MiRdOvWLf3x7777jsqVK2fKEh0dzdixY+nSpQu3ZTHzcO3ataxbt46rrrrK7luWEDJhwgRatmxpOsZF7d4NFSuaThH6QqXeEhh5qnd8vNXYFitmzci95prAhhNXU6MbKCkpkJDg/OtWrQqFCzv/uueoVKkSr776aobHGjVqRIUKFVi9enWGx9977z3y5cvHF198kaHJBXj11VeJiYnh66+/ztDont/kgnV2uVu3bowdO5bvv/8+y0a3T58+anLDyMSJE01HyJHdu+GOO0ynCH2hUm8JDNv1njsXWrWCf/4TZs6EsmUDG0xcT41uoCQkQM2azr9ufDzUqOH8656jevXqWS5rqFChAj/++GP6v6ekpLB+/XrKli3LO++8k+W+8ufPz8aNGzM8tn//foYPH853333Htm3bSElJyfD8rl27stxXVs2viEm6K5qIg778Ejp3ts7mxsYaPykkZqjRDZSqVa2m08TrGlaiRIksH7/kkktIS0tL//eDBw8CsG/fPl5//fVs93du03zo0CFq1apFUlISt99+Ox07dqRUqVJccsklHDx4kPfee4+TJ09muR9NiRC3OXzYuiva3/9uOomIh/n9MGgQ9O8PTz0FMTEaHxbGVPlAKVzY+JlVtytevDgANWrU4KeffsrRz3z++eckJSUxYMAA+vfvn+G5VatW8d5772X7s7rBhbjN2V8+6IyuSJCcOWPNyP3sM6vZffll0J8FYU1TF8QxRYoU4YYbbuDXX39NP7t7MVu2bAHgoYceyvTc0qVLA5pPQtu5t6x2K90VLXBCod4SODmq97Fj8MAD8MUX1rKFV15RkytqdMVZPXv25NSpU3Tq1InDhw9nev7gwYOsXbs2/d/PXoi2ePHiDNutXbuWIUOGBDeshJRQuFPW2UZX48XyLhTqLYFz0Xrv2wf168PixdbNIDp2dCSXuJ+WLkhQnT+SLDo6mvj4eD788EOuueYaGjduTIUKFThw4ACJiYksX76cTp068eGHHwLQvn17hg8fTo8ePVi8eDHXXnstmzdvZvbs2Tz00EPExsaaeFviQm3atDEd4aJ274bSpaFgQdNJQl8o1FsC54L13roVIiPhyBFYutTMheHiWmp0JUd8Pl+u17xm9zMxMTE0adKEjz/+mAULFnDo0CFKly5NxYoV6dOnD+3atUvftly5cixfvpyXXnqJuLg45s2bR7Vq1fjoo4+47777smx07WQVccKuXVq2IBJQq1dDs2ZQsqQ1IzeLcZQS3tToSiaVKlXKMC0B4Msvv+TLL7+86HbnOn+5wbmaNWtGs2bNcpSnWrVqTJ8+Pcvnsnr9rLKKuIFGi4kE0OzZ8MgjcPPN1i19y5QxnUhcSGt0RcQT4uLiTEe4KDW6gRMK9ZbAyVTvzz+HFi2gYUNYsEBNrmRLja6IeMKwYcNMR7goLV0InFCotwROer39fhgwwLoRxJNPwpQpuhGEXJCWLoiIJ7j9wsTUVNizRzeLCBS311sCKzY2Fk6fhqeftsaHDR4ML72k8WFyUWp0RcQTCrv8rM6ff1rNbvnyppN4g9vrLYFVOC3NWqowfz6MHQv/+pfpSBIi1OiKiDjg7F3RdEZXJJf++19o3hwSEuC776x1uSI5pEZXRMQBO3da39XoiuTC5s3WjNyUFFi+HKpXN51IQowuRhMRT+jdu7fpCBe0axfkywdly5pO4g1ur7cEwI8/Qp06kD8/ve+/X02u2KJGV0Q8oWLFiqYjXNDZiQsR+tQNCLfXW/Jo5kyoVw+uuw5WrKDijTeaTiQhSh+5IuIJXbt2NR3hgnbt0rKFQHJ7vSUPPvkEWra0lizMnw+lSqneYpsaXRERB6jRFbkIvx/69bNGiD37LEyaBIUKmU4lIU4Xo4mIOGDXLtBvX0Wycfq0dQOI0aNh6FDo3VszciUgdEZXRDwhISHBdIQL2rlTM3QDye31llw4cgTuvx++/tr66tMnU5OreotdanRFxBP69OljOkK2/voLjh7V0oVAcnO9JRf27oV774WVK2HOHGjbNsvNVG+xS42uhITRo0cTERHBmDFjMjxeqVIlKleubCiVuElMTIzpCNnSzSICz831lhzatAlq17aa3eXL4b77st1U9Ra71OhKSPGd9+ssn8+X6TEJT24eN3X2ZhEVKpjN4SVurrfkwKpV1ozcQoWsf7755gturnqLXboYTULaokWLTEcQuag//rC+X3ml2RwirjBtGrRpA7VqwfTpULKk6UTiYWp0JaRp2YKEgp074fLLIX9+00lEDPvoI+jSBR58EL76CgoWNJ1IPE5LFySTpKQkIiIiiI6OJikpiUcffZQyZcpQqFAhatWqxezZszP9zMmTJ3nrrbe48cYbufTSSylevDh169Zl0qRJAdl/drJao3vuet7Fixdz7733UqxYMYoXL07z5s2zvXo3JSWFIUOGUL16dYoUKULRokWpU6cOsbGxmbY9ffo0MTExNG3alKuuuoqCBQtSunRpGjZsyNy5cy+Y9ciRI/Ts2ZNKlSqRP39+Bg4cmOP3K9kbOnSo6QjZ+uMPLVsINDfXW7Lg98PLL1vzcbt2hYkTc9Xkqt5il87oSra2b9/O7bffzjXXXEOHDh3Yv38/EydOpEWLFixYsIB7770XgFOnTtG4cWOWLVtGtWrV6NKlC8eOHWPy5Mm0bt2adevW8eabb9re/8Vkt0Z31qxZTJ8+naZNm/LMM8/w22+/8d1337F69Wo2bNhA6dKl07c9dOgQ9evXZ926ddSsWZPHH3+ctLQ05s6dS9u2bfntt99444030rffv38/PXr04M4776Rx48aULVuW3bt3M3PmTJo2bcpnn33G448/ninnqVOnqFevHocOHSIyMpJixYpx9dVX5+h9yoWlpKSYjpAtjRYLPDfXW85z6hQ88YR1BnfECOjZM9e7UL3FNn+YiI+P9wP++Pj4gGznZYmJiX6fz+f3+Xz+119/PcNz8+bN8/t8Pn/Tpk3THxs8eLDf5/P5mzVr5k9NTU1//M8///RXqlTJ7/P5/CtXrrS9f7/f7//yyy/9Pp/PP2bMmAyPX3XVVf7KlStnuW2+fPn8ixYtyvBc3759/T6fzz9s2LAMj3fo0MHv8/n8w4cPz/D4iRMn/JGRkf6IiAj/unXr0h8/efKkf9euXf7zHT582P/Pf/7TX6pUKf/x48czZfX5fP6GDRv6U1JSMv1sdvT/ZOj75z/9/q5dTacQMeDwYb+/QQO/P18+v3/CBNNpJAQE+s88ndENkJQUMDHPumpVKFw4OPuuVKkSr776aobHGjVqRIUKFVi9enX6Y1988QURERG88847RET8/2qYsmXL0q9fP5544gk+//xzateubWv/dj366KPUq1cvw2NPPvkkb731Vob979+/n3HjxlGrVi169eqVYfsCBQrw1ltvMW/ePMaPH8/N/7syOH/+/FyZxZVFxYoVIzo6ml69erF69WruvvvuDM/7fD5GjBhBId3WMqz88YfO6EoY2rMHmjaFxESYNw/O+zwWcYIa3QBJSICaNZ1/3fh4qFEjOPuuXr16lssCKlSowI8//gjAkSNH2Lp1K+XLl6dKlSqZtq1fvz4A69ats7X/vLj11lszPVb+f93GwYMH0x9bvXo1aWlpAAwYMCDTz5w+fRqAjRs3Znj8t99+Y/jw4Sxbtoy9e/dy4sSJDM/v3r07074KFizIjboPbFg5cgQOH1ajK2Fm40Zo0gTOnLFm5OpzTwxRoxsgVataTaeJ1w2WEiVKZPn4JZdckt4YHj58GIBy5cplue0VV1wBWGtg7ew/L7La/yWXWP/Lp6ampj+2f/9+wGp4szuT7PP5OHbsWPq///DDD9SvX5+0tDTuu+8+WrZsSbFixYiIiGDt2rVMnz6dkydPZtrPZZddlqf3JNlLTk6mTJkypmNkohm6weHWeguwYoV1S98rr7TudhaA//lVb7FLjW6AFC4cvDOrbla8eHEA9u7dm+Xze/bsybCdG53N1rNnT95+++0c/cygQYM4ceIES5YsoW7duhmeGzJkCNOnT8/y53Rzi+Dp1KkTM2bMMB0jk7ONrs7oBpZb6x32vv3Wuo3vHXdY83KzOaGRW6q32KXxYpInRYsW5ZprrmHnzp1s2bIl0/OLFy8GoIaL/xZw++23ExERwbJly3L8M1u2bKF06dKZmlyApUuXBjKe5FBWy07c4OzNInT738Bya73DWkwMtGoFLVpYa3ID1OSC6i32qdGVPOvUqRN+v5/evXtnWHKQnJzMG2+8gc/no1OnTgYTXljZsmV57LHH+Omnnxg0aFCWyya2bt1KUlJS+r9XrlyZ/fv388svv2TYbtSoUXz//ffBjixZcOtfpnSziOBwa73DUloavPiiNR/3+edhwgQoUCCgL6F6i11auiC2+P3+9H/u1asXc+bMYfr06dx88800adKElJQUJk2aRHJyMn369KFOnTq29++EmJgYNm/eTP/+/fnqq6+48847ufzyy9m9ezcbN27kp59+IjY2lkqVKgHQo0cP5s2bx1133cUjjzxCsWLF+Omnn1ixYgWtWrVi8uTJjuYX99qxAypWNJ1CJEhOnYLoaKu5/fe/oUcP04lEMlCjK7nm8/kyrDXNly8f8+fP55133mH8+PHExMSQL18+br75Zt5//31at26dp/1n99jZx3Py8xdTtGhRli5dyqeffsr48eP59ttvOXHiBFdccQX/+Mc/ePfdd2nQoEH69o0bN2bmzJkMGjSIiRMncskll3DbbbexZMkStm7dypQpU3KUVbxPd0UTzzp8GB56yJqqMHEiPPyw6UQimfj8Tp86M2TNmjXUrFmT+Pj4C/4KJKfbiThF/0/mzKhRozLdjc4NqlWDyEjrZJcEjlvrHTZ27bJm5O7YYV10ds89QX051Tt8BPrPPK3RFRFPWLNmjekImfj9WroQLG6sd9jYsAFq14aDByEuLuhNLqjeYp8aXRHxhJEjR5qOkMnBg9ZdE7V0IfDcWO+wsGwZ3HmnNVFh1Sq44QZHXlb1FrvU6IqIBMmOHdZ3ndEVT5g0CRo2hFtusdblamaehAA1uiIiQaJGVzzjvfegdWvr4rM5c8DFNwESOZcaXRGRIPnjD8iXD3TnZwlZaWnQq5c1NqxXLxg3LuAzckWCSY2uiHhCVFSU6QiZ7Nhhrc+N0CdtwLmx3p5z8iQ89hi88w68/z4MG2bsf2bVW+zSHF0R8YQuXbqYjpCJJi4Ejxvr7SmHDsEDD1gXnE2aZC1ZMEj1FrvU6IqIJzRq1Mh0hEz++AOuvtp0Cm9yY709Y+dOaNLEmpW7YAHcdZfpRKq32KZfqImIBMn27RotJiHm11+tGbl//QUrVriiyRXJCzW6IiJBcPo07N4NlSqZTiKSQ0uWWI1t6dLWkoVq1UwnEskzLV3IxsaNG01HEAH0/2JOTZs2jZYtW5qOkW7nTuuC9auuMp3Em9xW75A3cSK0bw9168KUKVCsmOlEGajeYpca3Wy0a9fOdAQRyYUJEya46g/C7dut7zqjGxxuq3dIe+cdeOEFaNcORo2C/PlNJ8pE9Ra71Oiep2rVqsTHx5uOIZJJ1apVTUdwtYkTJ5qOkEFSkvVdUxeCw231DklpaVaD++678NJLMHgw+HymU2VJ9Ra71Oiep3DhwtSoUcN0DBEJcdu3w+WXQ8GCppOIZOHECWupwuTJMHIkPPus6UQiQaFGV0QkCLZv17IFcamDB6FlS/jPf6z1uA88YDqRSNCo0RURCYKkJF2IJi70xx8QGQl798LChVCnjulEIkGl8WLiWdHR0aYjiIPcVu/t29XoBpPb6h0Sfv4Z7rgDUlJg5cqQanJVb7FLja54lu6kE17cVO/UVOvEmZYuBI+b6h0SFi2Cu++2Fo6vWgXXXWc6Ua6o3mKXGl3xrDZt2piOIA5yU7337LFuGKEzusHjpnq73vjx1nKFO+6ApUvhiitMJ8o11VvsUqMrIhJgZ2foqtEVo/x+GD4cHnsM2rSBWbOgaFHTqUQcpUZXRCTAzs7QVaMrxqSmQvfu0KcPvPIKjB4N+fKZTiXiODW64llxcXGmI4iD3FTvxEQoXVonz4LJTfV2nePH4ZFHrPm4H30Egwa59kYQOaV6i11qdMWzhg0bZjqCOMhN9U5MhKuvNp3C29xUb1c5cAAaNYI5c2DqVHj6adOJAkL1Frs0R1c8KzY21nQEcZCb6p2YCJUrm07hbW6qt2ts3w5NmsCff1pTFu64w3SigFG9xS6d0RXPKly4sOkI4iA31XvbNjW6weamervCunVQu7Z1a9+VKz3V5ILqLfap0RURCaDTp60Zump0xTHz50PdunDlldaM3CpVTCcScQ01uiIiAfTHH5CWpjW64pCvvoKmTeGuu2DJEuuGECKSTo2ueFbv3r1NRxAHuaXeiYnWd53RDS631NsYvx/eegvat4d//QumT4ciRUynCpqwr7fYpkZXPKtixYqmI4iD3FLvxERrkpNL4niWW+ptRGoqdOkCfftC//4wapTnZ+SGdb0lTzR1QTyra9eupiOIg9xS78REKF8e8uc3ncTb3FJvxx0/Dm3bwsyZ8Omn0Lmz6USOCNt6S56p0RURCaBt27Q+V4Jk/364/35Yv95aqtCsmelEIq6nRldEJIASE6FaNdMpxHMSE60Zufv3w+LFcNttphOJhASt0RXPSkhIMB1BHOSWeutmEc5wS70dsWaNNSP3zBlrfFgYNrlhVW8JKDW64ll9+vQxHUEc5IZ6Hzli3ZTqmmtMJ/E+N9TbEfPmwT33WFc3rlwJ115rOpERYVNvCTg1uuJZMTExpiOIg9xQ761bre9h2os4yg31DroxY6B5c6vRXbwYLrvMdCJjwqLeEhRqdMWzNI4mvLih3lu2WN91Rjf43FDvoPH74c03oWNHiI6GadPg0ktNpzLK0/WWoFKjKyISIFu3QvHiULq06SQSss6cgWeegVdfhYED4ZNP4BJdNy5il44eEZEA2bLFOpvr85lOIiEpJQUefRS++866CUSnTqYTiYQ8ndEVzxo6dKjpCOIgN9R761atz3WKG+odUPv2Qf36sHChdTMINbkZeK7e4hid0RXPSklJMR1BHOSGem/ZAnfcYTpFeHBDvQNm2zaIjITDh2HpUrj1VtOJXMdT9RZH6YyueNbAgQNNRxAHma73iROwc6fO6DrFdL0D5qefrBm5YM3IVZObJc/UWxynRldEJAASE62L5TVxQXJszhy4917rDiMrVuje0SJBoEZXRCQANENXcuWLL+D+++G++2DRIihb1nQiEU9SoyuelZycbDqCOMh0vbdsgYIFoVw5ozHChul62+b3w+uvw+OPwxNPwJQpULiw6VSuF7L1FuPU6IpnddJVy2HFdL23bLF+8xyhT1VHmK63LWfOwFNPwWuvwaBB8NFHmpGbQyFZb3EFHWHiWQMGDDAdQRxkut6//w7XXWc0QlgxXe9cO3YMWreGefPgyy+tu55JjoVcvcU1dO5BPKtGjRqmI4iDTNf799+hShWjEcKK6Xrnyp9/Qr161uiwWbPU5NoQUvUWV9EZXRGRPDp+HHbs0BldycKWLdaM3KNHrUZXDZuIo3RGV0Qkj7Zuta4x0hldyeA//4E6dax1uKtWqckVMcB1je6xY8d47bXXiIyMpFSpUkRERDBmzJhM23Xs2JGIiIhMX9WqVTOQWtxo1KhRpiOIg0zWe9Mm67saXee4/vieNctarnDttdaM3MqVTScKaa6vt7iW6xrdffv28cYbb7Bp0yaqV68OgM/ny3LbAgUKMG7cuAxfb7/9tpNxxcXWrFljOoK8VaT7AAAgAElEQVQ4yGS9f/8dSpSAMmWMRQg7rj6+P/sMWrSARo1g4UIoXdp0opDn6nqLq7luje6VV17J3r17ueyyy4iPj6dWrVrZbpsvXz7atm3rYDoJJSNHjjQdQRxkst5nJy5k83dyCQJXHt9+PwwYYM3JffZZeP99+NvfTKfyBFfWW0KC687o5s+fn8suuwwAv99/wW39fj9paWn89ddfTkQTEcnSpk1athD2Tp+2bgDx+uswZAjExKjJFXEB1zW6uZGSkkKxYsUoUaIEpUuXpkuXLhw7dsx0LBEJMxotFuaOHrWWKowda3299JJO74u4hOuWLuTUlVdeyYsvvkiNGjVIS0tjzpw5fPjhh6xfv54lS5bwN/1NWkQcsH+/9aXRYmHqv/+FZs2sv+189x00bGg6kYicI2TP6A4ePJjBgwfTqlUrHnnkEb788kvefPNNVqxYweTJk03HExeIiooyHUEcZKremzdb33VG11muOL43b4batWHXLli2TE1uELmi3hKSQrbRzcrzzz9PREQECxcuNB1FXKBLly6mI4iDTNX77Gixa6818vJhy/jx/cMPVpNboID1z/+bEiTBYbzeErI81egWLFiQUqVKceDAgWy3adq0KVFRURm+ateuzbRp0zJs9/3332f5N8jnnnsu0zy/NWvWEBUVRXJycobHX3vtNYYOHZrhsR07dhAVFUVCQkKGxz/44AN69+6d4bGUlBSioqKIi4vL8PiECROIjo7OlK1169Z6H+e8j/NvGRmq78Mr9Qj2+2jUqJGR9/Hpp69RsuRQLr00MO/jrFCvR7DfR5kyZcy9jxkzoH59qFqV1lWrMm3tWtvvwyv1CPb7KFy4sCfeh1fqEaj3MWHChPRerHLlylSvXp3nn38+037ywue/2GgDg3766Sduu+02Ro8eTfv27S+6/ZEjRyhevDhPPfUUH330UYbn1qxZQ82aNYmPj9c9s0UkYO6/H9LSYPZs00nEEZ98Yo0Oa9kSxo2DQoVMJxLxlED3ayF5RvfkyZMcOXIk0+NvvPEGAJGRkU5HEpEwtXEj6IaMYcDvh1dfhaefhueeg2++UZMrEgJcOXUhJiaGQ4cOsXv3bgBmzJjBjh07AOjWrRsHDhzglltuoW3btlz3v0ud582bx5w5c2jSpAktWrQwll3cY9q0abRs2dJ0DHGIiXofPw7btsH11zv6soLD9T59Gjp3hjFjYNgw6NVL48Mcps9zscuVje6IESPYvn07YN3+d+rUqXz77bf4fD7at29PyZIluf/++5k/fz5jxowhNTWVf/zjHwwZMoRevXoZTi9uMWHCBH0whhET9d60yTrRp0bXeY7V+8gRaNUKFi+Gr78G3Y3TCH2ei12ubHQTExMvus3YsWMdSCKhbOLEiaYjiINM1HvDBuu7li44z5F6790LTZvC1q0wd651AZoYoc9zscuVja6ISCjYsAH+/ncoXtx0Egm4TZsgMhJOnYLly+Gmm0wnEhEbQvJiNBERN9CFaB61ciXUqQOFC8OqVWpyRUKYGl0REZs2bND6XM+ZNg3uuw/++U+Ii4OKFU0nEpE8UKMrnpXV8GrxLqfrfeqUdQdYNbpmBKXeH34IDz1kDUeeNw9Klgz8a4gt+jwXu9ToimedvVOWhAen6715M6SmqtE1JaD19vuhb19rPm63bhAbCwULBm7/kmf6PBe7dDGaeFabNm1MRxAHOV3v336zvmuNrhkBq/epU/D449ZdzkaMgJ49A7NfCSh9notdanRFRGz45RcoVw7KlDGdRGz76y9rqcKyZTBhAjz6qOlEIhJganRFRGz4+WddjB/Sdu+2ZuQmJVnrce+913QiEQkCrdEVz4qLizMdQRzkdL3V6JqVp3pv3Ai1a0NysjVZQU2u6+nzXOxSoyueNWzYMNMRxEFO1vvwYetEoBpdc2zXOy4O7rwTihWDH36wxoiJ6+nzXOxSoyueFRsbazqCOMjJev/6q/Vdja45tuo9ZQo0aAA332zd7ax8+cAHk6DQ57nYpUZXPKtw4cKmI4iDnKz3zz/DJZdA1aqOvaScJ9f1/uADePhhaNkS5s6FEiWCE0yCQp/nYpcaXRGRXPr5Z2usWP78ppPIRaWlQZ8+1nzcnj1h/HgoUMB0KhFxiKYuiIjkki5ECxEnT0KnTtbosHffhe7dTScSEYfpjK54Vu/evU1HEAc5Ve+0NGuGrhpdsy5a78OHrfFhU6bAxIlqckOcPs/FLp3RFc+qWLGi6QjiIKfqvX07HDmiRte0C9Z71y6ryd2xA77/HurWdS6YBIU+z8UuNbriWV27djUdQRzkVL1/+cX6rkbXrGzr/dtv0KSJ9c9xcXDDDc6FkqDR57nYpaULIiK5sG4dlC5t3f5XXGbpUrjrLihZElatUpMrImp0RURyIz4eatQAn890Esngm2+gUSOrOMuWwd//bjqRiLiAGl3xrISEBNMRxEFO1Ts+HmrWdOSl5AIy1Pvdd+HRR6FVK5gzB4oXNxdMgkKf52KXGl3xrD59+piOIA5yot5791rXOanRNa9Pnz7WCIwXXoDnn4feveGrrzTc2KP0eS526WI08ayYmBjTEcRBTtQ7Pt76fuutQX8puYiYESOgbVtrycIHH0CXLqYjSRDp81zsUqMrnqVxNOHFiXrHx0OpUnDVVUF/KbmQQ4eo2Lkz/PADTJ4MDz5oOpEEmT7PxS41uiIiOXR2fa4uRDPojz+s8WG7d8OCBdaUBRGRbGiNrohIDv30k9bnGvXLL1C7Nhw9CitWqMkVkYtSoyueNXToUNMRxEHBrvfevdZJRK3PNWTxYrj7bihbFlatYuiMGaYTiYP0eS52qdEVz0pJSTEdQRwU7HqfvRBNZ3QNiI2FyEi47TbrphDlyun4DjOqt9ilRlc8a+DAgaYjiIOCXe+fftKFaI7z+2HECGjTBlq3hlmzoFgxQMd3uFG9xS41uiIiOaAL0RyWlmbNx+3VC/r2hTFjNCNXRHJNja6IyEX4/bB6tdbnOubECesM7gcfwMiRMHiw/oYhIrao0RXPSk5ONh1BHBTMem/fbl2MVrt20F5Czjp4EBo1spYpTJkCzz6b5WY6vsOL6i12qdEVz+rUqZPpCOKgYNZ75UrruxrdINuxwxoZ9ttvsHAhtGyZ7aY6vsOL6i126YYR4lkDBgwwHUEcFMx6r1oFVapAmTJBewlZvx6aNrXW4a5cCdddd8HNdXyHF9Vb7NIZXfGsGjVqmI4gDgpmvVeu1NncoFq40JqRe8UV1t8qLtLkgo7vcKN6i11qdEVELuDoUetkY506ppN41Pjx1i1969SBJUusZldEJEDU6IqIXMDq1ZCaqkY34Px+GDYMHnsM2raFmTOhaFHTqUTEY9ToimeNGjXKdARxULDqvXKldY+C668Pyu7DU2oqdO8OL74Ir74KX34J+fLlahc6vsOL6i12qdEVz1qzZo3pCOKgYNX77PrcCH1aBsbx4/DII9Z83I8/hjfesDUjV8d3eFG9xS59dItnjRw50nQEcVAw6p2WZl0bpQvRAuTAAWjYEObMgWnT4KmnbO9Kx3d4Ub3FLo0XExHJxu+/W/cv0PrcAEhKsi46S06GxYvh9ttNJxKRMKAzuiIi2Vi50vqtunqyPFq3zjotfuqU9R9V/0FFxCFqdEVEshEXBzfeaF2MJjbNn2/NyC1f3mpy//EP04lEJIyo0RXPioqKMh1BHBSMei9eDPXqBXy34eOrr6y7nd19t/Uf8/LLA7ZrHd/hRfUWu9Toimd16dLFdARxUKDrnZhoLSutXz+guw0Pfj8MGQLt21tf06dDkSIBfQkd3+FF9Ra71OiKZzVq1Mh0BHFQoOu9aJE1Uqxu3YDu1vtSU+G55+Dll+G11+Dzz3M9IzcndHyHF9Vb7NLUBRGRLCxaBDVrQokSppOEkOPH//8uZ599Bk88YTqRiIQ5NboiIufx+61Gt0MH00lCSHIyREXB+vXWUoVmzUwnEhHR0gXxrmnTppmOIA4KZL0TEmDvXq3PzbHERLjzTtiyBZYscaTJ1fEdXlRvsUuNrnjWhAkTTEcQBwWy3osWWctK77wzYLv0rjVrrBm5aWnW+LBatRx5WR3f4UX1FrvU6IpnTZw40XQEcVAg671oEdxxB1x6acB26U3z5sE998BVV1lN7rXXOvbSOr7Di+otdqnRFRE5R1qaNfJVyxYuYvRoaN4c7r3X+ptB2bKmE4mIZKJGV0TkHOvXw8GDanSz5ffDoEEQHW19TZ2qU98i4lqauiAico5Fi6BQIbj9dtNJXOjMGejSBT75BF5/HV59FXw+06lERLKlM7riWdHR0aYjiIMCVe8FC+Cuu6BAgYDszjuOHYMHH7RuAPHFF9Cvn9EmV8d3eFG9xS41uuJZupNOeAlEvY8ds9bnagTsefbtg/vus053z5xpLVkwTMd3eFG9xS4tXRDPatOmjekI4qBA1HvhQjh50rrGSv5n61aIjIS//rJm5N56q+lEgI7vcKN6i106oysi8j+zZkHVqnDNNaaTuMTq1VCnjrVEYdUq1zS5IiI5pUZXRARrmMCsWTqbm+6776zRYVdfbc3Ivfpq04lERHJNja54VlxcnOkI4qC81nvtWtizR40uAKNGQVQUNGhgrecoU8Z0okx0fIcX1VvsUqMrnjVs2DDTEcRBea33rFlQooT1m/qw5ffDwIHwxBPQuTNMmQKFC5tOlSUd3+FF9Ra7dDGaeFZsbKzpCOKgvNZ75kxo3Bjy5QtQoFBz5gw884w1PuzNN6FvX1fPyNXxHV5Ub7FLja54VmGXnomS4MhLvffsgZ9+gu7dAxgolBw7Bq1bw7x51q19O3QwneiidHyHF9Vb7FKjKyJh77vvICLCmqIVdv7801qYvHEjzJ4NmlcqIh6iRldEwt6sWVC7tiuvuQquLVus7v7YMVi2DG65xXQiEZGA0sVo4lm9e/c2HUEcZLfeJ07A/PlhOG3hP/+xuvtLLrFm5IZYk6vjO7yo3mKXGl3xrIoVK5qOIA6yW+9Fi6wTmmHV6M6aZc3IrVIFVqyASpVMJ8o1Hd/hRfUWu9Toimd17drVdARxkN16T5wI1arBDTcEOJBbffoptGhhjZhYsABKlzadyBYd3+FF9Ra71OiKSNg6cQKmToVHH3X1JK3A8Puhf3946ilrjNjkyVCokOlUIiJBpYvRRCRsffcdHDliTdbytNOnrQb3yy/hrbegT58w6OxFRHRGVzwsISHBdARxkJ16x8Za12Bdd10QArnF0aPW7XzHjYOvvoIXX/REk6vjO7yo3mKXGl3xrD59+piOIA7Kbb2PHLGuyXr00SAFcoP//te66GzFCuv0dbt2phMFjI7v8KJ6i11auiCeFRMTYzqCOCi39Z45E44f9/Cyhd9/t2bknjhhzcitXt10ooDS8R1eVG+xS2d0xbM0jia85LbesbFQpw5cdVWQApn0ww/WmytY0JqR67EmF3R8hxvVW+xSoysiYefgQZg716PLFmbMgPr1rZlpcXEe7eRFRHJGja6IhJ2pUyE1FR5+2HSSAPv4Y3jgAWja1LrdW6lSphOJiBilRlc8a+jQoaYjiINyU+8JE6xrtK64Inh5HOX3wyuvWPNxu3Sx7oJRsKDpVEGl4zu8qN5ily5GE89KSUkxHUEclNN6//e/1m1/P/44yIGccvo0PPEEjB0Lw4fDCy94YnzYxej4Di+qt9jl8/v9/pxsGB8fjy8PH55VqlShSJEitn8+r9asWUPNmjWJj4+nRo0axnKIiFnvvQe9esHevSF799v/d+QIPPQQLFkCo0dD27amE4mI5Emg+7Ucn9GtVauW7Rfx+XzMnz+f+vXr296HiEhe+f3w6afWMtaQb3L37IFmzWDrVuvKOn2+iohkkqulC08++SR33HFHrl7g6NGjdOvWLVc/IyISDCtXwoYN1lndkJaQYM3IPXMGli+Hm24ynUhExJVy1ejWrVuXtrn81VhycrIaXTEiOTmZMmXKmI4hDslJvT/5BK65JsRPfq5YYd3St1w5mDMHKlQwncgIHd/hRfUWu3I8deHbb7/lnnvuyfULFC9enG+//ZabdMZBHNapUyfTEcRBF6v3gQPwzTfQuTNEhOq8malToUEDuPFG60xumDa5oOM73KjeYleOz+i2bNnS1gvky5fP9s+K5MWAAQNMRxAHXaze48ZZs3M7dnQkTuCNHAldu1rDf8eOhQIFTCcySsd3eFG9xa5QPa8hclGarhFeLlRvv99attCyJVx+uYOhAiEtDV56yZqP26OHNQQ4zJtc0PEdblRvsStXa3RHjBiR6xFjPXv2zNX2IiKBFrIXoZ06BY8/bp2OfucdeP5504lEREJKrhrd3r1752rnPp9Pja6IGPfpp3D11SF2Edpff1kzcpctg9hYaN3adCIRkZCTq6UL27Ztu+BXYmIi3377LdWrVwegWLFiQQktkhOjRo0yHUEclF29Dx4MwYvQdu+GunVh9WqYN09NbhZ0fIcX1VvsytXHfqVKlbL9OnToED169ODBBx8kKSmJ1157jaSkpCDFFrm4NWvWmI4gDsqu3l99ZY2bjY52OJBdGzZA7dqwfz/ExcG995pO5Eo6vsOL6i125fgWwNlZu3YtAwcOZMaMGZQoUYIePXrQo0cP153N1S2ARcKP329N4qpWDSZNMp0mB5Yvt2bkli9vzcgtX950IhERRwW6X7P9i7y1a9fSokULatasyfLlyxk4cCBJSUn079/fdU2uiISnZcvgt9/gySdNJ8mBKVOgYUOoXt1qeNXkiojkWa4b3TVr1qQ3uHFxcQwcOJDExET69eunBldEXGXoUOuMboMGppNcxPvvW/NxH3gA5s6FEiVMJxIR8YRcTV2Iiopi1qxZlCxZktdff53u3btTtGjRYGUTEbHt55+t3/5/9RXkciqic87OyB0+HHr1sjrzkLliTkTE/XLV6M6aNQuAokWLMnHiRL755ptst/X7/fh8Pn7++ee8JRSxKSoqihkzZpiOIQ45v97Dh0PFii4eWHDypHWFXGwsvPsudO9uOlFI0fEdXlRvsStXjW7dunXx+Xzk9Pq13N5cQiSQunTpYjqCOOjcem/fbt1AbMQIyJfPYKjsHD5sLVNYudKafdaqlelEIUfHd3hRvcWuXDW6S5YsCVIMkcBr1KiR6QjioHPr/c47ULw4PPGEwUDZ2bkTmjaFP/6A+fPh7rtNJwpJOr7Di+otduWq0RURcbv9++Hzz60lr5deajrNeX77DSIjrUXDK1bA9debTiQi4mm66kFEPGXkSGt+rut+07l0Kdx1F5QqBT/8oCZXRMQBOW50ixYtyiQbE9f3799PsWLFWLZsWa5/ViQvpk2bZjqCOGjatGmkpMAHH0CnTlC2rOlE5/jmG2jUCGrWtIb7Xnml6UQhT8d3eFG9xa4cN7rHjh3j9OnTtl7k6NGjnDlzxtbPitg1YcIE0xHEQRMmTODLL+HAAXjhBdNpzvHvf1ujHx5+GL77zlo8LHmm4zu8qN5iV66WLrRr146IiIhcfZV11WkVCScTJ040HUEc9PXXE3n7bXjkEahc2XQarBm5PXtaXy++CGPHQv78plN5ho7v8KJ6i105vhitf//+tl/E5/NR2RV/8oiIV02aBElJ8O23ppNgzcht394K9cEHLlwwLCISHnLc6A4YMCCIMURE7EtNhcGDrWWwt9xiOMzBg9aM3B9/hMmT4cEHDQcSEQlfGi8mIiFv/Hj49VdrrJhRf/wBTZrA7t2wYAHceafhQCIi4U3jxcSzoqOjTUcQB5w8Cf36QcWK0dx+u8Egv/wCtWvD0aPWHc/U5AaVju/wonqLXWp0xbN0J53w8Mkn1onUrl0N1nvxYmtGbtmysGoVVK1qLkuY0PEdXlRvsUuNrnhWmzZtTEeQIDtyBAYNgo4doVcvQ/WOjYXGjeH2262bQpQrZyZHmNHxHV5Ub7HLdY3usWPHeO2114iMjKRUqVJEREQwZsyYLLfduHEjkZGRFC1alNKlS9O+fXuSk5MdTiwiprzzDvz1Fxi5Vtbvh7ffhjZtrK9Zs6BYMQNBREQkO65rdPft28cbb7zBpk2bqF69OmCNJzvfzp07qVu3Ltu2bWPIkCH06tWL2bNn07BhQ9s3thCR0PHnn1af2aULVKjg8IunpkKPHtC7N7z8MowerRm5IiIulKdGd/v27Tz11FNUqVKFkiVLpt/mNzk5ma5du7J27dpc7/PKK69k7969JCYmMnz48Gy3Gzx4MMePH2fRokV06dKFvn378s0337B+/XpGjx5t9y2Jh8TFxZmOIEE0eDBEREDfvta/O1bvEyfg0UchJgY+/BDefBOy+Mu4BJeO7/CieotdthvdDRs2cMstt/DNN99QuXJlDh8+nH6b3zJlyhAXF0dMTEyu95s/f34uu+wyAPx+f7bbTZkyhebNm1O+fPn0x+677z6qVKnCN998k+vXFe8ZNmyY6QgSJElJ8NFH0KcPlC5tPeZIvQ8csIb1zppl3ZnimWeC/5qSJR3f4UX1Frtsz9Ht06cPJUqU4IcffiAiIiK9OT2rWbNmQWs4d+3axb59+7j11lszPVerVi3mzJkTlNeV0BIbG2s6ggTJa69ByZLW6oGzgl7vHTsgMtJaM7FokTVKTIzR8R1eVG+xy/YZ3WXLlvHMM89kanDPqlixIjt37rQd7EL27NkDQLksrm4uV64cBw4c0DpdoXDhwqYjSBD88gt89RX07w+XXvr/jwe13uvXwx13WMsWVq5Uk+sCOr7Di+otdtludNPS0rj03D9lzpOcnEyBAgXs7v6Cjh8/DpDl/gsWLJhhGxHxlldegauvhs6dHXrBhQvh7rutsWErV0KVKg69sIiI5JXtRveWW25h9uzZWT535swZYmNjueOOO2wHu5BChQoBcPLkyUzPnThxIsM2IuIdixfDzJnwxhuQL58DL/j119YtfevUgSVL4IorHHhREREJFNuN7ssvv8zcuXN5+umn+fXXXwHYu3cv8+fPp2HDhmzYsIGXXnopYEHPdXbJwtklDOfas2cPpUuXJl82fwo2bdqUqKioDF+1a9dm2rRpGbb7/vvviYqKyvTzzz33HKNGjcrw2Jo1a4iKiso0w/e1115j6NChGR7bsWMHUVFRJCQkZHj8gw8+oHfv3hkeS0lJISoqKtPVphMmTMjydoitW7fW+zjnfTz33HOeeB9eqUde30fPnr159lnrBmSPPpr5fZx9nwF5H88+y6hWraBdO3jsMZg5kzWbN6seLnof7du398T78Eo9gv0+zr9hRKi+D6/UI1DvY8KECem9WOXKlalevTrPP/98pv3kiT8Pxo4d6y9RooTf5/Nl+CpevLh//Pjxedm13+/3+1evXu33+Xz+MWPGZHrusssu8z/yyCOZHq9SpYq/QYMGmR6Pj4/3A/74+Pg855LQ8P7775uOIAE0ZIjf/7e/+f0//5z18wGr95kzfn+XLn4/+P39+vn9aWmB2a8ElI7v8KJ6h49A92u2py4A/Otf/+KBBx5g/vz5bN68mbS0NK699loaN25M0aJFA9OJZ+Ohhx5izJgx7Ny5M33E2MKFC9m8eTMvvPBCUF9bQkPXrl1NR5AA2bHDWq7QrRvceGPW2wSk3sePW2dwp0+HTz6BJ5/M+z4lKHR8hxfVW+zKU6MLUKRIER544IFAZEkXExPDoUOH2L17NwAzZsxgx44dAHTr1o1ixYrx8ssvM2nSJOrVq0f37t05cuQIw4cP56abbsry9LmIhK4ePaB48SDf6nf/foiKgrVrYdo0uP/+IL6YiIg4wXajO3/+fBYvXszgwYOzfP7ll1+mQYMG1K9fP9f7HjFiBNu3bwes2/9OnTqVb7/9Fp/PR/v27SlWrBjly5dn6dKl9OzZk5deeokCBQrQvHlzRowYke36XBEJPXPmwNSpEBsLxYoF6UWSkqwZufv3W1e83X57kF5IREScZPtitEGDBqWfZc3Krl27GDRokK19JyYmkpaWRlpaGqmpqaSmpqb/c8WKFdO3u/7665k7dy5Hjx5l//79jB07lrJly9p6TfGe8xfWS+g5cQK6doUGDeCRRy68re16r11rzcU9fdoaH6YmNyTo+A4vqrfYZbvR/eWXXy44PqxWrVqsX7/e7u5F8qxPnz6mI0geDR1qrc+NiQGf78Lb2qr3/PlQty6UL281uf/4h72g4jgd3+FF9Ra7bDe6J0+e5NSpUxd8PiUlxe7uRfIsJibGdATJg61bYcgQ6NULrrvu4tvnut5jx0LTplaju3gxXH65vaBihI7v8KJ6i122G90bbriBqVOnZvmc3+9n6tSpXH/99baDieTVuctcJLT4/daShcsvt+6ElhM5rrffD4MHQ4cO1tf06VCkiP2wYoSO7/Cieotdthvdbt26sWLFClq1asXPP//MmTNnOHPmDOvXr6dVq1asXLlS40BExJZp06yL0N5/Hy5wp/HcS02F556zuucBA+Czz+CSPA+fERERl7L9Cd+uXTu2bt3K66+/zrfffktEhNUzp6WlERERQb9+/ejYsWOgcopImDh2DLp3h2bNrGlfAZOSAm3bwqxZ8Pnn8PjjAdy5iIi4ke0zumDdSm7Tpk0MGzaMzp0707lzZ4YPH86mTZsYOHBgoDKK2HL+bQ4lNPTrB/v2WWdzL3YB2rkuWO/kZLjvPuvisxkz1OR6gI7v8KJ6i122zugeO3aMu+++myeffJKnn36aXr16BTqXSJ7pYsjQs3QpvPsuvP02XH117n4223pv22bNyD10CJYsgVq18pxTzNPxHV5Ub7HL1hndSy+9lKSkJHy5Od0i4jD9ViG0HD0K0dFw113W0oXcyrLe8fHWjFy/H1atUpPrITq+w4vqLXbZXroQGRnJvHnzAplFRMJY797w55/w5Zfwt78FYIdz58I990ClStaM3GuuCU+3C6gAACAASURBVMBORUQklNhudPv168fvv/9Ou3btiIuLY9euXRw4cCDTl4jIxcybBx9/bC1ZCEg/Ono0NG8O9erBokWgOyaKiIQl21MXbrjhBgA2bNjA+PHjs9zG5/ORmppq9yVE8iQ5OZkyZcqYjiEXceiQdW1Yw4bw1FP295OcnEyZ0qVh0CDo3x+efBJGjtT4MI/S8R1eVG+xy/afAP3797/oNlrDKyZ16tSJGTNmmI4hF9G9u7U+d9So3E1ZOF+n6GhmlCtnzcZ94w1rVq4+gzxLx3d4Ub3FLtuN7oABAwIYQyTw9P+o+02bZt2Jd8wYqFAhDzs6dowBhw9bd5n44gvrqjbxNB3f4UX1FrsC9ju948ePA1CoUKFA7VIkT2rUqGE6glzAvn3WUoUWLeBf/8rjjpo3p8Zvv1k3g4iMDFhGcS8d3+FF9Ra78nTDiO3bt9OxY0cuu+wyLr30Ui699FIuv/xyoqOj2b59e6AyiojH+P3wzDPWHXk/+SQPKwy2boU6dWD7dmsIr5pcERE5h+0zugkJCdx5550cPnyYhg0bUrVq1fTHx44dy8yZM1mxYgXXXXddwMKKiDfExsKUKTBpElx+uc2drF5t3Se4ZElrRm7lygHNKCIioc/2Gd2XXnqJv/3tb6xdu5Y5c+bw73//m3//+9/MmTOH9evXExERwYsvvhjIrCK5MmrUKNMRJAu7d8Nzz0GbNtCqlc2dfPcd3HuvNYtsxQqoXFn1DjOqd3hRvcUu243u0qVL6dq1KzfeeGOm5/75z3/StWtXlixZkpdsInmyZs0a0xHkPKmp0LEjFCgAMTE2dzJqFERFWfPIFi6E/40cUr3Di+odXlRvsct2o3v69GkKFy6c7fOFCxfm9OnTdncvkmcjR440HUHOM2QILFgA48ZBqVK5/GG/HwYMgCeesGbkTpkC53wGqd7hRfUOL6q32GW70b3lllv4/PPPOXToUKbnDh8+zKhRo3SVpIikW7IEXnvNupfDfffl8ofPnIHOnWHgQBg82LoRREDuEywiIl5m+2K0gQMH0rhxY6pVq0bHjh3TLzpLSEhgzJgx7N+/nxjbv5sUES/573+tNbn33AP9+uXyh48ehdat4fvvrYG77dsHJaOIiHiP7Ua3fv36zJkzh969ezN06NAMz1WvXp1x48ZRv379PAcUkdCWmgrt2lkrD8aPz+WJ2D//tCYrJCTA7NnQqFHQcoqIiPfkaY5ugwYNWLt2Lbt27WLlypWsXLmSXbt2sWbNGu7L9e8mRQIrKirKdATBWpe7cCF8/TVccUUufnDzZqhdG3buhGXLLtrkqt7hRfUOL6q32JWrM7p9+/alTZs23HTTTRkeL1euHOXKlQtoMJG86tKli+kIYc/2utwff4TmzaF0aWtGbqVKF/0R1Tu8qN7hRfUWu3J1Rnfo0KH8+uuv6f+enJxMREQEixYtCngwkbxqpF9zG2V7Xe6sWVCvHlSpYs3IzUGTC6p3uFG9w4vqLXblaemCiEhWzq7LTUvL5brcTz+FFi2sW/kuWGCd0RUREbFJja6IBFyu1+X6/db6hqeegmefte4NXKhQ0HOKiIi3qdEVz5o2bZrpCGHp7Lrcfv2gQYMc/MDp09CpE7zxBgwdCu+/b2tGruodXlTv8KJ6i125Hi+WmJiYfiu+szeL+P333ylRokSW2+umEWLKhAkTaNmypekYYeXsuty6da0TtBd19Ci0agWLFlm3S3vsMduvrXqHF9U7vKjeYpfP7/f7c7pxRETuTgD7fD5SU1NzHSoY1qxZQ82aNYmPj1fzLRIEp09bE8A2bIB16+Cig1j27rVm5G7eDFOn2rhdmoiIeE2g+7VcndH94osv8vyCIuJNzz8PcXHWydmLNrmbNlkXnJ08CcuXw803O5JRRETCS64a3Y4dOwYphoiEss8+g5Ej4eOP4e67L7LxqlVw//1w2WWwdClUrOhIRhERCT+6GE1E8iQuDp57Dp55xhqacEHTp0P9+nD99dYPqskVEZEgUqMrnhUdHW06guft2AEPPmjdqffddy+y8UcfWRs3bw7ffw+lSgU0i+odXlTv8KJ6i11qdMWzdCed4EpJgZYtoXBhmDwZ8ufPZkO/H15+2ZqP27UrTJwIBQsGPI/qHV5U7/CieotduR4vJhIq2rRpYzqCZ/n91ujbTZtg5UooWzabDU+dgs6dYexYePtt6NkTfL6gZFK9w4vqHV5Ub7FLja6I5Npbb1knZr/55gIDE44cgYcesu4gMX68NWBXRETEQWp0RSRXZs6EV16BV1+Fhx/OZqM9e6BpU9i2DebNg3r1HM0oIiICWqMrHhYXF2c6guds3GjdvCwqCgYOzGajhATr6rR9+6zJCg41uap3eFG9w4vqLXap0RXPGjZsmOkInnLwoNXgVqwIX30FWd4occUKuPNOKFLEmpd7442O5VO9w4vqHV5Ub7FLja54VmxsrOkInnHmDDz6KOzfb43CLVo0i42mToUGDazmNi4OKlRwNKPqHV5U7/CieotdanTFswoXLmw6gif4/dCjByxYYF18ds01WWwUE2NdeBYVZa3JLVHC8Zyqd3hRvcOL6i12qdEVkQt65x3r9r4ffWSdsM0gLQ1efNGaj9ujB0yYAAUKGMkpIiJyPk1dEJFsTZoEvXpB377w5JPnPXnqlDVMd/x4qxt+/nkjGUVERLKjM7riWb179zYdIaTFxcG//gVt28KgQec9efiwNT5s0iSIjXVFk6t6hxfVO7yo3mKXzuiKZ1WsWNF0hJC1aRO0aAF33AFffHHehIXdu6FJE9ixA77/Hu65x1jOc6ne4UX1Di+qt9jl8/v9ftMhnLBmzRpq1qxJfHw8NWrUMB1HxLX+/NNqcAsWtKaFlSx5zpMbNkBkpHWF2ty5cMMNxnKKiIj3BLpf09IFEUmXkgL33w/Hj8OcOec1ucuXWzNyS5SwZuSqyRUREZdToysiAKSmWutxf/sNZs2Cq64658nJk6FhQ7jlFqvhLV/eWE4REZGcUqMrnpWQkGA6Qsg4Oyt35kxrVm7Nmuc8+d578Mgj8OCD1mne4sWN5bwQ1Tu8qN7hRfUWu9Toimf16dPHdISQ8e9/W/d8+PBDa5gCYM3I7d3b6oB79YJx41w9I1f1Di+qd3hRvcUuTV0Qz4qJiTEdISRMmgQvvAAvvQRPPfW/B0+ehI4dYeJE64xut24mI+aI6h1eVO/wonqLXWp0xbM0jubiVqywZuW2aQNvvvm/Bw8dggcesC44mzTJurVvCPi/9u49Lsf7/wP467qLDkoUOecQZjOncopyyOggZS3CnGJjm/M2fM02h8z5sGHsa2LGlrOGnEIhpRBmxhxD5NCISAd1/f7o6/6tVdTl7r7u+7pez8fD4/Ht6r6v+9Ve32t7u/rcn5t9qwv7Vhf2TVJx0CVSqT/+yNthoU0bYPXq/+2Vm5SUt0furVtARATg5iZ3TCIiIsk46BKp0LVrQLdugIMD8Ntv/1t6+8cfeUOuRpP3sWhvvSV3TCIiotfCN6ORYs2ZM0fuCAbpzp28ncLKlQP27s3bFheHDgGuroCdXd6SBSMcctm3urBvdWHfJBUHXVKs9PR0uSMYnNRUwMMj7wMhIiKAKlWQ94azbt2Ali2Bw4eB6tXljikJ+1YX9q0u7Juk4kcAE6lEenrePHv+fN4827gxgIUL87Zc6N8fCAkBypaVOyYREakYPwKYiEosKwsICABOn877zIfGb+YC48b9/75iP//MIZeIiBSHb0YjUricHGDQIODAASA8HGjdNAPoOyhv67ClS4ERI+SOSEREVCp4R5cUKyUlRe4IshNFYNSovI/1DQ0F3nF+mLdId/t2YMsWRQ257Ftd2Le6sG+SioMuKdaQIUPkjiC7r78Gli8HVqwA/FvdzNsX948/8m7vvvuu3PF0in2rC/tWF/ZNUnHpAinW1KlT5Y4gq0WLgBkzgHnzgKGtfgfaeuWtwz16FGjUSO54Oqf2vtWGfasL+yapeEeXFEvNu2v89BPw6ad57zP73Olg3p3cKlXy9shV4JALqLtvNWLf6sK+SSoOukQKExYGDB0KDBsGzGwSCnh6Am3b5n0oRNWqcscjIiLSGw66RAqydy8QGAgEBIhYVm8+hPf7AX37Ajt3AtbWcscjIiLSKw66pFghISFyR9CrgweBnj2Bbl1FrK30KUz+Mx744ou8dQxlysgdr9SprW+1Y9/qwr5JKg66pFgJCQlyR9CbI0eAHj2Ajm452FymL8r+sDhvu4VvvgEEQe54eqGmvol9qw37Jqn4EcBERi42Nu+jfVu3yMbO556wOB0LrF8P+PrKHY2IiKhEdD2vcXsxIiN24kTee81avJWB7ffawyLlet4ahrZt5Y5GREQkOw66REbq9Om8O7mN6zxF+PXmKGeZA8TEAA0byh2NiIjIIHCNLpER+uMP4J13AMfKj7D76huwrmmTt4aBQy4REZEWB11SLF+FrlG9cAHo0gWoZfUAe680gI1rEyAqKu8DIVRMqX1T4di3urBvkoqDLinWyJEj5Y6gc5cuAe7uIuyFe4i43hC2A32A7dsBKyu5o8lOiX1T0di3urBvkoqDLilWt27d5I6gU9eu5Q25Nhl3sf9uE1T6egQQEqKKPXKLQ2l908uxb3Vh3yQV34xGZARu3ADcO+fC/OEdHEhvjSorZgAffih3LCIiIoPGQZfIwN26Bbh3yoFw5w4OohOqb/8B8PGROxYREZHB49IFUqywsDC5I7y2W7cAd7csZN28i4OWPVDr8C8ccoughL6p+Ni3urBvkoqDLilWaGio3BFey82bQMe2Gci4fg9R1d9HnbgNQOvWcscyWMbeN5UM+1YX9k1ScekCKdaGDRvkjiDZjRtA57bpyLlzH4eajEadiA2Avb3csQyaMfdNJce+1YV9k1QcdIkMTGIi0LlVGoSU+zjUeTpq7/gVKFdO7lhERERGh0sXiAzI1SsiOrZ4BJOUuzjU5wfU3reSQy4REZFEHHSJDMTlC8/RqflDmKXeRdSnO1Dr1zmAKX/pQkREJBUHXVKsoKAguSMU28Uzz9CpRSosn9xD1PyTqLlgHCAIcscyKsbUN70+9q0u7Juk4u0iUixj+SSdC7EP4d7xOWye/42Da2+hWv++ckcySsbSN+kG+1YX9k1ScdAlxerb1/AHxj/33oR7d3PY4W8c3J2BKh5d5I5ktIyhb9Id9q0u7Juk4tIFIpn8seEcOnuZw97kb0TGmKOKR3O5IxERESkKB10iGfy+LBqd+1RBNYuHOHimEuxb15E7EhERkeJw0CXFio6OljtCoU5P+w3uIxqhVvlHOHChJio1qiR3JEUw1L6pdLBvdWHfJBUHXVKsuXPnyh0hP1FE/LCVcJ/qhjqVnmL/pdqwq2UpdyrFMLi+qVSxb3Vh3yQVB11SrPXr18sd4f89f47DvvPxzo+90ajWU+y/6ABbe74XVJcMqm8qdexbXdg3ScVBlxTL0tJA7pY+fYp9rtPhuXMEWjZ6in1/1kKFitwjV9cMpm/SC/atLuybpOKgS1Sa7t1DWIup6BE3Ge6tnyA8oRqsrOQORUREpA4cdIlKy+XLCG0yEwGXZsG3Szq2HrGHhYXcoYiIiNSDgy4p1vjx4+V78fh4hLRYivfvLUT/954hdE9FlC0rXxw1kLVv0jv2rS7sm6TioEuK5eDgIM8Lh4djsetGfPDkW3wUlIlVG61hyvedlTrZ+iZZsG91Yd8klSCKoih3CH1ISEiAs7MzTp48CScnJ7njkFKtXImZwxIxWZyB8WOzMWdhGQh83xkREVGx6Hpe4x1dIl0QRYhTpuKLD+9hsjgD06bkcsglIiKSGX+hSvS6srMhDv8IY1c3xWKMwfx5Ij77nH+HJCIikhv/a0yKdeHChdJ/kSdPkOP7Lob91A5LhNH44Qfgs895G1cOeumbDAb7Vhf2TVJx0CXFmjBhQum+wN27yO7YBQP3D8AqYQjWrBEwfHjpviQVrdT7JoPCvtWFfZNUXLpAirV06dLSO/nFi3jm0ROByYuwB12xcaOA994rvZejVyvVvsngsG91Yd8kFQddUqxS247m2DE89u4D34yNOK5piR2/aeDhUTovRcXH7YfUhX2rC/smqTjoEpXEjh2433sEPE0icLVsA0Ts0qBdO7lDERERUWG4RpeouP77X9zwGwU301jctmqIQ4c55BIRERkyDrqkWHPmzNHNiUQR+OorXPhoEdqXO4WsytURfVRA06a6OT3phs76JqPAvtWFfZNURjvoRkVFQaPRFPonPj5e7nhkANLT01//JNnZQFAQTs7YBbdyCahQpwKiowU4Or7+qUm3dNI3GQ32rS7sm6Qy+jW6Y8aMQatWrfIdc+QUQgCmTZv2eidISwN69ULU/ufwNT+Gxk3KIDwcsLXVTT7Srdfum4wK+1YX9k1SGf2g6+bmBn9/f7ljkNLcuQN0747t5xugt+YXdHAzwdatgJWV3MGIiIiouIx26cILoigiLS0Nz58/lzsKKcVffwEuLvj5qiv8s0LRw88EO3ZwyCUiIjI2Rj/oBgUFwcbGBhYWFnB3d8fJkyfljkQGIiUlpeRPio0F2rXDd8+GYVDqdwgKErB+PWBmpvt8pFuS+iajxb7VhX2TVEY76JqZmSEgIACLFy/G9u3bMWPGDJw9exZubm44ffq03PHIAAwZMqRkTwgLg9jZHV+XW4Cxdydh4kRgxQrAxKR08pFulbhvMmrsW13YN0kliKIoyh1CV65cuYKmTZuiQ4cO2L17d77vJSQkwNnZGSdPnoSTk5NMCUmfEhISit/1smXIHTkao+vtxPdXPDF7NjBxYunmI90qUd9k9Ni3urBv9dD1vGa0d3QL4+joCD8/P0RGRkJB8ztJVKwLRBSBSZOQNWIs3m8Qj2VXPbBiBYdcY8T/CKoL+1YX9k1SKWrQBYCaNWsiKysLT58+LfT73t7e8PX1zffHxcUFYWFh+R63b98++Pr6Fnj+iBEjEBISku9YQkICfH19C6whmjJlSoFNrm/cuAFfX19cuHAh3/ElS5Zg/Pjx+Y6lp6fD19cX0dHR+Y6HhoYiKCioQLbAwED+HCX5ObKycCMgAN6z56FT7d3YmtgCmzYJ+PBDI/s5oJA++HPw5+DPwZ+DP4eqfo7Q0FDtLFa3bl00b94c48aNK3Ce16GopQsAEBAQgN27dxcYdLl0gfJ5/Bh47z3cO3Qe3WudwcUUO/z2G9Cpk9zBiIiI1ItLF/7n/v37BY6dOXMG27dvR7du3WRIRIbm33/D1bp9G+jYEdfi7qG9/UUkpdvh8GEOucauyL5Jkdi3urBvkspoB93AwED4+Pjgm2++wY8//ohx48ahXbt2sLKywuzZs+WORwYgISGh4MHz5wEXF5y5XQntzE4AFpY4ehRo1kz/+Ui3Cu2bFIt9qwv7JqmMdunCkiVL8Msvv+Dy5ct4/Pgx7O3t0aVLF0yZMgX16tUr8HguXSBERwO+vjhUwQ++KSFo0FCDXbsAe3u5gxERERGg+3nNaD8CeNSoURg1apTcMchYbN0K9OuHrY7j0e/KdLi6Cti2DbC2ljsYERERlRajXbpAVGxLlwIBAfjv20vQ68J09OwpIDycQy4REZHScdAl5crNBSZOhDhqFKa77MZHJz/EiBECfv2VH+lLRESkBhx0SZmysuDr4ICcuQswov1pTInxwDffAN99B2j4/3pFKmx/SVIu9q0u7JukMto1ukRFevQI8PfHsDsPENgmEdtia2LlSmDoULmDUWkaOXKk3BFIj9i3urBvkoqDLinLrVuAtzceX3+IBU1u49iZCti2DeDNAOXj/tnqwr7VhX2TVBx0STnOnQO8vHA7pwq8q/6F64kW2LcPcHOTOxgRERHJgasVSRkOHwZcXXHevAVchFikPLHAkSMccomIiNSMgy4Zv02bgK5dcbTeALS/vw3WFUwRGwtcvhwmdzLSo7Aw9q0m7Ftd2DdJxUGXjNu33wKBgQhrMwvv/PkdmjbTIDoaqFULCA0NlTsd6RH7Vhf2rS7sm6TioEvGKTcX+PxzYNw4LO+yGe8dHYcePQTs2QNUqJD3kA0bNsibkfSKfasL+1YX9k1S8c1oZHwyM4HBgyGu34Avu8Vj5r5WGDMGWLiQe+QSERHR/+OgS8YlNRV4911kxxzHBx2v4ed9tTFvHvDZZ4AgyB2OiIiIDAkHXTIeSUmAlxfSkh6hV/ObOBhTEb/8AvTrJ3cwIiIiMkT8RS8Zhz/+AFxccPdhWXSu8RdizlfEnj0vH3KDgoL0l49kx77VhX2rC/smqTjokuGLigJcXXHRqgVcTOJw+0HeHrnu7i9/Gj9JR13Yt7qwb3Vh3yQVB10ybOvXAx4eiGvQH+3vh8HMMm+P3GbNXv3Uvn37ln4+MhjsW13Yt7qwb5KKgy4ZroULgb59sbPdTHQ+twRvNNLg6FGgdm25gxEREZEx4KBLhic3F/j0U+Czz7DCcwv8Dn8KT08BERGAra3c4YiIiMhYcNAlw5KZCfTvD3HRt5jc7TiG7/HHJ58I2LQJsLAo2amio6NLJyMZJPatLuxbXdg3ScVBlwzH48eAtzeyNm/HALdrmLmvJebNAxYvBkxMSn66uXPn6j4jGSz2rS7sW13YN0nFfXTJMCQnA97eeHT1b/g3uYnouIpYvx4IDJR+yvXr1+suHxk89q0u7Ftd2DdJxUGX5PfXX4CnJ24+qwRv+4tIumqOiAigQ4fXO62lpaVu8pFRYN/qwr7VhX2TVFy6QPKKiwPat8cZoTlchFikZZsjJub1h1wiIiIiDrokn/BwoHNnRFTpD7f7W1Glet4euW++KXcwIiIiUgIOuiSPVasAPz+saTQL3hcXob2rgEOHgGrVdPcS48eP193JyOCxb3Vh3+rCvkkqDrqkX6IIzJgBcehQBLfYisGnxmDwYAHbtwNWVrp9KQcHB92ekAwa+1YX9q0u7JukEkRRFOUOoQ8JCQlwdnbGyZMn4eTkJHccdcrJAUaNQvbyH/GxUzxCElogOBiYPBkQBLnDERERkdx0Pa9x1wXSj2fPgP79kbZtP3q/fRX7f6+FNWuAgQPlDkZERERKxUGXSt/Dh4CfH5KPJ6F73Ru4fN0Gu3YBXbvKHYyIiIiUjGt0qXTdugV06IBzZ56jbYXzuJdpg+ho/Qy5Fy5cKP0XIYPBvtWFfasL+yapOOhS6fnrL6B9exy82xjtxSOoYG+GY8eApk318/ITJkzQzwuRQWDf6sK+1YV9k1QcdKl0HD8OuLri56w+8EwNReu2JjhyBKhZU38Rli5dqr8XI9mxb3Vh3+rCvkkqDrqkexEREDt1xnTzmRiUPBsDBggIDwfKl9dvDG5Hoy7sW13Yt7qwb5KKgy7p1oYNyPLuiSEVt2FK0ocIDgZWrgTKlJE7GBEREakNd10g3fn+ezwaORkBVeNw6F5jrF0L9O8vdygiIiJSK97RpdcnisCUKbg5cjZc7c7jREZj7NsnyD7kzpkzR94ApFfsW13Yt7qwb5KKd3Tp9eTkACNG4NR/49Dd+jzMrMshZpeAN9+UOxiQnp4udwTSI/atLuxbXdg3ScWPACbpMjOB99/H7q3P0LtsGBo1KYMdO4CqVeUORkRERMZI1/Maly6QNI8fA15eWPFbFfQQdqBztzKIiuKQS0RERIaDgy6V3N27yO3YGZNifDD8+ff4+BMNtm0DypWTOxgRERHR/+OgSyVz7Roy27vj/QtfYk7WOCxYACxeDJiYyB2soJSUFLkjkB6xb3Vh3+rCvkkqDrpUfPHx+Lu1F9659RPCRD9s2iTg008BQZA7WOGGDBkidwTSI/atLuxbXdg3ScVdF6h4tmzB5X5fw1uzD6lWNXFwuwYuLnKHermpU6fKHYH0iH2rC/tWF/ZNUvGOLr2cKAJz5yImYAHaIhYah5o4Fmf4Qy4A7q6hMuxbXdi3urBvkop3dKlo2dnAiBHY+GMqBpocQpu2pti2TYCtrdzBiIiIiF6Nd3SpcI8eQfTujjkhlRCIjQjoUwb79nHIJSIiIuPBQZcKSkxEtksHDD/cD//JnYmvvgLWrgXMzOQOVjIhISFyRyA9Yt/qwr7VhX2TVBx0Kb+4ODxu1QU9Epdgde4grFoFTJ9uuDsrvExCQoLcEUiP2Le6sG91Yd8kFT8CmP7fli24+f5/0N1kN66b1MPWbRp06SJ3KCIiIlILXc9rfDMa5e2sMG8eTk0MhY9FHEztKiJmt4DGjeUORkRERCQdB121y84GPvkEu1beQu8yx/Bm47LYsUNA1apyByMiIiJ6PVyjq2apqYC3N5atMkcPYSe6eJkhKopDLhERESkDB121SkxEbjtXfH7EFyNyl2DUaA22bgXKlZM7mO74+vrKHYH0iH2rC/tWF/ZNUnHpghqdPIl07wD0f/IDwrK64bvvgNGj5Q6leyNHjpQ7AukR+1YX9q0u7Juk4qCrNnv24O57n8BX2IE/xMYICxOg1L8od+vWTe4IpEfsW13Yt7qwb5KKg66arF6NPz9YiO5msciwrozD4QKcneUORURERFQ6uEZXDUQRmD4dB4esRTvTOFjVs0dcvIZDLhERESkaB12le/4cGD4cP025Cg9NBFp3tED0UQEODnIHK31hYWFyRyA9Yt/qwr7VhX2TVBx0lezpU4h+PfHVytoIwk8IGmqC8HABNjZyB9OP0NBQuSOQHrFvdWHf6sK+SSqu0VWqe/eQ6dUTQ34fg1/FQMyeDUyYAAiC3MH0Z8OGDXJHID1i3+rCvtWFfZNUHHSV6NIl/N21D3reXobjmlbYsAHo3VvuUERERET6xUFXaeLicNlrFLyfbMFD61o4uEODdu3kDkVERESkfxx0lWTHDkQHfIueuftgV9sax/aawNFR7lBERERE8uCb0ZTiv//Fer9QdHm+B43bWiM2nkNupqLGhgAAG41JREFUUFCQ3BFIj9i3urBvdWHfJBUHXWMnihAnf4mZH11HX/FX9O5rgn37TWBrK3cw+fGTdNSFfasL+1YX9k1ScdA1ZtnZyB70AT6YWReTMRNffyXi57UamJnJHcww9O3bV+4IpEfsW13Yt7qwb5KKa3SNVVoaUv0GISBqBA6bdMJPIcCgQSraO4yIiIjoFTjoGqM7d3D9naHofn4ebpVrgL3bTdC5s9yhiIiIiAwLly4Ym7/+wgmnYWh7fhXSqzkiJr4Mh9wiREdHyx2B9Ih9qwv7Vhf2TVJx0DUmMTH4rWUwOt5Zj9rNKuBYghnefFPuUIZr7ty5ckcgPWLf6sK+1YV9k1QcdI3Ftm34ruNWvPvkZ3j5mOJgtBns7eUOZdjWr18vdwTSI/atLuxbXdg3ScVB1wjkLFmG0f5JGPt8Pj4bm4uNYWVhaSl3KsNnyX9IqsK+1YV9qwv7Jqn4ZjRDJop48vlU9F3YErsFbyxfmouPPmFlRERERMXBqclQZWXhdt/P0GPrYFw0a4Id20zg5SV3KCIiIiLjwaULhujxY5ztOAJtt47HXds3ER1XlkOuBOPHj5c7AukR+1YX9q0u7Juk4qBraG7fxr4WE9H+2ALY1a+IuN8t0ayZ3KGMk4ODg9wRSI/Yt7qwb3Vh3ySVIIqiKHcIfUhISICzszNOnjwJJycnueMU7vx5/Oj6Ez5+8A083J5hfbg1rK3lDkVERESkH7qe13hH10DkHo7GJKc9GPZgDj4alIHfDnLIJSIiInodfDOaAcgI3YbB/bOxMXcMFsx4hnFfWEEQ5E5FREREZNx4R1dmKbN+xDv9KuM3oSc2r8/Bp5MtOOTqyIULF+SOQHrEvtWFfasL+yapOOjKJTcXlz+cg3ZfdMRFi+aIOmIK/8AycqdSlAkTJsgdgfSIfasL+1YX9k1ScemCHDIzEeMzE777R6GSvQbHYq1Qr57coZRn6dKlckcgPWLf6sK+1YV9k1S8o6tvjx5ho/NsuO+fhLffEhFz3pZDbinhdjTqwr7VhX2rC/smqTjo6pGYdAtzG61C4LkpCOj6CHsTKsPWVu5URERERMrEQVdPnp85h48bRWLinXH46qN7WLu3CszM5E5FREREpFwcdPUgbXc0erRMRsjTQKxa8BDTl9tzZwU9mDNnjtwRSI/Yt7qwb3Vh3yQV34xWypKW74DPCAdc07TAnt8y0cW3otyRVCM9PV3uCKRH7Ftd2Le6sG+Sih8BXIrOjF+H7vM7wcTSHOHRNni7BbcPIyIiIioKPwLYGOTmYo//CrjO90OVKsCxi7YccomIiIj0jIOurmVm4sc2K+GzbQg6v52CQ5droloN/mMmIiIi0jdOYDokPkzF5IYbMezEMHzkeR3bTteFlZXcqdQrJSVF7gikR+xbXdi3urBvkoqDro5kXr6J/nWjMfPGAMwfcQ1LdjnCxETuVOo2ZMgQuSOQHrFvdWHf6sK+SSruuqADD4+ew7vuj3Esqys2LkxCr3F15Y5EAKZOnSp3BNIj9q0u7Ftd2DdJxUH3NSWGxsK7f0XcExrhQFga2vvVlDsS/Y++dtcgw8C+1YV9qwv7Jqm4dOE1nPhmL9r2q4tMs/KIiS+D9n6V5I5ERERERP/DQVcKUcTOodvQ8UtX1Kn0BMcuVUJDJ77rjIiIiMiQcNAtqdxcLHPfDL9VvvBomIiDiY6oXKOs3KmoECEhIXJHID1i3+rCvtWFfZNUHHRLIDc9AxPe3IERUb0wqvM5bPqzMSzLCXLHoiIkJCTIHYH0iH2rC/tWF/ZNUvEjgIspI/khBjU9hU0pnbBo6DmMWdmkFFISERERqZeuPwKYuy4Uw99nkuDnchcnn7lgy6xLePc/HHKJiIiIDB0H3Ve4svM8vHuWxQOxDiLX30PbwDfkjkRERERExcA1ui8R//1xuPhWgmhqitjoXLQNrC13JCIiIiIqJg66RdjxWRQ6jWyM+jb3EXPBDvVdKssdiUrI19dX7gikR+xbXdi3urBvkoqD7r+JIv7rvxc9F7rB0+E8DlxvgEp1uEeuMRo5cqTcEUiP2Le6sG91Yd8kFQfdfxCf52Byq334aJsHPml1HJuuOMGifBm5Y5FE3bp1kzsC6RH7Vhf2rS7sm6Tim9H+J+vRMwxtEod1Nz0wzz8Wn212gcAtcomIiIiMltHe0c3MzMTEiRNRvXp1WFpaom3btti/f7+kcz269gDetc9h400XhH5+Ep9v4ZBLREREZOyMdtAdPHgwFi1ahAEDBmDx4sUwMTGBt7c3jh49WqLzJMXehFuj+zj5uD72LbuCPvOcSykx6VtYWJjcEUiP2Le6sG91Yd8klVEOuvHx8diwYQNmz56NOXPm4IMPPsDBgwdRu3ZtTJgwodjn+WPzBbi4apCaY43onY/Q8eO3SjE16ducOXPkjkB6xL7VhX2rC/smqYxy0N28eTNMTU0xbNgw7TEzMzMMHToUsbGxuHXr1ivPETnvBFx7VYWd2VMcO1kGjb25R67SVK7MLeHUhH2rC/tWF/ZNUhnloHvq1Ck0bNgQVlb5t/1q1aoVAOD06dMvff6vHx+Bx4SmaG13FYcvV0f1ZryAiIiIiJTGKAfd5ORkVKtWrcDxF8du375d5HN/Gnkc7//ghn71jyP8RhOUr849comIiIiUyCgH3WfPnsHMzKzAcXNzc+33i7IkthW+6nAIq/9qhzKW3COXiIiISKmMch9dCwsLZGZmFjiekZGh/f6/vRh+h7Tfhp6L3sWp06dKNyTJLj4+HgkJCXLHID1h3+rCvtWFfavH+fPnAbz8pmVJGOWgW61atUKXJyQnJwMAqlevXuB7iYmJAIBVR2dglfOMUs1HhsPZmdvFqQn7Vhf2rS7sW10SExPRvn371z6PUQ66LVq0QFRUFNLS0mBtba09HhcXBwBo3rx5ged4eHhg3bp1qFOnTqF3fImIiIhIXs+ePUNiYiI8PDx0cj5BFEVRJ2fSo/j4eLRt2xbz5s3DZ599BiDvk9LefvttVK5cGTExMTInJCIiIiK5GeUd3datW6NXr16YNGkS7t27B0dHR6xZswY3btzA6tWr5Y5HRERERAbAKO/oAnl3cL/66iusW7cODx8+RLNmzRAcHIyuXbvKHY2IiIiIDIDRDrpERERERC9jlPvoEhERERG9iuIH3czMTEycOBHVq1eHpaUl2rZti/3798sdi0pBVFQUNBpNoX/i4+Pljkev4enTp5gyZQo8PT1ha2sLjUaDNWvWFPrY8+fPw9PTE9bW1rCzs8PAgQORkpKi58T0Oorb9+DBgwu93t98800ZUpNUx48fx8iRI9G4cWNYWVmhdu3aCAwMxKVLlwo8lte38Stu37q6vo3yzWglMXjwYGzZsgXjxo1DgwYNsHr1anh7eyMyMlIn+7OR4RkzZgxatWqV75ijo6NMaUgX7t+/j+DgYNSuXRvNmzdHVFQUBEEo8LikpCR06NABFStWxKxZs5CWlob58+fj7NmziI+PR5ky/DREY1DcvgHAzMwMISEh+Y7Z2NjoIybpyJw5cxAbG4tevXqhadOmSE5OxtKlS+Hk5IRjx46hcePGAHh9K0Vx+wZ0dH2LChYXFycKgiAuWLBAeywjI0OsX7++2K5dOxmTUWmIjIwUBUEQt2zZIncU0rHMzEzx7t27oiiK4okTJ0RBEMQ1a9YUeNzHH38slitXTrx586b22P79+0VBEMQVK1boLS+9nuL2PWjQINHa2lrf8UjHYmJixOzs7HzHLl26JJqbm4v9+/fXHuP1rQzF7VtX17eily5s3rwZpqamGDZsmPaYmZkZhg4ditjYWNy6dUvGdFRaRFFEWloanj9/LncU0pGyZcvC3t4eQF6/RdmyZQt8fHxQs2ZN7bEuXbqgYcOG2LhxY6nnJN0obt8vvp+bm4vHjx/rIxqVAhcXF5ia5v8Fc/369fHWW2/hwoUL2mO8vpWhuH0Durm+FT3onjp1Cg0bNoSVlVW+4y9+rX369Gk5YlEpCwoKgo2NDSwsLODu7o6TJ0/KHYn04NatW7h//z5atmxZ4HutWrXCqVOnZEhFpS09PR3ly5dHhQoVYGdnh5EjR+Lp06dyx6LXJIoi7t69i0qVKgHg9a10/+77BV1c34peo5ucnIxq1aoVOP7i2O3bt/UdiUqRmZkZAgIC4O3tjUqVKuHcuXOYP38+3NzcEBMTU+hHQ5NyJCcnA0CR1/yDBw+QnZ3NdXwKUr16dUycOBFOTk7Izc3F7t27sWzZMpw5cwZRUVEwMTGROyJJ9Msvv+D27duYMWMGAF7fSvfvvgHdXd+KHnSfPXsGMzOzAsfNzc213yflcHFxgYuLi/ZrHx8fBAQEoGnTppg0aRJ2794tYzoqbS+u51dd8/wPoXLMnDkz39e9e/dGw4YNMXnyZGzevBmBgYEyJaPXceHCBYwYMQLt2rXDoEGDAPD6VrLC+gZ0d30reumChYUFMjMzCxzPyMjQfp+UzdHREX5+foiMjHzlWj8ybi+uZ17z6jZu3DhoNBocOHBA7igkwZ07d9C9e3dUrFgRmzdv1u62wetbmYrquyhSrm9F39GtVq1aocsTXvwKpHr16vqORDKoWbMmsrKy8PTp0wLrtUk5XvxK88X1/U/Jycmws7Pj3R4VMDc3h62tLR48eCB3FCqhR48ewcvLC48fP8aRI0dQtWpV7fd4fSvPy/ouipTrW9F3dFu0aIGLFy8iLS0t3/G4uDgA4JpNlbh69SosLCw45CpcjRo1ULlyZRw/frzA9+Lj43m9q0RaWhpSUlJQuXJluaNQCWRkZKBHjx64fPkydu7ciUaNGuX7Pq9vZXlV30WRcn0retANCAhATk4OVqxYoT2WmZmJ1atXo23btqhRo4aM6UjX7t+/X+DYmTNnsH37dnTr1k2GRKRv7733Hnbu3ImkpCTtsQMHDuDSpUvo1auXjMlI1zIzMwvcxACA4OBgAICnp6e+I5FEOTk5CAwMRFxcHDZt2oQ2bdoU+jhe38pQnL51eX0LosIXLgYGBmLbtm0YN24cHB0dsWbNGpw4cQIHDhyAq6ur3PFIh9zd3WFpaQkXFxfY29vjzz//xIoVK2BmZobY2Fi88cYbckek17B06VKkpqbi9u3b+OGHH+Dv76+9izN69GiUL18eSUlJaNGiBSpUqIAxY8YgLS0N8+bNg4ODA44fP85fbRqRV/X94MEDtGjRAv369dNe23v37sXu3bvh5eWF8PBwOeNTCYwdOxaLFy9Gjx49Ch1Y+/fvDwC8vhWiOH0nJibq7vp+7Y+cMHAZGRni+PHjxWrVqonm5uZimzZtxH379skdi0rB4sWLxTZt2oh2dnZimTJlxBo1aogDBw4Ur1y5Inc00oE6deqIgiCIgiCIGo1G1Gg02v99/fp17ePOnTsnenh4iOXKlRNtbW3FAQMGiPfu3ZMxOUnxqr5TU1PFAQMGiA0aNBDLlSsnmpubi02aNBFnz54tPn/+XO74VAKdOnXS9vvvPxqNJt9jeX0bv+L0rcvrW/F3dImIiIhInRS9RpeIiIiI1IuDLhEREREpEgddIiIiIlIkDrpEREREpEgcdImIiIhIkTjoEhEREZEicdAlIiIiIkXioEtEREREisRBl4iIiIgUiYMuERERESkSB10iUpypU6dCoym9f70NHjwYdevWLbXzq4lGo9H+Wbhwoc7Om5qamu/cCxYs0Nm5ich4cNAlIoP2008/5RtYLCwsUKNGDXh6emLJkiV48uRJgecIggBBEEo1V2mf31jt2rUL06ZNK9Fz/P39sW7dOnTv3l1nOaysrLBu3TosWrQIAPsiUisOukRkFIKDg7Fu3Tr88MMPGD16NABg7NixaNKkCc6ePZvvsV9++SWePXtWqnlEUSzV8xsrKYNu06ZN0a9fP7zxxhs6y2Fqaop+/frBz89PZ+ckIuNjKncAIqLi8PLygpOTk/briRMnIjIyEj4+PvD19cX58+dhbm4OADAxMYGJiYlcUVWPd0+JyFDwji4RGa3OnTvjq6++wvXr17Fu3Trt8cLW6EZERMDV1RUVK1aEtbU1GjVqhMmTJ2u/HxUVBY1Gg40bN+KLL75A1apVYWVlBT8/PyQlJb0yy/z589GuXTtUqlQJlpaWaNmyJbZs2VLoY9etW4fWrVujXLlysLW1RceOHREREZHvMbt374abmxusrKxQvnx5+Pj44M8//8z3mMGDB8Pa2hrXr19H9+7dYW1tjZo1a+L7778HAPz+++9wd3eHlZUV6tSpg9DQ0AJZUlNTMXbsWNSqVQvm5uZo0KAB5s6dm++OdWJionad64oVK+Do6Ahzc3O0bt0aJ06cyJdn2bJlEEUx33ITKV4sWYmOjsbIkSNRqVIlVKxYEcOHD0d2djYePnyIgQMHwtbWFnZ2dpg4caKk1yEiZeMdXSIyagMGDMAXX3yBiIgIfPDBB9rj/7yreO7cOfj4+KB58+YIDg6GmZkZLl26hJiYmALn++abb6DRaDBp0iTcvXsX3377Ld555x2cPn1ae8e4MIsXL4afnx8GDBiArKwshIaGolevXti5cye8vb21j5s2bRqmTZuG9u3bIzg4GGXLlsWxY8cQGRmJrl27AgDWrl2LwYMHw9PTE3PnzsXTp0+xfPlyuLq64tSpU6hdu7b2fDk5OfD29kanTp3Qo0cPrFu3DqNGjYKZmRm+/vprDBw4EAEBAVi+fDkGDhwIFxcX1KlTBwCQnp6Ojh07Ijk5GcOHD4eDgwOOHj2KSZMmITk5Wbu+9YVff/0VaWlp+PjjjwEAc+fOhb+/P65evQpTU1N89NFHSE5ORkRERL6/eLyOUaNGoXr16pgxYwZiY2Px448/wsbGBkeOHEH9+vUxa9YshIeHY968eXj77bcxYMAAnbwuESmESERkwFavXi0KgiCePHmyyMfY2NiIzs7O2q+nTJkiCoKg/XrRokWiIAji33//XeQ5IiMjRUEQxFq1aolPnjzRHt+0aZMoCIK4ePFi7bFBgwaJderUyff8jIyMfF9nZ2eLTZo0Ebt06aI9dunSJVGj0YjvvfdekTnS0tLEChUqiMOHD893/O7du2KFChXEYcOG5cshCII4e/Zs7bHU1FTR0tJSFARB3Lhxo/b4X3/9JQqCIE6dOlV7LDg4WLSyshIvX76c77UmTZokmpqaijdv3hRFURSvXbsmCoIgVq5cWUxNTdU+bvv27aIgCOLOnTu1x0aMGJHvn/2rCIIgTps2rcDxF717eXnlO96uXTtREATxk08+0R7LyckRa9WqJXbq1KnAeV5kX7BgQbEzEZFycOkCERk9KysrpKWlFfn9ihUrAgDCwsKQm5v70nMNHDgQ5cqV034dEBCAatWqYdeuXS99npmZmfZ/P3z4EKmpqXB1dUVCQoL2eFhYGERRxNdff13keSIiIvDo0SP06dMHKSkp2j8ajQatW7dGZGRkgef88062jY0NGjZsCGtra/Tq1Ut7vGHDhqhQoQKuXbumPbZp0yZ06NABFSpUyPdaXbp0QU5ODg4fPpzvdQIDA2FjY6P92tXVFQDynVPXhg4dmu/r1q1bFziu0Wjg7OyMq1evlloOIjJOXLpAREbvyZMnqFq1apHfDwwMxMqVK/HBBx/gP//5D7p06QJ/f38EBAQUeONUgwYNCjzf0dER169ff2mGnTt3YsaMGThz5gwyMzO1x/+5RvXKlSvQaDR46623ijzPpUuXAADu7u6Ffv+fgyYAWFhYwM7OrsBjatasWeC55cuXx8OHD/O91tmzZ1G5cuUCjxUEAffv3893zMHBId/XL/4C8c9z6tq/X/PFz1+rVq0Cx0szBxEZJw66RGTUkpKS8PjxY9SvX7/Ix5ibm+Pw4cOIjIxEeHg49uzZgw0bNsDd3R379u177Q+XOHLkCHx9fdGpUycsX74c1apVQ5kyZbBq1Sr8+uuvJTrXizvO69atK3R4NzXN/6/torIXteuE+I83mYmiiG7dumHChAmFPvbfQ39xzqlrRb1mYcdLMwcRGScOukRk1NauXQsA8PDweOnjBEGAu7s73N3dsWDBAsyaNQuTJ09GZGQkunTpon3cxYsX8z1PFEVcvnwZzZs3L/LcW7ZsgaWlJfbu3YsyZcpoj4eEhOR7XP369ZGbm4tz586hWbNmhZ7rxcBeuXLlIu/q6oqjoyPS0tJ0+jrcWoyIDAnX6BKR0Tp48CCCg4NRr149vP/++0U+rrBfab8YNLOysvId//nnn/N92trmzZtx584deHl5FXn+F3cXnz9/rj2WmJiIsLCwfINfz549odFoMH369CLvPnp4eKB8+fKYOXNmvvO9kJKSku/r1xkse/fujdjYWOzbt6/A91JTU5GTk1Pic75Y3/zo0SPJuaTikE1E/8Y7ukRkFHbt2oU///wTz58/x927d3Hw4EHs378fderUwfbt21G2bNkinztt2jQcOXIE3bt3h4ODA+7du4dly5ahVq1a2jdUvWBnZwdXV1cEBQXhzp07+Pbbb9GgQQN8+OGHRZ7fx8cHixYtgqenJ/r27as9f4MGDfD7779rH+fo6IjJkycjODgYbm5uePfdd2FmZobjx4+jRo0amDlzJqytrbF8+XIMGDAATk5O6NOnDypVqoQbN24gPDwcrq6uWLJkifacRQ3Mxfk1/vjx47F9+3b4+Phg8ODBcHJywtOnT3H27Fls2bIF169fh62t7SvP808tW7YEAIwePRrdunWDiYkJ+vTpU6JzSMWlC0T0bxx0icigvbhL92KngrJly8LW1hZNmzbFd999h6CgoHy7JLx4zj/v7vn5+eH69etYtWoVUlJSULlyZXTq1AnTpk2DtbV1vud+8cUXOHPmDGbNmoUnT56ga9euWLZsWb49dP99/s6dOyMkJASzZ8/GuHHjUK9ePcydOxfXrl0r8PHE06ZNQ926dbFkyRJ8+eWXsLS0RLNmzTBo0CDtY/r27Yvq1atj9uzZmDdvHjIzM1GzZk24ublhyJAhReYozvF/srCwwKFDhzBz5kxs2rQJP//8M8qXL4833ngD06dPR/ny5Quc41X8/f0xatQorF+/XruXrtRBt6Q/G+/oEtG/CSL/CkxEhKioKLi7u2Pz5s3w9/eXO45qaDQajB8/HuPHj4eVldVLP5SjpFJSUnDz5k04Oztj/vz5+PTTT3V2biIyDlyjS0REspo3bx7s7e2xbNkynZ0zNTUV9vb2cHZ25p1eIhXj0gUiIpJNRESEdhAtbA9jqaytrbF//37t17o8NxEZDw66RET/wzt/+vfPrd10ycTEpNS3ZyMiw8c1ukRERESkSFyjS0RERESKxEGXiIiIiBSJgy4RERERKRIHXSIiIiJSJA66RERERKRIHHSJiIiISJE46BIRERGRInHQJSIiIiJF+j9uLD6kEO5elgAAAABJRU5ErkJggg=="
|
||
],
|
||
"text/plain": [
|
||
"Figure(PyObject <matplotlib.figure.Figure object at 0x114260a90>)"
|
||
]
|
||
},
|
||
"metadata": {},
|
||
"output_type": "display_data"
|
||
},
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"PyObject <matplotlib.text.Text object at 0x1158e5c90>"
|
||
]
|
||
},
|
||
"execution_count": 14,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"PyPlot.plot(-reshape(uall[2, 3, 1, :], N), -t*F[2, 3], color=\"red\", linewidth=1.0, label=\"linear\")\n",
|
||
"PyPlot.plot(-reshape(uall[2, 3, 2, :], N), -t*F[2, 3], color=\"blue\", linewidth=1.0, label=\"nonlinear\")\n",
|
||
"PyPlot.legend(loc=\"best\")\n",
|
||
"PyPlot.grid()\n",
|
||
"PyPlot.xlabel(\"Displacement [m]\")\n",
|
||
"PyPlot.ylabel(\"Force [N]\")"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Last remarks. For some unknown reason this notebooks takes occasionally a very long time to complete for 40 different solutions. Sometimes NLsolve takes a very long time when using Newton algorithm. Without Newton all passes well. Actually my implementation of Newton method does not have any problems with speed."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 15,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"elapsed time: 57.255576466 seconds\n"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"57.255576466"
|
||
]
|
||
},
|
||
"execution_count": 15,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"toc()"
|
||
]
|
||
}
|
||
],
|
||
"metadata": {
|
||
"kernelspec": {
|
||
"display_name": "Julia 0.3.8",
|
||
"language": "julia",
|
||
"name": "julia-0.3"
|
||
},
|
||
"language_info": {
|
||
"name": "julia",
|
||
"version": "0.3.8"
|
||
}
|
||
},
|
||
"nbformat": 4,
|
||
"nbformat_minor": 0
|
||
}
|