Commit 53f479f5 authored by Jeremy BLEYER's avatar Jeremy BLEYER

Started contact penalty example

parent 06f459d0
......@@ -12,6 +12,6 @@ Contents:
demo/viscoelasticity/linear_viscoelasticity.ipynb
demo/2D_plasticity/vonMises_plasticity.py.rst
demo/contact/penalty.ipynb
......@@ -24,7 +24,7 @@
<script type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
<link rel="index" title="Index" href="../../genindex.html" />
<link rel="search" title="Search" href="../../search.html" />
<link rel="next" title="Beams and plates" href="../../beams_and_plates.html" />
<link rel="next" title="Hertzian contact with a rigid indenter using a penalty approach" href="../contact/penalty.html" />
<link rel="prev" title="Linear viscoelasticity" href="../viscoelasticity/linear_viscoelasticity.html" />
</head>
<body>
......@@ -35,7 +35,7 @@
<div class="rel" role="navigation" aria-label="related navigation">
<a href="../viscoelasticity/linear_viscoelasticity.html" title="Linear viscoelasticity"
accesskey="P">previous</a> |
<a href="../../beams_and_plates.html" title="Beams and plates"
<a href="../contact/penalty.html" title="Hertzian contact with a rigid indenter using a penalty approach"
accesskey="N">next</a> |
<a href="../../genindex.html" title="General Index"
accesskey="I">index</a>
......@@ -431,6 +431,7 @@ when considering a zero hardening modulus.</p>
<li class="toctree-l1 current"><a class="reference internal" href="../../nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul class="current">
<li class="toctree-l2"><a class="reference internal" href="../viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2 current"><a class="current reference internal" href="#">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="../contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../../beams_and_plates.html">Beams and plates</a><ul>
......@@ -465,7 +466,7 @@ when considering a zero hardening modulus.</p>
<div role="navigation" aria-label="related navigaton">
<a href="../viscoelasticity/linear_viscoelasticity.html" title="Linear viscoelasticity"
>previous</a> |
<a href="../../beams_and_plates.html" title="Beams and plates"
<a href="../contact/penalty.html" title="Hertzian contact with a rigid indenter using a penalty approach"
>next</a> |
<a href="../../genindex.html" title="General Index"
>index</a>
......
......@@ -241,6 +241,7 @@ writing/reading. Prefered output format is now .xdmf:</p>
<li class="toctree-l1"><a class="reference internal" href="../../nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="../viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="../2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="../contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../../beams_and_plates.html">Beams and plates</a><ul>
......
......@@ -2826,6 +2826,7 @@ if (IPython.notebook.kernel != null) {
<li class="toctree-l1"><a class="reference internal" href="../../nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="../viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="../2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="../contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../../beams_and_plates.html">Beams and plates</a><ul>
......
......@@ -219,6 +219,7 @@ that:</p>
<li class="toctree-l1"><a class="reference internal" href="../../nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="../viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="../2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="../contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../../beams_and_plates.html">Beams and plates</a><ul>
......
......@@ -282,6 +282,7 @@ and the beam theory eigenfrequencies :</p>
<li class="toctree-l1"><a class="reference internal" href="../../nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="../viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="../2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="../contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../../beams_and_plates.html">Beams and plates</a><ul>
......
......@@ -2162,6 +2162,7 @@ if (IPython.notebook.kernel != null) {
<li class="toctree-l1"><a class="reference internal" href="../../nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="../viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="../2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="../contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../../beams_and_plates.html">Beams and plates</a><ul>
......
......@@ -216,6 +216,7 @@ so that no term arise in the linear functional):</p>
<li class="toctree-l1"><a class="reference internal" href="../../nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="../viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="../2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="../contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1 current"><a class="reference internal" href="../../beams_and_plates.html">Beams and plates</a><ul class="current">
......
......@@ -209,6 +209,7 @@ the thin plate solution.</p>
<li class="toctree-l1"><a class="reference internal" href="../../nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="../viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="../2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="../contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1 current"><a class="reference internal" href="../../beams_and_plates.html">Beams and plates</a><ul class="current">
......
......@@ -4348,6 +4348,7 @@ if (IPython.notebook.kernel != null) {
<li class="toctree-l1"><a class="reference internal" href="../../nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="../viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="../2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="../contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../../beams_and_plates.html">Beams and plates</a><ul>
......
......@@ -2912,6 +2912,7 @@ Mechanics and Engineering, 85(3), 349-365.
<li class="toctree-l1"><a class="reference internal" href="../../nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="../viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="../2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="../contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../../beams_and_plates.html">Beams and plates</a><ul>
......
......@@ -84,6 +84,7 @@
<li class="toctree-l1"><a class="reference internal" href="nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="demo/viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="beams_and_plates.html">Beams and plates</a><ul>
......
......@@ -92,6 +92,7 @@
<li class="toctree-l1"><a class="reference internal" href="nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="demo/viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="beams_and_plates.html">Beams and plates</a><ul>
......
......@@ -76,6 +76,7 @@
<li class="toctree-l1"><a class="reference internal" href="nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="demo/viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="beams_and_plates.html">Beams and plates</a><ul>
......@@ -123,6 +124,7 @@
<li class="toctree-l1"><a class="reference internal" href="nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="demo/viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="beams_and_plates.html">Beams and plates</a><ul>
......
......@@ -98,6 +98,7 @@
<li class="toctree-l1"><a class="reference internal" href="nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="demo/viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="beams_and_plates.html">Beams and plates</a><ul>
......
......@@ -58,6 +58,7 @@
<ul>
<li class="toctree-l1"><a class="reference internal" href="demo/viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l1"><a class="reference internal" href="demo/2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l1"><a class="reference internal" href="demo/contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</div>
</div>
......@@ -93,6 +94,7 @@
<li class="toctree-l1 current"><a class="current reference internal" href="#">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="demo/viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="beams_and_plates.html">Beams and plates</a><ul>
......
......@@ -107,6 +107,7 @@
<li class="toctree-l1"><a class="reference internal" href="nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="demo/viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/2D_plasticity/vonMises_plasticity.py.html">Elasto-plastic analysis of a 2D von Mises material</a></li>
<li class="toctree-l2"><a class="reference internal" href="demo/contact/penalty.html">Hertzian contact with a rigid indenter using a penalty approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="beams_and_plates.html">Beams and plates</a><ul>
......
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -12,6 +12,6 @@ Contents:
demo/viscoelasticity/linear_viscoelasticity.ipynb
demo/2D_plasticity/vonMises_plasticity.py.rst
demo/contact/penalty.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Contact of an elastic membrane with a rigid indenter: an augmented-Lagrangian approach"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Residual: 7.55796079282e-06\n",
"Residual: 6.04663326501e-06\n",
"Residual: 4.88782492198e-06\n",
"Residual: 4.05065021525e-06\n",
"Residual: 3.32603080899e-06\n",
"Residual: 2.74757531832e-06\n",
"Residual: 2.28568557151e-06\n",
"Residual: 1.91675781292e-06\n"
]
}
],
"source": [
"from dolfin import *\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"parameters[\"form_compiler\"][\"representation\"] = 'quadrature'\n",
"import warnings\n",
"from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning\n",
"warnings.simplefilter(\"once\", QuadratureRepresentationDeprecationWarning)\n",
"\n",
"ppos = lambda x: (x-abs(x))/2.\n",
"Heav = lambda x: (1-x/abs(x))/2.\n",
"\n",
"N = 40\n",
"mesh = UnitSquareMesh(N, N, \"crossed\")\n",
"\n",
"obstacle = Expression(\"y+pow(x[0]-0.5,2)+pow(x[1]-0.5, 2)\", y=-0.001, degree=2)\n",
"\n",
"alpha = 100.\n",
"Nitermax = 10\n",
"\n",
"V = FunctionSpace(mesh, \"CG\", 2)\n",
"#P = FunctionSpace(mesh, \"Quadrature\", quadrature_scheme=\"default\")\n",
"Pe = FiniteElement(\"Quadrature\", mesh.ufl_cell(), degree=2, quad_scheme='default')\n",
"P = FunctionSpace(mesh, Pe)\n",
"\n",
"u = Function(V)\n",
"u_old = Function(V)\n",
"du = TrialFunction(V)\n",
"u_ = TestFunction(V)\n",
"h = Function(V)\n",
"\n",
"p = Function(P)\n",
"p_old = Function(P)\n",
"gap = Function(P)\n",
"\n",
"corr = Constant(0.)\n",
"rho = Function(P)\n",
"dxm = dx(metadata={\"quadrature_degree\": 2, \"quadrature_scheme\": \"default\"})\n",
"lagrangian = 0.5*inner(grad(u), grad(u))*dx + p*ppos(h-u)*dx + alpha/2.*ppos(h-u)**2*dx\n",
"form = inner(grad(u), grad(u_))*dx - p*Heav(h-u)*u_*dxm\n",
"J = derivative(form, u, du)\n",
"\n",
"def local_project(v, V, u=None):\n",
" dv = TrialFunction(V)\n",
" v_ = TestFunction(V)\n",
" a_proj = inner(dv, v_)*dxm\n",
" b_proj = inner(v, v_)*dxm\n",
" solver = LocalSolver(a_proj, b_proj)\n",
" solver.factorize()\n",
" if u is None:\n",
" u = Function(V)\n",
" solver.solve_local_rhs(u)\n",
" return u\n",
" else:\n",
" solver.solve_local_rhs(u)\n",
" return\n",
"\n",
"def border(x, on_boundary):\n",
" return on_boundary\n",
"bc = DirichletBC(V, Constant(0), border)\n",
"\n",
"def solve_Uzawa(tol=2e-6):\n",
" t = 1.\n",
" for i in range(Nitermax):\n",
" u_old.assign(u)\n",
" p_old.assign(p)\n",
" told = t\n",
"\n",
" local_project(p + alpha*ppos(h-u), P, p)\n",
" solve(form == 0, u, bc, J=J)\n",
" #gap.assign(project(h-u, P))\n",
" #p.vector()[:] = np.minimum(0*p_old.vector().get_local(), rho.vector().get_local() + \n",
" # alpha*gap.vector().get_local())\n",
" #direction = assemble((h-u)*(p-p_old)*dx)\n",
" #t = (1+sqrt(1+4*told**2))/2.\n",
" #if direction >= 0:\n",
" # corr=(t-1.)/told\n",
" #else:\n",
" # print \" Restart...\"\n",
" # corr = 0\n",
" # t = 1.\n",
" #rho.vector()[:] = (1+corr)*p.vector().get_local() - corr*p_old.vector().get_local()\n",
" res = max(0*errornorm(u, u_old), assemble((p- p_old)**2*dxm))\n",
" print \"Residual:\", res\n",
" if res < tol:\n",
" break\n",
"\n",
"h.interpolate(obstacle) \n",
"solve_Uzawa()\n"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAAEACAYAAADyRL7nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvXmQHMd95/vN6nvOBgacwUFcAxI8RclDciVZ1q6eBFpe\n0WtF2IC0lFaUQxKBcMjxQtSLAOXYPzZi40X4gS/W1IuQvYGhJVvSo/AoUM9BhSmLAihLlrQLWgR0\nAQRJEY0BQADEAHNhjr4r94+qmqmuriOzKuuayU9Ex0x3ZmX9uiv727/65S8zCaUUEolEIgkPJW4D\nJBKJZLUjhVYikUhCRgqtRCKRhIwUWolEIgkZKbQSiUQSMlJoJRKJJGSYhJYQcoyhzigh5CAhZI/+\ntxzcPIlEIkk/xC2PlhCyB8AogMOUUuLaECEnKaX36/+XATxNKd0n0liJRJJ+CCGjAPYCOAVgDMA4\npXSWt27UZYGglHo+tGqu5WMAjllem2FpWz7kQz7W1gPASdP/ZQBH/dSNuizIQ1SMdhSAVfWnCSFj\ngtqXSCSrAF0Tpo3nVPMW9/DWjbosKKKEdr2gdiQSyeqGxylzqxt1WSCyQRvQmYbmZptxFF9CyH4A\n+wGgVMzdv3XrRgBAPrcARWmiVl8HAMhk6ijmZ7FYHdGPo+gpTqJaXw9VzQEASoUptNolNFs9ehvz\nUEgbtYZmTjZTQyF3E4u1Yb0NFT3F66jWh6Cq2tvvKd5As9Wz3EYhfxMARb0xqLWRrSKfXcBidQSE\nqCCkjZ7iDSzVNoDSjN7GdTRafWi1SnobcwAI6o0BAEAuu4RcdglLtQ0AAEVpoVSYwlLtFlCq/d71\nFidRbw6g1S4CAIr5Wag0g0azf7mNbKaKan1Ib6OJUmEaS7VhGCH03tI11BpltNsFUKqgVJyCqubQ\naPbpn88iMkod1fp6/TNuoJifWf6MV9pYh3Y7r3/G02irBTSavaFdJ0JU9BRuCLlOS7Vb9DbCu05L\n9Q3Lz4NeJwAoFmZCu06/Pr1wg1J6CwLw0P9WolPTKlPdX/y6cQZAzfTSOKV0XP+fxylzqxt1WSBE\nCW0FNkZSSk/ZVdY/9HEAeOc7+umJH2T0kkFLzR79YWaz5fmw5blV73v1h0HGpo0Ry/N1lud9APow\nt7ANg30XHdrYaHlu/Tj69YebHZssz4d8tLHyXLPXOoY5YHle0h/2bWhYv6Pir9Pcwk4M9uVMr/i/\nTu5tiLlOK30BCHqdNDZYnou7Tj2bFy4gIFPTKn72fev57OnZPFGjlD7gUMzjlLnVjbosEL5DB3o6\nVxnoFlR95O44SzuGx5MGFqrWL3qySZO9abIVSJ+9CYLHKXOrG3VZIFyFlhAyRgg5qP9/SE/3MjgE\n4GOm548ZebTQ0iMeYzHAuFVKA5euvT9uE7hIk71pshVIn71JwcspY3Xgoi4LimvoQD/xKQBP2pTt\nc6gLUcZJJJJVyWO6A2fkqpqdskMAjkEPLXrUjbrMN6JitL7J5xbRHY9KJrcO/8+4TeAiTfamyVYg\nffYmCTenjMeBi7osCLELbUapx20CM4N9E3GbwEWa7E2TrUD67A1KGxRzanq+q0kj9kVljPSVNHCm\n8kjcJnCRJnvTZCuQPnsl8RK7RyuJlyi8lEGlEPo5JJIkE7vQZjINdOcHJpOBvktxm7AMi0AWeycS\ncbsnytYkCXaS+oIk+cQutMX8DNIitHfvOBr6OUQK487tR4S1FTYstrJ+NlEIchR9QbJ6iD1Gu5ii\nxO8TZx4X0s6cWnd8iOT02YNC2wsTkbZG8fmK6gtpoU2BOZXtIekmdo82VbgvydtF7LftNPbfUXYi\nstXpmnB7wZx9QbK2kULLA3FeJB1IgLBaISlyL2K21XrtPIXXoy9IJGZcd1iIgrF3FijrYhVJI3HC\nKgmNJA3E8dKzeeKkyyIvTNx7X55+5wXrwjf23LntauDzrTZiv7esNawrMCWXVyf2hRZPDYPzF9KT\n65l0W63X/dUJuUuThJ3YQwfGWppJxvhyzSzcii0x28LD4uJ21/IoBy4GPX7SvWxNEnNqHTMLt2JO\nrafa05VER+xCm1TS4LHaYRbPNpIzCuxlh9lWL1FOEkY/We2C2wLBlFqM24zUErvQlgrT6F6wOD7c\nBHZ0xzcitMQZVvHcuD0Z9rJgttXt/SVFhK19Ya0IrsQfsXfbtpqMjskSd11YHI3ImhWC5CnWYrDX\nL6y2JiVv06kvpCV+L4mW2IXW2N8oTli/GJPXfy9kO8SKyOyNcO0VSRBb4xBer74gBVdiJvbQQZwk\n4YuQlBjqasL6mcYZbpAhBQmQAKHN5xbQvdlf+PgR2ZHhHwk4b+AmmFl3y4+iO1lAwrTV/JmLEl3e\nvpD2DIU2VTCrWjfglLASu9AqSjPS8wXxYkvFqz7P6fuUgcj7tDcOorJVlOj66QtpF1uJf2IX2lo9\nugkLQUMFExcfwb13/18c5wt0usBcu/QIdtzZaW8SUnSGlFrXa3a2hk0Q0eXtCyvnlKEEVvTNEfdi\nZf+ucUrpLG/dkMoOATgIYBbAKwAOUEorTu8ldqGNiqjisXGLK7Aipi2qJEJYrdjZZNhqJ8JREHUO\nr/RumThKKb0fAAghrwB4GoDTlDy3umGUnaOUfWWh2IU2k6kDCDf2I0pk+/vedDmHkFNw4yakuT7H\nH9jEYdhq936iFF9WwXXrC+znkmLrBCFkDMC08ZxSOksI2cNbN4wyP8QutMX8LMIUWpGe7Lat/79N\n+8KaZ4LHQ+3d/N0QLRGLm63W9xyF8HoJrl1f8HeeVSm2G3QP0GCcUjruWNueUWi35WamCSFj+k61\nTHXDKNPPXyaE7NXrPATgL53CGkAChDbMhb9FhwvOnD3YEZeLQmSD3PrPvvEFrLvjrwRaEx48tkYp\nvE6Ca+0Lq50WMphq9bFWvyFg9S6eXVvd6oZRBnTGa6cBvATgfqfKsQttGglbYJMYV00y5s8rLNEN\nM4a7Sr1aWwgh+wHscqlyjFJ6HNpte9lS5iR+bnXDKIPZe6WUniKEjBFCyk5ebexCS0JaQDmMwS9V\nqYcmsmGIK1Hin5DBiihbjc8xTMEdVABF8Ge7VsSWI4RQgY2w2oQNXOsSQhBC2RiAp42BMlNZckMH\nPcVJAGIX/g5DZOdUYPvup4S3G6b3Wr79r13Lo0xALytL7uUetvISppc7pwJbQugLa0VsWTCJHYDl\nVKvjlufTlNJZt7phlEET9sOmsj0AnnN7P7ELbbXOE4rxJiyRBYArE49is6AVvKIID9y88AjUrc+H\nfh4WnETdEOCbFx7BQEi79obh5V6ZeBTY8Y3ErCa2SnmMEHIQK3msj5nKDgE4BmCcoa7QMj0DoaKH\nQQAtFGI+rovYhVZVc3Gb4Ig1TNCoBfe8wxZYs6C1alviv8AeGPa2aluW//fyfv0iUnCNvmCEEkSR\nVK+2TRXMtaOdgquHCYxQwXFL2T6OumGUdTz3IunfQy5EerOiY7FhCexqnH9ufU+ihVe0hytabCWr\nj9iFtlSYAjActxkdOIns5h1f5W4rDIFlFdfMtmeFnzss3Gw1v1+RohtEcK19QWRWQlK9Wol/Yv8d\nbrVLQtoR5c26ebILN+9lbmdKLQoV2Vm1Z/nBijq/W9j5w4bVVj+fgxd+rpVTX0jCFGxJ8ohdaJut\n5Nz6en1Jbk6/m6kdUQIbVFTozO8EOv9Uq4/5ERQ/toYhuKy49QURYpuEtZIl4og9dCACEZ1SxJdD\npMCGjQhx5GlvKLsg9HxmRIYWws7DlaxNPIWWc6myMQDG1LsygOfclg4DgHxuHt0TMKKFVWTXDx9z\nLBMhsqIFVrnlJwDEi6ofnGwwBNiwNSiiMhe8VhJz6wuAmAGyJMVq48g6WE2weLQ8S5XtoZQ+aTwh\nhBwGcMCtcYW0GU21J8pbrGxuzvb1oCIrWmANUcuQFtoJEFk3rLaK8nxFCK6b2Dr1BTMyG0Fi4NoN\n7JYKA+C2VNgBQgiXe1prpMObBYDJy3s7ngcd8BIZY7SLlxavfUhI21Fg2Coy7gsE/xFzur7WvuBE\n0JCUjNWuDrw8Wp6lygBttsZ5QsgT+vMnbOokhiBfgiR4sUkICYSJ+f0F8XSDerdB47bSs5V4CS3X\n/FhK6bju0RrhguPoFmpjBZ/9ADAyUsaJ0/8VALB15CfoK13D2QnNWyj3V7B723fxr2e+AADIZBp4\n8K6v4HTlE1hY2og2KG4b/Spm5+7FjSltFHjTxmPI5eZw8ZLWxsDAWWzZ9H2cff1x7Q1nF3Dn7q/g\n1XOfRbN+CwBgy+hh3Jx+APOz2hoRQxu/B0Vp4vqVjwIA+gZ/g/KGH6PZKGPitS9Bzd7E4OjXMFf5\nDNSm5pEP7jqM6o3fQ2PuHgBA76YXQGkWS29/GABQKP8ShXWnMHP+c9qHkJ9CdvuzaFU+DbS1Ldcz\no38LdfIDoAu3AQCUTd8DWv1Qr79f+9zW/QI3S2+hdFmzq1i4jtqWF9Az8UkQfYbd4s6vo3Dtg8gu\nbQVp9kOpbkKmvgH5ae29Ndb9Eu3SZZSuPAwAaJeuorbpRfRWPg2AAKBYHP06ilc/jEx1EwCguvkF\nZKpbkJ95l9bG+pNoF26gdFV7b62eS6iP/BC95z+t9QOliaUdz6B4+WFk9M+4uuV5ZBd2ITenpUXV\nh14GzS4se7IUbUDNoXfik9rzTBVL259F6a2PQmmsQw3A0tbvYHD+NtC5d2ifz/APAaUF9e3f1z6f\ngdegDJ1A+/yfap9xbg7ZHc+gNfFJoDmIGwAyO/8e+en7Pa/TzfOf0a5J4ToGdnwTs+f2Y6bVhyxR\nse32pwC0MfHalwAAw1ueQ6s5iOnJh7Q+t/5l9A2cxpWJz2qXungFg6PfwKuvPQ5V1WKu99z1JC5e\n+mPM69d6x7YjqNY24drkB7Q2b/kp+norqEw8igwIBvou4e4dR3HizOMAJQCheM89T+HViX24ubBV\na3P0COYWduCtyfcCWPk+AROQxAuh1Hn1LH1h2wOU0odMr80A+JCdR0sIOWjEaHUxfYJS6rYkGsbu\nK9KfvbjJl/FBbqv8eLNqu4AZwrx7RRdBvFhf3quaAxg2vwxzkGMww+hFMtpqJoiXGyR2O6TUoLYL\nUDJ8/S+IVxtkUKxn88TJoOvDbrxnPf1P33rIuyKA//aubwc+32rD69IzL1Wmr2BzylRnHMBzepzX\nkcVa9LPC/IYMKm/8H76OCxKLDRKvNLxDM3Ptnq5HmLCez85WL4J8NoGuiVrExd8+zn2cnMywdnEN\nHfAsVQZt0GwP7BdmEI5fb9ZvZ/cbkw0isCJIYkqO1SZmr9cB47Py4+HOqj2+vNsWlUFXCTss6V1M\nS5XpojxqWjqsDMBzsj0hKoAMn9URY4gs4fwi+xFZEQJrCFkuJQt/z7V7kFPqmGv3BBLdKZ/pYX7E\nlrcvGPgdGEtSTq2EH0+h5VyqzHXxWzt6itcheuFvJ/x4s2ZPtryLbXH4OLxYO8/1RkLWomXBsNX8\nPvyIrl/vlldsy7vGY90eXZIuYr//qdaH4jaBmZsTn/Ks49eL9SuybnHPoSt/4KtNu/bdHiKwszVI\n+34+T55rZ/QFPyGlNMZq21TBzWaR6SHpJva1DlSV3wQ/8dmg3iwAtPVUJSeiChWwik+WYTKICKF0\na4PVK3Wz1Wif18P1492yerbmviA923DgnP7vWNdvmantY+bMK17bgAQIbVLh9VR4RTZMgQ27jaDn\n8xuHDSK4YYhtxzk4xdZPrHYNxml5pv+71fVVpmdSjcJ+NiyPbfELbU/xBoCRUM/B6806iezAzq/Z\nvp5Ukb2x5R8Tn3VgiOaNLf/IfTyP4PJ6t15i69QXJGKwm/6vCx9XXb9l+vPjer3ljRh5bTOIPUab\npPVovajPdKcEhy2yfuOUc+0eKHPsC5XHhfH+/Njq57Ph+fzdrq1dX+C9C0pjrDZCHKf/c9b1WybK\nNgApFFre+KwobxYA6rPv4mvM2jbHl9yPiFgHqAYXRrmOj5PBhdFAPyo8iBBbp74Qxe7GKWADIeQV\n02O/9yFd8Ez/d6vrt8zv+WyJPXSQZni8WV6R5UFUeGCm1SukHQBYl130fayfFC/ecAJP3NbvpAYW\neGO1ccVp25Rgrsm87dQNpym4uui6Tcs/pt+yT6N7oWongXOr67fMDe7jYhfaQv4mgHVxmwHA2xvp\n2fji8v9hiKxogZ1c7zwpT6SosrbvJr5OtvIKKM+khyBia+4LXe3KLARH9Kn5LDBP/3erq89s5S4T\naBuABAgt4LyoTVB4wgYst3yEtADEL7KsdallUfWwxdUL8/mtomu11QqP4PLU9Su2Rl+QhAPP9H+3\nun7LgthmR+wx2npjMG4TmFm8+nBqRBYARqYexEyrd/mRJKx2jUw9yHRcGJ+VnyyQxasPu7fJEauV\ng2KOPEYIOaiP6O9F9/T/jzHW9VVGCBnTlx8AIeSQJbPArc0uEuDRshPWavNhDGCIFlkegTHEKy0L\nn8y0erls5fVuRXq2YcZrJZ1wTv93qxu07ElYcDvOjtiFNputAhC/U0AYXoI68AbT8jdxiazVa53u\nucx0nIHI6ZMDOb4Y5XTP5WX7WQfSWEU0DLHND57xbo8jVpv0XRjaVMFCMx+3GakldqHNZxcQhtCy\nwurNzqo9UIZOeLcnUGT9CqzB1f7fOh4T9px0a/tewmu2lUdwWb1b0WJb2vBTzzphsQZniKWe2H9D\nl2ru6wckieUtUhxIksgCwD3XPtDxPM6FP7zObbUVAFdsWeRnysL0uT9jqifzaiVAAoQ2DFjDBjze\nrChECQKPCCVtVSVewWd9r6I+W9YfTNHbxEtWL7GHDghpI+kLfy+Tm3MsYvlyihRZL242i1hU6r4F\nliM5fZnBXJX7GMO+OsNt/Uyr1zOcwBIiYKnjGUJw6QtdbTHGapMep5X4J3ah1RaViWbhbz+YvZbs\njmds60QlsjweLACcZIgpA/5ElbUdVvF9ef3PgWbRM5Ybpdi6YfQFmYUgYSH238+l2gameqypXaLD\nBmZaPjYQBKITWevt+P1T73E+X7O0/AgT1vMYtrKEFFhCCSI+c7cfUN6+kPZYrUoJFpt5poekm9iF\nltLkhg26YnDN7skVXt5sFCLrJE7Fdqe4RSWuTrid32orq+C6nk/A4Jfj9TX1BaExfDl5YVUSe+gg\nzYjcSNEJFpH1PIcPYeXNmezLNbjqGzZ5hRZueoQTvEIJXiGCoCGEuJApXukido9W25wxWlhu4+y8\nlMzOv+c6j5eIBhFZFo/vpcFTzCK70Mx3PHixHs/ahuHh/nzoZ451wvZs/YQQrH2BxatNe/hA4p/Y\nhbYhwCs0CPu2SzXFPIOGDIKKrGvbunjdXt3mWi+IsLLA0/7I/B2uPwpePyxhhxGs11t1iX9LJFZi\nF9pWK9p4oV9vFgDozTuF2BCWyFrjn1vqw111whZXJ7zOa9jq5YGHJba8QiyqL6QFVSVYbOSZHpJu\nZIzWB0G82TBF1g0eYfU7ctzLGKc1bHGK63rFb93itizpX054xWu9cmtZUr1YcmplPu3qI3ahLeTn\n4H9HCfG4xdqUjT/wPD7ILapokf1V7xtMAisqJceuHTfxNQvur3rf6Cqfa5aEi62owS+WviCRGCTg\nd5N41ghreURu1GygTAM3EfYjsm6pWgvNPJotdwGNIu+R5RwLzTwyDl3RK27rhNvnGeSOY/n6q/Y+\nipyWK7EjdqGtNwYiO1fQUV918oOu5X6/wH5F1gnDUxyr7ewq400sFxWT8zrv3Qu3O3rffsU2TLz6\ngkRiJnahFYWIjAMvbySMhbTDElkrLOLqW0Q5j3Gzw01snd630+cUllcbFJnmtfaIXWhz2fQkizcH\nzjqWif7i8oqs3Yh+JX8NgLuwhTFazNKmVfgNWwH3gbskiO1Uqw9k8DeO5SLCByyOQ5QhNZUS1Jo5\npoekGym0HDQHXxXanpMI+BFZO35Nph1FllVcg36pvETXENw3TUIL2P9w+CWM/dKU8q+FtynphBAy\nauzLpf+1bvHNVNdvmanOMZvXDhFCKCFkhhByTN+g0ZHYsw5YF5UJitftmpcXMtXqQ++lP8Hi6N93\nlUVxG8orsovNPD5Ruxvf6vll5+se4srrkdjVL+aa9jbp5+7Nd2ci7Jl/J54f/Neu1xea+a40MKds\nBK/punb4zUKoT3wKxd1f4T5OwsVRSun9AEAIeQXA0wD2+ajrq0zfeHEUgHlTRoNzlFLvkXyd2D3a\ntQqvN2uHm8h2vebiVYq+7fNqz827tcPufYoMITgR5EfS84dbxmldIYSMAZg2nlNKZ2EveK51/Zbp\nz49TSsdFvJ/YhVZRWnGb4ImR0qPmZ7rK/HizIkIGduJjjXvO6onxXgLrRr2R9Xy44XQOq/Av2xqi\n2DrhR1Dt+oJEKKMAZi2vTeviyFPXb5kXZULIXj3kcMgtrAEkIHRQKkwh6MLfUS0tV731+WhOZIJH\nZK0cVSqAjci6iauXcLIcU8h3/3ga57SGFRYbefTmG/he8bWV15p55llmbpMarPiZNeYUWqje+jyq\njBs5rgYoJTx9Y4N+G24w7sMz5JnF5FbXb5kX47oHDELINICXANzvVNnzk9ODvHuh7WE+Zj6BQ/29\n5ueU0ufc2k/C5oyso8Q9Fz6Ope3PMtUV4c0GEdnFRh6PNO/Akdzry685CawfcXXD3J5VdGvNnK3Y\n/qf27fiH0soW3nZiaxevdYI3Xssbq+XpC2uQG5TSB+wKCCH7AexyOfYYpfQ4tFt6q5foJIxudf2W\nuWLWQErpKULIGCGk7KSNLN8w5oA0IeQggAql9DndlX4JgKvQ0hByU60EiYeZZ4IR60Lagga6RCfd\nG7fkJdPltRNZN4FtNfgXZM/m247nMAuunXebV/PL3q0Bq9iG7dXaYe0LdgTd5mY1rnnA4dlWYCN6\nlNJTPHUJIfBT5maYHlp42tBF03GODqjrZeQMSJcB/IXhwVJKZ62GrBV4vVnbNnx6s3bxWKvIOsVW\nW43M8sMPbsfbnc8pdtvx3OY9s8ZrRf2AeeXU+kUOiDljFTv9zvq4+bkRF3Wr67fMgwqAw6bj9sDD\nofTyaB2DxTaq/wCAih46mIUWZniOUlpxO0FvcRLAJg8zksGiw+aMQRAZMrCK1DezZ21F1oqTsKpN\nNsFVct2erNGm2ct18m6LuSa+mXWeDMITs2XByavlCR+E0RckXTym3yUbYcvHTGWHABwDMM5Q11eZ\n7mgaGQqHoIc1KKWzhJCKHgYBtFCIuc0uvISWJ1g8qhtqGPIKgJOwicfoBu4HgJGRMk6c/iIAYOvI\nT9BXuoazE1qYt9xfwe5t38XpV78EAFCUOu6+8ymcO/8oqlVtAO220a9ievZe3Jx+t2bw8DFkc3OY\nvKy10dN/FsrwP2P2zT/Xzp1dQHnXOG5OfArt+i1oQ0Fm+zNQZ+8DnXuHdp7hHwJKC+rbv49eqqDZ\n/yaa606i99znQHPzUHPzeHvzP2Ho8h8i29Q8muu3Po++2ftQWtiJDVBwbejnIDSD4WltAPPtnrew\n2HsBd0++HwBQyy7gteGf4cEb70Ne1QTyxIZ/wW3zd2JdbSMA4FTfWZTUAu5aGoVKCd7IX8Wl3BQ+\ntHgvVEpwQ1nCD4pvYN/SfVD0vdf+LnsGe9rbsJX2o58W8P/iHIZRxHtwC1SV4GUyjQt0CR/L3ApQ\ngotqFc/hKr5Q3AkFBCoo/mr+IvaVRrCtoP0IfGvpbezIFvG7eS2c9S/1GVxTG9hXGgEAVFpVPF+b\nxON920EIUIeKv65N4JH8ZmxS9DbUC7hD6cf9ZB3QBn6mXMMCmvj32AI0AAoVz+Rew6dadwFNoKE0\n8A+lM/hI7U6Udc/veP+vcFtjBKMN7byv9v0Wbah45+JuAMDlwiSu9b+OB6fep33GmSpODp3Au6cf\nREH3Ss+M/Aib5m/HcFX7cbdep3r/61gaeB0bLv8hAKCVn8XU5u+j58LHl8MFizueQc/EfwSB9pmr\nm/8RaPVDva5dW7LuF1D630D74scxAxWZ4lUMbD+C2d9+HlTVtp8p7/4yFq/8EeYXdwAARrYeQaO2\nCTPXP6CVb/gpir0VXLrwKACgt/cCdm4/gtNnDwJUAYiKe+96Eq9O7MPNha0AgHtGj2BuYQfemnxv\nx/cJmLB+Bbmh1F84Kdg56SloAghYPE1K6T6OukHLnrSxjcXzXYZQSp0LNe/0AKX0IdNrMwA+ZON2\n7wFwmFK6y/QaBbDLzau9Y/c2+ssfO19AlmmGXlkHbrdobgNh1tvC3sqfLk9YcLqd5Akb+PVmnQa/\nzNSaOfwZvQP/nWiDYWZP1u4L4+i91hmChAX7C2D1dK0xXLNn+2f0Dvxd/nRHuXVig51Xa43X2sVq\nnQbF7LxaN4/WXGbuC26ZB14xWq+1ab1itCz7hvVsnjjpNDjFSmF0C93yf36eqe75T/7nwOdbbXh9\ni3gD0lYcg8MScdiJrBk3kVWbmW6RrSsrDxYc6lvbtZ7bGsaw2s0Sr7USNFbrZ4DTLU4rl02UAB5C\nyxmQrgCYNZ7rfyteMdpiPlwtFjngUBt5yddxQebae3mzbiL7T7jsKbId2IglaSieD7c2rEJuHSwz\n7PsnXO6yn4Wot+Qx8NsXJGsTlvQunoD0PgB/QQg5By026zQveRmVBov7RDVZAQCI7rnwhg3sYA0b\n8GAVqWKrAEC7NXUVWRtx5cFcn+bVzjb1sILazHSEElqNzHIood7Ioi/nPF3XK+XLCmu6V5BULyJw\nU1HJ6sdTaDkD0hUAT/AY0Gj281SPlcLUu9EadB4dFw2vN2um3sji32Vuwa/acx0ia+vF6lgFVmmw\nrZmh5lfi/KShrIit0b6D2Jp5nzqC0xltWqt1UoOX2LJMZOCZwOCUfWB+XVRfYNlDTJJ+Yp+CGyc8\nA2F+YQ0biPRmvQa+tErOXqyTwJpfN4ur9XWjrQ7v1kZszV6t3ftxWgVMEgOUMKf7SbqJfd5Jutaj\nPe1diQGWwZkg3qzBSdq58MnyF8XixRrCqDRIh2gaz62vu5WZn3d4yJa4rYHxQ3CSzrjOVONdmJz1\nh8tv/NyjnPqvAAAgAElEQVTaF0T9MFuJMjQmCY/YhTabYZs6mQRafeeExGdF4+TNnqmv/IjZeSNO\nXqxZLJWG98N6jLm9jgEzF7F9XZ3vst9tYMz6w8MyKCZyqnOr7xxz3TAzDxKzcanEldiFtlofitsE\nZkqXPxpKu7xhA1bv7tHCrQDsB77sRNZOYA0yje7HyvH2gmsbgnBIGftEZhvTe4oKrx/OsPqCZHUS\nu9CuZvzelrJMUDDgis3aiJxVFA3RtBNVqxdrV+7Url0YoeMHwLRYvZNXy5tXGzTubdtmjHcukvQS\n+2CYojQBpCPI3i5cD9xG0NtXFm/WENmras0xLgvYiyxgL5wGmQbQzneWqflOD9d4TWkQqHm6/Lcr\nIwErg2NX1eAj7zzLKAZFRF9IFSrYJ7BIuoj9kysVpr0r+UT06ki1LS8IbS8oTt6swZHGla7XeEXW\n7LXmFuy9XauXaz7e0bO1fGmPNK5EMpfe7ofOz51H0vqCJNnELrRLteG4TejCaQS5Z+KTws/ldXvL\nMu3Uic8Xdmr/OHgidiJrDR0oHa9RZBoUuQVqeX2lrtGOuV0mW4s7Op6zhg/iIoy+IFm9xC60HBtJ\nxg5R7UfBRcbtgk4pNXuFeXR+tlZvVvtf+2u+9c/YiGt+niJTB/I39b/zdLnMKtTW9ry8WrWZQcFn\nVwzyQxQEp74gkdgRu9DGRRoX+xDpzZnF1uyBZiwCm6lDf1AUbraRravI1Kn+MI63F1ujTR6Chg+s\nP1RhDIg5EVYurST9xC60vaVrcZvAzJVtRwMdLzKP0ys+qzYzeGrhQlc6l503a/xvHdAyBDZTp8jW\nVSj6I3+ztSy4hodrFtvl4228Wif+av6Ca3kcuN2pLO78eoSWSNJO7EJba7ju0psoyvrCzmnho0X+\n+PeyN6uLbFYX1txcE5laG/m5BjK1dpfgmsXWzat1Ch/Y2ZrkOG3h2gfjNiFaKGFayY13QaK1Quzp\nXe2298LFSaGwxL4tepClEUUxmrW/bbYLG9hheLGZWhuZehukpi3Una21kKln0S7ot/mF4F+u0WwJ\nsF/2IJFkl7ZCzsmSsBK70ErE4Te+2RWb1b1ZAMsim5ntXJNChZb93C5mkK2r0G6OCNoWZ9PIu5U4\nI1fwskdf/3ovVpZoHXfaadatboCyMWh7IZYBPAjgCWN9bR7bgAQIbbEwA2BD3GYwMTPyz3GbwMXR\nqr/4tzkmC2DZk0WtDlRrQKkIBbrYWrzaTIPCEFyloU1cMDAmLjjamqK7zuqmF+M2YS1w1NhJW9+D\n8Gk4r3HtVpe7TN+44AFje3R9q65jWNkDkce2+Lu2mqI0mWyDZ6/K4ARNXRpROo93GggzkzHdDxve\nLIBlkaXVmia2y/W1ciNWu/w6Z7aB2dawJy6IGJTM1NPhHKQV3Ztcns2ke4t7eOv6LYO22ax5be1X\nAIwSQso8thnELrSNZnpSYvpn3hnp+YJusf1vC+u4jzGHzNtFk+AVtQJSKgKloql+xvT/ipDzhgvM\ntjqtUZsk8tP3x21CtNDupTGdHgA2EEJeMT32ezVvwyi69xyc1kWOp66vMn3Dg4dMrz8AYFYXVR7b\nACQgdCCJDmPdAe1/Z69WE0wFWWhCmjVCB+sGNc9WF11a1LqPWlDQMg2IGeEC1SK2TmGDqGHdaUHi\nmxsCdsHluX10q+u3DJb9Dg9gZRsv7lvb2IU2n1sEMBC3GUwslMUs/B0V/6PBtvGlIYh2wtsuZpCp\nZ7G8FExRq2yIbLuYQaugoF0gaBeAdt7fTD9WW5NCY90v4zYhleje7S6XKscopceh3Zpbcz+dBM6t\nrt8yq83PUkqf4znOTOxCm1HSkyTTKF1lrrsuuxh5ilc23+6Ib0607D03s2fb+ToAEABU/6sgX1fR\nGMwjU8isxGuhebrtYgaqR2oXawhholVLyyJuAIB26bKwttZSxoExuMRABTbiZd2Z26suIQR+yoz/\n9UGwii7+fmwDkIAYbbUe7QBTENZffci7UoL4RM9GzzpOQtguaCGExkAWakFBu5hBYzCPdkH72xzM\nLYcMzN6sNVzgib6XmJetSds/rHTl4bhNWNVYRUtPpzpufq5nBrjW9VumPx8DMG2ILCFkL8txdsTu\n0a4lBnI1odNwmSioQF1bB5Y0FMc4badAEmQaVB8Y0zxbFBRk6+qyB2sILNA5gAZo4m1ur/N/LU5r\nXZeW2EQcCvmW7Vsy74gbJnY74Uoi5TFCyEGs5Ko+Zio7BC3dapyhLneZLp4n9f+NuhUAzzG02UXs\nQpvJNABEt/CHQVlZ4l9Ypkfc7aIfevON5emnxVxzeVpqId/qWu9AybVxse18S2oW3LZp4e4VD9cq\ntoD5BsjwYrVjVoTbLLJtG4F14qKq7R0XNOMgsoW/bcJIQ9mFSM4dB4Sybz8vCt1zNLzH45ayfRx1\nucv0gTDHN+zWph2xhw6K+RnvSgmh5pCkLtLzCSoUZqFimbBgFkXr/+080R4FzWttDBA9nLASKrCK\nrEHbJvOgy5strHi1zzXY499mvFLgBnPhbP7p1BckEjtiF9rF6kjcJjDTW/m08Da9hCBILu0X+/UN\nD3VBMwTOELyVvyvHWMW2Q3DzBI1+0iWwZk/WLLgraV7u3qySa+MLxZ0dr5nDBkmLzwLh9AXJ6iV2\noQ0Tv6O5zreAyVqk3CxAdvFMxcZeJ7EFNIG0CmY7DzT7Vl4zP5p9nccYbWjtwvY8dt6sYStL2CCq\n+Kw3yeoLkmQTe4w2XcSfcG+O0zphpHmpZnv1QTEnzLm0RszWK4PAWu4lsnYoOU1c1Yg+W7vJCuuy\niz5air8vSNJD7B5tmhb+XhzlW+zZ7gvMMitJ1IDOl2vnl4XMjJ1Xaw0jtPPWwazuh4G5rlu4wMmb\nBYCvqG8u/88aNrCGVayfW1jxWYC/L6Qe2rl/nNtD0k3sQltr8M/Hj4vi1Q+H0m6QOK1T+CCbb2Nv\nfhOAFa/RLHBWsTX/r+a7BdftYT3GTsDtRNb8I/DHyhbXzwCINmzgNMBpvB5WX5CsTmIX2naKFivN\nVDd5fgGjwE1wzGK7LWOTs+sgtk7erZ0H61Rmbcd6HjNmkc3m27iVlLrs5xkEiyqtyyBT3RTp+STp\nJnahldhjFQ5Wr9aJjhCCjdgC3d6t9fbfTnSNenZCbW3fOK9dOMML64+Ln2wMcfHZbtxyaMuK/x/h\nQY9v6KCSnh1K1jKxC22pMO1dKSFUN7/AfQxrnJY3nsji1X67/VbHSL6b2Fq9W6twOj0M7AS2I1xg\nE5cFVvJ+v91+S2hKV5jxWcBfX5CsXTyFVp9TfJAQskf/y7SbIiHkMEu9thrfL7Kbp2HnoWSq3nHE\nMOH1arcTbeYbi9gC3bf4TqLqVtYhsDbnsIYMDHblnPuBlzcrOmzAEgaKuy9EDVG1TBSWh6QbFo/2\nKKX0SX1hhXFoWza4oq94w7TYb6MZbIUrr1srkeRn3uVaLjpO6yUgVgGyDoy9W1lZsMdVbG28W7u4\nqpPomo/rwDLw5SSyhXwLD5q2MwrDmxWxBq35+nr1BVbW0spdaxlXmfKzZYPu8U6jewXyVQOvoIoM\nH1i9OTexVZROQbSKrZvgAp2iaxZSp9ed2rLGZN0mJlhFNkxvVlR8ViLxwssf5N6yAcAet3UZreRz\n4S7EIdJjaKw/KawtVoIIyQlc75oxZhW5roEpQyQdYqp2nq7TcV1ibnN+w74TuM4tsnaE5c1aMfeF\n1byYjEQMXkLLtVisHjLwXMmmwwDF/TYxzlFV6xeoXbghtP2wvdpJaO3bia2rd2tgFk+vhwkngbWG\nC8x2zfkQqzC8Wda7Fda+ECTjQLJ68JqCy7xlg75+47Tb3uamuvuhx3CHh4dw4vQXAQBbR36CvtI1\nnJ3YCwAo91ewe9t3cfrVLwEAFKWOu+98CufOP4pqdTMA4LbRr2J69l7cnH63ZtzwMWRzc5i8rLXR\n038WyvA/Y/bNP9fOnV1Aedc4bk58Cu36LZpB249Anb0PdO4d2nmGfwgoLahv/z56qYJm/5torjuJ\n3olPQC1MQc3NA1u/g+zFjyGrby55/dbn0Td7H0YWdqINBdeGfg5CMxie1pz/ub4KLpSu4u7J9wMA\natkFvDb8Mzx4433I6/lSJzb8C26bvxPratoi2Kf6zqKkFnDX0ihUSvBG/iou5abwocV7AQCTpIof\nFN/AvqX7kIOCNiX4u+wZ7Glvw1b0YwNK+Cr9LYZRxHtyt0BVCV5Wp3GBLuFjmVuBEnCxXcNzjav4\nQnEnlCKBCoov187jT7Kbl/Nwv7X0NnZki/jdvNYV/qU+g2tqA/tK2oJAlVYVz9cm8cX+7QCAOlT8\ndW0Cj+Q3Y5NSBAjFt9oXcYfSj/vJOigKxU/oJBbQxL/HFhBC0d/K429yv8KnWncBABpKA/+AM/hI\n7U6U1SIUQvGDzK9xW2MEow3tvK/2/RZtqHjn4m4AwFTpMi4oFTw49T7tM85UcXLoBO669n4U2trA\n4JmRH2HT/O0Y1vNgrdep3v86lgZex4bLfwgAaOVnMbX5+9j01n8AaWu5vos7nkHPxb2guXkAgLr5\nH4FWP9Tr2rUl634Bpf8NtC9+HDNQkSlexcD2I5j97edB9cHf8u4vY/HKH2F+cQcAYGTrETRqmzBz\n/QNa+YafothbwaULjwIAensvYOf2Izh99iBAFYCouPeuJ/HqxD7cXNgKALhn9AjmFnbgrcn3dnyf\ngAlI4oVQ6jxnWw8RPG3sX66/NkMp7ZrOpa8+bhbhw9A2NDtu2eSsgzt2b6O//LH7HiZzqvt2N3MO\nd7MGU6r7Yttu69JOtVZ26e2t/CkWR/9+5bxt++OcXnfa2sZuMfC5ZvcavQuW7cfttiM3r4PwmeY9\n+Bu80VXHunYtEM4W33axWKt3bXjfn2neg6/lzgBgCxnYebM8YQNej9b6urkvBMmh9Qpticij7dk8\ncTLoZomljVvprke/yFT3zP/9xcDnW224erSmfXUA2G8nAd2LNW1cZpQdZtkfKJOpA+BcgDsmWj2X\nmOoNZpZsxZZnH7HBXNVWbM305hpdYmtedOYSmUcxqwmZsUg4sCJ2ZsE1i2IQ0XUa6LJbXcwck71E\nNO+QZZpt1CJrh9EXgsRnZcbB2oElOeoxI48WwF50byfxMXNlQkhZ3+IB+nGjbo0X88lOTjB/keoj\nP+woE5XOxTpYYycwdt6eIVbHMxeXX7NLmXLaKsaIp1rjqn7qWmOxTvYcz1y0Fdkg6/GKwu46W/uC\nHTI+GwyeHH63ugHKxggh+/XXj5q1jBByiBBCCSEzhJBjXjrnuUwiz3YS+muzAJ7UH55EsfD3kFLz\nDB+w0Hv+0x2hAz8E9Wr7co2uEIKTZ/vxxd9Zvh0HOre/MTCLoF1IAfC3vYyTiDvlyH6ufTe+hc4t\nvJMQMnBCRF/wIsoc8YRy1AhbEkJegZbD36U5DHW5y3TBfcC4K9cdzWNY2Sr9HKWUeVFieSkRjufB\n+8XlSUFi9WwzpDv+Xsw1HcXO8D6dRNILr+PtztubbzB7sqwiKxK5QWM88OTwu9X1WwYttfUJ02le\nATBKGGfGWol94W9iIwa8DCreA2JBGMouYKrVB2qTiuYUj3XDyau12yXXKVbL4tk2oS6LmHWxcEP0\nrB6ugV+xtcNJ2M0C28TKBQwqsmF7swBAleaayp8llGt67QbdOzQYZxmvseCYw2+Tp++W7++rTB+f\nesj0+gMAZk1ZVWU9AWAWwEMA/tIt4yp2oe0pTgLYHLcZTCzteIarvh8Rtm3Hp9ge7fn1yusOOzOY\nRdBJdP3gumC3jQd7tOfXjvFYnnxZXpH1y9KOZ1z3bg6abZBybgjIOuDJ4Xer67cMlmypA+gcnxo3\nhJUQMg3gJQD3w4HYhbZa55oTERos248XLz+M2pbuVZvC9moBfrEFgPfN34sfFFfSu5y8W4Mgosu6\nPoFTVsFHGrvw49zZrtedRFZUyMDNm3Utu/phYOt3hNiwltBz6He5VDmmr6vCnMPvUddvmdXmZ82Z\nVWbvVfd+xwghZSevNnahVVVxXpQbQQfEhrILqBkTHDhwE+EwxRYAhqm9z+UluIDYnWe9UrZ6cw2s\nr/Z1vc4rslF5swCAWnp2b04SHCGECmxEz2F6v2NdPT2Vu8z4Xx8Eq+jib7zWNb9AP84xdCAHwwTh\n1zNywkk0nETGSZR6cw3HW3JjICqsLWK82nazLQqRFX3NDEQMrrJkHKzmRb+tgmqXw28MTLnV9Vum\nPx+DNk/guP58r15UgTYhy6i3B0DHPAIrsXu0pcIUgOG4zQDgHT4obD8C0b6RW7pXUM/2pd7Ty/8b\ngmY3mwzo9jq9dtplacOxno2IGra6xWN5RTYshrILoNuejfSca5TH9Jz8UwDG0J3Dfwza0q1edbnL\ndNE9qf9v1K0AeI5SOksIqeghBUALhZjb7CJ2oW213Wc/scKSeRA0fKDO7wbKv3G2wSVM4CeE4Iab\n2ALadN2tzSGcybzVUe4luMv1QvBy3SYfbG0O4ULRcaa2L5EN05tV53cjU/ifnvXsSOVAGI1+h1ue\nHH6Putxl+kCYY56sOZTAQuyhg2bLexApKbdIdOZ3Ah3v9gV2EgU3IXEbEOrLNbC74byBoNttu2i8\nztWXa+DO5kbH8ihF1gsjpcupL8jZYBI7YhfapOGZluOROxlGgrtfsVUI9d6lQRdBkcLL06ZXqEC0\nyHohJyhIwiD20EE+N4/uDItkotzyk8Bt+A0hGMLiFLMFulf8Otuj3YqbwwleuAmjOdwQRJTtxNWw\n1cDtBySIyAYJGZh/ZP32hVSGDSSBid2jVQj/PPogsHR0R682q60wFdSr9RNCMODxbqtK5/KSfbnG\n8sMPQT1ft3Mbtrp5sUCwgS+h3qreF8yIChvINQ5WH7Ff0lpDnDcbdgdVr35k+f8ki60hVGMLdznW\nCyK4PLCK+9jCXZ6TELxENkhclsebBTr7wlpAm4JLmR6SbmIPHUjs8cpEcEr9MhjMVZEh3gtA2Akg\nS4iBpz0vWGwNK1zAUs4CizcrwwZrl9iFNpupAQi25TgvLGledjm1pO/Nznb0xWac8Jqa61XOIraA\nfdwWAG4UJh3jt25E4ekCnaGOG4VJ2zpBvFhAjIja3b1Y+4JE4kbsoYNC7mbcJjCjDP+I+5ig3tS6\n7KLvUMKb/a+tnMcj9hkVhh1WW8y2GkQhsrwhAwNzXxCZ0sUa/kpKyqOEjdiFdrEmdlaYyDit9QvU\nrnyuqw7LUnkibl1ZxNYqTO+58W+7z2USuqiEl+V8Zlvt3osVEesXBPF27fqCGzJssLaJPXTAyqBS\n8NykkQdRuy4A3iEEIHgYAWCbQeYVTug6r0X8eEIMrG2ywppREIUnC7D9iMoJChIWYhdaQlQA4ndg\nFUVHrDYTwipQJljFFnDeUddgIFeDmuEXvDjCC6y2snixkYksR19YFd6sCmTE+TlrjthDBz3F67Gd\nm/ULYHgt2dGvO7clIIRg1BERSgCAMxt/xHQbHgeGXYZtZzb+yLV+lCLLQnb068K9WZk/u3qJ/dJW\n60PC2wyrw7YufNy1XJTYstbzGii7c/J9y/9bhS0unGww22qGZTAQECuyLNdRufgnTG1JJEACQgeq\nGq8JrLHasrKEGw3vHwUR8Vreek7hhKKDHVahY43n+oFV2K228gx2RS2yANBmXAR+VYQNJIGJXWjT\nRAZsO0CKFlsAgQTXCycx5BFgUZ6yaIHlqccqsmVlCTNMNdnhuQuTqV3pI3ah7SneAMC2LQhP5gHP\nzrisXu3Azq+he4a7f3hElGdfMkOsXh0OtghOlGGGK5teTIXIGgzs/Jp3m6vImyUUyNTl9Fq/xB6j\nZVmPNinUZ8aYB0B4vrg8osEzmLO9uok5xhkHhm3rsosYXHDbr28Fns8gDJE1rn99Zoz5GIlECq0O\ni/dRn30XAPbcyTDElqfu4MLo8v9mUYsTJzvMttrB+yMTpsgCK33BsV0Ob1ZmG9ij7wt2kBCyR//r\nuAKVW90AZWP663sJIYf17W24bQMSEDoIE57wQViwxGsNeMIDPGEHK1aR443pBjkXL7zpWDz1/Yqs\nZ7urKGQQM0eNnWYJIa8AeBrAPh91/Za9BGCnvkfYegBHAdzPcFwXsQttIX8TwLq4zQDgHavt2fji\n8v9eGzl2tMsptgC7gLrVn1xvtzNzN25iyCLCIrxkq61+8l2jFFlzXwgCrze7VgbCjB1ojee62O3h\nreu3TGenaQvx5Xo8thnELrQAX4Cddyour1frJraEtDqe84otgFC8W6O+gXEcFbCoelShBsPWsAUW\n4B/4ssPaF5bblt4sAGzQvTyDcUrpuGNte0YBzFpemyaEjFm3CXer67eMUnrKJLIAcADAEz5sA5AA\noa03BuM2gZnFqw8jP/B6x2s8YguE691aj9sw9SAWei5zHRsHg5kljEzdj2v9v/V1LA+8IusUMrDr\nC6tZZIlKka0zeyw3KKUPBDzlekF1/ZYBWN52fC+AY6adb3lsA5AAoY0CkV6tHWGKLeBfcDNQbT3d\nJBB0Kqyf40WJrCjW4iAYIWQ/ALcUE0PQptG9maCTwLnV9VsGYHnb8ScJIfsJIccopQ9x2gYgAUKb\nzVYBsItOVNiJbX7wjGP9sMUW4Bfcat952+N52hCFlzBabfXbjh1+QgVeImvtC6vZmxUJRwihAhvx\ncrg1d6xLCIHPslEAeymlT+ovfxuAkXnAYxuABAhtPruAKIRWRAZCacNPXcv9iC3AHrc1YBXchfKv\nPduwI4gI+/VU3WwN0m4YIgt09gU/IuvHm10rA2FAhxACWL6FP255Pk0pnXWr67cMWhzWPOd+FMCs\n7uHCzTY7YhfapRrbnHEzotemdcLq1c6dO4B1d/yV6zG8Ygv4824Bb8G95a2P4tqO/893u1HiZGuU\nAguwhwtY+oIkMI8RQg4COAVgDMBjprJDAI4BGGeoy11GKT1OCCnroQ4AeAjAhxjb7CJ2oY0SP16t\nnwXC/YotwO/dAvGGBMIgqNCHLbId54rIm40dCijsg2FiTqndihu348ctZfs46vote870dNxS5nic\nHZ5Caxp1M5R73JL2YK47BuABaIHiBwE8Ybjazu23keSFv4EVsVVytm/bFj9iCwQTXKDTy23lgqcx\nRQXJz6VKYJXcbKQiu5bCBqsRFo+WaQaEPgXtASPYrSfwHoP7CKO+qMxmTrP9hw/8xmqHlBow6r2Q\niBnjC+xXcP2KLaAJbmvbt2EkzyXR0zULa3Xrd3y3EyQv1m9mwehtf+P7nJK1h+vvq90MCABOMyBG\nsZLQCwCvABj1mgO8VNvAZmkCqJ7/tK/j/H6Zh7ILgUSkdGllcWpjrQDzI0q8zm+2lZWgn4/v66LU\n8Na5A9zHpTJkIBGCl0fLPANCH8F7yPTSA9BG6VzvtymNPmzg16ttNdf53tQxqHdrwOPlKs1+13LR\nmQdBxNvLVjNBZ3YFyY81wgWtZjKmjUvSgZfQcs2AsMRjD8BhJE4fydsPAMPDQzhx+r8AALaO/AR9\npWs4O7EXAFDur2D3tu/iX898AQCQyTTw4F1fwenKJ7CwtBEAsHP0bzE7dy9uTL0bALBp4zHkcnO4\neElrY2DgLLZs+j7Ovv649oazC7hz91cwef6zWNJXyd8yehg3px/A/Ky2XsTQxu9BUZq4fuWjAIC+\nwd+gvOHHqNc2YuK1LyGbm0Fp59cxV/kM1KbmsA/uOozqjd9DY+4eAEDvphdAaRZLb38YAFAo/xKF\ndadAz38GbShAfgrZ7c+iVfk00NbWE8iM/i3UyQ+ALtwGAFA2fQ9o9UO9/n7tc1v3C6zvfwPtix9H\niypoF66jtuUF9Ex8EkTNAQAWd34dhWsfRHZpK5T6EJTqJmTqG5Cf1t5bY90v0S5dRunKwwCAdukq\napteRG/l0wAIAIrF0a9jZPL9yFQ3AQCqm19AproF+RltxarG+pNoF26gdFV7b62eS6iP/BC9lT/V\n+oHSxNKOZ1C8/DAy+mdc3fI8sgu7kJu7FwBQH3oZNLuA4jVtIJc0+wE1h96JT2ptZKpY2v4sSm99\nFEpDE7XCjm9Cnb0Prbl3aJ/P8A8BpQX17d/X2hh4DcrQCbTPa3YgN4fsjmfQmvgk0NSCKOt3/Xcs\nTn7Y8zrdPP8Z7ZoUrmNgxzcxe24/Mu0ezAPYdvtTaDbKmHjtSwCA4S3PodUcxPSk5mcMrH8ZfQOn\ncWXiswCAfPEK7hr9Bl597XGoqhZrveeuJ3Hx0h9jXr/WO7YdQbW2CdcmP6C1ectP0ddbQWXiUWRA\nMNB3CXfvOIoTZx4HKAEIxXvueQqvTuzDzYWtWpujRzC3sANvTb4XwMr3CZhAUIhKkakFn9K9ViGU\nOq81QAjZC+CAPhvCeG0GwIfcknN1IZ22jNrZ8jv3lej/eHEjn9UmgqZ58Xi2rWYfsqYBpqDblfvx\nbu1w8nJJqwSajX5XWz842SpiXQIg+Cwv68CXtS+4ETRkEHQgrGfzxMmgU2IH+rfQB+//PFPdH/74\nPwc+32rDqwtwz4DQB8EqLCILAI0AAz5RM3vj33U8DzobSNQUTyNWaRWl3Mz9DkckD7OtTu/HD2Vl\nSbjIAt19wYm4RVaSDFy7gVVQ7WZnWBfKhebJHtef7/UyoNUq8drcQdCOyPNFWNBvWc2IEFuRc+rN\nIpWbv01Yu2EylF1AaWFUmLgC4gTW6fra9QWJxAmW9C6m2Rm6CJ8EOqanVQAwebZxEnR6rvFlDBJK\nCDJY5kSWqF3CFSRlTBSixNQOYXcJAtYukFkGEgNPoWWdnaEPhBFwUsjPwceqYx2ImJLLIra3bH7e\ntdxvRoIZs1AEFV1l4w+6XnMTOZEizCumdrbyIPSugEFkvfqCCJGVYYPVQwKm4HJrc2h4ia2qj+y7\nIcK7NQjs5ap8lzdMT9MTTlsNohZYA7e+sCo9WQpk6jLrwC+xd4l6Y0BIO6J+/d2+JFNvf4S5HZHL\n5n/qvTEAAAsRSURBVBnxRl5RUSc/KMyGsOGx1e/n4Qbv9XLqC6JEVnqzq4sEeLTJQ9SmjiK9WwOR\noYU0EdYi3CJ/EFelJysRQuxCm8suAWCfFeSGyOUT7cS2v3zSV1thCC7gLbpk8DdCzxcmVlvD3t0g\nqMBa+4JIkZXe7OpjVQmtaKxiO7D+FefKDIQluIC9MM14LKadJMrrX0YmZHEFxHmw5r4gPVmJF7F3\nEdGLyoj2BgaVlS/S5Qr/QiJ2uOVnCuXCIx3xzLC9RFasNpWVpeUpr2Eh+jO/XDnQ0TdEkVRvlqgU\npNZieki6id2jDYMwdmAIw2sJ08N1gkVsg8R+kyLmQLj7eIWxFFJSRVYSnNiFVlFaSPrC3wY9heuh\ntGsWBJGim/Fpbxxi6ddWO8IUWOMHd1JwX5Ai2w3npgOOdQOUjUFL8i9D28rmkGnPsEMADkJb3fAV\naGvCOG5yELvQlgpT8LPwtxdheLW37/rq8v8ishLsECm6Azu+GdScyAhqaxShGPNdjbkvSEKDadMB\nhrp+y14CsJNSOksIWQ/gKABjUY5zlFLmSQAJiNHyb87Iimgv4bU3/tzUdviDIEZc0a+IzJ7b710p\nIfDaav5swhZZu2tt7gvB25ferBWeTQfc6vot09lp8qCnEYDYPVpKw1UrkZ5ty2aKqvEFDMvDNbCK\nCYu3SxOwrgErXrZGMnhowe2H1K4v+DvHqhTZDbp3aDBubHHFAfOmA251/ZZRSk9ZwhQH0LmDTFlf\nNGsWWljhL902OYhdaKMgiu3JoxJcAzvhiXJQLWziEFaDqNK1UiWylEKpNVhr3xCwHi3PAihudf2W\nAeiI4R4zViXUMcdyp6GFGRzXJY1daHuLkwA2hX4eEWJ71x1PMZxH+xuV4JqxitO63f8NilJPtAAb\nNhu2xgmPwLL0BfdzpUhkBaJvCuC2YashaNPQBqHMOAmjW12/ZQCWF8t6khCynxByzNgEwey96tt4\njRFCyk5ebexCW2+KWeuAhaBie/nqH2Dbre6rNq2cS/sbh+Aa3Hj7DzC85Xkm71C0GPN6pIatUePX\ne+XpC93nXJsiCwAcIQSeTQcc6+pLtvopGwWwl1L6pP7ytwEc1l8vA3jaGEQzHZfc0EGrHa23FURs\nb968CwDfl8v8RY5adJfm2e2N81Yd4LM1KCJCA376gnbutSuyPJiEEID9pgPQNhmYdavrtwxa/HbI\nZNIotM1mK/pmB4dNx+2Bx7rbsQttHBidPey4bfd5V/6P09Ndi8Q9TVYKrC+YNh1gqMtdRik9Tggp\n66EOQBvw+pBeNksIqZjKdsFhI1oD180Zo+Cd9w7QE8eGvCuGBI/Y3py/DQP9b4Zkh/g2l+ZvQ09I\n9oomDFvDFFeevhC3yIrYnHGwMEJ/d/Mnmep+f+IpuTmjhdg9WpXGOyuMx7tt6ltWh2NH53MRwtsK\n0V7RiLA1Sq+VpS/ELbCS5BD7hIVGMxkrdw0qBc8vxtW3H3ItF4mRJG9+8DI9GZ29QeG1VcTnEwS3\nvsDSlyRri9g92qQRV/yWBScxWc3x3rhjqzxIcZU4EbvQJnU9WjvB3TD0clzmuOIkRiNDLycizcwJ\ns91mW9OAuS9IgZV4EbvQZjNVJFFoDcxfotrg6Rgt4adssjfpIlZO2Wd767o30CsFVsJI7F+/aj2+\njANezlc+l6r425uVz8ZtAjNpsNW49oNKAb9589G4zYkWlQLVGttD0kXsHm1aMYttEuO5EjGk5UdV\nkmxiF1pFaSItC3/39bxt+7rdlzEJ4lsqXYnbBGaSYCuPqDr1BYnEjtiFtlSYRhgLf4fBvaPfYq6b\nBPHdtfMbkZ4vCFHbGtRT5ekLEknsMdql2nDcJjDz87PBFns2x/jMj7B49bXHQ2tbNGHZGtZnHrQv\nSNYWsXu0HLtBxE67nQ+lXZYvvh9vWFXTE1/0Y2uc8dOw+kJiUVVQOdDlm9iFVsLGah+UyYCs+vco\nWbvELrS9pWtIS4z239zz5bhN4OK99/w/UEg6xCttn23a7JXES+wx2lrDusB5cnnj4h/FbQIXabI3\nTbYC6bNXEi+xC227nQ6PCwBm50fjNoGLNNmbJluB9NkriZfYQwcSiST5UEqhysEw33gKrWkXSGMF\n8nGnvXF46hoUCzMANnCaHQ937XDdrSJxpMneNNkKpM/eNCJKe/yWWdo/TCk94Mc2gM2jPWpsQqbv\n1f40gH0C6gIAVDXHYEIyWKiOYLDvYtxmMJMme9NkK5A+e1OKKO3xWwb99T0A9gM4YHqZS+tcY7SE\nkDFoW/ICWN7lcU/QumYazT6vKonh0rX3x20CF2myN022AumzN22I0h6/Zaa2y3qdWdNr3FrnNRg2\naj6BzrR+oiB1JRKJxA1R2uO3zGCPzRbn3FrnFTro2vNcRF1990hjB8l6z2akZDHS/30DgBtxW8FO\nmuxNk61Ayuy9I2gD83T6xWPNI6yDKUX9dtpgnFI67ljbHlHa47fMCBkctynisQ2At9BOA7Amujqd\nhLmu/qGPA1p8Iy07ZqbJViBd9qbJViBd9lpEzxeU0j8QZMt+aNtzO3GMUnoc4rTHV5k+2DXtMMDF\nYxsAb6Gt2DVg40rz1pVIJGsQDs9WiPYQQuCzbC+A9YQQ48e0rP9IHOe0DYCH0JqMAbCs8sctz6cp\npbNedSUSiYQVUdoToKwjf09P7xo3PXe0zQ6W9K7HCCEHsZIv9pip7BCAY9DDAB51neCN3cRJmmwF\n0mVvmmwF0mVvmmw1I0p7/JYZWQf79f8PAniOUlrxOs4KoZSyvmmJRCKR+CD2tQ4kEolktSOFNsUQ\nQkYJIQcJIXv0v0xLoRFCDodtmyRaCCHHGOr46i+S4IQeOgh7rYQYbR0D8AC0NI8HATyhx24igxBy\n0jQNsAzgaUqp65RnPTfwGI14awvea6uP+i5jHZwIG599AdD6w3NR9QX9eo4COOx1Tf30F4kgKKWh\nPgCcNP1fhjZHOHDdOG3Vy/abnu8BcC5iW8egCab5tRmPY8r6ca71EtAPDgLYa6p7MkzbRNhreX44\nBnup6P4iH+IeoYYOolgrQRSc5x8F8ITp+SsARiO+FfMz5dluOmHocPaDMoC/oLoHS7X0nfsjMXTF\nBt6+eCAFt+FyinyMhB2jTdNaCczn18XqIdNLDwCYpRGGOcA5DdBlOmEU8FzbBwBUCCF7TbHEqFfZ\n5u2LhwCcJ4Ts15Pan3CoFyfc00Yl4gh74e9Q1koICa7z084Y3AGw5QyLhHkaoMd0wijg+WxHod3m\nHqeUzurTR0/CfdqmaHj7wrju0RrL6B1Ht1DHDfe0UYk4wvZoQ1krISR8nV/3YJ6lEQ/WgG8a4BiA\nB0weV1n/PypPkeezrQCoGD8K+t/RiL1arr5ACDlIKX1SD3EchpZInzTkFPkYCdujTdNaCdzn12/H\nK1RbBCNSKN8URdfphBHA2w+sRO0dMtur94FTpjrjhJBdhJCxuEVMTpFPDqF6tNaOZicGxiCCV92w\n4bFVfz4GrRMf1593pCNFxGNGXiS0VCTrFMWPmSsTQsr6tEFEGfvk7AcVALPGc/1vhUaYOsfZF6ah\n3TG4thEWhJAx0zU9pPcFA2sfcOsvkhCJIo92DNqIrd1+PUehpZyMe9WNAlZb9S/eOcvhFUpplHHE\nVMHZD0ahxTvPQYvNHo5SaH3YuxcrHnAZWnxZ3pJLlpFrHUgkEknIyCm4EolEEjJSaCUSiSRkpNBK\nJBJJyEihlUgkkpCRQiuRSCQhI4VWIpFIQkYKrUQikYSMFFqJRCIJGSm0EolEEjL/CwxhSM1Rlh+s\nAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f4019511890>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAUYAAAEACAYAAADP1t+BAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztndtvXFeWn3+7LrxfihRJkZQluim3ux0nxlgtB2i4g3Zg\nCh60EQdtqD1IXvJkCUGCgTEDSOgAecjTQPoPRD/mIfBIsREBdiKImrgBGwbGkg0YceLx2GxdWiIp\nSlWkSIqXuuw8nH1KR8Wqc/au2zqLXB9QkKpq16mPp4qL+3bWUlprCIIgCE9IUAsIgiDEDQmMgiAI\nFUhgFARBqEACoyAIQgUSGAVBECqQwCgIglCBVWBUSl21aDOtlDqjlJox/2Ya1xMEQWg/Kmwfo1Jq\nBsA0gAtaaxV6IKVuaK1/Yf6fAfC+1vp3zZQVBEFoB6GBsdxIKR0WGJVSxwCc01qfCDyW01oPNUdT\nEAShfTRrjnEawErFY1kTMAVBEFjRrMA43KTjCIIgkJNq0nGyACoXW2oGS6XUKQCnAKCzq+MXhw6P\nAwCKyU1oVUSq0AcAKCXyKKbWkd55MiLPd+SQyg9A6aR3P/0IiVIHksWukGNsIL3j6WloFDpWKo6x\nikSxC8lSJwCgkHoMQCNV6DXH2EExuel5qBK0KqGQXkUqPwilE+VjJIvdSJQ6zDE2ACikCj2eV2Ib\npeQW0vlBz0MVUUg/QmonAwVlfrYVJAu9SJTS5hjrUDqJZLHb/GxbKCV2kM4PPHWMyvOTLPR5x9AJ\nFNKPdh1DJ/JI5fu9n00VUEyv7T5Gvh8J7X09Cuk1qFI64hw39jl5n8tqcz6n8jlu3eeU2hks32/4\nc6pxjGZ9Trf/X/aB1noUDXDiX3brh9mSVduvv9m5orX+80bej5pmBcZ5VAmEWuuvqjXWWs8CmAWA\nqX9yQP/lhz9tkkZr6dg8iJ3uJWoNazj5cnIFePm+98K1W40e42G2hM//16RV257JmyONvh81dQ+l\nzfacDLA7ACqlpgHMWR3H9AY4EPxrzQFOvpxcAX6+ghuhgVEpdUwpdcb8/5zZvuNzDsA7gfvv+vsY\nAZwE8K6NgD904MBA7iVqBSc4+XJyBfj5Cm6EDqVNT/ArAOerPPe7Gm0By96iIAhCHCG/JLCY3KJW\nsGYt8y21ghOcfDm5Avx8BTeatfhSNzqRp1awZrt7kVrBCU6+nFwBfr6NUoTGammbWqNtkPcY/e0I\nHBhZeJ1awQlOvpxcAX6+ghvkgVEQBCFukAfGkipQK1iz3XWfWsEJTr6cXAF+voIb5IGxmF6jVrAm\nO/4ptYITnHw5uQL8fAU3yAMjp42yE7d4ZVHj5MvJFeDnK7hBvirNivCUlPGDky8nV4Cfb4MUNZAt\n7Z+fWQKjCyo6d2Ws4OTLyRXg5xsjzCXDJ+FdEHIMwKzWujJtYWRbk9bwuGmaAXBJaz3fDEfywJjv\nyFErWLMwdZFawQlOvpxcAX6+MeNiINv/dQDvA6g1NxHWdkZrXb4qTyl1AcDpZgiSzzEmGe1jHF58\njVrBCU6+nFwBfr5xwfTysv590/ubqbPt6VbVliIPjH4uOQ50bo1RKzjByZeTK8DPN0a4ZPuPansO\nwB+VUqdMjtezzZLkE5UEQSCjAIUVkyDYghEz7PWZNTlYAbds/6FttdazpsfoD5/nsDuQ1gV5YCww\n2sf4YOIatYITnHw5uQL8fNvMA6318RrPuWT7D22rlDpj5hjPmx7jVQBH6/DdBflQWpnU7hzo3Byn\nVnCCky8nV4Cfb4xwyfZfs63J+/pV4LFZAJeaVYCPPDD69Sk40L/yIrWCE5x8ObkC/HzjQlS2f4fK\nAFl423dCj18v5ENpQRD2He+aygD+3sRgtv9z8IbEs2FtTa9x2gyhAW/I/UGzBMkDYzG5Sa1gzaOh\nb6gVnODky8kV4OcbJ8Ky/btUBtBaX2qVI3lg1KpIrWANp83oAC9fTq4AP99GKeokHhb7qDXaBvkc\no19XmAMHln5NreAEJ19OrgA/X8EN8sAoCIIQN8gDY4lRzZetnnvUCk5w8uXkCvDzFdwgD4zF1Dq1\ngjW50c+pFZzg5MvJFeDnK7hBHhglUW3r4OTLyRWw980Velts0h4KSCBb7LO67QXIA6Mg7FX8oLhX\nguN+QgKjAzrBp3AXwMuXkysQ7VsZDCU48oI8MHLaD7Z45ENqBSc4+XJyBcJ9awVBCY58IA+MqfwA\ntYI1IwtV82nGFk6+nFyB2r5RwU+CIw/IA6PSSWoFa9LbLqnk6OHky8kVqO07lNoIfV3U80I8IL8k\nUBD2GkOpjao9Q85BsagTWCn2UGu0DfLAmE8/olawZnnyCrWCE5x8ObkC0b6VwZFzUGw2TawSaH0c\nV8gDY6LUQa1gTffGFNY6+GRV4eTLyRWw8/WDowTFXTSrSqDLcZwgn2PklKi2b/Xn1ApOcPLl5ArY\n+0pQfJpmVQl0OU49kAdGQRD2Fc2qEuhyHGcih9KO8wHHAPhFcDIALmmt58OOzylR7erw19QKTnDy\n5eQK8PNtlIJOuGw1akeVwJZuY7CZY3QZx8+Yql0w7S/gSWnDqnBKVFtkNizi5MvJFeDn22baUSXQ\n5TjOhA6l6xjHn/YL2djCKVHt8P1fUSs4wcmXkyvAzzdGNKVKoONxnInqMdYcx9cQOAfgj0qps+b+\n2SptBEHYp5giVuX71aoEAshqrVfC2kYdp1GiAqNT11RrPWt6jP7weQ67AytMZa9TADAyNoSJm+8A\n8AoM5Tty5bTxWz33kBv9vJziSScKWDzyIUYWZspXHixPXkH3xlR5lXB1+GsUUxvlv+ibvXeweuA6\nxm//FgBQSm5h6fBljN57A6mdQQDA/UOfoPfR8+hdew4AsDLyJbQqYGj5lwCAx303sZb5BslCHyZu\nvoNiegP3D32MsbtvIpn35l2WnrmM/pWX0LP+LAAgN/oFlE4h8+AVAMBG/w/YGPgeY3d/AwAodKxi\nefIKDt55CwmzMr945CMMPjyO7o3DAIDs2GdIFnoxmH0ZALA++B02e29h9N4bAIB8ZxYPJuYwfvtt\nqJL3US5MXcTQ8qvoejyJZKEPHZsHkd4ZwkDuJQDAWuZbbHcvYmThdQDAdtd9ZMc/9c6xVoDSWJi6\niOHF19C5NQbAKy7fuTleLhnais9JqxJUKd2Uz+ngn94CgJZ+TlqVyt/bRj8nAHh48A8t+5xiSMNV\nAi2eawilta79pFInAZzWWp8IPJYD8Hq1HqNS6ow/x2iC31mt9dEwgSP/dFD/9aVX6vVvK6qUhmaU\ncZyTLydXgJfvey9cuxEy52fF5IsZfeoDuzo3/+WfXW74/aiJ2q5jPY5XSs3gSZlDmFWoS1HL5+kd\npylJUvzeDBc4+XJyBfj5NkpRJ7Ba6La67QVCA2NlAKw2HxBYbMnC686GHkMQBCHu2GzXsZoPMJOh\n02YIDXhL6R9EHVyj9lA+bpSSW9QKTnDy5eQK8PMV3IgMjKbH5/f65iqe+13F/UuuAoWOplzz3RaW\nDl+mVnCCky8nV4Cfr+AG+SWBnBLV+quMXODky8kV4OcruEEeGDklqvW3jXCBky8nV4Cfr+AGedox\nQRDiT1EnsJLfGyvONpD3GPPpVWoFa+4f+oRawQlOvpxcAX6+ghvkgTHBKB9j76PnqRWc4OTLyRXg\n5yu4QR4Yk6VOagVr/EvRuMDJl5MrwM9XcIM8MAr7CykfKnCAfPGlkHpMrWDNysiX1ApOxM3XD4rV\n6qDEzTUKbr6NUtIJrOf5THs1Sgx6jHyufNGqQK3gRJx8K3uKlffj5GoDN1/BDfLAmGI0tPLTW3Eh\nLr61hs9PlReNiast3Hw5Yi4xPqOUmjH/1sw4E9ZWKXVMKXXK3M6YnA+hkA+lhb1N1JyilBcVQmhW\nmVXnkivkPcZSYodawZrHfTepFZyIg29U0POfj4OrC9x8udGsMqsG55Ir5IGRU5XAtQyfgvBAfHxr\nBcfg47Vcc4Xe8i1OxOXcxpQRpdT1wO1U9Et20awyq8CTkiun/ATaUW9OHhjTeT7XnPop87kQJ9/K\n4Fh5P06uNnDzbZSiVniU77S6wVQJDNxmo45fhWaVWfWTZv8NvOHzaZtjyxyj0DaGUhsNzSnKfCR/\nTI8trNzJVa31HJpXZjVYcuW8ef+rEQ70gVGrErWCNcU0r1/KOPrWCmxxdA2Dm29ccOg9NqXMarWS\nK0qpoyGVTgHEYChdYJVE4mNqBSc4+dZyHUptPHWLC5zOLUdcyqpEtK2r5Ap5YEwxmmMcu/smtYIT\nnHw5uQL8fJnyrr83EcBJ7C6r8k5UWxMA54P7GGFRcoV8KK00eWy2xq9NzAVOvpxcAX6+jVLUCus7\n7U344lhWJaytc8kVPlFJEAShTZAHRk6Japee4VUAiZMvJ1eAn6/gBnlgTBb5pEvvX3mJWsEJTr6c\nXAF+voIb5IExUeqgVrCmZ/1ZagUnOPlycgX4+QpukAdGQRCEuEG+Kl2I0d60KHKjX1ArOMHJl5Mr\nwM+3UUpaYTOfptZoGzHoMSpqAWuUJv874gQnX06uAD9fwQ3ywJgq9FArWJN58Aq1ghOcfDm5Avx8\nBTfIA6MgCELcIA+MxcQ2tYI1G/0/UCs4wcmXkyvAz1dwgzwwlpJb1ArWbAx8T63gBCdfTq4AP1/B\nDfLAyClR7djd31ArOMHJl5MrwM+3UbRW2NxOW932ArK0JghCLDHpw07CSw5xDMCsqecS9pqrWusT\nVR4/GbwflViCPDBqVaRWsKbQwee6boCXLydXgJ8vU6yrBJp0Y9OoUjDLpBqb11pfMjkcrwGId2As\npB9RK1izPHmFWsGJRn2fqvvc4o34++3cCuFUq/xngl9VTDkEvzRq8DgZAL/XWg/5xwHwi6j3j5xj\ndCl6bdqfDN6ijp/acapqSMrBO7wKIDXT17VKn2tlv/18boWquFQJDOM4vES1JwMxbDrqRTY9Rpfu\nrHOXVTG68iVR7KJWcILKNxgQbQtYybmNN1or5PNJ2+YjJlb4zNZRKdClSmAY0/DmJ+dMr/M6gBto\npBiWS3e23i6rwAebIFetlyjV/fYdD7TWx6s90aIqgWHMw+usrQDlGDatlJrWWs/XelFUj7Fmd7ZK\nMZlyl9W85hiAS2FvDgD5jtBFplixeOQjagUnGvV1DWZhQ+eo4Ljfzu1+pUVVAqOOU0lk0IkKjC4R\n2rrLav5qnAKAkbEhTNz0ato8GvoG+Y4cDiz9GgCw1XMPudHPMXHLG7nrRAGLRz7EyMIM0tue2vLk\nFXRvTKFv9ecAgNXhr1FMbWD4/q8AAJu9d7B64DrGb/8WgLehfOnwZYzeewOpHW8P5f1Dn6D30fPo\nXXsOALAy8iW0KmBo+ZcAgMd9N7GW+QaTN/8tiql1FNMbuH/oY4zdfbNc+2PpmcvoX3mpnKcvN/oF\nlE6Vr6nd6P8BGwPfl/e/FTpWsTx5BQfvvFUeli0e+QiDD4+je+MwACA79hmShV4MZl8GAKwPfofN\n3lsYvfcGACDfmcWDiTmM334bquR9lAtTFzG0/Cq6Hk8iWejD/UMfI70zhIGcl1h1LfMttrsXMbLw\nOgBgu+s+suOfeudYK0BpLExdxPDia+jcGgMAPJi4hs7NcfSvvBj6Of3MFKEvqjx+OPQ/ceT+v0D3\nzhAA4ObYp+hfO1rzc9KqhKXD/6Mpn9NB36OFn9PBO/+6XK+o0c8JAB4e/EPLPieOmNKn5fvVqgQC\nyEZt39FazyulVpRSGROXMvB6kKEdNqW1rv2k1/s7HdwXpJTKAXi9SsnCGQAXtNZHA49pAEfDJI4+\nP6X/8vJPwxxjw8TNd7Dw7N9Sa1hD4Vur1xjV+5Rz2zree+HajVpDW1u6jh7Sz/zNv7dq++Nf/OeG\n3w8oT+X5daGf2seolLoIb9g9W9H2HIDzeDIk94PoaQA/wuuoXYgKjFE9Rtei15XwGScLTWEotbEr\nOMr8olAPdVYJPF/lOPMAzrq8d2hgdOnO1ttlLaTWXXxJyY59Rq3gBJVvMDjaBkU5t/FGlxQK2+Tb\nntuGzU/6rtmG43dnK4teXwXgT6j+DsDvlVJ+l7Xqtp4gSltvASAn6biXjxpKX9deopxbIU5EBkbH\n7qxzl5VTlcDB7Mt4PPCP1BrWcPLl5Arw8xXc2D99Y0FgSrUFLZm3bS3kaceKjPIxrg9+R63ghI2v\n66V+rWIvnluBL+Q9xlJih1rBms3eW9QKTkT5xiUoAnvv3O45NKDz5P2otkH+k6bzA9QK1vgbdrkQ\n5hunoAjsrXMr8Ie8xyi0n7gFRSEcmU9sP+SBkVOi2nxnNrpRjKjlW7kJOw4bsvfKuRX2BuRDaU6J\nah9MzEU3ihFhvrWCX71BsdFe6F46t0Gkd84T8sCYNkkGODB++21qBSeifCuDYKM9xUaCwF47t8CT\n8yHBkR/kQ2lO+JlRuGDj6w+rmzF8buQYe+3cVgZD9jkptYLaJu9HtY3985MKNWlGT5H1L32TqdVD\nlJ6jGy5lVZRSx5RSp0y7i7XKF1TWhKkF+Z9pTvniFqYuUis4wcmXkytQ2zcq+MkfESesyqqYgHk8\nkIJsBl4Oh6MV7Wbg5YE9HfXG5IExWeijVrBmaPlV5ByyqlCv9tr61vplvbk1gr9fnirf7+vYBgCs\n73Q+1e7EhP1VILXeq5prnINIrXNbLe1a5fNCNI5VAqfh5Wjwk9lcBzDtZ/oyx8uY41mlQiQPjIlS\nmlrBGj/bMhdsfIMLBMFf2ptbI1jJd2Mp14/CQy9zte4sAcBTc02Hnluu+nqX96rmGtY2DoSd21rB\nMY4/R4yxLqti0iOeCDx0HMBKRXbvGVOkz+rNyQOjQEe1BQIAWC10YyXfjfV8FwrbKXQtJpFeA0od\nT1LEJXaAjWdL2Mynn3q9/8vfyJByLyxcVAZHbv670EAib13Rs+1VAivyvp5GID2i6Wk67QcjD4yc\nEtU+PPgHagUnwnxrBa5gULy3MVC+Prb7oUbno1K53cpPqufRdFlgCAY83zVs4SJOwcXmu9DMFX9m\nkFUJNMf/QGt9ydy3qg1TCXlg5JSoNr0zhJ3uJev21L8QYb7VhnvBoPgo34nNfBpqO4GkN7WIrqUt\npBayWP8zbxhZSuvy6wZTm3U5+oEjvTOEpXT4H8k4BRnb70JcfONCK6sEmp7hvF/rxXAMwLBSyg/U\nGRM858KqC5Bv1+GUqNav4MaFKN+h1Eb5F3cotYHB1CYy6U30pbcwkN7GaO+TQJXa8oJgabgfyS2v\n5+gwtAp18F2jgkicggy37wI3qhTb21VWJbh9x1+sCRTAOmmOc0lrPevfzGOzjRbDEmKEax0VV/zj\nDqY2sZr2/mD9Q984AK9Xv3G4G4D3eMc6sLmtyu1tjgsgct5NFi6EAFZlVUzQvAEAgcWVeQCX/Dsm\niJ4y/z+DiJr35IGRU6Latcy3ZO9dmfTBJlC4+PrH/Emnt8qMTmCl2IPuvm3k+zqxdih6cFEtANpe\nkx105bBwQfld2C/YllUxAS50+GLmGM+jShXBapAHRp3IUytYs929SPK+1XpQNsHRxjcYwDLJxxhO\nruNA0htCPyz24chwDrdfBB7lupF+6H1dSp3am1/sK5SPU+26a5eeX6Vr3BcuqL4LVKgSkNhufOqE\nC+RzjKl8P7WCNSMLr7f9PcNWeaNWgKN8awXFTGIb06ktHE0/xH88/Hf4i59+hX/+8z/i4J8tAVOP\nkXkuh0PPLeOVo7fw68kfmpKpp5prXIMiQPNdENoHeY9RoCEsqA4nNAYTnQC2sZJcx1BqA8/05LyV\n6r40RnvXMdmzikx6E892PWiftCC0CfLAWFKF6EYxYbvrftvfM+wSs8ihdIhv8LjVjj+d2kK25A2d\nMsnHAICB9DbQi5YERYpz2wjcfAU3yANjMb1GrWBNdvxTq3bNXj2uFhxtjm3rC9TuQT4seteyD6Y2\ngR7gT4+HWtJTdHGNA9x8BTfI5xg5JaqduLUrsccuwkoGNEIwENoG3CjfyuPkCr1YKfYgW+zDQ3PL\nFvuwUuwp73N8pifXkuGzzbmNE9x8BTfIe4ys0OGrcvWuHtvifJwIX/+YYcF8pdjzVNuWLYhYuMYK\nbr6NovfXqrQERheUrvlU1OoxyQpriG+Qam7BgNgWLF1jAzdfwQnywCiJaltHo742UwH1BPyqx91n\n51aIN+RzjElG+xiHF1+r+VxYgKDajxfm2wya+XO12rXZcPMV3CAPjAlN3mm1pnNrLPT5Wtf+UhHl\nGyc4uQL8fAU3+EQlJgQXM+J85YYgOKEBRlfvNgx5YCww2sf4YOKaVbu4BERb33qpd1Gp2mta7dps\nuPkKbkQOpV1KGFa8zqpMoWJU86Vzc5xawYlGfdsZ4PfbuRWiqaN86oxS6qRS6kKwfKptadUgNnOM\nF7XW500CyFl4JQyjfiC/TGEkyWKXTbNY0L/yIrWCE83wbdei0n48t0IkLrHnGoDrpqTBDQAXgadL\nq2qtzwO4AC+PYyihgbFaCUMAtUoY+q9xKlMoxJ+4LSoJe586Ys9PAnVdsoHH/dKqPuXSqmHvH9Vj\nrFnCMOQ1M2F1GSopJuurFULBo6FvqBWcaKZvPZckurCfz+0eZEQpdT1wsxo9VuAUeyqKXZ2GCYYm\nFkWVVt1F1OKLUwnDesoUalV0aU4Kp83oQPN9W9lL3O/nNu4ojXJRNAtqVgl0wCn2AOW6MCfxpNIg\ngPDSqrWICozWJQxdyhSavyCnAGBk9AAmbr4DwPsrnO/I4cDSrwEAWz33kBv9vHzBvk4UsHjkQ4ws\nzCC97WksT15B98YU+lZ/DgBYHf4axdQGhu//CgCw2XsHqweuY/z2bwEApeQWlg5fxui9N5DaGQQA\n3D/0CXofPY/etecAACsjX0KrAoaWfwkAeNx3E2uZbzB+57fId2RRTG/g/qGPMXb3TSTz3tacpWcu\no3/lJfSsPwsAyI1+AaVTyDx4BQCw0f8DNga+x9jd3wAACh2rWJ68goN33kLCzLMuHvkIgw+Po3vj\nsHfyxz5DstCLwezLAID1we+w2XsLo/feAADkO7N4MDGH8dtvQ5W8j3Jh6iKGll9F1+NJpHeGsXj4\nI6R3hsrFm9Yy32K7e7GcaHW76z6y459651grQGksTF3E8OJr5b16DyauoXNzvDyv1orPKVnow93p\n/9qUz+ngn94CgJZ+TmN330TRlP5t9HMCvHKsrfqc4kQry6eaAHjeLLRc1VoHe4q7SquGempd+5pP\n0219X2v9i8BjOa31rpQ4pipXUPwCvOgcWqbw6PNT+i8v/zTKMxZM3HwHC8/+LbWGNWM3/w3+4ZnL\nVm2p5wy5nVtOvu+9cO1Goz247onDevrf/ZVV2/977q8afj/H2DMN4KRZXPHXOXIAjvqxx4xmUVFa\ntSahPUat9VeBqltVSxjC9BIro7BS6oJNDdkSo12jWz33qBWcWO+yr0tCXV+F27nl5ssNl9gDbz7y\nQODl0/DmEf2g6JdW/crcPxnVa7TZ4G1VwjAg7FSm0B+OcCA3+jm1ghN3D3zp1J4yOHI7t9x8mWIV\ne7TWc0qpTGCR5wSA14FyAA0trVqNyMBoW8Iw8JhTmUJuiWq5DJ8A4Gd3/5X1UNqHKjhyO7fcfBum\nBCR22vuWLrGnogc4G3g8srRqNcgvCRRaRxLFqkEuKp0Y9bBaEKiRwOiATvAp3AXU9g0rsOU/D9iX\nZmhGEG3VuXWpbe0Ct++C4AZ52rG4bScIY/HIh9QKToT5NqMWdDPZS+dW4A95YEzlB6gVrBlZCL0a\nMnZE+VYGQcrh8147twJvyAOj0klqBWv8zcpcsPH1gyH1nOJePLcCX2SOUSAPikL8URpItnlVmhLy\nwJhPP6JWsGZ58gq1ghON+rYzYLbq3LbqZ+D2XRDcIB9KJ0od1ArWdG9MUSs4wcmXkyvAz1dwgzww\nckpU6ydA4AInX06uAD9fwQ3ywCgIghA3yOcYOSWqXR3+mlrBCU6+nFwBfr6Nst8WX8h7jJwS1RaZ\nrd5y8uXkCvDzFdwgD4ypQh+1gjV+UlUucPLl5Arw8+VIsyqU1nMc8qG0IAhCDS76iWqVUtfhVQnc\nldErSKBC6elGjkPeY+SUqHaz9w61ghOcfDm5Avx8udGsCqX1HAeIQWDkNFezeuA6tYITnHw5uQL8\nfNtM26sEGqpVKK3nOPSBMb1jNW0QC/xCTVzg5MvJFeDn2zAaSG5rqxtMlcDALbLESRWaVaG0rova\nZY5REIS20YoqgREVSp2rDQIxCIwatasUxo1ScotawQlOvpxcAX6+ccGh9ziPKgGsylAZ8OrBDCul\n/MqEfv2XOcfjlCEPjIWOyDLUsWHpsFv9FGo4+XJyBfj5cqOZFUrDjlML8jlGTolq/QLqXODky8kV\n4OfLlHf9/YcATmJ3lcB3go1NpcAz5v9nTBCMOk5VyHuMnBLVpnYGqRWc4OTLyRXg58uRZlUoDTtO\nLcgDoyAI8UeVgNQWn/WARiEfSufTq9QK1tw/9Am1ghOcfDm5Avx8BTfIA2OCUT7G3kfPUys4EUff\nXKG3aknTOLqGwc1XcIM8MCZLndQK1vSuPUet4ETcfIMBsTI4xs01Cm6+ghvkgVHYH1TrJVZ7TBDi\nAPniSyH1mFrBmpWRL6kVnIiLb1gAzBV6MZTaiI2rLdx8G8VLVCuLL22Ez8nWqkCt4AQnX06uAD9f\nwQ3ywJhiNJwaWv4ltYITcfENK2HqPxcXV1u4+QpukAdGYX9QLTi2s261ILhAHhhLCT4Vdh733aRW\ncCJuvsFAWBkU4+YaBTdfwQ3yxRdOVQLXMt9QKzgRR99avcQ4uobBzVdwg7zHmM7zueb04J/eolZw\ngpMvJ1eAn2/DlDSSWyWr216AvMcoCIJQDZMd5yS8BBDHAMzWSEZb+boLWuvTgfvHAJRzNQK4pLWe\nDztGZGB0kQsIZAC8AuBslIBWfP7CFNO8Fgs4+XJyBfj5MqVZVQJntNbnA20uVDy/C5seo5WcqdBV\nru9gBK8iPI05CqySSHxMreAEJ19OrgA/X25Uq+5nYkrYa3ZVCTScVkpZ9TZ9QucYHUsPTgM4G7h/\nHcB0VHHs2oZrAAALQUlEQVTrFKM5xrG7b1IrOMHJl5MrwM+XIc2qEgh4SW3/qJQ6ZUoenK3S5imi\neow15SoFTCryE4GHjgNYiYrSSpOv/1iTzPPZjA7w8uXkCvDzbRSlgeR20bb5iBld+szWUSmwWVUC\nobWeNR00f/g8h91x7SmiAqOTXMV84mnUSCFuovYpABgZPYCJm16G8kdD3yDfkcOBpV8DALZ67iE3\n+jkmbnkjd50oYPHIhxhZmEF621NbnryC7o0p9K3+HACwOvw1iqkNDN//FQCvMPrqgevlcpel5BaW\nDl/G6L03ylmY7x/6BL2Pni9nTFkZ+RJaFcpXNzzuu4m1zDdI7wxj4uY7KKY3cP/Qxxi7+2b5F2Tp\nmcvoX3kJPevPAgByo19A6RQyD14BAGz0/4CNge8xdvc3AIBCxyqWJ6/g4J23yqnXFo98hMGHx9G9\ncRgAkB37DMlCLwazLwMA1ge/w2bvrXJa/XxnFg8m5jB++22okvdRLkxdxNDyq+h6PIn0zjA6Ng8i\nvTOEgdxLAIC1zLfY7l7EyMLrAIDtrvvIjn/qnWOtAKWxMHURw4uvoXNrDADwYOIaOjfH0b/yYss+\np2ShD6qUbsrn5K8Yt/JzShb6yt/bRj8nAHh48A8t+5wIeKC1Pl7tCYIqgVBKnTFzjOfN+0dO8Smt\na1+rrJQ6CeC01vpE4LEcgNfDqmyZN89WFqmpxpEXM/qv/3vVcxg7EoUulFJ8qsNx8uXkCvDyfe+F\nazdqBSpbBgae0a8c/w9Wbf/uf/+nht/PDJnf99c3zGM5rfVQlbYn8XTQ9BdX5uCNemGCrd/+HIAP\nwmJY1DjWufSg6dLO2wRFAEgWu22axYL+lZeoFZzg5MvJFeDny43KGFOtSqC/fqG1vqS1nvVv5rFZ\nM4LNwttNE3r8SkIDo4ucuX8MXk9xztw/GXZ8AEiUOqKaxAZ/+MUFTr6cXAF+vkxpuEqgiWHz/sKL\nef6DqDe22a7zrjmYv4+xUu4qgFkTNG8YKf/5eQBWPUdBEIQgTawS6ByDIgOjrZzptio4UmCUYSU3\n+gW1ghOcfDm5Avx8G6akkdiyXpVmTwz2yjjHUjKU5nUFJSdfTq4AP1/BDfLAmCr0UCtY42/p4AIn\nX06uAD9fwQ3ywCgIghA3yANjMbFNrWDNRv8P1ApOcPLl5Arw8xXcIJ8oKSV5bJIFgI2B76kVnODk\ny8kV4OfbKEprJLbz1Bptg7zHyClRrX+ZGBdcfHOF3qdurlS+3vUYcm6FOEEeGAVBEOIGeWDUis/e\nqEIHn9yRAC9fTq4AP1/BDfLAWEg/olawZnnyCrWCE5x8ObkC/HwFN8gDY2onNI9trDh4h1cBJE6+\nnFwBfr6CG+Sr0orRlS9+Pj4uuPjWKmvartfLuY05JY3EpqxKC4Ig7FvIe4z5Duv6NOQsHvmIWsEJ\nTr6cXAF+vhxxrFB6DsAZeCULrsNLsD0feP6pFIhRGXfIe4xJRvu6Bh/yyDTuw8mXkyvAz5cpF7XW\n501+11l4FUpr8aPWWmmth7TWJyqC4hmgHAznAPw+6o3JA2OilKZWsMavxcIFTr6cXAF+vtxwrFAa\ndpwMgN/7PUSt9UqwXEItyAOjIAh7jhGl1PXA7VQdx3Atn5pRSp1USs0opc4FKgsch5fB23/ujBmi\nh0I+x1hIrVMrWJMd+4xawQlOvpxcAX6+DVMqAZvWeQ1qVgl0wKlCKQLzj0qpLIBrAH4BL8AeAzCn\ntV4xZV1vIKJKIHlgVDpJrWANp/lQgJcvJ1eAn29caEX5VKA81Pb//5VS6pjpNc7DK8634rcztaqm\nK8o9PwV5YORUJXAw+zIeD/wjtYY1nHw5uQL8fOOCX8XPAusKpdVKrZq2K0qpasEvciuMzDEKghA7\nHCuUzsOrJe0/NwNThM/0Clf8tn4vMqy3CMSgx1hklI9xffA7agUnOPlycgX4+TLFqkKp3zMMLPIc\nrWj7OwC/V0r9aJ7bVWGwEvLAWErsUCtYs9l7i1rBCU6+nFwBfr4NU9LAVnuz7buUT/Vr2dc4zjyA\nsy7vTT6UTucHqBWsGb33BrWCE5x8ObkC/HwFN8gDoyAIQtwgD4ycEtXmO7PRjWIEJ19OrgA/X8EN\n8sDIKVHtg4ma0xixhJMvJ1eAn6/gBnlgTO8MUStYM377bWoFJzj5cnIF+PkKbpCvSnNClXidLk6+\nnFwBfr4NUypBP96ktmgb5D1GQRCEuEEeGPMdOWoFaxamLlIrOMHJl5MrwM9XcIM8MCYLfdQK1gwt\nv0qt4AQnX06uAD9fwQ3ywMgpUW3X40lqBSc4+XJyBfj5Cm7ssxlkQRDqQWuNkn0+RvZEBkbHgjTW\nbX04Jap9ePAP1ApOcPLl5Arw8+WIazyxKXillLqgtT4d9d42PcaLfp4zk/32fdTOTuHS1hNllKg2\nvTOEne4lag1rOPlycgX4+TLFOp6YLDzzWutLJrXYNZjUY4E2MwBOAYgMjKFzjC4FaeotXsMpUe1A\n7iVqBSc4+XJyBfj5csMx9kQWvDJtsrBIUgtEL764FKRxLV4jCIJQC5d4YlPwaqZa9u9aRA2lXQrS\nWLc1CSX9pJLb771w+/84vA8h10YAPKC2sIeTLydXgJnvzxo9wJrOXrma/28jls27zNDXZ9ahpIGP\nS+wJLXhlhtBOF7dHBUaXgjTWbc1JmgW8uYMmVBRrC5xcAV6+nFwBXr4VQaoutNZ/3iSXVhTDqlnw\nyjyfjVoEriQqMFoXpHFsKwjCPqQVxbBM20r8QHgMwLBSyv8jljHBea7uKoGmDGH5frWCNDDROKqt\nIAiCLY6xZ14ptaKUypjeYrDg1VPBz2zXiQzONtt1rArSWLSthevcAyWcXAFevpxcAV6+nFyDuMSe\n0IJXJlieMv8/A+BSWI9Raa2b9UMIgiDsCcivlRYEQYgbEhgZY1bezgT2blWu4tV63YXoVgInlFJX\nLdrU9X3Zj7R8KN3qa60JXY/B21iaAfAKgLNhcxatQCl1I3DJVAbA+5X1dqu8ZgbelggV1q7ZtOK6\n11ZS53cB8L4PofNXTfacgbeP70LUZ1rP92XforVu6Q3AjcD/M/Cuf2y4LaWree5U4P4MgB/b7HoM\nXoALPpaLeE3GvC60XQy+B2cAnAy0vdFKt2b4Vty/QOCrm/192c+3lg6l23GtdbNwfP9pAGcD968D\nmG7z0KSeSzCdLotqFs2+7rXV1PFdPM1gWCqX7DrQ6jlGTtdaW7+/CS4nAg8dB7Ci2zjsh9slU3Vd\nFtVEmn3da6tx/S6eA/BHpdQps3n4bI12lDh9X/Y7rU5U25JrrVuE0/vrp+eQTsNuz2Yzsb5kKrgZ\ntuVW1Wnada9twvW7MGt6jH46qzlYZnFpIy6X2O17Wt1jbMm11i2irvc3PYQPdJsXB+B2ydQxAMcD\nPZqM+X+7emINXfcKb5qinb1Gp++CUuqM1vq8GfJfgLfxOG7IJbsOtLrHyOlaa+f3N8PTee1d9N5W\ntNslU5UJO60ui2oizbrutV1Y+5rvwFeBNrNKqaNKqWPUQUcu2a2flvYYK78Y1X55/UnrqLatxsXV\n3D8G70s3Z+4/tb2kTbzr70uDt7Wk8pKpd4KNlVIZczkU2jl35/g9mAfgX+/qL8bM6zZuhXL8LmTh\n9chDj9EqlFLHAp/pOfNd8Kn8DoR9X4QA7djHeAzeit6u/WBKqYvwthDMRrVtB7au5hflx4qXz2ut\n2zkPxgrH78E0vPk6/7rXC+0MjHX4nsSTHmYG3vyoDFEZI9dKC4IgVCCXBAqCIFQggVEQBKECCYyC\nIAgVSGAUBEGoQAKjIAhCBRIYBUEQKpDAKAiCUIEERkEQhAokMAqCIFTw/wFuDu+E63QqZAAAAABJ\nRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f401cad1110>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"pl = plot(u)\n",
"plt.colorbar(pl)\n",
"plt.show()\n",
"\n",
"pl = plot(local_project(p, V))\n",
"plt.colorbar(pl)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"ename": "RuntimeError",
"evalue": "\n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n*** fenics-support@googlegroups.com\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a *minimal* running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error: Unable to solve nonlinear system with NewtonSolver.\n*** Reason: Newton solver did not converge because maximum number of iterations reached.\n*** Where: This error was encountered inside NewtonSolver.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2017.2.0\n*** Git changeset: unknown\n*** -------------------------------------------------------------------------\n",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mRuntimeError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-14-c6af918a3cf7>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mobstacle\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0md\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mh\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minterpolate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobstacle\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0msolve_Uzawa\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 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0mpl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m<ipython-input-12-d2d4361d090f>\u001b[0m in \u001b[0;36msolve_Uzawa\u001b[0;34m(tol)\u001b[0m\n\u001b[1;32m 43\u001b[0m \u001b[0mtold\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 44\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 45\u001b[0;31m \u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mform\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 46\u001b[0m \u001b[0;31m#gap.assign(project(h-u, P))\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 47\u001b[0m \u001b[0;31m#p.vector()[:] = np.minimum(0*p_old.vector().get_local(), rho.vector().get_local() +\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/usr/lib/python2.7/dist-packages/dolfin/fem/solving.pyc\u001b[0m in \u001b[0;36msolve\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 298\u001b[0m \u001b[0;31m# tolerance)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 299\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0misinstance\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[0mufl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclasses\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mEquation\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 300\u001b[0;31m \u001b[0m_solve_varproblem\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[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 301\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 302\u001b[0m \u001b[0;31m# Default case, just call the wrapped C++ solve function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/usr/lib/python2.7/dist-packages/dolfin/fem/solving.pyc\u001b[0m in \u001b[0;36m_solve_varproblem\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 347\u001b[0m \u001b[0msolver\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mNonlinearVariationalSolver\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 348\u001b[0m \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparameters\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msolver_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 349\u001b[0;31m \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\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 350\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 351\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mRuntimeError\u001b[0m: \n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n*** fenics-support@googlegroups.com\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a *minimal* running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error: Unable to solve nonlinear system with NewtonSolver.\n*** Reason: Newton solver did not converge because maximum number of iterations reached.\n*** Where: This error was encountered inside NewtonSolver.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2017.2.0\n*** Git changeset: unknown\n*** -------------------------------------------------------------------------\n"
]
}
],
"source": [
"for d in np.linspace(1e-4, 1e-3, 11):\n",
" obstacle.y = d\n",
" h.interpolate(obstacle)\n",
" solve_Uzawa()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Calling FFC just-in-time (JIT) compiler, this may take some time.\n"
]
},
{
"data": {
"text/plain": [
"<dolfin.cpp.la.Vector; proxy of <Swig Object of type 'std::shared_ptr< dolfin::Vector > *' at 0x7f401e7570c0> >"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"assemble(form)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.12"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Hertzian contact with a rigid indenter using a penalty approach\n",
"\n",
"In this numerical tour, we explore the formulation of frictionless contact between a rigid surface (the indenter) and an elastic domain, representing an infinite half-space for the present case. Contact will be solved using a penalty formulation, that is small interpenetration between the solid and the identer will be authorized. Prerequisites for this tour are the formulation of an linear elasticity problem and the resolution of a general nonlinear variational problem."
]
},
{
"cell_type": "raw",
"metadata": {
"raw_mimetype": "text/restructuredtext"
},
"source": [
".. figure:: contact_problem.png\n",
" :scale: 20%\n",
" :align: center\n",
" \n",
" Hertzian contact problem of a rigid spherical indenter on a semi-infinite domain"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The elastic (isotropic $E,\\nu$) solid will be represented by a 3D cubic domain of unit dimension and contact will take place on its top surface $z=0$, centered at $x=y=0$. Symmetry conditions will be applied on the $x=0$ and $y=0$ surfaces whereas the bottom surface $z=-1$ will be fully fixed. If contact appears on a small region of extent $a\\ll 1$, the problem can then be considered to be a good approximation of contact on a semi-infinite domain.\n",
"\n",
"The rigid indenter will not be explictly modeled (in terms of mesh in particular) but instead its distance with respect to the solid top surface will be given as an ``Expression``. We will also consider that the indenter radius $R$ is sufficiently large with respect to the contact region characteristic size $a$ so that the spherical surface can be approximated by a parabola. In this case, the distance between such an indenter and the top surface can be written as:\n",
"\n",
"\\begin{equation}\n",
"h(x,y) = h_0 + \\dfrac{1}{2R}(x^2+y^2)\n",
"\\end{equation}\n",
"\n",
"where $h_0$ is the initial gap between both surfaces at $x=y=0$. Obviously, if $h_0>0$ there is no contact between both surfaces, contact appears only if $h_0=-d<0$ where $d$ will be the indenter depth inside the surface. This classical problem admits the following known analytical solution **[REF?]**:\n",
"\n",
"* the contact area is of circular shape and radius $a=\\sqrt{Rd}$\n",
"* the force exerted by the indenter onto the surface is $F=\\dfrac{4}{3}\\dfrac{E}{1-\\nu^2}ad$\n",
"* the pressure distribution on the contact region is given by $p(r) = p_0\\sqrt{1-(r/a)^2}$ where $p_0=3F/(2\\pi a^2)$ is the maximal pressure\n"
]
},
{
"cell_type": "raw",
"metadata": {
"raw_mimetype": "text/restructuredtext"
},
"source": [
".. figure:: contact_solution.png\n",
" :scale: 20%\n",
" :align: center\n",
" \n",
" Indent at depth :math:`d` (left) and contact solution (right)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Contact problem formulation and penalty approach\n",
"\n",
"The unilateral contact (Signorini) condition on the top surface $\\Gamma$ writes as:\n",
"\n",
"\\begin{equation}\n",
"g \\geq 0, p\\geq 0, g\\cdot p =0 \\text{ on }\\Gamma\n",
"\\end{equation}\n",
"\n",
"where $g=h-u$ is the gap between the obstacle surface and the solid surface and $p$ is the pressure. One of the most simple way to solve approximately this contact condition consists in replacing the previous complementary conditions by the following penalized condition:\n",
"\n",
"\\begin{equation}\n",
"p = k\\ppos(-g)\n",
"\\end{equation}\n",
"\n",
"where $k$ is a large penalizing stiffness coefficient. With the previous relation, the pressure will be positive but a small negative gap will be authorized. "
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
"Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
"Maximum pressure FE: 1.42763078681 Hertz: 1.39916433487:\n",
"Applied force FE: 0.0320642723007 Hertz: 0.029304029304:\n"
]
}
],
"source": [
"from dolfin import *\n",
"import numpy as np\n",
"\n",
"N = 30\n",
"mesh = UnitCubeMesh.create(N, N, N/2, CellType.Type_hexahedron)\n",
"mesh.coordinates()[:, :2] = mesh.coordinates()[:, :2]**2\n",
"mesh.coordinates()[:, 2] = -mesh.coordinates()[:, 2]**2\n",
"\n",
"class Top(SubDomain):\n",
" def inside(self, x, on_boundary):\n",
" return near(x[2], 0.) and on_boundary\n",