{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import theme\n", "theme.load_style()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "import sympy as sp\n", "sp.init_printing()\n", "%pylab inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Lesson 16: The Method of Manufactured Solutions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\"Creative\n", "\n", "This lecture by Tim Fuller is licensed under the\n", "Creative Commons Attribution 4.0 International License. All code examples are also licensed under the [MIT license](http://opensource.org/licenses/MIT)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Validation and Verification" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Validation Question\n", "How do we know if the equations we have been solving are the right equations?\n", "\n", "* Compare with data or with more sophisticated models.\n", "\n", "## Verification Question\n", "How do we know that we are solving the equations right?\n", "\n", "* Compare the approximate solution with an exact solution whenever possible. Otherwise, perform a variety of \"sanity checks\". One such sanity check is the method of manufactured solutions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Method of Manufactured Solutions (MMS)\n", "\n", "The method of manufactured solutions is a straightforward, yet powerful, method to verify that the finite element equations are being solved correctly. In the MMS, we start with a \"manufactured solution\" $u$ and determine the form of the sourcing functions required to produce the solution by applying the governing equations to it. The sourcing functions are then applied to the code and solution checked.\n", "\n", "For example, consider the following boundary value problem\n", "\n", "$$\n", "\\frac{d\\sigma}{dx} + b = 0\n", "$$\n", "$$u(0) = 0$$\n", "$$\\sigma(L)=0$$\n", "\n", "Let $\\sigma = E\\epsilon$ (linear elasticity) and suppose $u = f(x)$. Applying the governing equations to $u$ we get\n", "\n", "$$\n", "\\frac{d\\sigma}{dx} + b = E\\frac{d^2f}{dx^2} + b = 0 \n", "\\quad \\longrightarrow \\quad b = -E\\frac{d^2f}{dx^2}\n", "$$\n", "\n", "Enforcing boundary conditions supplies additional constraints on $f$:\n", "\n", "$$f(0) = 0$$\n", "$$E\\frac{df}{dx}=0$$\n", "\n", "Other constraints (as needed) can be imposed on $f$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example\n", "\n", "Let $u = u_0 + u_1 x + u_2 x^2 + u_2 x^3$. The boundary conditions require that\n", "\n", "$$u(0) = u_0 = 0$$\n", "$$\\sigma(L) = E\\frac{du}{dx} = \n", "E\\frac{d}{dx}\\left(u_0 + u_1 x + u_2 x^2 + u_2 x^3\\right) = 0$$\n", "\n", "We have two equations, but three unknowns. We choose the third equation by arbitrarily scaling the displacement at L to 1\n", "\n", "$$u(L) = u_0 + u_1 L + u_2 L^2 + u_3 L^3 = 1$$\n", "\n", "We solve the preceding equations below using sympy" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAAyBAMAAACpAxLYAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMna7q2ZEmYki\n71TRS9i1AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAK70lEQVRoBeVZfYhcVxU/b752Pt7MDgkxaA07\n3VprgrprLVVqbAb8Ryy420BoSawZ/KgiSrYxsirVPCKotIYdFEH/KB38aGpC7DS4RBN0n1ZLEWU3\nsVbFGka0/wRDNh+m2pau55z78e59774xiTOJsAf23XN/59xzzz333nPvnQVIoZ+k4KsIXt93rJVO\nX/GqEJYa/YZ5tJ9wtcj29Rlo7YQS3qiY61ZW5n90nfoudNM73tCWsqfPpStdG4m3Eb5zbXpK9OKd\nSkAa+ITkXnPgukco04X7tV/XmDmT2l9NxyWvuVTlIQvmwiF30Md8NnWbZTuqWSxCX1X4UErn8fqp\nRFfFIAENFvBmpL3aK2mGo3mzI5RvpLUYCO48Xs/ftb0N3rYHvzmv+tiqmKGVX1eW71BMvLxdA3aE\nHqlrwVCYe5JWvfMBnIIb4EvTt0ph6URSa8BIXi2iuZQB5/6te7Qi5D2r8eEwmSBh11upw+H6k7BQ\n/4eUbWgmlAYOvElaHJt2my4vadyK0EhH44IpxOquamXHAy5YYv6KQeA7zoULAHPtOrxRG/m85jQz\nfzRlrrXGf2V+ZfixAovSXqbhblid0bgVoam4H2fjgG4XMY/DY1ElwY20LWjWqnHlJEYoAPinkuSS\nUSwH5ZYSX235OavhiIxAPiVVTwVa3YrQJg0Lxlu4jAgdgrlYM7P6C7MCkJgDADzL5tqQuwieUM10\n7CZYK3eLywnwygC/YemrafD1xFhi2BnquhkhPx7Q/NhlRAj6Xol/qHp6eDvZykyrui5xnR6Gv+aX\n4W4BLYZaFDHlRsRfFVfuqWbCkZtl9TkF26VxBTEjVJ601eD9lxWhfvndUyaLTa+F1otLsS4ACl3v\nVOnVwrIfCNF3ExoI/M/Ze62yKh3ZL+tGKJQGlVHgyl+48FktqcYnuHs5EfK2xptpg7g/erJSCKGB\nrB+dolIAcOxA0/vM3TvUZfX3WhAx6+SUvzeCUrhi0y3Qa7kgHFlsC72FMNLP77o0O7ulgYD3UoSa\n3GITa79+7CP4/crKLQEUw0SE5DvcO7QddYt7z3N445naP3g/WSKiqWNb5dv8GQL+TJ9+5F+0pNwY\nI90i1Mev6pTqkoTLAK+b/x7AHoVaJa9l05FqV8gxBUY02kD7LayXUtITZVy8kC9Qo/P4tw7iEVLv\n8EzP24wKZBFpqseF/qyHiloqPyaQbHlbDrF8N3/7fOI5mRoD5F6lb6aOH9kp1QUpl+GtsDaEjIKt\nshJQ1XBkpEEA+j4tSv4+ihWvg2yNuzMkkqUMMFKH6hIG8WXkTx984nEper0o1Tscd/ERRB7tosVn\nYTGQWrLA+9FfBMtTx7Zg34WQsJ1CACAtqqouY8mQG29o+v8ihR30oU4FeSGX0mXKcZUG5JVYSoUq\npyHTETwZmMY6ouTvQhPHM41s3l7IWoXeItVJGEGxXASLdSl8QZTqPfcxgLMo4tW2GfYrLan8O4C9\nTebzARZsqxK85wRBZ0P6IkmLomJ8sy2jIhtnwyI3fp5E3Cnr1Ggw2mUozAD+ZuEtMYjrQEhF7V4q\nTEfUm3V0RqpTMQHQhRCZ8itUTRKlp8w5jlBmmcTFXV+UWnI8KvO/iFPZY4sAd+14UCqpYqJOEfo2\nVrMEsa01gJsASe/ItAiNdkhNEzfOzR8NESlNEozDkCRjIF0GHGwJdwd3gxoqQqU2Vril6YjYtrhn\nW9IYFZeApzQ9QjI9Zc8lN7scD7/DcSViLGnPihxhdKHZLfX8SaxsI4ATx9MAH6TKWEBfpNQIoV2D\nzKxT6ZAg6lTFAGcCXQbcMJRhj5ESkpLe18Bdw0ZNRzzetrgAUaqodH52S8iVEbKXJNWINnq02YWe\nGI94h+OZgllqqiuSVdIOIrjjf76zLaeObWVCeIBUR1UA0iKkQygMm44UAsREMmGhioF0d7ELOXRs\nn2ioIuSvw1OD17IYlHJELofqpFTHonJR5fmR5Qg1OJENASgdRZtdKMgIrdA7XEdIn1iGEcFmWzil\nLSh2qcq2vEM7QqqMMoRMWoRiad90JNvDdkanUYTIZZAR+hvxSErqXQJ4mBHLEYwlUXYJ1ryd6A2c\nDsqQI1RGyHrsIiy3Jl89aLPLdxK1UOO5wG8o8NA87jLe1xgwJu8m6ufWJlfohoTZLstCy5ZOMypC\ncS84UUXGzMbVNpo1OlUxkLelMVzVuDKm2IMoQvBECG9jzLQFlHSRMEKaMCd5eMVBKrt3WU5sTf6V\nC+MOX1NNj4yPv3t8fBNWMbfMBVii+bM9znLFEKtx4h+o/Jfk1Fm2RIQii/GmMBVYkNmYI0SplTqt\njI/fdHJ8vEPa4oc5NE33GI6QKV3siLWMidgYlIrQJBkQRKkVniS+4j7LxC6rdeA34k7ZIV1Jcsbx\nLONL6O1Y1jlZl5WGWZ4GP8QTv85Tx/fTjhJf4S6zGvMuo2HITtUaYpdxOeC2xrn/huxKSWHknFjL\nli2Qu6xqRIiuQ7kZam7fh3J7ngkIBY+T1/sAPkA5Sx18LFK7DC9Bhwm4jx/0ZJEvGkJHf3MdyIc4\nmT/rEGTbGg0IQ1K7TNRQbT5gdowmUpPVmDO10amKAbuMDi9BuZHI1Djol0+zPcsWOM6yLaj2SJN0\ni+pNQBUoreDthonC6r/r+IFJyC5jgq1LmAo5HnqHw1gDsj160KPFCk5agp46fuzTCGZOsgXb1mJP\nqscjtNFfYoleZFyzGldmEDM6lRFilz28BNzJPwEclB2o+AG85cakLR/3IxFtWkHe7pXZ2Q+Lu0TN\njtDfVSj+gKoZzJyT6ydenJ1QceP2ajz4DofsDHjz9zTJ4u4VcjpOEysrmNKhxls9ZmuqLbWVRVnN\nz/B1Hy8opkW7Md4YrU5lDNhleAfADdv2oLWN0mIUobMdgmxb6qf6MZaR3CTekATwY7ekfIYPmUo2\nHxuPLUyrvdkhoMcKU2SRvaj2JF5oSMZRbLaxKAYGLpeicZZBJjTkklX/U52aTsrElVjgtKZK398u\nvf6pS1lg30oXXZlE92FYJC9GX5BeVJbSDT5ki0qBXedaRQ3ZKdUNVDd4JLvopAT5fuqFtZao6/l1\ntRkQhm/eOLEXix3x6qLHZyoVwlSRFqzRXF+mMCPEC6FTbZdE5f3Uk6f/WODUHigojxDTJnsx1gOR\nHdXbx9RQvN9SXHr5yXSRKVFn6l4TjPidclvx/bQQejJzq7hGigPnco6rGHsx2lbnr2OVaTee0Vwa\n47zBOpSnQgE+55AhNNUTOD92s6H6P5/jZ3Z3+6tHK41kW/YiE8g1BPuTGkNAjgibaf8NUicqP5wr\n9UpHqHvyjjAEh5TJaqC4qGQvSq1aQ0C4364B/VH0kfYfxfKykIuH8/FDyqOtihlaecZhWXjx0DG5\n9fMNh86godqSsFiYdFsuycQzgWLzBV9tuvUHh77TYSruBd7Wh07ZQHQxOp3S1R0Ctx67CBUbKfqD\ngp0dxL1QP4ENqlOXHdUHvr/dJAT2Y5c0P+pWHxj62zBpKuFFIW1ek22vFvGfly1vSbNQ7ZDEfuwS\nsrZJ36GRd6fDdMIL7zaH1mChtT1hr+S4fAiJuLhaD2cW+JuFfEhf/vkibjvpxVNhXGfQ9Y9Lg4VO\nquVNKLEfu0L1l6ktBiF4rcOIwwu/6dAbJOS1pbUzaWkIf/UethODHNDQbHl/Sjddm0yXrR5Jn00G\ncGD1xCF9pF9OF+F/hIZ/pPbr/v9CVmr0deMHfaWrQrie8/R/APZc60Kg3k3CAAAAAElFTkSuQmCC\n", "text/latex": [ "$$- E \\left(\\frac{1}{L^{5}} \\left(2 L^{4} + 6.0 L^{3} - 2 \\sqrt{L^{8}}\\right) + \\frac{6 x}{L^{6}} \\left(- L^{3} \\left(L + 2.0\\right) + \\sqrt{L^{8}}\\right)\\right)$$" ], "text/plain": [ " ⎛ ⎛ ____⎞ ⎛ ____⎞⎞\n", " ⎜ ⎜ 4 3 ╱ 8 ⎟ ⎜ 3 ╱ 8 ⎟⎟\n", " ⎜2⋅⎝L + 3.0⋅L - ╲╱ L ⎠ 6⋅x⋅⎝- L ⋅(L + 2.0) + ╲╱ L ⎠⎟\n", "-E⋅⎜───────────────────────── + ──────────────────────────────⎟\n", " ⎜ 5 6 ⎟\n", " ⎝ L L ⎠" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Manufactured solution\n", "uL = 1\n", "x, L, b, s, E = sp.symbols(\"x L b sigma E\")\n", "u0, u1, u2 = sp.symbols(\"u0 u1 u2\")\n", "a = u0 + u1 * x ** 2 + u2 * x ** 3\n", "eps = sp.diff(a, x)\n", "sig = E * eps * (eps + 2.) / (2. * eps + 2.)\n", "coeffs = sp.solve([a.subs(x, 0), sig.subs(x, L), a.subs(x, L) - uL], \n", " (u0, u1, u2))\n", "coeffs = dict(zip((u0, u1, u2), coeffs[0]))\n", "ua = a.subs(coeffs)\n", "eps = ua.diff(x)\n", "sig = E * eps\n", "b = -sig.diff(x)\n", "b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Yuck! We now use the $b$ previously determined as an input to the simple linear FE program below to verify that we get back the manufactured solution $u$." ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true }, "outputs": [], "source": [ "ROOT3 = sqrt(3.)\n", "def wundee(xa, xb, num_elem, E, q):\n", " \"\"\"Solves\n", " \n", " ds du\n", " -- + q = 0, where s = Ee, e = --\n", " dx dx\n", " \n", " subject to\n", " \n", " u(0) = 0, s(L) = 0\n", " \n", " Parameters\n", " ----------\n", " xa, xb : real\n", " Origin and end of bar\n", " num_elem : int\n", " Number of elements\n", " E : real\n", " Young's modulus\n", " q : callable\n", " Body force per unit area\n", "\n", " Returns\n", " -------\n", " X, u : array\n", " Nodal coordinates and displacements\n", "\n", " \"\"\"\n", " # Linear 1D mesh\n", " n = num_elem + 1\n", " X = linspace(xa, xb, n)\n", " connect = array([[i, i+1] for i in range(num_elem)])\n", "\n", " # Global stiffness and force\n", " K = zeros((n, n))\n", " F = zeros(n)\n", "\n", " # Shape functions and derivatives, Jacobian\n", " N = lambda xi: array([1 - xi, 1 + xi]) / 2.\n", " dN = lambda xi: array([-1., 1.]) / 2.\n", " Jac = lambda xi, xc: dot(dN(xi), xc)\n", " x = lambda xi, xc: xc[0] + (xc[1] - xc[0]) / 2. * (1. + xi)\n", "\n", " # assemble global stiffness\n", " for (e, nodes) in enumerate(connect):\n", " np = len(nodes)\n", " xc = X[nodes]\n", " ke = zeros((2, 2))\n", " fe = zeros(2)\n", " for xi in (-1./ROOT3, 1./ROOT3):\n", " # Convert shape function derivatives to derivatives\n", " # wrt global coords\n", " J = Jac(xi, xc) \n", " dNdx = dN(xi) / J\n", " ke += E * J * outer(dNdx, dNdx)\n", " fe += q(x(xi, xc)) * N(xi) * J\n", "\n", " # add contributions to global matrices\n", " for a in range(np):\n", " I = nodes[a]\n", " F[I] += fe[a]\n", " for b in range(np):\n", " J = nodes[b]\n", " K[I, J] += ke[a, b]\n", "\n", " # Apply boundary conditions\n", " K[0,0] = 1e30\n", " F[0] = 0.\n", " u = linalg.solve(K, F)\n", "\n", " return X, u" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEACAYAAAD1KqK3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xlc1NX+x/HXYV8FBFFABSnEJZcUl9y3FHMtc8slzWzT\nbC/Le1tu12uL3azs9iszNTM1l3JfKjX3BcsdVBQUATdUlB2G8/sDJVCEQYHvAJ/n4+HDmfmeGT5z\nHjDvOd/v+Z6v0lojhBBClDUrowsQQghROUkACSGEMIQEkBBCCENIAAkhhDCEBJAQQghDSAAJIYQw\nRJEBpJT6Til1Til1sJA2nyuljiul9iul7i/ZEoUQQlRE5oyAZgGht9uolHoIuFdrHQQ8BXxVQrUJ\nIYSowIoMIK31FuByIU36AnOut90FuCulqpdMeUIIISqqkjgG5AfE5Ll/BqhZAq8rhBCiArMpoddR\nN93Pt75PYmKirPcjhBAVnJub281ZUKiSGAHFArXy3K95/TEhhBDitkoigJYDIwGUUq2BK1rrcyXw\nukIIISqwInfBKaXmAx0BL6VUDPAOYAugtf5aa71aKfWQUioSSAZGF/Z6bm5ud191BRYWFgZASEiI\nwZVYPukr80lfFU9566+U9CROxB7hZNwRTsSGE3P+BKbsrHxtqjh74Ovpj49nbXyu/1/Dsxb2tg4F\nvmZ2tqZWf4hPAF8v2Pw/CPS7dQ9bYmLiHdddZABprYea0Wb8HVcghBCiWLTWnL8cy6GoMA5H7eFk\nXDjZOjt3u0LhV60O9/g2INC3PoG+9XF38bzt6z0xWfPKUGgY+HfAWFkp/veaxtcLmgfn3C9pJTUJ\nQQghRCnKMmVyIvYIh6PCOBwVxoXE+NxtVsqKe3wbcI9fQwJ961PHJxhHe+dbXkNrTUYm2NvlDxNX\nZ/hlCzQMzN++X/uSD528JICEEMJCmbJNHIs5QFjEHxw8uZu0jJTcbU4OrjQIaMZ9dVpQz78pTvYu\nt32dQyc1M5bD8i0wIhT+NTb/9onDwc62tN7F7RkeQFprMjIykCuz5vD39wcgLS3N4Eosm1Kl+81M\nCKNorTl17jhhEX/w17GtXEv9+xiLj2dtGtZpwX11QgioURcrK2uzXvPoKfhiUc7tbQdu3e7jZczf\nk+EBlJGRgY2NDdbW5nVkRefgUPABQZGfyWTC29ub8+fPG12KECXi/OVYwiI2s/fo5ny717zdfWle\nryPN67bH28O30NeIOaf5YjF8NC5/oPRoBW8Mh77toVWDUin/jhgeQFprCR9RbNbW1jg6OhpdhhB3\nJcuUyYETu9h6cC2RZw7lPl7FyYNmwe0JCe5ALe97zB7xV3OHb5bB68M0Xu5/P8fFSTHl2RIv/64Z\nHkBCCFHZXLp6ge2H1rPj8K9cS7kCgJ2NPfcHtSWkXkeCat53291rW/drft4MK7fBxi/At9rfQeNg\nr1jyH42DXZm8jbsmASSEEGUgW2cTceovth5Yy+Hovejr06ZrVK1Fu8Y9aVGvY4Ez1272zrew8c+c\n26t3wJN982/vGlJ+jo9KABnk3Xff5cSJE8ydO7fYz92yZQtjx44lIiKiFCorWKdOnRgxYgRjxoy5\no+e7urpy8OBBAgICSrYwISxcRlY6e8I3sWHvL7nHdqytbGgS1JZ2jUO5x7dBgbvY1u/SONpD+6b5\nt43pA82CoV97eOC+MnkLpUYCyAydOnXiwIEDnD17Fju7khnbFmcWl5WVFZGRkQQG5kzSb9++/R2F\nz5UrV3j55ZdZs2YNycnJ+Pj48MQTT/DGG2+YVa+5NRcUVteuXSt2vUKUZylpSWw5sIbN+1bmzmSr\n6lqNNo160LpBN6o4uxf6/Kh42LwP2jfN//hj3RWPdS+tqsuWBFARoqOj2b17N7Vr12b58uU8+uij\nJfK6xZ12XhLT1F966SVSU1OJiIjAzc2No0ePcujQoaKfWEwyRVpUZpeuXmDTvhVsP7SejMyc0ylq\nVgukW8gjNLn3AayvH9sxmTS7jsAvmyEzCz59If/fTf8O4FLB59mUxGKkFdr3339Pt27dGDFiBHPm\nzMm3bdSoUYwbN47evXtTpUoVWrduzcmTJ3O3v/DCC9SuXRs3NzdCQkLYunVrvuff+KDu1asX06dP\nz7etcePG/PLLL3Ts2BGAJk2a4OrqyqJFi9i0aRO1av29AHlMTAyPPPII3t7eeHl58fzzzxf4XsLC\nwhg6dGjuenzBwcEMGDAgd/v27dtp0aIF7u7utGzZkh07dhT4Ou+++y4jRozIvR8dHY2VlRUmk4lJ\nkyaxZcsWxo8fj6urKxMmTAByRnE3+iYxMZGRI0fi7e1NQEAAkydPzg3Y2bNn065dO1577TWqVq1K\nYGAga9euLbAOISzJuUtnmLtuGv+a8wyb/lpORmYawbWbMO7h93ht6Cc0q9suN3wAjp+Bds/A1B9h\nxnJITc//JbN6VcWwHhX7y5zFj4AmfNb/rl/j8xd+uePnfv/997z33nu0bNmS9957j/Pnz+Pt7Z27\nfeHChaxdu5b777+fxx9/nEmTJjF//nwAWrZsybvvvoubmxvTpk1j4MCBnDp1Knc33o0P3VGjRvHJ\nJ58wfnzOknr79+8nLi6O3r17079/f6ysrDhw4EDuLrhNmzbl/nyTyUTv3r3p1q0b8+bNw8rKKnch\nxZu1bt2aSZMmcfnyZdq2bUtQUFDutkuXLuUG4dChQ/npp5/o1asXJ06cwMPDI9/r3G6Eo5Ri8uTJ\nbN++nREjRvDEE08U2O7555/n2rVrREVFcfHiRbp37567OxBg9+7djB49moSEBL7++mvGjBlDbKxc\n4UNYpgtX4lmzawF7j25B62yUsqJ53fZ0af4wtbwD0Vrzwzp47EGNtfXffzv1/BXdW2rqB+Qcz7Gz\n+E/jkicjoEJs3bqV2NhY+vbtS1BQEA0aNODHH3/M3a6U4pFHHiEkJARra2uGDRvGvn37crcPGzYM\nDw8PrKysePnll0lPT+fo0aO3/Jw+ffpw7NgxTpw4AcDcuXMZMmQINjZF/0bu3r2b+Ph4Pv74Yxwd\nHbG3t6dt27YFtv3iiy8YNmwY06dPp2HDhgQFBeWOLlatWkVwcDDDhg3DysqKIUOGUK9ePZYvX37L\n65izO/B2bUwmEwsXLmTKlCk4Ozvj7+/PK6+8km8yhr+/P2PGjEEpxciRI4mPj5cTToXFSbh6jh9/\n/YLJ348jLOIPlFK0va8Hbz/+FY/3fIVa3jlfGJVSTP0Rdh6+9TXWfqr49AVFp2YqXzhVFhafuXcz\nerlbc+bMoXv37ri6ugIwcOBA5syZw4svvpjbpnr16rm3HR0dSUpKyr0/depUvvvuO+Li4lBKcfXq\nVS5evHjLz3FwcGDQoEHMnTuXd955hwULFrBkyRKzaoyJicHf3x8rq6K/Szg4OPDmm2/y5ptvcu3a\nNT744AMGDhzI6dOniYuLo3bt2vna+/v7ExcXZ1YdN7vdKOnixYtkZmbmLjkEULt27XwjnBo1auTe\ndnJyAiApKSnfyFMIo1y+doH1uxez48hvZGebsFJWtKzXFW+Px9j0Z1UCaoDnTVedefUxkEOjt7L4\nADJKamoqP/30E9nZ2fj4+ACQnp7OlStXOHDgAI0bNy70+Vu2bOHjjz9mw4YNNGzYEICqVavedmTw\n+OOPM3LkSNq2bYuTkxOtWrUyq85atWpx+vRpTCZTsVaUcHV15c0332TKlClER0fj5+fH0qVL87U5\ndeoUPXv2vOW5Li4upKT8vSji2bNn820vbBKCl5cXtra2REdHU79+fQBOnz5NzZo1za5dCCNcS0lk\n/Z5FbD24FpMpC6WsaFGvE6GtBjP1xxp8PC+nXXIatGqY/7kjQiV9CiK74G7jl19+wcbGhvDwcPbv\n38/+/fsJDw+nffv2fP/990Dhu6KuXbuGjY0NXl5eZGRk8K9//YurV6/etv0DDzyAUopXX32VkSNH\n5ttWvXr13N1zN2vZsiU+Pj5MnDiRlJQU0tLS2L59e4Ft33//fcLCwsjIyCAtLY3PPvsMDw8PgoOD\n6dmzJ8eOHWP+/PlkZWWxcOFCIiIi6N279y2v07RpUzZv3kxMTAyJiYlMmTLF7Hqtra0ZNGgQkyZN\nIikpiVOnTvHpp58yfPjw2/aNEEbKzMrg17ClvD/nWdbv3sD5hJo0q9uON4d/xogeL1LN3YcHW0Dt\n6jD+URjUxeiKyw8JoNv4/vvveeKJJ6hZsybe3t54e3tTvXp1xo8fz48//ojJZCrw3Jgb90NDQwkN\nDaVu3boEBATg6OiYbxdXQc8dOXIkBw8evOXD+N133+Xxxx/Hw8ODxYsX53uutbU1K1asIDIyktq1\na1OrVi1++umnAt+TlZUVo0ePplq1avj5+fH777+zatUqnJyc8PT0ZOXKlXzyySd4eXkxdepUVq5c\nSdWqVW95nW7dujF48GAaN25MixYt6NOnT7738sILL7B48WKqVq2ab3flDV988QXOzs4EBgbSvn17\nhg0bxujRo2/bLzKtWxhBa01YxB/8+/txrNj2PWkZKbg6Psj+4x8yquer1Kj690zULs0hagl8/pKi\nbWP5fTWXKovLICQmJub+kJsvyZ2WliYrQF83d+5cZsyYwebNm40upVy4cOECp06dKjeXTTZSebvE\ntNE+/+E3Vu+9Qt2AHwDw9Qqgf7tRBNVswuB/wg/v5Ky7JvJfktvNza1YnSLHgCxESkoKX375Ze5U\nbCFE2Yu9EMfQd86xdX9XAOrU3Mvjod1oWb9T7uKgi/9jZIUVi+yCswDr1q3D29sbHx8fHnvsMaPL\nEaLSSc9MY8W2uXw8fwLxCSlYW2XQscl+Jgx4h9YNu5p94TdRPDICsgA9evTIN31bCFE2tNbsj9zB\nz5u/43JSzikSrw07hG2GLY0CrKgfIIcHSpMEkBCiUjp/OZb5v33Hibi9ANT0DmRgp6ep4xN829VE\nRMmSABJCVCrpmWms372IDX8uY+nGSTSt68arjwXR9r7usqutjEkACSEqjQMndrFk04zc3W1j+x1j\n81/P0q6RjUz3N4AEkBCiwruafJnFm2awLzLnJO28u9v0SC3hYxAJICFEhaW1Zsfh31i2dTbHY+rg\n4RrAyNButG/cM3d3m4SPcWQatoXLex2d4tqyZQv16tUr4YqEKB/OX47ji6X/5Mdfv2LT3r4s2/Qv\n9oZ/TOsGveRYj4WQACpEQEAATk5OuLq64urqSpUqVTh79mzuBdhuPH7j36JFiwp8ncOHD9O9e3c8\nPT3x8PAgJCSENWvWlHi9N4fVnV66W4jyzGTK4tc9S/hg3gtEnjlERmYw+449CkrRt50N1vKpZzFk\nF1whlFKsXLmSLl3yry4YHR0N5CxBYc5lEPr06cO4ceNYvXo1Wmv27NlTIpfYLkhZLK0khKU6c+Ek\n89Z/TuzFaABa1u/Mw+1H80BDKwJqQOfmsrvNksh3gVJ28eJFoqOjGTt2LDY2Ntja2tKmTZt8F42b\nMWMGQUFBeHp60q9fP+Lj4wt8rU6dOjFz5szc+7Nnz6Z9+/YAdOjQASj80t3h4eF06tQJDw8P7rvv\nPlasWJG7rajLiwthyUymLNbuWsjUBa8RezEazyrVea7/uwzv/gLOjlUY3UtJ+Fggiw8gq7Y699/t\ntt/J88xV2IjCnNGGp6cn9957L8OGDWPZsmWcO3cu3/YNGzbw1ltvsWjRIuLj4/H392fIkCEFvlZB\nK0XfcGMB0wMHDnDt2jUGDhyYb3tmZiZ9+vQhNDSUCxcu5F4d9dixY7ltFi5cyLvvvsvly5e59957\nmTRpUpHvTwijxSfE8OlPE1m9cz6pafYcPP4hLw78jHr+TY0uTRTB4gPISFpr+vfvj4eHBx4eHjzy\nyCP5tnt5eeVu8/DwKPBy20opNm7cSEBAAK+88gq+vr507NiRyMhIAObNm8eYMWNo2rQpdnZ2TJky\nhR07dnD69OkSfS87d+4kOTmZiRMnYmNjQ+fOnenduzfz58/PbVPY5cWFsDTZ2SZ+3/szH89/mdPn\nI/FwrcZLg97Ay60u/11gb3R5wgwWfwwoe1vhw+bbbS/qeeZQSrFs2bJbjgHdkJCQYNYxID8/P774\n4gsAzpw5w1NPPcXIkSPZvn078fHx+ZbId3Z2xtPTk9jY2FsukX034uLi8u2Og/yX3FZKFXp5cSEs\nyfnLccz79XOi4nMm2bRu2I2H2z+Bo70TMyZqufx1OWHxAVTR1KxZk+eeey531WtfX9/cSQ0AycnJ\nJCQk4Ofnd8tznZ2dSU5Ozr1/86WwC+Pr60tMTAxa/33S3alTp2SatihXtNZsObCGZVtnk5mVgaIO\nT/UdRsM6f3+Jc3eV9CkvZBfcXTDnGNCVK1d45513OHHiBNnZ2Vy8eJHvvvuOBx54AIChQ4cya9Ys\n9u/fT3p6Om+99RatW7cucPTTtGlTli5dSmpqKpGRkfkmJEDhl8Ju1aoVTk5OfPTRR2RmZrJp0yZW\nrlyZe7xJZs8JS3c1+QpfL3ufxZu+ISMzg6tJ45nxyycciWpudGniDkkA3QV3d/d85wFNmzbtljZ2\ndnacOnWKbt264ebmRqNGjXB0dGT27NkAdO3alffff58BAwbg6+tLVFQUCxYsyH1+3kkHL730EnZ2\ndlSvXp3Ro0czfPjwfNsLu3S3nZ0dK1asYM2aNVSrVo3x48czd+5c6tatm/tz5FLYwlIdjgrjg3kv\ncOTUnzjZu+Bg+1++X92VtAzFjkNGVyfuVJGX5FZKhQLTAGvgW631hzdt9wJ+AGqQs0tvqtZ6dt42\nckluURrkktzmK6+X5M7ISmfZljlsObAagLo1GzGs+wtYKU86jYM3RsCQbiX/Ram89pcRSu2S3Eop\na2A60A2IBfYopZZrrcPzNBsP/KW1fvN6GB1VSv2gtc4qTiFCCJHXmQsnmbP2v5y7dAZrKxt6txlO\n52Z9sVI5O272ztJYWckovTwrahJCSyBSax0NoJRaAPQD8gZQPND4+u0qQIKEjxDiTmXrbDb9tYIV\n2+diMmXh7FCPpvdOoGtz33ztJHzKv6ICyA+IyXP/DNDqpjYzgA1KqTjAFRhU2AvefKVBf39/2QUn\n7opcvdJ8lt5XqRnJbDu+jLgrOatw1K3RHFebvjz3kRcOmQfw88oo03osvb8sQVBQ0B0/t6gAMmdq\n1FvAPq11J6XUPcCvSqkmWutrd1yVEKLSOXslmi3HfiE1Mwl7G0fa3NuHWp51gUwmDjqNrY3M1Kxo\nigqgWCDv2Yu1yBkF5dUGmAygtT6hlIoCgoECvzrcfFAvLS2tGOUKcSs5UFw0Sz6onp1tYu3un/j1\n8E9oNIE+DXi858t4uHrltinrsi25vyxN3kkIxVXUNOwwIEgpFaCUsgMGA8tvahNBziQFlFLVyQkf\ns1exVEphMpnMr1gIwGQykZqaanQZ4i4lJl1i+s/vsHbXQgDcnF5i5dZ/42jnaXBloiwUOgLSWmcp\npcYD68iZhj1Tax2ulHr6+vavgf8As5RS+8kJtNe11pfMLcDOzo6MjAwyMzPv+E1UJNeu5ey5dHV1\nNbgSy6aU4vz580aXIe5C+Km/mLtuGkmpibg6uXPu4kd88VM1AGathmcfNrhAUeqKXIpHa70GWHPT\nY1/nuX0R6HOnBSilsLeXhQNvOHQo56w6GfqLisqUbWLVjh/5LWwJAMG1mjCix0t8v8YNWxv44Fl4\nup/BRYoyIWvBCSHKTGLyJWav+YQTsYdRyoperYfSrcUArJQV4wZoureEurVlenVlIQEkhCgTx88c\nYvaaqVxLuYKDXQ1GhU6gQZ0GuduVUtQtuQXgRTkga8EJIUqV1prfwpby5dK3uZZyhaCajUhNncaX\nS+sbXZowmIyAhBClJiU9iXnrP+fgyd0AdAsZQK8HHiMpxYqR/4KkFI2Lk+xyq6wkgIQQpeLMhZN8\nt+ojLiaexdHOieE9XqRRYEsA3Fxg2UcGFygMJwEkhChxu478zk8bvibTlIHJ1JlGgSNpFOhhdFnC\nwkgACSFKTJYpk6V/zGTrwbVorbhw+XWWbGzND2sUbRtp6vjK7jbxN5mEIIQoEYnJl5i+5G22HlyL\ntbUN/do9z97wBzCZFE/1A79qRlcoLI2MgIQQdy0qPoKZqz7kavJl3Fw8ebLXG/jXqMv89zTnL8ND\nbWTkI24lASSEuCvbDq5j8aYZmLKzuMevIaN7vkYVZ3cAQupL8IjbkwASQtyRzKxMFm/6hh2Hf0Vr\n2H1oGs/1q0kVZ/lYEeaR3xQhRLFdSUpg5qoPOXX2GLbWdgzu+izBNf2Z+BX88qHR1YnyQgJICFEs\nUfFHmbnyA66mXMbDtRpP9p5ILe97aFFPMyLU6OpEeSIBJIQw264jG1iw4X8kXKlGk6CmPNXnJVyd\n3ICctdzcXAwuUJQrEkBCiCKZsk0s3zqHjX8tJzy6M1v3PYuLvTWuTtZGlybKMQkgIUShUtKSmL1m\nKhGn93HmfFN+3z0BgKvJYDJprK1lppu4MxJAQojbOnsphhkrpnDhShzOjlX48NmBVHGADvfD6F45\nu92EuFMSQEKIAh2OCmPO2v+SlpGCn1cAY/u8RdUq3sz+p9GViYpCAkgIkY/Wmo1/LeOn3xdxPKYt\nw3okM6z7BOxtHYwuTVQwEkBCiFxZpkx+2vB/7DzyOyhHIqJH4OLgir2t7GoTJU8CSAgBQFLqVWau\n+pATsYextbFj9EPP82w/V6o4G12ZqKgkgIQQxCfE8H/LJnP52lncnKsyts9b1K5+r9FliQpOLscg\nRCV3JPpPJkz7gc8X/hNHu5a8MuRjCR9RJmQEJEQlpbVm8/5VvD3jHFv2TQQgOeV13F3kY0GUDflN\nE6ISMpmyWPzHt2w7uBafaoHY2Wbxj1HWvDlCVjYQZUcCSIhKJjU9mVmrPybi9D5srG15dUg/Pn/B\nBm8PmekmypYEkBCVyJkL55i+ZDop6QdxcXRjbJ83qeNTz+iyRCUlkxCEqCQuXItl3Cc/M2vF43i6\nBfDK4I8kfIShZAQkRCUQffEI244vp3aNLJJSGjGw02Q83eQEH2EsCSAhKjCtNb+GLWHz0aUAtLnv\nQaY93wpra/nTF8aT30IhKqj0jEz++e1q0jJ/AKCZf1eGdH1OVrAWFkMCSIgK6GJiEt2eP8WBE315\nsFUcT/Vwxd+znoSPsCgSQEJUMAlXz/HImwc4cKIbDnbJjOnVG3+3c0aXJcQtJICEqEBOn4vk6+X/\npkFgGgmJ1fh2Yi0eaFSLsDAJIGF5JICEqCAOntzNnDWfkJGVzn11GvHf5+/Fyd7F6LKEuK0izwNS\nSoUqpSKUUseVUm/cpk0npdRfSqlDSqlNJV6lEKJQK7at49sVU8jISqdV/S480/9tCR9h8QoNIKWU\nNTAdCAUaAEOVUvVvauMOfAn00VrfBzxaSrUKIW6SnW3iu5XzGPp2CBeu+PNQ66E89uDz2FjbGl2a\nEEUqahdcSyBSax0NoJRaAPQDwvO0eQxYorU+A6C1vlgKdQohbpKRmc736z7lwImddGx2Fl/PcYS2\nCjK6LCHMVlQA+QExee6fAVrd1CYIsFVKbQRcgc+01nNLrkQhxM2upSTyzYrJnDp7DEd7Z6a90J2g\nmhI+onxRWuvbb1RqABCqtR57/f5woJXW+vk8baYDzYCugBOwA+iltT5+o01iYmLuDzl+PPdhIcQd\niE24yvYT80nNvICzvRtdGwzF3cnL6LJEJRUU9PcXHzc3t2KdaFbUJIRYoFae+7XIGQXlFQOs11qn\naq0TgM1Ak+IUIYQwz+HTCYyZFsQvm0fg4eTDQ41HS/iIcquoEZANcJSc0U0csBsYqrUOz9OmHjkT\nFXoA9sAuYLDW+siNNnlHQG5ubiX8FiqWsLAwAEJCQgyuxPJVtr76PSyMfm8EkpLmga/XefbMrIKP\nl6NZz61sfXW3pL/Ml5iYmHu7uCOgQo8Baa2zlFLjgXWANTBTax2ulHr6+vavtdYRSqm1wAEgG5iR\nN3yEEHdv018rWLb1OwL9xqKzG/HH/3ypWkWuXirKtyJPRNVarwHW3PTY1zfdnwpMLdnShBDZOptf\nNs9i074VoODD55Lp0qwm9nayppso/2QlBCEsVHpmOj+sn8b+yB1YW9nw2IPjaVGvk9FlCVFiJICE\nsEDJqVd5+M0wsrJdCKnvxJO9J1K3VmOjyxKiRMkluYWwMAlXz/Hpojep4fUT+48P4PEeH0j4iApJ\nAkgICxJz/iSfLpzI+cux1A+w48g8GxoG1ja6LCFKheyCE8JCHIn+i1mrPyQ9M42gmo14svdEHO2d\njS5LiFIjIyAhLMCXS/6k6/POXEmypXlwB57p97aEj6jwZAQkhIG01kz8aiefzG9OdrYtycmvMKJH\nY6yUfDcUFZ8EkBAGMWWbWLTx/zgYdY3s7JYM6hrFD+80wUrJOT6icpAAEsIA6ZlpzF4zlcNRYdT3\nt2NUz8MM6iIz3UTlIgEkRBnLeykFJwdXnuoziUDfekaXJUSZkx3NQpShiFNn6fP6VqLjj1PVtRov\nDfpAwkdUWjICEqKMnD4XyazVHxAZ8yIOdmP49/ttcHOuanRZQhhGAkiIMnAkei/frf6YjMw0Jgxe\nyeAuE3BzdjK6LCEMJQEkRCnbdeR35v/2Jdk6mxb1OjG02zhsrG2NLksIw0kACVFKsrM1z32yj4uJ\nv+FbLZsHQwbQu81wlEyzFgKQABKiVGRkmug/8TBrdzbFwS6QRf/ZRa8HuhtdlhAWRWbBCVHCMjLT\n+XTRl/zxVy2srTJ5Z8w5CR8hCiAjICFKUFLqVb5ZPpnYC0d5tMsFQluNYeiDdY0uSwiLJAEkRAlJ\nSDzHV7+8x/krcXi4VuOtEU9Ro2oto8sSwmJJAAlRAlZuO8Nvez8kW8fh5xXAM/3exs1FzvERojBy\nDEiIuxR+6i+mzN3Bwl/HE1CjORMenSzhI4QZZAQkxF3YHb6RH3+bTrN6JuoHuPJs/4k42ss5PkKY\nQwJIiDugtebXPYtZuWMeAN1C+tOnbXe5jo8QxSABJEQxpaaZeGfmetKy5qFQPNJxDB2b9ja6LCHK\nHQkgIYouY/aBAAAWvklEQVTh3OV0Ojx7luMxoYQ+cJT/PN2SpkFtjC5LiHJJAkgIMyWnXqX/xCMc\nj2mFi+Mlnnu4N02D7jW6LCHKLQkgIcyQkHiOr5b9i8b3XuZq0kvM+acfIfUkfIS4GxJAQhQh5vxJ\nvl72PldTLlPHN4APnw3E3cXT6LKEKPckgIQoxNYDB1m2dTLpmWnUrdmIMb0n4mjvbHRZQlQIMmdU\niNv4+Y9t9HixJvEXvWge3IGn+70t4SNECZIRkBA30Vqzfs8iNu77kbZNOuHtPpgRPdrKOT5ClDAJ\nICHyMGWbWLTxa7YfWo9C8Y9RdenQpL3RZQlRIUkACXHdqbNprNr+XyJidmNrbcfI0Jdpcm9ro8sS\nosKSABIC2HX4KqEvm/D1aknv9uE83XcSgb71jC5LiApNAkhUenvCz9FpvAvpGa5UcfbnmX4fUsfH\n1+iyhKjwJIBEpRYVH8GSzf8h0HckNtZe/PqZP94eHkaXJUSlUOS0HqVUqFIqQil1XCn1RiHtWiil\nspRSj5RsiUKUjv2RO5m+5G1S0q7y3IBd7PwmWMJHiDJU6AhIKWUNTAe6AbHAHqXUcq11eAHtPgTW\nAqqUahWiRGit2bx/FUv/mIlG0+a+BxnY+RmsrayNLk2ISqWoXXAtgUitdTSAUmoB0A8Iv6nd88Bi\noEVJFyhEScrW2Tz0ylGydST1AjS9HxjGgy0eRSn53iREWSsqgPyAmDz3zwCt8jZQSvmRE0pdyAkg\nXZIFClFSMrLSmbtuGl7uMazd8Qb/GNWMDk06GF2WEJWW0vr2eaGUGgCEaq3HXr8/HGiltX4+T5tF\nwFSt9S6l1GxghdZ6Sd7XSUxMzP0hx48fL9l3IIQZ0jJT2Bj+ExeuncHW2p72QQOp6RlgdFlClHtB\nQUG5t93c3Iq1K6GoEVAsUCvP/VrkjILyag4suL4LwwvoqZTK1FovL04hQpSWxORLbIxYwNW0SzjZ\nudK1wVA8nL2NLkuISq+oAAoDgpRSAUAcMBgYmreB1jrwxm2l1CxyRkC3DZ+QkJA7rbVSCAsLA6Sf\nzGFOX3255AyT51SlT/s0An0DeLrfPyvlpRTk96p4pL/Ml5iYeMfPLTSAtNZZSqnxwDrAGpiptQ5X\nSj19ffvXd/yThShlr//vJFPn+QNWXLk2ghcGdsLBztHosoQQ1xV5IqrWeg2w5qbHCgwerfXoEqpL\niDumtWbDn8vYH3kSeIEh3Xby/dvdsbGWadZCWBJZCUFUKKZsE0s2zWDrwbUE+8OgLo14olc3mWYt\nhAWSABIVRnpGKrPWTOVI9F5srG0Z3v0FmtVtZ3RZQojbkAASFcLhk5eZMC2MBoF7cXF0ZWyftwj0\nrW90WUKIQkgAiXIv9kIUP/z6EUeiXsTRfgjfvdUBbw9ZzVoISycBJMq1M5eOs3D3VNIz0xg3YBFP\n9J6At0cVo8sSQphBAkiUW+FxewiLWo9GExLckaHdxmNrY2t0WUIIM0kAiXLn/GUTfV8/Rc3qp/Gt\npunZagihrQbLTDchypkirwckhCXZfjCVekOvsftIHTb/+RRt7+1Hz9ZDJHyEKIckgES5cenqBZZv\n+4DUdCt8qx3nwycjuKd6I6PLEkLcIQkgUS6cOnuc/y58nZT0/TzZ73N2f+tKo9rVjC5LCHEX5BiQ\nsGhXkzWRsduYt/5zMk0ZBNVsxJheL+Lk4EJcdKzR5Qkh7oIEkLBYR09pOoxLpm+Hr3G0z6DNfQ/y\naKensLGWmW5CVAQSQMIiZWSls+PIF9SqXoe4i/fx6tD6dGraRyYbCFGBSAAJi5OYfIkZK6Zw+txx\nOofsZVToKzSsI9dlEaKikQASFiE1XfPtcujbPoqZqyZzJSmBqlW8earPJHy9/I0uTwhRCiSAhOGO\nRGmGvA2HTsKSP7bROCiBQJ/6jOn9Bq5O7kaXJ4QoJRJAwnA/rMsJH3fXM/h47aVl/c4M7vKcLKsj\nRAUnASQMlZ6ZRoDvdFrdV4OmQat4pNNAujV/WCYbCFEJSAAJw1y6ep4ZK/5D7MVo2jd14nGZbCBE\npSIBJMrUog2a0+egf4cjzFz1EUmpiVRz82Fs37eoUbWW0eUJIcqQBJAoU60awDMfpXM4+gtcnBKp\nV7spo3q+ipODi9GlCSHKmASQKDMmUxY7w2fycOftODkk0vn+vvRt9zjWVtZGlyaEMIAEkCg1JlPO\n7rY6voprKVeYtfpjImMP4+psw5AuE2jVoIvRJQohDCQBJEpF7AXN8PfgRCws/+gESzf/hytJCVRx\n8mBM74nU8Qk2ukQhhMEkgESJM5k03SbA0dNQtUo6H//4DZ7uCdTxqccTvV7Hzbmq0SUKISyABJAo\ncdbWiinPZvHed3E0r/9PnBwSadcolEc6jpGVrIUQuSSARIm7mnyZ6LMf0a5pODY2NgzqPJ4HGnYz\nuiwhhIWRABJ3bfFGTcemUM1DERUfwcxVH3I1+TLurp482esN/GvUNbpEIYQFkgASd21POMxaBa8P\nW80vW7/DlJ3FPX4NGd3zNao4y2KiQoiCWRldgCj//jk6nQCfVSz+4xtM2Vl0bNqb8Q+/J+EjhCiU\njIBEsWitAXIXCz136QwzV32ItU0MdrYODO06jubB7Y0sUQhRTkgACbMlJmme/hDuD4Y3hsOfx7Yy\n/7fppGemUb1qTcb0ekPWcxNCmE0CSJglOj7n3J6TcfDrHo2P11zCji4FoFnd9gzt+hz2do4GVymE\nKE8kgIRZ/KpB9arg6pRFaJtPCTu6HWsrG/q3H0WHJr3k+j1CiGKTABJmsbVRTHluPyu2f0Zy2iXc\nXTwZ/dDrsqSOEOKOSQCJAsVe0PhVyxnVmExZrNwxj9/3/gxAcO0mjOzxMq5ObkaWKIQo58yahq2U\nClVKRSiljiul3ihg+zCl1H6l1AGl1DalVOOSL1WUleRUTasnYdsBzaWr5/ls8SR+3/szVsqK3m2G\n82z/dyR8hBB3rcgRkFLKGpgOdANigT1KqeVa6/A8zU4CHbTWiUqpUOAboHVpFCxKn7Oj4qvXNMu2\nnGLF9n+Qkp6Eu4snj4e+wj1+DYwuTwhRQZizC64lEKm1jgZQSi0A+gG5AaS13pGn/S6gZgnWKMpY\nZlYmGaY5pGWthCxoWCeE4Q9OwNmxitGlCSEqEHXjxMLbNlDqUaCH1nrs9fvDgVZa6+dv0/5VoK7W\n+qkbjyUmJub+kOPHj5dE3aKEZGQpZq2vweAO53F3MXEt9RKbj/5MQnI8VsqKZv5dqe/bUma5CSEK\nFBQUlHvbzc2tWB8U5oyACk+oPJRSnYEngLbFKUIYI+aCPZNm1yHijDORcY483edndp9cR1Z2Bi72\n7nQIfhgvVz+jyxRCVFDmBFAskPf09lrAmZsbXZ94MAMI1Vpfvt2LhYSEFLfGSiUsLAwom346v10T\ncQb8a2TTtvkqtkeuAOD+oLYM7vosTvYupV7D3SjLvirvpK+KR/rLfImJiXf8XHMCKAwIUkoFAHHA\nYGBo3gZKqdrAUmC41jryjqsRZeqhNoqPx5/m3OWPuJpyBntbBwZ2fpoW9TrJLjchRKkrMoC01llK\nqfHAOsAamKm1DldKPX19+9fA24AH8NX1D65MrXXL0itb3AmtdW6wZJkyWb1jPqfP/4xGE1AjmBE9\nXqSau4/BVQohKguzTkTVWq8B1tz02Nd5bj8JPFmypYmStOuw5uN58NO/NRcS45iz9hPOnD+JUlaE\nthxEj5aDsLayNrpMIUQlIishVBLNgiEuQfPJgjDOXf6EjKx0qlbxZmSPlwj0rW90eUKISkgCqJK4\nlnKBRzv/H9Hn/gQgJLgjAzs/haO9s8GVCSEqKwmgCig7W7MnHFo1VGit2XH4N37e8h3pGak4O1Zh\ncOdnaBrUxugyhRCVnARQBXP+smbUv+HXPbDy40Si4j/jyKmcUU+Te1ozqMszuDrJpbKFEMaTAKpg\nxk2FtTvBzSWTmav+Rw3PP3Gyd+HRTmNpHtxBplcLISyGBFAF8+6YK0TGnuf+4I9wcUqgQUBzhnYd\nh5tLVaNLE0KIfCSAKohsnc2OQ7+yfOsc2t+fgr2dI490GE/rBl1l1COEsEgSQOXY6u2aKs4QVDuO\nBb//jxOxh4Gc1asHdX4aD9dqBlcohBC3JwFUjmVlmxjwVjoDur6FjXUiro5uDOg0lvuD2sqoRwhh\n8SSAyqlTZ49xJPpL2jV1wdrqKq0adKV/+1E4O7gaXZoQQphFAqgcSc/QmHQyq3fMZ8v+1Wg0jYOq\nM6TLuwTXbmJ0eUIIUSwSQOVAUopm/H81x2Mu0bHZqySlXcFKWdGlWX96thqCna290SUKIUSxSQBZ\nuJQ0TbPRGUSescPG2gV/XzdaNvBhYKen8asWYHR5QghxxySALFhKehKrd86niosXVas045HOX/Nk\nn/5yvR4hRIUgAWSBtNbsOrKB5VvncC01kQca2fLykCs83GGSLB4qhKgwJIAszKaD10jKWsSllDgA\nAn3ry+42IUSFJAFkIc5djmX51jnM+r0F7i6t6BLyB33ajpDdbUKICksCyAAx5zTLt4K3B4S2vsq6\n3QvZenAd2dkmGtZJQ5m68o/HB2Nv62B0qUIIUWokgMrYsi2ahyfm3G4QkMCeoxNIy0hBKSva3Pcg\nvo4NcLJzlvARQlR4EkClJCtLc+gkNK2bf/dZ6wYmnB2yqVXjIH7VN5OWkUJ9/2b0azcSX68AwsLC\nDKpYCCHKlgRQKUlJh47j4MwvGldnhSnbRFjEH6zdvZCRvS9ibWXC19Offu3fob7//UaXK4QQZU4C\n6C7FXdCs2AaDuoBHlb9HO1WcFY9110TFZ5Np2saanQs4fyVnZpuPpx89Ww3h/rptsVJWRpUuhBCG\nkgC6C2OmaGatzLnt4gjDevy9LVtn81T/nazauYD4hNMAeLpVp2erITQP7oC1lbUBFQshhOWQADJD\nVpYmKRXcXfMfzwn0BUd7eLAF+Hpdb2vKZO/RLWz485fc4KnqWo0eLQfRsn5nrK2ly4UQAiSAzPLl\nUjgcBd+8kf/x8QPgpcHg5KBITU9hw5/r2fjXChKTEgBwd/HkwRaP8kDDbthY2xpQuRBCWC4JoDzO\nXdKER0OnZvlHOr3awNJNt7Z3c1FcTb7M8m0r2XZgDakZKQD4eNama/OHaVa3nQSPEELchgQQcOGy\npv9E2HkYPFzh7AqNjc3fIXRvTcWmL3W+58ReiGbz/lXsjtiIyZSV086vIV2bP0yDgOayeoEQQhSh\n0gWQyaSxsiJfQHi5w+lzYGcLrRtCwlWoXjX/85RSZGZl8NfxbWw9uJbo+KM5j6Nock9ruoY8QkCN\numX5VoQQolyrdAHU/ln46jVoEvT3Y0opVnykubcmuDjdOnK5cCWebQfXsevI7ySnXQPAwc6JlvU7\n0aFJL7w9/MqqfCGEqDAqbABdvKLJ1uDtkT9QWjWE38LyBxDcumJBZlYGh6PC2H5oPRGn9+U+XrNa\nIO0ah9I8uIMslyOEEHehwgXQul2a/8yBbQfhjeEw+en826c8A/Z2BT83W2cTeeYwYUf/YP/x7bmT\nCmyt7bi/blvaNe6Jf/UgOb4jhBAloMIFUFoGbNkPtjZw+dqt2x3s84eH1prYi1GERWxm77EtuVOo\nAfyq1aFlvc60bNAZZwfX0i5dCCEqlXIbQEeiNO98C4sm5w+UB1vAgn9BaOuc5XAKkq2ziTl3gsNR\nYeyL3M7ZSzG526pW8SYkuCPNgzvg41mrVN+DEEJUZuUigK4m61vC5B6/nGM55y/rfMd5nBwUg7re\n+hppGakcPb2fw1F7OBy9l2spV3K3OTu4cn/ddoQEd6SOT7DsYhNCiDJg0QH01c+axRtydqmdXKyp\n6f13MNjbKf6aranmXvBzs7NNxCecJjL2MIejwjgeeyj3fB0AD9dqNKwTwn11WhBcq7EskSOEEGWs\nyE9dpVQoMA2wBr7VWn9YQJvPgZ5ACjBKa/1XSRT38ybY+CfYWENYBNT0zr89wOfvQMrMyuD0ueOc\niAvnZOwRouIjcicRQM75OnV86l0PnRB8PP1lpCOEEAYqNICUUtbAdKAbEAvsUUot11qH52nzEHCv\n1jpIKdUK+ApoXZwivlmm8fGEPu3yB8Irj8GoXtCzdf5LHWRkpnP2UgzxCaeJTzhFdPwxTp0/nm+E\nAznHcwJ961OvdlMaBDTHxbFKccoSQghRiooaAbUEIrXW0QBKqQVAPyA8T5u+wBwArfUupZS7Uqq6\n1vpcQS+otb5l5GGl4Mf10Kfd34+Zsk20aniFxKQEos6eZ/vhUzmBc/E0FxPPosm/NI5C4esVQKBv\nfe7xbUCgb308XL3M6AIhhBBGKCqA/ICYPPfPAK3MaFMTKDCAer58iSnPRZCZlUFmVgYZWem4V8mm\nabCJb1ce40rSJRKTEriacgWtswssysrKmhoefvh41sbHszY1qwVSx6ceTg4uRbwdIYQQlkJprW+/\nUakBQKjWeuz1+8OBVlrr5/O0WQF8oLXedv3+b8DrWus/b7RJTEy8/Q8RQghRIbi5uRXrwHpR14OO\nBfKeDFOLnBFOYW1qXn9MCCGEuK2iAigMCFJKBSil7IDBwPKb2iwHRgIopVoDV253/EcIIYS4odBj\nQFrrLKXUeGAdOdOwZ2qtw5VST1/f/rXWerVS6iGlVCSQDIy++XWKOywTQghR8RV6DEgIIYQoLUXt\ngis2pVSoUipCKXVcKfXGbdp8fn37fqXU/SVdQ3lRVF8ppYZd76MDSqltSqnGRtRpCcz5vbreroVS\nKksp9UhZ1mdJzPwb7KSU+kspdUgptamMS7QYZvwNeiml1iql9l3vq1EGlGkRlFLfKaXOKaUOFtKm\neJ/tWusS+0fObrpIIACwBfYB9W9q8xCw+vrtVsDOkqyhvPwzs68eANyu3w6Vvrp9X+VptwFYCQww\num5L7SvAHTgM1Lx+38voui24r94FptzoJyABsDG6doP6qz1wP3DwNtuL/dle0iOg3BNXtdaZwI0T\nV/PKd+Iq4K6Uql7CdZQHRfaV1nqH1jrx+t1d5MwwrIzM+b0CeB5YDFwoy+IsjDl99RiwRGt9BkBr\nfbGMa7QU5vRVPHBjCZUqQILWOotKSGu9BbhcSJNif7aXdAAVdFLqzdervt2Jq5WNOX2V1xhgdalW\nZLmK7CullB85Hx5fXX+osh7cNOf3KgioqpTaqJQKU0qNKLPqLIs5fTUDaKiUigP2Ay+UUW3lUbE/\n20t6CWhz/+hvnhVXGT8szH7PSqnOwBNA29Irx6KZ01fTgIlaa61y1nqqrDMvzekrW6AZ0BVwAnYo\npXZqrY+XamWWx5y+egvYp7XupJS6B/hVKdVEa13A5S4FxfxsL+kAkhNXzWdOX3F94sEMclakKGz4\nW5GZ01fNgQXX1xn0AnoqpTK11jeft1bRmdNXMcBFrXUqkKqU2gw0ASpbAJnTV22AyQBa6xNKqSgg\nmJxzJEV+xf5sL+ldcHLiqvmK7CulVG1gKTBcax1pQI2Wosi+0loHaq3raK3rkHMc6NlKGD5g3t/g\nMqCdUspaKeVEzgHjI2VcpyUwp68iyLkaANePZwQDJ8u0yvKj2J/tJToC0iV04mplYE5fAW8DHsBX\n17/ZZ2qtWxpVs1HM7CuB2X+DEUqptcABIBuYobWudAFk5u/Vf4BZSqn95Hxhf11rfcmwog2klJoP\ndAS8lFIxwDvk7M694892ORFVCCGEIUr8RFQhhBDCHBJAQgghDCEBJIQQwhASQEIIIQwhASSEEMIQ\nEkBCCCEMIQEkhBDCEBJAQgghDPH/M+b9g/fwPJYAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "xa, xb, n, E_ = 0., 1., 5, 10\n", "xc = linspace(xa, xb)\n", "plot(xc, sp.lambdify(x, ua.subs({E: E_, L: xb}))(xc), \n", " lw=2, label='Analytic Solution')\n", "\n", "q = sp.lambdify(x, b.subs({E: E_, L: xb}))\n", "xp, u = wundee(xa, xb, n, E_, q)\n", "plot(xp, u, '-.', lw=2, label='FE Solution')\n", "legend(loc='best');" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.10" } }, "nbformat": 4, "nbformat_minor": 0 }