{ "cells": [ { "cell_type": "markdown", "id": "organic-pierre", "metadata": {}, "source": [ "# Coupling to the interface\n", "\n", "This is an example of how to use Dune-MMesh and solve coupled problems on the bulk and interface grid." ] }, { "cell_type": "markdown", "id": "ambient-charge", "metadata": {}, "source": [ "## Grid creation\n", "\n", "We use the [horizontal](grids/horizontal.rst) grid file that contains an interface $\\Gamma = [0.25, 0.75] \\times {0.5}$ embedded in a domain $\\Omega = [0,1]^2$. Grid creation from a mesh file works as follow." ] }, { "cell_type": "code", "execution_count": 1, "id": "round-motel", "metadata": {}, "outputs": [], "source": [ "from dune.grid import reader\n", "from dune.mmesh import mmesh\n", "dim = 2\n", "file = \"grids/horizontal.msh\"\n", "gridView = mmesh((reader.gmsh, file), dim)\n", "igridView = gridView.hierarchicalGrid.interfaceGrid" ] }, { "cell_type": "markdown", "id": "equivalent-spanking", "metadata": {}, "source": [ "## Solve a problem on the bulk grid" ] }, { "cell_type": "markdown", "id": "general-persian", "metadata": {}, "source": [ "Let us solve the Poisson equation\n", "\n", "\\begin{align}\n", "-\\Delta u = f & \\qquad \\text{in}\\ \\Omega\n", "\\end{align}\n", "\n", "on the bulk grid. We use the manufactured solution $\\hat u = \\sin(4 \\pi x y)$ and therefore apply the source term $f = -\\operatorname{div}( \\nabla \\hat u )$. The weak form of the problem above reads\n", "\\begin{align}\n", "\\int_\\Omega \\nabla u \\cdot \\nabla v~dx &= \\int_\\Omega f v~dx\n", "\\end{align}\n", "for all corresponding test functions $v$. This can be implemented as follows." ] }, { "cell_type": "code", "execution_count": 2, "id": "bearing-lighting", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5.628259763933402e-07" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from ufl import *\n", "from dune.ufl import DirichletBC\n", "from dune.fem.space import lagrange\n", "from dune.fem.scheme import galerkin\n", "from dune.fem.function import integrate\n", "\n", "space = lagrange(gridView, order=3)\n", "u = TrialFunction(space)\n", "v = TestFunction(space)\n", "\n", "x = SpatialCoordinate(space)\n", "exact = sin(x[0]*x[1]*4*pi)\n", "f = -div(grad(exact))\n", "\n", "a = inner(grad(u), grad(v)) * dx\n", "b = f * v * dx\n", "\n", "scheme = galerkin([a == b, DirichletBC(space, exact)], solver=(\"suitesparse\", \"umfpack\"))\n", "uh = space.interpolate(0, name=\"solution\")\n", "scheme.solve(target=uh)\n", "\n", "def L2(u1, u2):\n", " return sqrt(integrate(u1.grid, (u1-u2)**2, order=5))\n", "\n", "L2(uh, exact)" ] }, { "cell_type": "markdown", "id": "combined-hurricane", "metadata": {}, "source": [ "## Solve a problem on the interface\n", "\n", "We can solve similar problem on the interface $\\Gamma$ like\n", "\\begin{align}\n", "-\\Delta u_\\Gamma = f & \\qquad \\text{in}\\ \\Gamma\n", "\\end{align}\n", "with the weak form\n", "\\begin{align}\n", "\\int_\\Gamma \\nabla u_\\Gamma \\cdot \\nabla v_\\Gamma~dx &= \\int_\\Gamma f v_\\Gamma~dx\n", "\\end{align}\n", "for all corresponding test functions $v_\\Gamma$." ] }, { "cell_type": "code", "execution_count": 3, "id": "through-munich", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5.807742030532398e-08" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ispace = lagrange(igridView, order=3)\n", "iuh = ispace.interpolate(0, name=\"isolution\")\n", "\n", "iu = TrialFunction(ispace)\n", "iv = TestFunction(ispace)\n", "\n", "ix = SpatialCoordinate(ispace)\n", "iexact = sin(0.5*ix[dim-2]*4*pi)\n", "iF = -div(grad(iexact))\n", " \n", "ia = inner(grad(iu), grad(iv)) * dx\n", "ib = iF * iv * dx\n", "\n", "ischeme = galerkin([ia == ib, DirichletBC(ispace, iexact)])\n", "ischeme.solve(target=iuh)\n", "L2(iuh, iexact)" ] }, { "cell_type": "markdown", "id": "joint-workstation", "metadata": {}, "source": [ "We can use the `plotPointData` function to visualize the solution of both grids." ] }, { "cell_type": "code", "execution_count": 4, "id": "experimental-hopkins", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAANAAAACYCAYAAACVr2d1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABTTElEQVR4nO29fbgt21XW+Rtzzlpr33PvDTeQAMEgiIL4ASIi4mO3fLbQ2BIUxcDjRxREBBqFhhYebKGl7Y6ioggdtBEQmhYwmhAViFFAUCESWyURWggQYhK+QpJ77zl7r1VVc47+Y4xZNWvttfbHOfvcfe7NHudZZ61VVat21Vr11jvGO8ccQ1SVG7uxG7s7C9d9ADd2Y09nuwHQjd3YPdgNgG7sxu7BbgB0Yzd2D3YDoBu7sXuwGwDd2I3dg90A6MYeSBORbxSRXxKR1x1YLyLyNSLyehH5MRH50GbdHxeRn/LHH2+W/zYRea1/5mtERO71OM8F0L2cyI3d2D3YNwOfcMb6/x54f398FvASABF5V+DLgd8BfDjw5SLybP/MS4A/1XzurP1fyC7CQN98zh/aeyI3dmP3Yqr6g8DbztjkBcC3qNmPAI+JyPOAjwdepapvU9W3A68CPsHXPUtVf0Qte+BbgE++1+M8F0D3cCI3dmP3034V8F+b92/yZWctf9Oe5fdk6V53wOED/vndDUXkszCWQtar37b6Vc8hhEIXCl3I9pD5kShEMZQHBGF2WRVQlIKSVRkJjBoYNDJosucSGUogl0DJAYogGXuU9lntURRKgfqsiqqCpzuJCIhAEJAAIUAUNPgjChpBA5QIBNAIxEIISoplOsdVyCTJdNj7iBJFmP/N54mfZwGyMp1r356nRsZs56pZ9p5r8HOlKJKX54kWtGj9oRDBz1HsPEOYzzH4efq5agQV6N/0preq6nPb3/zjP/qWvvVt5dRF8//+2PY/A5tm0d9V1b97sUvuwbGrANCFzb+gvwuwft/n6/P/j8/mWY9seO7Dt3nPh57keUeP8/zV23iP7nHeMz7Ou8UTHgvKI5JYS6KT+XALhUEzxzrwZMk8XhK/kB/ll8Zn8fPDs3nz9jF+YfMu/Pzxo7z9zi3u3F5TnlwRnwx0Twrdk7C6DavbhdWThfTkQLrdE+5s4c4JbLfoyQbte8owghYkdUiXkPUaeegIHjpCHz4iP7xmfLRjeCSyfVagf1QYHoHxERieVdBHR9aPbnn2wyc899ZtnvfQE7zX0Tt4/uptPDc9yXvGx3ksbHmXALfOOddjzbwtJ36l3OIXxsd4y/AYv9A/xptP3oVfOnmEX779CHfuHJGf7IhPxulcu9uwuq2snsx0T47E2wPhzhY5PoHjE3SzRfsBHQcAJEZktULWazhawyO3KLdW5EfWDI8k+mdF+kcCw6MwPAJ5DT/1F77w53Z/819+W+bffO97nboWbr3XGzaq+mH3cDm9GXjv5v3zfdmbgY/aWf4Dvvz5e7a/J7sKAB06kQtbUaGo3XOze5WFQEbIFIrfhwuFsON1RsRYStTu7M5e6zCyCiNdyMRYkKhoKmgMaAJNUBKUKJQkaBfQZA9JCYYBYoQYkZzREkD9rj2OkDOMIzJkwpAJfSQMgdgrsRfKFsoKwlbIq8CwShx3HU+mNY90ax4fH+KReIu1DNySLUcystaRTgqRTCQszrWTyJF/F4+GkZ4tQ3ySTekYusg2J/qS2ObEmAMnWchZkMq8BaQIUgJSEihIKYgz7MTtWtCc7bHZEsDWiRzw921pGPb/tooyaL7M5XBRewXweSLy7Zhg8Liq/ryIvBL43xvh4PcAX6qqbxORJ0TkI4BXA38M+Nv3ehBXAaC9J3KRD2oDHKhAChQ18LRWUGLz3i4u+2Gib9uRCVLoZDQXMGRWMU8gylHRpJOrVSYQGZA0BgNNig6eYO5MjGgZzcXJ2ZaPIzImGEZkSIQ+E/pA3ApxpeSVEHqIW6GsAmUd2fQdd9KKJ7ojHk49j8db3Ao9D5ctRzLQSWGlhY7AQGYtYTrXQiER6SjcEmUjA0PY8G7pNoNGtquObUmMJTDmSC6BbRbG3CFZGKsrVxxUOSJlBaoEVcjFvkVV6Hs0ZwNT39t3ABAcRCLo4ucJZwAIBk67cOeZiPwDjEmeIyJvwpS1zg5Rvx74buATgdcDx8Cf8HVvE5GvBH7Ud/WXVLXG8J+DiWIPAd/jj3uycwF0tydyGSsqFITid7OJhVTIClmUQoEFhCwusgesKEQpEwNZrFEs1vAYhKBoVANMwkEklM4YqHSBkAKSIoSAOANpLkgQA5AWyAWGEdJobDRkZCjEvlAOsdAmMnaRTdfx5LDmVv8QD4WBW3HLw2HLURh4WHvu6OgxkTDouHDlogTWJDLKoyGT6Rn0mE3s2GjHxkHUe+ynReiLMJaElACFho0iojUmguDxjxRFi4GpDCNaRoIEY6kgUJTgF440M2HysH9IRYGtXh5Aqvpp56xX4HMPrPtG4Bv3LH8N8JsvfTBn2LkAupcTuRurzJP1cmO8QSCiE/vMblymi8ZAIRVKchZKHggnB1HE3buARlm6cWG0gJqMFkXczTEWysjEQpGwLcRuZqG4hdIJZS2UbWLTFe50K55Ia26lIx4ejYWOZOBh6enkmDsl04VCkDC5rZWFqisHMEhmCBv6GOk1MnSRUaMJDM5COQfGLIxZkNy6chioNNGpIqpIsf2KFnRjzlkZxikugs0pJpp+t6P9v5eq0j+D55w9pSLCPqvfrbGNuXDZGSkjdtP0GCirXVTVogSLgRCiFr9zmwu3DgPrOJKqupcym2gMpJGJhUrcYaLqvqUIKSFjtvc5U115zdnYKRdnoAHpE9JFYh9mFtoKpcNcuY0wdoGxi5x0HbfTmifSwEOxNzcubDkqA0cysAqFjeaZYRtXDnBXTnk0KIOO5HjCgKuPHg+NR6bOlSLcLkIpwlhMLpMCouIAMtcOheixkNQfpQ8TiKggkq2ByJkIB10Ylt7B9PsiDOxnp2eCXTuAgEUcBJxy5c6z6sYZA81CQnK5eBEHpSYOasWEziTakoTQRXPjksdBMYIEc+NcUJCcXUiwxyEWKtWV64SwFnIf6fu0YKG3R3Pl1jJwJL0xqA5nCArJZG11V64MDOGYISUKwnbVMWigL8ZCYw6cFCEXgRIcPMZEaDBXzkEzgQiX0Esxeb+ybm/HINN378/9/ktJgUFvAHR/bCEghAWQ8gQiIauSUdYhnFLiAoFYx4tE97pxSQopFmIsFgc5iEoUZyIXFbodNc5jIEQMSDlMjr9OABqRMRqI/HXLQsVZqPRQNoGSlDElNqlwZ2KhgbeHhzmSgaMwcCSmxnVOeYFAJ7CryhUKt4BeRnKweKiPkb5Lxj4aGIt9r6UENkXImhA1Zc4G08TeqzGI6HpS3qRxvbQfDERjsTgQi5fs+ECH/UqbAeiZm3J57QykKugCODLFP6bKMd3uspZFUL2rxFXwRFFWu0pcdDEhFUpUExEq+0wPmdQ4jXEC0ASkMIIaC4G5dTMLjcg2IikR0sxCpVNyL4QtxA5KF8irSL81WfvxdMRRHF1Q6B1A7sphg8xbRgICDYgCgUQEgYdDoJRMCSf0GikxsO06MsGUOfUBVmCrwlgAFxEqiCoTUSCpuWfSytsA42juHBl6Hy+yHxFZr/b+vgWhZ797d5aJyCcAf8u+aL5BVV+8s/6rgY/2t7eAd1fVx0Tko4Gvbjb9QOCFqvpyEflY4KswzN8GXqSqr7/0wTV27QCq1krYYCJCFRSqhlNjodZaJS6yZKCjMJwaDwqxERL2uXGuxkkX0WSAIEUY3ZUrCqI+ap8nEMkwWMzUjxYLbYPJ100spNv2b0Q22447KfNEWnMUH2IdR26F3l25YXLlohSCjtyS1YKBO7HxnFsCWQZyyDzGiV20Gi2G1PmRVUyZ0xXjxDyVbVqWWJFgZiJfuhAW+t6OooJsHPf+roqwKZe7zEQkAl8H/HdYZsuPisgrVPXHp/2qfkGz/f8I/FZf/v3Ah/jyd8XU4X/um74EeIGq/oSIfA7wF4AXXergduyBARDM7LOUtM0bP3MwVUxI6MSEhJoeUx+LOCgVxmgDqotxoOrGJTGFLgVIwVy3lCAOLhzk0yw0jgaucUSGaLHQEInbQumElNSUuFRZSMhdZIiJk9TxRHQWigMPhYfpZOQoDKw8xScwEoKw1YG1dIvvIEqgqHJLEoURQm9gSYGC3ZAGV+aKs/2TKgwKo6aZeRQgNNL0DohccWtBVPqeoMVk79V+BrJxoEsz0IcDr1fVn8H+9rdjOZc/fmD7T8OGV3btDwLfo6rHzeE8y1+/C/CWyx7Yrl07gHQaA2rcuB0WyhTyHvbZtYmFpsHUkXUYJzl7lTIx1KyEWc4+5calQEiCVjEhWnaC5jIpckgwFhI9yEKhCxYPJZ2AU7YQk6AOou224zhlHncQrYOBpwKok5FObYA1IgQyncTmnC0+AriFkhl5NPT2vcWaGxgnhq8x0W2FUYVR6zdXf5D5tWg3fa8tE1V1rmYsCNh3c+D3HfTSANqXX/k79m0oIu8D/Brg+/asfiHwN5r3nwl8t4icAE8AH3HZA9u1awXQ7vBAmeRrk7IvYvUOHDylp5NCx6zETXGQy9k1DhqTUrqlG5c7ITpbVCFBO2MUkjNNk5mAFosf2lhoGJAUCX1Eu0DYBmKSeWA1YdK2Cwo5JU7SilXMvCMdsQoj62gyfCfZVblCpz2BWdrejYcqiDIKPsia9ZiSwnRDqkCyuBNuO4BGmNw5fI9mCRUheRqP1IyEOpTgGQua83JUtf2NEXrde5k9R0Re07y/22TSFwIvVV3mC/mMgA8CXtks/gLgE1X11SLyxRi4PvMu/uZk185AsJMLV9mH6oJYNgLC3myEatN4kDNQwN25ykBxpIs2HjTJ2S4mTONBEwOJZSW0YsLezISGhdyNk2lcKBK6SEyBksQEBY+zygYHbSAnZYiJO3HFKmVjoTiylgqgYR7fCgPRsxTw72SpSApHEi1fTjIlbCkE+hTJhAlAYMG9qnAHA9Dge6DZW2utOwctGzUg2mOmwu39zd56RjLpZfIrX8j+gfxPBV6mqgOAiDwX+C2q+mpf/x3A9x7Y54Xt2gG0zIVr1beliHBoMLWO0VchITgLrdyFO3JBYVWzEoLJ2RIKGuNpN67bAVEKs5iQM8QRYkZ0DwtNsVC2HLk+ElIgpoB2tr+ytb+lW6a/UVJkmzpuJ2PJVbDjXYfRATS6OFKgZAgjj4glty4GWQUigawKAYuJOFnkFU4KZ7Psjsp+EEmYRxqCg0hkmtYhwad39IMNtO7x4hRh2M9AZ9mPAu8vIr8GA84LgU/f3UhEPhB4NvDDe/bxacCXNu/fDryLiHyAqv4kJlD8xGUPbNeuHUAwJ5UWV42mMaAGRFm1ue3ttyiWkRCdfWoMMcVBMZMmOVsZU5kGULWNgTpP9XG1LHYRhghDQLoE44jWgdWyEwuNIzIESBHpI5ICoQuErZCSUJI2ggWTYFFduRQLXcysPB7qwtgIIiOrsLGsC0bWkhaCQk33Wbs6Z4tHMieAj7V17Q2rARFHSxBNq2zfKkzuXBRZMlEI0ItR2Y4VhI3HUhc1VR1F5PMw9ysC36iq/1lE/hLwGlV9hW/6QuDbdac+tYi8L8Zg/2pnn38K+EciUjBA/clLHdgeu+aB1MOr5oHUJYjSARBFMbcmMMdBrZxdE0vXcTQ3LmVyigs3bs7QFnInBL/4tSaYVjk7JRuh38dCOaPjiAzu0vXBXLkY0C6QNjr9jbjBM8MDOZordxxWdKEsWGgCjzNRkC3R4yFgMcgaCJ5Dp+CJuDk4E8XdZN2lm9aCSCvLi7+ukwnBcuFECDHASZixdswpu0sRAVX9bixRuV32F3fef8WBz76BPbNNVfVlwMsufTBn2LUz0G4uXHbptebD2To8BtovZe8KCVGrEjcLCeto40GrmEnBsrMlFhcL9qhxnbj0bABoxQQZs83QPMBCkqMlok4slAlxFhRSdeOmMajZletjx51YSNGl9+Y8VpIJnq5EGAmM3PLvYC1LMFSlztaPEGYmOsvuyJqRSibmwin4FAZ36cQvnGAunblzB5JJ4ZCI8IywB+LMducETYzjIoKBx9fvzAtqrQoJnegUBx25mrWWwUGUWSVjoSGlSY2zi9gv6Brw+5hQ6QKSZveNFCFbTLRPkdNhMOl7GCHERSxUkliuqmd9a6jsF9AAOSqb2BF9+neSYgmxDqIgNm0jckwMI1GyM68QJbBPmTtCKWTLanXbzTMMjYp2B8gCo2CKW9WwRdBg7zUIKTp4RBA5BKC7Y6Cni10IQBdIq/jVwN8HHvNtvsQp+EI2xUCNS1EIzTiQjQXZ8n3zgiylZ5mRUCbXZx0GD87HSc6OsSzcuMNjQuaqaBfRIZuYEMdZkSt6ioWkzhfyugmSoo2dpECKDsqtTgO45sqZO5dTZAjKSSx0Yc2qiYei3xRWPlgcRAklE4KAjqxJB+VtsyWIAOKeyF9EuSPKKB2DJDTMzKMiJi4E8QwgIYbAAQK6AdBF0iqwlIjvVNWXiMhvxHzX973IAbR5cMA0er5vPlBGORSOVk88ihBUF+NBqzq9weXso2Su3KZx4+aMhIaFOrV0nEEJQzRJO0XoXJE7wEKascTTnGdZOxqLhRiI2+Kxj0yxlxXoqPGQZSncCSYqJLFH9GnrASUkYyJCT9TMLWHK3N4VFloQBTJF+r2jAUF8agKKCNwRGARGiVOBEQ1+nAE0RIjmIRD3B6dFhW25nIjwdLKLMNBF0iruKkWiBc+kwGmdByRLFlIosl/KhnluUNDdOGgeU6l5cQs3LiZy59nZXavGLV25KT8uR5vKXVkoJWOhGNFc5wphrp5XspmSUnsbS4rRklZt2lFzYcZGVAiJbYDbQYmhkEJ16byKj7tyAJEtMdSLP9DOIarLFkwUMpTestOZWShKMVbD/mYQ5bYovXSMwZloF0QxoCGhBwB0xjjQmXYBr+dFWGJoHR/6WlX9Bl+Xgdf68jeq6if5cgH+N+APYVnIL1HVr7n0wTV2EQBdJK3iK4B/7kl9DwMfd5mDaAf3WjMhYY6DbJkpcfty4qapDVNmtrFQFJ0GJB+Kw+TGTWqcK2Ta5KzZeJA/r4QyuiI37LJQ9HGhmX32CgohICEQvBRWdeXmCX5Nqag4M9FJ6AjBWSgUu7gpBB8sjqm4y7q1ulUM3JLu1HTwXRCVkAk6YhGPsU8Fj71Xgigiyu2gbAXy5M4ZeKjlvILdEPaZIgzlcgC6oNcD8B2q+nl7dnGiqh+yZ/mLMHn7A1W1iMi7X+rA9thViQifBnyzqv51EfmdwLeKyG9WXU6Gb+vCxXd9DJhZqDTP1YVrJWyriXb2WFCbmd2JErXM+WSTnG0xxcpVrjrVOydFPbVnmr/TCbmD0AysShfR0Sv3pGzToHO2XDBx4OwKClWlimF25YLYAGuAksISQMGC9RxgFOUkKjHYxRzQ6WK3eVB1mRK0BwqRkU4sNecQiMzlNXUO7kxsBkwKZgoWZ0VRnhRlE5QckyXYhmCDtc5GJV6piHDZZNKL2p8BPr1el6r6S/e4vwsB6CJpFZ+Bl/9V1R8WkSPgOcDiAHfrws3L7XkXPNWqkAB4gcHTSlww737KzK5xUKBM7LNupjccxZE7ntozpEhJhRIDpVPC4CAa2mkOPtmuC+iQkFTmcSFno30sJLk0eXLN2FDwDIUoaFB/VPdoBlEJJircEXNJU6iqnF3ksCMEhB7IJl97ZsLuBMQ2+TSQiV5Sp4ovK3ERuwFnDIUngnISlBySCwihKSB5gIH0rly4iyaTfoqI/G7gJ4EvUNX6mSPPsxuBF6vqy335rwX+sIj8fuCXgc9X1Z+67MG1dhEAXSSt4o3AxwLfLCK/ATjyAzzbmoHUsice2gXTRWu71PGglTPQJGeXvMjOngdVC2NU6AolRUqCUN23zsCUO0G64I8mFqos1HVWxbQsWUgzpgcHGzOqrhwhEKIQp1jCGIgdEOHx0BCU27L2wqjOODXob15HirOKgciyt+v3shQW5iHQDGEgNj9IHW+KoiSX02Owx3FcMYaOIUaffGgze/f/xF6L4bTdazLpPwH+gapuReRPYyrwx/i691HVN4vI+wHfJyKvVdWfBtbARlU/TET+AFa557+9xN88ZecC6IJpFf8T8H+JyBdgsHjRbnrFwf37827+G9SE0gomH1XHSlxllb1CQlGd4iDLi5vTYOqU6XUYOYoDq7jy9J7CkAo5BbQZEyoDDh4IuyxUAZQLuBsndcKdFpO1q6BQMxSGHVcuVCaqJXPdpYstiExCziGxFZOYK4DAp7GLpS+RcFfsxAda5wTPOgV8dwB6TfIbTgZGghxPibgmKhTSBCSLxUJQjoPSBxcXookL+39fYbt/Qt09JZOq6q80b78B+KvNujf788+IyA9gk+1+GmOyf+ybvQz4pgN//8J2oRjovLQKD+5+190exJxxUJW4ME1paGsj0Eyss8+dFhJgWWSkPuaKpVaDYB3NjasZ2kOXKYMnl3ZKGWSpxu2yUK7jQhHGOM8XamVtGkFhGBERdPCJeVIBhLGOtKrWrHJZCrQDS2Ar8GQFEnPQf9oMRKUB0b4KP5Ow4LlzoWRiPJlUvimrvSnSUpnoTlC20eIiTQdiIL2rmgjnej0i8rymgOcn4YmhXpH02JnpOdh1WcH1cmwa+M8CH4m5fvdk15eJsGcuUPt6d2p3lbKRs3Pi9o4H1cxsd+OOwjDJ2WsXFLaxEFImewpPSMs4qGWhsoqEUdFVguz11Erxgu1lIWvvunISgqlyIki0QDyEQJxYR6ZBypoBUAcsswSyRDbSNW7czEa7ltlOMVFGz1TngmdxEEaiFmBrN59UFnOr1mGc6ky8I2aeiEdsos1r2v8zH3ThDtoFvZ7PF5FPwuKctzFPzf4NwN/xhNGAxUBVfHgx8G3uKd3mHucCwbVPqBPaoiLmrs3q2zSte0rBn+3MlB6x8aCVCNtmULW6cfVCOIoDR7HjeCEm1DGhRkwYoYwNC42uyOWIDBlGTzTNcZa1YXbliiKhzKoczFMCRHbiIQeShIWogGcAZLE0mxOXmGUHPEH2RIrhfHXOCZ5bggsLhch2qjke6hSRMHqy61y08vGYOY7r/b8xMN5FVZ4LeD1fynK6Ql3+b7GJdPv2+Q7g9176YM6wByIXDpZTGqpZfYRmYp1Yiaty4I5bL4QaB1U5e6U1I6EWGhkWUxxWcaRLiZgKORZL7FyZGyfjDgutxMC0CkgOyMozs0tn9aWntij9aVeObDljwTO2Y4S+R2KYkjKXeWcslmkDolHg5JCk35QnmG4+0pPDyC0iRw6W3dy5qtqB5RWGkAnasyJPbnAVF9Zey7sO8L4jHpB49PIM9HSyBwJApwor1go9zJkJpyv07BcSqllgbG6c3UXHeTxoR0w4Sh2b0cSEsbPaZ5Zmo7Mr5yyURwhdsHltg80XCqPL2l3HVIgw58aVyxMbkTMEQXrmqdFh8IRMsfEhcVFhSuQUz4iu/5k8UAvu7rIQQOl83/UXjpDLACFTUG45iHZz5wDW0lGnj8NIDCYuzN+hM3nIjSBzeEbq3TDQ08WuH0B1DKhJ5ZlamzQStmVmz0rc4Zw4m1ZW3bgqZ3dSptmdEwvFkVUe52nU3cgw2piQdgFdKWVsWMiBlFcQRmtbIlnR0WXtXMydi3meMwQzePBUn+CA6Pt5VmdwAMX6HgdPmCazzbPYDFxFDER7puFMlgk2voVQwgnQe8bCyBqlI54q2rgcK4IomU4zke2pmhPT9xj2l7UqCH2+YaCnzKa5QA6eWpapunBViTsrpadaO71hpXmhxh15jtxDceA4ZNZpZDVaLJRTJHfFXLg2FsoeC43mzkkOlKymyOW4FBS8744yIBqphemBuZKNT4cWMXlNHDjB3bdYS0mFmYkAn2YtjAQsecdA1GZ1tIom7FR6pedRyQZSv4EdAlEVF2pc1LlLV1m8fpcP7WR5T6Y3DHT/zH88G3+cf3TgwA8/K3FwWEiYJ9gpnQh948ZZlZs1a68EelxW5oKEzhNME0OXyUOATtHR8uNkZAbOKIQVlCzkMSCjIjl5qxDPkdO5ZWQdYJVQptraNR5SXJGDWVQQm/EJgCRX5OpFuGSiFkTT17pbMrlN1K1VeuKGUjJZ1OcLFdZSy1jtj4s6AgETEur405H0DqRDhRVhLJcH0AWSST8bKyaSMUXts1T1x0Xkw/FsF/uy+ApVfZlnx/wgNpiasEo++2rJXcqunYGWP/bpvLg6FtQqcWcJCdWqnN1mJdT0lyPpWclqEhMeij1HccXGJe0+ZnIXyVkmFhJnHhl9cHXE2oVkoeTgTasSoRQkm6BAKqDJ2Kio1ZX2YiA1HhKwGguhMpCBaQZQlbTn70mb16iDSH0qti4FmVFnd3jQZO5xdGEmnJBDpsrcwFQueF9cNJB5VISoo2V5sFl8p3t/X4R8SQBdMJn0//H+VLic/TewdLLXAR/mUvjzgP8kIv8E2AIfo6q3RaQD/rWIfI83xr5ru3YAAXul7FofDpyNhMmN2xUS2HE/qkWxC6lmJazUJqIdhYG1txKxeKhzSTtxkkZWXWTMmTJ6YZHsAMoOoAw5i80KGIMvjxYPta6c993RXN266so5iJq4W1wokBAmJgrTeUxnZN9XFROme0iNF2FUOHYAtRMVJ5eugilGZ/UTiozcasSFonrQpcsq3BLYkK1BtG5ZYd/r/t/2rly4c5NJVfWJZvuH8W+jqUIKllJWlyvGVGAN4jqab/Bu7RoBdFqDnavyeInfVsZmN83nbCGhytlQ6ETonIFWHhAfhYHjkrkVeo7DiofiwCZ2rONIH5OxkMdCOoplJ4yCrAxAkisLOROtgnV9Kzq7ckUnEMEketGCaDHIutgGpDayWlx/s9MqftOx0ryuzilkhY0XT2xd4yme7OYbVNbAEDZWiFEsj27tf2NfDl2QMCWoDmQTGLybxKHf+bIMxAWTSUXkc4EvxIT7j2mW/w4sz+19gD+qqqMvj8C/B34d8HVNjbi7tmtmIFm4HK21QJql7FmJq0KCpfXvt3l6wzzJrhbnqGNCG7VY6CSvpoHV7Q4LlS5QRmMhyRYLSbYYqPj7kmsDX2OgMMVAFg+RM1LiTjzElC+nw+jfCE28U0XrXautSKZP+HsBrSASNlrZ3eKQUQOD35QGjfSaGGJiwNko9Kfiol2Xzr7X/QLDPlOFMe8F0D1XJlXVrwO+TkQ+HZsV/cd9+auB3+SJzX/fXbWNVy/9EBF5DHiZT7l53WX+5q49EC5ctdbdqApce6dsB1PPmhfUWnXjOlfjOsokJhxJ1+TH9ZyUjltpYJMTfUwMKZO7gGZjoTKqy9fOQJV9Mg4enIEiWtW40hlgus7YpcZDmSlre5GpIGLydj0BCcuYqNqUKAdT6ZwCYEDORcgFNqXJL/TvcCyRYTXXzc4qBiQ9NjYKIwMjR4SFSwf7BYZAOOzCMQtCO3ZVlUkBvh3rvLD829aF4TbWF/U1zfJ3iMj3M8dMd20PBICW/noLotaVmwdUd4WEfVO8WzcuoEQROtRZaC44ciQDW+lMTAgD25i4laxBVZ8jYxcszunEVDcHTekqCzmQii0nh8mFC6pQmFvJq7VFaaFQ8+UmZQ53zMfRxojcWiZa6Ce1RYmaCzi5cw6qokKvoN6hLpfQsJG1gRw6aw+Zo32/g/ZWoF4KWQZuedpPW9S+fse1b2s4eEfzzniXs4skk75/M5fn9wI/5ct/DfBfXUR4H6w/0Bu8tO/g4HkIEyj+ymUPbNeuHUCnZqQy10OY8+LCKRYqME1tsInchzOzp0HVZkyoiglHWt24LduSOCkdR3Fk48VHcg7OQgJZ0BzItVlvNrFN/G4v0+tgrFNMWECNhcDApDAPsnIARB43SdFJgZOpL+ncJVvU2c+BU/ud1uVZoZRE731SS7FmW32JjCUyroyFhi6y0Y5BE5tojbqGsGUImcLImsKRu3SnUoDgoEegCuWSMdAFk0k/T0Q+DhMf3467b8B/A3yJiAz2TfE5qvpWEflgzJ3zKYB8p6r+00sd2B67krJWvs2nAl+B3UD/k6qeqmV8ns1jQEtpO7sLUsSTTD0Osm3PyswOUzms6FMcOnQpJrgrt8tCRzGxTSN9jpZk2hVKHTgd1WIeZ56JlTLkXF05v/gLFg/l5KJB59kJPsgKM3h2QdSbNDydnhroQjvVSusj+LMxj1hBcaRYl25rMiwce8/UXGRiom1JBqKUPDaKU2zUs/U8OlPpavaCTXPYSUg9YHfBQBdJJv2zBz73rcC37ln+Y3gTrqu0KylrJSLvj2XG/i5VffuFizU0A6mtYgRMoJm6NDQpPbaeyYU7OzM7LNy4OLlxNpq+kW4SE448FtqWxCZ29GlgqG5cEfrRi8RnndW33DJRZSBvKa8RKc4mjRo3DbLCTqbCHhDV7O36RfnrUPenan9H8fb1YWpjj9px4ayYs4GoZPuehxInV7Uvke0qsS2JberY6m02oWOIdxhCoKfnSAoPy+hTwYWCuvt2GDx3Mw70dLKrKmv1pzBZ8O3A5Yo1HEok1VpsPiyUuDalp8ZB540HVTeu08CATmLCykE0SORIOm6FniFGv4gSfYn0aWDMftde2Z1bs1A85pnYRu1Oax2tmxbyWUGjsUZth9gAgVL2ytunYqLG5aMBT2RmG0rtuB0mMFeWFDVQ5wKahW2xAeCcA6MDqC+JfpXYaseg0W4iGk1Q0Q2Phg2DZHoZrS+rzkpdPJDUi8dfz1S7qrJWHwAgIv8Gc/O+QlUv3Huldqlb5nLNNeGqEjdoNBXOhYQaB51Xtm83tSejrLQw7GGhtaZpqsOt1NuFlRNDydao192hkgOyUmsVP7lLFUTOSC4pVzdrISqAg2JlPXamhr4HQLQDuuk1db86MR6VjdpnZ0pxJipZGPw8xhwYSmDIxkYnuWPbJTZdx0Y7ttpxHNZsYnIQDQxkHpbsLp3aZLwDVt7JAXTR/bw/8FGY5PiDIvJBPoFpsn1lrVqb6h8ws02Z5NfTQkLm/PGgqhQtUnsaSXuXhbbSccvduEEDt1Kiz9FdIGtAXHJBV0LxGIcGNPidv8Yh9jo4a6RpPoYBobp0aWcQ9QwmarY30HiclSsr4S6iz1fK8/FNjywmjmRhHG1+U86BYYxsc5puGiel46SsOI5r3jXdptfIJq7YhBN6tgwycCSZW1LoDrhxqkyZ6M9Eu6qyVm8CXu3dwH5WRH4SA9SPthstylq9z3srOudu7c0kpk7nbsDTCAltUinsl7Nbm8aEPMF0hbHQlKEdBgbdTu7LoObW3Eou93oAXla147Ula853e2YgqUxKGMyKWXt00jzrzvNBELkMrqWNgVwir+9LagBj7lzYidmkLF26vlR3LtBni4c2ObHNieOuY1sSx2nNcTlmEzs28YSNbHg49Aw+z+qQLasDXswukEz6u4G/CXww1sb+pc2678X6n/5rVf0fmuWfB/w5rLzVc1X1rZc/sqVdVVmrl2PFFb/JCzl8APAzlzmQCTg1F27K36ptHpt5QvVZfW6Q1GkN++Xs3TGhgE33LhgLHclAL5FBIkOIJmuHnjF6XBATfRpsDKXzbtfZaidoURMVmrt8lZHzYkwmTBRSVbQKAFhmIJwFIkShB+l0HpidXDj1ZFb1R/TnQPB4rT6LJ8aGMZAHIY+VjSwmqmy0yYnjvOKk6zguK47Tio123Clr7sQTHvPB18MAkkvHQBdMJn0jVgfhi/bs4quwqUx/emf5vwH+KfADlzqgM+yqylq9Evg9IvLjWHr5F++UHTr7b0x/a0dQICyUuEnObsaCahwUznDjqlUxoai3g4SJhY5kYAg2oHgkA0M0JtqGxJCsTfzoYyilBHTlMYRXtarqV413ZvCoNfCl+HM9Fpq4Zzbx+dg6dYFbJp9KZE5CrQmr9uWZzJ07cCBRumksqnXpJAthNFaasilGyz7XUehzIA8OotFYaDMmNrkzIOUVz+7ucFxWbGLHcVwfzMbm7ly4iySTvsHXneI3Vf2XIvJRe5b/B//MZY/noF1VWSvFkvq+8NJH0EjZsGSi3DDSDKSlGlfjoGpZy141ro4J2RCg0onFRpWFVpI5YmAIiVtawWOy9qBhAaDBU2XKyu6uqh4P7RvMVE89cpduN7NtyrpuflRhZVMcwGoo5DyBZuqIrQWhQ7U3MWKlBiivyzCxUEmTS1fduNAASfIeNhoC45HFRuNoQNqsOo7Hjs26405ac5JXPCudTG7dURg4aJd34S7c5v667dozEVo7BZ5pLKgBUmUl8QFVPCNB5qyEQxZ9Ls6hWKgOrPYhMmjPoIkh2kj96HLu2BmYpopCxY5vjoeYBjPzRK3t1AMsqPd2DMZNzC3kwaqX1o/WwVYqeNzNKwqj1Zibdp3jYZduLJQxImMkZAiDWAOJ0dlobOY4DZ4kOwSGMZDHQD9GtmNkM3Y8a20sdLtbcTuveTIdcStu93/ph2Xsq2pzf612zTNS7Ytty1st8+HmOCjrrpBQgaUXduN2WaiNhY7I5DBSSqAXY6IqKAwhMoY4u3LJmCiXOu+melNiKZVVQGgk7BrVTCzkxRRp5v1MW7WDpyLeKsW6YqsXK1nMJyoOmM6EBKkHVDKSV5YdPqZpqoXFRmFio+BAqmwU1pDHQF7bFI4yBLZDZFxH+iGxHRN3Vj3Hq47jbs2dbs1DZzLQXgBdZTLptdm1M9Bu3AMzEy0EA8IkJPQaOZpcuIu7cbCfhTLqY0OZwRNMq6Aw6NYnn9lxjFNGs4GoAn8o9SYQKFpMNDhlu+faxERTA1+mad2IIMNwjkI3T42wL8/9pew1u4sio73WUpDRpp6HMRDWViAyrA1EMoqBZ4AweOGUQSi9PY+DkFeRcYxsh8R2SJysOm53K26lAwCqaUWXswu1uX8Q7NoBBBjNN+wDzZjQIjPBQSRzHOQZ/IushLPcuH0stOwnVBaCwuDuXNbAGC0BswKozrXJLiqMWskgsNDZF+asO/tr0/sI4LG4eA24+lqDuVaTuCC6Py5Sm3Mktf1kVTnUK6eOhTCWaQatZHfpRiG4EhcGK981s5EgQ6T0gbIuDGOwJNR1ZDN0HK9WrNL+mggmf14uaL+IcCUivx2rb/1s4PeJyP+qqr/Jv68fwrKwHxGRNwGfoaqvFJHPB/5n4D2BHxOR71bVe6pO+mAAyE13gFRnpNb5/HPPoDkOsvR7JcppN+5QhnYnyUvcxum2XjPuIr1V5dTifVbnNvNdsEo+tbvDURpZJyuHdbIyharvE3kVyetI2QbKNpBXQlxD2UJeQ+yF0AvjkRJ7iNtA7ANxm4jbROgL0mdCP8KwRvoBGQYYszHSMCLjaMBxl05znh6SPPu7qnRep24SGEpBS0JGY6QyRGQdjJWGQF4zxUdlsMKmYY0xUR/Jo6CrQD8YkLbrRDpQF84mBN2XZNIfxVy7fZ/d23FBrRvd11z6YM6way7t6zfOvcwzF1aclDfEgeQulkY6KdNU79o7qObGHczPeoZbLeRY2aj2c53YaByRsYNVRxwLMiTCEMlDJKzcvRuEPHgNvLUJDXlwVhqscms+Esoq0veZIR320+7ChXva2APFQFOXhipTN0pc0UCvibUO8/iQzG5cj82VqTMn92VoD8McJw1qYBtUGeqzKoPaulGVsdRHIedCLoU82qPkQhnz9NAx+EPQEXQEsqJZ0exNt3YfRX0MSZ0ZFCnZZ7LOG8rOB1XNJZufdXoNPsZUB16LzmyUs7HRGIyNcoExEXJGc0JyMqUuR0IvhCEYeAZT7cLK2KisIa+MRfNaKb2gSTlk5xRQelrbAwGgZkD+YEbCHPM0cVDjxp1W484WE57RVktnNdPGJzZStS55uWGj3CFjsXYtQ0KGRFkFwhiJg7mgeRB7Xhkzhc6AFAahbG2G7l67Sxfu6WLXDqCqwu2XsttiInPsM2gi40zkokCrxlU3jh05oWu77FYBQgWwdJusBuDO2S+V+gjEEIglWDdwseo0gUhQe4jandsq8/gjeyOtKEjk9COotRoNSgjizxCCeFd5BQGZEtsAKagoVvLXfGD1RzWJEVRmubsEtIxzLORlh8kZWXW2bEzIaECSIROGRBkKZeXx0RAIKyGuhNwL49rjo60BaX8PLT+eGxfuqbGyA57qyo0aJzdu0OhTGxJZBgYNrKhTvJ2B5HwW6pp5/lkLo2SOtJBRHtWRjQ48FrZsNHJnygFb8UQ+4risebIccTvb44nhiDt5xe1hze1hxfGw4mTo2PaJvk+MfWQYArK12CJsAqGvYgKEXohDXZaIvRK3SuyV0BdiXyZhQYaM9CMyjMiYbcxoHGG0jg/WQmX0G0MG4nLcSMVculishnfO0CVk6uMakdXK/o6zUegDZRWJq0BeO5B6Y6MKpEMAqkm2V2ki8q7AdwDvC7wB+NQ6F21nu0sllYrIB2Jd6z4U+DJV/WvnHcsDUtp3zsZubR5ElbmCzJQTZ8s68ik3rmWhsyRtWMranaU5T2NDOLMdhXFKR8mhjsSEeYA3zZnk+2wQJUv0mQ8BlTKN81gLEyoJMrczsebDMYAGIUYrPB+iWKXFKMhgA7J4jyEJYm0kxQZcRbzNZDg9ZRxAtKC5TMUfLSYKXlU1ImNHGDJ0cRIajIkCZRWI3QykMxno6l24LwH+paq+WES+xN//+T3bXTap9G3A5wOffNEDuX4GagLMVsYey1yNp3XnisrkwtUk01aNsxiIhoXs7ntWLFQHV2nGhjoERDmiYO3grThGDnNy60UseBOsQayVmiVnWxXRtv/P1EjLAaNBCVHQGIhRIcrUvc5cvEAII1PT4jEvX4+jsczojJPz0q0bCyrB3b0CMdo2MVpc1HUum0foHEh9JHSRsIqUVbQ4aTAglQNduu9yIPU8ewE29wysufAPsAdAl00q9ZnUvyQiF27Cdf0AAvthmcHTzg3a7dQwNx62u39PpJNxUuMydoO+LAv5JBvriKDz+1r5p5ApMtQPTLaPddoGwG37RRFlDEoJ0duTBAdP21SYmXWaZfN7B1EqxkIpQApIP9oHh2Ad8mKwJl4hWF/WIAaQJouhunWameOjGL2nUUFSts57DZDoEnGIhD5ab6TegHQQQHA3yaTn2Xs0/VF/AXiPK/8LF7QHA0CN2QRLmR5jieQ4y9mDx0O9RjqNrGS09xKdPQrRWajzDIWLsFA1q8SZDZAyTS6guj3gU9suQEAVPLUVYx9MFBh9cL5MgAhTX9SJhQQ0MvdOjVCSsVFMgvbWhCvEgCRjHnHAkLI1Ph4qqwQYRstikOCum53P3vioMlH2/VQgDSPSJUgJ6RKyqvFRQg+VRuKgjH1mMqmI/AssY2DXvqx9o6oq+zqMPUX2ACSTzuV9981KrSCylJo4CQmzAhcYiKw0O+hmFhqwqcb7FLldqyyUtY2HmPyP0oCoDtxmCXuBFEUJUgzQbTNgDEhbUUax0lwaggHJXTSCgVMjaG+AMlcOQnQ2ikKItjxGS/EJwdwpYoYhG5BChMH7D0UDADVb4az4aB+QcjJ3bxxNxRsiMnRIikg/ovHAt3vYhTsrmRRV/bhD60TkF2uXbu/AcPEiNldsV1YXzrf7FOClwG9X1dfs2+Ysq4DZlbJ346BaYGTQxOBFEs2Vy5baAxMLXUSRq9a6cjUe6gjNFaBkCkeM9QM4Nu2t3whjc8Us3DmUGJyNRBlCJEtEgwNJggOlcd/i/BwdSBohOBtpEkIKaCzT68mtGypzeBpQCCZfO5DIPqC6w0i1n+vCtRtGNHrM1CUkJdtPZ0qdhP3fqQBnzPa+W3sFVkjxxf78XVf+Fy5oV1IXzrd7FPizwOUr3nvK/+lcuFmFq/HPUBLF3baqxhWRSUyIKJG8NxayW/vFQdRJZFDrFdp5rt0phz7gbeGbRZV9KM4+heSP6s7FoISgDEEZh0gJ5tJJVdmC90mNDRs5uEqEmKAkMak7NW5dKrNblwIy2FiUCQtxBtLoEnbDSHbjmDvpQU0LMklc1GIjY6bR4qXRXDrkgAt3f0SEFwPfKSKfAfwc8KkAIvJhwGfXBNHLJpWKyHtiNbSfBRQR+XPAb9xppbKwq6oLB/CVWK3hL77INzBZI2XDHJSfAhE2HmT5cHHBQL3miX3qYx8LIecLCtXaaQ9ocHdQ3KGfr4hpOltlISwJtS6vbeInt06UFIolrEphGxLbYOJClgjRAnINgZJmNirRphjEAJowtpnYyN26ZOk1IfnAZ4pI52NHQ7aLfRhnII3jgpGsU0RE2jlHEyNhg7FBjLlisVl41aULh2OgqxYRvFzAx+5Z/hrgM5v3l0oqVdVf4ECC6iG7krpwIvKhwHur6j8TkYMAOlTWapmNMMdDk5igLiaEkaFEcpBmenfYYR89yEI4mC7CQoXiIMJ8kEaZa0FUgRJRYyNa960QswGqunDJW8OnYIxUmagPhSEoeQiUGCc20uRiwSAGHAdSqDFREkKygcw4QEiRMBQDUudActdLcnZWiqcZqQWSiI0JqUypQG1Wg8VJZe5APoxnAugmE+EME5GAtdd70Xnb7pa1arWTChzgVNwzNuNBu25creUcSOeyEHApV65QptdWo/c0iJoJ1QDOLGqPUsFT6IqxUArZWEgqgOojEYPSx8g4KCUENAW0t7EnTaCDq3bJWKhU4CSx3DSfCBeHQEhqk+UmIBXC0DCSZ2EbkGzcSFKy56JTVsNUZ2FPnKS5hzqWdCjOuT8u3ANjV1EX7lGs/8oP+MDUewKvEJFPupCQME159re6v7NazTxY67hw4zrNDJpYSd7LQgMzCwGnXLmzuny38VBwebuCqP6NVuIOYjUIAkoIdT6RAeZYsjOQuZt3wpoUMikUVjHThY5NKMSY6EOxyXtDpARFk1JScDaCMrgKl5jAFEbMhWtBNeoeILlrNzobjRHGMrNSjJYNHj0jIRu4tJQ5ThKd6zP4oCxnTB15ZwfQmdNrVfVx4Dn1vYj8APBFl1XhtB1MbZS4sQRKMBeuhHHuaSONCqejjQtJWDBQ7zzSe6p3V4HUuHJRwoVAlLVMY0RTTCR1GyVINnC6MleVuKB1op73JCojnazpQiFJdiZaTWyUYra4KCaGWBhitP5EKVKSoCnYzNTBRIOSDEgGIMuSDoNMbBQ6IQ5qyyuQfOKcjObuyVjAx4xk7OZMBp8KTmWjCqhSkBznGCkfltnuRy7cg2RXVRfuSqzNRFhI2E1VnqliqYZG0k5Wz0DTFMDXx4AQRRcpPrB05c6beDcxETTgczDWq8OZLTIQQ2Gj5pKtdLRjEjU3zo+vk0ySPNXhvpNXrGJmMyaOQ0eKhT5ZJaB+jDYPKcYFG5XORANNgnpC5y6Qwji7dmGEmLwOwhiQoRAGtSneQ5wysSW74FDdO4+fFiCqokOdPp7zQbHgOpJJvbnWy7BbWgf87aar9/cCz8Ou/x8CPldVs4h8JSaQFWxs6UWq+pazjuVK6sLtLP+oi+xz8n7qAGqTElNBtMtCgwaSRjotjRJn2Qi9Rr/bWyzUa5wC90F9VHQSEMLkyl0kHoJdYWGWuOs4UaBR5MhElE2VsoPSqYNn6o63pgsNG4XMKoysQje5ddsxsQmFLlkRjyEqeQyUaPFRGR1Io6tvg1DGCh4DkozGQmE0QOWVEgdx4JibJ0MwUNU4KZcpTjrFSrlMYJJc5lgpBmt1tc+unoEukkz688DvVNWtiDwCvM6HX96CAe4JsZjjpcAfwtpEfpWq/i8ALnX/ReCzzzqQByaVRxdjQUziQdfWR2jL/Pqy2iw3eFWdIq7IiVriKcEkZCz+qfHQQBUF7g1EhbKj7s3sF8SBo15jQYq7cBYHrdV6tHY5e62FFSdeb2GTOzYx0cXMNqeJnYbRymvlMVCGPUAaxGKkhM0gHZyBKrBqnYNaE25QB5M6YwVkVMIYEXfvJDcy+JiniXhWY6GbXbw7+37Y60kmVdW2VOqaJmekGddJWIdv3VkO8HBdfpY9GADamdbQytmtEmePyFjKXBNB814WqgOZUJmhTME/uoyHLgsiYJGx0EkVJ2Zxoc3Lq4+OTBesgGMnI10ZWXnRkm3p6Mrs0p2EFavQOZgSm5hYpcRm6OhTPA2krj57YUQHU0ggq9mFk7EFT/N6ZAZSy0q1gs+oOy5eN7t4XTLZ+4DdBwBdKJlURN4b+GdYW/svbt0xEXklNsb5PRgL1eV/GfhjwOPAR593INcOoKnh9u6AKjtxUE3l0VqJZxYTgpaJhaLq5CpFD9574rRuABCLh6htOfx1IF54GviuQtfGReY6LsHTstFKMysdWUlmU4yNNtKxDgPHZcVR8HaTMXGSV5zkzhgpJ9ZxtBYko7Vd2Q5WLTWPgZwimgNlEGQME4hknOOhCiBbBqHWyV6ASf29F6avcdKos4s3FlPvsrt64+GqPAcAdN+TSVX1vwIfLCLvBbxcRF6qqr/o6z5eRI6AbwM+BniVL/8y4MtE5EuBzwO+fP+JmT0AyaRKjdDn9jYzE2WvST1KYJTIKIXR44iaE1dZKFLLUaUZSCgwEhd5bjUemkWF+joglwcRLFy689iol0IoTadwHbgja450MDaS0cCU185IK7ZxWABpExND9i5yuTJSsb4/7taRDUQMQskOpgqkvMNKGSus2Lp2uQWTMZKMLZjKAkz77IxcuKcsmVRV3yIirwP+Wxq2UdWNiHwX5hK+audj34bF/Q8wgBo7mI19QJEbqpjg7lzwOm7RZePWlYMmwBd8/VJUAEwUqK3cLwMi32916eqg6y4bdSiDKJ0Wc8000ZHZ0NHFkU3pZjbSjiMZ2WjiVuw5zquJkbYlzYzkjZD7bGpdn63j9jhG8hgpXsOtjIKOYXLvpJbznVy6hp3y7NZJZaTpeT+YDrpwytQP9grt3GRSEXk+8CuqeiIiz8a6d3+1CwqPOvgS8HsxJQ4ReX9V/SnfxQuA/++8A7l+AB0YSJ0fO/URSmQQG0MZSiBFYyEDkLlJJiSkhSvXqwGjStt2wZuoMOyAKIpw0ZioWmWjCGQVQq2SqrW5sU3O2+fWrdRUwy5Yh4iNdGx0NQGpMtK2dDzkTLR1MPUlssnGQpuUFu5dzlYUvuQwx0orK73FWIHUAqWC5xwwOTPJDpgO2TUlk/4G4K+7eyfAX1PV14rIe2AD/VVY+H7g6+t+ReTXY7rhz3GOAgcPAoCquaTNbi6cv65uXJIyiQlBdWKhOlVg0DhlR7euHO4t1liFaZmBqGWifcICXJ6NamzUKnWToODgGSj0Ug4CadDEJnRsSscmdGyn58Q29Rznjr54Dx9npb6zvj5DiTMrjZFx5aLDGLxWnYGpAumiYJpBpd5ryJjpkF01gC6STKqqr8K61+1u84vAbz+w30+57LE8GACaQKM7SpyDRwNRy6TGJbUW7ZWFQjAwRS1stJvm6FRXzpTK3hhCZzC1IKruXBFZCAvtYOtZGQut7cZGrVJX3bpBi81acHXuEJAGIkdqYOo1sq1gUgPTQxVM0VrUV1bqO2uK1XdxYqWhWB/UYYzeI7WCKUwgKvkAmMqOAJGZGnVVEB36bd+pMxHuq1XWaXy4KiRMqlsRYpAGPIWx2J29slB9RC2TKxe8nWP0NJ7gp1pFhcMgwhW5WhvBAFB8XtDdsNEhty47Ew26ZKQsmUEDR4wMGthoZ32LNLHRjlshsSmdt6N0ViozK22zgykaO5mbZ26dsdIOmEo4zEwToGb2mZgpN6A6UEBeuC8x0ANjDwYDwSIWmhJKi1BCW963MJZg2czOQqHGFKILV66KCpWN2mkGZqdBVBxExceJBoz1ArqIiy7NRnAQSJWRSuPazYAqDAQ6KQxeUPJIBwbi5N61rDRoYqOJbWwAVdLk4p0HpjFbf9SSxfLvsgGqjNZLKI/eZ3VwUFUgOTsd+l3Pcu+e7nbtAJKGhSYlbseVyyUQQwuemYWCBnvsuHITgKgzRt2NW9gMoox1qgP1mafmXh2Mi+DCbLTYZg+QOmrxk9py0urb9RSKFJs0KIFBR5fuA4MYWHpN3gTM2GnQOLl4g8YlO+0BU58jfYkM2Vy9MQcbX2vdvBwolZkcSFRmKiCDnFn77caFe0rMr+SdPLjaAS4XExBqbpxNHXAGEt3ryk021S1Ygsg6fWePe/Ic/O+IC0WE4IJGZaNwF2wEy3SgCqSalBqJU5+j7GNHBSb3Lkth8PSkXgpHOtBrZCBOKU3VxatAqm5eC6ZR4+TmbSuQSprAVNlpzPbamCnMrl6eXb2SHTzjWVO6r9aFu0Rl0gy81t++UVU/yZd/M/CRWLYBWNLof/R1HwX8TSwB9a2q+pFnHcv1x0AssxFqkXn15dWNq2KCFB9bqc8aZiA1rhzgrpuPA02A2mGiZQjmQoOecukOsdFlY6NT2zSMBOyNk/axUiYbmDwrvU4sbMFUm4RVAPWa2BZnqVhBZSAaNJ7JToM/jyVMgJrcPQfUPhPuCwNdtDLpiap+yIF9fLGqvrRdICKPAf8n8Amq+kYReffzDuTBYSBX4kR1Ma27unBWEmpmoYmJOO3K2fRMv+jrtVraMlQrMiOWQ91P08JtekRhRaZOeKgu3VlsdLdu3WK7RrWbUoP8KFpW6kTJ3pKlQ93Fk0l46LEOeqfA5DUk+tp5T5fMVBmrZadBgzNTYtRwLqD2/656P0SEc5NJ79I+HfjHqvpGoFYqPdOupKyViHwhpr+PwC8Df1JVf+5Ch6zNoy5q4qCiQlCmSXYLFlJlbBlI1FPnR0KIcxp9+9sugASwopOx+fs24zWTrVywB/RFvY8qBsxDbt3dMFK7XXCXsJY+qbGSHfoMJgOP+nvodclMxkphjplccCjIAlAtO9XpIfsANWhgWxJjdfcahhqLxU2Hft+zBlnv0i5amfTI8+1G4MWq+vJm3V8Wkb8I/EvgS1R1C3wA0Pmk0EeBv6Wq33LWgVxVWav/AHyYqh6LyJ8B/irwh8/b9ymbWMgZiCUb5WKVNYNY46vggKkMRE4QPTO6xj0NiDKNm1GwBl3iPYakFimxeKhloyowtLFRdeuumpHabZex0nJgtgoPu8xUgCyFrDBIoSj0Eig6nGKnrMFZysA1NCCqgNqWjiHFKXl3qx1DcTD5shpDHbLrSiYF3kdV3ywi7wd8n4i8VlV/GvhSDHgrrEbHnwf+EoaH34YN0j4E/LCI/Iiq/uShc7uSslaq+v3N9j8C/JEL7Nc/LEsGmlw4T+MpgSxYPbW4dOlGMRctFCUQ/bXfCSt4/LlYTg7AVNHU8upqhdPRX5vStSKTp4FVq3raMc1jpeN8ILWMlOs8IurhXS5WWjKTZY3vMhM+RSPvsJOxkj0XhR6bV9UCqnY/r4CyGhSngTW5f85Sg8aJpQ7ZARfuvieTquqb/flnnFV+K/DTDXttReSbgC/y92/C8ufuAHdE5AeB3wLcE4DOLWu1Y5+BzbE4ZYuyVs9+9jJ4h504SKz+eZxZaI6FZlcOPNbZZxU8IUxMVLyouz2MebKz0C4bBXSKjSJK9pmmlwWSbTuzUlXu4N6ZCczNS1LLd51mp1aEyKKsNU/AagFVxQh73u/yzVVhl27foPtdOLm+ZNJnA8c+I/U5wO/CPCMa8AnWyuR1/rHvAr7Wk0xX2HX+1WcdyJWKCCLyR4APwyTCU7Yoa/Wr33v6VkV3lbjqxulCSBDJZA2TS9C6cO1za0UEa0/SmQxbARPmMsGH2Mguc2kAtARSoHZ/8PJZ/twh/nchqBBp3TtT90L17Rqv8tLCQ32/h50qoMAZyoEEBqLCkqGyf3U21iRk+oMsVUE1sZYD6ZBJvnIAXTSZ9O+I1ME7XtyEHd8mIs/Fvv3/iCeNqupPeL2EH8O+jm9Q1ddxhl1FWSv84D8O808/0gOyi5nuvK7DQVNyqU3NFlVKCYhYB258LIgAY3G1zRqyzWyjMr0vWsjBywCLJ6qewUaDJGubQqDXYl0gdoBUXbsqobdigyWNyh5WMoHilOgAdwWmvdv7fnZjJ+AUQwF7XD6LoVqWmmrtTa1mToNq/++rRs9XaBdMJv23wAcd+PzHnLHvr8Iac13I7rmsFYCI/Fbg72D6+d1Vyp8yEmpamg2qliIEYQKP5VwFQsjksnOROYhKrZIZgJwcONaAq4hQQjiTjaxhVz4IpF4jKy9csvJ6b5HiKThlwUb12cCEA0e8hLHl3LXMZId9dzHTvm2DzG6f2cxQtnzJUnVZZSlwsKEOrGyuIDIBq45FHfpdDyaaPgPsqspafRXwCPAPvbjiNOp79s4r3ey8ZpayaWKh4uCBYlJ0mVIMZmuUN2OZsvMsB9lokDgV/DgEpLpsIHoWtc2IjQ6cHmVFBZVeAkycipmuip1OfabZV8tSu24f7GeqeXllq3NGSt/Zk0nPK2t1lmJy/s7bfXpunFTw4Kk9ihaxetEeF+UiEE6DqKgnoPqFWZNRS7CmWfvYaMgOnJBPAckAE70clYNHi4FJEsFZyTIHbDrCIIHojb0uCiZwQHnMdB472fb3zlAwsxQsmQqY/patO8xWZ5mUcwD2NLYHIxOhyUCocY9O8nbNjwPUXLmFPr0Lop2xn5IbQO2w0YjXdpN8CkgDyeu2WcHGTkZ6SQYW1EGlEytFUTpGBiyhdSWZoGUBJptyvnTz4DA7AQcBtctQ9+LytXZIoLDjXLIVzKA6ZKJ6P0SEB8YeDADtMwVt5phMkJn/YxdERePEOjWLIYVi+XQNG41iGQxpD5AGMZfMpltHb01i7BPE6i9YmV5jpUBHJ6MVk29eb3aYKarXRpAKotwAy+tsKwSpBfFNHo9UwMwMBfP0jH0sZcvnu/5u5dV7ESiCnP5s4TwX7moZ6IKVST+apQT9gcALVfXl3jfoUV/+7sC/U9VPFpF3Af5v4Fdj2PhrqvpNZx3L9QOoTeVxd808gjY+mtfhrlwLIlWlBCGGAiWY9N24cUG0AU8wUO0AaVCr+hOkMEokhcwokUGsomgXzBXrwuj1DKyqjnVf6CZmitJNyzfSUWti19oM+wDVTSW4KogcVKozoGhbqZx2+2x5mRhpl6ns9dWy1YU+r/dFxj43mdQH9z8EJsC9Hvjnvm7qGyQi/4h5HOlzgR9X1d/nMvd/EZFv2ynSuLBrBZAobUXfaTxoEhRkdt0sEA0QygJEqtZfBwqqAQ12gbZsFLxKafu6AqkCKIXMiDpwlFASKWQi/qzWWSGW1DTNqt0XClu6Baiid6mb2cnanEx1smVcAoolQ7UuHziotBYjmcG0j6mggkywMsN13QyuUXNz4TdiRf38Hqa5K6CpenfwK7UXcLlk0j8IfI+qHrcLReRZWE24P+GLFHjUB1gfAd4GnDkd8PoZCGYVTnRmIWSa26DF03D2gEgACqgGogOpSCSGYhPAVBBnoKyBKOUUeEYJJLUsh6RhboalVha4LpvAJMqWuXXJBJaGnaKDJzgL1W3iBLJuAagZPPN2h0CFOis1TAUs2IpDwDrgBtZ11Zbu4H7mmtefDaz7wECXbXP/QqyH1a59MsZkT/j7r8WyHN6CuXh/WPVsifH6ATQBByMdZgGudeV2QaRe7xobU0V8XKV4UXk1AvNMAZ2AlMVmt06zXD05NZUyTdJLwXv8OAsFnHVCJpQ0tWqcAWQXeufstAuWCqgKmIldmu1Cw0zAaWBxAFi0DFUmYAETY9n6WtbY3cFpsLWJsRZMlXfWy+K9LQsEhOFgdy37TdlfM+4paXPvuXIfhA3D7NqnAd/QvP94LDPhY4BfC7xKRH5I77FH6v23KQNBGr+ued4HIh8vETHwSFC02M+s4TSQsrtKIszgqcwiM5AWYKqAcsD0Ze74kEKmF6/BfQBQFShnsdRFQNW6f8BBYAHngmt+r/P2FUxaGqby5z2uob1fgmzXBVz8uPtFhKeqMumnAi9T1UXvCM+P+3Dg9zeL/wSW8qPA60XkZzHx4d8d2vn1z0iV5rU/G1waV24XRNJ8TgQVEx5CcBXcgVRECVVEECWLTG3mbWqEX0zBumdXcPQVSGHubToBagGy3fWVNWZQRQfbBKqyOgWqKGUHUPVCPx9cwKn4aiqqfwbA5vUNY7EEWd2+ZTFqHOZffzyEm/Y3PqPw/F3aucmkjX0aNn1h1/4g8E9VddMseyOWIvRDXoDx1wM/c9aBXDsDmfcm1V+bWEew2EdFrTN0gzZtgSdOMSJkVW8jb+6aBKyqjNTlUEp9be/NnZvBVIEkogzFwCByGkDGQmXxvoLoFKjKaVAFv6DPY6wWXLb9fuYCTgHMlp3NYMAEsrqP6W8ttlkCrS6r53DYFMqViwgXbXP/vlge57/as48X+n5a+0rgm0XktdgV9udV9a1nHci1A2iyxW/QgAixDtFiF/kEGEAz1lFa1Ly/YNuLtV4gOKCogIIFmIJYelD7Hub2jPtAVcFxCFT7lgGHwbXYZgmwCi5gL8Bs+dkgA84EWt33LpvV1+3fPw9wB3/XQ50b7tIukkzq79+ATcfZt4+P2rPsLcDvucyxPDgAgjNBBOocJLPCIDuMtPNcc8xoQGOAYgJUZSOYATNKWIBKpotI94IKOAisdt3Zy8vBZa1bWLdpQTZv07xvX+8CZg/Y5tft31gCbn69ZLezf9P7ImM/MHb9AGpintmNqyt9jeKxDqjUD+xhpFpfrmEl8NS6ymBwGlD+pyqo7DN6ipnaZdAAiSUAW1BVQAELgNX1YWf9xdaVg8v3MRqw1220z58GzC7DARPL1X1Nf695fdBucuHuv4nOTLMwlWmmXQskWzWDZwEmNcDUfU2AAgdKBSynQFX3tQssf9lscxpILNYvgbcLsmn5zufCYtsloKb1BwC3b/28rCy32wEecAp87fa7AKzb277OAIie3cX76W4PBoDc5ZpBVBe228jknlU3DN0BExVsFVD+2QYwtufT4JiBOK9rFp8C2P5tlwy2fLblu2x2aB2wAN3u+n3AW37mNLD2PnO3250G5EG7ehXugbGrKmu1Br4Fq2jyK9gI7hsudSQtiBwcU6xzKsapV6xfuzvgqR+ddugv7ViXy1tgMX1u+Rn73GmwnN6mvtgPnjP30+xqF3ztPnaBs+/vhD3rdsF4aF+HgLncx+ntD9o7ewx0wbJWnwG8XVV/nYi8EPgr3GVZK1hckzMjtWCaDq55IadfntpZ81oXV7XubNO4hdP2y210sf3OCewsl1P7P3+dLI77NCB2tzm8n+UFvtjtBfa7C5CLfGZhquhwtdXlReQPAV+B1T34cFff9m2398bvs6u/HXg34N8Df1RV+7shgispa+Xvv8JfvxSrbCI+onvQJrbZNV1uM3t188ZTksLuDk/9kT1/4BTI9m1zauenLv6dQzp4DPvOsf7Jxda7n9333ey9UPed976/ef7+92bFnHH855nC/YiBXgf8AayMwF4758b/V4CvVtVvF5GvxwjgJdwFEVxVWatpG58C/jiG7jMHoWB5zewFE+y/PvZCcz9YLmrLv3+BDzYxz8HPXeLvX2r7y+x3clfPt4O/wV1udz9EBFX9CQA5G8V7b/wi8hNYrlut6/H3sZv/S7gLInhKRYS2Lhyw/dk/+0Vnlgx6hthzuMCN5Bliv353wZO8/ZX/onznc/ZsW8vuVlskk16BHbrxvxvwDlUdm+W/avczFyWCqyprVbd5kxelexfMh1xYWxdORF5zVjLhM8XeWc4T7Fx3l6nqJ9zlvg5mY6vqd93NPu+HXUlZK+bkvh/GkvS+77z458Zu7Cy7p0I1Zodu/L8CPCYiyVmoJYQLEUFr504x9D9Sy1r9BPCdtayViNTSVX8PeDcReT3whdgU2xu7seu06cYvIivsxv8Kv7F/P3ajh2U2dyUCuCgRWPncp/4BfNZ1/e2b83x6nys2h+dNwBb4ReCVvvy9gO9utvtErDD8T2OuX13+ftgcn9cD/xBY+/Ijf/96X/9+5x2L+Adv7MZu7C7s3sqx3NiNvZPbfQeQiHyCiPwXEXm9lyDaXb8Wke/w9a/2SVBPO7vAeb5IRH5ZRP6jPz5z334edBORbxSRXxKRvUMQYvY1/j38mIh86FN9jE+p3WdfNWL+5/th/Vb+E/Abd7b5HODr/fULge+4bl/+Pp3ni4Cvve5jvYJz/d3AhwKvO7D+E7H+UAJ8BPDq6z7m+/m43ww0jQarFaeraUCtvQAbDQYb/f1YOWeI+QG0i5znM8JU9QexemmH7AXAt6jZj2CS8fOemqN76u1+A2jfaPDuFNvF6C9QR3+fTnaR8wT4FHdrXioi771n/TPBLvpdPCPsRkR46uyfAO+rqh8MvIqZdW/saWz3G0CXSQPioqO/D6Cde56q+is6d+77Bixl/ploF+po+Eyx+w2gvaPBO9tcfvT3wbNzz3MnDvgkLKvjmWivAP6Yq3EfATyucxneZ5zd12xsvVh3u78HfKunAb0Nu/ieVnbB8/x8T30asfN80bUd8D2YiPwD4KOw0rxvAr4c6ABU9euxRmyfiI3mHzMXbn9G2k0mwo3d2D3YjYhwYzd2D3YDoBu7sXuwGwDd2I3dg90A6MZu7B7sBkA3dmP3YDcAurEbuwe7AdCN3dg92A2AbuzG7sH+f9jUMnNbP9sjAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "from dune.fem.plotting import plotPointData as plot\n", "\n", "figure = plt.figure(figsize=(3,3))\n", "plot(uh, figure=figure, gridLines=None)\n", "plot(iuh, figure=figure, linewidth=0.04, colorbar=None)" ] }, { "cell_type": "markdown", "id": "electric-taiwan", "metadata": {}, "source": [ "## Couple bulk to surface\n", "\n", "Dune-MMesh makes it possible to compute traces of discrete functions on $\\Omega$ along $\\Gamma$." ] }, { "cell_type": "code", "execution_count": 5, "id": "union-independence", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4.266679479547976e-08" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from dune.mmesh import trace\n", "tr = avg(trace(uh))\n", "ib = inner(grad(tr), grad(iv)) * dx\n", "\n", "iuh.interpolate(0)\n", "ischeme = galerkin([ia == ib, DirichletBC(ispace, avg(trace(uh)))])\n", "ischeme.solve(target=iuh)\n", "L2(iuh, iexact)" ] }, { "cell_type": "markdown", "id": "invisible-certificate", "metadata": {}, "source": [ "## Couple surface to bulk\n", "\n", "Similarly, we can evaluate a discrete function on $\\Gamma$ at the skeleton of the triangulation of $\\Omega$." ] }, { "cell_type": "code", "execution_count": 6, "id": "inside-workshop", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAANAAAACSCAYAAAA0NEQTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAz/UlEQVR4nO2debwuR1nnv0/1e865uQkhgQgKASMCOuCCEAVHURCIwAxEBDG4IIpmEBhZBMXR0YhxJoIjysAHjIAsnw+bYcQgmxBBkDEMSYAMQSAsEQMMCGHJzc255327n/mjqrqrqquXdzn3nHvv+/t8+ry9VFdV96lfP0s9/bSoKmusscZiMHvdgTXWOJaxJtAaayyBNYHWWGMJrAm0xhpLYE2gNdZYAmsCrbHGElgTaI09hYg8WEQ+LiKfFJFnZY5vicjr3PH3i8hZbv9ZInKziHzILS92+w+KyJtF5GMico2IXBTU9QQR+b+u/D+KyN2WvgBV7V2AlwFfAj7ScVyA5wOfBK4G7jlU53pZL6oKUACfAu4EbAIfBu6WlHki8GK3fh7wOrd+Vm5MAgeB+7v1TeC9wEPc9qlBuYcDb1v2GsZIoJcDD+45/hDgLm45H3jRiDrXWAPgB4BPquqnVXUHeC1wblLmXOAVbv0S4AEiIl0VquphVX2XW98BrgLOdNvfCIqeDCwdRTBIIFV9D3BDT5FzgVeqxeXAaSLyLct2bI0TArcH/jXYvt7ty5ZR1RnwdeDW7ti3icgHReQfROS+aeUichrwMOCyYN+TRORTwHOAX1v2AibLVkD3TfhCWlBEzsdKKWRz814bt73NCppf42hj51+v/7KqflO475z7H9Sv3FBG5a66eucaYDvYdbGqXryibnwBuKOqfkVE7gW8UUTu7qWMiEyA1wDPV9VP+5NU9YXAC0XkZ4DfAX5hmU6sgkCj4W7exQBbd7yD3v4ZTx0+p1NY7x3kBA8f/MxTnvEv6b4v31Dy3rd9c7TvlNt9dltVz+6p6nPAHYLtM92+XJnrHSluCXxFrSFzBEBVr3RS5a7AFe68i4FrVfVPO9p+LSswN1bhhRtzE0ZDJV72I1bWv0r2z7LsPUGZahUtI/AB4C4i8m0isol1ElyalLmURko8Cvh7VVUR+SYRKQBE5E5YG/zTbvtCLNGeGlYkIncJNv8DcO0815jDKiTQpcCTReS1wL2Br6tqS30bg/1KmD74Po+WSisYrLuCtF9mPjGrwJRRpGnOUZ2JyJOBt2M9ci9T1WtE5NnAFap6KfBS4FUi8kmsLX6eO/1HgGeLyBSogCeo6g0icibw28DHgKucv+EFqvoS7Dh9IDAFvsqS6huMIJCIvAa4H3CGiFwP/B6w4W7Ai4G3AA/FurEPA7+4SEeORfLMhbHE2W31cOx99v0dSaQKZXuc1Imgqm/BjqFw3+8G69vAT2XOewPwhsz+6+m4SlV9ytwdHMAggVT1MQPHFXjSMp1YmjyrJt8Cg1ilRwoNkedo2lRpWyu6d6owPQFtw6PqRFgZdltahfXPMSh6SZQ9YaAbKxqQvQ8o30ZXmUpGSSFFmB73akQbe06gue55T9lV/e9ag3ZBMtXISZ+OenbLu5fWm71XfUQaQSIFtrVYoHfHNvacQKPQQY7deOCFdWbJtOwgz5zfS5xl28vcI99eJ5EWuK8Vwg4nHoH2fzBp8s/sdSGLLr9k2urrz1wYSx4NlmWhdNaXueTmnAWamaqJljFYNJg0OH5HETkkIs8YqlNE3hsEn35eRN44/5XG2FMJNChBMuSJj4/4T4/R/7vqdA22XNU9kmisHdQq03XOqiVQRlXLSqRUEg2ocdYGmm84uXmcFwIPwkawfEBELlXVjwbFHg98VVXvLCLnAX8E/HRw/E+At46pU1XvG5R7A/A3c3U4g/0rgYJ/XksSdD06JbMMtdFXPmkn7sNA3RDbP0F3B8nTI4FGCs/hujok0qJQFbZ1Ei0jsFQwqYj8BPAZ4Jp56hSRU4EfA944xyVmsT8JlJAnPpb+1+kY/HMsufM62pybREPQZL1DzeolSUfZLFHHSr6hY62iwo5OogU7d3hFsJyfnLZwMKmInAL8JvD7C9T5E8BlSXT2Qth/ToQu8uSI07e9YJstFafe1kilW/RpLSlhMuvd0mTERSYnd6pnEF9joKYu4pyxNlDLifDlgVi4ZXAB8DxVPdTzdkMXHgO8ZBWd2H8EcugkzxjiLEKmnBs3tAOUPInGeObc8Sx5elW7RUZyco6rdNDO6SLRSK9chbCtm/P2duFgUmzY2KNE5DnAaUAlItvAlX11isgZWDXvEfN2Nof9RaAuz1ru+CgijRATtYcg3BfUma4v+oiO2kx+ScmVek8WaKMmQOwBaRGpg0TzwjoR5nZj18Gk2EF+HvAzSRkfTPpPBMGkQOgQuAA4pKovcCTrq/NRwN+6EKGlsb8I5FD/c8eQp0/NG4PcLGPf4ApINJcUGtOFEcTpu8Ssmgauf20i9UmZeZ8TqvMTaMlg0rnqDIqcB1yUPXkB7EsCtSDJb2t9wD4aC08OyBMpQ6KxaEgS/2bJM8Ye6msjrTKSog2RsiRaUAopeMfBfOctGEyalL9gqM7g2P3m7mQP9g+B3D+tJX36yNNjG83z9IzmdyBPpByJ/OF5PVld5OlU6YavIUIsbDrUNWmTKK1mjmdEhXCk2pizo8c+RrmxR8wW31FE3uXeT79aRB66kt6l5AldzCHBAvJ1zhn1TJ60zutqJ/wNj6X7u5AjSEKeyP2s7XNGBVJofG5nna7tLumY63f3pVkVLlxOBIx5H2jMbPHvAK9X1Re5XFtvwaYdmgtZ26fuSLgeD+ped3d6btRgu7wSDKpUXcs5FRiWQt2u6zZ5sttpnzsbis/RpL+RRKqvL5FEfSpcTzTCIjbQ8YAxEmjMbLECp7r1WwKfn6sXfa7orDSx+2OpoTGxwsUXEbtEbaTlXD113fXSlnh1+TFIpU+OPInUqI8nxzrj20ZKnrbEyfwD5lQbFWG72oiWMVgiseKDRORKlyjxShH5seCcn3aa0DUi8kfB/h8RkatEZCYij5rvCvMYQ6AxM7sXAD/n3lh9C/CfcxWJyPl+Vro8dFO+tZxalLF/eolDQ5ZwafrRcTyjusVqXdK3eY1tJSuNJDfoBwjTqbJ1EKZTuuXUSvr3dVza3CpcoN08BLgb8JhMttA6Fg54HjYWDuDLwMNU9buxbu5XuTpvDTwXeICq3h34ZhF5gDvns8DjgFePu6phrCqU5zHAy1X1TOzr3a8SkVbdqnqxqp6tqmcXp5wcH+v0sPlfzZPHl2mRYYTBkCwtIgXSqGlHGSLOoOGt0kmeFhk6ulvfntyx3BLcrvzEbtDpOaWPPUWYaREtI7BwLJyqflBVvaZzDXCSiGxhs5xeq6r/5o69E3gkgKpep6pXw5zJG3owhkBjZosfD7weQFX/CTgAnDF3b7qM8i7yJMSJRxILLBkiubbajon4PA3LQx1IOmgbZYiUlRxzICdx+iTPKrx/qjCtTLSwi7FwSZlHAlep6hFsbo7vEJs7e4KNe7sDu4Qxbuwxs8WfBR4AvFxE/h2WQP/GImjZKM1/sUUeAtWsz1tmzx5orF1WkOAsRRF7xtiBlahTte2jHeTJSYZ52rOdblY1f7r3HaDxOpGnIa5raH7IurFbw2k3Y+EAEJG7Y9W6cwBU9asi8qvA67CS5n8D375b7Y9JKjJmtvjXgb8Qkadhb/XjXLjFYmgZ6j3k6SRO0nz2n69xsVZQYk0ZV0yTPbEjYBGE5OkiTi+hpPtQOB8Ue94C4oTHyW+Pg4xV20IsEwuHS2H118BjVfVT/gRVfRPwJlfmfKBklzBqInXEbPFHgR9aqAfS889K7Q23LRkytYgTOQ369CjQcKax5kWwj1AaORK1BnH/qMtJn07yZLaz6347vNa0aECSaJ/WP8mlCiljxxDKqnBHLxZObN7rNwPPUtX3Rf0VuY2qfklETsd+3eHR83ZsLPbP+0ChXZFAA+O9RZ76HD8SceXULmbYHhKjdfmovsTOabvANdvf9gW0t7PkyRCLKthfNfZNtARlUudD6uqWju3sMyYp24cKYacqomUIzqbx2s0/Y+cSrxGRZ4vIw12xl2Lf//kk8HTAu7qfDNwZ+N3gNW2fbP3PROSjwPuAi1T1EwAi8v3OU/xTwJ+LSBgjtxD2TyhPiGBwam6gtsjmyRVLn5bkyQ12bY6pij1HrfKG+naDR7TY53YkhdyjPDuhOiBZQvK03NIk+3oQ3opIHQtPD1W21olxQUHnU+NUmM0vgZZJrHghcGFHndlchqr6AdynTlaF/UmgDkhIptCa9+RJiZM6GrKV2h9VIvLU6+qJJDV/vCrnydSKUsghmDhNpUBOWrScDEHZXP9D4kimaKquhSSKzKO5bZ+mrdnIRCLHE/YfgTqkT5Y8IWFaJPL1NcMoZwt5+0f86BEsYTyRnLRBFYxA5faGEickT4ZIkhIhUbFa6lwqjcLf1gUQkQLp7krL5lkhFJhVawLtHTrsn3yZtsomKdFojvdJIF9G1Z6nKtbPFhCJKlTpvAUe/NYjNtafcpIjJU0necaocRnpU5eXiFttUkmPFCIoMBKqMsruOd6wPx4ZOcmQkz7hKf6YaRwFvpyIYow1+o1pnAl9iwkWMRW10yBwMkSOB2n3MwuFUH3LGvSKdRC4Bb/u9qfOhMjBUAXnh+U13s4RtePWt/s3Al4ChcsYLBELd2v3BsAhEXlBUP4WgVPhQyLyZRH506TOR4qIisjSc1R7LoGyWXe69iXesVTqSOApy9lEfWM8dCKIOuHj9DStxJK0Sp7OofQZgAQDN5IwweDPec8GnQiBGeZFiL9NPjqiLV2Sm6FRVe0bpZl9rSJCOacKt2ReuG3gvwLf5RbbD9UbgXsEbVwJ/K9g+xbAU4D3z9XZDuwPCeQR/pNabmsIVbcu8kTSRhQTLqbKL0m5UCqFrvBaEhG3Hbm7c8gQQlLy5CTQvIs25w/ZVn1SaA7NrblEtfNA4TICy8TC3aSq/0j8Ccn4OkTuCtwG+6Vujz/AknAlORH2lkCSX8/GnbkyoSrVRR4j9q2ehjSWDEXHMilsmaKICSWRWhe0W6/bX+3qO/nBGKpVNXkSEoVk8sfTJVTt6vJd6lpO8oUqWqquzU0qK4HCZQRWFQvXhfOA1/moGBG5J3AHVX3zyPMHsecqXAuhMZxKH7deq2qBvVPvp5FQIpZEjTOhexQYrAqnKlROfVMVKmObNjQhvFo1KpwIaE3yhD3BEz6rtgXkCO2eukxSR9+9gsDHESxq2qpb6DxI64j8IHMkwFOFsv0lijNE5Ipge5UfGR6D84CfB3BvB/wJ9nWGlWF/EKgmSjji0mOx6pYjj8kQJ9yGbhJ54oBiFKrK2G1n91QmIJGzCZoIbOlW4ZL5n6wUCCVOQJzUG5e9bZ4oTiL7CMTa7qliEnnyCA0/NCRVaDBBTKgeKFC254GGgkmXioXrg4h8LzBR1SvdrltgbaV3u0SM3wxcKiIPV9UrOqoZxP4gUBekGU2p6tZIoZg8hWkTR0QxxKE4nlBVEFGpakNSrDOhQirBiNgna2WojBvUxo4y1T7iJOuBquRtlZa9E0inIW9Z+I6SJ0/YF1UQkyGRBhIpVN8kXh9LnAbCrJzbIlgmL9wQHoP9zD0Aqvp1gldsROTdwDOWIQ/sBwLVUsf+1OM5HJj1vlB1ayRSSh5jNCJOYRrpY5KRGJq6Va3C2YnSUoxTS5zscWKookLVIKJorcZR9625GLKSJ3RN5xwAkcGvzS0KJ228BInUtZBIBrRyUsk0JIrqCKVP5nceWBVuPgItmxdORK7DphLYFJto/pzAg/do7Mudu4pRBBKRBwN/hr3Il6hqKzGdiDwa+2q3Ah9W1fRJMqKh5BeIfKyR08B5y6AmjzEVIlCYqiZOSDDjCNWFCqlJVFYG0QoRoRQQEWaAYpy9ZHUerRrJmBdFweWFalqHJ63TrZ2gTlvnVTfPc4mPhWqbryvdl6pt0a3vOpa7fwt8gXyZvHCqelZPvXcaaPd+8/SzCyvJyiMidwF+C/gh90LTbfK1zYHADgpVMrtN/RQ2pgrUNkuewpGlMBUiykSqyKUNtCRR5SSPJ1BhKsrKYCqDQZlSUKgVGyUGo0rp5ocaVa490nI2z5A7OpZCocUf1CugIrHkgVj6eEL5dQ1UuFBVy6hwnYTqgH/onGgYI4FqXz2AiHhffTjZ9SvAC1X1qwCq+qW5ehH+oyS2VdI5FgklUKDKWTe1JU9hKiuFRCmkYhLO9dBW4zwqlVoKzRxxrNOhoHkny1gPnXNte4+cSNseSqWIdxikkQZSpiTSWBJBPHhr0mjjPDBO3JigjO+OcSTRWJULVbXY+xZcxzyeuAUk0LGOMQTK+ervnZS5K4CIvA872i5Q1belFbm3A88HKE4/PTmY/JOc8yAtE0qfWqpkyDMR9+vWPYEmxjqjc6pcSJ6JVMyMQcoCI8q0bKwlK6UsiTQgkYpkNZ3QORDZOmVAnrIhT8upkHEghHaXJYSCkYg0XuKEjoj6Ve7wWECYRewfXBWLqHDHOlblRJgAdwHuh3VFvkdEvltVvxYWcnMAFwNs3fEOrRGc+8elkiaVPvXEJ23ybBYlE0lJVHXaQjWBpLAviDny7ITkQdDCDr5KpY5O0Fq/jC8okiCB48D439LuM6Umkkhjr1xYbWDreMkjApVjQtZkcWpdl5TUhEj+2Ggy6YkpgVaVled64FJVnarqZ4BPYAk1P7wt0VLj2tKncE6D2ong7J2QPJtFyaYp2SxmHCimHCimnFzscFIx5eTJkXo5ye0/UEw5OPG/Uw4UMw4UM7aKGRtFyWZR1hEM9VxTGJng++sRSp9QsjipQwVmppiZJZNdV8wUt2izb6ZIqZipIjPcOYopGzKaUiP3eKRC+lsbSp4cOhwXQ/84LeNl1FlLfGRYRH7L7f+4iPx4sP86l3DxQ+FErojcQ0Qu9/tF5AfmvcoUq8rK80as3/0vxX7A6K7Ap4cqHny6JU9LaEska+doTShPnk1TMjEVm2bGpinZMCUTKdmQqpZCRTBKSudA2HJfmJ5pwTQoy6zpQ1kZqsLODxkjLoOVEwOBKoTfHdg/Xl0zTn0zM3XEAamlkLakUHhPVAIpVAlqlGrS3FBxzgwVRyJJ1LWAVJpud/1PhjxxC0igZYJJXQLG84C7A7cD3ikid1VVb6zeX1W/nDT5HOD3VfWtYvO3PwerNS2MVWXleTtwjnsPvQSeqaqDs8XzoI5CoCFORCJT1Q4Dr7YdKKZMpOKkYocNqdgqZkykpKBiw7QTtZRqqNTmeD5STZiago1qgimd565U62iY2HkiSyITvVLRmgeC2HEQkaghj5k66RIQCP9r/xF474oa7HxUQT34DUI1sWUkJXTV+AI0kEBzC5khzG88jXFQnYudHgEbTPoCsaEE5wKvdbngPuPmiX4AO+Ha2UOWSUGdwaqy8ig24cPT5+5BRsrUhxL3Nc0DviV9vIPAk2fTzGrybJkZJxVTNqTkgLG/RiqKZAhNtaBEmFYTtsyMI9WEI1Jazx0KbNau7rIwTCtL3MoYK4USN7aXPC3XdS19GvLU6lkJMlNHInUu7aafKoIWVuqYUtAJaCFUqCURjUewJV1c7+rq+tS48P8wxhZScmrbUCzcGAdVFEwqIl/HBpPeHrg8OdcHoirwd2KfuH8etPlU4O0i8sdY8+XfD1zVIPY+EsEjsXPS/2w7t3Vz3LurQxJtmrImzkGzw8Fihy0z5YBYAm1ISRFY56Ua6zhQK3msFNrgcLlZOx88Zmoo1TArDaUxzLxLO/cwCMmjIXm0IY8n0EytDTNTpKwaSRQSqBCYWRJVE+PeW7I3sJKmLR/iRIZMTSyPu59djoSMY6EXbRVu1xMrduCHVfVzbj7yHSLyMVV9D/CrwNNU9Q1u4v+lwAOXaWj/EKgPAVnCcJxWlEEtiazN46XPKZNtDpodTjZH2JIpB8yUAjfh6mKsKwylClOdsKMFR3SDbdmwNpBT4SpHMv/G5bQoKKrKDv6UPMGTvyZRGTgKyoY8xVSRaeXUOFdfWUGpSFU1BHI5h3VirPetUnTDqp0ixh12k7tGIhunXvd9C38zt1vnIQ7NNc6JZYJJO89VVf/7JRH5a6xq9x5sTN1TXPm/YgVf6j42COQgEZGoJzrTSdKJVEzESqCDZoeDZofTisMcNEc41dxsJRBWAnk1rsSSp0TYVvt5jptkK5JW3ung54qmVWFd3aayr4FjmjAaYvKkDoTIszatKHYqZFZhphUyLZFZBZ48To1TESgEmRp0YpCqsBHjmwadKWK8BJImgFRdUpSgT60spKv57+Qk0BCWSax4KfBqEfkTrBPhLsD/EZGTAaOqN7r1c4Bnu7o+D/wo8G7gx4Br5+1wiv1PoIwXLoWf0/HkMaJsSMWGlBwsjnALs82tikOcara5hdnmgJRsiEY2UIl9nWFbC7Z1wrZscMBMOVBNKcqAaCpM1bAzmdgEgmbCtKgCNU7jGX3vfavVN6e67Tjpc6TC7FSYnRKzUyLTEqYlMistgcrKTuCos22KAoyBrQmUaudLjWDchG6kLnpyJHYQJMRZBYMUGOm6rk9ZIpjUlXs91uEwA56kqqWI3Bb4a/fKwgR4dTCp/yvYpIsT7BupabL7ubH/CTQScZxbVe+z0ucmbl0c4puKbU4WOCgFBqEQwWCoqChVqVC2tWRbZxzWHQ7qjpVUVDVJSzXMqoJZVbBtJmxOZmzPJnEG1ABZx4Gf29lx5DlSYnZmyJEpslPCdAazGZSlWwLdqDAwmSCzTWSrsgHiRtCJIKXCRGJCZMgh+d1LYwEVbtlg0j8E/jDZ92ngezvK/yNwr/l72Y1jhkC9+a0zME4CbcrMqWFVTZ4tmbAhhYurtqikYlpPIZSWVMw4IFOmpuCATjlc2bompqwdFl5t9JO7NYKnvvfCNW5rbZwH08qSZ2eGHJnBzo4l0M4UdSTSsoRKwQhSFI5AM+AgslHYpdR4AtQnruu7pyt2Zc8/+Xrs45ghkM3X1v0fqpLBUqlhqgXb1QaHqy0OV5t8TabsyMypcIYiEBclXgJVbKtwuJpwk25wWDe5sTyJw9Umh6stjlQTbi43mKlL3+TmjnwUd6MjEU16amHnjaqJUm4KUhmnainoBOMcBDIpLIE2N2o1TspgzmoycSrcJrq5gW5O0MJQTUzdnm1fWtKwvpf+d5WRNwoyO/FCeY4ZAnmothWQyo2U5pUE45wCBdu6wY3VATbKGTtacKrZ5oDM2JCKIvhQWYk4B8IW29WEbd3gsG5xU7XFjeUBbqwO8PXZSXxjdhI3l5scnm2yXU6YlQVlaewsfIbfahoCWRKBTMRGDmwJ9St9IkghyMTYZTaBWYl7U82WMWLJUxh0o0A3J1SbE6rNwpKnkLYaeTTH9AIq3LGO/U+gHrWgCkZHTZ46otp+ZvBQeaB2AEzVEmNTZhwQG5sTubFr0m1aL5yTXIcdiQ6VWxwqN9kuN9guJy59k31rVdWFsvgAUokXFUegiSCVUm7SlKVARTCFIMYgRYGUJVJOGk8cxG7sorAubL8U4sjqgl2NazN4zaIze9AKSGbV1NWxVURuhf1I1lnAdcCj/esySblfwH4lHuBCVX2F238v4OXASVgb6ynOe3cB1pngPwD3X1T1LSKygXVr3xPLi1eq6n8f6uf+J1AAdfEolUr9tXVVaZPHxbMdKScYNtkQqwLtOAJZm2iWSCBDpYYdLWqiHak2LIHKTQ6VW9xUbrFdbnB4tsH2bMLObGInUyuDVjY4LTKDJJBAgRRqiNTYDSoGLWxsnUwEmRmr3pVaRyKo9d2jxqATsarbhj1PN3yEgicOgT897lNre0VYxInQg2cBl6nqRS7I9FnAb0btWZL9HnA29lF7pYul+yrwIixR3o8l0IOBt7pTn6eqf5y091PAlqp+t4gcBD4qIq9R1ev6OnlsEMi9f2PDUKwK16Sgkmh9VhlmUjCTiiMywUjFjeUBplpw0EUXbNUTqfF/vHR2kyXQhCPuc+03lxvcVG5xc7nRkj5VZdU3DSKfa/gnf0KiqrBSqJoE8zN1WYMpBZmoi4Vr4uFsAKlYKeM8b1UR/NbSR+I2Q0INSZ6usoP/I+tlXCHOpQn0fAV27uY3kzI/DrxDVW8AEJF3AA92CUNOVdXL3f5XYr+V+la6ocDJzsV9ErADfGOok/uDQHWYiSOK4kaUtsv5sBT8i23Nm6SzymBEmalhp5oEmXesijOtCrbMjMPVZh0P5+HtptD5cKSacKSasOMcB97umZaFtX0qG87jCdypbgpWcoQOhcpfiB2lZmqvWUq1k6KFNk6GgOfezqlJ4slTWLvK7ydR3yJCZKRQNgyJjn1dl7laCXRbVf2CW/9/wG0zZboSM97eraf7PZ4sIo8FrgB+3UmsS7Ck/QJwEBvyc8NQJ/eeQD0hI/Wg1GYuXd1+dWTzaaiqSAJZ93RIkBJhQ0qmal+SK6iiV7srFafGWTtoVrmIbDVO6mywUxZszwLbp3QeuEpq+yfUiTTYVatwWs+LOhsuuAGCjSYwLppAARdGVEcQCPbN07pOK3mqCY0KFyxI0o+gb2OlzDwv1SXoDSYVkXdi87Ol+O2oWmu7rMpJ/iJsel91v/8D+CVsuE+JjWo4HXiviLzTR4p3YWVZeVy5R2KZ/P3L5ttKoU7a2PFm7SAJVLfSSR8q2BHr2dqRiZU+KkzUqnZTLSjQaLLVu8C9FJpVhX0nqCrYqWzEwXY5Yad0eZ/Lwtk+4tQ3OzLrD9rVN4Ra+kQDuXDTOuqJZAurUcxMECON6lbE8zV24Ae2jlfbHHkq70AwXlLRkkIL2T1D5+RVuN5gUlXtDOQUkS+KyLeo6hdE5FuAXJ6NzxG/z3MmVtX7HPGX6MI4uS8GbfwF8Ldu82eAt6nqFPiSS09wNgPvtQ2+kRq89PQQ4G7AY9zLTGm5ubPe9z5TNF2CJ7sjU00ebVJRqZdCatgpCzf4rQS5udzgptkWN1feptm02+Wm3a7ccVf+cGlVNk+eI85tPXO2T1mZRvoghDZQ6vHy7/Dk7CHrVLAkqDaEasM6GcoNodowlBtSL9UkWFxZLYSqiMmTs38U2kTKEWsBgglB1IWPCF8OPgYO9/s3mTL+PbTTxX5Q+Bzg7U71+4aI3Me9O/RYf74jo8cjgI+49c9i4+NwMXT3AT421MlVZeWBJuv9M0fU2Q1Plgy7/MtgKp48xNlDg+fBDlAZp9oZZyO5qIGdqshm5vFzSH5ydKcsaiJOK2v3TJ3dMyu988C0SZ5zJIQkcpKlwkqhCsGgwesI4tJQ+TdSkxFdD/6EKAF5QgcGGXUu6leGMHNLKV25DXQR8HoReTzwL7gvbYv9ps8TVPWXVfUGEfkDbFAqwLMDu+WJNG7st9I4EJ4jIvewPeY64D+5/S/EvlF9Dfbu/KWqXj3UyZVk5ZEg672IdBKoNyuPios4js4I3mNR52Sw23aXOskjdoKxshHRIkpVJeSpTJ3eCvoIZMv6c6aORCF5yiqwfbSx1dqvXweZetxgFid50GAdqERsvJw0UdX1PFHaVc9TE0gXTx7vsEhJ48/rUuMy9tC8JFolgdwbzQ/I7L8C+OVg+2XAyzrKfVdm/893tHeIjpi7PiztRJA5st4PZeVpytHOUZh46SIp5H9tahqgoDCVHUilJYX3zvVlJ/XzSJ5EpTaf6mgkTxN5oJXJOg/iG0SUclcNNjrbJLyoSeHmh0Sjd3jS+aWYRBKQiNZ6S9qEqlpCnDaxRtruq3djHxMYQ6Chl56Wz3ofkiVU3zQ+rjWJqKVQVUHz3QRj80W5fZXaDKNVJa13h1J4R0LppI93TJSVfdGudHZPLdmqUDoGpE8Geih9qByRTEOMrnvhiQS4d3ri+xSm0aqljsQEqomU2DrhebXWGRKnS/IMSKQVq3DHBJbOyrPSrPcdRKpVJPFqnH06W9vBqnI2rax/qU1tEKb4Ly3Eb65CnFgxjKULJ2fLSpwUEiongSqVSPrUMXA5W8Vfht/tiRMSTalNN3U2EO63HtzajnCwN5vIGRCqczV5oEWklgs7R54FbKB1LFwGI196Whw5h4EnkkI9oRpIHS+FBGszVZXBuFzWULknsrHHnWqnjkB1gt5kDshea/Nlhgrv1cN53KQmTyR9+izx4BIjCWCICeSEZ00ItZctKdHC+lzTYf0tAkEjmXJkC+oI9+WODUGwwv9Ew0qy8iT777dQTxxpGvum4Y497KQQtkBVWb+Bd2WFJLLkqBrioJQqtdQRaTJd2z7b33BS1jsHykBt8+ShlkLSVt9CInXZGwGJ6m9S+cFnknsQPjhSJATKSqOAUNn5oH7uz4cVEmgFwaR/iHVfn66qpwTlnwfc320eBG6jqqcFx0/FepjfqKpPHurn4DzQUUX6CGw9at1PYAvVUQlQ2yhe5fLr3nNW+mw6yefYS7X5DUovaUKVLZU8lWni3sK5nw5PWboezs1Ec0MFrXkiwgDUzBxS9rwe8kT7B1S4rG3U+79j1fNAPpj0LsBlbjtCEEx6b+x0y++5+SCAN7l9cTdVn6aq91DVewD/k+AL3g5/gE1AMgr7g0C5x2uqugRk8U/9ekLVRwTQlCkrcbFqRU0mHz2QLp4oZWlq4szceqy2BaStbZ9Q98lch1cVw0GbesqieZyYYJoQjJB0OeJkQnlykinqEz2EmgMrJtC5NF/nfgU2GDRFHUzqpNM7sFHXqOrlQSxdF6Kv2LlXIG4L/N3YTu59LFyK0NMUqnH+A6Bq1Tmvyml4ks/ACYhzHohY9U2CxHI5+8e3F0o0T9aaPGHoDg2R7cmZUZcMztqccwO9TrtrnMkXSLIoIUhPE+mgzzkKWrZRD1lib9xIF7brZ8aNvcxHhpcJJh2EiHwr8G3A37ttg42L+znmyBW3twSKvG7NujiSRN/jAGdgW/IoYsP9TUwif8xXJ25EhLGIFZZcGpFH6loaKRdLnNjjBn6GMhxnLX+IJ04wmIVm3RNJPaGg/n5PaBdG9wyigR9qulkJkyFTVvp0qJ1RGx0QslKnNxZuj4JJPc4DLglyaT8ReIuqXh8+bIewvyRQSCi3XTviAMegWgJ5EtVvaorPl01NEA3Wg1pqKRM1HxAHaBwGAWla0kcTUg0hkELiw3oCEnlChtKIHIF8XX53D3H6pFKufy2yjBlPik1FPAd2MZh0DM4DnhRs/yBwXxF5InAK9rurh1S1ZXuF2HMCpYn+vMQRgv2RFGqTyP56+QOgbv6HNiklGeehiuTJgidT+Kq2J5YnTFv69F2kjxbw11t3zfcnQ57aA5eR0uH1RNcWqmYpeTLn5STSosjk618GPpj0IvqDSf9b4Dg4B/up0V6IyHdiX1moE9Gr6s8Gxx8HnD1EHtg3ToSO/1zgMIie9AQDOvqlmeTUZr6mcQDEdoy6/VXVLLZM85ZpRJ5QEmnQ7/pY5jpStSiRFL1SInUGmMwSlpOgXEedo4kiHbm+u7B6L9xFwINE5FqsTXIRgIicLSIvAXCBoz6Y9AMEwaQi8hwRuR44KCLXu1wIHudhv+ywtFq45xKoRnQpMuLRnpdEdVVOAvk8CoO3yttAidcvJk/y6xvrULE0vIxQdQt/Q9unKdbcguC8vJevvT7KJoK29OlT7QYgrIQ0NVYQTPobwG901H3BQNsvx0ZyD2L/ECiF+rFj/3oVrVHl3J+QRE43iogENZmyzWQGf2oLxXaO1NInetWgh6ARaeq+O8Ik+6Nk8IEK2gquzbRRI6PK1WUy66vCvDbQ8YD9S6AQ2kEiiGwiICISSu19G5JAmhChyXGQIw9R2VEIjL2aML6OQKBJsDva7pJAYSES4vjuJ/uyZFtC+vgOy2y42PGG/UegepAE8iccIBqQqH5CS82QkEi2vAwPisSREPUhRx4NgjsDQmW1znDgS6ymhapcJJWC8zQ5P+cmb7UX/LYIFZbrsofmtX9cx09ECbQvnAi1KpS7/9o+nncf+xFBrX6FjoOorC8TOBKickldWfKkToM+y7xzoDbHUrslNfxzjoWuGLdsHeT39fZv6FhadIVOBBG5lYi8Q0Sudb+nd5T7BVfmWhcXh4gcFJE3i8jHROQaEbkoKP+tInKZiFwtIu8WkTODY28Tka+JyN/m2sph7wnU9dDKebVyJKrL+rrCEdTsb3nfcpOiLeLQkKev/2MfvKJ5lSknLQIS9BGlRa4OkgzZRM2xxaSIqLrk+c2yJJaNhftjVf1O4PuAHxKRh/j92Kyj34P9blCYffS5QPaN1S6MIpAMf4r86SLyUcfqy1yYxGLIDUgN/q9ZSZSXMKNGXZY0xOTpkj7LjpGMpEntkxaZepboksLb10We8PiQMB0hhfZLLJyqHlbVdwGo6g5wFU2WnrvhwneAd7l2cGUvA26cp5ODBBqZleeD2Imn78GmtXrOPJ2o0ZI4wb4ciRTqb4T2EWnuRWLyuDbnJk8ysNU/4ROi5FzNWdtmaEnQRayWm7tuI7iokWpb0xg2h124uFi4YJnng1YriYUTkdOAh2GlGMCHgZ90648AbiEit56jXxFWkpXHs93hcmxA3lxoRyTQ/BNrSzsoV/uGbbk6dwBpJSMRNB5/xVrqdnOSsSWZPIz2f/JQgkt015SWDqsbnZoglWDJepY8XdJnThJlnAh7Ggvn0vS+Bnh+kCDxGcALXLTBe7DhQAvHUKwkK0+Cx9Ofg3gYEXnceuCOajbDQsREio6PaDL894SSjoQ8I6RP37/aEqV5WtTE8SSqr5foHsxxKYPEyW0vJX3c6fPaPUchFu5i4FpV/dOgzc/jJJCInAI8UlW/NlfHA6zUiSAiP4fN5vjcjuPne3FeHrqpORANXmnvTwa0V+cksz9V7cYuORUOesiTk0hDSNWlcDAnts485tuQWjeWPDmSzkNcKTValoSPhYM5EysCiMiF2C96PzXqo8gZ7tUFsHFzrSiGeTCGQGM+RY6IPBAreh+uqkdyFanqxap6tqqeXZxycn+rLfL4X4m2ayJF1jQj7JxkAVLitMiT9s3tHxdQGp3S2p8b5C3jv8PW6SJNq45MO539nNsGUhsZHy7LYeFYOOea/m2szX6ViHxIRHz4z/2Aj4vIJ7B2Vf2NVRF5L/BXwANc/NyPD3Vy6aw8ruHvA/4c6wHJidpR6LJv2r/Bcey+OhNWl+o2OAOZKdZB2tZ6bjtpKo2JqxvzOlyqwtHsa/WvY3BnHQ+546kkzJ2bg+m7SPvt11VhmVg4Vb2ejrukqpdgHV25Y/edt5+rysrzXOw7FH/lXkb6rKo+fN7O9HeENomgk0j1aZKuxGhJjzHqZE4qjUFAnoYsbRLV/QqrDx8QY9oJT81IvKaRnEpHfv8QTsBIhJVk5ekzBkfDDZhOKRSU6fLQxW4r9zP2f9oi0oDUCW2kPoQSJ9mXJRGQSp6hQNIceqVR0Om5SdIDqU68vFb7LxYuxBCJ8OvpaKM9aMe2F20nv4vWS6KeBb8RieqCcVtZNW6g603hdDtDnr7yIyG6EsfBMYe9D+XJQDJP+m5JkG5Le8mhr1ykpmXaTaXPouMmVZV8pUEmn5zjYPASc+dFbssetS53fAwUmFXxsgSWiYVLjl8qIh8JtrP1isgznbPhQyLyEREpXahQL/YXgYKBmFWN0oGdI1K6H4ZHXNYj17W+uK6TN+LjLjb74wE/r+s6qiMhztBka4Q5LleqKlqWxLKxcIjITwKHxtSrqs8N8sX9FvAPYz7xuL8IBPmneU5ChOtZ0syx5M7LttX0YS7pMzQww0vrIlJKqGw73WWzxOnpw9yw787Hy3I4lyXywrlJ0qcDFy5Qb5Qvrg/72gbKJhxJbZyMpyqupKeBvvHY5VRgDsdErtrgElq2XNBuzqybt/HxdtGIc0ZAllTbEiwbC+e/f3p4nnrFfuL+wcBgWl/YrwQKnAaSDihNDwTn5f758wz2LukXNrFKhwK0vXQJkTJdWAwj1LOheaReqELZItCefGRY7Bfovl1VnyYiZ3V3OVvvw4D3jVHfYL8SCCISQYc08gfCc1bWfnv0jCFPVCYMKE2I0iJRWueQZB2DsRKIESTtm0T1aKtte/WR4R8EzhaR67Bj/DYi8m734YOhes9jpPoGe2wDDT5TkuNZM2CMx20sOupptbuiAd3pPcudt8jSVU+ClbizVaEs42U5LBwLp6ovUtXbqepZwA8Dnwi+GtJZr4jcEvjRjray2H9OhBSZwdprT+c8bmOXTBu9kQqLYIhEvkwfEeZpq6eezmfOQk4EYFbGy3JYKi/cvPU6PAL4O1W9KXtmBvtXhQuRcxjQHtzLCqBeibhK9TCjzvX2YRU2UKatLBZuS1fheWtqWzIvXHD8OoKPDXfV6469nJH54Dz2nEAt26YPAw6DZbxjo9ocQLb9rhfrUueBb27gITEvRt/bvnJj7B9lFWrbMYc9J9DCSP+nq3pK7wYJ+0g00OZKPHBdGFP3GPIAqKKzEy8x3L4g0FxSqAu7MfDnwFKSYln3+6JtDGEseTzabuzjHqvKyrMlIq9zx9/f53vvbGOPCbAMRvXd6HwDclHP21iP3Cr7CrvhhTsmsKqsPI8HvqqqdwaeB/zRIp0ZE62y3zB3f/3gXGSQ7hZW0idFyzJalsGywaRikyZ+PAgQvY3bn03BJjbhon979RoRecKYfo6RQHVWHpdjy2flCRHGF12CfSV2YaUsDenaL6TalX6lg3cvlhVAFXQ6i5YlsXQwKfCzPkA0eFP6g+RTsH0B+EEXTHpv4FkicruhTo4h0JjvUNZlVHUGfB1YONdWDrnBe7SXNXqgq5VALBlM2t1NfZeq+vi4y3EJF1V1J8jlscVI8+aoOhFcYj2fXO/IZ57yjI/0lT+GcQbw5b3uxC7iO9IdN/LVt7+zev0Zye4DfbFwA1hFYsW/FJESeANwYeaDWlEKNhG5A/Bm4M7AM10KrF6MIdCYrDy+zPUumd0tga+kFbmbd7Hr7BV9cVLHMo7nawN7fek+Ve198nfUs5uJFX9WVT8nIrfAEujngVcGbfsUbD8atPOvwPc41e2NInKJqn6xr5ExYqrOyiMim9hgu0uTMmF80aOAv1/F5/PWOL6hqg9U1e/KLH+DC/oEGAgmzT7cVdX/3gi8Gmsj4errTcHmJM9HgMEsPYMEcjaNz8rzz8DrfVYeEfGZd14K3FpEPol9iWnw46xrrDGAhYNJRWQiImcAiMgG8B+xhAhTsD08TMEmImeKyElu/XRsEOrHB3upqnuyAOfvVdvra9v/14d1Ql0GXAu8E7iV23828JKg3C8Bn3TLL7p9JwNXAlcD1wB/BhTu2DuBLwIfcsulbv+DXPkPu99R1yju5DXWWGMB7P/XGdZYYx9j1wl0NMKA9gojru1xIvJvwWz4L+fq2Y8QkZeJyJckSAmVHBcReb679qtF5J5Hu4/7ArusxxbAp4A7AZtY/fJuSZknAi926+cBr9trHX+F1/Y44AV73dcFr+9HgHsCH+k4/lDsHIoA9wHev9d93otltyXQUQ8DOooYc23HLFT1PUDf253nYr81qqp6OXCadzufSNhtAu2LMKBdwphrA3ikU3EucTPdxwvGXv9xjbUTYXfxJuAstYGL76CRtGscJ9htAs0TBuS/aZkNA9qHGLw2Vf2KNjPdLwHudZT6djQw6sNrxzt2m0DHcxjQ4LUlNsHDsZEcxwsuBR7rvHH3Ab6uTfDnCYNdjcbWcR/neinwKhcGdAN2IO57jLy2X3PhTjPstT1uzzo8J0TkNdikhWeIyPXY9242AFT1xdjvRT0UGwFwGPjFvenp3mIdibDGGktg7URYY40lsCbQGmssgTWB1lhjCawJtMYaS2BNoDXWWAJrAq2xxhJYE2iNNZbAmkBrrLEE/j/jvn2TuT6SrQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from dune.mmesh import skeleton\n", "\n", "sk = skeleton(iuh)\n", "b = avg(sk) * avg(v) * dS\n", "\n", "uh.interpolate(0)\n", "scheme = galerkin([a == b, DirichletBC(space, 0)])\n", "scheme.solve(target=uh)\n", "\n", "figure = plt.figure(figsize=(3,3))\n", "plot(uh, figure=figure, gridLines=None)" ] }, { "cell_type": "markdown", "id": "exposed-shadow", "metadata": {}, "source": [ "## Compute jump of gradient of traces\n", "\n", "We can use trace and skeleton within UFL expressions. This is useful when using jump and average operators, but also to compute gradients and more. Together with the utility function `normals`, which returns the positive normal of the underlying bulk facet, we can determine the orientation of the restricted values." ] }, { "cell_type": "code", "execution_count": 7, "id": "analyzed-change", "metadata": {}, "outputs": [], "source": [ "from dune.mmesh import normals\n", "inormal = normals(igridView)\n", "jmp = jump(grad(trace(uh)), inormal)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.5" } }, "nbformat": 4, "nbformat_minor": 5 }