NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1.6 Error estimation & adaptive refinement\n",
    "\n",
    "\n",
    "In this tutorial, we apply a Zienkiewicz-Zhu type error estimator and run an adaptive loop with these steps:\n",
    "$$\n",
    "\\text{SOLVE}\\rightarrow\n",
    "\\text{ESIMATE}\\rightarrow\n",
    "\\text{MARK}\\rightarrow\n",
    "\\text{REFINE}\\rightarrow\n",
    "\\text{SOLVE} \\rightarrow \\ldots\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometry\n",
    "\n",
    "The following geometry represents a heated chip embedded in another material that conducts away the heat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7fef7f27c7ec45778f0709fc66984142",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 172,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def MakeGeometryOCC():\n",
    "    square = Circle( (0,0), 2).Face() # create a rectangle with boundary conditions\n",
    "    square.edges.name = \"outer\"\n",
    "    square.faces.name = \"air\"\n",
    "\n",
    "    el1 = MoveTo(-0.4, 0.2).Rectangle( 0.8, 0.1).Face() # create a rectangle with boundary conditions\n",
    "    el1.edges.name = \"el1\"\n",
    "    el1.vertices.name = \"el1\"\n",
    "\n",
    "    el2 = MoveTo(-0.4, -0.2).Rectangle( 0.8, 0.1).Face() # create a rectangle with boundary conditions\n",
    "    el2.edges.name = \"el2\"\n",
    "    el2.vertices.name = \"el2\"\n",
    "\n",
    "    dielec = MoveTo(-0.9, -0.1).Rectangle( 1.8, 0.3).Face()\n",
    "    dielec.faces.name = \"dielec\"\n",
    "\n",
    "    air = square  # subtract the rectangles from the air rectangle\n",
    "    shape = Glue([air - dielec, dielec])\n",
    "    shape =shape - el1 - el2\n",
    "\n",
    "    ### adding extra specifications on the shape\n",
    "    #predefined mesh size for the shape\n",
    "\n",
    "\n",
    "    return OCCGeometry(shape, dim=2)\n",
    "\n",
    "mesh = Mesh(MakeGeometryOCC().GenerateMesh(maxh=0.5))\n",
    "mesh.Curve(3)\n",
    "Draw(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spaces & forms\n",
    "\n",
    "The problem is to find $u$ in $H_{0,D}^1$ satisfying \n",
    "\n",
    "$$\n",
    "\\int_\\Omega \\varepsilon \\nabla u \\cdot \\nabla v = \\int_\\Omega f v \n",
    "$$\n",
    "\n",
    "for all $v$ in $H_{0,D}^1$. We expect the solution to have singularities due to the nonconvex re-enrant angles and discontinuities in $\\lambda$.\n",
    "\n",
    "**On the notes adaptivity is explained using $\\lambda$ instead of $\\varepsilon$.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=3, dirichlet=\"el.*\")\n",
    "u, v = fes.TnT()\n",
    "\n",
    "# one heat conductivity coefficient per sub-domain\n",
    "eps = CoefficientFunction([1, 2])\n",
    "a = BilinearForm(eps*grad(u)*grad(v)*dx)\n",
    "\n",
    "# heat-source in inner subdomain\n",
    "f = LinearForm(fes)\n",
    "\n",
    "c = Preconditioner(a, type=\"multigrid\", inverse=\"sparsecholesky\")\n",
    "\n",
    "electrode = mesh.BoundaryCF({\"el1\":1, \"el2\":-1}, default = 0 ) # define the boundary conditions\n",
    "gfu = GridFunction(fes)\n",
    "gfu.Set(electrode, definedon=mesh.Boundaries(\"el.*\")) # set the boundary conditions\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the linear system is not yet assembled above.\n",
    "\n",
    "### Solve \n",
    "\n",
    "Since we must solve multiple times, we define a function to solve the boundary value problem, where assembly, update, and solve occurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveBVP(gfu):\n",
    "    fes.Update()\n",
    "    gfu.Update()\n",
    "    a.Assemble()\n",
    "    f.Assemble()\n",
    "    inv = CGSolver(a.mat, c.mat)\n",
    "    gfu.vec.data +=inv* (f.vec - a.mat * gfu.vec) # solve the system\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "70b23b7eb62f46b8945c7ed40cc99689",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "SolveBVP(gfu)\n",
    "Draw(gfu);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate\n",
    "\n",
    "We implement a gradient-recovery-type error estimator. For this, we need an H(div) space for flux recovery. We must compute the flux  of the computed solution and interpolate it into this H(div) space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [],
   "source": [
    "space_flux = HDiv(mesh, order=2)\n",
    "#space_flux = VectorH1(mesh, order=2)\n",
    "\n",
    "gf_flux = GridFunction(space_flux, \"flux\")\n",
    "\n",
    "flux = eps * grad(gfu)\n",
    "gf_flux.Set(flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Element-wise error estimator:** On each element $T$, set \n",
    "\n",
    "$$\n",
    "\\eta_T^2 = \\int_T \\frac{1}{\\varepsilon} \n",
    "|\\varepsilon \\nabla u_h - I_h(\\varepsilon \\nabla u_h) |^2\n",
    "$$\n",
    "\n",
    "where $u_h$ is the computed solution `gfu` and $I_h$ is the interpolation performed by `Set` in NGSolve.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f92a80dd8eb747e7b47313bec5dda428",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 177,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "err = 1/eps*(flux-gf_flux)*(flux-gf_flux)\n",
    "Draw(err, mesh, 'error_representation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[4.21842882e-08 3.51994462e-08 1.15858338e-07 9.16675035e-07\n",
      " 8.56871323e-08 4.69599389e-07 1.35455371e-07 1.08785321e-06\n",
      " 3.81165502e-08 2.91867924e-07] ...\n"
     ]
    }
   ],
   "source": [
    "eta2 = Integrate(err, mesh, VOL, element_wise=True)\n",
    "\n",
    "# print only the first 10 values\n",
    "print(np.array(eta2)[:10], \"...\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above values, one per element, lead us to identify elements which might have large error.\n",
    "\n",
    "\n",
    "### Mark \n",
    "\n",
    "We mark elements with large error estimator for refinement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "maxerr =  0.17618067298119255\n"
     ]
    }
   ],
   "source": [
    "maxerr = max(eta2)\n",
    "print (\"maxerr = \", maxerr)\n",
    "\n",
    "for el in mesh.Elements():\n",
    "    mesh.SetRefinementFlag(el, eta2[el.nr] > 0.25*maxerr)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Refine & solve again \n",
    "\n",
    "Refine marked elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7d3615ee88764461a39972a00d6a1e3e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 180,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mesh.Refine()\n",
    "SolveBVP(gfu)\n",
    "Draw(gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automate the above steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = []    # l = list of estimated total error\n",
    "\n",
    "def CalcError():\n",
    "\n",
    "    # compute the flux:\n",
    "    space_flux.Update()      \n",
    "    gf_flux.Update()\n",
    "    flux = eps * grad(gfu)        \n",
    "    gf_flux.Set(flux) \n",
    "    \n",
    "    # compute estimator:\n",
    "    err = 1/eps*(flux-gf_flux)*(flux-gf_flux)\n",
    "    eta2 = Integrate(err, mesh, VOL, element_wise=True)\n",
    "    maxerr = max(eta2)\n",
    "    l.append ((fes.ndof, sqrt(sum(eta2))))\n",
    "    print(\"ndof =\", fes.ndof, \" maxerr =\", maxerr)\n",
    "    \n",
    "    # mark for refinement (vectorized alternative)\n",
    "    mesh.ngmesh.Elements2D().NumPy()[\"refine\"] = eta2.NumPy() > 0.25*maxerr\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ndof = 1652  maxerr = 0.09112441447041292\n"
     ]
    }
   ],
   "source": [
    "CalcError()\n",
    "mesh.Refine()\n",
    "mesh.Curve(4);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the adaptive loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ndof = 2261  maxerr = 0.040022812393294735\n",
      "ndof = 3026  maxerr = 0.01686238424438765\n",
      "ndof = 3884  maxerr = 0.00716470226976139\n",
      "ndof = 4763  maxerr = 0.0030597661863823703\n",
      "adaptive step # 5\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1b5d44ee657c49d2b82dff207e0f47fc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ndof = 5519  maxerr = 0.0013107177162579886\n",
      "ndof = 6428  maxerr = 0.0005624494603309932\n",
      "ndof = 7277  maxerr = 0.00024166036673661242\n",
      "ndof = 8507  maxerr = 0.00010387514474365285\n",
      "ndof = 9851  maxerr = 4.4691827734222984e-05\n",
      "adaptive step # 10\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "57481c0c1d004b478f7da7d74de55832",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ndof = 11477  maxerr = 1.9224574256791775e-05\n",
      "ndof = 12926  maxerr = 8.273945685310195e-06\n",
      "ndof = 14885  maxerr = 3.55981350896467e-06\n",
      "ndof = 17162  maxerr = 1.5322545647250393e-06\n",
      "ndof = 20009  maxerr = 9.539191899191035e-07\n",
      "adaptive step # 15\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "57547913d49a4142b9a6f38afaaad8fe",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ndof = 21641  maxerr = 3.8866562381930035e-07\n",
      "ndof = 26072  maxerr = 1.6733131295093797e-07\n",
      "ndof = 31460  maxerr = 6.298646211685861e-08\n",
      "ndof = 38234  maxerr = 2.516914093622385e-08\n",
      "ndof = 45122  maxerr = 1.049478136862962e-08\n",
      "adaptive step # 20\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "18bd38e0b75941cb89929066da90ef4e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ndof = 53699  maxerr = 4.390017091991425e-09\n",
      "ndof = 64940  maxerr = 1.8902728268024055e-09\n",
      "ndof = 78116  maxerr = 8.385455776480657e-10\n",
      "ndof = 93812  maxerr = 4.85783536567953e-10\n",
      "ndof = 105272  maxerr = 1.7098027943176962e-10\n"
     ]
    }
   ],
   "source": [
    "level = 0 \n",
    "while fes.ndof < 100000:  \n",
    "    SolveBVP(gfu)\n",
    "    level = level + 1\n",
    "    if level%5 == 0:\n",
    "        print('adaptive step #', level)\n",
    "        Draw(gfu)\n",
    "    CalcError()\n",
    "    mesh.Refine()\n",
    "    mesh.Curve(4)\n",
    "\n",
    "SolveBVP(gfu)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot history of adaptive convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAG1CAYAAAAV2Js8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABjdUlEQVR4nO3deVxU9f7H8deZYRNQcEVB3JdEEUzRTDMty9RMbbOyMivb1Bazrv663TZNW6+p3BZb1Faz1DbTzH3fENTcd3DfANlhZn5/TGIoCijDGeD9fDx4XOfMmXPe5kU+fs/3+/0YDofDgYiIiEg5ZDE7gIiIiIhZVAiJiIhIuaVCSERERMotFUIiIiJSbqkQEhERkXJLhZCIiIiUWyqEREREpNxSISQiIiLllofZAdyd3W7n0KFDVKxYEcMwzI4jIiIiheBwODhz5gzBwcFYLBcf91EhVIBDhw4RGhpqdgwRERG5DPHx8dSuXfui76sQKkDFihUB53/ISpUqmZxGRERECiM5OZnQ0NDcn+MXo0KoAGcfh1WqVEmFkIiISClT0LQWTZYWERGRckuFkIiIiJRbejQmIiLlhs1mIzs72+wYUgw8PT2xWq1XfB0VQiIiUuY5HA6OHDlCYmKi2VGkGAUGBlKzZs0r2t5GhZCIiJR5Z4ugGjVq4Ovrq33hSjmHw0FaWhrHjh0DoFatWpd9LRVCIiJSptlsttwiqGrVqmbHkWJSoUIFAI4dO0aNGjUu+zGZJkuLiEiZdnZOkK+vr8lJpLid/TO9knlfKoRERKRc0OOwsqc4/kzLRSH066+/0rRpUxo3bsynn35qdhwRERFxE2V+jlBOTg7Dhg1j4cKFBAQE0Lp1a/r27avnxCIiIlL2R4TWrFlD8+bNCQkJwd/fn+7du/PHH3+YHUtERKRE7Nu3D8MwiI2NNS2DYRjMmjXLtPtfitsXQkuWLKFXr14EBwdf9D9kdHQ09erVw8fHh3bt2rFmzZrc9w4dOkRISEju65CQEA4ePFgS0UVERC7bmDFjiIqKomLFitSoUYM+ffqwffv2S37moYceok+fPiUTsIxw+0IoNTWViIgIoqOj831/2rRpDBs2jFdeeYWYmBgiIiLo1q1b7t4CRZWZmUlycnKeL5fZuwSSElx3fRERKbUWL17M4MGDWbVqFfPmzSM7O5ubb76Z1NRUs6PlKysrq1Re2+0Loe7duzNq1Cj69u2b7/vvv/8+gwYNYuDAgYSFhfHRRx/h6+vL559/DkBwcHCeEaCDBw8SHBx80fuNGTOGgICA3K/Q0NDi/Q2dlX4apg+EiVGw7L+Q47o/ZBERyUdW6sW/sjOKcG564c4tojlz5vDQQw/RvHlzIiIimDx5MgcOHGD9+vX5nv/qq68yZcoUfvrpJwzDwDAMFi1alPv+nj176NKlC76+vkRERLBy5cpL3v/AgQP07t0bf39/KlWqxN13383Ro0fz3C8yMpJPP/2U+vXr4+PjA8DOnTvp1KkTPj4+hIWFMW/evAuuHR8fz913301gYCBVqlShd+/e7Nu3L/f9syNbo0ePJjg4mKZNmxbhv1zRlOrJ0llZWaxfv56RI0fmHrNYLHTt2jX3D7ht27Zs3ryZgwcPEhAQwO+//87LL7980WuOHDmSYcOG5b5OTk52TTGUkQxVG0H8KvjzVYj9Fnq+C/U7Ff+9RETkQm9e/B/FNL4Z+k8/9/qdRpCdlv+5dTvCwN/OvR4XDmknLzzv1aTLy/m3pCTn56tUqZLv+8OHD2fr1q0kJyfzxRdf5J576NAhAF566SXeffddGjduzEsvvcS9997Lrl278PC4sBSw2+25RdDixYvJyclh8ODB9OvXL09xtWvXLn788UdmzJiB1WrFbrdz++23ExQUxOrVq0lKSuLZZ5/Nc+3s7Gy6detG+/btWbp0KR4eHowaNYpbbrmFjRs34uXlBcD8+fOpVKlSvoVUcSrVhdCJEyew2WwEBQXlOR4UFMS2bdsA8PDw4L333qNLly7Y7XZefPHFS64Y8/b2xtvb26W5AahcFwb+DnHfwrz/wIntMKUXhN8FN4+CijVdn0FEREoFu93Os88+S4cOHWjRokW+5/j7+1OhQgUyMzOpWfPCnyHDhw+nZ8+eALz22ms0b96cXbt2cdVVV11w7vz589m0aRN79+7NHQyYOnUqzZs3Z+3atURFRQHOAYmpU6dSvXp1AP744w+2bdvG3Llzc5++vPnmm3Tv3j332tOmTcNut/Ppp5/m7gP0xRdfEBgYyKJFi7j55psB8PPz49NPP80tjFylVBdChXXbbbdx2223mR0jj40JiYyZvY2RPXrScmgPWDAK1n4Gm6bDjj/g6Q3gpyX+IiIu83+HLv6ecV67hhd2XeLc82aZPLvp8jNdxODBg9m8eTPLli277Gu0bNky99dne3MdO3Ys30Jo69athIaG5nkiEhYWRmBgIFu3bs0thOrWrZtbBP3zc/+cgtK+ffs8146Li2PXrl1UrFgxz/GMjAx2796d+zo8PNzlRRCU8kKoWrVqWK3WPM8sAY4ePZpvNexOZsQcZOWek8yIOUjL25pDz/eg1f3w2/NQo5mKIBERV/PyM//cQhgyZAi//vorS5YsoXbt2pd9HU9Pz9xfnx2JsdvtV5TNz6/ov9eUlBRat27N119/fcF7/yyqLufal8PtJ0tfipeXF61bt2b+/Pm5x+x2O/Pnz7+gAi2q6OhowsLCcqve4pBwOo1NCUlsPpjE9PXxAPwSd4jNB5PYlJBEQoWm8Mif0P3tcx86tRd+eQZS83neLCIiZZbD4WDIkCHMnDmTBQsWUL9+/QI/4+Xlhc1mu+J7N2vWjPj4eOLj43OPbdmyhcTERMLCwgr83OHDh3OPrVq1Ks85V199NTt37qRGjRo0atQoz1dAQMAVZy8qty+EUlJSiI2Nzd0Iau/evcTGxnLgwAEAhg0bxqRJk5gyZQpbt27lySefJDU1lYEDB17RfQcPHsyWLVtYu3btlf4WcnV8ayG9Ji7j1gnLSM10/h/1ZGoWt05YRq+Jy+j41kKwWPL+a2LOCFg/GSZcDes+hyus3kVEpHQYPHgwX331Fd988w0VK1bkyJEjHDlyhPT09It+pl69emzcuJHt27dz4sSJy25G2rVrV8LDw+nfvz8xMTGsWbOGBx98kOuvv542bdpc8nNNmjRhwIABxMXFsXTpUl566aU85/Tv359q1arRu3dvli5dyt69e1m0aBFPP/00CQklv6WM2xdC69ato1WrVrRq1QpwFj6tWrXiP//5DwD9+vXj3Xff5T//+Q+RkZHExsYyZ86cCyZQu4Nx/SLxsFy8QVxkaCDr95/C4XCcO9hxGAS1gIxE+PU5+PRGOLTB9WFFRMRUH374IUlJSXTu3JlatWrlfk2bNu2inxk0aBBNmzalTZs2VK9eneXLl1/WvQ3D4KeffqJy5cp06tSJrl270qBBg0veG5wrt2fOnEl6ejpt27bl0UcfZfTo0XnO8fX1ZcmSJdSpU4fbb7+dZs2a8cgjj5CRkUGlSpUuK++VMBx5furK+ZKTkwkICCApKalY/oA2H0zi1gmXnuzWIqQSA9rXo1dEMD6eVrDlwNpJsGA0ZJ0BDGjzMNz4MlSofMWZRETKsoyMDPbu3ZtnrxspGy71Z1vYn99uPyJUVv09Ty33f8ff04p+bULx9rCw+WAyL/ywkWvHLuCduds4nJIN1zwJQ9c5l9fjgHWfwbovTMsvIiJSFqgQughXTJYGqOrvRXV/b8JDAhjdtwXhIQFU9/cmqn5l3rqzJatG3si/brmKkMAKnErNInrhbjq+tZDBX8ew5oQXjtsnwYM/Q5Nb4Jqnzl3YdnnPgUVERMozPRorQHE/GgPIzLHhZbVgGAYOh4Msmx1vj7x7VuTY7Py59RiTV+xl1Z5Tuceb1arEwGvrcVvk34/NwPno7LOuUOda6DISvPPuzSAiUp7p0VjZVRyPxlQIFcAVhVBRbTuSzJQV+5i54SAZ2c5VY5V9PbmnbR3uv6YuIUcXw7f9nCdXrAXdRkPz2889dxMRKcdUCJVdKoRKgDsUQmclpmXx/bp4pqzYz8FE5/JJiwE3h9Xkmbr7uWrD6xin9zpPbtAZerwL1RqbF1hExA2c/WFZr149KlSoYHYcKUbp6ens27dPk6XLi0BfLx7r1JAlL3bhkwda06FRVewOmPPXEbrP9uY2+7tsbPwUDqs37FkE/2sP819XZ3sRKdfO7qiclnaRpqlSap39M/3nrtlFVapbbJRXVovBzc1rcnPzmuw4eoYpK/YxI+Ygm45mctvRjjT3acYHAd/SKGkF7FsOXfTHLCLll9VqJTAwkGPHjgHOfWwMTR0o1RwOB2lpaRw7dozAwECsVmvBH7oIPRq7iOjoaKKjo7HZbOzYscMtHo1dSlJaNtPXxzN15X4OnEoDHHSzriO4fhg3db6B9g2rYmSlQOoJqFLwNu0iImWJw+HgyJEjJCYmmh1FilFgYCA1a9bMt7DVHKFi4k5zhArDZnewcNsxpqzcx9KdJ3KPNwnyZ1yVH2gW/z1Gx2HQ4Rnw9GFjQiJjZm9jZI+raFk70LzgIiIlwGazXXbbCXEvnp6elxwJKuzPbz0zKWOsFoOuYUF0DQti17EzTFmxnx9jEth5NJnjp2IJs2bAojfJ3vAtnre+y4yttVi55yQzYg6qEBKRMs9qtV7RYxQpezQiVIDSNiKUn+SMbH5Yl8CUFXtpmTiff3t+RZCRCMCftOXljAfI8qvFlIfb4nBAZT9Palf2NTe0iIjIFdCjsWJSFgqhs+x2B4t3HGfo5MU86/EjD1nn4mHYSXV4MzR7KAvsV+eeu29sTxOTioiIXBktn5cLWCwGXa6qwah+1zLG/iC3Zr3JWnsT7FjYZHdOoLYY8N+7I0xOKiIiUjI0R6gc6tMqhEY1/Ll1goO7s/5DA+Mwx3F2sbc74Oif41luPMi1kc21xFRERMo0jQhdhKuarrodw8JuRwhny50bPTfxRNrHhM/qymfv/Ytl24+gp6ciIlJWaY5QAcrSHKF/OpyUzm0TllMr0Id+UaFMWxvP4cQMpvX2w2fOMIJTtwDwl70u31R/htt69qFdg6ompxYRESkcTZYuJmW1EALIzLHhZbVgGAYOh4Msmx1vDyvYbSSv+AzPha9TwXYGgGk5nVkYOpjHukdxdZ3KJicXERG5NBVCxaQsF0IFSj1B2uyX8P3rO8A5OtQz601uuCqIYTc1oUVIgMkBRURE8qdVY3Ll/Krhe9fH8PBcsqqFsa7eIKwWCwu2HePWCct4/Mt1bDuSbHZKERGRy6YRoQKU6xGhf7LbwLCw92Qa4+fvxGvjVzQ39vGe7S6uC2/Ms12b0KiGv9kpRUREAD0aKzYqhPKRmYLtvy2wZpzmhKMSY7LvY6bjOvpE1uaZro2pW9XP7IQiIlLO6dGYuI63P9Z+U6FaU6oZybzn9RHfeb7OX7ErueG9xYz4cSMJp9PMTikiIlIgjQhdRHR0NNHR0dhsNnbs2KERofzkZMGqaFj8NmSnYcPC5zm3MC7nDrKsvtwTVYfBXRpRM8DH7KQiIlLO6NFYMdGjsUJIjIe5I2HrLzgweLn6B3wVXw0ALw8L97ery5OdG1K9orfJQUVEpLxQIVRMVAgVwc4/4dAGuP4FVu4+yfvztrN9XwLJ+FHB08qD19bliU4NqeznlfuRjQmJjJm9jZE9rqJl7UDzsouISJmiOUJS8hp3hetfAKB9w6p8f2cNYvyf5Z3Ks3Bkp/Hx4j1c9/ZC3v9jO0np2QDMiDnIyj0nmRFz0MzkIiJSTmlEqAAaEboCC8fA4rEApPuG8I4xkM9PNgMM/Lys9G0VwuzNRziVmkVVPy+mPNwWhwMq+3lSu7KvudlFRKRU06OxYqJC6Ao4HLDtN5gzApLiAThW83ru2N+HeEfQJT+6b2zPkkgoIiJllB6NifkMA5rdCoNXQ8fnwOJJjSOLWVRhBI97/JLvRzwsBuP6RZZsThERKbdUCInreflB11fhyRVQvxNWeyYD24fme+r1TavTsXG1ks0nIiLllgohKTnVm8CDP0O/rzkZ/hjgHDS6yjhALU4CMH/rMW54dxFfrtqPza6ntiIi4loqhKRk/f24rEqAP9X9vYkM9uO7qpNY4PMCz/jM5qrqPiRn5PDyrM30/d9yNiYkmp1YRETKME2WvgjtLO16mTk2vNKPY3w/AOJXAeCodhVz6g7nxXWVOJOZg2FA/3Z1eOHmqwjw9TQ5sYiIlBZaNVZMtGqsBNjtEPctzPsPpJ0AIKPZHYzOvo8vN2cCUNXPi5E9mnHH1SEYhmFmWhERKQVUCBUTFUIlKP00zH8D1n0OOMC7EjG3zODFhWnsOpYCQFS9yrzRpwVX1dSfhYiIXJwKoWKiQsgEB9fDb8+Dpx889CtZNgefL9/LB3/uJD3bhtViMPDaejx7UxP8vT3MTisiIm5IhVAxUSFkErsN0hPBr6rzdUYSKXNH8crp7vy4LR2AoErevHxrGD3Da+lxmYiI5KENFaV0s1jPFUEAC9/Ef8MnvHf0YeZct5t6VXw4mpzJkG828MBna9h9PMW8rCIiUmqpEJLSoXlfCGoB6ae5au3LLAgczehrbHh5WFi26wS3jFvCu3O3k55lMzupiIiUIno0VgA9GnMjthxYOwkWjIasM2BYSG7xICMSezN7p/NxWUhgBV67rTldwy7dy0xERMo2PRqTssfqAdc8CUPWQos7wWGn0qbJRNeczUf3tyY4wIeDiek8OnUdj05ZS/ypNLMTi4iIm1MhJKVPpVpw52fOdh2122Jc/yK3tKjJn89fzxOdGuBhMfhz6zG6vr+YiQt2kpnjfFy2MSGRez9Zpd2qRUQklwohKb0aXA+PzgP/GgD4enkwImUMa6MW0bmeD5k5dt79Ywfdxy1l6c7jzIg5yMo9J5kRc9Dk4CIi4i40R6gAmiNUihxcD5NuAMBRsRbrrhrOY+vqcDo9GwAvq4Usm52qfl5MebgtDgdU9vOkdmVfM1OLiIgLaB+hK6ReY6XUznkw+wU4vReApbYWvJLzEHscwRf9yL6xPUsqnYiIlBAVQsVEI0KlUHYGLB8HS98HWyZZDiuf2G5lQk5fMvHKPc3DYvDuXRH0aRViXlYREXEJrRqT8svTBzqPgMGroNFNeBk2eltWYJC35u/eoia3tqxlUkgREXEHatQkZVeVBtB/OvtXTOel3/aQaXiDAyzYqcVJftkIienZTLzvagIqeJqdVkRETKARISnbDAOvFr3Y7tuG8JAARvdtwQtVlzPfZzjPec5k9c7D9P3fcvaoRYeISLmkOUIF0ByhsiEzx4aX1YJhGDimPYCx9WcA4o1avJT5ILFerflf/9Z0bFzN5KQiIlIcNEdI5B+8Pay5HeqNu6fCHZ+Bf01CHYeZ6vUWY2zv8q8vfmfqyn3mBhURkRKlQkjKH8OA8DudrTqueQqHYaWndQ1/eD7P2l8m8e9Zm8i22c1OKSIiJUCFkJRfPpXgljEYjy/BEXoNFYxs9hDMV6sO8OBnazidmmV2QhERcTEVQiI1W2AM/B3Lo/N47v478fOysnLPST4Y/xZ79+02O52IiLiQCiERAIsFareha1gQPz51LR0CTvBSxn+p9kUHdv7yLthyzE4oIiIuoEJI5DxX1axE9D0R7PdqSEUjncbr3+Dkf6/FEb/G7GgiIlLMVAiJ5COwfiR1XljBjODnSXT4UTVlO8ZnN2GbNQRST5odT0REiokKIZGL8PLypO+gl/nt+l+YbrseAGvsl9g+6gQ5mkgtIlIWqBASuQTDMOh/Q2uq3f8pDzheZ6u9Dh+l38i2ExlmRxMRkWKgnaULoJ2l5axdx87w2OTV7D+Vho+XFx/c04quPts4se5H/nXqNp65tQ0taweaHVNERNDO0iLFrlGNivw4uBPtGtYgNcvG41+u4fQPT1NtyxTGHh7I7nmfgv5dISJSqqgQEimCyn5eTHm4Lb0jg7E5LAxOvI/djmCqG8n03f8GqR93Y8fGNSScTjM7qoiIFIIKoYuIjo4mLCyMqKgos6OIm/G0Wvgp9hAAK+wtuCVzLG9l30Oawxu/I6tp8GM3Zr/3KGSeMTmpiIgURHOECqA5QpKfWRsOMnx6HDn2c98+wZzgP55fcot1rfPA3VMhrLdJCUVEyjfNERJxoT6tQpg1uEOeY4eoxhPZz7Hzpi/g6gHQ7LZzb2ZrlZmIiDtSISRyhQwj7+uRG4PI6jHu3Btpp2DC1TD/DcjS3CEREXeiQkjkMlX196K6vzfhIQGM7tuCBtX9AFi3P5FnvttAjs3uPHHTD5B8EJa+C/9rB9t/NzG1iIj8k+YIFUBzhORSMnNseFktGIaBw+FgwbZjPPlVDFk2O30ig3nv7kisBrD1F5gzEpITnB9s0h26j4XK9cyMLyJSZmmOkEgJ8PawYvz9CMwwDG5sFkR0/6vxsBjMij3ESzM34QAIuw2GrIEOz4LFA3b8DtHtYMm72ntIRMREKoREitlNYUH8t18kFgO+WxvPa79sweFwgJcf3PQaPLkC6l0HORlwet+Fk4xERKTEeJgdQKQs6hURTGaOneHT45i8Yh8+nlb+dUtT5+hR9aYw4BfY/CM06HzuQ0kJgAEBIWbFFhEpdzQiJOIid7auzag+LQD4aPFuJizYde5Nw4DwO8GvmvO1wwG/PgcTo2D5eLBlm5BYRKT8USEk4kL3X1OXf/dsBsD783Ywacme/E/MSoGMJMhOhXkvw0fXwb5lJZhURKR8UiEk4mKPXteA4Tc3AWD07K18uXLfhSd5V4SBc6B3NPhWheNbYXJPmPEYnDlasoFFRMoRFUIiJWDIDY0Z3KUhAC//9BfT18VfeJLFAq3uhyHroM3DgAEbp8HENpCwvmQDi4iUEyqERErI8Jub8nCH+gD868eN/Bx3KP8TfavArf+FQfOhVqRzhCioeckFFREpR1QIiZQQwzB4+dZm3NeuDnYHPDctlrl/HWFjQiL3frKKjQmJeT8Q0hoGLYABP4Onj/OYLQcWjnG27RARkSumQkikBBmGwajeLbi9VQg2u4Oh32xg/PydrNxzkhkxBy/8gMUKgXXOvV47CRaPdfYuWz8F7PaSCy8iUgapEBIpYRaLwTM3NqZjo2pk2ez8ufUYAL/EHWLzwSQ2JSSRcPoizVmDr4YazSH9NPzyNHx+MxyOK8H0IiJli3qNFUC9xsQV6o34rcBz9o3tmf8bthxY87HzEVnWGTAsEPUodHkJKgQWb1ARkVJKvcZE3Ni4fpF4WPJvreFhMRjXL/LiH7Z6QPvBMGQttLgDHHZY8wn8NNg1YUVEyjAVQiIm6NMqhFmDO+T73lOdG9KnVSHabFSqBXd+Dg/+BDXCoMv/FXNKEZGyT73GRExmGHkb0I9fsAt/Hw8e69SwcBdo0NnZyPWfzVvnvQL2HOg8wrlZo4iI5EsjQiImqervRXV/b8JDAhjdtwUtQwKo4GkF4M3Z23h7zjYKPYXvn0XQ6f2wYgKsnAgT28JfM/NWWiIikkuTpQugydLiSpk5NrysFgzDwOFwkGWz8/myfbw1ZxsA/dvV4fXeLbBeZD7RRe34A35/AU7vc75u0AV6vAvVGhXvb0BExE1psrRIKeDtYcX4ezTHMAy8Paw82bkhb/YNxzDg69UHeOa7DWTlFHG/oCY3w1Or4PoRYPWGPQvhw/Yw/w3IusjSfBGRckiFkIgbuq9dHSbc2wpPq8GvGw/z2JfrSM+yFe0inhWgy0h4aiU06gq2LFj9MWQmuya0iEgpVC4Kob59+1K5cmXuvPNOs6OIFNqtLYOZ9GAbfDwtLNp+nAc+W01SenbRL1S1IfT/Afp9BT3ehoo1z72Xcrz4AouIlELlohB65plnmDp1qtkxRIqsc9MafPVIOyr5eLBu/2nu+WQVx89kFv1ChgHNekHkfeeO7fwTxoXDkncg5zKuKSJSBpSLQqhz585UrKglxFI6talXhWmPt6eavzdbDydz10criD9VDPN8tsyEnHRYMAo+vBZ2L7jya4qIlDKmF0JLliyhV69eBAcHYxgGs2bNuuCc6Oho6tWrh4+PD+3atWPNmjUlH1TERM1qVeKHJ9pTu3IF9p1M486PVrDz6Jkru+htE+H2T8E/CE7ugi/7wvcDICmf5q8iImWU6YVQamoqERERREdH5/v+tGnTGDZsGK+88goxMTFERETQrVs3jh07lntOZGQkLVq0uODr0KFDRc6TmZlJcnJyni8Rd1Cvmh8/PHEtjWv4czQ5k7s+XklsfCIAGxMSufeTVWxMSCz8BQ0DWt7lbNXR7klnz7Its2BiFMR86YrfgoiI2zG9EOrevTujRo2ib9+++b7//vvvM2jQIAYOHEhYWBgfffQRvr6+fP7557nnxMbGsnnz5gu+goODi5xnzJgxBAQE5H6FhoZe9u9NpLjVDPDh+8fbExEaSGJaNv0nrWLFrhPMiDnIyj0nmRFzGaM5PgHQfSw8vgRC20F2qpq3iki5YXohdClZWVmsX7+erl275h6zWCx07dqVlStXuuSeI0eOJCkpKfcrPj7eJfcRuVyV/bz4+tF2XF0nkNQsGw9+voYfYxIA+CXuEJsPJrEpIYmE00WcR1QzHAbOgf4/wlW3nju+dymkHLv450RESjG37jV24sQJbDYbQUFBeY4HBQWxbdu2Ql+na9euxMXFkZqaSu3atZk+fTrt27fP91xvb2+8vb2vKLeIq/l7exBzIBGAHLuDMxk5AJxKzeLWCctyz9s3tmfRLmyxQONz//Ag9SR8/wDY7XDDvyHqEbBYrzS+iIjbcOsRoeLy559/cvz4cdLS0khISLhoESRSmozrF4nHea03zvbL8bAYjOsXeeU3ST8NgXUhM8nZsuOTzhC/9sqvKyLiJty6EKpWrRpWq5WjR4/mOX706FFq1qx5kU+JlA99WoUwa3CHfN+bNbgDfVqFXPlNqjWCQQug53vOuURHNsJnXeHnoZB26sqvLyJiMrcuhLy8vGjdujXz58/PPWa325k/f77LR3Wio6MJCwsjKirKpfcRKQ7GeT1ZZ28+XHwXt1gh6lEYsh4i/t6QMWaqc3VZemLx3UdExASmF0IpKSnExsYSGxsLwN69e4mNjeXAgQMADBs2jEmTJjFlyhS2bt3Kk08+SWpqKgMHDnRprsGDB7NlyxbWrtVjAHFfVf29qO7vTXhIAKP7tKCqvxcA/1u4my9X7ivem/lXh74fwsDfoUYYhN2m1WUiUuoZDofDUfBprrNo0SK6dOlywfEBAwYwefJkACZOnMg777zDkSNHiIyMZPz48bRr165E8iUnJxMQEEBSUhKVKlUqkXuKFEVmjg0vqwXDMLDb7bzx61a+WLEPgNdua86Aa+sV/01t2c62HN7+ztcndsHaSdB5pIojEXELhf35fVmF0NKlS/n444/ZvXs3P/zwAyEhIXz55ZfUr1+fjh07XlFwd6NCSEobh8PB2Dnb+HjxHgD+c2sYD3es79qbftnX2aLDrwbcPApa3n3h8zoRkRJU2J/fRX409uOPP9KtWzcqVKjAhg0byMx0NmtMSkrizTffvPzEIlIsDMNgxC1X8VTnhgC8/usWPl26x7U37fAMVG0Mqcdg5mMw+VY4ttW19xQRKQZFLoRGjRrFRx99xKRJk/D09Mw93qFDB2JiYoo1nJk0WVpKM8MweKFbU4be0AiAUb9t5ePFu113wwad4cnlcON/wKMC7F8GH3WEP/4NmSmuu6+IyBUqciG0fft2OnXqdMHxgIAAEhMTiyOTW9BkaSntDMNg2E1NeObGxgCM+X0b0Qt3ue6GHt5w3fMweLVzZ2p7DqyYADFTXHdPEZErVORCqGbNmuzadeFfpsuWLaNBgwbFEkpEiodhGDx3UxOG3dQEgHfmbmfC/J2uvWnlunDP13Df99DoJogadO49u8219xYRKaIiF0KDBg3imWeeYfXq1RiGwaFDh/j6668ZPnw4Tz75pCsyisgVevrGxrzQrSkA783bwbg/d1xex/qiaNIN7v8BPJxL+rFlw6QusGAUZBWxD5qIiIsUudfYiBEjsNvt3HjjjaSlpdGpUye8vb0ZPnw4Q4cOdUVGESkGg7s0wmIYvDVnG+P+3MnCbceIS0hiRsxBWtYOdH2Arb/A4Tjn18Zp0P1taNrd9fcVEbmEy95HKCsri127dpGSkkJYWBj+/v7Fnc0taPm8lCUJp9P4bNlevli+L/dYVT8vpjzcFocDKvt5Uruyr2tu7nA4i6E5IyD5oPNYk+7QfSxUrueae4pIueWyfYQefvhhPvjgAypWrJjneGpqKkOHDuXzzz+/vMRuJjo6mujoaGw2Gzt27FAhJGVCvRG/FXhOkTvWF1VmCix5G1ZGOydUe/hAp+HQ4TmwFnmQWkQkXy4rhKxWK4cPH6ZGjRp5jp84cYKaNWuSk5NzeYndlEaEpCyZteEgw6fHkWO/8NveYsD7d0cWT7PWwji2DWYPh31Lod51MOAXbcIoIsWmsD+/C/3Pr+TkZBwOBw6HgzNnzuDj45P7ns1mY/bs2RcURyLiXvq0CqFRDX9unbDsgvfsDthyOJke4bXw8iiBNoQ1rnIWP5t+gJrh54qgjCTnqFFACRVkIlKuFboQCgwMxDAMDMOgSZMmF7xvGAavvfZasYYTEdcxDOe0nbP/C/DJkj2s2nOS8fe0ol41v5IJ0fKuvMcWjIYNX0HnEXDNk2D1zP+zIiLFoNCF0MKFC3E4HNxwww38+OOPVKlSJfc9Ly8v6tatS3BwsEtCikjxOduxvlagD/2iQpm2Np7DiRk8f3MTxvy+jY0JSfQcv5TRfcNL7jHZWbYcOLIJslNh3ssQ+w30fA/qdSjZHCJSbhR5jtD+/fsJDQ3FYimBoXM3oDlCUhb9s2O9w+Egy2bH28PKocR0np0Wy5q9pwC4/eoQXu/dAn/vEpzEbLdD7Nfw5yuQdtJ5rOU9cPMb4K/H7yJSOC7tPg+QlpbGgQMHyMrKynO8ZcuWl3M5t6VCSMobm93BxAW7+GD+DuwOqF/Nj/H3tCK8dkDJBkk7BfNfh/WTAQd4B8DdU6Bhl5LNISKlkssKoePHjzNw4EB+//33fN+32crGFvpaPi/l3Zq9p3j2uw0cSsrA02rwr1uu4uEO9bFYSnhlV8J6+O05OL0PhqzTqJCIFEphC6EiP9969tlnSUxMZPXq1VSoUIE5c+YwZcoUGjduzM8//3xFod2Jmq5Kede2fhVmP3Md3ZoHkW1zMOq3rTw8ZS0nUjJLNkjt1jBoIQyck7cIWjPJOWokInIFilwILViwgPfff582bdpgsVioW7cu999/P2+//TZjxoxxRUYRMUmgrxcf3d+aUX1a4O1hYdH243T/YCnLdp4AcH2/srMsVggKO/d6+xznHkQTWsP6Kc55RSIil6HIhVBqamrufkGVK1fm+PHjAISHhxMTE1O86UTEdIZhcP81dfl5SEca1/Dn+JlMHvh8NWN/38b09Qms3HOSGTEHSzaUb1Wo0RzST8EvT8PnN8PhjSWbQUTKhCIXQk2bNmX79u0ARERE8PHHH3Pw4EE++ugjatWqVewBRcQ9NK1ZkZ+HdKR3ZDAOB3y0eDffrj4AwC9xh9h8MIlNCUkknC6BzvKhUfD4Yuj2Jnj5Q8Ja+OR6mP2ic0NGEZFCKvJk6a+++oqcnBweeugh1q9fzy233MKpU6fw8vJi8uTJ9OvXz1VZTaFVYyJ5uUW/sn9KPgRzX4K/Zjhf124Lj84rufuLiFty+fL5s9LS0ti2bRt16tShWrVqV3Ipt6RCSCSvS/Ur87AYvHtXRMlvxAiweyHMfgFuHgVNbyn5+4uIWymxQqisUyEkcqHNB5Py7Vd2T1QoY24PxzCreaotO29LjnVfwKndcP0I8PY3J5OImKLYm66e5XA4+OGHH1i4cCHHjh3Dft5qjRkzZhQ9rRv65z5CIpK/f/YpA/hubTwZ2TbeurMl3h7Wkg/0zyIoPdG5O3VGEmz6EW55E8L6qMO9iORxWfsIPfDAA+zduxd/f38CAgLyfJUV2kdI5OLO9isLDwlgdN8WtKwdgL+3BxYDZsUe4oHP1pCYllXwhVypQiD0/QQC68KZQzD9Ifjqdji529xcIuJWivxorEqVKnz11Vf06NHDVZncih6NieQvv35la/ae4qmvYjiTmUODan58MTCKulVLoIv9pWSnw7L/Or9sWWD1gg7PwHXPg2cFc7OJiMu4bGfpgIAAGjRocEXhRKT08/aw5s4FMgwDbw8r1zWuzg9PXktIYAX2nEil7/9WsH6/ybs/e1aALv8HT62Chjc6i6Gl78GJnebmEhG3UORC6NVXX+W1114jPT3dFXlEpJRrWrMiM5+6lvCQAE6lZnHvpNX8uvGQ2bGgakO4/0e4e6qzMKr1jwbRmWfMyyUipiryo7H09HT69u3L8uXLqVevHp6ennneL2u7S+vRmMjlScvK4elvY/lz61EAXrylKU9e39C8FWUXc2wbfHYTdHgarn0aPLzNTiQixcBlq8YGDBjA+vXruf/++wkKCnK/v9RExC34ennw8QOtGfXbFr5Yvo+352znwMk03ujTAk9rkQejXSf2a8hMhgWjIO476PEONLzB7FQiUkKKPCLk5+fH3Llz6dixo6syuRWNCIlcucnL9/L6r1uwO+C6xtWI7n81lXw8C/5gSXA4YNN05+7Uqcecx8L6ONt3BJiwMaSIFAuXTZYODQ1VQSAiRfJQh/p88kAbKnhaWbrzBHd9uJKDiekl173+UgwDWt4NQ9dBuyfAsMCWWTAxCtZMMi+XiJSIIhdC7733Hi+++CL79u1zQRz3ER0dTVhYGFFRUWZHESkTuoYFMf2J9tSo6M32o2foE72cT5bsMad7fX58AqD7W/DYYme/suxU59J7ESnTivxorHLlyqSlpZGTk4Ovr+8Fk6VPnTJ5qWwx06MxkeK1fv8pnv8+jn0nz3Wpr+rnxZSH2+JwQGU/T2pX9jUxIWC3O5u4hvU+t1v1kc3gX8P5JSJuz2WTpceNG3cluUSknLvjw5UXHDuZmpWnd1mJdq/Pj8UC4Xeee52T5dyZOuUY3PgytHkYLCa0EBGRYndZq8ZERC7XuH6RBXavdzupx8DLFzKTYPZw2PAl9HwfarcxO5mIXKFCzRFKTk7O8+tLfYmIXEqfViHMGtwh3/da1QnkprCgEk5UCAG1YdBC6PEueAfA4Tj4tCv88gykla3pACLlTaEKocqVK3PsmHNZaWBgIJUrV77g6+xxEZHCOrsN2dndyNbuO83dH6/kaHKGaZkuymKFtoOcq8si7gUcsH4yTGgNiQfMTicil6lQj8YWLFhAlSpVAFi4cKFLA4lI2Xe2e32tQB/6RYUybW08B06lgQP+OpRMn+jlfDYgirBgN1yg4F8D+n4EVz8Ivz0P/kEQEGp2KhG5TEVeNXbgwAFCQ0Mv2FHa4XAQHx9PnTp1ijWg2bRqTMQ18utefyw5k4e+WMPu46n4eVmZ2P9qujR141VatmzISAK/as7Xaadg+ThnZ3ufAFOjiZR3LttQsX79+hw/fvyC46dOnaJ+/fpFvZyIlFP5da8PreLLjCc70L5BVVKzbDw6ZR1frtpvctJLsHqeK4IAFrwByz+ACW1g4/fOXatFxK0VuRByOBz59hdLSUnBx8enWEKJSPkV4OvJlIfbcmfr2tjsDl6etZlRv27Bls8qM7fT7Dao2si5ymzGIJjSy9nUVUTcVqEfjQ0bNgyADz74gEGDBuHre27DM5vNxurVq7FarSxfvtw1SU2iR2Mi5nA4HPxv0W7embsdgG7NgxjXrxUVvNx8/56cTFgxAZa8CznpYPGAa56C6/8F3v5mpxMpN4p9Q8UNGzYAzr+cNm3ahJeXV+57Xl5eREREMHz48CuILCJyjmEYDO7SiNqVK/DC9I3M/eso93yykkkD2lCjohuPPnt4Q6fhEH4XzBkJ23+DFePBYYduo81OJyLnKfJk6YEDB/LBBx+Um9ERjQiJmG/tvlM8NnUdp9OyCQmswOcPRdG0ZkWzYxXO9jmweCz0/xH8qjqPORzn9g4QEZdw2WTpL774Is8Fk5OTmTVrFtu2la3n4Gq6KuI+oupVYeZTHahfzY+Dienc+eEKlu50Ltpwiw72l9L0FudmjP8sgqbdDwtGq6mriBso8ojQ3XffTadOnRgyZAjp6elERESwb98+HA4H3333HXfccYersppCI0Ii7uN0ahaPf7meNftO4WExGNWnBduOnGHyin08dG09Xr2tudkRC7ZvOUzu4fx1YF3o/razWBKRYuWyEaElS5Zw3XXXATBz5kwcDgeJiYmMHz+eUaNGXX5iEZECVPbz4stH23JzWA1y7A5GzNjEtLXxAPwSd4jNB5PYlJBEwum0Aq5korrXwt1ToVIIJO6Hb/vBt/dpd2oRkxR5RKhChQrs2LGD0NBQHnzwQYKDgxk7diwHDhwgLCyMlJQUV2U1hUaERNxPvRG/FXiO6R3sC5KZAovfglX/A3sOeFRwTrK+9mnw8Cr48yJySS4bEQoNDWXlypWkpqYyZ84cbr75ZgBOnz6tfYREpESM6xeJhyX/ycYeFoNx/SJLNtDl8PaHm9+AJ5ZB3Y7OpfaxXztXl4lIiSn08vmznn32Wfr374+/vz916tShc+fOgPORWXh4eHHnExG5QJ9WITSq4c+tE5Zd8N7Uh9tybaNq+XzKTdVoBg/9CpumO/uWef79D0pbDqQeh0q1zM0nUsYVeUToqaeeYuXKlXz++ecsX74ci8V5iQYNGmiOkIiUuPNXob/440YOJZay1ViGAS3vhgbXnzu2dhJMbAMrJjp7momISxS5EAJo06YNPXv25ODBg+Tk5ADQs2dPOnToUKzhREQu5mwH+/CQAEb3bUGTIH8sBiScdi6v33WsFM9XdDhg15+QlQJ/vAQfd4L9K8xOJVImFXmydFpaGkOHDmXKlCkA7NixgwYNGjB06FBCQkIYMWKES4KaRZOlRdzX+R3s951M45Epa9lzPJXKvp58MbAtkaGBZse8PHY7xH4F816B9FPOYxH3wk2vg38Nc7OJlAIumyw9cuRI4uLiWLRoUZ7J0V27dmXatGmXl1ZE5DKc38G+fjU/fnjiWiJqB3A6LZv7Jq3K3Xix1LFY4OoHYeh6aP0QYEDct87O9lt/NTudSJlR5EJo1qxZTJw4kY4dO+bpQt+8eXN2795drOFERIqqip8X3wy6husaVyMty8bDk9fyS9whs2NdPt8q0OsDeHQ+1IpwPi6rUt/sVCJlRpELoePHj1OjxoXDsqmpqXkKIxERs/h5e/DpgDb0bFmLbJuDp7/bwJcr95kd68rUbu1s1TFwNgT9Ywftv2ZB2inTYomUdkUuhNq0acNvv53bzOxs8fPpp5/Svn374ksmInIFvD2sjL+nFQ9cUxeHA17+6S/+O28HcfGn3bs32aVYrFDnmnOvj/4FPzwME1pDzFTnvCIRKZIi7yP05ptv0r17d7Zs2UJOTg4ffPABW7ZsYcWKFSxevNgVGUVELovVYvB67+ZU9fdi3J87+WD+Tub8dYTtR84wI+YgLWsHmh3xytiyoFoTOL4Vfh4KMV9Cz/egVkuzk4mUGkUeEerYsSOxsbHk5OQQHh7OH3/8QY0aNVi5ciWtW7d2RUYRkctmGAZ3tq7NE9c3AGD7kTNAKepNdinBreCJpXDzKPDyh4Q18Mn1MPtFyEgyO51IqVDk5fPljZbPi5R+ZaI3WUGSD8Hc/4O/ZjpfV20Eg9c4H6eJlEMuWz7/Tz179uTw4cNXcgkREZe7VG8yi0Hp6E1WkErBcNdkeGCmswiKelRFkEghXFEhtGTJEtLTS9lW9iJS7vRpFcKswfnvfG93wI6jZ8ixlZGJxg1vgCdXQNSgc8d2zYd5/3F2vBeRPK6oECrLoqOjCQsLIyoqyuwoIlKMzu7y8c/dPv63aDf9P13NsTMZ5oQqbh7eYP17LUxOFsweDss/gOi2sOUnZwsPEQGKWAjl5OQwdepUjh49CkDdunXx9PR0STCzDR48mC1btrB27Vqzo4hIMTi/N1l4SADV/b15o09z/LysrN57ip7jl7Fqz0mzoxYvDy/o9iYE1oHkg/D9g/D1nXBSG+CKwGVMlvb19WXr1q3UrVvXVZnciiZLi5Qd5/cmy7LZ8fawsvt4Ck99FcP2o2ewGDC8W1Oe6NQQy0XmFZVKWWmw7H3nyJAtC6ze0PFZ6PgceFYwO51IsXPZZOm2bdsSGxt7JdlERExxfm8ybw/nZOKG1f2ZNbgDt18dgt0Bb8/ZzqCp60hMyzIzbvHy8oUb/g1PrXLOI7JlwuK3YI/2f5PyrcgjQt9//z0jR47kueeeo3Xr1vj5+eV5v2XLsrWRl0aERMoPh8PBtLXx/Ofnv8jKsRMSWIEP77+69G+8eD6HwzlXaPd86DX+3ISpnEzn/CKRMqCwP7+LXAhZLBcOIp0dZjYMA5vNVvS0bkyFkEj589ehJJ76Oob9J9Pwslp4uVcY97erg2EYbExIZMzsbYzscVXZKpBST8LH10Gbh+HaoSqIpNQr7M/vIrfY2Lt37xUFExFxd82DA/hlaEdemB7H3L+O8vKszazde4oxt4czI+YgK/ecLBstOv5pw5fOydQL3oC476Dnu9Cgs9mpRFxOO0sXQCNCIuWXw+Hgs2V7GTN7KzYH1K5cgZSMHBLTs6nq58WUh9vicEBlP09qV/Y1O+6VcThg03SY+xKkHnMea347dBvt3KxRpJRx2aMxgN27dzNu3Di2bt0KQFhYGM888wwNGza8/MRuSoWQiOTXosMA/vmXZ6lv0XFWeiIsfBPWTgKH3dnD7IZ/wzVPmp1MpEhctmps7ty5hIWFsWbNGlq2bEnLli1ZvXo1zZs3Z968eVcUWkTEHY3rF4n1vKX0Z4sgD4tRNlp0nFUhEHq8DY8tgpA2kJUCJ3aYnUrEZYo8ItSqVSu6devG2LFj8xwfMWIEf/zxBzExMcUa0GwaERIRgM0Hk7h1wrILjv88uAMtQwNLPlBJsNsh7hto2gN8qziPJcaDhw/4Vzc3m0gBXDYitHXrVh555JELjj/88MNs2bKlqJcTESlVjPP2WHx77nYyssvWatlcFgu0uv9cEeRwwM9DYGJrWDMJ7GX09y3lSpELoerVq+e7oWJsbCw1atQojkwiIm7n/BYddao4d2NetusED36+hqS0bJMTloD005B2CjKSnP3LJt0ACevNTiVyRYr8aOz111/nv//9LyNGjODaa68FYPny5bz11lsMGzaMl19+2SVBzaJHYyJy1vktOpbsPM6QrzdwJjOHJkH+TB7YluDAMt6uwm6DdZ/D/DcgMwkwoPUAuPGVcyNHIm7AZavGHA4H48aN47333uPQoUMABAcH88ILL/D000/nbl9fVqgQEpFL2Xo4mYe+WMPR5ExqVvJh8sNRXFWzHPxdkXIM/ngZNn7nfO1bFR78CWqGm5tL5G8uKYRycnL45ptv6NatG0FBQZw5cwaAihUrXnliN6VCSEQKcjAxnYc+X8POYylU9PHgkwfa0L5hVbNjlYx9y+G358FhgyeWO7vdi7gBl40Iqfu8iMiFktKyGTR1HWv2ncLLauG9uyPoFVFONiK0ZUPyIahc99zr5eOg7WPgE2BqNCm/XNp9fsOGDVcUTkSkrAnw9WTqI23pEV6TLJudod9u4NOle8yOVTKsnueKIIBVH8KCUTAxCjZ+71xtJuKmitxr7KmnnuL5558nISGhXHSfFxEpLB9PKxPuvZoaFbcwecU+Rv22lSNJGfxfj2ZYLGVr/uQl1WoJVRrCqd0wYxDETIUe70KNq8xOJnIBdZ8vgB6NiUhRORwOPlmyhzG/bwPg1pa1eO/uCLYfOVM2O9fnJycTVoyHJe9BTjpYPKD9YOj0Inj7m51OygF1nxcRMYlhGDx+fUOCKvnwwg9x/LrxMCdSMqlfzb9sdq7Pj4c3dHoBwu+GOSNg+2xY/oFzLtEdn5qdTiRXkUaEsrOzueqqq/j1119p1qyZK3O5DY0IiciVmBVzkJEzN5GebcNqGNgcjrLXub4wtv8Of/wb+n2tR2RSIlwyIuTp6UlGRsYVhxMRKS+e/T4299e2v//deTI1K0/fsjLTuf5SmnaHxt2cbTvO+uPf4OkLHZ8DzzK+EaW4rSKvGhs8eDBvvfUWOTk5rsgjIlKmjOsXicdFJkqXuc71BflnEXRiJ6yMhsVvQXQ72DHXvFxSrhV5snTfvn2ZP38+/v7+hIeHX7BqbMaMGcUa0Gx6NCYiV+pinetvDgti/L2t8PG0mpDKZA4HbPkJ5oyEM84uBTTtCd3HQmAdc7NJmeCyfYQCAwO544476NatG8HBwQQEBOT5cjfx8fF07tyZsLAwWrZsyfTp082OJCLl1PkdiP7YcpS7PlpJwuk0cwKZyTCgeR8YsgauHepcVbb9N5jYFpa+BzlZZieUcqLII0KlzeHDhzl69CiRkZEcOXKE1q1bs2PHjgtGsi5GI0IicqUOJ6Vz24Tl1Ar0oV9UKNPWxrP/ZBo4HCRl5FDZ15Px97biusbVzY5qnmNbna069i937kY9ZD34l+P/HnLFXNZiA5w9xxYtWsTu3bu57777qFixIocOHaJSpUr4+7v3/hARERH8+uuvhIaGFup8FUIiUhzO71yfZbNz/EwmT30dw8aEJAwDht/clCevb1i+Nl/8J4fDuRM1Doi459zxtFPqbC9F5rJHY/v37yc8PJzevXszePBgjh8/DsBbb73F8OHDixx0yZIl9OrVi+DgYAzDYNasWRecEx0dTb169fDx8aFdu3asWbOmyPcBWL9+PTabrdBFkIhIcfH2sGL8/WzMMAy8PazUruzL94+35962oTgc8M7c7Tz+1XqSM7JNTmsSw4CIfnmLoO1zYFxL58RqmxbpSPErciH0zDPP0KZNG06fPk2FCueWO56dRF1UqampREREEB0dne/706ZNY9iwYbzyyivExMQQERFBt27dOHbsWO45kZGRtGjR4oKvQ4cO5Z5z6tQpHnzwQT755JNL5snMzCQ5OTnPl4iIq/h4Whlze0veuiMcLw8L87Yc5bYJy9h2RH/3ALDpe8g6A3P/Dz7uBPtXmp1IypgiPxqrWrUqK1asoGnTplSsWJG4uDgaNGjAvn37CAsLIy3t8if9GYbBzJkz6dOnT+6xdu3aERUVxcSJEwGw2+2EhoYydOhQRowYUajrZmZmctNNNzFo0CAeeOCBS5776quv8tprr11wXI/GRMTVNiYk8uRXMRxMTKeCp5Wxd4TTOzLE7Fjmstthw5fw5yuQftp5LOI+uOl1zSGSS3LZozG73Z5vP7GEhAQqVqxY1MtdUlZWFuvXr6dr1665xywWC127dmXlysL9q8DhcPDQQw9xww03FFgEAYwcOZKkpKTcr/j4+MvOLyJSFC1rB/LL0I5c17ga6dk2nvkulld//otsm93saOaxWKD1ABgaA1cPcB6L+wYmtoZNP5ibTcqEIhdCN998M+PGjct9bRgGKSkpvPLKK/To0aM4s3HixAlsNhtBQUF5jgcFBXHkyJFCXWP58uVMmzaNWbNmERkZSWRkJJs2bbro+d7e3lSqVCnPl4hISani58XkgW0Z0qURAJNX7OPeT1ZxLNm5q//GhETu/WQVGxMSTUxpAt8qcNt4eORPqNkSMpLAu3j/8S3lU5Gbrr733nt069aNsLAwMjIyuO+++9i5cyfVqlXj22+/dUXGK9KxY0fs9nL8rykRKXWsFoPh3ZoSERrIsGmxrNt/mp4TlhF939XM3nS4/DRuzU9oFDy2CHb+AU26nTu+dykENdfqMimyIhdCtWvXJi4ujmnTphEXF0dKSgqPPPII/fv3zzN5ujhUq1YNq9XK0aNH8xw/evQoNWvWLNZ7iYi4m5vCgvh5aEcembyWPSdSufeTlXj/vQv1L3GHuLN17fLVuPUsi9XZu+yslGPwXX+wejjnDkXcl7edh8glFLkQAvDw8KB///7079+/uPPk4eXlRevWrZk/f37uBGq73c78+fMZMmSIS+8dHR1NdHR0vvOhRERKSv1qfuw5kQqAzQFpWc6/k8pl49aLST0BlYLh+Fb4aTDETIWe70HNcLOTSSlgesmckpJCbGwssbGxAOzdu5fY2FgOHDgAwLBhw5g0aRJTpkxh69atPPnkk6SmpjJw4ECX5ho8eDBbtmxh7dq1Lr2PiEhB1Li1AEFh8MRSuOkN8PSD+NXOpfa/j4AMbUMgl2Z6i41FixbRpUuXC44PGDCAyZMnAzBx4kTeeecdjhw5QmRkJOPHj6ddu3Ylkk87S4uIO7hY49a37ginX5SalOZKOujcc2jLLOfrSiEwZC14Fa6tkpQdLm2xUZ6oEBIRd3C2EDIMZyeKszws8EafcO5tq2Ioj13zYfYL0Kgr9Hjb7DRiApftIyQiIiWvqr8X1f29CQ8JYHTfFrQIroS3h4UcO4ycsYmXZ20m22Yvv8vrz9foRnhqJdz4n3PHjm2Dea9AVqp5ucTtaEToIv45WXrHjh0aERIR053fuDUzx8Zny/bx7h/bcTigXf0q1Kvmx7S18Tx0bT1eva252ZHdh8MBU3rBvqVQqTbcMgaa9XL2N5MyqcQfjcXFxXH11VeXuVVWejQmIu5u2poDvPrLFtKzbVgMsDugqp8XUx5uWz6X11/M9t/h9xch0bkYh0Y3Qfe3oGpDc3OJS5hSCLVq1arMbV6oQkhE3F29Eb8VeE65Xl7/T1lpsOx9WP4B2LLA6g0dn4OOz4Jn8e6FJ+Yq7M/vQu8jdPvtt1/y/aSkJAwNMYqIlLhx/SIZPj2OHPuF/661WgzeuyvChFRuyssXbvg3tLwHZg+HPQth8VhnA9eoR81OJyYo9GTpX375hYyMDAICAvL98vf3d2VOERG5iD6tQpg1uEO+73lbLWTl2NF00PNUawQPzIS7JkPDG881dAVnx3spNwo9ItSsWTPuuOMOHnnkkXzfj42N5ddffy22YCIiUnRnl9cbgANIy7bx4o8b+X3zYcbe0ZKgSj5mR3QfhgHN+zq/zsrJhM9ugrA+0H4IeHiZFk9KRqFHhFq3bk1MTMxF3/f29qZOnbKzj0V0dDRhYWFERUWZHUVEpEDnL68Prx1ANX9vBnduiJeHhYXbj3PT+4uZuSFBo0OXsvlHOBwH81+DjzrAnsVmJxIXK/Rk6czMTGw2G76+5WvlgSZLi0hpcf7y+iybHW8PKzuPnmH49DjiEpIAZzPXN/uGU72it8mJ3ZDDARunwR//htTjzmMt7oCbR0OlWuZmkyLRztLFRIWQiJQFOTY7Hy/Zw7g/d5Btc1DZ15PXe7egV0Sw2dHcU3oiLBwNaz8Fhx28KkKXkdDuSXW2LyWKvRBKTi5c47qyViyoEBKRsmTr4WSe/z6OLYedf6f3DK/F672bU9Vfo0P5OhQLvz0PB9dBgy7OCdZaIV0qFHshZLFYLrk83uFwYBiGNlQUEXFz2TY7ExfsInrhLnLsDqr6eTG6bwtuaaFHP/my22HDVKjb0bnaDCD9NNhynMvuxS0VeyG0ePG5CWMOh4MePXrw6aefEhISkue866+//jIjuycVQiJSVm0+mMTz38ex/egZAHpHBvPabc0J9NVKqQL98iz8NcPZy6z1QLBYzU4k53H5HKGKFSsSFxdHgwYNLjtkaaBCSETKsswcG+Pn7+TDRbuxO6B6RW/G3h7Ojc2CzI7mvnIy4bOb4XCs83VwK+j5HoS0NjWW5KXu81dIy+dFpDzw9rDyQrermPFUBxpW9+P4mUwembKO57+PIyk9G0Ad7c/n4Q2DFkD3d8A7AA5tgEk3OkeJ0k6ZnU6KSCNCBdCIkIiUFxnZNt6ft4NJS/fgcEDNSj6MvSOcRduPM3nFPnW0z0/KMfjjZdj4nfO1b1Xo9zXUbW9uLimZESH1FhMRKTt8PK38X49mTH+8PbUDfTiSnMFDX6zl2zXObu2/xB1i88EkNiUkkXA6zeS0bsK/Btz+MTz0G1S/yrkPUfWmZqeSIrjspqsZGRk88cQT+Pn55Tk+Y8aM4kkmIiKmaFOvCgmJGbmvM3OcvbdOpmZx64RlucfV0f4f6nWEJ5bB8e3gW8V5zOGAdZ9D+F3goycK7qrQI0LnN1m9//77CQ4OvuC4iIiUfuP6ReJhyX/U38NiMK5fZMkGKg2snlCzxbnXW3+B34bBxCjY9IOzMBK3U+gRoS+++MKVOURExI30aRVCoxr+eUaAzgoPCeDahlVNSFXKVKgMVRrCqd3w4yOwfrJzdZkenbkVrRoTEZFLOjsd9Oz40Ib4RLp/sJSF24+ZlqlUqH8dPLUSuvwbPHxg31L4sAPMewWyUs1OJ39TISQiIvnKr6N9ZV9PGlb342RqFgO/WMsbv24hM6dsdRQoVh7ecP0LMHg1NOkO9mxYPg6+vdfsZPI3NV29iOjoaKKjo7HZbOzYsUPL50WkXMqvo73DAWN/38bkFfsAaB5cifH3tqJhdX9zw5YG23+H31+EW8dBoxvNTlOmqft8MdE+QiIi+ftzy1Fe+CGO02nZVPC08lrv5tzVura2VilIThZ4/KONydpPIeU4dHwOPH3My1XGaGdpERFxqa5hQfz+TCfaN6hKeraNF3/YyNPfxZKckW12NPf2zyIo9STMexUWj4X/tYMdf5gWq7xSISQiIpetZoAPXz3ajhe6NcVqMfgl7hA9PlhKzIHTZkcrHXyrQO8JULEWnN4H39wF3/WHxANmJys39GisAHo0JiJSODEHTvP0txtIOJ2O1WIw7KYmPHF9Q6wX2Y9I/iHzDCwaC6s+BIcNPCo4J1m3H5p3BEkKTY/GRESkRF1dpzKzn7mOXhHB2OwO3pm7nQc+W83R5IyCP1zeeVeEbqPhiaVQ51rISYcFo517EIlLqRASEZFiU8nHk/H3RPLOnS2p4Gllxe6T3DJuCX9uOZp7jrrZX0JQcxg4G/p8BF1GQo1m597T3kMuoUJIRESKlWEY3NUmlF+f7kjz4EqcTsvm0anrePXnv8jItjEj5iAr95xkRsxBs6O6J8OAyHuh0wvnjh3ZDO+Hwcr/gS3HvGxlkOYIFUBzhERELl9mjo2352zns2V7AahX1ZfEtGwS07Op6ufFlIfb4nBAZT9Palf2NTmtG/ttOKyd5Px1UAtnq44615ibyc1pH6ErpA0VRUSKT70Rv11wzAD++QNI3ewvwW6HDVPhz1ch/e8VeZH94abXwa+aqdHclQqhYqIRIRGRKzdrw0Genx6HzX7hjxwPi8G7d0XQp1WICclKmdSTMP9ViJnqfO0TCN3ehFb9zUzllrRqTERE3EafViH8NLhDvu99OqCNiqDC8qsKt02AR+ZBzXDISHQuvZfLpkJIRERK1PkdOJ79LpY1e0+ZE6a0Cm0LgxY5V5dFPXru+JHNkKb/lkWhQkhERErE+d3sm9asiNVikJiezX2TVjF15T40W6MIrB7O1WVWD+fr7AyYdj9MbAMbvnLOK5ICaY5QATRHSESk+JzfzT4pPZuXf/qLX+IOAXBn69qM6tMCH0+ryUlLoVN74dt74Pg25+vQa5yry2q2MDeXSTRHSERE3I63hzW3O71hGAT6ejH+nkhe6tEMiwE/rE/g7o9XcjAx3eSkpVCV+vDEMudKMk8/iF8FH3eCOSMhI9nsdG5LI0IF0IiQiEjJWLbzBEO/jeF0WjZV/LyYeF8rrm2opeGXJekgzB0JW35yvvavCY8vhoo1zc1VgjQiJCIipUrHxtX4eYhzN+pTqVk88NkaPlu2V/OGLkdACNw9Fe7/Eao0gFotwT/I7FRuSSNCBdCIkIhIycrItjFyxiZmbnC24OgdGczY21tSwUvzhi5LdoZzib1/defr1BOw+mPo+Cx4+ZkazZU0IiQiIqWSj6eV9++O4JVeYVgtBj/FHuKOD1cQfyrN7Gilk6fPuSII4M9XYMnbEN0Otv4C5Xw8RIWQiIi4HcMwGNihPl8/2o6qfl5sOZxMr4nLWLrzuNnRSr+mPSGgDiTFO5fbf30XnNpjdirTqBASERG3dU2DqvwytCMRtQNITMtmwOdr+Gjxbs0buhJX9YDBq+G64WD1gl3zIPoaWDTW+RitnFEhdBHR0dGEhYURFRVldhQRkXItOLAC0x5vz12ta2N3wNjftzHk2w2kZuaYHa308vKFG1+GJ1dCgy5gy4RFY2DJO2YnK3GaLF0ATZYWEXEPDoeDr1Yf4PVf/iLb5qBpUEU+fqA19ar5sTEhkTGztzGyx1W0rB1odtTSxeGALbNg6fsw4GeoUPnc8fP7oZQimiwtIiJlimEYPHBNXb4ddA3VK3qz/egZbpu4jIXbjzEj5iAr95xkRsxBs2OWPoYBzfvC40vyFkHf9XcWRzlZ5uZzMY0IFUAjQiIi7udocgYPT17LX4ecOyb7ellJy7JR1c+LKQ+3xeGAyn6e1K7sa3LSUmrXfPjqduevqzWBHu9Cg+vNzVREhf35rUKoACqERETcU70RvxV4zr6xPUsgSRnkcMDGafDHvyH175V6Le6EbqNLze7UejQmIiJl2rh+kXhY8p/D4mExGNcvsmQDlSWGARH3wJB1EDUIDAts/gEmtIGV/wNb2ZmorkJIRERKpT6tQpg1uEO+7z3coR69I4NLOFEZVCEQer4LgxZCSGvIOgPrvwCH3exkxcbD7AAiIiJXyjDybpD8ydK9HDiVzlt3tiSggqd5wcqK4Eh45E/YMBWqNgIPL+dxWzZkJIFf6W2OqxEhEREptar6e1Hd35vwkABG921By5AA/L2teFhgzl9HuHXCUuLiE82OWTZYLND6IajX8dyxVf+DCa1h7Wdgt5kW7UposnQBNFlaRMS9ZebY8LJaMAwDh8NBls3OtsNnGPJtDPGn0vG0Gozs3oyBHephlOJ9cdyOwwFfdIcDK52vg1tBz/ch5Gpzc/1Nk6VFRKRc8Paw5hY4hmHg7WElIjSQX4deR/cWNcm2OXj91y08/uV6ktKyTU5bhhgGDPgVbnkLvCvBoQ0w6Qb4dRiknzY7XaFpRKgAGhESESm9HA4HX67az6hft5JlsxMSWIHo/lcTGRpodrSy5cwR+ONl2PS987VvVbj9E2jU1bRIGhESEZFyzzAMHmxfjx+fvJY6VXw5mJjOnR+u4NOle9S4tThVrAl3THKOEFW/CjKSISDU7FSFohGhAmhESESkbEjOyGbkj5v4bdNhALo2C+Ldu1oS6OtlcrIyxpYN8Wug3j+2NtjyMzToDD4l93NUI0IiIiL/UMnHk4n3teKN3s3xslr4c+tReo5fRsyB0jOfpVSweuYtgg7FwvQBMDEKNv2Qd58DN6BCSEREyg3DMHigfT1mPHUtdas6H5Xd/dFKJi1xPirbmJDIvZ+sYmNCotlRy46cDKhcD1KOwI+PwNTb4PgOs1Pl0qOxAujRmIhI2XQmI5sRMzbx20bno7Ibr6pB9YrefLc2noeurcertzU3OWEZkp0BK8bD0vechZHFE64dAp1eAC8/l9xSTVeLiQohEZGyy+FwMHHhLj6Yv5McmwOLAXYH6mLvKqf2wu//gp1zna+DwuHxJc7NGotZYX9+q8XGRURHRxMdHY3NVjp3yhQRkYIZhsF7f5x7TGP/e2jgZGoWt05YlntcXeyLSZX6cN802D4bfh8BbR5ySRFUFBoRKoBGhEREyrZZGw4yfHocOfYLfxxaLQbv3RVBn1YhJiQr47LSwMMbLFaXXF6rxkRERArhUl3srRZISs/Gnk+RJFfIy9dlRVBRqBASERH529lWZGc7kmXlOHjl578Y8MUajiRlmJZLXEeFkIiIlHvnd7EPrx1AdX9vnuvaGG8PC0t3nqDbuCX8EnfI7KhSzDRHqACaIyQiUj7k18Xe28PKrmMpDPs+lo0JSQDcFhHMG71bEODraXJiuRTNERIRESmC/LrYAzSq4c+PT17L0zc2xmox+DnuEN3GLWHpzuNmxpViokJIRESkAJ5WC8NuasKPT15L/Wp+HEnO4IHP1vDqz3+RnqVtVkozFUIiIiKFFBkayG9Pd+SBa+oCMHnFPnpOWKqWHKWYCiEREZEi8PXy4I0+LZjycFtqVPRmz/FUbv/fCj74cyc5NrvZ8aSIVAiJiIhchuubVOeP5zrRs2UtcuwO/vvnDu74aCV7jqcAqIFrKaFCSERE5DIF+nox8d5WfHBPJJV8PIiLT6TH+KV8uWo/P65PYOWek8yIOWh2TLkELZ8vgJbPi4hIYRxOSmfINzGs358IgKfVINvmUANXk6jpqoiISAmqFVAhtwgCyLY5xxnUwNW96dGYiIhIMRnXLxIPi5Hve1aLwbh+kSUbSAqkQkhERKSYXKqBa0UfDwIqaDdqd6NCSERExAXOb+CamJbNwMlrGTljIymZOablkrxUCImIiBSj/Bq4VvP34p6oUAwDvl0TT/cPlrB6z0mzowpaNVYgrRoTEZGiulgD11V7TjJ8ehwJp9MxDHikQ32Gd2uKj6fV7MhljpquioiImORiDVyvaVCVOc924p6oUBwO+HTZXnpNWMamvzvbS8lTISQiIlKC/L09GHtHSz4b0IZq/t7sPJZC3/8tZ9yfO8hWi44Sp0JIRETEBDc2C3K26Ah3tugY9+dO7vhwBbuOnTE7WrmiQkhERMQkVfy8mHifs0VHQAVPNiYk0WP8Mj5duge7XVN4S4IKIRERERMZhkHvyBD+eK4T1zepTlaOnVG/beXeSauIP5VmdrwyT4WQiIiIGwiq5MPkgVGM7tsCXy8rq/ee4pZxS5i29gBnF3iro33xK/OFUGJiIm3atCEyMpIWLVowadIksyOJiIjkyzAM+rery+/PXEebupVJzbLxrx838ciUdRxLzmBGzEF1tC9mZX4fIZvNRmZmJr6+vqSmptKiRQvWrVtH1apVC/V57SMkIiJmsNkdfLp0D+/O3U623YG/jwc4ICUzRx3tC0Hd5/9mtVrx9XX+HyQzMxOHw0EZr/1ERKQMsFoMHr++IWN+3wZASsa5thzqaF98TH80tmTJEnr16kVwcDCGYTBr1qwLzomOjqZevXr4+PjQrl071qxZU6R7JCYmEhERQe3atXnhhReoVq1aMaUXERFxrXH9IrFepKO9hzraXzHTC6HU1FQiIiKIjo7O9/1p06YxbNgwXnnlFWJiYoiIiKBbt24cO3Ys95yz83/O/zp06BAAgYGBxMXFsXfvXr755huOHj160TyZmZkkJyfn+RIRETFLn1Yh/HSRjvYdG1Wlc9PqJZyobHGrOUKGYTBz5kz69OmTe6xdu3ZERUUxceJEAOx2O6GhoQwdOpQRI0YU+R5PPfUUN9xwA3feeWe+77/66qu89tprFxzXHCERETHL5oNJ3DphGYYB5//Url7RmzF9w+kaFmROODdVJnqNZWVlsX79erp27Zp7zGKx0LVrV1auXFmoaxw9epQzZ5y7dCYlJbFkyRKaNm160fNHjhxJUlJS7ld8fPyV/SZERESu0Pkd7VvWDiDQ15O6VSpw/Ewmj05dx/Pfx5GUnm121FLHrSdLnzhxApvNRlBQ3io3KCiIbdu2Feoa+/fv57HHHsudJD106FDCw8Mver63tzfe3t5XlFtERKQ41QqowLIRXXI72t/Xtg5ZNjsOB7w/bweTlu7hx5gElu06ztjbW9LlqhpmRy413LoQKg5t27YlNjbW7BgiIiJX5GwHe8jb0f7/ejSjW/Mghk/fyN4TqQycvJa729Tm37eGUcnH06y4pYZbPxqrVq0aVqv1gsnNR48epWbNmialEhERcS+t61Zh9tPX8UjH+hgGfL8ugW7/XcKSHcfNjub23LoQ8vLyonXr1syfPz/3mN1uZ/78+bRv396l946OjiYsLIyoqCiX3kdERKQ4VPCy8vKtYXz/eHvqVvXlcFIGD36+hpEzNnImQ3OHLsb0VWMpKSns2rULgFatWvH+++/TpUsXqlSpQp06dZg2bRoDBgzg448/pm3btowbN47vv/+ebdu2XTB3yBW0s7SIiJQ2aVk5vD1nO5NX7AMgJLACb93Rko6Ny88+eoX9+W16IbRo0SK6dOlywfEBAwYwefJkACZOnMg777zDkSNHiIyMZPz48bRr165E8qkQEhGR0mrVnpO8+MNGDvzdxb5/uzqM7NEMf2/nFOGNCYmMmb2NkT2uomXtQBOTFr9SUwi5OxVCIiJSmqVm5vDWnG1MXbkfgNqVK/D2nS25tmE1Xv35Lyav2MdD19bj1duam5y0eKkQKiYqhEREpCxYsesEL/ywkYOJ6QDcGl6TFXtOcSo1q0w2cVUhdIWio6OJjo7GZrOxY8cOFUIiIlLqpWTm0OKVuRccN4B/FgNloYlrmdhZ2kyDBw9my5YtrF271uwoIiIixcLf24Nx/SI5v4fr2SKoPDZxVSEkIiJSjvRpFcLPQzrm+96Y21vQp1VICScylwohERGRcso4b2ToxR82MWb2VjKybeYEMoEKIRERkXLm/CauzYMr4e1hwQF8vGQPPccvJTY+0eyYJUKTpQugVWMiIlIWZebYcpu4OhwOsmx2lu44wciZmzh+JhOLAU92bsjTNzbO0+estNBk6SukFhsiIlKWeXtYMf5+Nna2iWvXsCDmPdeJ3pHB2B0QvXA3t01YzuaDSSandR2NCBVAI0IiIlIezdl8mJdmbuZkahYeFoPBXRoxuEsjvDxKxxiKRoRERETkst3SohZ/PNeJHuE1ybE7+GD+TvpEL2fr4WSzoxUrFUIiIiKSr6r+3vyvf2sm3teKyr6ebDmczG0TlzFxwU5ybHaz4xULFUIiIiJySbe2DOaP567n5rAgsm0O3v1jB7d/uIKdR8+YHe2KqRASERGRAlWv6M3HD7RmXL9IKvl4sDEhiZ7jl/HR4t3Y7M7pxhsTErn3k1VsTEg0N2wRqBASERGRQjEMgz6tQpg37Hq6NK1Ols3O2N+3cedHK9h9PIUZMQdZueckM2IOmh210LRq7CLUdFVEROTiHA4H09cn8NrPf5GaZcPTauBptZCWZXOLbvbqPl9MtHxeRETk4uqN+K3Ac8zoZq/l8yIiIuJy4/pF4nF+O/u/lYZu9iqERERE5LL1aRXCrMEd8n0vPCSAtvWrlHCiolEhJCIiIsXi/G72G+IT6fbfJUxfF4+7zsRRISQiIiJX5Pxu9i1rB1DZ15PmwZU4k5nDCz9sZNDUdRw7k2F21AtosnQBNFlaRESkYPl1s/ewWPhkyR7+O28HWTY7gb6evNG7Bb0iggHnvkNjZm9jZI+raFk7sFjzaLK0iIiIlJj8utlbLQZPdm7Iz0M70Dy4Eolp2Qz9dgODv47hVGqWW+w7pBGhi9A+QiIiIsUn22Zn4oJdTFiwE7sDAit4YHPAmYwcl+w7pH2EiokejYmIiBSfktp3SI/GRERExO2M6xeJ1Y32HVIhJCIiIiWmT6sQfrrIvkOzBnegT6uQEs2jQkhERERMcXbfofP3HypJHubdWkRERMqjs/sO1Qr0oV9UKNPWxnM4MYOq/l4lnkWTpQugydIiIiLFL799h7w9rMV2/cL+/NaIkIiIiJS4fxY9Z/cdMoPmCImIiEi5pULoIqKjowkLCyMqKsrsKCIiIuIimiNUAM0REhERKX20oaKIiIhIAVQIiYiISLmlQkhERETKLRVCIiIiUm6pEBIREZFyS4WQiIiIlFsqhERERKTcUouNApzdZik5OdnkJCIiIlJYZ39uF7RdogqhApw5cwaA0NBQk5OIiIhIUZ05c4aAgICLvq+dpQtgt9s5dOgQFStWxDAMs+MUKCoqirVr15odo9DMzOvqexf39Yvjepd7jcv5XFE+k5ycTGhoKPHx8drBvRD0fe4+99b3eeE+Y8b3uMPh4MyZMwQHB2OxXHwmkEaECmCxWKhdu7bZMQrNarWWqh8kZuZ19b2L+/rFcb3LvcblfO5yPlOpUqVS9f9fs+j73H3ure/zon2mpL/HLzUSdJYmS5cxgwcPNjtCkZiZ19X3Lu7rF8f1Lvcal/O50vb/xdKktP231fd5yV5P3+dFo0djImI6NTcWKdvc+XtcI0IiYjpvb29eeeUVvL29zY4iIi7gzt/jGhESERGRcksjQiIiIlJuqRASERGRckuFkIiIiJRbKoRERESk3FIhJCIiIuWWCiERcXtpaWnUrVuX4cOHmx1FRFygXr16tGzZksjISLp06VKi91aLDRFxe6NHj+aaa64xO4aIuNCKFSvw9/cv8ftqREhE3NrOnTvZtm0b3bt3NzuKiJRBKoRExGWWLFlCr169CA4OxjAMZs2adcE50dHR1KtXDx8fH9q1a8eaNWvyvD98+HDGjBlTQolFpKiK4/vcMAyuv/56oqKi+Prrr0souZMKIRFxmdTUVCIiIoiOjs73/WnTpjFs2DBeeeUVYmJiiIiIoFu3bhw7dgyAn376iSZNmtCkSZOSjC0iRXCl3+cAy5YtY/369fz888+8+eabbNy4saTiq8WGiJQMwzCYOXMmffr0yT3Wrl07oqKimDhxIgB2u53Q0FCGDh3KiBEjGDlyJF999RVWq5WUlBSys7N5/vnn+c9//mPS70JELuVyvs/P98ILL9C8eXMeeuihEsmsESERMUVWVhbr16+na9euuccsFgtdu3Zl5cqVAIwZM4b4+Hj27dvHu+++y6BBg1QEiZQihfk+T01N5cyZMwCkpKSwYMECmjdvXmIZtWpMRExx4sQJbDYbQUFBeY4HBQWxbds2k1KJSHEqzPf50aNH6du3LwA2m41BgwYRFRVVYhlVCIlIqVBSw+QiUrIaNGhAXFycaffXozERMUW1atWwWq0cPXo0z/GjR49Ss2ZNk1KJSHEqDd/nKoRExBReXl60bt2a+fPn5x6z2+3Mnz+f9u3bm5hMRIpLafg+16MxEXGZlJQUdu3alft67969xMbGUqVKFerUqcOwYcMYMGAAbdq0oW3btowbN47U1FQGDhxoYmoRKYrS/n2u5fMi4jKLFi3Kt2/QgAEDmDx5MgATJ07knXfe4ciRI0RGRjJ+/HjatWtXwklF5HKV9u9zFUIiIiJSbmmOkIiIiJRbKoRERESk3FIhJCIiIuWWCiEREREpt1QIiYiISLmlQkhERETKLRVCIiIiUm6pEBIREZFyS4WQiJRLhmEwa9asIn3m1VdfJSgo6LI+KyLuSb3GREQKYevWrbz22mvMnDmTa665hsqVK5sdSUSKgQohEZFC2L17NwC9e/fGMAyT04hIcdGjMREp9Tp37szTTz/Niy++SJUqVahZsyavvvpq7vs7d+6kU6dO+Pj4EBYWxrx58y64xqZNm7jhhhuoUKECVatW5bHHHiMlJQVwPhLr1asXABaLRYWQSBmiQkhEyoQpU6bg5+fH6tWrefvtt3n99deZN28edrud22+/HS8vL1avXs1HH33Ev/71rzyfTU1NpVu3blSuXJm1a9cyffp0/vzzT4YMGQLA8OHD+eKLLwA4fPgwhw8fLvHfn4i4hh6NiUiZ0LJlS1555RUAGjduzMSJE5k/fz4Oh4Nt27Yxd+5cgoODAXjzzTfp3r177me/+eYbMjIymDp1Kn5+fgBMnDiRXr168dZbbxEUFERgYCAANWvWLNnfmIi4lEaERKRMaNmyZZ7XtWrV4tixY2zdupXQ0NDcIgigffv2ec7dunUrERERuUUQQIcOHbDb7Wzfvt21wUXEVCqERKRM8PT0zPPaMAzsdrtJaUSktFAhJCJlWrNmzYiPj88zr2fVqlUXnBMXF0dqamruseXLl2OxWGjatGmJZRWRkqdCSETKtK5du9KkSRMGDBhAXFwcS5cu5aWXXspzTv/+/fHx8WHAgAFs3ryZhQsXMnToUB544AGCgoJMSi4iJUGFkIiUaRaLhZkzZ5Kenk7btm159NFHGT16dJ5zfH19mTt3LqdOnSIqKoo777yTG2+8kYkTJ5qUWkRKiuFwOBxmhxARERExg0aEREREpNxSISQiIiLllgohERERKbdUCImIiEi5pUJIREREyi0VQiIiIlJuqRASERGRckuFkIiIiJRbKoRERESk3FIhJCIiIuWWCiEREREpt1QIiYiISLn1/2lo+VYbi1JdAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.xlabel(\"ndof\")\n",
    "plt.ylabel(\"H1 error-estimate\")\n",
    "ndof,err = zip(*l)\n",
    "plt.plot(ndof,err, \"-*\")\n",
    "\n",
    "plt.plot(ndof, [(ndof[0]/n)**(1.5)for n in ndof], \"--\", label=\"2 th order\")\n",
    "\n",
    "# plot the convergence rate line\n",
    "\n",
    "plt.legend()\n",
    "plt.ion()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5deda0bb82924648939f44d6bcabb034",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1e6825d2a9e74e2ab4d4d62a56298f5d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {'Objects': {'Surface': False…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 185,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 10\n",
    "# fac = 0 if mesh.dim == 2 else 1\n",
    "# p = [(-1+4*i/N,-2+4*j/N,fac * 2*k/N) for i in range(1,2*N) for j in range(1,2*N) for k in range(1,N)]\n",
    "p = [(i / N, -1 + (2*j) / N, 0) for i in range(N+1) for j in range(N+1)]\n",
    "\n",
    "fieldlines = flux._BuildFieldLines(\n",
    "    mesh, p, num_fieldlines=N**3//5, randomized=False, length=0.3)\n",
    "clipping = {\"clipping\": {\"y\": 0, \"z\": -1}}\n",
    "\n",
    "\n",
    "Draw(-eps*grad(gfu), mesh, \"X\" , min = 0, max = 10)\n",
    "Draw(-eps*grad(gfu), mesh,  \"X\", draw_vol=True, draw_surf=True, objects=[fieldlines],\n",
    "     autoscale=True, min=0, max=10, settings={\"Objects\": {\"Surface\": False}}, **clipping)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}