From 636e13b2c8562ab07499dffc29d63e45a6c31cd2 Mon Sep 17 00:00:00 2001 From: ratnania Date: Sat, 7 Sep 2019 18:45:49 +0200 Subject: [PATCH 1/2] add test for anidiff 2d --- gelato/tests/test_gelatize_2d.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/gelato/tests/test_gelatize_2d.py b/gelato/tests/test_gelatize_2d.py index fe9ed27..59a3644 100644 --- a/gelato/tests/test_gelatize_2d.py +++ b/gelato/tests/test_gelatize_2d.py @@ -195,6 +195,33 @@ def test_gelatize_2d_7(): # print('> input >>> {0}'.format(expr)) # print('> gelatized >>> {0}'.format(gelatize(expr, degrees=degrees))) +#============================================================================== +def test_gelatize_2d_8(): + + domain = Domain('Omega', dim=2) + + V = ScalarFunctionSpace('V', domain) + u,v = elements_of(V, names='u,v') + + bx = Constant('bx') + by = Constant('by') + b = Tuple(bx, by) + + expr = integral(domain, dot(b, grad(v)) * dot(b, grad(u))) + expr = BilinearForm((u,v), expr) + + nx, ny = symbols('nx ny', integer=True) + px, py = symbols('px py', integer=True) + tx, ty = symbols('tx ty') + + expected = bx**2*nx*Mass(py,ty)*Stiffness(px,tx)/ny + 2*bx*by*Advection(px,tx)*Advection(py,ty) + by**2*ny*Mass(px,tx)*Stiffness(py,ty)/nx + assert(gelatize(expr) == expected) + + + degrees = [3, 3] + print('> gelatized >>> {0}'.format(gelatize(expr, degrees=degrees))) + + #============================================================================== # CLEAN UP SYMPY NAMESPACE #============================================================================== From 7a43a00564aec9e1f486b2cdffa64d02869f7df5 Mon Sep 17 00:00:00 2001 From: ratnania Date: Sat, 7 Sep 2019 18:46:02 +0200 Subject: [PATCH 2/2] add new notebooks --- ...nisotropic_Diffusion_2d_non_coercive.ipynb | 197 +++++++++++++----- notebooks/psydac_example_1.ipynb | 184 ++++++++++++++++ notebooks/psydac_example_2.ipynb | 167 +++++++++++++++ 3 files changed, 495 insertions(+), 53 deletions(-) create mode 100644 notebooks/psydac_example_1.ipynb create mode 100644 notebooks/psydac_example_2.ipynb diff --git a/notebooks/Anisotropic_Diffusion_2d_non_coercive.ipynb b/notebooks/Anisotropic_Diffusion_2d_non_coercive.ipynb index 0b19458..2b2fa85 100644 --- a/notebooks/Anisotropic_Diffusion_2d_non_coercive.ipynb +++ b/notebooks/Anisotropic_Diffusion_2d_non_coercive.ipynb @@ -35,11 +35,13 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "from sympy.core.containers import Tuple\n", + "from sympy import Function\n", + "from IPython.display import Math\n", "\n", "from sympde.core import Constant\n", "from sympde.calculus import grad, dot\n", @@ -49,19 +51,26 @@ "from sympde.expr import BilinearForm\n", "from sympde.expr import integral\n", "\n", - "DIM = 1\n", + "DIM = 2\n", "domain = Domain('\\Omega', dim=DIM)\n", + "x,y = domain.coordinates\n", "\n", "V = ScalarFunctionSpace('V', domain)\n", "\n", "u,v = elements_of(V, names='u,v')\n", "\n", - "bx = Constant('bx')\n", - "by = Constant('by')\n", - "b = Tuple(bx, by)\n", + "#bx = Constant('b_x')\n", + "#by = Constant('b_y')\n", + "#b = Tuple(bx, by)\n", + "\n", + "bx = Function('b_x')\n", + "by = Function('b_y')\n", + "b = Tuple(bx(x,y), by(x,y))\n", "\n", "expr = integral(domain, dot(b, grad(v)) * dot(b, grad(u))) \n", - "a = BilinearForm((u,v), expr)" + "a = BilinearForm((u,v), expr)\n", + "\n", + "mass = BilinearForm((u,v), integral(domain, u*v))" ] }, { @@ -73,7 +82,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 20, "metadata": {}, "outputs": [], "source": [ @@ -91,83 +100,165 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "glt = gelatize(a, degrees=[3,3])" + ] + }, + { + "cell_type": "code", + "execution_count": 22, "metadata": {}, "outputs": [ { - "ename": "TypeError", - "evalue": "Can't multiply sequence by non-integer of type 'dx'", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m/usr/local/lib/python3.6/dist-packages/sympy/core/compatibility.py\u001b[0m in \u001b[0;36mas_int\u001b[0;34m(n)\u001b[0m\n\u001b[1;32m 372\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 373\u001b[0;31m \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 374\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mresult\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/usr/local/lib/python3.6/dist-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m__int__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 223\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_number\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 224\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"can't convert symbols to int\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 225\u001b[0m \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mround\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mTypeError\u001b[0m: can't convert symbols to int", - "\nDuring handling of the above exception, another exception occurred:\n", - "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m/usr/local/lib/python3.6/dist-packages/sympy/core/containers.py\u001b[0m in \u001b[0;36m__mul__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m 87\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 88\u001b[0;31m \u001b[0mn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mas_int\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 89\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/usr/local/lib/python3.6/dist-packages/sympy/core/compatibility.py\u001b[0m in \u001b[0;36mas_int\u001b[0;34m(n)\u001b[0m\n\u001b[1;32m 376\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 377\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'%s is not an integer'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 378\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mValueError\u001b[0m: dx(v) is not an integer", - "\nDuring handling of the above exception, another exception occurred:\n", - "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mglt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgelatize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", - "\u001b[0;32m~/projects/GeLaTo/gelato/expr.py\u001b[0m in \u001b[0;36mgelatize\u001b[0;34m(a, degrees, n_elements, evaluate, mapping, human)\u001b[0m\n\u001b[1;32m 38\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 39\u001b[0m \u001b[0;31m# ... compute tensor form\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 40\u001b[0;31m \u001b[0mexpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mTensorExpr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmapping\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmapping\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 41\u001b[0m \u001b[0;31m# ...\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 42\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/projects/sympde/sympde/expr/evaluation.py\u001b[0m in \u001b[0;36m__new__\u001b[0;34m(cls, *args, **options)\u001b[0m\n\u001b[1;32m 1075\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1076\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0moptions\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'evaluate'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1077\u001b[0;31m \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1078\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1079\u001b[0m \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/projects/sympde/sympde/expr/evaluation.py\u001b[0m in \u001b[0;36meval\u001b[0;34m(cls, *_args, **kwargs)\u001b[0m\n\u001b[1;32m 1167\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1168\u001b[0m \u001b[0;31m# ... # TODO improve\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1169\u001b[0;31m \u001b[0mterminal_expr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mTerminalExpr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1170\u001b[0m \u001b[0;31m# ...\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1171\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/projects/sympde/sympde/expr/evaluation.py\u001b[0m in \u001b[0;36m__new__\u001b[0;34m(cls, *args, **options)\u001b[0m\n\u001b[1;32m 546\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 547\u001b[0m \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_annotate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 548\u001b[0;31m \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 549\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 550\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/projects/sympde/sympde/expr/evaluation.py\u001b[0m in \u001b[0;36meval\u001b[0;34m(cls, *_args, **kwargs)\u001b[0m\n\u001b[1;32m 656\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 657\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 658\u001b[0;31m \u001b[0mnewexpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdim\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdim\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 659\u001b[0m \u001b[0mnewexpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mexpand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnewexpr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 660\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/projects/sympde/sympde/expr/evaluation.py\u001b[0m in \u001b[0;36meval\u001b[0;34m(cls, *_args, **kwargs)\u001b[0m\n\u001b[1;32m 808\u001b[0m \u001b[0mdim\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdomain\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdim\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 809\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 810\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_args\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdim\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdim\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 811\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 812\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mNormalVector\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/projects/sympde/sympde/expr/evaluation.py\u001b[0m in \u001b[0;36meval\u001b[0;34m(cls, *_args, **kwargs)\u001b[0m\n\u001b[1;32m 607\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 608\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mMul\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 609\u001b[0;31m \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdim\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdim\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0ma\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 610\u001b[0m \u001b[0mo\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 611\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/projects/sympde/sympde/expr/evaluation.py\u001b[0m in \u001b[0;36m\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m 607\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 608\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mMul\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 609\u001b[0;31m \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdim\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdim\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0ma\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 610\u001b[0m \u001b[0mo\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 611\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/projects/sympde/sympde/expr/evaluation.py\u001b[0m in \u001b[0;36meval\u001b[0;34m(cls, *_args, **kwargs)\u001b[0m\n\u001b[1;32m 826\u001b[0m \u001b[0mnew\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'{0}_{1}d'\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mop\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdim\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 827\u001b[0m \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdim\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdim\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 828\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mnew\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 829\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 830\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mTrace\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/projects/sympde/sympde/core/algebra.py\u001b[0m in \u001b[0;36m__new__\u001b[0;34m(cls, *args, **options)\u001b[0m\n\u001b[1;32m 108\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 109\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0moptions\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'evaluate'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 110\u001b[0;31m \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 111\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 112\u001b[0m \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/projects/sympde/sympde/core/algebra.py\u001b[0m in \u001b[0;36meval\u001b[0;34m(cls, *_args)\u001b[0m\n\u001b[1;32m 132\u001b[0m \u001b[0mv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_args\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 133\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 134\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mu\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 135\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 136\u001b[0m \u001b[0;32mclass\u001b[0m \u001b[0mDot_2d\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mDotBasic\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/usr/local/lib/python3.6/dist-packages/sympy/core/containers.py\u001b[0m in \u001b[0;36m__mul__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m 88\u001b[0m \u001b[0mn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mas_int\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 89\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 90\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Can't multiply sequence by non-integer of type '%s'\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mtype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 91\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 92\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mTypeError\u001b[0m: Can't multiply sequence by non-integer of type 'dx'" - ] + "data": { + "text/latex": [ + "$\\displaystyle \\frac{n_{x} \\left(- \\frac{\\cos{\\left (\\theta_x \\right )}}{4} - \\frac{2 \\cos{\\left (2 \\theta_x \\right )}}{5} - \\frac{\\cos{\\left (3 \\theta_x \\right )}}{60} + \\frac{2}{3}\\right) \\left(\\frac{397 \\cos{\\left (\\theta_y \\right )}}{840} + \\frac{\\cos{\\left (2 \\theta_y \\right )}}{21} + \\frac{\\cos{\\left (3 \\theta_y \\right )}}{2520} + \\frac{151}{315}\\right) \\operatorname{b_{x}}^{2}{\\left (x,y \\right )}}{n_{y}} + 2 \\left(- \\frac{49 \\sin{\\left (\\theta_x \\right )}}{72} - \\frac{7 \\sin{\\left (2 \\theta_x \\right )}}{45} - \\frac{\\sin{\\left (3 \\theta_x \\right )}}{360}\\right) \\left(- \\frac{49 \\sin{\\left (\\theta_y \\right )}}{72} - \\frac{7 \\sin{\\left (2 \\theta_y \\right )}}{45} - \\frac{\\sin{\\left (3 \\theta_y \\right )}}{360}\\right) \\operatorname{b_{x}}{\\left (x,y \\right )} \\operatorname{b_{y}}{\\left (x,y \\right )} + \\frac{n_{y} \\left(\\frac{397 \\cos{\\left (\\theta_x \\right )}}{840} + \\frac{\\cos{\\left (2 \\theta_x \\right )}}{21} + \\frac{\\cos{\\left (3 \\theta_x \\right )}}{2520} + \\frac{151}{315}\\right) \\left(- \\frac{\\cos{\\left (\\theta_y \\right )}}{4} - \\frac{2 \\cos{\\left (2 \\theta_y \\right )}}{5} - \\frac{\\cos{\\left (3 \\theta_y \\right )}}{60} + \\frac{2}{3}\\right) \\operatorname{b_{y}}^{2}{\\left (x,y \\right )}}{n_{x}}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "p = 3\n", - "\n", - "glt = gelatize(a, p)" + "Math(latex(glt))" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 23, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "# imports from sympy\n", + "from sympy import Symbol\n", + "from sympy import ratsimp\n", + "from sympy import expand\n", + "\n", + "tx = Symbol(\"tx\")\n", + "ty = Symbol(\"ty\")\n", + "nx = Symbol(\"nx\", integer=True)\n", + "ny = Symbol(\"ny\", integer=True)\n", + "N = Symbol(\"N\", integer=True)\n", + "wx = Symbol('\\omega_x')\n", + "wy = Symbol('\\omega_y')" + ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 24, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "p = 3\n", + "degrees = [p,p]\n", + "\n", + "sa = gelatize(a, degrees=degrees)\n", + "sm = gelatize(mass, degrees=degrees)\n", + "\n", + "expr = sa/sm" + ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 25, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "expr = expr.subs({nx: N, ny: N})" + ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 26, "metadata": {}, - "outputs": [], - "source": [] + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\frac{N^{2} \\left(2 \\left(- \\frac{49 \\sin{\\left (\\theta_x \\right )}}{72} - \\frac{7 \\sin{\\left (2 \\theta_x \\right )}}{45} - \\frac{\\sin{\\left (3 \\theta_x \\right )}}{360}\\right) \\left(- \\frac{49 \\sin{\\left (\\theta_y \\right )}}{72} - \\frac{7 \\sin{\\left (2 \\theta_y \\right )}}{45} - \\frac{\\sin{\\left (3 \\theta_y \\right )}}{360}\\right) \\operatorname{b_{x}}{\\left (x,y \\right )} \\operatorname{b_{y}}{\\left (x,y \\right )} + \\left(- \\frac{\\cos{\\left (\\theta_x \\right )}}{4} - \\frac{2 \\cos{\\left (2 \\theta_x \\right )}}{5} - \\frac{\\cos{\\left (3 \\theta_x \\right )}}{60} + \\frac{2}{3}\\right) \\left(\\frac{397 \\cos{\\left (\\theta_y \\right )}}{840} + \\frac{\\cos{\\left (2 \\theta_y \\right )}}{21} + \\frac{\\cos{\\left (3 \\theta_y \\right )}}{2520} + \\frac{151}{315}\\right) \\operatorname{b_{x}}^{2}{\\left (x,y \\right )} + \\left(\\frac{397 \\cos{\\left (\\theta_x \\right )}}{840} + \\frac{\\cos{\\left (2 \\theta_x \\right )}}{21} + \\frac{\\cos{\\left (3 \\theta_x \\right )}}{2520} + \\frac{151}{315}\\right) \\left(- \\frac{\\cos{\\left (\\theta_y \\right )}}{4} - \\frac{2 \\cos{\\left (2 \\theta_y \\right )}}{5} - \\frac{\\cos{\\left (3 \\theta_y \\right )}}{60} + \\frac{2}{3}\\right) \\operatorname{b_{y}}^{2}{\\left (x,y \\right )}\\right)}{\\left(\\frac{397 \\cos{\\left (\\theta_x \\right )}}{840} + \\frac{\\cos{\\left (2 \\theta_x \\right )}}{21} + \\frac{\\cos{\\left (3 \\theta_x \\right )}}{2520} + \\frac{151}{315}\\right) \\left(\\frac{397 \\cos{\\left (\\theta_y \\right )}}{840} + \\frac{\\cos{\\left (2 \\theta_y \\right )}}{21} + \\frac{\\cos{\\left (3 \\theta_y \\right )}}{2520} + \\frac{151}{315}\\right)}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Math(latex(expr))" + ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle - \\frac{42 \\operatorname{b_{x}}^{2}{\\left (x,y \\right )}}{5} - \\frac{42 \\operatorname{b_{y}}^{2}{\\left (x,y \\right )}}{5} - \\frac{- 55824552 \\operatorname{b_{x}}^{2}{\\left (x,y \\right )} \\cos{\\left (\\theta_x \\right )} \\cos{\\left (\\theta_y \\right )} - 5624640 \\operatorname{b_{x}}^{2}{\\left (x,y \\right )} \\cos{\\left (\\theta_x \\right )} \\cos{\\left (2 \\theta_y \\right )} - 46872 \\operatorname{b_{x}}^{2}{\\left (x,y \\right )} \\cos{\\left (\\theta_x \\right )} \\cos{\\left (3 \\theta_y \\right )} - 56621376 \\operatorname{b_{x}}^{2}{\\left (x,y \\right )} \\cos{\\left (\\theta_x \\right )} + 200088 \\operatorname{b_{x}}^{2}{\\left (x,y \\right )} \\cos{\\left (3 \\theta_x \\right )} \\cos{\\left (\\theta_y \\right )} + 20160 \\operatorname{b_{x}}^{2}{\\left (x,y \\right )} \\cos{\\left (3 \\theta_x \\right )} \\cos{\\left (2 \\theta_y \\right )} + 168 \\operatorname{b_{x}}^{2}{\\left (x,y \\right )} \\cos{\\left (3 \\theta_x \\right )} \\cos{\\left (3 \\theta_y \\right )} + 202944 \\operatorname{b_{x}}^{2}{\\left (x,y \\right )} \\cos{\\left (3 \\theta_x \\right )} - 70430976 \\operatorname{b_{x}}^{2}{\\left (x,y \\right )} \\cos{\\left (\\theta_y \\right )} - 7096320 \\operatorname{b_{x}}^{2}{\\left (x,y \\right )} \\cos{\\left (2 \\theta_y \\right )} - 59136 \\operatorname{b_{x}}^{2}{\\left (x,y \\right )} \\cos{\\left (3 \\theta_y \\right )} - 71436288 \\operatorname{b_{x}}^{2}{\\left (x,y \\right )} - 29412250 \\operatorname{b_{x}}{\\left (x,y \\right )} \\operatorname{b_{y}}{\\left (x,y \\right )} \\sin{\\left (\\theta_x \\right )} \\sin{\\left (\\theta_y \\right )} - 6722800 \\operatorname{b_{x}}{\\left (x,y \\right )} \\operatorname{b_{y}}{\\left (x,y \\right )} \\sin{\\left (\\theta_x \\right )} \\sin{\\left (2 \\theta_y \\right )} - 120050 \\operatorname{b_{x}}{\\left (x,y \\right )} \\operatorname{b_{y}}{\\left (x,y \\right )} \\sin{\\left (\\theta_x \\right )} \\sin{\\left (3 \\theta_y \\right )} - 6722800 \\operatorname{b_{x}}{\\left (x,y \\right )} \\operatorname{b_{y}}{\\left (x,y \\right )} \\sin{\\left (2 \\theta_x \\right )} \\sin{\\left (\\theta_y \\right )} - 1536640 \\operatorname{b_{x}}{\\left (x,y \\right )} \\operatorname{b_{y}}{\\left (x,y \\right )} \\sin{\\left (2 \\theta_x \\right )} \\sin{\\left (2 \\theta_y \\right )} - 27440 \\operatorname{b_{x}}{\\left (x,y \\right )} \\operatorname{b_{y}}{\\left (x,y \\right )} \\sin{\\left (2 \\theta_x \\right )} \\sin{\\left (3 \\theta_y \\right )} - 120050 \\operatorname{b_{x}}{\\left (x,y \\right )} \\operatorname{b_{y}}{\\left (x,y \\right )} \\sin{\\left (3 \\theta_x \\right )} \\sin{\\left (\\theta_y \\right )} - 27440 \\operatorname{b_{x}}{\\left (x,y \\right )} \\operatorname{b_{y}}{\\left (x,y \\right )} \\sin{\\left (3 \\theta_x \\right )} \\sin{\\left (2 \\theta_y \\right )} - 490 \\operatorname{b_{x}}{\\left (x,y \\right )} \\operatorname{b_{y}}{\\left (x,y \\right )} \\sin{\\left (3 \\theta_x \\right )} \\sin{\\left (3 \\theta_y \\right )} - 55824552 \\operatorname{b_{y}}^{2}{\\left (x,y \\right )} \\cos{\\left (\\theta_x \\right )} \\cos{\\left (\\theta_y \\right )} + 200088 \\operatorname{b_{y}}^{2}{\\left (x,y \\right )} \\cos{\\left (\\theta_x \\right )} \\cos{\\left (3 \\theta_y \\right )} - 70430976 \\operatorname{b_{y}}^{2}{\\left (x,y \\right )} \\cos{\\left (\\theta_x \\right )} - 5624640 \\operatorname{b_{y}}^{2}{\\left (x,y \\right )} \\cos{\\left (2 \\theta_x \\right )} \\cos{\\left (\\theta_y \\right )} + 20160 \\operatorname{b_{y}}^{2}{\\left (x,y \\right )} \\cos{\\left (2 \\theta_x \\right )} \\cos{\\left (3 \\theta_y \\right )} - 7096320 \\operatorname{b_{y}}^{2}{\\left (x,y \\right )} \\cos{\\left (2 \\theta_x \\right )} - 46872 \\operatorname{b_{y}}^{2}{\\left (x,y \\right )} \\cos{\\left (3 \\theta_x \\right )} \\cos{\\left (\\theta_y \\right )} + 168 \\operatorname{b_{y}}^{2}{\\left (x,y \\right )} \\cos{\\left (3 \\theta_x \\right )} \\cos{\\left (3 \\theta_y \\right )} - 59136 \\operatorname{b_{y}}^{2}{\\left (x,y \\right )} \\cos{\\left (3 \\theta_x \\right )} - 56621376 \\operatorname{b_{y}}^{2}{\\left (x,y \\right )} \\cos{\\left (\\theta_y \\right )} + 202944 \\operatorname{b_{y}}^{2}{\\left (x,y \\right )} \\cos{\\left (3 \\theta_y \\right )} - 71436288 \\operatorname{b_{y}}^{2}{\\left (x,y \\right )}}{7092405 \\cos{\\left (\\theta_x \\right )} \\cos{\\left (\\theta_y \\right )} + 714600 \\cos{\\left (\\theta_x \\right )} \\cos{\\left (2 \\theta_y \\right )} + 5955 \\cos{\\left (\\theta_x \\right )} \\cos{\\left (3 \\theta_y \\right )} + 7193640 \\cos{\\left (\\theta_x \\right )} + 714600 \\cos{\\left (2 \\theta_x \\right )} \\cos{\\left (\\theta_y \\right )} + 72000 \\cos{\\left (2 \\theta_x \\right )} \\cos{\\left (2 \\theta_y \\right )} + 600 \\cos{\\left (2 \\theta_x \\right )} \\cos{\\left (3 \\theta_y \\right )} + 724800 \\cos{\\left (2 \\theta_x \\right )} + 5955 \\cos{\\left (3 \\theta_x \\right )} \\cos{\\left (\\theta_y \\right )} + 600 \\cos{\\left (3 \\theta_x \\right )} \\cos{\\left (2 \\theta_y \\right )} + 5 \\cos{\\left (3 \\theta_x \\right )} \\cos{\\left (3 \\theta_y \\right )} + 6040 \\cos{\\left (3 \\theta_x \\right )} + 7193640 \\cos{\\left (\\theta_y \\right )} + 724800 \\cos{\\left (2 \\theta_y \\right )} + 6040 \\cos{\\left (3 \\theta_y \\right )} + 7296320}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Math(latex(ratsimp(expr/N**2)))" + ] + }, + { + "cell_type": "code", + "execution_count": 28, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "expr = expr.subs({tx: wx / N, ty: wy / N})\n", + "expr = expr.subs({nx: N-p, ny: N-p}) #N = nx + p\n", + "expr = expand(expr)" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\frac{\\frac{\\omega_x^{8} \\operatorname{b_{x}}^{2}{\\left (x,y \\right )}}{30240} + \\frac{\\omega_y^{8} \\operatorname{b_{y}}^{2}{\\left (x,y \\right )}}{30240}}{N^{6}} + O\\left(\\frac{1}{N^{8}}; N\\rightarrow \\infty\\right)$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "e_exact = (wx*bx(x,y) + wy*by(x,y))**2\n", + "\n", + "# here we consider an expansion up to the forth order\n", + "order = 8\n", + "e = expr - e_exact.expand()\n", + "e = e.series(1/N, 0, order)\n", + "Math(latex(e))" + ] }, { "cell_type": "code", diff --git a/notebooks/psydac_example_1.ipynb b/notebooks/psydac_example_1.ipynb new file mode 100644 index 0000000..9228b6b --- /dev/null +++ b/notebooks/psydac_example_1.ipynb @@ -0,0 +1,184 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Convection-Diffusion " + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "from sympde.core import Constant\n", + "from sympde.topology import dx\n", + "from sympde.topology import ScalarFunctionSpace\n", + "from sympde.topology import Line\n", + "from sympde.topology import element_of\n", + "from sympde.expr import BilinearForm\n", + "from sympde.expr import integral\n", + "\n", + "from gelato import GltExpr\n", + "from psydac.api.discretization import discretize\n", + "from scipy.linalg import eig as eig_solver\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [], + "source": [ + "# abstract model\n", + "domain = Line()\n", + "\n", + "V = ScalarFunctionSpace('V', domain)\n", + "\n", + "u = element_of(V, name='u')\n", + "v = element_of(V, name='v')\n", + "\n", + "alpha = Constant('alpha', real=True)\n", + "beta = Constant('beta', real=True)\n", + "\n", + "expr = beta*dx(v)*dx(u) + alpha*dx(u)*v\n", + "a = BilinearForm((v,u), integral(domain, expr))\n", + "\n", + "glt_a = GltExpr(a)" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "# define numer of cells\n", + "ne = 64\n", + "\n", + "# spline degree\n", + "degree = 3\n", + "\n", + "# Peclet number\n", + "# [0.1, 0.25, 0.9, 5]\n", + "Pe = 0.25\n", + "\n", + "alpha = 1.\n", + "h = 1./ne\n", + "beta = alpha * h / (2. * Pe)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "# create the computational domain from a topological domain\n", + "domain_h = discretize(domain, ncells=[ne])\n", + "\n", + "# discrete spaces\n", + "Vh = discretize(V, domain_h, degree=[degree])\n", + "\n", + "# dsicretize the equation using Dirichlet bc\n", + "ah = discretize(a, domain_h, [Vh, Vh])\n", + "\n", + "# dsicretize the glt symbol\n", + "glt_ah = discretize(glt_a, domain_h, [Vh, Vh])" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "# approximate the eigenvalues using the GLT symbol\n", + "eigh = glt_ah.eig(alpha=alpha, beta=beta)\n", + "eigh_r = eigh.real\n", + "eigh_i = eigh.imag" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "# compute eigenvalues using LAPACK eigenvalue solver\n", + "M = ah.assemble(alpha=alpha, beta=beta).tosparse().todense()\n", + "w, v = eig_solver(M)\n", + "eig_r = w.real\n", + "eig_i = w.imag" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEYCAYAAABC0LFYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de5QU9bXo8e/mOYogMRCPiAgqElCQp0dFFBS8QFDQowafaOIhYghGJUSOAQY0F4OCiUvFy9EIYowiEUWDQYXBJMbHDAoqIIiIEUHFFzIgILDvH79qpqann9PdU13d+7NWra5XV/2qYGr371miqhhjjDHx1As6AcYYY/KbBQpjjDEJWaAwxhiTkAUKY4wxCVmgMMYYk5AFCmOMMQlZoDDGGJOQBQpjjDEJWaAoQCKySkT6Bp2OeERko4j0DzgNs0XkNm8+q/fLf7xsX2uu/m1FpIOIrBCR7SIyJtvHN+FmgSKkvAfQtyJS6ZvuAVDVE1R1WcBJDI1U71eqD/1s3f9Y58vhv+04oExVm6rq3ZkeTES+JyLq/b/cKSKbReSXWUhnvPMdJiILRGSHiHwoIpfG2a+xiDzo7bPdC46DovZZJiK7fH9Xa3OV7rCwQBFu56rqIb5pdNAJKmYi0iDoNGTgaGBVbb4Y57q7Alu9/5cHA6OAu0SkdQZpTOReYA9wOHAZMFNEToixXwPgI+BM4FDgN8A8EWkbtd9o399VhxylOTQsUBQg/y9REekuIm96v56eEJHHI0Uu3vZWIvIXEdkqIh/4ix2844wVkbdEZJv33RJv269FZH7Uef8gInd78zeLyPveeVeLyPkJ0qsicpxv2V8sFDd9vnR87J1nrYicHecc3UTkDW+/x4GSWPcr3jFFZC7QBnjG+5U5zvfdX4vIW8AOEWkQIyfQy7sHX4nIQ5F7mOjak5wv8m/b0fv1+7VXJHVe1DXH/feL2m8p0A+4xzvX8Skeu9p1Rx22K1DuW37N+2wU698nEyLSBPgvYIKqVqrqP4GFwBXR+6rqDlUtVdWNqrpfVZ8FPgB6ZDtdBUVVbQrhBGwE+ifahvuj/BC4HmgIXID71XWbt189YDkw0dv3GGAD8H98x3kdaAUcBqwBrvW2HQ3sBJp6y/WBLcAp3vJF3vfqAT8GdgBHxEo7oMBxvuXZwG0ppK8D7tdhK2+5LXBsjPsRuQ83ePfhQuA73304kJ5Ex4x1z711K4CjgINiHG8j8I63/TDg5ch5E117kvP1965jPfA/3vWdBWwHOkTtG/PfL8Y9WgZc482neuxq1x11vIeBSd58c+BBoAKQFP5vPwt8HWd6Nsb+3YCdUevGAs+kcK7DgV3AD6PuxVbgc+/fq2/Qf+9BT5ajCLenvF98kem/o7afgstq362q36nqk7gHR0QvoKWqTlHVPaq6AfhfYLhvn7tVdbOqfgk8g/uliKp+CLwBRHIKZ+H+WF/1tj/hfW+/qj4OvAecnOb1JUvfPqAx0ElEGqr7lfh+jOOcgnv4/d67D/Op/mvXL9Vj+t2tqh+p6rdxtt/jbf8S+C1wSZLjpeIU4BDgdu/eLMU9YKOPHfPfL4vHjnfdXYFficiXuGCvuKLSpMNVq+oQVW0eZxoS4yuHAN9ErdsGNE10HhFpCPwJmKOq7/o2/Rr3o+RIYBYuV3dssnQXMgsU4TYs6o/of6O2twI+jvrj/Mg3fzTQyh9scL8iD/ft84lvfifujzLiUaoeHpd6ywCIyJVeRWHkuCcCLdK8voTpU9X1wC+BUuAzEXlMRFrFOE6s+/BhrBOmcUy/j9LY/qGXnky1Aj5S1f1Rxz4yar9E/36ZHjvmdYtIY6Aj7lf6Yap6rKpeo6pbUjh3bVQCzaLWNcPlgmISkXrAXFwOu1rdnqq+pqrbVXW3qs7B5SoGZzfJ4WKBorBtAY4UEfGtO8o3/xHwQVSwaaqqqf5RPAH0FVdBeT5eoBCRo3G//EcD31fV5rjiF4lznJ3Awb7l/0g1far6qKqejgsqCvwuxvFj3Yc28S4qwTHj/RpO9ivZf8/bAJt9y/GuPdlxNwNHeQ88/7E/TpKWVKR67HjpOxHYoaqb4p1ARK4SkTIRKReRM6O2PSfVW/P5p+diHG4d0EBE2vvWnUScynnv/8GDuB8c/6Wq38VLp0eJ/3+3KFigKGyv4IpSRnuVrEOpXvzzOrDdq5Q8SETqi8iJItIrlYOr6lZcee5DuAf6Gm9TE9wf11YAEbka9/CIZwVwqXf+gbgWKUnTJ67t/1neL9hdwLfA/hjHfwXYC4wRkYYicgFxisGSHPNTXJFEun4uIq1F5DDgFuDxFK492flewwWZcd419QXOBR6rRfqyfexuJGhBJSI9cb/QzwIG4OoTDlDVQVq9NZ9/GhR9PFXdATwJTBGRJiLSGxiKyzHEMhOX4zk3uthMRJqLyP8RkRLvb+Yy4Azgbylee0GyQBFukRYxkWmBf6Oq7sFVYP8UVxF4Oa6sebe3fR8wBFee/AGu8u4BXLPBVD2Kq1w9UOykqquB6bgH9KdAZ1z2PZ7rcQ+ir3FNG59KMX2Ngdu99Z8APwDGRx/cdx+uAr7EVa4/GSctiY45FfiNVww2Ns73Y3kUeB5XEf8+rqI+4bUnO593TecCg7y03gdcGVXWXitZOHZXXA4ynv8CjgfKcNf7de1Te8B1wEHAZ8CfgVGqugoO5FD+x5s/GviZl8ZPfH87l3nHaYj794lUZv8CV8S7LgtpDC1JoW7JFBAReQ24X1UfCjotpjiJyHRggbpmrIhIA1XdG3CyTAIWKAqcV/67Fvfr6DLgfuCYHFYsGpOQiHQA/ohrorwbl1v5NNhUmUTC3JPUpKYDMA9Xb7ABuNCChAmSqq4FegedDpM6y1EYY4xJKNDKbBEZKG6IhPUicnOM7W28JnRvihuGoKjbMhtjTBACy1GISH1c++cBwCZcT9lLvBYzkX1mAW+q6kwR6QQsUtW2iY7bokULbds24S7GGGOiLF++/HNVbRlrW5B1FCcD671hGRCRx3Btn1f79lGqelweSvWOSjG1bduWioqKLCfVGGMKm4jEHK0Agi16OpLqQwBsouYQAaXA5SKyCViEa9Ncg4iMFJEKEanYunVrLtJqjDFFK9873F0CzFbV1rienHOjhhUAQFVnqWpPVe3ZsmXMnJMxxphaCjJQfEz1MXBaU3MsmZ/imnaiqq/g3iGQ7sByxhhjMhBkHUU50F5E2uECxHDcCKR+/wbOBmaLSEdcoEi7bOm7775j06ZN7Nq1K8Mkm9ooKSmhdevWNGzYMOikGGNqIbBAoap7RWQ0sBj30ps/quoqEZkCVKjqQuAm4H9F5AZcxfZVqYxnH23Tpk00bdqUtm3bUn0AUZNrqsoXX3zBpk2baNeuXdDJMcbUQqA9s1V1Ea6S2r9uom9+NVnowblr1y4LEgEREb7//e/jb2RQWuomY0w45HtldtZYkKg7m6MaMUff+8mTa37HAocx+atoAoXJnuhAkO72WGIFD2NMfrBAEbCrrrqK+fPnA/D73/+enTt31sl5N27cyIknJnqXUE19+/aloqKiVoHg669BxE1QNZ8oJ2G5DGPygwWKBOr6QVWXgSLbNm+Gigo3QdV8JKg0bw6qbgL3OWmSy0nECx6WyzAmP1igSCCbD6pbb72VDh06cPrpp3PJJZdw5513Vtt+9913s3nzZvr160e/fv1qfP/mm2+mU6dOdOnShbFjx7J9+3batWvHd9+51/1+8803B5b79u3LDTfcQM+ePenYsSPl5eVccMEFtG/fnt/85jcHjrl3714uu+wyOnbsyIUXXnggSD3++BK6detG586d+clPfsLu3bvZvBm2b4fV3gAr0YGgVSvo2dNNUDXfqlX8e1JaWjN4qFpOwph8Y4GiDpSXl/OXv/yFlStX8txzz8Uci2rMmDG0atWKsrIyysrKqm374osvWLBgAatWreKtt97iN7/5DU2bNqVv37789a9/BeCxxx7jggsuONBXoVGjRlRUVHDttdcydOhQ7r33Xt555x1mz57NF198AcDatWu57rrrWLNmDc2aNeO+++5j165djBlzFY8//jhvv/02e/fuZebMmbRqBU2bQqdOLk2pBIJ4Jk2Kv620NHYRVd++6Z/HGJMdFiiixHtQZfIr9+WXX2bo0KGUlJTQtGlTzj333LS+f+ihh1JSUsJPf/pTnnzySQ4++GAArrnmGh56yL3R9KGHHuLqq68+8J3zzjsPgM6dO3PCCSdwxBFH0LhxY4455hg++sgNsXXUUUfRu7drfXz55Zfzz3/+k7Vr13Lkke04/vjjARgxYgR///vf00pvsuAR615Ggke8XMZLL6WVBGNMFlmgiJKPxSENGjTg9ddf58ILL+TZZ59l4MCBAPTu3ZuNGzeybNky9u3bV61yunHjxgDUq1fvwHxkee9e93pif7PVzz+HbduEVavccqRo6fPPa6YnWSCoTS7DipuMyV8WKOpA7969eeaZZ9i1axeVlZU8++yzMfdr2rQp27dvr7G+srKSbdu2MXjwYO666y5Wrlx5YNuVV17JpZdeWi03kciePVXz//73v3nllVcAeP75Rxky5HQuuKADmzdvpHnz9fTsCc89N5czzzyz2jFqEwhq48wzs5+7M8akzwJFAonK0tPRq1cvzjvvPLp06cKgQYPo3Lkzhx56aI39Ro4cycCBA2tUZm/fvp0hQ4bQpUsXTj/9dGbMmHFg22WXXcZXX33FJZdcklJadu+umu/QoQP33nsvHTt25KuvvmLUqFGUlJQwceJDXHTRRXTu3Jl69epx7bXX1u7CM7RsWeyWUhYojKlbBffO7J49e2p0ZfGaNWvo2LFjQClyKisrOeSQQ9i5cydnnHEGs2bNonv37hkfd/78+Tz99NPMnTs3pf0rKqpaJsWzeXP2cw2Z/huIuEAR+TTGZJeILFfVmE+HQMd6KiYjR45k9erV7Nq1ixEjRmQlSPziF7/gueeeY9GiRQn327y5em/pSBxt1Sp2QKiroqV0ZCt3Z4xJn+UoikwqOYpcyOTfoLQ0dp8WK4YyJnsS5SisjsLkt8GDKW02o3pdxfQZ6KDBFiSMqSMWKApEquMv5WOxUlzTpkHbtjB2LHgV+PdwHXrTTdC/f7BpM6aIWKAoEAUXKKZNgwYN4Ikn4NprXbA45hiuYyb3MQpuvDHoFBpTNCxQmPwybRqUlUGvXjB1KowfD3/6Exx8MHzwAdKuHaO5L+hUGlNULFAEbPDgwXz99de1+m6yEVtTVZshx3OmVy+4+GI3P2+eq8WurIQdO/iEluz/YCO/ZIZ1vjOmDgUaKERkoIisFZH1InJznH0uFpHVIrJKRB7NeaIiv2j9ysrc+hxYtGgRzZs3r9V3azNia16L3ON581ywKCtzQWL/fmjXjv9oodQbdS3TGesqtG2kWWPqRGCBQkTqA/cCg4BOwCUi0ilqn/bAeKC3qp4A/DLnCYv8oo0Ei7Iyt9yrV8aHfuSRRzj55JPp2rUrP/vZz9i3bx9t27blc29ApXhDkb///vsMHDiQHj160KdPH959913AvfRozJgxnHbaaQwdesyBFyANHz78wKiykf3mz5/Pxo0b6dOnD927d6d79+7861//qpHG2bNnM3r06APLQ4YMYdmyZQA8//zznHrqqXTv3p2LLrqIyspKoOYQ6GnxB+bIvX/zTTjpJLj1Vhckjj/ejXE+fjw75jzBTK6FF19M7zzGmNpT1UAm4FRgsW95PDA+ap9pwDXpHLdHjx4abfXq1TXWJbR0qWqLFqoTJrjPpUvT+34Mq1ev1iFDhuiePXtUVXXUqFE6Z84cPfroo3Xr1q36+uuv60knnaTffvutfvPNN3rcccfpHXfcoaqqZ511lq5bt05VVV999VXt16+fqqqOGDFCL7zwQt23b5+Wla3SY489VlVVn3zySb3yyitVVXX37t3aunVr3blzp+7YsUO//fZbVVVdt26dRu7VBx98oCeccIKqqj700EP685///EC6f/SjH2lZWZlu3bpV+/Tpo5WVlaqqevvtt+vkyZP1888/1+OPP17379+vqqpfffVV3OuPKXKvI/d4+nRVEdWGDV2L2JISt83b75dM1+f7/y7Nu2+MSQao0DjP1SB7Zh8JfORb3gT8Z9Q+xwOIyMtAfaBUVf8WfSARGQmMBGjTpk3mKevXD0aNcr9oJ0xwyxlasmQJy5cvp5eXM/n222/5wQ9+cGC7fyjykpKSA0ORV1ZW8q9//YuLLrrowL67fQM2DRs2jHr16tG3byc+/fRTAAYNGsT111/P7t27+dvf/sYZZ5zBQQcdxLZt2xg9ejQrVqygfv36rFu3LuX0v/rqq6xevfrAsOR79uzh1FNPrTYE+pAhQxgyZEhqB5w2zeUg+vWrKmoaNMi1cmrY0I1eeMUVcPXVbtu8eTBvHg3PKmfAC+NSTrcxJnP5PoRHA6A90BdoDfxdRDqrarXaX1WdBcwC1zM747OWlcHMmS5IzJzpHmYZBgtVZcSIEUydOrXa+tmzZyf83v79+2nSpDkrVqyIud0/hLh6PdJKSkro27cvixcv5vHHH2f48OEA3HXXXRx++OGsXLmS/fv3U1JSUuN4DRo0YP/+/QeWd+3adeDYAwYM4M9//nON77z++ussWbKE+fPnc88997B06dKE1wRUFTPNm+fu7aBBMHcutG8Pn30Gv/61u/dXX83swfNYfVY5dzAO6Mcd3miy1jPbmLoRZGX2x8BRvuXW3jq/TcBCVf1OVT8A1uECR+5E6iTmzYMpU6pXrGbg7LPPZv78+Xz22WcAfPnll3z44YcHtscbirxZs2a0atWOJ554AnAPbP8w4/H8+Mc/5qGHHuIf//jHgfdXbNu2jSOOOIJ69eoxd+5c9u3bV+N7bdu2ZcWKFezfv5+PPvqI119/HYBTTjmFl19+mfXr1wOwY8cO1q1bl3AI9IT8OYkrr4RHHoEBA2D9epg4sdq9v+oqmKbj8uodIcYUkyADRTnQXkTaiUgjYDiwMGqfp3C5CUSkBa4oakNuU1Ve9SsXqh5o5eUZHbZTp07cdtttnHPOOXTp0oUBAwawZcuWA9ujhyJv374zX311KBUVcOutf2LGjAc5/viT6NDhBJ5++umk5zvnnHN46aWX6N+/P40aNQLguuuuY86cOZx00km8++67NGnSpMb3evfuTbt27ejUqRNjxow5MHhhy5YtmT17NpdccgldunTh1FNP5d133004BHpM/sprf06iWzfX2/rOO13/ibKyrN17Y0yG4lVe1MUEDMblEt4HbvHWTQHO8+YFmAGsBt4Ghic7ZlYqswOyfft2VVXdsWOH9ujRQ5cvX66qquXlQaYqOw78G/grr5cudZXVjRurNmtWVaG9dKnq76pXWE+aVP3TGJNdJKjMttFj88ill15abSjy8ePHA8GN+Jo1n3zCmi1b6Nitm1suK4Nzz3VvUWrSBBYscOv9dRZR7D0UxuSWvY8iJB59NHZ/wtB2oIs4+GDYurWqOAlcq6a9e2HMmKp1kWKmLLQyM8ZkT9EM4RHmnFPYA4U2bQotW7ocw8SJcP75cNBBVa3K/HUW46qavpaW2juzjckHRREoSkpK+OKLL/IyWKQ7LlNofPIJfPMNqsoXX3xBSfPmruL61ltdbuKpp5K2KistrfnObGvtZEzdK4qip9atW7Np0ya2bt0adFJq+PBD2LYt6FTkwK5drripZUtKmjen9bPPuiaw3brBBl/DNX/LJityMiYvFUWgaNiwIe3atQs6GTF16lTAlbRlZXDOOS4n8cgjrunrjTdW76sS6cwYI0iUllblHuyd2cYEpyiKnvJNwZa9R4+86+8ncfnlVS8bSrF/hP892aG/N8aEWFE0j81nBdXsMzqnMGOGezPd5ZfDc8/FbfoaT0HdG2PyXKLmsZajMNkTPSzH2LGuuOnhh1MeCqVgc1vGhJgFioCFuuw91kuewL1LopbFTdbSyZj8Y4EiYKF+AMZ6ydP557tgMGGCK26KrrMYZ0OEGxM2FihM7fmLmiId6VRT6iMRS3TQDHVuy5gCYoHCZMb/kqeePV2QqOXIu/5WThDy3JYxBcQCRY4V3MMuul6irAz+8Ac4+2yI9S4KK24yJvQsUORY9K/k0PPXS5SVwbBhrlnSLbfU6iVP1srJmPxn/ShyrCD7AkT6S5x0khsDfcGCquKmsjJX1FSLXERB3itjQsL6UdSxgvuVHK/H9ZIl1YcJj2yzoiZjCooFihwouL4A0c1gZ8xwYzddcUX1YcLTEOteWCsnY/JToIFCRAaKyFoRWS8iNyfY779EREUkzO95C68s9LiOFqvuJrSB1JgCF1igEJH6wL3AIKATcImIdIqxX1PgeuC1uk1hdhTMr+RIM9ha9rg2xoRXkDmKk4H1qrpBVfcAjwFDY+x3K/A7YFddJi5bQvkrOdbQHDNmwPTpGfW4Lri6G2OKRJCB4kjgI9/yJm/dASLSHThKVf+a6EAiMlJEKkSkIh9fThQ6seokxo51nepq0eM6ouDqbowpEnlbmS0i9YAZwE3J9lXVWaraU1V7tmzZMveJK3TRQ3NMnFj10iH/dituMqYoBBkoPgaO8i239tZFNAVOBJaJyEbgFGBhPlZoF+QvYv/QHDfeWBUk/NszaAZbMHU3xhSBIANFOdBeRNqJSCNgOLAwslFVt6lqC1Vtq6ptgVeB81Q1f3rTeQqi93WioTlq2QQ2IlYgLcjgakyBCixQqOpeYDSwGFgDzFPVVSIyRUTOCypdRSvLQ3P4FUQgNaaIBVpHoaqLVPV4VT1WVX/rrZuoqgtj7Ns3n3ITBdeCx18v8dvfuouJDM1hdRLGFLW8rczOdwXZgidSL5GFoTkKLpAaU8RsUMAsCO1gdtOmuSIn/4B+w4a5dStXulyEP1jUUmjvjzFFJNGggA3qOjGFKLQteCL1EvPmuWV/vQRUbctCsDDGhJcVPWVBaItTslwvEe8+hDaQGmMAK3oy4DrU3XqrG55jypRaH8aKmIwJL3sfhYmvrMz1k5gwIeP+EsaYwmSBotj4O9ZF3lQ3fjwccoi9ytQYE5MFimLj71hXXu6CxNSpVa2falEvUXDNhI0x1VigSENBPPz8FdiVlS5I+Fs22atMjTFRLFCkoWCGovAP+DdqVK2av9qrTI0pHtbqKQ2hbdWTg451ob0XxpiYrNVTBgqisjaHA/4ZYwqfBYokCqKyNksd6woiaBpj0pa06ElEfqeqv062Ll9Y0VMCWepYBwVwL4wx1WRa9DQgxrpBmSUpnEJdWZtBxzrLMRhT3OIGChEZJSJvAz8Ukbd80wfAW3WXxPwR2gdmpGPdvHkuJ5FmvUSs1l6hDprGmLQkylE8CpwLPO19RqYeqnp5HaTNZMLfA7u8vGqE2GnTsvIiotAGTWNM2uIGClXdBnwEdFPVD33Tl3WXPFNr/pZOkQ50F1/s1kPSjnVWcW2MiUhYR6Gq+4C1ItImFycXkYEislZE1ovIzTG23ygiq70iryUicnQu0lGQ/C2dJk5M+90SBdHayxiTFam8uOh7wCoReR3YEVmpqudlcmIRqQ/ci6ss3wSUi8hCVV3t2+1NoKeq7hSRUcA04MeZnLeo+HtgT5hgLyAyxtRKKq2eJgBDgCnAdN+UqZOB9aq6QVX3AI8BQ/07qGqZqu70Fl8FWmfhvCkL/a/nWrR0sqE5jDHRAhvCQ0QuBAaq6jXe8hXAf6rq6Dj73wN8oqq3xdg2EhgJ0KZNmx4ffvhhltIYwr4CkeE6oPprTh97DJ58MmnxUyiv2RiTsYz6UYjIKSJSLiKVIrJHRPaJyDfZT2bCNFwO9ATuiLVdVWepak9V7dmyZcu6TFr+iVRiP/ZYVZC4+GIYPjzjlk7GmOKUStHTPcAlwHvAQcA1uLqFTH0MHOVbbu2tq0ZE+gO3AOep6u4snDeh0Lf2iVRiP/lk9f4TkeE6YrR0Cv01G2NyKpUhPCpUtaeIvKWqXbx1b6pqt4xOLNIAWAecjQsQ5cClqrrKt083YD6uiOq9VI6bzSE8Ql0MU8vhOkJ9zcaYWst0CI+dItIIWCEi00TkhhS/l5Cq7gVGA4uBNcA8VV0lIlNEJNKi6g7gEOAJEVkhIgszPW9RSKMS23INxphkUslRHA18CjQCbgAOBe5T1fW5T176spmjKC0N4YM0urgpejlKdA4ilNdsjMlYohxFSq2evBzFDwEF1nrNWfNSLkePzWuR1k7l5VUvKYq8FzuyPkb9hBU1GWMg81ZPPwLeB+7GVWyvF5GiHD02r0VaO/mDhH/ZFySs8toYk45Uip7eBYZEippE5Fjgr6r6wzpIX9qKNkcBVcFh1ChXN5HCkB2WozDGQOaV2duj6iM2ANuzkjKTXf4hO0aNihkkLNdgjElXKoGiQkQWichVIjICeAY3LtMFInJBjtNn0pFCa6fod0vY8BzGmGRSCRQluFZPZwJ9ga24jnfn4saAKjih/NVdy5cThfJajTF1KrCxnnIlG3UUoSy3j7R68hc3ea2eSneOi/uWOgsUxhjIQvPYMCm6QJEgQFhzWGNMqjKtzC4KoW0y6n+THVRvFuuT99dhjMlbqTSPre+96S4Uii5HASk1i41ck/W8NsbEkmmO4j0RuUNEOmU5XSZbUmgWG2FBwhiTrlQCxUm4UV4fEJFXRWSkiDTLcboCFbomo3GaxYa2OM0Yk1fSqswWkTOBR4HmuOG/b823wQGLrmd2ioMAhq44zRhTpzId66m+iJwnIguA3+Pel30MruPdoqym1KSvvLx6UIi8uMh7k53lHowxmUqlMnsDUAY8qKr/itp2t6qOyWH60lZ0OYokrBLbGJOKWucoRKQ+MFtVfxodJADyLUgUnWnTava8Litz66NYkDDG1FbCQOE1iy3IYToKQpw+FLNX9bJKbGNM1qRS9HQX0BB4HNgRWa+qb+Q2abVTdEVPSfpQWCW2MSYViYqeGqTw/a7e5xTfOgXOykLCBgJ/AOoDD6jq7VHbGwMPAz2AL4Afq+rGTM9bg28YjANl+WVlcMcdcNRRVfsNHw6PPQabN8Nnn8EPfgB9+sD771dtv+MO+NWv3LJ/3qtcpkED2FOKVYoAABcaSURBVLvXDa8xbVriZYCf/cx9/r//V5UO/zp/H4qzz076/gljTB1Ic2idvD+XqgYy4YLD+7gWVI2AlUCnqH2uA+735ocDjyc7bo8ePTRtS5eqtmihunSpgm95+nTVZs1UmzSpmho3VgXVhg1VRVRHjarap1kz951DD605v3SpWxZxn6rJl5cudd899FA3H2vd0qVu/qCDqu/ngfRvhzEmQ75nSszlPDwXUKHxntfxNmj1B/aPgHHAxMiUyveSHPNUYLFveTwwPmqfxcCp3nwD4HO84rJ4U60CheqBmzuZCTVverNmVQHCHyQGDHCfV1xR9fCeMCH+fCT4tGiR+nIkGMRa16yZ6sEHVwWiGP9BLFAYE5BYf7d5fK6MAgVwP6745yNgEvA2rqlspoHiQlxxU2T5CuCeqH3eAVr7lt8HWsQ41kigAqho06ZN2jdo0iR3JyYzQRV0MhMU3HpVdTc/EiQiU58+1T8nTKjaL9G8/3ipLsdbd9xxNddNn66vfG9QjeSC6plnpn1rjDGZiPV3m6fnyjRQvBX1eQjwj2TfS+G4WQsU/slyFJajMCYvFFmO4jXv81WgFdAYWJ/seykcN3+KnqyOwhiTTQVWR5FKq6dnRaQ5cAfwBq7F0wMpfC+ZcqC9iLQDPsZVVl8atc9CYATwCi4HstS7oOzyDYMxaRJVw2DccYdryRQRq9VT27ZV+0RaPS1Y4Jb985FWT3fe6Vo1gftMtNyvX9Wx/UN0RK8bM8a1ejrtNGv1ZEw+SDS0Trb/RuvgXOkOCtgYKFHVbVk5uchg3PhR9YE/qupvRWQKLrItFJESYC7QDfgSGK6qGxId0/pRWD8KY0z6Mu1HgYicBrSN7C8iqOrDmSZMVRcRNbCgqk70ze8CLsr0PAUreqTYfv1cz+zB87j64erBAuwd2caY2kkaKERkLnAssAKIvOlOcS2hTJDiZDmvKi/nKnXrLEdhjMlUKkN4rMF1hAvF46boip6SsNFjjTGpyPRVqO8A/5HdJJmsSTKCbORtfZMn13G6jDEFI5VA0QJYLSKLRWRhZMp1wkyK4owgS69egOUijDGZSyVQlALDgP+Le7tdZCpYoXq4RprCXXwxTJxYrXLb3pltjMmGtJrHhkE26ihCWQE8caLrSzFhAkyZUmNzKK/JGFNnalVHISL/9D63i8g3vmm7iHyTq8SaWigrc30oJkxwn9F1Fj6WmzDGpCtuoFDV073PpqrazDc1VdVmdZfEuhHaYhp/X4opU6qKoaKChVVqG2NqK5XmsYfFWL1dVb/LTZIyU3RFT2m+tCRU12aMqTOZNo99A9gKrAPe8+Y3isgbItIje8k0tTJunAsS/may/fq59V4z2dDmlowxeSGVQPECMFhVW6jq94FBwLO4t8/dl8vEBSVSTBMqCZrJlpZWvZUCquYtUBhjUpFKoDhFVRdHFlT1edzQ36/ihhwvOKF8gCZoJptIKK/VGFOnUgkUW0Tk1yJytDeNAz4VkfrA/hynz6SjXz83iuytt7rPGEEiOrdkldvGmGRSCRSXAq2Bp7ypjbeuPnBx7pJm0pZCM1nLQRhj0pU0UKjq56r6C1Xt5k2jVXWrqu5R1fV1kUiTghSbyUKImwIbYwKRyjDjLYFxwAlASWS9qp6Vw3SZdEWGHI+8Sc//lqvIdq+5rH8kWWsua4xJJpWipz8B7wLtgMnARtxrTAteqH5hR5rJ+ls/+Ze9QQKTCdU1G2PqRCod7parag8ReUtVu3jrylU1tSdPHcvm+yhC+2s7yetR/aLfUxHaazbGZCTTDneRHthbRORHItINiNVbO50EHSYiL4jIe97n92Ls01VEXhGRVSLyloj8OJNzFpUUWj9FWA7CGJNMKoHiNhE5FLgJGAs8ANyQ4XlvBpaoantgibccbSdwpaqeAAwEfi8izTM8b1IFUdGbxiCBUCDXbIzJmUCGGReRtUBfVd0iIkcAy1S1Q5LvrAQuVNX3Eu1X9EVP/tZP/frVXE4ilNdsjMlYoqKnVFo9tQN+AbT176+q52WQpsNVdYs3/wlweJI0nAw0At6Ps30kMBKgTZs2GSSrAKTR+skYY1KRNFDgOtk9CDxDGj2xReRFYr9r+xb/gqqqiMT9DevlOOYCI1Q15vlVdRYwC1yOItU0JhPKMZ/8QSA6JxFZTiDWNUdXeBtjiksqrZ5eU9X/zOpJUyx6EpFmwDLg/6rq/FSOnc2ip9BLo/VTIlYcZUzhy7TV0x9EZJKInCoi3SNThmlaCIzw5kcAT0fvICKNgAXAw6kGCRMljdZPxhgTTyqBojPw38DtwHRvujPD894ODBCR94D+3jIi0lNEHvD2uRg4A7hKRFZ4U9cMz1tc0mz95GctoYwxEakUPa0HOqnqnrpJUmas6Mnjb+1UXg4NGsDUqdVbQ6VYsW1FT8YUvoxaPQHvAM2Bz7KaKpNbkdZP/ors8eOrWj+lULFtjDGQWtFTc+BdEVksIgsjU64Tlo9CVewSGfsJqprITp0KlZVp9auA+C2hjDHFIZWipzNjrVfVl3KSogzlsugp9EUwEye6iu0JE9xQ5BkI/b0wxlSTUdFTvgYEk6ZIxfbZZ8Pdd7vcRCRHkUZ9hTGm+MQtehKR7SLyTYxpu4h8U5eJDFJBtP7xV2zfcovLCpx/vlsf2ZbCMOQFcS+MMWkLZKynXLKipximTXOBwJ+DGDbMrVu5slYd8UJ7L4wxMWXa4c6Enb9iG9z89dfDkiVZ7YhnOQtjCpMFijSEcuynWDLoiBcR615MnpyFtBlj8o4FijQUxC9mf33FIYe4vhWRV6dGtk+blvQwBXEvjDEpsUBRbPwd8Xr1cn0rIh3x0qjYjrAKbmMKn1VmF7ssjTALVsFtTJhZZbaJz0aYNcYkYYGi2EV3xPNXbKdYXxFRMJX9xphqLFBkQWjL47PUES8i3n0I7f0xxgBWR5EVoS2bz0FHvFhCe3+MKSJWR2Fiq6OOeMaYcLNAUUsF2Sw0i/UVBXl/jClSVvSUBQVRtOKvrwBXBCUCCxa45TTfYeFXEPfHmAKXd0VPInKYiLwgIu95n99LsG8zEdkkIvfUZRqLjr8jXr9+8NRT7un+299mFCSMMeEXVNHTzcASVW0PLPGW47kV+HudpKqWCqJZaA7rK+wNecaEW1CBYigwx5ufAwyLtZOI9AAOB56vo3TVSkE+9PwDB86Y4abo7SnWWcS6PzaAoDHhkfQNdzlyuKpu8eY/wQWDakSkHjAduBzon+hgIjISGAnQpk2b7Ka0GPnrK/r1g+bNYexYt+3GG2vWZxhjClrOchQi8qKIvBNjGurfT11teqyqzuuARaq6Kdm5VHWWqvZU1Z4tW7bM0hUUMX99BbjgcOedLncxcWKt6yysJZQx4RRIqycRWQv0VdUtInIEsExVO0Tt8yegD7AfOARoBNynqonqM/JuUMDS0gJ6EE6c6MaEmjABpkzJ6FDWEsqY/JJ3rZ6AhcAIb34E8HT0Dqp6maq2UdW2wFjg4WRBIh8VTFl8FvtYGGPCJahAcTswQETew9U/3A4gIj1F5IGA0mTiyfKYUGAtoYwJE+twlwOlpbFzEpMmhfRhaGNCGVPwEhU9WaDIsYJ9+EXqK664Ah5+uGp9WZmrDB83Lu1DFuy9MiYE8rGOwoRZpL7iiivgkUeq+ljYq1SNKUhB9aMoGgXRa9svuo9F166uj8WKFfDcc2kXQ/lbhVmOwpj8ZDmKHCu4X8ax+lhcfjnMnQsnnVRzf2sRZUzoWaAw6YkeE6qszOUkJkyAigpXyR1pOptmUVR07qvggqwxIWWV2ab2oouhyspcs1lVN6DgzJkZtYiyoihj6o5VZuexUP9qji6G6tfPvb+iVy/XImrQoJq5DyuGMiZ0LFAELNQ9t6OLoSJWrqx1iyhrBWVM/rFWTyZ7stAiylpBGZN/LEcRgIL91ZyoRZQVQxkTWhYoAlBa6n4pR34tR+ZDHyjitYiqZTGUvxVU6O+NMSFmrZ4CVrDFK9HFUDNmuGKorl1hwwZX6e0fOyrJsB8Fe5+MyRPW6imPFVzP7Yh4xVBvvgl79lTtV8vRZ40xdccCRcAKtkglUce8Ro1cx7wkb8sr2LocY0LGAkUIhP7B6C+GmjLFFTt9+63razFqVPUiKF8Fd8HW5RgTMhYoQiDUfS2gZjEUuFxFgwZVb8uzIihj8pYFCpN7/mKoSEB45hl4/nnYvds1nR02rPpQIFFNZyN1OZabMKbuBRIoROQwEXlBRN7zPr8XZ782IvK8iKwRkdUi0rZuUxqcgi2f9+cu+vWDiy5yweK446qCRIycReS6Q5+7MiaEAmkeKyLTgC9V9XYRuRn4nqr+OsZ+y4DfquoLInIIsF9VdyY6dtiax6aiYJuGRoLCoEGun0UKTWcL9l4YE7B8bB47FJjjzc8BhkXvICKdgAaq+gKAqlYmCxImRPwV3A8/XNV0trLSffr2+duSBoyTaYWXuzImJIIKFIer6hZv/hPg8Bj7HA98LSJPisibInKHiNSPdTARGSkiFSJSsXXr1lylOTAF2dfCXwTl78Fdrx7cdBP06OGGLB8/noFvTGXa0l7o0jJ+xTRr/WRMHctZoBCRF0XknRjTUP9+6sq+YhUmNAD6AGOBXsAxwFWxzqWqs1S1p6r2bNmyZXYvJA+k8kAM3UMzUsEdnbNYvBgaN4Y33oDt211fi3nz3HcuvpjvaMALA2yMKGPqUs4Char2V9UTY0xPA5+KyBEA3udnMQ6xCVihqhtUdS/wFNA9V+kNu9BW8sZqOtu4MbRtC/v3w44d8KtfuWAyfjy3HTyVdS9uhMGDg0qxMUUnqKKnhcAIb34E8HSMfcqB5iISySKcBayug7SZuhSr6eykSa6u4oor3Prly6FpU5g6lSYjLmIU90P//sGl2ZgiE1SguB0YICLvAf29ZUSkp4g8AKCq+3DFTktE5G1AgP8NKL15qeCa0JaXw/jxMHWqy2VcfTUcfLC7qA8+4JPPhf0z7+cm7kRuujHc12pMiNjosQWiYJqNTptW1YfCK25i8mTYt88VQ7Vrh3ywoTCu1Zg8ko/NY42JLVIU5c9dXHYZ7NwJ7dqhH3zAPVwXdCqNKSoWKApEqk1oQ1NUM24c7N3rem7ffz/ceSds2MB9jOI6Zla9BMkYk3NW9FRkQldENXiwq7i+8UbAS//0GfDii7BoUcCJM6ZwWNGTCa9Fiyj95sbqlfY33Yg8tyg8uSNjQs4CRREIe+uo0lJXtGbvpTAmGFb0VGRSKXoqLc2/h3Ak3aErOjMmJKzoyaQln3t5F+S4V8bkuQZBJ8DUrTA9aPv2hZdeqlqOFJ1B/uV4jClklqMoMvEesOnUY9TVQ/qll+yd2cbkA6ujMDUkqwfIRT1BrHoR/3msbsKY3LI6ClOnkv3ij7U9Ui8SL2dz5plZTKAxJi0WKEwNseox0imaSlYZnmh7aWns4qZly5Kn2xiTGxYoTA3x6iWyXV8Q9v4dxhQLCxQmK5I99GNtnzw5cUe6MLXQMqaQWWW2SVuyDnm1qQy3ympjgmWV2SarclE0ZLkHY/KXBQqTdcke+vEqy40x+cmKnowxxuRf0ZOIHCYiL4jIe97n9+LsN01EVonIGhG5W8Q/iIMxxpi6EFTR083AElVtDyzxlqsRkdOA3kAX4ESgF2Ddrowxpo4FFSiGAnO8+TnAsBj7KFACNAIaAw2BT+skdcYYYw4IKlAcrqpbvPlPgMOjd1DVV4AyYIs3LVbVNbEOJiIjRaRCRCq2bt2aqzQbY0xRytkw4yLyIvAfMTbd4l9QVRWRGjXqInIc0BFo7a16QUT6qOo/ovdV1VnALHCV2Zmm3RhjTJWcBQpV7R9vm4h8KiJHqOoWETkC+CzGbucDr6pqpfed54BTgRqBwm/58uWfi8iHGSS9BfB5Bt8PWtjTD3YN+SDs6Qe7hnQdHW9DUC8uWgiMAG73Pp+Osc+/gf8WkamA4Cqyf5/swKraMpOEiUhFvCZiYRD29INdQz4Ie/rBriGbgqqjuB0YICLvAf29ZUSkp4g84O0zH3gfeBtYCaxU1WeCSKwxxhSzQHIUqvoFcHaM9RXANd78PuBndZw0Y4wxUWwIj5pmBZ2ADIU9/WDXkA/Cnn6wa8iaghvCwxhjTHZZjsIYY0xCFiiMMcYkZIHCIyIDRWStiKwXkRpjT+U7EfmjiHwmIu8EnZbaEpGjRKRMRFZ7g0FeH3Sa0iEiJSLyuois9NKf5O3h+UtE6ovImyLybNBpqQ0R2Sgib4vIChEJ3XDSItJcROaLyLveoKinBpoeq6NwfxTAOmAAsAkoBy5R1dWBJiwNInIGUAk8rKonBp2e2vA6Xx6hqm+ISFNgOTAsLP8O3ujGTVS1UkQaAv8ErlfVVwNOWtpE5EagJ9BMVYcEnZ50ichGoKeqhrLDnYjMAf6hqg+ISCPgYFX9Oqj0WI7CORlYr6obVHUP8Bhu4MLQUNW/A18GnY5MqOoWVX3Dm98OrAGODDZVqVOn0lts6E2h+yUmIq2BHwEPJNvXZJ+IHAqcATwIoKp7ggwSYIEi4kjgI9/yJkL0gCpEItIW6Aa8FmxK0uMV2azADUvzgqqGKv2e3wPjgP1BJyQDCjwvIstFZGTQiUlTO2Ar8JBX/PeAiDQJMkEWKEzeEZFDgL8Av1TVb4JOTzpUdZ+qdsUNZnmyiISqGFBEhgCfqeryoNOSodNVtTswCPi5VzQbFg2A7sBMVe0G7CDGO3vqkgUK52PgKN9ya2+dqWNe2f5fgD+p6pNBp6e2vKKCMmBg0GlJU2/gPK+M/zHgLBF5JNgkpU9VP/Y+PwMW4IqXw2ITsMmXG52PCxyBsUDhlAPtRaSdV3E0HDdwoalDXmXwg8AaVZ0RdHrSJSItRaS5N38QrnHEu8GmKj2qOl5VW6tqW9zfwVJVvTzgZKVFRJp4jSHwimzOAULTGlBVPwE+EpEO3qqzgUAbdAQ1emxeUdW9IjIaWAzUB/6oqqsCTlZaROTPQF+ghYhsAiap6oPBpiptvYErgLe9cn6A/1HVRQGmKR1HAHO8VnT1gHmqGsrmpSF3OLDA/e6gAfCoqv4t2CSl7RfAn7wfrhuAq4NMjDWPNcYYk5AVPRljjEnIAoUxxpiELFAYY4xJyAKFMcaYhCxQGGOMScgChTExiMg+b+TRd0TkmUj/iFoea6OItMhm+rzj9hWR07J9XGOiWaAwJrZvVbWrNxLvl8DPg06Qn4g0wPWbsUBhcs4ChTHJvYJvkEgR+ZWIlIvIW/53TojIU94gdKtSGYhORCpF5C5v/yUi0tJb/9/e8VeKyF9E5GBv/WwRuV9EXgPmAdcCN3g5nz7ZvmhjIixQGJOA18v6bLwhXUTkHKA9buygrkAP34BzP1HVHrj3OIwRke8nOXwToEJVTwBeAiZ5659U1V6qehJuqPWf+r7TGjhNVS8A7gfu8nI+/8j0Wo2JxwKFMbEd5A0j8gluSIgXvPXneNObwBvAD3GBA1xwWAm8ihtksj2J7Qce9+YfAU735k8UkX+IyNvAZcAJvu88oar7an1VxtSCBQpjYvvWGy78aECoqqMQYKr3K76rqh6nqg+KSF+gP3CqlxN4EyhJ85yR8XRmA6NVtTMwOeo4O2p1NcZkwAKFMQmo6k5gDHCTV4G8GPiJ984MRORIEfkBcCjwlaruFJEfAqekcPh6wIXe/KW4V6cCNAW2eEOuX5bg+9u9fY3JKQsUxiShqm8Cb+Heo/488Cjwilc0NB/3sP4b0EBE1gC344qfktmBe7nRO8BZwBRv/QTcm/1eJvEw5c8A51tltsk1Gz3WmICISKWqHhJ0OoxJxnIUxhhjErIchTHGmIQsR2GMMSYhCxTGGGMSskBhjDEmIQsUxhhjErJAYYwxJqH/D9OKy6eTpTGqAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "title = 'Eigenvalues distribution for $P_e = {Pe}$'\n", + "\n", + "plt.plot(eigh_r, eigh_i, \"+b\", label=\"glt symbol\")\n", + "plt.plot(eig_r, eig_i, \"xr\", label=\"eigenvalues\")\n", + "plt.legend(loc=2);\n", + "\n", + "plt.title(title.format(Pe=Pe))\n", + "plt.xlabel('Real part')\n", + "plt.ylabel('Imaginary part');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.8" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/notebooks/psydac_example_2.ipynb b/notebooks/psydac_example_2.ipynb new file mode 100644 index 0000000..9bb7d4b --- /dev/null +++ b/notebooks/psydac_example_2.ipynb @@ -0,0 +1,167 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 2D Poisson" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "from sympy.core.containers import Tuple\n", + "from scipy.linalg import eig as eig_solver\n", + "import matplotlib.pyplot as plt\n", + "from IPython.display import Math\n", + "\n", + "from sympde.core import Constant\n", + "from sympde.calculus import grad, dot\n", + "from sympde.topology import ScalarFunctionSpace\n", + "from sympde.topology import Square\n", + "from sympde.topology import element_of\n", + "from sympde.expr import BilinearForm\n", + "from sympde.expr import integral\n", + "from gelato import GltExpr\n", + "from psydac.api.discretization import discretize" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "# abstract model\n", + "domain = Square()\n", + "\n", + "V = ScalarFunctionSpace('V', domain)\n", + "\n", + "u = element_of(V, name='u')\n", + "v = element_of(V, name='v')\n", + "\n", + "k_para = Constant('k_para', real=True)\n", + "k_perp = Constant('k_perp', real=True)\n", + "\n", + "bx = Constant('b_x', real=True)\n", + "by = Constant('b_y', real=True)\n", + "b = Tuple(bx, by)\n", + "\n", + "expr = k_para * dot(b, grad(v)) * dot(b, grad(u)) + k_perp * dot(grad(v), grad(u))\n", + "a = BilinearForm((v,u), integral(domain, expr))\n", + "\n", + "glt_a = GltExpr(a)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "ncells = [16,16]\n", + "degrees = [3,3]" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "# create the computational domain from a topological domain\n", + "domain_h = discretize(domain, ncells=ncells)\n", + "\n", + "# discrete spaces\n", + "Vh = discretize(V, domain_h, degree=degrees)\n", + "\n", + "# dsicretize the equation using Dirichlet bc\n", + "ah = discretize(a, domain_h, [Vh, Vh])\n", + "\n", + "# dsicretize the glt symbol\n", + "glt_ah = discretize(glt_a, domain_h, [Vh, Vh])" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "# ...\n", + "eigh = glt_ah.eig(k_perp=1., k_para=1.e4, b_x=1., b_y=1.)\n", + "eigh = eigh.ravel()\n", + "eigh.sort()\n", + "# ..." + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "# use eigenvalue solver\n", + "M = ah.assemble(k_perp=1., k_para=1.e4, b_x=1., b_y=1.).tosparse().todense()\n", + "w, v = eig_solver(M)\n", + "eig = w.real\n", + "eig.sort()" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de3hU1bn48e87k5Ag4hWkXNQgIgcUL1y8FBVi1UKKgKeCKLVobUGsR89Ra6UeSQAtGCtWTxWx/XnBYhFsUbBYRQ1q7UECR7xSFAQroIKKXITcZt7fH2vvzM4kE5IwyUyS9/M888zM2ntmVkbc76z1rouoKsYYY1q3UKorYIwxJvUsGBhjjLFgYIwxxoKBMcYYLBgYY4wBMlJdgYbq0KGD5uTkpLoaxhjTrKxevfpLVe0YX95sg0FOTg6rVq1KdTWMMaZZEZFPaiq3biJjjDEWDIwxxlgwMMYYQzPOGdSkvLyczZs3U1JSkuqqtDrZ2dl069aNzMzMVFfFGNMA+w0GIpINvAZkeec/rar5ItIdmA8cCawGrlDVMhHJAuYC/YGvgEtVdZP3XpOBq4EIcL2qvuCVDwXuA8LAH1R1ZkP+mM2bN9O+fXtycnIQkYa8hWkAVeWrr75i8+bNdO/ePdXVMcY0QF26iUqB81T1FOBUYKiInAncBdyrqscDO3AXebz7HV75vd55iEgfYCxwIjAUeFBEwiISBh4AhgF9gMu8c+utpKSEI4880gJBExMRjjzySGuRGZNshYVwxhkwfDgUFbmyoiK46CJXXliYtI/abzBQZ4/3NNO7KXAe8LRX/jgwyns80nuOd/x74q7OI4H5qlqqqhuB9cDp3m29qn6sqmW41sbIhv5BFghSw753YxrBwIHw3nvw17+6ADBrlgsMzz0H77/vjidJnRLI3i/4NcA2YBmwAfhGVSu8UzYDXb3HXYFPAbzjO3FdSZXlca9JVF5TPSaIyCoRWbV9+/a6VN0YY5qv3Fx34T/oIPj2W7jpJti7F9q1gyVL3PEkqVMwUNWIqp4KdMP9kv+3pNWgHlT1YVUdoKoDOnasNoEurV155ZU8/bRrSP32t79l7969TfK5mzZt4qSTTqrXa4YMGWIT+oxJF7m5LggE3XhjUgMB1HNoqap+AxQBZwGHiYifgO4GbPEebwGOBvCOH4pLJFeWx70mUXmTKShoyk9r2mBgjGnmiorgnnuqls2aFcshJMl+g4GIdBSRw7zHbYELgLW4oHCJd9p44Fnv8WLvOd7xV9Rtp7YYGCsiWd5IpJ7ASqAY6Cki3UWkDS7JvDgZf1xdTZ2avPeaPn06vXr14uyzz+ayyy7jN7/5TZXj999/P1u3biU3N5fcGiL7rbfeSp8+fTj55JO5+eab2b17N927d6e8vByAXbt2VT4fMmQI//Vf/8WAAQPo3bs3xcXF/Pu//zs9e/bkv//7vyvfs6KignHjxtG7d28uueSSykD08ssvc9ppp9G3b19+8pOfUFpamrwvwhhTf4WFMHFi7DZrFuTlua6hjAw480zIzHRdRhddlNyAoKq13oCTgbeAd4D3gCle+XG4i/l6YCGQ5ZVne8/Xe8ePC7zXbbh8wzpgWKA8D/jQO3bb/uqkqvTv31/jffDBB9XK6qIyTX6AVq5cqaeccoru27dPd+3apccff7zefffdqqo6fvx4XbhwoaqqHnvssbp9+/Zqr//yyy/1hBNO0Gg0qqqqO3bsUFXVK6+8UhctWqSqqnPmzNEbb7xRVVUHDx6st9xyi6qq/va3v9XOnTvr1q1btaSkRLt27apffvmlbty4UQH9+9//rqqqV111ld599926b98+7datm65bt05VVa+44gq99957K9+3uLi43n9/Q79/Y4znlVdUDzlEtV07d8vMVBVRDYVUs7NVDz1U9Z57VIcPVz39dNW77qr3RwCrtIZral1GE72jqqep6smqepKqTvPKP1bV01X1eFUdraqlXnmJ9/x47/jHgfe6U1V7qGovVX0+UL5UVU/wjt3ZsLBWPwUFIOJuEHt8IF1Gb7zxBiNHjiQ7O5v27dtz0UUX1ev1hx56KNnZ2Vx99dX85S9/4aCDDgLgpz/9KY8++igAjz76KFdddVXla0aMGAFA3759OfHEE+ncuTNZWVkcd9xxfPqpy8sfffTRDBo0CIAf/ehH/P3vf2fdunV0796dE044AYDx48fz2muvNfyPN8YcuNxceOYZCIehogLKy0HVtQaysmDRIpcvWLIE3nwTbrklaR/dapejKChw37FrmMQeN3X+ICgjI4OVK1dyySWX8NxzzzF06FAABg0axKZNm1i+fDmRSKRKQjgrKwuAUChU+dh/XlHhBnvFD/u0YaDGpLHcXLjhBgh225aWwvXXJz1pHNRqg0FjGDRoEEuWLKGkpIQ9e/bw3HPP1Xhe+/bt2b17d7XyPXv2sHPnTvLy8rj33nt5++23K4/9+Mc/5vLLL6/SKqirf/3rX/zv//4vAE8++SRnn302vXr1YtOmTaxfvx6AJ554gsGDB9f7vY0xSVZUBPfd51oCvqwsuP/+pCeNgywYAPn5yXmfgQMHMmLECE4++WSGDRtG3759OfTQQ6udN2HCBIYOHVotgbx7926GDx/OySefzNlnn82sWbMqj40bN44dO3Zw2WWX1btevXr14oEHHqB3797s2LGDSZMmkZ2dzaOPPsro0aPp27cvoVCIa665pv5/tDEmeYqKYNQoiERcwrhdOzfHICPDdRtdfHGjBQRRv5+kmRkwYIDGj4Vfu3YtvXv3TlGNnD179nDwwQezd+9ezj33XB5++GH69et3wO/79NNP8+yzz/LEE08koZaNIx2+f2OatcJC2LAh9nzsWHc/f36srEePA8oViMhqVR0QX96iVi1NBxMmTOCDDz6gpKSE8ePHJyUQ/Md//AfPP/88S5cuTUINjTFpyQ8EY8fGcgNFRS4QHGAAqAsLBkn25JNPJv09/+d//ifp72mMSRN+EOjVy134n3oKpkxxgeCll2KjiBqZBQNjjEmlgQPhzjvd2Pb8fBcI/OUn2rVzgaARRxH5LIFsjDGp5M8tUIXbb686pLQR1iBKxIKBMcakmj+3YO9eN2ooIwPatm304aRBFgyMMSbV4hejy8qCO+5wrYVGHE4aZMGgCeTl5fHNN9+ktA4NWcraGNME/LkF0Sj84AcuKITDMG2ayyFceikUFzd6NVpvArmw0CVugv1xRUXuS0/yEC4bEmqMSai42A0nDQ4pPe00N7KoogLmzGmSarTelsHAgTBmTNV9RceMOeBt5P74xz9y+umnc+qppzJx4kQikQg5OTl8+eWXQOIlrjds2MDQoUPp378/55xzDv/85z8BtynO9ddfz3e/+12OO+64yg1yxo4dy1//+tfKz/U3z9m0aRPnnHMO/fr1o1+/fvzjH/+oVsfHHnuM6667rvL58OHDWb58OQAvvvgiZ511Fv369WP06NHs2eN2PI1fWtsYkyS33OIu+MEfprm5rqyR5xZUUdNSps3hlpQlrF95RbVDB9Xbb3f3r7xSv9fX8PnDhw/XsrIyVVWdNGmSPv7445VLVte2xPV5552nH374oaqqrlixQnNzc1XVLX19ySWXaCQS0ffff1979Oihqqp/+ctf9Mc//rGqqpaWlmq3bt107969+u233+q+fftUVfXDDz9U/3vauHGjnnjiiaqq+uijj+rPf/7zynr/4Ac/0KKiIt2+fbuec845umfPHlVVnTlzpk6dOjXh0to1/f3GmHq4667q151XXmnQ0tR1RYIlrFtvNxG46DtpEkyf7oZ0HeAQrpdffpnVq1cz0Gtd7Nu3j6OOOqryeHCJ6+zs7Molrvfs2cM//vEPRo8eXXlucKOZUaNGEQqF6NOnD1988QUAw4YN44YbbqC0tJS//e1vnHvuubRt25adO3dy3XXXsWbNGsLhMB9++GGd679ixQo++OCDyuWuy8rKOOuss6osrT18+HCGDx/e8C/JGOMUFrpRQ2PGwIIF7voza5abZ7BkSZNXp3UHg6IimD3bBYLZs91/jAMICKrK+PHjmTFjRpXyxx57rNbXRaNRDjvsMNasWVPj8eDS1OqtJZWdnc2QIUN44YUXeOqppxjrrWFy77330qlTJ95++22i0SjZ2dnV3i8jI4NoNFr5vKSkpPK9L7jgAv70pz9Ve83KlSt5+eWXefrpp/nd737HK6+8UuvfZIypRV4e5OTAwoUwebILCMceC6tXuwRyE80tCGq9OQM/R7BggcvaL1hQNYfQAN/73vd4+umn2bZtGwBff/01n3zySeXxREtcH3LIIXTv3p2FCxcC7qIcXL46kUsvvZRHH32U119/vXLvg507d9K5c2dCoRBPPPEEkUik2utycnJYs2YN0WiUTz/9lJUrVwJw5pln8sYbb1Qua/3tt9/y4Ycf1rq0tjGmAc4/Hx56CEaPhhkzoH17Fwj693cTzVKg9bYMiotjTTNw9wsWuPIGRuU+ffpwxx13cOGFFxKNRsnMzOSBBx6oPB5c4rpTp05VlrieN28ekyZN4o477qC8vJyxY8dyyimn1Pp5F154IVdccQUjR46kTZs2AFx77bX88Ic/ZO7cuQwdOpR27dpVe92gQYPo3r07ffr0oXfv3pWL6XXs2JHHHnuMyy67rLKb6o477qB9+/aMHDmSkpISVLXK0trGmAbwL/g33wwdOsDGjdC9O3zyiftBmoKWgS1h3cQaa4nrdNAcvn9j0spxx7lA0LGjm2A2ebJrKQR/qCZZoiWsW283UYpMmDCBU089lX79+vHDH/6wxQQCY0w9XXttrEXw5ZexLqPJk5tkklm81ttNlCKNscS1MaaZmTXLDVqZNAkefNA9v/lmuOYaN9GsKecXeFpcMFBV2/A9BZprd6MxKfHSS27UkJ878O9feskFhxRoUTmDjRs30r59e4488kgLCE1IVfnqq6/YvXs33bt3T3V1jDG1aBXbXnbr1o3Nmzezffv2VFel1cnOzqZbt26proYxpoH2GwxE5GhgLtAJUOBhVb1PRAqAnwH+lfdXqrrUe81k4GogAlyvqi945UOB+4Aw8AdVnemVdwfmA0cCq4ErVLWsvn9MZmam/TI1xpgGqMtoogrgJlXtA5wJ/FxE+njH7lXVU72bHwj6AGOBE4GhwIMiEhaRMPAAMAzoA1wWeJ+7vPc6HtiBCyTGGNOyFBbCxIlVJ7cWFbmywsLU1Ys6tAxU9TPgM+/xbhFZC3St5SUjgfmqWgpsFJH1wOnesfWq+jGAiMwHRnrvdx5wuXfO40ABMLv+f44xxqQxf7/jp56KbXI/apTb/7gJNr2vTb3mGYhIDnAa8KZXdJ2IvCMij4jI4V5ZV+DTwMs2e2WJyo8EvlHVirjymj5/goisEpFVlhcwxjQr/i9/f7/j4cNh2DAoK2uyTe9rU+dgICIHA38G/lNVd+F+ufcATsW1HO6p5eVJoaoPq+oAVR3QsWPHxv44Y4xJjuAKpRDb77i01E02S3EggDqOJhKRTFwgmKeqfwFQ1S8Cx38PPOc93QIcHXh5N6+MBOVfAYeJSIbXOgieb4wxzZ+/mdbkya5baN8+V56ZCYsXp2w9oqD9tgzEDdj/f8BaVZ0VKO8cOO1i4D3v8WJgrIhkeaOEegIrgWKgp4h0F5E2uCTzYm+zhSLgEu/144FnD+zPMsaYNOIvhDl1KuzZA+XlcMEF8MILTbrpfW3q0k00CLgCOE9E1ni3PKBQRN4VkXeAXOC/AFT1fWAB8AHwN+DnqhrxfvVfB7wArAUWeOcC/BK40Us2H4kLPsYY03Lk5sLxx7uN7/v2hbfecuXPPNNkm97XpkXNQDbGmLTlrz/0ox/B8883yQqlNbFVS40xJlWKityOir/5Dcyd6wJAClcorYkFA2OMaWzFxfDcc7EF6fwcQopWKK2JdRMZY0wrYt1ExhhjErJgYIwxjSkvzyWPg2bNcuVpxIKBMcY0pvPPd6OI/IDgjyo6//zU1itOi9rPwBhj0o6fNL75Zjen4O9/d6OK/PI0YS0DY4xpLIWFbljpjTfC2WfD66/DSSe5UURpxloGxhjTGIKL040e7VoE3bvDu+/Cpk2prl011jIwxpjGMHCgm1j23e/C7NnQr58LAiNGwEMPVU8qp5gFA2OMaQz+xLLnn3ctgtWr3VIUzz7rcgYvvZTqGlZhwcAYYxpLbi6MHQsbN8I557jA4OcQli5Nde2qsGBgjDGNZdYs+OMf4YorYO1atxbRmDEpX666JpZANsaYZPOTx/7idDfe6ALDlCkwbZpbqygNdjcLsmBgjDHJNnCg2+N4+nQXCIqKXDJ52rS0WpwuyIKBMcYkW26uW6V0zBj45hs3mqiJ9y2oL8sZGGNMY8jNhUmTXOtg0qQDDgQFBTBkiLtvDLaEtTHGJFNhoesmAtcymDQJ7r/fbW05Z06D31Yk9vhALtuJlrC2biJjjEmmgQNh1Ch39V60yJXddx889ZQbZlrPFkJjtQTiWTeRMcYkS2Ghux871v18LyqCiy92zxctqtcWl3630NSp7hYk4m7JDBTWMjDGmGTIy4OcHLj7bpcs7tTJ5QtCoViLoB6tgvgAENQYvfsWDIwxJhn8fQuuucZ1E+3b58rbtKn3WzVV11CQBQNjjEkGf3+Cm25yrYFoFC64IDbruA5DSwsKErcIBg9290OGJK3GVdhoImOMSabOneHzz+Hoo13rYMECV15cXK/JZsHRQ5C8rqFEo4n2m0AWkaNFpEhEPhCR90XkBq/8CBFZJiIfefeHe+UiIveLyHoReUdE+gXea7x3/kciMj5Q3l9E3vVec79I/NdgjDFpzN/E5tprXSDo2xc2b4bjj3etAmjwrOP8fHdrbHXpJqoAblLV/xOR9sBqEVkGXAm8rKozReRW4Fbgl8AwoKd3OwOYDZwhIkcA+cAAQL33WayqO7xzfga8CSwFhgLPJ+/PNMaYRuQvP7F3r5tXMHq0e/7mmy6H0IC1iPwA0FT5g/0GA1X9DPjMe7xbRNYCXYGRwBDvtMeB5bhgMBKYq67/aYWIHCYinb1zl6nq1wBeQBkqIsuBQ1R1hVc+FxiFBQNjTHORm+u2s3zvPejQwbUGnnsO3nrL7Vvw4IP1fsumTiLXa56BiOQAp+F+wXfyAgXA50An73FX4NPAyzZ7ZbWVb66hvKbPnyAiq0Rk1fbt2+tTdWOMaRx+F9Gbb7rk8fTpMGyYaw2k4b4FidQ5GIjIwcCfgf9U1V3BY14roNEz0ar6sKoOUNUBHTt2bOyPM8aYxAoLYeLE2D7Hs2a5ZSd69oQnnnDlzUidaisimbhAME9V/+IVfyEinVX1M68baJtXvgU4OvDybl7ZFmLdSn75cq+8Ww3nG2NMevL3K5g/3w37ufzy2JDSXbtc3mDGDDjttHrlCuK7hpq0q0hVa70BAswFfhtXfjdwq/f4VqDQe/wDXH+/AGcCK73yI4CNwOHebSNwhHdspXeueK/N21+9+vfvr8YYkxKvvKLaoYPqPfeoHnKIajis6kZ/ql5xReycu+6q19v6b+HfGgOwSmu4ptalZTAIuAJ4V0TWeGW/AmYCC0TkauATwBs/xVIgD1gP7AWu8oLO1yIyHfAX55imXjIZuBZ4DGjrBQNLHhtj0ldxsZtMNmMGnHIKvP66K8/Jie1zXMflJ4YMabyJZPVhk86MMaa+iopcnuC002DZslh5u3ZuN7MZM+q8mU1dZlXl5yevy6jBk86MMcbEyc11cwmWLYtdzSdNgnDYBYPJk+u1QmkifodRU+QOmle62xhjUs1PHj/+uBs59NFHbg2iuXNdIFi3br/7HA8ZAq++2nRVrgsLBsYYU1d+ILj9dhg/HhYudIFg2TLXMqioqNNuZsuXxx4Hu4maYtmJRCwYGGNMXQ0c6HIF48fDQw+5ZatfeskFgoULY4vSNVAqlq72WTAwxpj9KSyEDRvcJjULFrh1h3JyXIugXz+33MTo0Q1ag2jw4PQYTWTBwBhj9mfgQLjzTreP8aJFMGiQCwThsAsScUNJE/3CX768+oU/2GWUSja01Bhj6qKoKLaDWXk5ZGbCQQfBlCnVhpLWZxH+pr4E29BSY4xpCH8hutxcGDnSBQKAc891rYQZM5I2lDSVLBgYY0xt/KTxrFnu139mpitfscLdL1gAFRUU7L0Fkfq1CoDK16QyeQzWTWSMMYkVFrpg8NZbbiG6rCw3tDQ3F157zV3FFy2qljS2biJjjGkpCgvdmkMXX+wmkvXrB6WlUFICXbrAM8/ApZc2++4hn40mMsaYIH8Yaa9e7td/JOJmF/u5gkjEHatlIbpEk8dqGk2ULiwYGGNMkD+MVMRd1W+7zbUGwOULZs7c714Fqe7/bwjrJjLGGF9hobt/5hnXmX/77a5ryJeZ6YLAggUtpnvIZ8HAGGN8/sghgBtugL17XVAIhaBtW9dauPhid7yWheiaIwsGxhgDsVbBggVuctnMmbFjbdvCHXe4GccVFW67yxbGgoExxoBLGl98sRtGWlrqEsYicOaZsX0K8vNh3Djo0SPVtU06SyAbY0xhoRshpAq33grRqAsAkYhbgO6001xroI5LVDdH1jIwxpiBA90IoXHjXIsgEnG3SZNcObgg0MLyBEEWDIwxrVswV/DIIy5ZDG6m8ejRLXLkUE0sGBhjWjd/BNFbb7nn0agLBFlZLXbkUE0sGBhjWqfCQpg40T1esMDlCkpLXcsgI8MljFVb5MihmlgwMMa0Pv5exvPnx0YQRSLuWGamm4E8Y4YbPdQCRw7VxIKBMab18RPG+fnu1/8vfuG6hzIzITs7Nsu4oqLeXUQFBc1zOYr9LmEtIo8Aw4FtqnqSV1YA/AzY7p32K1Vd6h2bDFwNRIDrVfUFr3wocB8QBv6gqjO98u7AfOBIYDVwhaqW7a/itoS1MaZB/GWpweUKOneGd991z2+/3a03NGZMlZ3LahN/4Z861d2n6+4AiZawrss8g8eA3wFz48rvVdXfxH1IH2AscCLQBXhJRE7wDj8AXABsBopFZLGqfgDc5b3XfBF5CBdIZtf5LzPGmPrYsAF+/Wu3D8GwYfDEE648HIb773cBwB9BFAgGQ4bUvOKof/Fv7vYbDFT1NRHJqeP7jQTmq2opsFFE1gOne8fWq+rHACIyHxgpImuB84DLvXMeBwqwYGCMaSxjx7pcQV5ebDVSgAkTYN48l0NYtKha99Crr7pbXfkb3OTnN49uowOZgXydiPwYWAXcpKo7gK7AisA5m70ygE/jys/AdQ19o6oVNZxfjYhMACYAHHPMMQdQdWNMq+LvUQAuGOTnu53LwI0eysuDhQtd+bp1UFxMwav77yKqTXMJAr6GBoPZwHRAvft7gJ8kq1KJqOrDwMPgcgaN/XnGmBYgOHIoEoGnnnI7lfmiUdcddOONrmvIW25iaj33Mo7XnAIBNHA0kap+oaoRVY0CvyfWFbQFODpwajevLFH5V8BhIpIRV26MMQemsBDOOAM2bYqNHALYuRPWrnWPs7KgXTs3pwCSMrksPz/xTmfprEEtAxHprKqfeU8vBt7zHi8GnhSRWbgEck9gJSBAT2/k0BZckvlyVVURKQIuwY0oGg8829A/xhhjKg0c6LK7xcVwzTXugr93b9VzbrnFtQpGjYL58yl4NbdBCeHgxb+5tQh8dRla+idgCNAB+ALI956fiusm2gRM9IODiNyG6zKqAP5TVZ/3yvOA3+KGlj6iqnd65cfhAsERwFvAj7wEdK1saKkxJqHg8NHhw6sHAXCtguxslywGFzS8lsGQIdWTxYMHJ96/uDkFgERDS/cbDNKVBQNjTI38HMGMGW6IaFERTJ8eOy7iWgrz5rkcQkYGj41cxKbusYRxTa2DZnqprOZA5hkYY0z680cM9erlAsHkyXDRRfDtt1XPu+iiqiOHgA8eLuZuqo4eao79/gfCWgbGmJahqMj1/YvAlCnu5geCbt1g7lzXZbRvn2sZ5ORUdgtJLSOHBg+G5csbvfZNJlHLwNYmMsa0DLm58Mwzrj/n9ttjeYLvfCc2uey55+AHP4BNmyjYewsitQeC/PyWFQhqY8HAGNNy5ObCDTe4QKAKffu6xeYmT3brDQEsWQJLl1JQ4E5J1Dmi2rwSwwfKgoExpnnz5xPMmuW6iu65J3Zs7VoXCPwcQh13LBs8uJHqmsYsZ2CMad6KimLDR7OzY/sSlJe751lZLn8Qtxx1QYHrAqppuGhLbhHY0FJjTMsTnE/w/e+7ABAKuSUm7rnH7Uswf77boCZudnEwV9BML4MNYkNLjTEtiz+fwN974NZb3XyCaNTlCm680Z3nLUPdkn/tJ4MFA2NM85JoPoE/eigUgvfeczkEPyBQ+74DzW256cZgwcAY07wMHOj2KK5pPsEJJ8BDD7kcws03u7JAQEikNXUTJWKjiYwxzUdhocsB+HsXB+cTHHwwfP21e+zNJ/hw9kv7nUtgHGsZGGPSW3Bjml69YvsSdOkCH33kykMht22lP59gwQJYsoQTcKtpQvWA4E8oS7T4XGtjwcAYk978biFvUTnGjYPZs2OBAKBtW9ddFJxPsJ/N7FtrbiARCwbGmPTlDx195hm37lBJiQsEQVdcAc8+6/YrqGE+ga+1LTxXXzbPwBiTnvLy3GJyCxfWvBR1KOQmlLVp44LAunU1zicwVdk8A2NM81FY6ALBQw+5FUZHjYI9e6qe07ataw1MneruFy2q1jUU3xVkXUOJWcvAGJN+iopcInj06KrdQiJuHaJ33nGPMzJqbRXEJ42b6eUuqaxlYIxJb8FRQ2PHuq6hiy92F/yKCnd/0EHw61+7c+bPd/cVFTBnTmrq3IJYMDDGpJ6fH/CHjT71lFs6dOfO2DmRiBtJ5A8dTRAACgoSzza2mcaJ2aQzY0xqBfMD48a5+QK7d8Pixe54KASHHupyBw895LqOalmKevlyd7GvafSQv3+BBYLqrGVgjEmtgQPh7rtjF/uOHWHXrtjx225zieExY9w5mzbBgw8mfLtXX3U3Uz8WDIwxqZWb67p9xoyBDh1g27bYsawsuP/+2DnFxfDgg5W/7Pc3g9jmFtSdjSYyxqSOPyCTNRQAABlxSURBVKksNxcGDIDVq2PHRoxwV3t/5nFg6Gh91xqyHEGMjSYyxqSP4DLUY8bAd7/rAkFWFpSWuudLlrhuIX/nsjosMRGvmf7WTYn9JpBF5BER2SYi7wXKjhCRZSLykXd/uFcuInK/iKwXkXdEpF/gNeO98z8SkfGB8v4i8q73mvtFbH1BY1q0vDzX7z9/vpssNnp0LFlcVgaTJsGHH8byA3PmUNB5DkOW3mIrkDai/XYTici5wB5grqqe5JUVAl+r6kwRuRU4XFV/KSJ5wH8AecAZwH2qeoaIHAGsAgbgFhFcDfRX1R0ishK4HngTWArcr6rP76/i1k1kTDMSnEMQDrtE8UUXwfPPu3kC/nXoiitg7lw36ay4uHISWX0CwODBsTyCdQ1V1+BuIlV9TURy4opHAkO8x48Dy4FfeuVz1UWYFSJymIh09s5dpqpfe5VZBgwVkeXAIaq6wiufC4wC9hsMjDHNSPzKoxddFGsN+M45xwWHoiLIzaXg1VwoqP9HLV+ejAq3Pg3NGXRS1c+8x58DnbzHXYFPA+dt9spqK99cQ7kxprnLy3NzBLp2dTOKn3nGle3dWz0QtGnjlpiYMqVyUtnUqfvPDxx7LFx5ZeNUv7U54ASyqqqINEmaRkQmABMAjjnmmKb4SGNMQ51/Ptx0E2RnuxnFl1/ulqAO8C8cxSdfTd935yG/msYfuk6h16+LgcTBwBLDydfQGchfeN0/ePf+wOAtwNGB87p5ZbWVd6uhvEaq+rCqDlDVAR07dmxg1Y0xjaqwECZOhNNOg3vucQFg9+7q+xCI8C3teJBJDFj1EI+UjuPx0kv518cVXPiSLUPd1BraMlgMjAdmevfPBsqvE5H5uATyTlX9TEReAH7tjzoCLgQmq+rXIrJLRM7EJZB/DPxPA+tkjEm1wkKXE5g/37UGFi2C7t1h48bYOSIuZ7B8OeyKMI4nmc015LCJ4SxN+NY2gaxx1WU00Z9wCeAOwBdAPvAMsAA4BvgEGONd2AX4HTAU2AtcpaqrvPf5CfAr723vVNVHvfIBwGNAW1zi+D+0DjPhbDSRMWkkOG/A33pyyhSXHwj+7yxCWcZBfL98CQBjcSuPbqAHd1N7a8C6hpIj0Wgim4FsjDkwhYXw+utuFJC/v8CUKfDtt7FzsrLgJz+B2bMpD2cTychi/qWLuGpu4rxAfEvAhokmh81ANsYkV7A18Nprrqy01M0N8GcNg5tXkJHhJpcdfzyrb3qKtyOnsmFu7Uliu/g3LQsGxpj6CQaB+fNdDmDcOPj9790EMl8oBAcf7H7iT53qNqpZtIizuLHK28W3APa3+JxpHNZNZIypn6IityexSM1dQr527dxyEzNm8Ld+k/nkxXW15gZsMbmmYTkDY8yBKSyEP/8ZLr3UDRsdNcp1C5WWVj0vI8PlCPzHU6ZARQXyy5qDQDO9BDVbiYKB7XRmjKmbgQPhvffg5pvhrbdg5MjqgSA7G+66y+UJwC0xUVFRbaN6k36sZWCMqbuiIhg+3A0ZDcrIgJ/9DObNi3UfrVvHso97JJxA5i8oZ11DTctGExlj6ie40ii49YXArTX00Ufucdu2cMcdLkH85JOu43/dOtcamDOHNwogf5BLCvtbUTbT358tngUDY0zNNmxwv/R98+a5/QbKy2Nl+/a5+2eecSOLKioYsm4OQzoDBS5GmObBgoExJibYGujVy/X9l5W5WzQaOy872y1Jffvt6E0388KFv2HFWXNgb+IN6W05ifRmwcCY1i4+AMyfH7v49+gBa9dWf80vfgE33ginncZz580i9OJLTH3xxurnBfithOXLbc+BdGTBwJjWLtgdlJHhJpD5K4zWFAiysuD++91+xLm5jKhlFnGQ5QrSmwUDY1qjRN1B+/ZVX2oa3GziaNR1D4XD7NtdQel5F3Mxi6htSQnTfFgwMKY18XcfO+881x0Uibhb+/awa1fNrxFxwWLiRNeCOPdc5j7XBYCBFLO8lmDg5wmsWyj9WTAwprUoLIScHPfLf/lyt1TErbe60UFxO5BVGjECli+nbF+E6CNPsnxwPuGSCq6pZblp25C+ebJgYExLFt8dtHChu8AvXgw334yqIoHTowhRCSEI0VAGkRdeZfngfD55cR0AG16sSLi2kE0ia94sGBjTEsWvLOovKX3MMbHN6L1AEMGtSxNF2MdB3K7TyGcqr0fOZWukixcA5iT8KEsMtwwWDIxpifwRQhkZruPeX1k0bnSQ4gLBs4xgCMvJIEI+05hKPpkkbgWYlseCgTEtSbBFEA67ZSF+9atqC8pFEcrJ4COOpw9rOZ+XuZ1p/BuuO2h/gSA/3/YdaGksGBjTEiTacGbOnKozh3GBoCLzIH5/9DQu/XgGs5lEf1aTSQXX1NIdFGR5gZbHgoExzZm/x0D//tWDQE3zBUaMYNfi5WSUR/jRx9OYyhQyqeA6HtzvR9lyEi2bBQNjmpP4lUR79YJ33oGVK6F3b7ea6OzZKFSOEtoZPpx2ugeVMJEXXmUq+ZXdQXnnV/DGoFvwr/O1df1Ya6Bls2BgTDoLXvzXrHEtgHnzYjkA1cpuIPWSw+LdFCgjA41E+QUzvRFC51TtDnrJ3QYPtolhrZ0FA2PSVWGhGw0UnCm8cmXl0hDBEZ0qYULqho/65SVk8yvuJJ+pCUcI2bBQ47NgYEw6iO/+WbMGvvMdePFFOPVUWL0aKirchT4aRRFCgXAgGql8psBib6ioHwT+jXU2VNTUyoKBManiJ387doytFVRWRqS0DNEoAkQlTGjFisqLv58HEJRIXECIEGIOExnP44xgCbO5hgwi1UYJ2XpBpiYHtAeyiGwCduMmMVao6gAROQJ4CsgBNgFjVHWHiAhwH5AH7AWuVNX/895nPPDf3tveoaqP7++zbQ9k0+zUlPy97TYoKSESyuBfbXqSU1J9yWh/Ylj8xR8gQhghSgUZlNOGCjKYyhTOo4gQEYaztPr7WddQq9aYeyDnquqXgee3Ai+r6kwRudV7/ktgGNDTu50BzAbO8IJHPjAA9+9+tYgsVtUdSaibMemhsJB1f3id7puKANf9nxGKEKooIwSEohWVgUDiXipAFNcaANcCUEIIUUJEmM0kMohUnp9JBavzlwDufyybHGbqojG6iUYCQ7zHjwPLccFgJDBXXVNkhYgcJiKdvXOXqerXACKyDBgK/KkR6mZM0wpMBuv00WuUAWEitKGMkDcKKEqIMLGJYeWEyAg893/Ib+MoDmYPtzO9cmhoF7aSwyaGs7TKrGAbBmrq60CDgQIviogCc1T1YaCTqn7mHf8c6OQ97gp8GnjtZq8sUXk1IjIBmABwzDHHHGDVjWkEXh5g3Y6O7G7flRXf9GLCp/MIR0r5hs504TMyqajy6z9M1ZFBGUQpJ0wmEcoJU0Y2YSIcxXZmc01lDiA/H/z/0fKxAGAOzIEGg7NVdYuIHAUsE5F/Bg+qqnqBIim8YPMwuJxBst7XmAZLMAnshJISIoQ4xVsPNESUY73fPH4giP8H/D696cWHhIkQRnmWEXzBdyqPd2ErFxy/iSfHPcjg5XbxN8l1QMFAVbd499tEZBFwOvCFiHRW1c+8bqBt3ulbgKMDL+/mlW0h1q3kly8/kHoZ06iCAWDrVigqgkiEipIKJAREo4TA6/qp3hUUDAJRBEUopw3d2MotFHJNzyIO+nYb7foM4rNBt1R2/XwGXFQABU30Z5rWpcHBQETaASFV3e09vhCYBiwGxgMzvftnvZcsBq4Tkfm4BPJOL2C8APxaRA73zrsQmNzQehnTKLwAsGoVrPimF1dvmU9GRRkhrQCNEtIoYfCv/dXyAMGuoAoyKGYAJ/MOIBSRS+d+Xdi6FfJOqqDXMpf8PRq4oAn/RNO6HUjLoBOwyI0YJQN4UlX/JiLFwAIRuRr4BBjjnb8UN6x0PW5o6VUAqvq1iEwHir3zpvnJZGNSJr77Z+tWeOklTi6poC8hPqIHJ1J1GGhteQCAqGSwuv/P6PvuPE6uWMvfvzedDl+uI/uIHgxYZpPBTGod0DyDVLJ5BiaZCgpiI3EGvVFIJJTB4FenQnmEcBhCkTLC0fJqr4sfBhr/f1MFGd4Q0Cih7GzIynIbzaxbBz16wC0WBEzTasx5BsY0W34Sdu/UQi5jA7wKJWzlHF7nEcbxU/5AOBohFOjyCQaACoRwXAiIIm4paUJUSBavnj+N7huL6HX4Nre0REWFW2LamDRiwcC0anunFtKDDfyTXlzKfDKIECZCJmVcy+xqff8+//IfRr1JYG6t0ArvHQ7+QS506UJbYGiPCrhlSVP+WcbUmwUD02wFh1YmmmUbLB/0RiGHf72BrVuhSxeXBhjF2wxgNUqI1fRjAKvIpKLy9TX1/UNsFjBAKVnczjRyv7OOLl3c8QHnWBeQaV4sZ2CaLYnvsK/Fc+SxiRzGMa9y6YY2lFa58Fe+b9xzxV38w0SJSogK2iAZYQA25uSyu30XdhzRgwssCWyaAcsZmFbjF7iun6ByMrmW2bxPb3qyngwqqi365geB4IJwfgsgijCHifz8ZxE3hNTTy5LApoWwYGCalYICmDq19nN6sIFxzKt8HiZCNiVEEE5kbeUqoDVxx5RIKJNIOJNoRAiHXaA4o2/EEr+mxbJgYJqdwYNdHsBP/gZ1YQtH8hWZlCHeb3o/ARxGKwOBWwU0RgBEkLD73Z+R1YaMadPcEFDPgB49GvGvMia1LBiYZmXqVNcN1PnVDXRhK7m4JaH9EUD+hb+MzGoLwkFwOWgoJ4M2fs4gMxPatImdmJtrQ0BNq2LBwKQdf5RQcCTQ5fPyUAnxEF0rh4EeHNpHKBpBBUJadfhnG2ITxKrsFYwQkQy+PvJ4jvpyLfTu7SZ/+cOAfJYLMK2MBQOTFoYMiV34/V///iQwgGXkMInZ9CBEhAy28h0Oje4CQLT6CKBKoRDi7RtAOIxkZxMCOpVvhUmTYNMmWGJzAIyxYGCaTHDJh6BBbxRy2avuwn8Ka+hPR/7FMVzJo5X9/oIb2ZNBlDBl5PAvRKT2PRxFIByG00+Hd95x+wwHWwA5OfDgg0n/O41pjiwYmEYTv96+Pwro9Fdd4vcU1rCdjizlPKYzjzaUEvYGc0Ks398XgsqN4YHaA0EwB7B2LUyf7nIA1vVjTI0sGJikCgYA/+Lvj/t/yCv/J70Yxzyy2EcGUYbxN6JAZmAfX4j1+1cd9aOVCeDK8t69Yf36WHAYMABOPrlqxSwQGFMrCwam3uJ/8Qe7foLDPf2Lf19vyQefu4hHK/f5dTt7BY9VFdwT2C3/Bhx7LHzyiWsBbN0KM2fGhoFa8teYerNgYOot/hf/ZVCZ6I2/8AfH+ceraRG4+M3gwe0JrOI6h8LhEGRkwDffuATw6tW2EqgxSWDBwOxfXh6EQtC1K+B+8X+fv5HNPg7nG8Trww976/bXJn4DmPgloDO8ZDHgPhMgGkWy2sCdd1aZBGYJYGOSx4JBK1Xbip+D3ihk33sbKgfeREM5DFw1G7eqf4ifBi765YSr9fX7Eg33DKZ9wyhRCRMO4fr8o1EXBCZOhIj3vlu3wrZt9uvfmEZkwaCFqWnCli9R377fzXMOr3EQe9lKF85lNfJ5bFhnReWwzqq//DOJVB3hEydYKsH7UGx1oHDbbPCXfvAv/Dk51u9vTBOyJaxbmPhlnWtawbMLW+jLO3Ths8qy+L79MjKrzOIFar3oRxIcq1KdUMjdMjLc+P/c3Ni4f0v6GtMkbAnrZmB/m7XUVjbojUJ6f/BnFtORrbi+/S5sYSDFdORLIt46nbUldIMSDetMdNEPoURDGYSj3rwA/8Lvq2m4pwUAY9KGtQyaQF27bhItzRz/6/4U1hAlxGF8U1kWRejNP93yy96Fvy4XfahlKQdqHt0jGRmu/x5is3yjUXeztX6MSWvWMmgiNS25ELzIv/pq7HH8+jv/W8NF/jt8ThmZcaN2IjX+Oi8nTIhIwiBQ20U/0U+CDKJEvU/zt3nMFIGjjoKsLBg2LHby1q0u6Wtr/RjT7Fgw2I/aftXXVF7Tkgv9va6bLmyhBxv4Dp8Dyj/pTX9WB9bfiSbcdKW2UTu+TCJECFdu61iTurYDI4QQbzXQaLgNL3/vTjp8uY6tW2HECOyXvjEtjHUTJeAHgZq6bp4jjyihyr55cOPuD2Enn9MZgG84jAGsqvwVHyVU4xj8mhK1idSWwK163v5VGcsfGM9Phw5uYTe/m8cf3fPDH9rF35gWIO27iURkKHAfEAb+oKozk/oBhYWwYQOsWcO6HR3Z3b4rpZu2cIJsIKN8Lwft+4q9bY+kIvMgvt0LY/a5Lprr0coLPMBB7CWbfXRie+W4++Bkq8PZWXmuErvoJpqMVVOiNpHaErixc7z3ysyE8kCQiU/oduniErrWt2+MIU2CgYiEgQeAC4DNQLGILFbVD5L2IQMHuhmsZWX0LCmhvIZf5G3KvwWgQ9xLgxd4n99tU1uS1i27XHu3Tez9qidqaxJCiUqo2mYu4AWf4EX/zDPdL/uDDoJzz616sl30jTEBaREMgNOB9ar6MYCIzAdGAskLBrm58MwzMGoU5SWROnfNJLK/yVa+DCL7/TXvzou6/XcjCQJHYNROOBp1C7W1bRs7fthhbhOXU0+NldkF3xhTR+kSDLoCnwaebwbOiD9JRCYAEwCOOeaYen1AQQFMnZrLVG5gCtMbXtNgferQbQNuyQWCwzGDgv31qm5oZjCP8/nn0L59zaN2li5Nwl9hjDHpEwzqRFUfBh4Gl0Cuz2sLCqBgcBGMuo/SXXVP2u6PP9kq5E22qtb3LxK7nXmmW23z88/dBf+ccyxRa4xJC+kSDLYARweed/PKkqeoCEaNgrIyMimv8yiehIldv9sGCItAh6OgtBQ6e8nm449399u2xbpurNvGGJOm0iUYFAM9RaQ7LgiMBS5P7icUw9ixsGYNH9VhNFHOscR+wXeOjSZi715LyBpjWpy0CAaqWiEi1wEv4IaWPqKq7yf1QwIX6l4JTsny7jsm9YONMSb9pUUwAFDVpYBlRI0xJgUSrX5gjDGmFbFgYIwxxoKBMcYYCwbGGGNoxquWish24JMGvrwD8GUSq9NYrJ7J11zqavVMvuZS18au57GqWm3QZLMNBgdCRFbVtIRrurF6Jl9zqavVM/maS11TVU/rJjLGGGPBwBhjTOsNBg+nugJ1ZPVMvuZSV6tn8jWXuqaknq0yZ2CMMaaq1toyMMYYE2DBwBhjTOsKBiIyVETWich6Ebk11fWJJyKbRORdEVkjIqu8siNEZJmIfOTdH56Cej0iIttE5L1AWY31Eud+7zt+R0T6pbieBSKyxftO14hIXuDYZK+e60Tk+01Yz6NFpEhEPhCR90XkBq88Hb/TRHVNq+9VRLJFZKWIvO3Vc6pX3l1E3vTq85SItPHKs7zn673jOSmu52MisjHwfZ7qlTfdf3tVbRU33NLYG4DjgDbA20CfVNcrro6bgA5xZYXArd7jW4G7UlCvc4F+wHv7qxeQBzyP2xfoTODNFNezALi5hnP7eP8GsoDu3r+NcBPVszPQz3vcHvjQq086fqeJ6ppW36v33RzsPc4E3vS+qwXAWK/8IWCS9/ha4CHv8VjgqSb6PhPV8zHgkhrOb7L/9q2pZXA6sF5VP1bVMmA+MDLFdaqLkcDj3uPHgVFNXQFVfQ34Oq44Ub1GAnPVWQEcJiKdaQIJ6pnISGC+qpaq6kZgPe7fSKNT1c9U9f+8x7uBtbh9wNPxO01U10RS8r16380e72mmd1PgPOBprzz+O/W/66eB74lIwo0Nm6CeiTTZf/vWFAy6Ap8Gnm+m9n/UqaDAiyKyWkQmeGWdVPUz7/HnQKfUVK2aRPVKx+/5Oq+J/Uigmy0t6ul1T5yG+4WY1t9pXF0hzb5XEQmLyBpgG7AM1yr5RlUraqhLZT294zuBI1NRT1X1v887ve/zXhHx99pqsu+zNQWD5uBsVe0HDAN+LiJV9tZU125Mu7HA6Vovz2ygB3Aq8BlwT2qrEyMiBwN/Bv5TVXcFj6Xbd1pDXdPue1XViKqeittD/XTg31JcpRrF11NETgIm4+o7EDgC+GVT16s1BYMtwNGB5928srShqlu8+23AItw/6C/8ZqF3vy11NawiUb3S6ntW1S+8//miwO+JdVmktJ4ikom7uM5T1b94xWn5ndZU13T9Xr26fQMUAWfhulX8HR2Ddamsp3f8UOCrFNVzqNcdp6paCjxKCr7P1hQMioGe3uiCNrik0eIU16mSiLQTkfb+Y+BC4D1cHcd7p40Hnk1NDatJVK/FwI+9URBnAjsDXR9NLq5/9WLcdwqunmO9USXdgZ7AyiaqkwD/D1irqrMCh9LuO01U13T7XkWko4gc5j1uC1yAy28UAZd4p8V/p/53fQnwitcaS0U9/xn4ESC4vEbw+2ya//aNlZlOxxsuM/8hri/xtlTXJ65ux+FGYbwNvO/XD9eP+TLwEfAScEQK6vYnXFdAOa7P8upE9cKNenjA+47fBQakuJ5PePV4B/c/VufA+bd59VwHDGvCep6N6wJ6B1jj3fLS9DtNVNe0+l6Bk4G3vPq8B0zxyo/DBaP1wEIgyyvP9p6v944fl+J6vuJ9n+8BfyQ24qjJ/tvbchTGGGNaVTeRMcaYBCwYGGOMsWBgjDHGgoExxhgsGBhjjMGCgTHGGCwYGGOMAf4/3EDf5b55ygoAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(eigh, \"+b\", label=\"glt symbol\")\n", + "plt.plot(eig, \"xr\", label=\"eigenvalues\")\n", + "plt.legend(loc=2);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.8" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}