diff --git a/README.md b/README.md index a29c9e9..bfcac3a 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,4 @@ +![JuliaFEMLogo](https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/geometry/JuliaFEMLogo_256x256.png) # JuliaFEM [![Join the chat at https://gitter.im/JuliaFEM/JuliaFEM](https://badges.gitter.im/Join%20Chat.svg)](https://gitter.im/JuliaFEM/JuliaFEM?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge) @@ -6,7 +7,7 @@ Build Status: [![Build Status](https://travis-ci.org/JuliaFEM/JuliaFEM.jl.svg?br Code Coverage: [![Coverage Status](https://coveralls.io/repos/JuliaFEM/JuliaFEM.jl/badge.svg?branch=master)](https://coveralls.io/r/JuliaFEM/JuliaFEM.jl?branch=master) -http://juliafem.readthedocs.org/en/latest/ +Documentation: [![Documentation Status](https://readthedocs.org/projects/juliafem/badge/?version=latest)](https://readthedocs.org/projects/juliafem/?badge=latest) The JuliaFEM project develops open-source software for reliable, scalable, distributed Finite Element Method. diff --git a/REQUIRE b/REQUIRE index c09da18..15e0e0d 100644 --- a/REQUIRE +++ b/REQUIRE @@ -2,5 +2,6 @@ julia 0.4- Logging LightXML Docile +Lexicon FactCheck HDF5 diff --git a/docs/2015-05-23-one-element-solution.ipynb b/docs/2015-05-23-one-element-solution.ipynb deleted file mode 100644 index b1599e6..0000000 --- a/docs/2015-05-23-one-element-solution.ipynb +++ /dev/null @@ -1,920 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#One element solution with and without nonlinear geometrics.\n", - "\n", - "Author: Jukka Aho\n", - "\n", - "Email: " - ] - }, - { - "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 )" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/plain": [ - "PyObject " - ] - }, - "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 -} diff --git a/docs/README.md b/docs/README.md index 0f7f226..a29c9e9 100644 --- a/docs/README.md +++ b/docs/README.md @@ -1,3 +1,41 @@ -[FEATURES / TODO](https://github.com/ovainola/JuliaFEM/blob/master/docs/features.md) +# JuliaFEM -Use `git mv` to move files in repository to make file tracking work properly. +[![Join the chat at https://gitter.im/JuliaFEM/JuliaFEM](https://badges.gitter.im/Join%20Chat.svg)](https://gitter.im/JuliaFEM/JuliaFEM?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge) + +Build Status: [![Build Status](https://travis-ci.org/JuliaFEM/JuliaFEM.jl.svg?branch=master)](https://travis-ci.org/JuliaFEM/JuliaFEM.jl) + +Code Coverage: [![Coverage Status](https://coveralls.io/repos/JuliaFEM/JuliaFEM.jl/badge.svg?branch=master)](https://coveralls.io/r/JuliaFEM/JuliaFEM.jl?branch=master) + +http://juliafem.readthedocs.org/en/latest/ + +The JuliaFEM project develops open-source software for reliable, scalable, distributed Finite Element Method. + +The JuliaFEM software library is a framework that allows for the distributed processing of large Finite Element Models across clusters of computers using simple programming models. It is designed to scale up from single servers to thousands of machines, each offering local computation and storage. The basic design principle is: everything is nonlinear. All physics models are nonlinear from which the linearization are made as a special cases. + +JuliaFEM current status: project planning + +Initial road map for JuliaFEM: + +version | number of degree of freedom | number of cores +----------|-----------------------------------------|---------------------- +1.0 | 100 000 000 | 1 000 +2.0 | 1 000 000 000 | 10 000 +3.0 | 10 000 000 000 | 100 000 + +We strongly believe in the test driven development as well as building on top of previous work. Thus all the new code in this project should be 100% tested. Also other people have wisdom in style as well: + +[The Zen of Python](https://www.python.org/dev/peps/pep-0020/) +``` + Beautiful is better than ugly. + Explicit is better than implicit. + Simple is better than complex. + Complex is better than complicated. + Flat is better than nested. + Sparse is better than dense. + Readability counts. + Errors should never pass silently. +``` + +Interested in participating? Please start by reading [CONTRIBUTING.md](https://github.com/JuliaFEM/JuliaFEM/blob/master/CONTRIBUTING.md). + +Contributors: see [contributors](https://github.com/JuliaFEM/JuliaFEM/blob/master/contributors) diff --git a/docs/api/JuliaFEM.elasticity_solver.md b/docs/api/JuliaFEM.elasticity_solver.md new file mode 100644 index 0000000..d466b29 --- /dev/null +++ b/docs/api/JuliaFEM.elasticity_solver.md @@ -0,0 +1,154 @@ +# JuliaFEM.elasticity_solver + +## Internal + +--- + + +#### assemble!(fe, eldofs_, I, V) [¶](#method__assemble.1) +Assemble global RHS to I,V ready for sparse format + +Parameters +---------- +fe : local vector +eldofs_ : Array + degrees of freedom +I,V : Arrays for sparse matrix + +Notes +----- +eldofs can also be node ids for convenience. In that case dimension +is calculated and eldofs are "extended" to problem dimension. + + +*source:* +[JuliaFEM/src/elasticity_solver.jl:171](https://github.com/JuliaFEM/JuliaFEM.jl/tree/92c5e6c15a1ffaea4c153cba2af3a62ba3b42ebe/src/elasticity_solver.jl#L171) + +--- + + +#### assemble!(ke, eldofs_, I, J, V) [¶](#method__assemble.2) +Assemble global stiffness matrix to I,J,V ready for sparse format + +Parameters +---------- +ke : local matrix +eldofs_ : Array + degrees of freedom +I,J,V : Arrays for sparse matrix + +Notes +----- +eldofs can also be node ids for convenience. In that case dimension +is calculated and eldofs are "extended" to problem dimension. + + +*source:* +[JuliaFEM/src/elasticity_solver.jl:130](https://github.com/JuliaFEM/JuliaFEM.jl/tree/92c5e6c15a1ffaea4c153cba2af3a62ba3b42ebe/src/elasticity_solver.jl#L130) + +--- + + +#### calc_local_matrices!(X, u, R, Kt, N, dNdchi, lambda_, mu_, ipoints, iweights) [¶](#method__calc_local_matrices.1) +Calculate local tangent stiffness matrix and residual force vector R = T - F + + +*source:* +[JuliaFEM/src/elasticity_solver.jl:68](https://github.com/JuliaFEM/JuliaFEM.jl/tree/92c5e6c15a1ffaea4c153cba2af3a62ba3b42ebe/src/elasticity_solver.jl#L68) + +--- + + +#### eliminate_boundary_conditions(dirichletbc, I, J, V) [¶](#method__eliminate_boundary_conditions.1) +Eliminate Dirichlet boundary conditions from matrix + +Parameters +---------- +dirichletbc : array [dim x nnodes] +I, J, V : sparse matrix arrays + +Returns +------- +I, J, V : boundary conditions removed + +Notes +----- +pros: +- matrix assembly remains positive definite +cons: +- maybe inefficient because of extra sparse matrix operations. (It's hard to remove stuff from sparse matrix.) +- if u != 0 in dirichlet boundary requires extra care + +Raises +------ +Exception, if displacement boundary conditions given, i.e. +DX=2 for some node, for example. + + + +*source:* +[JuliaFEM/src/elasticity_solver.jl:218](https://github.com/JuliaFEM/JuliaFEM.jl/tree/92c5e6c15a1ffaea4c153cba2af3a62ba3b42ebe/src/elasticity_solver.jl#L218) + +--- + + +#### eliminate_boundary_conditions(dirichletbc, I, V) [¶](#method__eliminate_boundary_conditions.2) +Eliminate Dirichlet boundary conditions from vector + +Parameters +---------- +dirichletbc : array [dim x nnodes] +I, V : sparse vector arrays + +Returns +------- +I, V : boundary conditions removed + +Notes +----- +pros: +- matrix assembly remains positive definite +cons: +- maybe inefficient because of extra sparse matrix operations. (It's hard to remove stuff from sparse matrix.) +- if u != 0 in dirichlet boundary requires extra care + +Raises +------ +Exception, if displacement boundary conditions given, i.e. +DX=2 for some node, for example. + + +*source:* +[JuliaFEM/src/elasticity_solver.jl:257](https://github.com/JuliaFEM/JuliaFEM.jl/tree/92c5e6c15a1ffaea4c153cba2af3a62ba3b42ebe/src/elasticity_solver.jl#L257) + +--- + + +#### interpolate{T<:Real}(field::Array{T<:Real, 1}, basis::Function, ip) [¶](#method__interpolate.1) +Interpolate field variable using basis functions f for point ip. +This function tries to be as general as possible and allows interpolating +lot of different fields. + +Parameters +---------- +field :: Array{Number, dim} + Field variable +basis :: Function + Basis functions +ip :: Array{Number, 1} + Point to interpolate + + +*source:* +[JuliaFEM/src/elasticity_solver.jl:30](https://github.com/JuliaFEM/JuliaFEM.jl/tree/92c5e6c15a1ffaea4c153cba2af3a62ba3b42ebe/src/elasticity_solver.jl#L30) + +--- + + +#### solve_elasticity_increment!(X, u, du, elmap, nodalloads, dirichletbc, lambda, mu, N, dNdchi, ipoints, iweights) [¶](#method__solve_elasticity_increment.1) +Solve one increment of elasticity problem + + +*source:* +[JuliaFEM/src/elasticity_solver.jl:278](https://github.com/JuliaFEM/JuliaFEM.jl/tree/92c5e6c15a1ffaea4c153cba2af3a62ba3b42ebe/src/elasticity_solver.jl#L278) + diff --git a/docs/api/JuliaFEM.md b/docs/api/JuliaFEM.md new file mode 100644 index 0000000..0b7a06a --- /dev/null +++ b/docs/api/JuliaFEM.md @@ -0,0 +1,104 @@ +# JuliaFEM + +## Exported + +--- + + +#### add_elements(model, elements) [¶](#method__add_elements.1) +Add new elements to model. +Parameters +---------- +list of dicts, dict = {element_type => elcode, elids => [node ids..]} + +Examples +-------- +Create two tet4 element and add them: + +>>> m = new_model() +>>> const TET4 = 0x6 +>>> el1 = Dict("element_type" => TET4, "node_ids" => [1, 2, 3, 4]) +>>> el2 = Dict("element_type" => TET4, "node_ids" => [4, 3, 2, 1]) + +In dict key means element id + +>>> elements = Dict(1 => el1, 2 => el2) +>>> add_elements(m, elements) + + + +*source:* +[JuliaFEM/src/JuliaFEM.jl:102](https://github.com/JuliaFEM/JuliaFEM.jl/tree/92c5e6c15a1ffaea4c153cba2af3a62ba3b42ebe/src/JuliaFEM.jl#L102) + +--- + + +#### get_elements(model, element_ids) [¶](#method__get_elements.1) +Get subset of elements from model. +Parameters +---------- +element_ids : list of ints + Element id numbers + +Returns +------- +Dict +{element_type = XXX, node_ids = [a, b, c, d, e, ..., n]} + + +*source:* +[JuliaFEM/src/JuliaFEM.jl:126](https://github.com/JuliaFEM/JuliaFEM.jl/tree/92c5e6c15a1ffaea4c153cba2af3a62ba3b42ebe/src/JuliaFEM.jl#L126) + +--- + + +#### get_field(field_type, field_name) [¶](#method__get_field.1) +Get field from model. + +Parameters +---------- +field_type : Dict() + Target topology (model.model, model.nodes, model.elements, + model.element_nodes, model.element_gauss +field_name : str + Field name + +create_if_doesnt_exist : bool, optional + If field doesn't exists, create one and return empty field + +Raises +------ +Error, if field not found and create_if_doesnt_exist == false + + +*source:* +[JuliaFEM/src/JuliaFEM.jl:67](https://github.com/JuliaFEM/JuliaFEM.jl/tree/92c5e6c15a1ffaea4c153cba2af3a62ba3b42ebe/src/JuliaFEM.jl#L67) + +--- + + +#### new_model() [¶](#method__new_model.1) +Initialize empty model. + +Parameters +---------- +None + +Returns +------- +New model struct + + +*source:* +[JuliaFEM/src/JuliaFEM.jl:41](https://github.com/JuliaFEM/JuliaFEM.jl/tree/92c5e6c15a1ffaea4c153cba2af3a62ba3b42ebe/src/JuliaFEM.jl#L41) + +--- + + +#### JuliaFEM.Model [¶](#type__model.1) +Basic model + + +*source:* +[JuliaFEM/src/JuliaFEM.jl:19](https://github.com/JuliaFEM/JuliaFEM.jl/tree/92c5e6c15a1ffaea4c153cba2af3a62ba3b42ebe/src/JuliaFEM.jl#L19) + diff --git a/docs/api/README.md b/docs/api/README.md new file mode 100644 index 0000000..f11064e --- /dev/null +++ b/docs/api/README.md @@ -0,0 +1,2 @@ +Files in this directory are generated using the `build.jl` script. Make +all changes to the originating docstrings/files rather than these ones. diff --git a/docs/api/README_new.md b/docs/api/README_new.md new file mode 100644 index 0000000..f11064e --- /dev/null +++ b/docs/api/README_new.md @@ -0,0 +1,2 @@ +Files in this directory are generated using the `build.jl` script. Make +all changes to the originating docstrings/files rather than these ones. diff --git a/docs/api/empty.txt b/docs/api/empty.txt deleted file mode 100644 index 8b13789..0000000 --- a/docs/api/empty.txt +++ /dev/null @@ -1 +0,0 @@ - diff --git a/docs/api/index.md b/docs/api/index.md new file mode 100644 index 0000000..bcc500b --- /dev/null +++ b/docs/api/index.md @@ -0,0 +1,43 @@ +# API-INDEX + + +## MODULE: JuliaFEM + +--- + +## Methods [Exported] + +[add_elements(model, elements)](JuliaFEM.md#method__add_elements.1) Add new elements to model. + +[get_elements(model, element_ids)](JuliaFEM.md#method__get_elements.1) Get subset of elements from model. + +[get_field(field_type, field_name)](JuliaFEM.md#method__get_field.1) Get field from model. + +[new_model()](JuliaFEM.md#method__new_model.1) Initialize empty model. + +--- + +## Types [Exported] + +[JuliaFEM.Model](JuliaFEM.md#type__model.1) Basic model + +## MODULE: JuliaFEM.elasticity_solver + +--- + +## Methods [Internal] + +[assemble!(fe, eldofs_, I, V)](JuliaFEM.elasticity_solver.md#method__assemble.1) Assemble global RHS to I,V ready for sparse format + +[assemble!(ke, eldofs_, I, J, V)](JuliaFEM.elasticity_solver.md#method__assemble.2) Assemble global stiffness matrix to I,J,V ready for sparse format + +[calc_local_matrices!(X, u, R, Kt, N, dNdchi, lambda_, mu_, ipoints, iweights)](JuliaFEM.elasticity_solver.md#method__calc_local_matrices.1) Calculate local tangent stiffness matrix and residual force vector R = T - F + +[eliminate_boundary_conditions(dirichletbc, I, J, V)](JuliaFEM.elasticity_solver.md#method__eliminate_boundary_conditions.1) Eliminate Dirichlet boundary conditions from matrix + +[eliminate_boundary_conditions(dirichletbc, I, V)](JuliaFEM.elasticity_solver.md#method__eliminate_boundary_conditions.2) Eliminate Dirichlet boundary conditions from vector + +[interpolate{T<:Real}(field::Array{T<:Real, 1}, basis::Function, ip)](JuliaFEM.elasticity_solver.md#method__interpolate.1) Interpolate field variable using basis functions f for point ip. + +[solve_elasticity_increment!(X, u, du, elmap, nodalloads, dirichletbc, lambda, mu, N, dNdchi, ipoints, iweights)](JuliaFEM.elasticity_solver.md#method__solve_elasticity_increment.1) Solve one increment of elasticity problem + diff --git a/docs/build.jl b/docs/build.jl index 56cc80c..94e8c4d 100644 --- a/docs/build.jl +++ b/docs/build.jl @@ -3,6 +3,12 @@ using Docile, Lexicon, JuliaFEM const api_directory = "api" const modules = [JuliaFEM, JuliaFEM.elasticity_solver] +main_folder = dirname(dirname(@__FILE__)) +this_folder = dirname(@__FILE__) + +file_ = "README.md" +run(`cp $main_folder/$file_ $this_folder`) + cd(dirname(@__FILE__)) do # Generate and save the contents of docstrings as markdown files. index = Index() @@ -12,7 +18,7 @@ cd(dirname(@__FILE__)) do save(joinpath(api_directory, "index.md"), index; md_subheader = :category) # Add a reminder not to edit the generated files. - open(joinpath(api_directory, "README_new.md"), "w") do f + open(joinpath(api_directory, "README.md"), "w") do f print(f, """ Files in this directory are generated using the `build.jl` script. Make all changes to the originating docstrings/files rather than these ones. diff --git a/docs/doc.conf.py b/docs/doc.conf.py deleted file mode 100644 index 9eb45cc..0000000 --- a/docs/doc.conf.py +++ /dev/null @@ -1,15 +0,0 @@ -import os -import sys -import re - -import juliadoc - -extensions = ['sphinx.ext.mathjax', - 'juliadoc.julia', - 'juliadoc.jldoctest', - 'juliadoc.jlhelp'] - -master_doc = 'index' - -html_theme_path = [juliadoc.get_theme_dir()] -html_sidebars = juliadoc.default_sidebars() diff --git a/docs/features.md b/docs/features.md deleted file mode 100644 index bfb80d1..0000000 --- a/docs/features.md +++ /dev/null @@ -1,29 +0,0 @@ -=============== -FEATURES / TODO -=============== - -- Parallel design -- Sparse matrices -- Discontinuous Galerkin? -- 100% Tested - - Automatic coverage (Travis CI) -- Doctest examples and tutorials - - Intuitive to use -- Marketing - - Engineering porn -- Field functions -- Multiphysics platform - - Elasticity - - Thermal implemented -- Transient and steady solvers - - Implicit dynamics -- Nonlinearities - - Geometrical - - Material - - Contact - - Mortar -- Modular design - - E.g. contact formulation can be altered by user -- Mesh and results format, Xdmf? -- No scalars but field variables, e.g. no constant 210GPa for steel, we interpolate variable from nodes. - - (Constant is a special case of field variable). diff --git a/docs/requirements.txt b/docs/requirements.txt deleted file mode 100644 index 0cfb48a..0000000 --- a/docs/requirements.txt +++ /dev/null @@ -1 +0,0 @@ --e git+https://github.com/JuliaLang/JuliaDoc.git#egg=JuliaDoc diff --git a/mkdocs.yml b/mkdocs.yml new file mode 100644 index 0000000..3ad4b15 --- /dev/null +++ b/mkdocs.yml @@ -0,0 +1,14 @@ +# This is automatically generated by docs/build.jl. Edit that file if you +# wish to make any permenant changes. +site_name: JuliaFEM.jl +site_description: JuliaFEM.jl, open-source software for reliable, scalable, distributed Finite Element Method. +repo_name: GitHub +docs_dir: 'docs' +site_dir: 'site' +repo_url: https://github.com/JuliaFEM/JuliaFEM.jl +pages: +- Home: 'README.md' +- Overview: 'api/index.md' +- API Docs: + - JuliaFEM: 'api/JuliaFEM.md' + - JuliaFEM.elasticity_solver: 'api/JuliaFEM.elasticity_solver.md' diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 96f249a..9d3123a 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -2,7 +2,9 @@ # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md module JuliaFEM + VERSION < v"0.4-" && using Docile +using Lexicon include("elasticity_solver.jl") include("xdmf.jl") @@ -11,9 +13,9 @@ include("abaqus_reader.jl") export Model, new_model, new_field, get_field, add_nodes, get_nodes, add_elements, get_elements -@doc """ +""" Basic model -""" -> +""" type Model model # For global variables nodes # For nodes @@ -25,7 +27,7 @@ end -@doc """ +""" Initialize empty model. Parameters @@ -35,7 +37,7 @@ None Returns ------- New model struct -""" -> +""" function new_model() return Model(Dict(), Dict(), Dict(), Dict(), Dict()) end @@ -45,7 +47,7 @@ end -@doc """Get field from model. +"""Get field from model. Parameters ---------- @@ -61,7 +63,7 @@ create_if_doesnt_exist : bool, optional Raises ------ Error, if field not found and create_if_doesnt_exist == false -""" -> +""" function get_field(field_type, field_name; create_if_doesnt_exist=false) if !(field_name in keys(field_type)) if create_if_doesnt_exist @@ -77,7 +79,7 @@ end -@doc """Add new elements to model. +"""Add new elements to model. Parameters ---------- list of dicts, dict = {element_type => elcode, elids => [node ids..]} @@ -96,7 +98,7 @@ In dict key means element id >>> elements = Dict(1 => el1, 2 => el2) >>> add_elements(m, elements) -""" -> +""" function add_elements(model, elements) elfield = get_field(model.elements, "connectivity"; create_if_doesnt_exist=true) eltyfield = get_field(model.elements, "element_type"; create_if_doesnt_exist=true) @@ -110,7 +112,7 @@ end -@doc """Get subset of elements from model. +"""Get subset of elements from model. Parameters ---------- element_ids : list of ints @@ -120,7 +122,7 @@ Returns ------- Dict {element_type = XXX, node_ids = [a, b, c, d, e, ..., n]} -""" -> +""" function get_elements(model, element_ids) eltyfield = get_field(model.elements, "element_type") elfield = get_field(model.elements, "connectivity") diff --git a/src/abaqus_reader.jl b/src/abaqus_reader.jl index ce292d9..4791762 100644 --- a/src/abaqus_reader.jl +++ b/src/abaqus_reader.jl @@ -12,9 +12,9 @@ eldims = Dict({"C3D10" => 10}) global handlers = Dict() -@doc """ +""" Register new handler for parser -""" -> +""" function add_handler(section, function_name) handlers[section] = function_name end diff --git a/src/elasticity_solver.jl b/src/elasticity_solver.jl index 9881cdd..47245f0 100644 --- a/src/elasticity_solver.jl +++ b/src/elasticity_solver.jl @@ -13,7 +13,7 @@ VERSION < v"0.4-" && using Docile # solver. -@doc """ +""" Interpolate field variable using basis functions f for point ip. This function tries to be as general as possible and allows interpolating lot of different fields. @@ -26,7 +26,7 @@ basis :: Function Basis functions ip :: Array{Number, 1} Point to interpolate -""" -> +""" function interpolate{T<:Real}(field::Array{T,1}, basis::Function, ip) result = dot(field, basis(ip)) return result @@ -62,9 +62,9 @@ end -@doc """ +""" Calculate local tangent stiffness matrix and residual force vector R = T - F -""" -> +""" function calc_local_matrices!(X, u, R, Kt, N, dNdchi, lambda_, mu_, ipoints, iweights) dim, nnodes = size(X) I = eye(dim) @@ -112,7 +112,7 @@ function calc_local_matrices!(X, u, R, Kt, N, dNdchi, lambda_, mu_, ipoints, iwe end -@doc """ +""" Assemble global stiffness matrix to I,J,V ready for sparse format Parameters @@ -126,7 +126,7 @@ Notes ----- eldofs can also be node ids for convenience. In that case dimension is calculated and eldofs are "extended" to problem dimension. -""" -> +""" function assemble!(ke, eldofs_, I, J, V) n, m = size(ke) dim = round(Int, n/length(eldofs_)) @@ -153,7 +153,7 @@ function assemble!(ke, eldofs_, I, J, V) end -@doc """ +""" Assemble global RHS to I,V ready for sparse format Parameters @@ -167,7 +167,7 @@ Notes ----- eldofs can also be node ids for convenience. In that case dimension is calculated and eldofs are "extended" to problem dimension. -""" -> +""" function assemble!(fe, eldofs_, I, V) n = length(fe) dim = round(Int, n/length(eldofs_)) @@ -189,7 +189,7 @@ function assemble!(fe, eldofs_, I, V) end -@doc """ +""" Eliminate Dirichlet boundary conditions from matrix Parameters @@ -214,7 +214,7 @@ Raises Exception, if displacement boundary conditions given, i.e. DX=2 for some node, for example. -""" -> +""" function eliminate_boundary_conditions(dirichletbc, I, J, V) if any(dirichletbc .> 0) throw("displacement boundary condition not supported") @@ -229,7 +229,7 @@ function eliminate_boundary_conditions(dirichletbc, I, J, V) return findnz(A) end -@doc """ +""" Eliminate Dirichlet boundary conditions from vector Parameters @@ -253,7 +253,7 @@ Raises ------ Exception, if displacement boundary conditions given, i.e. DX=2 for some node, for example. -""" -> +""" function eliminate_boundary_conditions(dirichletbc, I, V) if any(dirichletbc .> 0) throw("displacement boundary condition not supported") @@ -272,9 +272,9 @@ end -@doc """ +""" Solve one increment of elasticity problem -""" -> +""" function solve_elasticity_increment!(X, u, du, elmap, nodalloads, dirichletbc, lambda, mu, N, dNdchi, ipoints, iweights) diff --git a/src/interfaces.jl b/src/interfaces.jl index 642fa39..d5dd11c 100644 --- a/src/interfaces.jl +++ b/src/interfaces.jl @@ -11,14 +11,14 @@ VERSION < v"0.4-" && using Docile #using JuliaFEM.elasticity_solver export solve_elasticity_interface! -@doc """ +""" This is generic interface that reads data from data model, solves elasticity problem and updates model. Parameters ---------- model : to be defined -""" -> +""" function solve_elasticity_interface!() return 0 end diff --git a/src/shape_functions.jl b/src/shape_functions.jl index 968f7f3..e49c6c3 100644 --- a/src/shape_functions.jl +++ b/src/shape_functions.jl @@ -1,13 +1,3 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -function C2D4(xx) - theta = 1 - chi = xx[1] - eta = xx[2] - return = [(1 - chi) * (1 - eta), - chi * (1 - eta), - chi * eta, - (1 - chi) * eta] -end -