JuliaFEM.jl/docs/tutorials/2015-09-15-2d-segmentation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating mortar projection matrices\n",
    "\n",
    "Author(s): Jukka Aho\n",
    "\n",
    "**Abstract**: Evaluate mortar projection matrices in 2d and 3d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using JuliaFEM\n",
    "using JuliaFEM: Seg2, Basis, Field, FieldSet, dinterpolate, interpolate, get_connectivity\n",
    "using PyPlot\n",
    "using ForwardDiff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create some test boundaries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rlinspace (generic function with 1 method)"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function rlinspace(X1, X2, R, n, p0)\n",
    "    fx(ϕ, xc) = R*cos(ϕ) + xc\n",
    "    fy(ϕ, yc) = R*sin(ϕ) + yc\n",
    "    function F(p)\n",
    "        ϕ1, ϕ2, xc, yc = p\n",
    "        return [\n",
    "            fx(ϕ1, xc) - X1[1]\n",
    "            fy(ϕ1, yc) - X1[2]\n",
    "            fx(ϕ2, xc) - X2[1]\n",
    "            fy(ϕ2, yc) - X2[2]\n",
    "        ]\n",
    "    end\n",
    "    J = ForwardDiff.jacobian(F)\n",
    "    p = copy(p0)\n",
    "    for i=1:10\n",
    "        dp = -J(p) \\ F(p)\n",
    "        p += dp\n",
    "        if norm(dp) < 1.0e-9\n",
    "            println(\"Converged. p = $p\")\n",
    "            break\n",
    "        end\n",
    "    end\n",
    "    ϕ1, ϕ2, xc, yc = p\n",
    "    ϕ = linspace(ϕ2, ϕ1, n)\n",
    "    return fx(ϕ, xc), fy(ϕ, yc)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "interpolate (generic function with 11 methods)"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function FEM.interpolate(basis::Basis, field::Field, xis::Array{Vector,1})\n",
    "    [FEM.interpolate(basis, field, xi) for xi in xis]\n",
    "end\n",
    "#FEM.interpolate(Γ₁[1], \"geometry\", Vector[[-1.0], [1.0]], 0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converged. p = [4.2905787641613236,5.624156522861784,0.04706336239279979,6.561746550428803]"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzkAAAEkCAYAAAAMxf91AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlclXX+///HQRYRARXc0HADjFwiCDVsDHdJcAsXPlnqVPqZFiczTdvNqWb8VTbm1/rUlEU6TEmacrLcnawxLbDRTEMyJBVUVAhjO8D5/cF4RmI9bAcOz/vtxk17X9f7ul6XZ0Z88l4ug9lsNiMiIiIiImInHGxdgIiIiIiISH1SyBEREREREbuikCMiIiIiInZFIUdEREREROyKQo6IiIiIiNgVhRwREREREbErCjkiIiIiImJXFHJERERERMSuKOSIiIiIiIhdUcgRERERERG7UqeQ8/zzz+Pg4MCAAQNqdH5WVhZz586lY8eOtG3blhEjRnDo0KG6lCAiIiIiIlKGwWw2m2vT8fTp0/Tt2xcHBwd69erF4cOHqzy/pKSE3/3udxw+fJjFixfj5eXFmjVr+Pnnn0lMTMTPz69WDyAiIiIiInKtWoecGTNmcPHiRYqKisjMzOTIkSNVnv/hhx8yY8YM4uPjmTJlCgCZmZkEBAQQERHB+vXra1OGiIiIiIhIGbWarvb555/z0Ucf8eqrr2I2mzEYDNX2iY+Pp0uXLpaAA+Dt7c20adPYvHkzJpOpNqWIiIiIiIiUYXXIKS4u5qGHHuK+++6jX79+Ne536NAhgoODy7WHhoaSm5tLcnKytaWIiIiIiIiUY3XIeeONN0hLS2P58uVW9UtPT6dr167l2q+2nT171tpSREREREREynG05uSLFy/y9NNP8/TTT+Pl5WXVjfLz83FxcSnX3rp1awDy8vLKHcvMzGTbtm307NkTV1dXq+4nIiIiIiL2Iy8vj9TUVMaOHYu3t3eV51oVcp588km8vb156KGHrC7K1dWVgoKCcu35+fmW47+1bds2Zs6cafW9RERERETEPq1bt44777yzynNqHHJOnDjBW2+9xauvvsrp06ct7fn5+RQWFnLq1Ck8PDxo3759hf27du1a4ZS09PR0AHx8fMod69mzp+VBAgMDa1qq1JMFCxawcuVKW5ch9Uifqf3RZ2p/9JnaH32m9kefqW0cO3aMmTNnWjJCVWoccs6cOUNJSQnz589n/vz55Y736tWLhx9+mFdeeaXC/kFBQezbt6/cbmwHDhzAzc2NgICAcn2uju4EBgZWuGmBNCxPT0/9udsZfab2R5+p/dFnan/0mdoffaa2VZNlLDUOOQMGDGDTpk1lAorZbObJJ5/kypUr/PWvf6VPnz5A6ehMdnY2fn5+ODqW3iI6Opr4+Hg2btzIHXfcAZSuudmwYQNRUVE4OTlZ9XAiIiIiIiIVqXHI8fLyYuLEieXarw7VTZgwwdK2dOlSYmNjSU1NxdfXFygNOUOGDGHOnDl8//33eHl5sWbNGsxmM8uWLavrc4iIiIiIiABWbjxQEYPBUO5loBW1OTg4sHXrVhYtWsSqVavIy8tj0KBBxMbG4u/vX9cyREREREREgHoIOXv27CnXtnbtWtauXVuuvV27drz11lu89dZbdb2tNIKYmBhblyD1TJ+p/dFnan/0mdoffab2R59p02cwm81mWxdRmaSkJEJCQkhMTNTiLhERERGRFsyabODQSDWJiIiIiIg0CoUcERERERGxKwo5IiIiIiJiVxRyRERERETErijkiIiIiIiIXVHIERERERERu6KQIyIiIiIidkUhR0RERERE7IpCjoiIiIiI2BWFHBERERERsSsKOSIiIiIiYlcUckRERERExK4o5IiIiIiIiF1RyBEREREREbtiVcg5evQoU6dOpU+fPri5ueHl5UVYWBjr16+vtu+7776Lg4NDhV/nz5+v9QOIiIiIiIhcy9Gak9PS0rhy5QqzZ8/Gx8eH3Nxc4uPjueuuu0hNTeWJJ56o9hrLly+nV69eZdo8PT2tq1pERERERKQSVoWciIgIIiIiyrQ98MADhISE8Oabb9Yo5ERERBAcHGxdlSIiIiIiIjVU5zU5Dg4OdO/eHScnpxqdbzabycnJobi4uK63FhERERERKadWISc3N5fMzEx+/PFHVq5cybZt21i8eHGN+g4fPhxPT0/c3NyYOHEiKSkptSlBREREpMk6duwY/v7+uLi44OzsjIuLC/7+/hw7dszWpYm0CFZNV7vqkUce4c033yy9gKMjq1atYu7cuVX2cXNzY86cOQwfPhwPDw+++eYbXnnlFcLCwkhKSqJ79+61KUVERESkSTl69ChBQUEUFRWVaU9JSWHgwIEcPnyYwMBAG1Un0jIYzGaz2dpOP/zwA2fOnOHs2bOsX7+eHTt28PbbbzNr1iyrrvPll18ybNgw5s6dy+uvv17ueFJSEiEhISQmJmodj4iIiDQL/v7+Vc5U8fPz48SJE41YkYh9sCYb1Gokp2/fvvTt2xeAmTNnMnbsWB5++GGmTZuGq6trja8zdOhQBg8ezM6dO2tThoiIiEiTk5aWVqfjIlJ3tQo5v3XHHXewY8cOfvjhB4KCgqzq2717d5KTk6s8Z8GCBeW2mY6JiSEmJsbqWkVEREQaUnWTZGoxiUakxYmLiyMuLq5MW3Z2do3710vIycvLA0p3WrPWyZMn6dixY5XnrFy5UtPVREREpFkwGAx1Oi4iFQ9oXJ2uVhNWpZILFy6UazOZTMTGxuLl5UW/fv0ASE9P5/jx42UW3FXUd+vWrSQlJTFu3DhryhARERFpksxmcHX1rfIcHx+fRqpGpOWyaiRn7ty55OTkMGzYMHx8fMjIyGD9+vUkJyezdu1aWrVqBcDSpUuJjY0lNTUVX9/S/6OHhYURHBxMSEgInp6eJCUl8c477+Dr68vjjz9e/08mIiIi0ohKSuDBByE7ewsODgMpKSmq8Dyz2czZs2cVdkQakFUjOTNmzMDBwYHXX3+d+++/n1dffZXevXuzfft27r77bst5BoOh3FDsjBkzOHHiBC+++CLz589n+/btzJs3j6+//rra6WoiIiIiTVlJCfzhD/DGG/C3vwXy3XeH8fPzw9nZGScnJ5ydnfHz82Pbtm2UlJQwfPhw0tPTbV22iN2q1RbSjUVbSIuIiEhTV1IC990Ha9fCO+/A7NlVn5+SkkJ4eDju7u7s3buXzp07N0qdIs2dNdnA+p0CRERERASA4mL4/e/h3XchNrb6gAOl78nZs2cP2dnZjBgxgvPnzzd0mSItjkKOiIiISC0UF5eGmvffL/2aObPmff39/dmzZw8XL15k1KhRZGZmNlidIi2RQo6IiIiIlYqK4K67IC4O/v53+J//sf4affv2Zffu3Zw7d45Ro0Zx8eLF+i9UpIVSyBERERGxgskEd94JGzbAP/4B06fX/lo33HADu3bt4syZM4wePZrLly/XX6EiLZhCjoiIiEgNmUwQEwMbN8KHH0J0dN2v2b9/f3bt2kVaWhqjR48mKyur7hcVaeEUckRERERqoLAQpk2DLVvgo49g8uT6u/bAgQPZuXMnP/30E2PHjiU7O7v+Li7SAinkiIiIiFSjoKB01GbrVti0CSZMqP97BAUFsWPHDk6dOsWhQ4fq/wYiLYijrQsQERERacry80sDzs6d8PHHEBHRcPcKDg7mxx9/xM3NreFuItICKOSIiIiIVCI/v3Ra2t69pdPUxoxp+Hsq4IjUnUKOiIiISAXy8mDiRPjiC0hIgFGjbF2RiNSUQo6IiIjIb+TmQlQUfPUVfPIJDB9u64pExBoKOSIiIiLX+PVXiIyEr7+GTz+FYcNsXZGIWEshR0REROQ/rlyB22+HQ4fgs8/g1lttXZGI1IZCjoiIiAiQk1O6c9rhw7BtG4SF2boiEakthRwRERFp8bKzSwPO0aOwYwcMHmzrikSkLhRyREREpEXLyoKxYyE5ufRdOKGhtq5IROrKwZqTjx49ytSpU+nTpw9ubm54eXkRFhbG+vXra9Q/KyuLuXPn0rFjR9q2bcuIESP0Rl8RERGxmcuXYfRoOHECdu1SwBGxF1aN5KSlpXHlyhVmz56Nj48Pubm5xMfHc9ddd5GamsoTTzxRad+SkhLGjx/P4cOHWbx4MV5eXqxZs4bw8HASExPx8/Or88OIiIiI1NSlS6UB59Qp2L0bgoJsXZGI1BeD2Ww21+UCJSUlhISEcOnSJU6dOlXpeR9++CEzZswgPj6eKVOmAJCZmUlAQAAREREVjgYlJSUREhJCYmIiwcHBdSlTRERExCIzszTgnD5dOoIzcKCtKxKR6liTDayarlbhBRwc6N69O05OTlWeFx8fT5cuXSwBB8Db25tp06axefNmTCZTXUsRERERqdaFCzBiBJw5A3v2KOCI2KNahZzc3FwyMzP58ccfWblyJdu2bWPx4sVV9jl06FCFiSs0NJTc3FySk5NrU4qIiIhIjX333Tn69JlAevrP7N0L/fvbuiIRaQi1CjmPPPIInTp1wt/fn8cee4xVq1Yxd+7cKvukp6fTtWvXcu1X286ePVubUkRERERqJCMDJk8uIDf3CG3aDMfD47StSxKRBlKrkLNgwQJ27txJbGwsI0eO5MEHH+S9996rsk9+fj4uLi7l2lu3bg1AXl5ebUoRERERqdbZsxAeDrm5vmzfvgcwMWLECP2QVcRO1eo9OX379qVv374AzJw5k7Fjx/Lwww8zbdo0XF1dK+zj6upKQUFBufb8/HzL8cosWLAAT0/PMm0xMTHExMTUpnwRERFpQc6cgeHDIS8P/vlP8PPryZ49e7jtttsYMWIEe/bsqXC2iYjYTlxcHHFxcWXasrOza9y/Xl4Gescdd7Bjxw5++OEHgirZf7Fr164V/rQkPT0dAB8fn0qvv3LlSu2uJiIiIlb7+efSgGMylQac3r1L23v37s2ePXsIDw9n5MiR7Nmzh86dO9u2WBGxqGhA4+ruajVR593V4L9TzRwcKr9cUFAQSUlJ/HbH6gMHDuDm5kZAQEB9lCIiIiIClL7/5rbboLi4bMC5ys/Pj927d5OVlcXIkSM5f/68bQoVkXpnVci5cOFCuTaTyURsbCxeXl7069cPKB2dOX78OEVFRZbzoqOjOXfuHBs3brS0ZWZmsmHDBqKioqrdglpERESkpn76qTTgAOzdCz17VnxeQEAAu3fvJjMzk1GjRpGZmdlYJYpIA7JqutrcuXPJyclh2LBh+Pj4kJGRwfr160lOTmbt2rW0atUKgKVLlxIbG0tqaiq+vr5AacgZMmQIc+bM4fvvv8fLy4s1a9ZgNptZtmxZ/T+ZiIiItEgnT5ZOUXNyKn0PznXXVX3+9ddfz+7duxk+fDijRo1i165deHl5NU6xItIgrBrJmTFjBg4ODrz++uvcf//9vPrqq/Tu3Zvt27dz9913W84zGAwYDIayN3JwYOvWrUyfPp1Vq1axePFiOnXqxO7du/H396+fpxEREZEWLSWldATHxaV0ilp1AeeqG264gV27dnHmzBlGjx7N5cuXG7ZQEWlQBvNvF8k0IVcXFyUmJmrjAREREalScnLpCI67e+kITm02TDt8+DAjRoygZ8+e7Ny5k3bt2tV/oSJSK9Zkg3rZeEBERETElo4fLx3B8fQsXYNT2x2hBw4cyM6dO/npp58YO3asVVvWikjToZAjIiIizdr335e+6NPbuzTgdOlSt+sFBQWxY8cOkpOTGTduHL/88kt9lCkijUghR0RERJqt774rDTidO8Pu3dCpU/1cNzg4mB07dnDs2DEiIiLIycmpnwuLSKNQyBEREZFm6fDh0jU43brBrl3QsWP9Xv/mm29m+/btfPfddzz66KP1e3ERaVBWbSEtIiIi0hR8+y2MHFn6/psdO6BDh4a5z6BBg9i1axe9evVqmBuISINQyBEREZFmJSkJRo2CPn1g+3Zo375h73fzzTc37A1EpN5pupqIiIg0G19/XTqCExBQOoLT0AFHRJonhRwRERFpFg4cKB3BCQwsHcHRK2xEpDIKOSIiItLk/etfMHo0DBwI27aBh4etKxKRpkwhR0RERJq0L76AsWPhppvg00/B3d3WFYlIU6eQIyIiIk3W55/DuHEQGgpbt0LbtrauSESaA4UcERERaZL27IGICLjlFjAawc3N1hWJSHOhkCMiIiJNzq5dMH483HorbNkCbdrYuiIRaU4UckRERKRJ2b4dIiPhtttg82ZwdbV1RSLS3CjkiIiIiE0YjUacnJwwGAyWL0dHJ26/3cjIkfDxx9C6ta2rFJHmyKqQ8/XXX/Pggw/Sr18/2rZtS48ePZg+fTonTpyotu+7776Lg4NDhV/nz5+v9QOIiIhI87N582aioqIoKioq015cXERxcRT33GPExcVGxYlIs+dozcl/+ctf2L9/P1OnTmXgwIGkp6ezevVqgoOD+eqrr+jXr1+111i+fDm9evUq0+bp6Wld1SIiItKsRUdHV3l82rTJmEymRqpGROyNVSFn4cKFhIaG4uj4327Tp09nwIAB/PnPf+b999+v9hoREREEBwdbX6mIiIjYjd+O4Fh7XESkKlZNV7vlllvKBBwAPz8/brjhBo4fP16ja5jNZnJyciguLrbm1iIiIiIiIjVS540HzGYz586dw9vbu0bnDx8+HE9PT9zc3Jg4cSIpKSl1LUFERESakYu5F8Fg6ypExJ7VOeSsX7+es2fPMn369CrPc3NzY86cOaxZs4aPP/6YxYsXs2vXLsLCwjh9+nRdyxAREZEmrqCogFf2v4Lfa35QzXLc384cERGxRp3+Bjl+/DgPPPAAYWFhzJo1q8pzp06dytSpUy3/PWHCBMaOHcuwYcN4/vnnef311+tSioiIiDRRZrOZjcc28tjOx/gp6yfmhcwj7O0w7rrjrkr7vPfee41YoYjYm1qHnIyMDMaPH0/79u2Jj4/HYLB+3Hno0KEMHjyYnTt3VnneggULyu3AFhMTQ0xMjNX3FBERkcbz9ZmveWT7I3yR9gURfhFsnrGZfp1Kd2Ntl9COyZMnl9tkwMPDgxUrVjBu3Dg6dOhgi7JFxMbi4uKIi4sr05adnV3j/gaz2Wy29qbZ2dmEh4dz+vRp9u3bx/XXX2/tJSymTZvG7t27yczMLHcsKSmJkJAQEhMTtSObiIhIM/Jz9s88vvtx1h1eR/9O/Xl5zMuM6TOmyj6pqalcf/31zJkzh/j4eHx9fdm5cyft27dvpKpFpCmzJhtYvSYnPz+fqKgoUlJSMBqNdQo4ACdPnqRjx451uoaIiIg0DTkFOTy5+0kCVgew48cdvBn5Jt/O+7bagAPQs2dPFixYwHvvvcff//53Tp06xZgxY8jKymqEykXEnlgVcoqLi5k+fToHDhxgw4YNDB48uMLzMjIyOH78eJnh5wsXLpQ7b+vWrSQlJTFu3DgryxYREZGmpLikmLcS38L/NX9e3v8yC29ZyImHTnBfyH20cmhV4+ssXbqUtm3bsm7dOnbu3MmPP/7I2LFjrZqmIiJi9ctAExISiIqKIjMzk3Xr1pU5PnPmTACWLFlCbGwsqamp+Pr6AhAWFkZwcDAhISF4enqSlJTEO++8g6+vL48//ng9PY6IiIg0th0/7mDh9oUcOX+EOwfcyQsjX8DX07dW1/Lw8OC5557jD3/4A/Pnz2fnzp2MHDmSiIgItm3bhru7ez1XLyL2yKqQ8+9//xuDwUBCQgIJCQlljhkMBkvIMRgM5TYimDFjBp988gnbt28nNzcXHx8f5s2bxzPPPKPpaiIiIs3Q9xe+Z9GORWw9sZVbfW/l4L0HCe0WWufr3nvvvbz22mssXLiQPXv2sGPHDkaNGkVERASfffYZbdu2rYfqRcSe1WrjgcaijQdERESanvO/nufZvc/yZuKb9GjXgxWjVjAlcEqtdlqtzGeffUZERASbNm1i0qRJHDhwgDFjxhAUFMTWrVtxc3Ort3uJSPPQoBsPiIiISMuUX5TPii9X4P+aP38/8nf+MuovfH//99xxwx31GnAAxo0bx5gxY1i0aBGFhYUMHjyYzz77jKSkJKKiosjNza3X+4mIfVHIERERkSqZzWY++O4DAv9fII/vepxZN84iZX4KC8MW4uLo0mD3femllzh58qTlheG33HILn376KQcPHmTChAmYTKYGu7eING8KOSIiIlKpr05/xdB3hjLjoxkM6DSA7+7/jlURq/Bu493g9x4wYAD33nsvy5Yt49KlSwDceuutbN26leHDh+PoWOt3mouInVPIERERkXJSs1KJ+SiGW96+hbyiPHbdvYstMVu43rtu78ez1nPPPYfJZGL58uWWtmHDhvHEE0/U+xQ5EbEfCjkiIiJikZ2fzZKdS7h+9fX8M/WfvDPhHb657xtG9Bphk3o6d+7M0qVLWb16NSkpKTapQUSaH43zioiICEUlRfwt6W88vedprhReYcmtS3g07FHaOtt+u+YFCxbg7e1Njx49bF2KiDQTCjkiIiItmNls5rOUz3h0x6Mcu3CMu2+8m+dHPE83j262Ls3C1dWVuXPn2roMEWlGFHJERERaqI/2fsTKgyv5Mu9LbutxG+/PfZ/grnovnYg0fwo5IiIiLUzGlQye2v0Uf7vvb7i2c2XTx5uYeP1ELeQXEbuhjQdERERaiDxTHi/sewH/1/z56NhH3LfoPvKS83D80VEBR0TsikZyRERE7FyJuYR/fPcPluxcQsaVDB4c9CBPDnuS9q3bc/Kzkzz66KOMHTsWJycnW5cqIlIvNJIjIiJix75I+4IhfxvCnRvv5Gafmzl6/1FeGfsKHVw7YDAYePnll0lOTub//u//bF2qiEi9UcgRERGxQz9e+pGpG6byu7W/o8Rcwj9n/5ON0zfi7+Vf5rwbb7yROXPm8Oyzz5KVlWWjakVE6pdCjoiIiB3Jys/i0e2PEvj/Atn/835iJ8Vy8L6DDOsxrNI+y5cvJz8/n+eff74RKxURaTgKOSIiInbAVGxi9cHV+K3y441v3uDp254m+aFk7rrxLhwMVX+79/Hx4bHHHmPVqlWcPHmykSoWEWk4NQ45X3/9NQ8++CD9+vWjbdu29OjRg+nTp3PixIka9c/KymLu3Ll07NiRtm3bMmLECA4dOlTrwkVERKT0ZZ4JPyQw4PUBzP90PpOun8SJh07w5LAnaePUpsbXWbhwIR07dmTJkiUNWK2ISOOo8e5qf/nLX9i/fz9Tp05l4MCBpKens3r1aoKDg/nqq6/o169fpX1LSkoYP348hw8fZvHixXh5ebFmzRrCw8NJTEzEz8+vXh5GRESkJfk241sWbl/I7p92M7LXSD6I/oAbu9xYq2u1adOGF154gVmzZvHll18ydOjQeq5WRKTxGMxms7kmJ+7fv5/Q0FAcHf+bi1JSUhgwYADR0dG8//77lfb98MMPmTFjBvHx8UyZMgWAzMxMAgICiIiIYP369RX2S0pKIiQkhMTERIKD9QZmERERgLM5Z3lq91Os/XYtfb378tLol7jd//Y6v+umpKTE8r1+//79ODhoVruINB3WZIMa/+11yy23lAk4AH5+ftxwww0cP368yr7x8fF06dLFEnAAvL29mTZtGps3b8ZkMtW0DBERkRbr18Jfee6fz+H/mj+bf9jMaxGvcfh/DzM+YHy9vMzTwcGBV155hYMHD/KPf/yjHioWEbGNOv2Ixmw2c+7cOby9vas879ChQxWmrdDQUHJzc0lOTq5LGSJShWPHjuHv74+LiwvOzs64uLjg7+/PsWPHbF2aiNRQibmE2H/H0nd1X57f9zwPhD5AyvwUHhj0AE6t6vcFnrfddhuTJk1i6dKl5OXl1eu1RUQaS51Czvr16zl79izTp0+v8rz09HS6du1arv1q29mzZ+tShohU4ujRowwcOJCUlBQKCwsxmUwUFhaSkpLCwIEDFXREmoG9qXsJfSuUWR/PIuy6MI49cIwVo1fQrnW7BrvnihUrOHv2LK+++mqD3UNEpCHVOuQcP36cBx54gLCwMGbNmlXlufn5+bi4uJRrb926NYB+UiTSQCZNmkRRUVGFx4qKipgwYUIjVyQiNXXi4gkmfzCZ4e8Nx9HBkS/mfMGHUz+kd/veDX5vf39/HnzwQV544QXOnTvX4PcTEalvtQo5GRkZjB8/nvbt2xMfH1/tPGBXV1cKCgrKtefn51uOi0j9S0tLq9NxEWl8l/Iu8fBnD3PDmhtISk/i71P+zv579jPUt3F3O3vqqadwcnLi6aefbtT7iojUhxpvIX1VdnY2ERER/PLLL+zbt48uXbpU26dr164VTklLT08HSl9CVpUFCxbg6elZpi0mJoaYmBgrKhdpearbPLGksITkB5LxivSi3fB2tGrdqpEqE2mZjEYjkydPLjPC6ujoyKZNmxgTMYY1X6/huX8+R1FJEcuHL+ePg/+Iq5NtfhDYoUMHnn76aRYuXMhDDz1E//79bVKHiLRMcXFxxMXFlWnLzs6ucX+rQk5+fj5RUVGkpKSwc+dOrr/++hr1CwoKYt++fZjN5jKjPgcOHMDNzY2AgIAq+69cuVJbSIvUQnWjrAYHA5c+vcTZNWdxaONA+1Ht8Yr0wmu8Fy4+5aeYikjtbd68mUmTJpVrLyoqIioqii73deF89/PcF3wfy8KX0bltZxtUWdb9999PQUEBvr6+ti5FRFqYigY0rm4hXRM1DjnFxcVMnz6dAwcOsHnzZgYPHlzheRkZGWRlZeHn52fZcjo6Opr4+Hg2btzIHXfcAZS+J2fDhg1ERUXh5FS/O8OISClfX19SUlIqPd6jdw8GJw8m93guF40XuWi8SPIfkqEY2ga3LQ08kV64h7hjcKj79rQiLVl0dHSVxzPezuBI+hH6d2o6IybOzs489thjti5DRMRqNQ45CxcuJCEhgaioKDIzM1m3bl2Z4zNnzgRgyZIlxMbGkpqaavnJT3R0NEOGDGHOnDl8//33eHl5sWbNGsxmM8uWLavHxxGRa23ZsoUBAwZQXFxc7pijoyNbtmzBYDDgFuiGW6Abvot8MV0ycWnbJS4aL3LmtTOceu4UTp2d8BpfOsLTfnQA3fm0AAAgAElEQVR7HN2tnukq0uJVtgmIRQlNKuCIiDRnNf6Xyr///W8MBgMJCQkkJCSUOWYwGCwhx2AwlJsi4+DgwNatW1m0aBGrVq0iLy+PQYMGERsbi7+/fz08hohUJDAwkPj4eCZPnoyjo6Pl/5++vr5s2bKFwMDAcn2cOjjROaYznWM6U1JUwi/7f7GM8mS8k4HByUC78HaWUR7X3to4RKRKZjOkpmIAql4lJyIi9cVgrm5lsg1dnXeXmJioNTkitfTCCy/w4osvkpmZWeFW7tbIO5nHxU9KA0/W3izMhWbaBLaxBB6PMA8cHOv0+i2R5u8/oYa9e//7lZaGC1BYbdcm+y1ZRMTmrMkGmnMiYueMRiNjxoypc8ABcO3tSveHutP9oe4U5RRxeefl0hGe2Ax+/v9+xrGdIx0iOuAV6UWHcR1w6qD1dtICVBJqMBjgppsgOhrCw2HKFKhiytrVdawiIlJ3+htVxI5duHCBr776irfffrver+3o7kjHyR3pOLkj5hIzOYk5lmlt5+POgwN4DvW0jPK0CWxT7W5vIs1CTUPN734H7dpZun20aRNRUVGVXnbTpk0NXLiISMuhkCNixz799FPMZjO33357g97H4GDAI9QDj1APei3rRcGZAi5uLQ08qc+mcvKxk7Tu1dqyPXW78HY4uGhamzQTtQw1vxUZGUlCQkK59+QArFu3jsjIyAZ8CBGRlkUhR8SOGY1GBg0aROfOjfu+DZduLvjc54PPfT4U5xWTtTeLi8aLZG7O5MxrZ3Bwc6DD6A60HtUap1ud6HFjj0atT6RK9RRqKhIZGYnJZALg2Wef5a9//SslJSXs37+fO++8s54fRESk5dKPUkXsVGFhIdu2bbP5T4dbubbCK8KLgP8XwJDUIdx8+GZ6PNGDwvOFvP3g2/QO6s3um3aTuiyVnMQczCVaeC2NzGyGn36CtWth1izo2RN694Z77oHvvisNNVu2wMWLkJgIL78MUVFWB5zfCg8PJysri3vuuYc33niDY8eO1c/ziIiIRnJE7NUXX3zBL7/8YvOQcy2DwUDbAW1pO6AtPZb24IWJL3DjdzfSuW9nfl75M6nPpuLc1bn0nTyRXrQf1Z5Wbq1sXbbYm5qO1Nx6K7Rv32BlDBkyBBcXF7p164avry+LFi3CaDQ22P1ERFoShRwRO2U0GvHx8SEoKMjWpVSosLCQHXt38Oijj9LvqX6UmErI/jLbsnlB+t/SMbgYaD+8felubeM74NpT7+SRWmgioea3WrduzZAhQ/jiiy9YsWIFU6dOZceOHYwePbrRahARsVcKOSJ2ymg0EhkZ2WR3NPvtSJODkwPtw9vTPrw9fi/5kXsi1/JOnpSHUzA/aMatv9t/38kzxANDq6b5bGJj14Sat996i48OHWJrfr7NQ01FwsPDWbVqFfHx8QwdOpSFCxdy6NAhWrXSCKaISF0o5IjYoeTkZE6cOMHLL79s61IqVd1IUxv/NrR5uA3XPXwdRdlFXNpxyTLCk/bnNBy9HPGK+M+0trHtcWqnd/K0WFWM1Dj7+vJpfj6Z69bhffvtNg81vxUeHs6yZcv47rvveOWVVxg8eDDvvvsu99xzj61LExFp1hRyROyQ0WikdevWjBw50talVMqakSZHT0c6RXeiU3QnzMVmfvn6F8u0tnPrzkEraPe7dpZRHtcA1yY7giX1wIrpZ7fl5ECPHnzu6sqUJhZw4L/rcvbu3csf//hH/ud//ocnn3ySadOm4e7ubuvyRESaLe2uJmKH9u3bx4gRI2jTpo2tS6nQ1ZGm2myKYGhlwHOIJ73/1JvQb0MZkjYE/9X+OLg58NOTP3Hw+oMcDDhIyoIULu+6TElhSQM8gTSq2u5+1r49vr6+9O7dm71799r6KSp0dV3O1fpefPFFsrKyWLFihW0LExFp5jSSI2KHPvroIy5dumTrMipVnyNNra9rTbf/7Ua3/+1GcW4xl3df5qLxIuc3nOf0q6dp5d6KDmM7lG5eENEB507O9fAE0qDqeaOA8PDwJhty4L/rckpKSvD19eWRRx7hpZdeYu7cuVx33XW2Lk9EpFnSSI6IHXJwcMDb29vWZVTKaDQ2yEhTqzat8I70pu8bfbnl51sIORTCdYuvI//nfI7POc6/uvyLpFuSOPX8Ka78+wpmc9l38hiNRpycnDAYDJYvJycnbevb0OowUlMT4eHhHDlyhMzMzIZ9jloKDw/n8uXLHDlyBIAlS5bg4eHB448/buPKRESaL4UcEWlUWVlZ7Nu3r8Hf32MwGHAPcqfnkz0J+SqEsPQw+r7TF+duzqT9OY1vgr7hK9+v+OF/fyDTmMnHH35MVFQURUVFZa5TVFREVFSUgk4lahUMGzjU/NZtt90GwOeff16r/g3t2nU5AO7u7ixfvpx169bxzTff2LY4EZFmStPVRKRRbd++naKiIsaPH9+o93Xu7EzX2V3pOrsrJYUlZH2eVbp5QcJF0v8vnfu5v8r+kydPxmQyNVK1zcPmzZuZNGlSufarwTAhIaE0zNr4PTXXrsuZMmVKvV+/rq5dl/PHP/4RgN///vekpaXRrVs3G1cnItI8WR1yfv31V1asWMGBAwc4ePAgWVlZrF27llmzZlXb99133+X3v/99hccyMjLo1KmTteWISDNjNBoZOHAgvr6+NqvBwdmBDqM60GFUB/xW+pH7Qy7nAs9V2aeoqIjOX36Jk8GAs4ND2V8NBpwcHHCupO23fSpqq8l1SgoKOJ2cjL+/P54eHpVe2/E/IyoNLTo6usrjkydOxDRzZpN4+WZzWpfj4OCAo6Mjf/rTn2xdlohIs2V1yLlw4QLLly+nR48eBAUFsXfvXqu/mS5fvpxevXqVafP09LS2FBFpZoqLi9m6dSvz5s2zdSkWBoMBt+vdKKH6Xdj+2L07hSUlFJrNmMxmCktKyv56ze+vFBdjKioq03bteab/XOfatqLfrBEqJy2tdHrXq6/CjTdWearTNQGoqlBlbdC6tq2opOo/s6KSkv9OP7PxyzfDw8N55513yMzMbJLr1e677z5iYmK09bmISD2xOuT4+PhYRl0SExMJDQ21+qYREREEBwdb3U9EmrcDBw5w8eLFBl+P01Ae79GjQa9v/k94uhqIrg1DJrOZo23acAfwVr9+BAQFVRu0KgpVlV3b2nBmMptLR2aqk5jYoH9mNXXtupymOGVN09JEROqX1SHH2dnZMq3stzsT1ZTZbCYnJ4c2bdrQqlWrWl1DRJofo9GIt7c3gwYNsnUp5Tg6OpbbdOC3xxua4eroCuBWwd+Nv/6nhpvatyekXbsGr6c6huJiW5dQY019XY6IiNQvm+yuNnz4cDw9PXFzc2PixImkpKTYogwRaWRGo5Hbb7+9Sf5wY9OmTXU63hiubnzg7Nw03vVTXfBrjGBojaa+LkdEROpPo4YcNzc35syZw5o1a/j4449ZvHgxu3btIiwsjNOnTzdmKSLSyE6dOsWRI0ea7FS1yMhIpk6dWq7d0dHxv7uE2VhhYSEATk5ONq6kVHMIhtdq6u/LERGR+tOoIWfq1Km8/fbbzJw5kwkTJvDcc8+xbds2Ll68yPPPP9+YpYhII9u+fTuOjo6MGTPG1qVU6vjx45Ywc/DgwdI1MiZTkwg40PRGciIjI0lISCg3YtOqVasmEwyv1dTflyMiIvXH5i8DHTp0KIMHD2bnzp22LkVEGtC9997L0aNHm+xOildHmsLCwoCmEySudXUkpynVFhkZiclkYt26dQC4urry4osvNrmAA2XX5YiIiH1rEhOmu3fvTnJycqXHFyxYUO4fRjExMcTExDR0aSJSTwwGAwEBAbYuo1KffPIJjo6OhISEAE1nSti1ro7kNMXargawoUOHsnfvXhYtWmTjiiqmdTkiIs1DXFwccXFxZdqys7Nr3L9JhJyTJ0/SsWPHSo+vXLlSW06LSIMyGo0MGzbMMkrSlEZLrmqKIzlXmUwmDAYDw4cP589//jNFRUVNbuMBaPrvyxERkVIVDWgkJSVZfhhZnQabrpaRkcHx48fLbMl64cKFcudt3bqVpKQkxo0b11CliIhU6ddff2X37t2WqVfQNINEU9t44FqFhYU4OzsTHh5OTk4Ohw4dsnVJFRo+fDgzZswgJyfH1qWIiEgDqtWP2VavXk1WVhZnz54FYMuWLaSlpQEwf/58PDw8WLJkCbGxsaSmpuLr6wtAWFgYwcHBhISE4OnpSVJSEu+88w6+vr48/vjj9fRIIiLW2bVrFwUFBURGRlqmzjbFINGUA5jJZMLJyYmbb76ZNm3asHfv3lq9LLqhde/evdz0BxERsT+1Cjkvv/wyp06dAkrn2W/atImNGzdiMBi4++678fDwwGAwYPjN27BnzJjBJ598wvbt28nNzcXHx4d58+bxzDPPVDldTUSkIRmNRgICAvD39+fo0aNA0wwSzWEkx9nZucmvyxEREftXq5Dz008/VXvO2rVrWbt2bZm25cuXs3z58trcUkSkQZjNZoxGo2Xeb1MOEiaTCUdHx3I/QGoKro7kQOm6l6a8LkdEROyfzbeQFhGxpUOHDpGenm7Z8rgpTwm7OlrSFF1bW1NflyMiIvZPIUdEWjSj0YiHhwe33nor0PRHcppiXVA25Fy7LkdERMQWFHJEpEUzGo2MGzfOEh5MJhMODg60atXKxpWV15RHcq4NYNeuyxEREbEFhRwRabHOnTvH119/zfjx4y1thYWFTXa0xGQyNdmQ89sAFh4ezr59+8q8RkBERKSxKOSISIvVqVMnjhw5wsSJEy1tzs7O+Pj42LCqyjX1AHZtbVqXIyIitqSQIyItlsFgoH///nh6elra5s6dy8mTJ21YVeV8fX255ZZbbF1GhX47knN1Xc6//vUvG1YlIiItlfb2FBFpJubNm8e8efNsXUaFhg0bxvXXX2/5b2dnZ3744Qe6detmw6pERKSl0kiOiIjU2ezZs1m8eHGZtu7duzfJd/oYjUacnJwsL602GAw4OTlhNBptXZqIiNQThRwREWkxNm/eTFRUVLkNEYqKioiKilLQERGxEwo5IiLSYkRHR1d5fPLkyY1USe1pJEpEpHpakyMiInbPZIKff4aiouIqzysqKmL16tW0adMGNzc33NzcKv29q6tro79PafPmzUyaNKnCuqOiokhISCAyMrJRaxIRaYoUckREpNkrLITTpyE1teKvM2egpARKv+2ZqrzWwoULKSwsrNF9W7duXWUQqsvvKwpRNRmJMpmqfj4RkZZAIUdERJq8mocYMBjAxwd69iz9Gjbsv78fPbr6AFBQUEBRURG5ubn8+uuv5X6t7Pe/bbty5Qrnzp2r8NyCgoIaPfdvQ1R1L1fVy1dFREop5IiIiFWMRiOTJ08u8w9qR0dHNm3aVOupUvUVYnr2hOuuAxeXiu/j6OhYZRBwdHS0/Orh4YGHh0etnqc6RUVF5OXlVRmaKvr98ePHq712SsqjeHgMwt19EK1b92iSO9yJiDQ0hRwREamx2q4JubomxpoQ06MH9OplXYipzqZNm4iKiqryeGNwdHTE3d0dd3d3q/q98sor1Z6TmfkRp0+/DICTU0fc3Qfh4RGKu/sg3N1DcXb2rlXNIiLNiUKOiIjUWHVrQiZNmsybb5oabCSmriIjI0lISKj3kajGUpORqCFDfqKw8Dw5OV/zyy8Hyck5yOnTqygqugRA69a9LSM9Hh6DaNv2Jlq1atNYjyAi0iisCjm//vorK1as4MCBAxw8eJCsrCzWrl3LrFmzatQ/KyuLxYsXs2nTJvLy8hg0aBAvv/wyN910U62KFxGRxlXdmo/i4iLuvdd2IaYmIiMjm+3i/JqORDk7d8LLazxeXuMBMJvN5OeftISeX345SGbmx5SU5AOtaNt2gCX0uLsPws3tBgyGxt05TkSkPlkVci5cuMDy5cvp0aMHQUFB7N27t8ZzfUtKShg/fjyHDx9m8eLFeHl5sWbNGsLDw0lMTMTPz69WDyAiIg0v+WIyxmQjGABz1efm5dk2xNiz2o5EGQwGXF374Orah86dYwAoKTHx66/fWULPL7/sJz39b0AJDg5uuLuHWEKPu3uo1veISLNiVcjx8fEhIyODTp06kZiYSGhoaI37xsfHs3//fuLj45kyZQoA06ZNIyAggGeeeYb169dbV7mIiDQYU7GJL3/+koQfEjCeMJJ8MZnWjq1rFHIUcBpWfY1EOTg44e5+E+7uN+HjMw+AoqIrXLmS9J8RnwOcP7+Bn39+Cbh2fc8gyzofJyevOtchItIQrAo5zs7OdOrUCSgd+rZGfHw8Xbp0sQQcAG9vb6ZNm8a6deswmUw4OTlZdU0REak/l/Iu8emJTzGeMPJZymdk5WfRtW1XIgMieWn0S4zsPRLPZz0pKql+dzJpnhwd29Ku3TDatRtmaSssPMcvv3xtGfE5ffqv16zv6VPB+h7XKu9x7NgxJkyYQFpaGmazGYPBgK+vL1u2bCEwMLBBn09EWo5G+2506NAhgoODy7WHhoby5ptvkpycTL9+/RqrHBGRFs9sNnM88zgJyQkYk418+fOXlJhLCOkawsODHyYyIJLgrsFlpig1ld3JpPE4O3fG2zsSb+/SqXAVr+/ZdM36noFldnS7dn3P0aNHCQoKKre2KyUlhYEDB3L48GEFHRGpF40WctLT0wkPDy/X3rVrVwDOnj2rkCMi0sAKiwv5/NTnGJONJCQncPLySVwdXRndZzRvjH+D8QHj8XH3qbR/c9+dTOquZut7viQ9/S1+u74nMvL9SjevKCoqYsKECZw4caIRn0ZE7FWjhZz8/HxcKpio3bp1awDy8vIaqxQRkRblwq8X2HpiK8YTRralbCOnMIfuHt2J9I8kqm8Uw3sOx9Wp6ilG12rOu5NJw6h+fc9Bzp/fwOnT56q8TlpaWmOUKyItQKOFHFdXVwoKCsq15+fnW45XZsGCBXh6epZpi4mJISYmpn6LFBGxA2azme/Of2cZrfnq9FcADOo2iMVDFxMZEMmNnW/UTlnSoCpa32MwOAGVr+mydr2viNivuLg44uLiyrRlZ2fXuH+jhZyuXbty9uzZcu3p6elA6c5tlVm5cmWF63lERKRUQVEBe1P3WtbXnMo+hZuTG6P7jOZvE/7GeP/xdG7b2dZlSgtnMDhUc1zBW0RKVTSgkZSUREhISI36N1rICQoKYt++fZadVK46cOAAbm5uBAQENFYpIiJ24dyVc3xy4hOMyUa2/7idX02/4uvpS1RAFFEBUdzW87bSbZ9FmghfX19SUlKqPC4iUh8aJORkZGSQlZWFn5+fZTvR6Oho4uPj2bhxI3fccQcAmZmZbNiwgaioKG0fLSJSDbPZzL/P/dsyDe3gmYMYMHDLdbfwxO+eIDIgkv6d+uun4dJkbdmyhYEDB1a4+YCjoyNbtmyxQVUiYo+sDjmrV68mKyvLMvVsy5YtloWC8+fPx8PDgyVLlhAbG0tqaqrlpzLR0dEMGTKEOXPm8P333+Pl5cWaNWswm80sW7asHh9JRKTpMxqNNdqhLM+Ux57UPZaXcp7+5TTuzu6M9RvLA6EPEOEXQUe3jrZ4BBGrBQYGcvjwYb0nR0QanNUh5+WXX+bUqVNA6dzZTZs2sXHjRgwGA3fffTceHh4YDIZyP0l0cHBg69atLFq0iFWrVpGXl8egQYOIjY3F39+/fp5GRKQZ2Lx5M5MmTSrXXlRURFRUFLEfxpLfOx/jCSM7T+4k15RLr3a9mHL9FKL6RjGsxzCcWznboHKRugsMDNQ20SLS4AzmJryVydXFRYmJidp4QETshpOTU6XvCgHAAA7POjD0uqFEBkQSGRBJoHegpqGJiEiLZk02aLSNB0REpFSVAQfADOcfPY9XG6/GKUhERMTOVL2Xo4iI2IQCjoiISO0p5IiINJYzZ+CPf0SraURERBqWQo6ISENLS4P774feveH99/GoZm3N1a33RUREpHYUckREGsrJk3DffeDnBx9+CM8+C6mprK3mXSCbNm1qnPpERETslEKOiEh9S06G2bMhIAC2bIEXXoDUVFi6FDw8iIyMJCEhodyIjaOjIwkJCWXekyMiIiLW05wIEZH6cuwYPP88xMVB587w0kswdy60aVPu1MjISEwmkw2KFBERsX8KOSIidXXkCPzpT7BhA3TrBqtWwT33QOvWtq5MRESkRdJ0NRGR2jp0CKZMgYED4cABeOMNSEmBBx5QwBEREbEhhRwREWsdPAhRURAcDIcPw9tvw4kTpVPTXFxsXZ2IiEiLp5AjIlJT//oXRETA4MGlmwvExsLx4/D734OTk62rExERkf9QyBERqc4//wmjRsHQoaXvvImLg++/h7vuAr3TRkREpMlRyBERqYjZDLt2wW23QXg4XLhQurHAkSMwYwa0amXrCkVERKQSCjkiItcym+Gzz+DWW0tHb65cgY8/Lt1kIDoaHPTXpoiISFOn79YiIlAabozG0vU2ERFQXAyffALffAMTJyrciIiINCNWf9cuKCjgsccew8fHhzZt2jBkyBB27txZbb93330XBweHCr/Onz9fq+JFROqspAQ2bYKQkNId01xcYPt22L8fbr8dDAZbVygiIiJWsnrF7OzZs/noo49YsGAB/v7+rF27lttvv509e/YwdOjQavsvX76cXr16lWnz9PS0tgwRkbopLoaPPip9ieeRIzB8OOzZU7oGR8FGRESkWbMq5Bw8eJAPPviAl156iUceeQSAu+66i/79+7N48WK+/PLLaq8RERFBcHBw7aoVEamr4mL44IPScHPsGIwZA2vWlK7BEREREbtg1XS1+Ph4HB0dmTt3rqXNxcWFe+65h/3793PmzJlqr2E2m8nJyaG4uNj6akVEqnDs2DH8/f1xcXHB2dkZFxcX/P39OXbsGBQVwXvvQWAg3Hkn9OpVOiVt2zYFHBERETtjVcg5dOgQAQEBtG3btkx7aGgoAN9++2211xg+fDienp64ubkxceJEUlJSrClBRKRCR48eZeDAgaSkpFBYWIjJZKKwsJCUlBQG9u/PsV69YPZsuOGG0s0EPvkEhgyxddkiIiLSAKyarpaenk7Xrl3LtV9tO3v2bKV93dzcmDNnDsOHD8fDw4NvvvmGV155hbCwMJKSkujevbuVpYuI/NekSZMoKiqq8FhRSQkTLl/mxLffwo03NnJlIiIi0tisCjl5eXm4uLiUa2/durXleGWmTp3K1KlTLf89YcIExo4dy7Bhw3j++ed5/fXXrSlFRKSMtLS0qo+bTAo4IiIiLYRVIcfV1ZWCgoJy7fn5+Zbj1hg6dCiDBw+udgvqBQsWlNuBLSYmhpiYGKvuJyL2y2w21+m4iIiINB1xcXHExcWVacvOzq5xf6tCTteuXSuckpaeng6Aj4+PNZcDoHv37iQnJ1d5zsqVK7Ujm4hUyVDNts/VHRcREZGmo6IBjaSkJEJCQmrU36qNB2666SaSk5PJyckp037gwAEAgoKCrLkcACdPnqRjx45W9xMRuZavr2+djouIiIj9sCrkREdHU1xczJtvvmlpKygoYO3atQwZMoRu3boBkJGRwfHjx8ssAr5w4UK5623dupWkpCTGjRtX2/pFRADYsmULjo4VD047OjqyZcuWRq5IREREbMWq6WqDBg1i6tSpLF26lPPnz9OnTx/ee+890tLSWLt2reW8JUuWEBsbS2pqquWnp2FhYQQHBxMSEoKnpydJSUm88847+Pr68vjjj9fvU4lIixMYGMjhw4eZMGECaWlpmM1mDAYDvr6+bNmyhcDAQFuXKCIiIo3EqpADEBsby1NPPcX777/P5cuXufHGGzEajdx6zcv0DAZDufnvM2bM4JNPPmH79u3k5ubi4+PDvHnzeOaZZzRdTUTqRWBgICdOnLB1GSIiImJjBnMT3nLo6uKixMREbTwgIiIiItKCWZMNrFqTIyIiIiIi0tQp5IiIiIiIiF1RyBEREREREbuikCMiIiIiInZFIUdEREREROyKQo6IiIiIiNgVhRwREREREbErCjkiIiIiImJXFHJERERERMSuKOSIiIiIiIhdUcgRERERERG7opAjIiIiIiJ2RSFHRERERETsikKOiIiIiIjYFYUcERERERGxK1aFnIKCAh577DF8fHxo06YNQ4YMYefOnTXqm5WVxdy5c+nYsSNt27ZlxIgRHDp0qFZFi4iIiIiIVMaqkDN79mxWrlzJXXfdxapVq2jVqhW33347X375ZZX9SkpKGD9+PHFxccyfP58VK1Zw/vx5wsPDSUlJqdMDiIiIiIiIXMuxpicePHiQDz74gJdeeolHHnkEgLvuuov+/fuzePHiKoNOfHw8+/fvJz4+nilTpgAwbdo0AgICeOaZZ1i/fn0dH0NERERERKRUjUdy4uPjcXR0ZO7cuZY2FxcX7rnnHvbv38+ZM2eq7NulSxdLwAHw9vZm2rRpbN68GZPJVMvyRUREREREyqpxyDl06BABAQG0bdu2THtoaCgA3377bZV9g4ODy7WHhoaSm5tLcnJyTcsQERERERGpUo1DTnp6Ol27di3XfrXt7NmzDdJXRERERETEGjUOOXl5ebi4uJRrb926teV4ZfLz82vdV0RERERExBo1Djmurq4UFBSUa8/Pz7ccb4i+YjtxcXG2LkHqmT5T+6PP1P7oM7U/+kztjz7Tpq/GIadr164VTitLT08HwMfHp0H6AixYsIAJEyaU+dL/uBqe/oztjz5T+6PP1P7oM7U/+kztjz7ThhcXF1fu3/8LFiyocf8abyF90003sXfvXnJycnB3d7e0HzhwAICgoKBK+wYFBbFv3z7MZjMGg6FMXzc3NwICAqq898qVKyvcuEBEREREROxPTEwMMTExZdqSkpIICQmpUf8aj+RER0dTXFzMm2++aWkrKChg7dq1DBkyhG7dugGQkZHB8ePHKSoqKtP33LlzbCDlwnoAAAlzSURBVNy40dKWmZnJhg0biIqKwsnJqaZliIiIiIiIVKnGIzmDBg1i6tSpLF26lPPnz9OnTx/ee+890tLSWLt2reW8JUuWEBsbS2pqKr6+vkBpyBkyZAhz5szh+++/x8vLizVr1mA2m1m2bFn9P5WIiIiIiLRYNQ45ALGxsTz11FO8//77XL58mRtvvBGj0citt95qOcdgMJSZkgbg4ODA1q1bWbRoEatWrSIvL49BgwYRGxuLv79/pfe7uuvasWPHrClT6kl2djZJSUm2LkPqkT5T+6PP1P7oM7U/+kztjz5T27iaCWqyM7PBbDabG7qg2lq/fj0zZ860dRkiIiL/f3t3G1Nj2MAB/F/L0Jv1rr44R2/IidCWCjGiWZSVQ4yJ2kgpPrRVbCizMTvoQ8bEivlQjA+GnZSKIlPL1lKx3ot1Fr1HOc+HZ/VIPTmdytU5/r+tD133/eF/du+u/e/7OtdFREQzRGZmJvbs2TPuOTO65LS1teHp06eQSCRcZpqIiIiI6B/W29uL2tpabN68GdbW1uOeO6NLDhERERER0URpvLoaERERERGRLmDJISIiIiIivcKSQ0REREREeoUlh4iIiIiI9ApLDhERERER6RWWHJqwiIgIGBoaIjAwUHQU0lJOTg7Cw8Ph4uICExMTODo6IiIiAq2traKjkQb6+/sRHx8PBwcHGBsbw8vLC0qlUnQs0lJJSQmOHj0KNzc3mJqaYsGCBZDL5aiurhYdjaZQSkoKDA0NIZPJREehSXj37h22bdsGKysrmJiYQCaT4erVq6Jj0Ri4hDRNyNu3b+Ht7Q0jIyNs3LgRjx49Eh2JtLBq1Sp8/foVoaGhcHZ2xsePH5GamgpjY2OUlZXBzs5OdEQax+7du5GdnY24uDg4OzsjPT0dJSUlyM3NhY+Pj+h4NEEhISEoKipCaGgo3N3d0dLSgtTUVHR1daG4uBhubm6iI9IkNTY2wtXVFYaGhpBKpSgvLxcdibTw7NkzBAYGYuXKlZDL5TA1NUVNTQ3UajXOnz8vOh79hiWHNKZWq+Hj4wM3NzcolUrIZDKWHB1VWFgIX1/fEWMFBQVYt24dEhMTcfbsWUHJ6E/evHkDLy8vXLx4EcePHwfw3zc7S5cuha2tLV6+fCk4IU1UUVERPD09YWRkNDxWU1MDmUyGkJAQZGRkCExHU2HXrl1QqVQYGBhAW1sb3r9/LzoSTVBHRwdcXFzg6+uLrKws0XFIA5yuRhrLyMhARUUFkpOTwW6s234vOACwZs0aWFpaorKyUkAi0lRWVhaMjIwQGRk5PDZ79mwcPHgQRUVFaGpqEpiOtLF69eoRBQcAnJycsGTJEt6PeiA/Px/Z2dlQKBRQq9UwMDAQHYm0cPfuXXz58gUpKSkAgO7ubvz8+VNwKhoPSw5ppLOzE/Hx8UhISOBUJj3V1dWFzs5OWFtbi45C4ygtLYWLiwtMTU1HjHt6egIAysrKRMSiKaZWq/H582fejzpucHAQ0dHRiIiI4LRDHadUKmFubo6Ghga4urrCzMwM8+bNw5EjR9Df3y86Ho2BJYc0cubMGZiYmCAuLk50FJomCoUCP378gFwuFx2FxtHS0gJ7e/tR40Njzc3NfzsSTYM7d+6gubmZ96OOS0tLQ319PacA64Hq6moMDAwgKCgIAQEBuH//PsLDw5GWloYDBw6IjkdjMPrzKaRP1Gq1xk8c5syZAwCoqqrClStXcO/ePcyaNWs645EWtLmmv8vPz8fp06chl8vh5+c3heloqvX29mL27NmjxoeubW9v79+ORFOssrISUVFR8Pb2xv79+0XHIS2pVCqcOnUKp06dgpWVleg4NEldXV3o6enB4cOHoVAoAABBQUH4/v07rl27hjNnzsDJyUlwSvoV3+T8Y168eAFjY2ONfqqqqgAAx44dg4+PD4KDgwWnp7Foc01/VVlZieDgYLi7u+PGjRsCPgFNxNy5c8cstX19fcPHSXe1trZi69atsLCwQFZWFr+/ocOSkpJgbW2N6Oho0VFoCgz9bd29e/eI8aHfi4uL/3omGh/f5PxjFi9ejFu3bml07vz58/H8+XM8ffoU9+/fR21t7fCxgYEB9PT0oK6uDpaWljAzM5uewPRHE72mv2poaIC/vz8sLCzw+PFjmJiYTENCmkr29vZjTklraWkBADg4OPztSDRFvn37hoCAAHR0dKCgoGDU/Uq6o7q6GtevX4dCoUBjY+PweF9fH75//466ujqYm5vDwsJCYEqaCAcHB1RUVIz6XrKtrS0AoL29XUQsGgdLzj/Gzs4O+/bt0/j8+vp6AMCOHTtGHWtuboZUKoVCoUBMTMyUZaSJmeg1HaJSqeDv748fP34gNzeXC0roCA8PD+Tl5aGzs3PEw4XXr18DAJYvXy4qGk1CX18fAgMDUVNTA6VSiUWLFomORJPQ1NSEnz9/IiYmZsz/j1KpFLGxsbh06ZKAdKSNVatWQalUorGxEc7OzsPjQw+dbGxsREWj/4P75NC4GhoaUFpaOmJMrVYjMjISEokEiYmJWLp0KRYuXCgoIWmju7sbGzZswIcPH5CbmwsPDw/RkUhDQ/vkXLhwASdOnADwv31ybGxs8OrVK8EJaaIGBwexY8cOPHnyBA8fPsSWLVtER6JJUqlUKCwsHDHdUK1WIykpCV1dXbh8+TIcHR254poOKSsrw4oVKxAWFobMzMzh8bCwMGRnZ6Ouro5vX2cYlhzSikQigbu7OzcD1VFBQUF49OgRwsPDRy00YGZmhu3bt4sJRhqRy+V48OAB4uLi4OjoiNu3b+Pt27fIyckZcw8kmtliY2Nx5coVBAYGIjQ0dNTxvXv3CkhF08HPzw8qlYqbgeqoQ4cO4ebNm9i5cyfWrl2LvLw8ZGVlISEhAcnJyaLj0W9YckgrUqkUMpmMJUdHSaVS1NfXj7mpq0QiwadPnwSkIk319/fj5MmTyMzMRHt7O5YtW4azZ89i06ZNoqORFtavX4/8/Pwx70cDAwMMDg4KSEXTYf369VCpVCgvLxcdhbQwMDCAc+fOIT09Hc3NzZBIJIiKiuKU/RmKJYeIiIiIiPQKl5AmIiIiIiK9wpJDRERERER6hSWHiIiIiIj0CksOERERERHpFZYcIiIiIiLSKyw5RERERESkV1hyiIiIiIhIr7DkEBERERGRXmHJISIiIiIivcKSQ0REREREeoUlh4iIiIiI9ApLDhERERER6ZX/AEV9OZn+FOlUAAAAAElFTkSuQmCC",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x000000002F321F98>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(-0.1,4.0)"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Converged. p = [2.110067133637596,1.2802355090457196,1.5675520985912952,-3.7904167887303712]\n"
     ]
    }
   ],
   "source": [
    "srand(42)\n",
    "\n",
    "nsl = 6\n",
    "nm = 5\n",
    "\n",
    "ϕ11 = 5.4\n",
    "ϕ12 = 4.6\n",
    "R1 = 5\n",
    "xc1 = -0.8\n",
    "yc1 = 9.9\n",
    "X1 = [-2.0, 2.0]\n",
    "X2 = [4.0, 3.5]\n",
    "x1, y1 = rlinspace(X1, X2, R1, nsl, [ϕ11, ϕ12, xc1, yc1])\n",
    "\n",
    "ϕ11 = 2.2\n",
    "ϕ12 = 1.1\n",
    "R1 = 5\n",
    "xc1 = 1.9\n",
    "yc1 = -4.1\n",
    "X1 = [-1.0, 0.5]\n",
    "X2 = [3.0, 1.0]\n",
    "x2, y2 = rlinspace(X1, X2, R1, nsl, [ϕ11, ϕ12, xc1, yc1])\n",
    "\n",
    "#x1 = linspace(0, 3, nsl) + rand(nsl)*0.1\n",
    "#y1 = 0.5*rand(nsl)\n",
    "#x2 = linspace(0, 3, nm) + rand(nm)*0.1 + 0.4\n",
    "#y2 = 0.5*rand(nm) + 0.5\n",
    "\n",
    "\"\"\"\n",
    "Calculate element local normal field\n",
    "\"\"\"\n",
    "function calculate_normals!(element::Seg2, time, field_name=\"normal\")\n",
    "    normal_fieldset = FieldSet(field_name)\n",
    "    normal_field = Field(time, Vector[])\n",
    "    tangent(xi) = dinterpolate(element, \"geometry\", xi, time)\n",
    "    for xi in Vector[[-1.0], [1.0]]\n",
    "        t = tangent(xi)\n",
    "        n = [-t[2], t[1]]\n",
    "        push!(normal_field, n/norm(n))\n",
    "    end\n",
    "    push!(normal_fieldset, normal_field)\n",
    "    push!(element, normal_fieldset)\n",
    "end\n",
    "\n",
    "function create_elements(X, sid=0)\n",
    "    Γ = []\n",
    "    nnodes = size(X, 2)\n",
    "    nelements = nnodes-1\n",
    "    for i=1:nelements\n",
    "        con = sid+[i, i+1]\n",
    "        el = Seg2(con)\n",
    "        push!(el, FieldSet(\"geometry\", [Field(0.0, Vector[X[:, i], X[:, i+1]])]))\n",
    "        push!(Γ, el)\n",
    "    end\n",
    "    return Γ\n",
    "end\n",
    "\n",
    "function plot_element(el; plot_with_normal=false)\n",
    "    # create a array of vectors\n",
    "    xis = linspace([-1.0], [1.0], 3)\n",
    "    time = 0.0\n",
    "    coords = interpolate(el, \"geometry\", xis, time)\n",
    "    ncoords = interpolate(el, \"geometry\", xis, time)\n",
    "    normals = interpolate(el, \"normal\", xis, time)\n",
    "    xs = [X[1] for X in coords]\n",
    "    ys = [X[2] for X in coords]\n",
    "    plot(xs, ys, \"-\")\n",
    "    plot([xs[1], xs[end]], [ys[1], ys[end]], \"ko\")\n",
    "    plot([xs[1], xs[end]], [ys[1], ys[end]], \"ko\")\n",
    "    if plot_with_normal\n",
    "        for i=1:3\n",
    "            p0 = ncoords[i]\n",
    "            p1 = ncoords[i]+0.3*normals[i]\n",
    "            plot([p0[1], p1[1]], [p0[2], p1[2]], \"-k\")\n",
    "        end\n",
    "    end\n",
    "end\n",
    "\n",
    "Γ₁ = create_elements([x1 y1]', 0)\n",
    "Γ₂ = create_elements([x2 y2]', nsl);\n",
    "for el in [Γ₁; Γ₂]\n",
    "    calculate_normals!(el, 0.0)\n",
    "end\n",
    "\n",
    "figure(figsize=(10, 3))\n",
    "for el in Γ₁\n",
    "    plot_element(el; plot_with_normal=true)\n",
    "end\n",
    "for el in Γ₂\n",
    "    plot_element(el; plot_with_normal=false)\n",
    "end\n",
    "axis(\"equal\")\n",
    "ylim(-0.1, 4.0)\n",
    "#axis(\"off\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating surface normals\n",
    "\n",
    "To define so called \"continuous normal field\", normals must be unambiguous in nodes. This can be done by averaging normals of adjacent elements:\n",
    "\\begin{equation}\n",
    "\\mathbf{n}_{k}=\\frac{\\sum_{e=1}^{n_{k}^{\\mathrm{adj}}} \\mathbf{n}_{k}^{\\left(e\\right)}}{\\left\\Vert \\sum_{e=1}^{n_{k}^{\\mathrm{adj}}} \\mathbf{n}_{k}^{\\left(e\\right)}\\right\\Vert }.\n",
    "\\end{equation}\n",
    "\n",
    "In practice, we take average of all normal vectors connecting to some arbitrary node."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "JuliaFEM.Field{Array{Array{T,1},1}}(0.0,1,Array{T,1}[[0.7021480119880804,-0.7120310170639945],[0.7021480119880804,-0.7120310170639945]])"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FEM.interpolate(Γ₁[1][\"normal\"], 0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzgAAAEeCAYAAABG5rCAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlclXXe//HXOeyCoKKmqLgkmJlG4BY2htVkKriLUjYuUzbdmYWVltU4Zc6ULU7m6J39DEONu5EkBXE0nWyc3MWlTEXKrQAVDUUB2c7vjzMyIds5bAcO7+fjwUP7Xtf3uj7XzFF5810ug8lkMiEiIiIiImIHjLYuQEREREREpKYo4IiIiIiIiN1QwBEREREREbuhgCMiIiIiInZDAUdEREREROyGAo6IiIiIiNgNBRwREREREbEbCjgiIiIiImI3FHBERERERMRuKOCIiIiIiIjdcLTm5CNHjvCnP/2JpKQk0tPTcXV1pVu3bjz11FM88sgjFfZdsWIFU6dOLfNYeno6rVu3LtWekZHBpk2b6NSpE25ubtaUKiIiIiIidiQnJ4dTp04xePBgWrZsWe55VgWcM2fOcPXqVSZPnoyPjw/Z2dnExsby6KOPcurUKV5++eVKrzFv3jw6d+5cos3Ly6vMczdt2sTEiROtKVFEREREROzYqlWrKhxcMZhMJlN1blBUVERQUBCXLl3i9OnT5Z53YwRn3759BAYGWnTtb775hnvuuYdVq1bRvXv36pTZIEVGRrJw4UJblyE2ps+B6DMgoM+BmOlzINB4PwdHjx5l4sSJ/Pvf/2bAgAHlnmfVCE5ZjEYj7du3Jysry6LzTSYTWVlZNGnSBAcHhwrPvTEtrXv37haHInvi5eXVKJ9bStLnQPQZENDnQMz0ORDQ56CypStV2mQgOzubjIwMfvjhBxYuXMimTZuYNWuWRX0HDRqEl5cX7u7ujBgxgpSUlKqUICIiItIoHD16FD8/P1xcXEhISMDFxQU/Pz+OHj1q69JE6qUqjeDMnDmTZcuWmS/g6MiiRYuYNm1ahX3c3d2ZMmUKgwYNwtPTk3379vHee+8RHBxMUlIS7du3r0opIiIiInbryJEjBAQEUFBQUNyWl5dHSkoKvXr14vDhw41yGr9IRaoUcCIjIwkPDyc1NZXVq1czffp03NzcmDRpUrl9xo0bx7hx44r/e/jw4QwePJiBAwcyf/58li5dWpVSREREROzWyJEjS4SbXysoKGD48OGcOHGijqsSqd+qFHC6detGt27dAJg4cSKDBw/m2WefJTw83KrtnAcMGEC/fv3YsmVLVcqwexEREbYuQeoBfQ5EnwEBfQ4aqzNnzlTruNgn/X1QsWrvogawbNky/vCHP5CUlERAQIBVfcPDw/nnP/9JRkZGqWNJSUkEBQUxcODAUltJR0RE6P9cERERsWvOzs7k5+eXe9zJyYm8vLw6rEikbsTExBATE1Oi7fLly/zrX/9i//79FW6yUO1d1MD80h0w76hmrR9//JFWrVpVeM7ChQsb9U4RIiIi0jgZDIZqHRdpqMoazLgx+FEZqxLJhQsXSrXl5+cTHR2Nt7c3PXr0ACAtLY1jx46VmDNaVt/ExESSkpJ46KGHrClDRERExO7l5oKDg2+F5/j6VnxcpDGyagRn2rRpZGVlMXDgQHx8fEhPT2f16tUkJycTFRVV/F6bl156iejoaE6dOlX8By84OJjAwECCgoLw8vIiKSmJjz/+GF9fX+bMmVPzTyYiIiLSQOXkwIgRUFS0HgeHXhQWlr3RwGOPPVbHlYnUf1YFnAkTJrB8+XKWLl3KxYsX8fT0pF+/fixevJj777+/+DyDwVBqyHTChAls2LCBzZs3k52djY+PD0888QRz586tdIqaiIiISGORnQ1hYbBrF2zc2J02bQ4zfPhwzpw5g8lkwmAw0KFDB+68807mzJmDr6+v1iWL/EqNbDJQW27Ms6tsIZGIiIiIPbh2DUJDYe9eSEyEgQPLP7eoqIipU6eycuVKYmJiCA8Pr7tCRWzA0mxQI5sMiIiIiEj1ZGXBsGFw4AD84x9wzz0Vn280Glm+fDmFhYU8/PDDODg4MGbMmLopVqQes37bMxERERGpUVeuwJAhcPAgbN5cebi5wcHBgRUrVhAeHs6ECROIi4ur3UJFGgAFHBEREREbunwZHnoIvv0WvvwS7r7buv4ODg5ER0czevRowsPDWb9+fe0UKtJAKOCIiIiI2EhmJjz4IBw9Clu2QL9+VbuOo6Mjq1atYsSIEYwdO5YNGzbUbKEiDYgCjoiIiIgN/PIL/Pa3cOIEbN0KffpU73pOTk7ExMQwbNgwRo8ezT/+8Y+aKVSkgVHAEREREaljFy/C/ffDyZPwz39CTW0W6+TkxGeffcbgwYMZOXIkmzdvrpkLizQgCjgiIiIidSgjwxxuzp41h5uAgJq9vrOzM2vWrOGBBx5gxIgRbN26tWZvIFLPKeCIiIiI1JELF+C++yA1Fb76Cnr1qp37uLi4EBsbS0hICGFhYezZs6d2biRSDyngiIiIiNSBc+dg0CA4fx62bYM77qjd+7m6uhIXF0dkZCQ9evSo3ZuJ1CN60aeIiIhILUtPN4/cZGaaw81tt9XNfV1dXZk/f37d3EyknlDAEREREalFqanmcJOVZQ43/v62rkjEvingiIiIiNSSn382T0vLyYGvv4auXW1dkYj9U8ARERERqQVnz5rDTX6+Odx06WLrikQaBwUcERERkRp2+rQ53JhM5nDTqZOtKxJpPLSLmoiIiEgNOnkS7r3X/Ptt2xRuROqaAo6IiIhIDfnhBwgJAUdH88hNx462rkik8VHAEREREakBKSnmcOPiYg43HTrYuiKRxkkBR0RERKSakpPN09Lc3c3hpl07W1ck0nhZHHCOHDnCuHHjuPXWW3F3d8fb25vg4GBWr15tUf/MzEymTZtGq1at8PDw4L777uPAgQNVLlxERESkPjh2zBxuvLzMa27atrV1RSKNm8W7qJ05c4arV68yefJkfHx8yM7OJjY2lkcffZRTp07x8ssvl9u3qKiIYcOGcfjwYWbNmoW3tzdLliwhJCSE/fv301WbwouIiEgD9P335pd4tmoFW7dC69a2rkhEDCaTyVTVzkVFRQQFBXHp0iVOnz5d7nl///vfmTBhArGxsYwePRqAjIwM/P39GTJkSLmjQElJSQQFBbF//34CAwOrWqaIiIhIjfvuO3O4adsWtmwxhxwRqT2WZoNqrcExGo20b98eJyenCs+LjY2lTZs2xeEGoGXLloSHh7Nu3Try8/OrU4aIiIhInTp0yPyem3btzCM3Cjci9YfVASc7O5uMjAx++OEHFi5cyKZNm5g1a1aFfQ4cOFBmyurTpw/Z2dkkJydbW4aIiIiITRw4YB658fU1h5uWLW1dkYj8mtUBZ+bMmbRu3Ro/Pz9mz57NokWLmDZtWoV90tLSaFvGirsbbampqdaWISIiIlLnVq3azz33fEiXLuZpaS1a2LoiEbmZxZsM3BAZGUl4eDipqamsXr2a6dOn4+bmxqRJk8rtk5ubi4uLS6l2V1dXAHJycqwtQ0RERKRO7d0Ljz22nuvX5/H00y1p3nyMrUsSkTJYHXC6detGt27dAJg4cSKDBw/m2WefJTw8HDc3tzL7uLm5cf369VLtubm5xcdFRERE6qtdu2DwYLjrrrn4+BzniScm0q1be/r162fr0kTkJlYHnJuNGTOGL7/8kuPHjxMQEFDmOW3bti1zGlpaWhoAPj4+Fd4jMjISLy+vEm0RERFERERUsWoRERERy+zYAQ89BHfeCYmJRpycVvDAAw8QFhbG7t276dy5s61LFLE7MTExxMTElGi7fPmyRX2rHXBuTC8zGstfzhMQEMD27dsxmUwYDIbi9t27d+Pu7o6/v3+F91i4cKG2iRYREZE69+9/w5AhEBgIGzaAhweAK1988QV33303Q4cOZceOHTRv3tzWpYrYlbIGM25sE10ZizcZuHDhQqm2/Px8oqOj8fb2pkePHoB5VObYsWMUFBQUnzd27FjOnTvH2rVri9syMjJYs2YNYWFhlW4zLSIiIlLXvv7aPHLTpw8kJt4IN2YtW7YkMTGR8+fPM2bMGPLy8mxXqIiUYPEIzrRp08jKymLgwIH4+PiQnp7O6tWrSU5OJioqCgcHBwBeeukloqOjOXXqFL6+voA54PTv358pU6bw/fff4+3tzZIlSzCZTLz22mu182QiIiIiVfTVVxAaCsHBsG4dNGlS+hw/Pz+++OILHnjgAaZNm0ZUVFSJmSoiYhsWB5wJEyawfPlyli5dysWLF/H09KRfv34sXryY+++/v/g8g8FQ6g+30WgkMTGRF154gUWLFpGTk0Pfvn2Jjo7Gz8+v5p5GREREpJq2bIHhw+E3v4EvvoCK9kL6zW9+Q1RUFI888gi33norr776at0VKiJlMphMJpOtiyjPjXl2+/fv1xocERERqXWbN8OIERASAnFx8J83WlTqjTfe4NVXX2XVqlU88sgjtVqjSGNlaTao9iYDIiIiIvZg40YYNQoeeAA+/xzKeIVfuV5++WV++OEHpk6dSocOHRg4cGDtFSoiFbJ4kwERERERe5WQACNHmt91Y224AfMU/Q8//JB77rmHkSNHcvz48dopVEQqpYAjIiIijdq6dTB6NAwbBmvWWB9ubnB2dubzzz+nTZs2DB06tMwdaEWk9ingiIiISKMVFwdjx5rX3Xz2GTg7V+96zZo1IzExkatXrzJixAhyc3NrplARsZgCjoiIiDRKsbEwbhyMGQOffgo19Vq+Tp06ER8fz8GDB5k0aRJFRUU1c2ERsYgCjoiIiDQ6n30GEybA+PGwalXNhZsb+vbty+rVq1mzZg0vv/xyzV5cRCqkgCMiIiKNyqefwsMPm7+io8GxlvaUHTVqFO+88w5vvvkmH330Ue3cRERK0TbRIiIi0misXAmTJ8Pvfgf/7/+Bg0Pt3i8yMpJz587RqVOn2r2RiBRTwBEREZFGYcUKmDrV/LVsGRjrYB6LwWDgrbfeqv0biUgxTVETERERu7d8uTnYPP543YUbEbEN/fEWERERu7ZsGTz2GDz5JCxdqnAjYu/0R1xERETs1pIl8MQT8PTTsHixwo1IY6A/5iIiImKXPvgAnnoKnn0W3n8fDAZbVyQidUEBR0REROzOwoUwYwY8/zy8957CjUhjooAjIiIiduWdd2DmTJg9GxYsULgRaWwUcERERMRuvPkmvPACvPwy/OUvCjcijZECjoiIiNiF+fPhpZdg7lyYN0/hRqSxUsARERGRBu+11+CVV8y//ulPCjcijZkCjoiIiDQICQkJODk5YTAYir+cnJyYMCGBP/3JPILzxz/aukoRsTWLA87evXuZPn06PXr0wMPDg44dOzJ+/HhOnDhRad8VK1ZgNBrL/Dp//ny1HkBERETs37p16wgLC6OgoKBEe0FBAZ99FsbvfpfAnDk2Kk5E6hVHS09866232LlzJ+PGjaNXr16kpaWxePFiAgMD2bVrFz169Kj0GvPmzaNz584l2ry8vKyvWkRERBqVsWPHVnj8009H8ckn+XVUjYjUZxYHnOeee44+ffrg6PjfLuPHj6dnz568+eabrFy5stJrDBkyhMDAwKpVKiIiIo3WzSM31h4XkcbD4ilqd999d4lwA9C1a1duv/12jh07ZtE1TCYTWVlZFBYWWleliIiIiIiIBaq1yYDJZOLcuXO0bNnSovMHDRqEl5cX7u7ujBgxgpSUlOrcXkREREREpASLp6iVZfXq1aSmpvLGG29UeJ67uztTpkxh0KBBeHp6sm/fPt577z2Cg4NJSkqiffv21SlDRERE7NTl3Mv8efufoTnwS/nn3TzLREQaryr/bXDs2DGeeuopgoODmTRpUoXnjhs3jnHjxhX/9/Dhwxk8eDADBw5k/vz5LF26tKpliIiIiB0qKCpg2f5lzN02l+z8bCL+GEFMZEy558fFxdVhdSJSn1Vpilp6ejrDhg2jefPmxMbGYqjC27QGDBhAv3792LJlS1VKEBERETtkMplIPJFIr6W9mJ44nVD/UJKnJ/Pps58SHx9faqTmxvcgW7dutUW5IlIPWT2Cc/nyZYYMGcKVK1fYvn07bdq0qfLN27dvT3JycqXnRUZGltpOOiIigoiIiCrfW0REROqXw+cO8/zm5/nyxy8J6RTCqtGrCGz7391XQ0NDyc8vuRV0UVERnTp14q9//SuhoaHcf//9dV22iNSCmJgYYmJKjtpevnzZor4Gk8lksvRGubm5PPjggxw4cIAtW7bQr18/6yq9Se/evbl27RpHjx4t83hSUhJBQUHs379f20uLiIjYqfSr6bz6z1f5+ODH3Nr8Vt558B3C/MMsniHy008/0aVLFxwdHUlNTaVZs2a1XLGI2IKl2cDiKWqFhYWMHz+e3bt3s2bNmnLDTXp6OseOHSuxH/2FCxdKnZeYmEhSUhIPPfSQpSWIiIiIHcnJz2H+v+bj94Efa4+tZeHghXz3P98xvNtwq6a/t2/fnpUrV5KTk8OAAQNqsWIRaQisetFnfHw8YWFhZGRksGrVqhLHJ06cCMCLL75IdHQ0p06dwtfXF4Dg4GACAwMJCgrCy8uLpKQkPv74Y3x9fZkzZ04NPo6IiIjUd0WmIj799lPmbJ1D+tV0nu77NK8MfIXmbs2rfM3x48ezfv16Pv30U2bMmMGiRYtqsGIRaUgsDjiHDh3CYDAQHx9PfHx8iWMGg6E44BgMhlI/dZkwYQIbNmxg8+bNZGdn4+PjwxNPPMHcuXNp1apVDTyGiIiINATbT29n5uaZ7Evdx+juo3nrgbfo2qJrjVx75cqVfPPNN3zwwQcMGzaMwYMH18h1RaRhsWoNTl3TGhwRERH78MOlH5i9ZTafH/2coLZBvDf4PQZ2HFjj90lNTaVTp044ODjw888/06JFixq/h4jYRo2vwRERERGx1i85v/Dcpufo/rfu7P55NytHrWTP43tqJdwA+Pj48Omnn5Kbm6v1OCKNlAKOiIiI1Lj8wnw+2P0BXT/oyof7P+SP9/6R49OPM7HXRIyG2v32Y+zYsTz66KPFLyUXkcZFAUdERERqjMlkYv3x9dyx9A6e+cczjL5tNCeePsErA1+hiVOTOqtjxYoVdOrUiSVLlrBx48Y6u6+I2J4CjoiIiNSIA2kHuD/6fkb83wg6eHbgwBMH+Gj4R7Rt2rbOazEajezcuRNnZ2dGjx5NRkZGndcgIrahgCMiIiLVkpqVytR1UwlaFkTa1TQSIhL48tEvubPNnTatq02bNlqPI9IIWbxNtIiIiMivXcu7xjs73mHBjgU0cWrC4qGLeTzwcZwcnGxdWrExY8YwefJkvv76a3Jzc3F1dbV1SSJSyxRwRERExCpFpiJWHlrJnH/OISM7g2f6PcOc38yhmWszW5dWpuXLlwPmaWsiYv8UcERERMRi205t47nNz5GUlkR4j3DevP9NOjfvbOuyKqRgI9K4KOCIiIhIpZIvJjPry1msO76Ofu368c3UbwjuEGzrskRESlHAERERkXL98PMPvL3jbZYfW45PUx9ixsQwvsd4DAaDrUsTESmTAo6IiIiUkleYxwc7P2DWiFk4+Dsw7+15PNPvGdyc3GxdmohIhTQpVURERIqZTCbijsbRY0kPZv1zFn1H9KVobxEjvEco3IhIg6CAIyIiIgDsT91PyCchjP77aLq26MqhPxxi29Jt+Pr68sILL9i6PBERiyjgiIiINHI/XfmJ38X9jt4f9eZi9kU2PrKRjY9s5I7Wd+Di4sKCBQvYsGEDW7ZssXWpIiKV0hocERGRRupq3lUWfLOAd3a8Q1OXpnwY+iFT75qKo7HktwdjxoxhwIABPPfccyQlJeHg4GCjikVEKqcRHBERkUamsKiQ5UnL8fvAjwXfLODZ/s9y4ukTTAuaVircABgMBt577z0OHz7MihUr6r5gERErKOCIiIg0Ilt+3ELgskAei3+M+zrfx/Hpx/nz/X/G08Wzwn59+/bl4Ycf5pVXXiErK6uOqhURsZ4CjoiISCNw9MJRQj8N5bcrf4uHswe7fr+L1aNX07FZR4uv8Ze//IXMzEwWLFhQi5WKiFSPxQFn7969TJ8+nR49euDh4UHHjh0ZP348J06csKh/ZmYm06ZNo1WrVnh4eHDfffdx4MCBKhcuIiIilbtw7QLTE6fTc2lPvr/wPWvGreHfU/5Nv/b9rL6Wr68vM2fO5J133uHs2bO1UK2ISPVZHHDeeust4uLi+O1vf8uiRYuYNm0a//rXvwgMDOTIkSMV9i0qKmLYsGHExMQwY8YMFixYwPnz5wkJCSElJaXaDyEiIiIlXS+4ztvfvE3XD7qy6vAq3nzgTY4+dZSxt4/FYDBU+bovvvginp6evPzyyzVYrYhIzbF4F7XnnnuOPn364Oj43y7jx4+nZ8+evPnmm6xcubLcvrGxsezcuZPY2FhGjx4NQHh4OP7+/sydO5fVq1dX4xFERETkBpPJROz3sczeMpszl8/wZO8nmRsyl5ZNWtbI9Zs2bcq8efN44oknmDFjBr17966R64qI1BSLR3DuvvvuEuEGoGvXrtx+++0cO3aswr6xsbG0adOmONwAtGzZkvDwcNatW0d+fr6VZYuIiMjNdv+0m3ui7iE8NpwerXvw3f98xwdDP6ixcHPD1KlTueOOO5g5cyYmk6lGry0iUl3V2mTAZDJx7tw5Wras+C/OAwcOEBgYWKq9T58+ZGdnk5ycXJ0yROzO0aNH8fPzw8XFBWdnZ1xcXPDz8+Po0aO2Lk1E6qHTmad5+POH6b+8P1fzrvLlo18SHxHPbS1vq5X7OTo68u6777J9+3bi4uJq5R4iIlVVrYCzevVqUlNTGT9+fIXnpaWl0bZt21LtN9pSU1OrU4aIXTly5Ai9evUiJSWFvLw88vPzycvLIyUlhV69einkiEixK9evMGfrHLot7sZXp75i+fDlJE1L4oEuD9T6vR988EGGDBnCrFmzyMvLq/X7iYhYqsoB59ixYzz11FMEBwczadKkCs/Nzc3FxcWlVLurqysAOTk5VS1DxO6MHDmSgoKCMo8VFBQwfPjwOq5IROqbgqICPtz3IX4f+PHXXX9l1oBZnHj6BFPvmoqD0aHO6nj77bc5efIkf/vb3+rsniIilalSwElPT2fYsGE0b96c2NjYSndjcXNz4/r166Xac3Nzi4+LiNmZM2eqdVxE7NumlE0E/G8Af9jwBx7q+hDJTyfz+qDX8XD2qPNaevTowbRp03j99de5ePFind9fRKQsFu+idsPly5cZMmQIV65cYfv27bRp06bSPm3bti1zGlpaWhoAPj4+FfaPjIzEy8urRFtERAQRERFWVC7SMFS2YFcLekUap+/Of8fzm59n0w+bGNhxIHsf30tvH9vvYPbaa6+xevVqXn/9dd5//31blyMidiImJoaYmJgSbZcvX7aor1UBJzc3l7CwMFJSUtiyZQu33WbZ4sWAgAC2b9+OyWQqMdqze/du3N3d8ff3r7D/woULy9ykQMQeVTYiaso3kTQgCe8wb7xDvXHv4V6td1qISP127uo55m6by0dJH9GleRfixscxotuIevPnvnXr1rz88svMmzePP/3pTzRv3tzWJYmIHShrMCMpKYmgoKBK+1o8Ra2wsJDx48eze/du1qxZQ79+Zb8BOT09nWPHjpVYQzB27FjOnTvH2rVri9syMjJYs2YNYWFhODk5WVqGiN3z9fWt8Hj7Vu1xbu3M6TdOs6/nPnZ13kXyU8lc3HiRwtzCOqpSRGpKQkICTk5OGAyG4i8nJyfWrlvLX7b/Bb8P/Pj7kb/z7oPvcuR/jjDytpH1Jtzc8Mwzz/Dtt98q3IhIvWDViz7j4+MJCwsjIyODVatWlTg+ceJEwPyG4+joaE6dOlX8jdrYsWPp378/U6ZM4fvvv8fb25slS5ZgMpl47bXXavBxRBq+9evX07NnTwoLS4cVR0dHEr9OpHv37hRdLyLz60wuJlzkYvxFUpekYmxipPkDzfEO9cZ7mDcuPqU39xCR+mPdunWMHDmyVHtBQQFjRo6BR+DZic/y6r2v0sKthQ0qtIyrqyudO3e2dRkiIoAVAefQoUMYDAbi4+OJj48vccxgMBQHnBs/ffo1o9FIYmIiL7zwAosWLSInJ4e+ffsSHR2Nn59fDTyGiP3o3r07b7/9NjNnziwe3TQYDPj6+rJ+/Xq6d+8OgNHFSIsHW9DiwRZ0fb8r2UezzWEn4SLJf0iGIvAI9DCHnVBvmgY1xWCsXz/1FWnsxo4dW+Fxh88cWLhqYR1VIyJiHwymerxi+cY8u/3792sNjjQqkyZNIikpiW+//bZK/fMv5nNp0yUuJlzk0sZLFGQW4HSLE97DzGGn+QPNcWxq9R4jIlLDLJlqVo//mRYRqVOWZgN9hyNSzxQWFpKYmMhjjz1W5Ws4eTtxy8O3cMvDt1BUUMSVHVeKR3fSP07H4GygWUiz4qlsbl20VbtInbh0Cf71L9i2DbZtwwnIt3VNIiJ2RgFHpJ7Zu3cvGRkZhIaG1sj1jI5Gmg1sRrOBzbh1wa3k/JDDxQ3msPPDcz+QMiOFJrc3KZ7K5nm3J0bHKr8DWER+7aZAw+HDYDJBp04QEkLhoUM2LlBExP4o4IjUMwkJCbRo0YL+/fvXyvXdbnWj/Yz2tJ/RnoKsAn758hfzyM6KdM4uOItjc0daDGmB9zBvWjzUAqcW2uVQxGKVBBoiI+Hee83/DRhXraLoV7uO3szRUf9Mi4hYS39zitQzCQkJDB06FAcHh1q/l2NTR1qNbkWr0a0wFZnI2pdVPJXt/KfnwQheA7zMozth3jS5rUm9255WxKasDDQ3i4uLIywsrNzLx8XF1ULRIiL2TQFHpB45e/Yshw4d4qWXXqrzexuMBjz7euLZ15POr3fm+s/Xi6eynfrTKX6c/SOunV2Lp7I1u7cZRhdNZZNGppqB5mahoaHEx8czatSoEu+PA1i7dm2NTVUVEWlMFHBE6pENGzbg4ODA4MGDbV0KLu1c8Jnmg880HwpzCsncZn7nTsYXGfz8wc8Y3c1K3GpcAAAgAElEQVTbVF8KvET7Ye3pclcXW5csUvNqONCUJTQ0lPz8/2418MknnzB58mQOHz7MqFGjqvkAIiKNj7aJFqlHQkNDuXbtGl999ZWtSymXyWTi2nfXiqey/c+O/8GEiSV9lhSP7njc5aGpbNIwVRZoQkKqHWgqU1RUhIeHB82aNSM1NbXW7iMi0tBom2iRBiY7O5utW7fyxhtv2LqUChkMBjx6euDR04MOsztwvMVxHr//cVwdXTn77llOzT2Fs4/zf9+5c39zHNxrfz2RSJXUwQiNtYxGI8OHD+ezzz4jKSlJP+ATEbGSAo5IPfHPf/6T3NzcBjXn/vjx42RezuTBJx+kxwM9KMov4vK/LxeP7qR9lIbBxUDz+5oXv3PHtaOrrcuWxqweBpqyLFiwgM8++4zZs2fz5Zdf2rQWEZGGRgFHpJ5ISEiga9eu+Pv727oUi+3YsQOj0Ujfvn0BMDoZaT6oOc0HNafru13JTs4u3qgg5ZkUTjx1Avee7v99504/TwwOmsomtejSJY59+in/LyqKNwsKcPz223oZaG7m6+uLn58f27Zto6CgQNtFi4hYQX9jitQDJpOJhIQExo0b16DWruzYsYOePXvi6elZ5vEm/k1o4t+EDpEdKLhcwKXNl4pHds785QyO3o54DzWP7DQf3BynZnrnjlRTGSM0l0wm3gUmhIbSe+bMehloyjJr1iwef/xx3n77bZvsrCgi0lBpj1eReuDQoUP8/PPPDWp6GsDOnTsJDg626FxHL0daj2tN90+6E5wezF077sLnCR+uHrrK9xO+55uW33Bw0EHOvnuW7OPZtVy52I1Ll+CLL+DZZyEgAFq2hFGjYN06uOsuiIqi9/HjNGnShG333guTJjWIcAMwdepUXFxc+Nvf/mbrUkREGhQFHJF6ICEhgaZNm/Kb3/zG1qVY7NKlSxw9epS7777b6r4GBwNed3vRZX4X+hzqQ//T/fH7wA+ju5GTr5xkz2172O23m5TIFH7Z+gtFeUW18ATSIFkQaDh50vwVFQWTJuHs78+AAQPYtm2brau3itFoJDQ0lJ9//pnvvvvO1uWIiDQYmqImUg/cfvvtvPjiizg7O9u6FIvt2rULwOIRnIq4+rrS7sl2tHuyHYXZhfzyz1+4mHCR82vO89Nff8KhqQMtBrfAO9SbFkNa4Ny64fzvJNVUQ5sChISE8NZbbzW49Sxvv/02n3/+OS+88AIbN260dTkiIg1Cw/lbXsSOjR492tYlWG3nzp20bt2aLl1q9gWfDk0caBnakpahLTGZTFw9dLV4V7ZjU44B4NnPs3ijAvde7g1q3ZJUopZ2OQsJCeHll1/m4MGD9O7du+brriWdO3emS5cubNmypcGFMxERW9EUNRGpkh07dhAcHFyr4cJgMNA0oCmdXulE0K4ggtOC6fZxN5zbOXPmzTPsC9jHLt9dJD+ZzMUNFynMKSzRPyEhAScnJwwGQ/GXk5MTCQkJtVazWKkKU86qsoamd+/e5nU4DWyaGsBzzz1HQUEBCxcutHUpIiINgsFkMplsXUR5LH1bqYjUrYKCApo1a8Yf//hHZs2aZZMaiq4Xkbk90zy6E3+R3B9zMboZaX6/+Z07uxx3MeaxMeX2j4+Pb3CbOtiFykZoQkJqbZezBx98EGdn5wYXcIuKinBzc+OWW27hzJkzti5HRMRmLM0GGsEREat9++23XLt2rUbW31SV0cVIiwda4PdXP/ql9KPP0T50er0TBVkFJD+VzNOPPV1h/1GjRtVRpfbJ4tGxOhqhsURISAjbt2+noKCgVq5fW4xGI0OHDuXs2bMcPXrU1uWIiNR7mswrIlbbuXMnTk5OBAUF2boUwDyVzf02d9xvc8f3eV/yf8knrUVahX0KCgrotXcvTgYDzkaj+VeDASej0fzrf9rLaiuvT0XXublPeddxAhwNBozG+vvzp3Xr1jFy5MhS7QUFBYSFhRE/Zw6h167V6BqamtBQ1+EALFiwgPPnz9frz4WISH1hVcC5du0aCxYsYPfu3ezZs4fMzEyioqKYNGlSpX1XrFjB1KlTyzyWnp5O69atrSlFRGxox44dBAYG4ubmZutSyuTU3AkTlc++HdSsGXkmE3lFReT/+leTiauFheQVFJRoyy8qMv96U58bbQU1MeP3229hxgwcP/kEl06dyg1KVQlPlYWwioLZr9vGjh9f4SOM+vOfya8HgeZmv16H09ACjp+fH998842tyxARaRCsCjgXLlxg3rx5dOzYkYCAALZt22b1AuN58+bRuXPnEm1eXl5WXUNEbGvHjh2MGDHC1mVU2/t+fjV6vaL/hJxfh6KywtONoFRW2+GLF3kLeKVrV7w6dKgwUJV17azCQvLLCWblhTWrg1lhYYWHC8A85ayecXZ2Ln4fzvPPP2/rckREpJZYFXB8fHyKR1v2799Pnz59rL7hkCFDtGGASAOWnp7OyZMnbbr+xhKOjo4VrrWoje12jb8aQakqbw8P3gKmduhAhw4daq64ChT9J+hUNkJ1o+3eBraG5dca6vtwRETEclb9K+zs7Fw8layqm6+ZTCaysrIorOQngCJSP+3cuROAu+++28aVVCwuLq5ax20lPz8foE5f+mo0GHAxGvFwdKSFkxO3ODvTwdWVLm5u3ObuTi8PD4KaNuVuLy8GNmtWZ3XVhpCQEK5cucLBgwdtXYqIiNSSOl+tOGjQILy8vHB3d2fEiBGkpKTUdQkiUg07duzA19eX9u3b27qUCoWGhhIfH4+Dg0OJdkdHx3q9RXReXh4ATk5ONq6kfJWNfNTnkZGG/D4cERGxTJ0FHHd3d6ZMmcKSJUv44osvmDVrFlu3biU4OJiffvqprsoQkWrauXNnvR+9uSE0NJQFCxYAkJiYiMlkIj8/v96GG7DNCI61GuroGJRchyMiIvapzgLOuHHjWL58ORMnTmT48OG8/vrrbNq0iYsXLzJ//vy6KkNEqiEvL499+/bV+/U3v7Z3714A7r33XhtXYpmGMIJzY3Ts5pGa+j46dkNDfR+OiIhYxqYb6g8YMIB+/fqxZcsWW5YhIhZydnbm5MmTTJw40dalWOzo0aM4OzvTpEkTW5dikRsjOPU54IA55OTn5/Pggw9iMBgaxOjYDVqHIyJi32w+Ubp9+/YkJydXeE5kZGSpraQjIiKIiIiozdJEpAxt27a1dQlW+emnn/D29rZ1GRbLy8vDwcGhwbzQsUuXLphMJs6fP99g3mfWkN+HIyLSWMTExBATE1Oi7fLlyxb1tXnA+fHHH2nVqlWF5yxcuFBbS4tIlWRmZtK3b19bl2GxvLy8er3+5ma33XYbAPv372fIkCE2rsYyeh+OiEj9V9ZgRlJSEkFBQZX2rZUfEaanp3Ps2LES85svXLhQ6rzExESSkpJ46KGHaqMMEWnkzp8/T2FhIXfccYetS7FYfn5+gwo4PXv2BODbb7+1cSXWmTx5coNZlyUiItaxegRn8eLFZGZmkpqaCsD69es5c+YMADNmzMDT05MXX3yR6OhoTp06ha+vLwDBwcEEBgYSFBSEl5cXSUlJfPzxx/j6+jJnzpwafCQREbOtW7cC9f+dPb+Wl5dX79ff/NqNKV7Hjh2zcSXWefjhh21dgoiI1BKrA867777L6dOnATAYDMTFxbF27VoMBgO/+93v8PT0xGAwYDAYSvSbMGECGzZsYPPmzWRnZ+Pj48MTTzzB3LlzK52iJiJSFTt27ADg/vvvt3EllmtoIzg3/s6/8e+CiIiIrVkdcE6ePFnpOVFRUURFRZVomzdvHvPmzbP2diIiVXbo0CGMRmPxSHJD0NBGcABcXV2LR/VFRERsrWFs0yMiUgUnT57E09PT1mVYpaGN4AB4eXmRkZFh6zJEREQABRwRsWMZGRm0a9fO1mVYpSGO4LRq1YqsrCxblyEiIgIo4IiIncrLyyM3N5du3brZuhSrNLRtosH8PrPr16/bugwRERFAAUdE7NSNDQYa2oscG+IUtc6dOwNoHY6IiNQLCjgiYpe+/vprAEJCQmxbiJUa4hS17t27A+YXsImIiNiaAo6I2CUvLy9uueUW+vTpY+tSrNIQR3ACAgKAhveyTxERsU8KOCJil5599lnS09NxdLR6N3yb6t27N/fee6+ty7DKjYBz/PhxG1ciIiJShffgiIhI7YmMjLR1CVbz8PCgc+fODep9QyIiYr8UcEREpNp+/PFHW5cgIiICaIqaiIg0UgkJCTg5OWEwGIq/nJycSEhIsHVpIiJSDQo4IiLS6Kxbt46wsDAKCgpKtBcUFBAWFqaQIyLSgCngiIhIozN27NgKj48aNaqOKqlZGpUSEdEaHBERaQRycuD0aTh1yvxVUGCq8PyCggJGjx6Nu7s77u7uNGnSpNTvy2q7+fcODg518nxgHpUaOXJkmc8SFhZGfHw8oaGhdVaPiIitKOCIiEiDd3OAufnr3Ln/nmvOHAYLrplDRkYG2dnZXLt2jWvXrhX/Pjc316K6nJ2dLQ5DVfn9r7dBt2RUKj8/36K6RUQaMgUcERGp96wNML6+0KkT3H47DB1q/v2NLx8fcHIquPkWpWzcuLHcY4WFheTk5JQKPtb+/uLFi2Uer0qAunk90c0qOy4iYi8UcEREpEIJCQmMGjWqxDfIjo6OxMXF1diUJ2sDTIcO5rDSvTsMGVIywLRrB5W939XR0bHCb/gre0Gsg4MDHh4eeHh4VP5wVVBUVFQceG7+tbzfv/7665Ve99ixqTRt2hdPz764u/fEaHSqlfpFRGxJAUdERMpVU+s6qhJgOnaseoCpTFxcHGFhYRUetyWj0Wh1gLIk4Fy9eohz51ZiMhVgMLjQtOldxYGnadO+uLl1xWCofPqeiEh9poAjIiLlsnRdR12PwFRXaGgo8fHxtT4yVZcsGZXq3Xs/hYU5XL16kKysPVy5sodLlxL5+edF/zmnOU2b9ikOPE2b9sHFpU1dPYKISI2w6p+Qa9eusWDBAnbv3s2ePXvIzMwkKiqKSZMmWdQ/MzOTWbNmERcXR05ODn379uXdd9/lrrvuqlLxIiJSuyxZ19G2LaSn/7fNkjUwtR1gLBEaGmpXi+4tHZVycHDDy+tuvLzuLj6Wn3+JrKx9XLmyh6ysPaSmfkR+/hsAuLj4lgg8TZsG4ejYtHYfRkSkGqz6J+bChQvMmzePjh07EhAQwLZt2yweyi4qKmLYsGEcPnyYWbNm4e3tzZIlSwgJCWH//v107dq1Sg8gIiK29eST9TPANDbVGZVycmpBixYP0qLFgwCYTCauXz/7n8Czl6ysPZw+/TqFhVcBA02a3F4cerSeR0TqG6v+CfLx8SE9PZ3WrVuzf/9++vTpY3Hf2NhYdu7cSWxsLKNHjwYgPDwcf39/5s6dy+rVq62rXEREalx+YT47zu4gPjmehOQEcAAKK+7zxz/WSWligZoalTIYDLi6+uLq6kvr1uZpiiZTIdnZx4pHea5c2VO8nsdodMXD4+b1PLdqPY+I2IRVAcfZ2ZnWrVsD5p/uWCM2NpY2bdoUhxuAli1bEh4ezqpVq8jPz8fJST/9ERGpa7/k/MLGlI0kJCewMWUjmbmZtPVoS6h/KClOKRQWlp9wKtttTOyHweCAu3sP3N170LbtFIBfrefZ+6v1PO8DpdfzeHr2xdn5Fls+gog0EnX2L9OBAwcIDAws1d6nTx+WLVtGcnIyPXr0qKtyREQaLZPJxPGLx0lITiA+OZ5vznxDoamQoLZBPNPvGcL8w7ir7V0YDUaGrxler3cbE9uybD3PsjLX83h69sXDIwhHx4p3ijt69CjDhw/nzJkzmEwmDAYDvr6+rF+/nu7du9fq84lIw1RnASctLY2QkJBS7W3btgUgNTVVAUdEpJbkFebx7zP/Jv54PAknEki5lIKboxsPdHmApcOWMtRvKO0825XqZ4+7jUntqng9j3lq23/X8xhxd7/9V1Pb+pRYz3PkyBECAgJKbXaRkpJCr169OHz4sEKOiJRSZwEnNzcXFxeXUu2urq4A5OTk1FUpIiKNQkZ2BhtPbCQ+OZ5NP2ziyvUrtGvajjD/MP46+K/c1/k+3JzcKr2Ove02JnXL8vU80aXW84SG/l+5O/kVFBQwfPhwTpw4UZePIyINQJ0FHDc3N65fv16qPTc3t/i4iIhUnclk4vsL3xdvELDzp50UmYro264vz9/9PGHdwrjzlju18FtsruL1PDfez7OBn346V+F1zpw5UxflikgDU2cBp23btqSmppZqT0tLA8w7tJUnMjISLy+vEm0RERFERETUbJEiIg3M9YLrfH36axKSE0hITuBk5kncndz57a2/5aOwjxjqN5Q2HnpRo9R/Za3nMRicgfJHD63d8EhEGo6YmBhiYmJKtF2+fNmivnUWcAICAti+fXvxAsEbdu/ejbu7O/7+/uX2XbhwYZkbFIiINEbnr50n8UQi8cnxbP5hM1fzruLr5UuYfxih/qGEdArB1dHV1mWKVFtlo40ajRSxX2UNZiQlJREUFFRp31oJOOnp6WRmZtK1a9fiLUTHjh1LbGwsa9euZcyYMQBkZGSwZs0awsLCtEW0iEg5TCYTh88dNo/SnEhg90+7Aejfvj8v3fMSYf5h3NH6Dn2zJ3bH19eXlJSUCo+LiNzM6oCzePFiMjMzi6ebrV+/vngO7IwZM/D09OTFF18kOjqaU6dOFf/lM3bsWPr378+UKVP4/vvv8fb2ZsmSJZhMJl577bUafCQRkYYvtyCXr05+Vbye5uyVs3g4ezD41sFEjYhiiN8QWru3tnWZIrVq/fr19OrVq8yNBhwdHVm/fr0NqhKR+s7qgPPuu+9y+vRpwDw0HBcXx9q1azEYDPzud7/D09MTg8FQ6ieJRqORxMREXnjhBRYtWkROTg59+/YlOjoaPz+/mnkaEZEGLC0rjQ0nNpCQnMCXP35Jdn42nZt1ZtRtowj1D2Vgx4G4OJbejVLEXnXv3p3Dhw/rPTgiYhWDqR6v0Lsxz27//v1agyMiDUpCQkKl744xmUwcSD9Q/MLNfan7MBqMBHcIJtQvlLBuYXRv2V1Tz0RERLA8G9TZJgMiIo3FunXrGDlyZKn2goICwsLCeOV/X+G8z3kSTiSQmpWKp4snD3V9iGf6PcNDXR+iZZOWNqhaRETEPijgiIjUsLFjx1Z4/I0n36Dr+10Z32M8Yf5h3ON7D04O2mhFRESkJijgiIjUsPLevF7MBCee1tvXRUREaoPR1gWIiNiVoiK0YkZERMR2FHBERGpCYSH83/9Br15ospmIiIjtKOCIiFRHQQGsWgU9ekBEBHTogJODQ4VdbrwAWURERGqeAo6ISFXk50NUFNx2Gzz6KPj7w+7dsHEj//fFFxV2jYuLq6MiRUREGh8FHBERa+TlwbJl5kAzdSr06gVJSbB+PfTtC0BoaCjx8fGlRmocHR2Jj48vfg+OiIiI1DzNkxARsURuLixfDm+9BT/9BOPGmUNNz55lnh4aGkp+fn4dFykiIiIawRERqUh2Nrz/PnTpAjNmwMCBcOQIfPZZueFGREREbEcjOCIiZbl6Ff73f+Htt+HiRZg4EebMMU9NExERkXpLAUdE5NeysuBvf4N334XMTJg8GV56yTyCIyIiIvWeAo6ICJjDzAcfwMKFcO0a/P73MHs2dOxo68pERETECgo4ItK4XboEf/0rLFoE16/D44/DrFnQvr2tKxMREZEqUMARkcbpwgV47z3zdLSCAnjySXj+eWjb1taViYiISDUo4IhI43LuHLzzDixZAgYDPPUUPPcctG5t68pERESkBijgiEjjkJoKCxbAhx+CszNERsKzz0LLlrauTERERGqQAo6I2LczZ8wv51y+HNzc4MUXze+zad7c1pWJiIhILVDAERH7dPIkvPkmREVB06bwxz/C9Ong6WnrykRERKQWGa05+fr168yePRsfHx+aNGlC//792bJlS6X9VqxYgdFoLPPr/PnzVS5eRKSUlBSYOhX8/CAuDubPh9OnzS/pVLgRERGxe1aN4EyePJnPP/+cyMhI/Pz8iIqKYujQoXz11VcMGDCg0v7z5s2jc+fOJdq8vLysq1hEpCzHjpnDzKefmjcMePtteOIJaNLE1pWJiIhIHbI44OzZs4fPPvuMd955h5kzZwLw6KOPcscddzBr1iy++eabSq8xZMgQAgMDq16tiMjNvvsO3ngD/v53aNcO3n/f/JJONzdbVyYiIiI2YPEUtdjYWBwdHZk2bVpxm4uLC7///e/ZuXMnP//8c6XXMJlMZGVlUVhYWLVqRURuOHgQxoyBnj1h1y5YutQ8PW36dIUbERGRRszigHPgwAH8/f3x8PAo0d6nTx8ADh48WOk1Bg0ahJeXF+7u7owYMYKUlBQryxWRxuLo0aP4+fnh4uKCs7MzLi4u+Pn5cXTNGhgxAu66yxxyli+HEyfM09FcXGxdtoiIiNiYxVPU0tLSaFvGG75vtKWmppbb193dnSlTpjBo0CA8PT3Zt28f7733HsHBwSQlJdG+ffsqlC4i9urIkSMEBARQUFBQoj0lJYVe4eEc7tiR7tHREBEBjtoMUkRERP7L4u8McnJycCnjp6Ourq7Fx8szbtw4xo0bV/zfw4cPZ/DgwQwcOJD58+ezdOlSa2oWETs3cuTIUuHmhgJguKMjJx59tG6LEhERkQbB4ilqbm5uXL9+vVR7bm5u8XFrDBgwgH79+lm0zbSINC5nzpyp+PjZs3VUiYiIiDQ0Fo/gtG3btsxpaGlpaQD4+PhYffP27duTnJxc6XmRkZGltpOOiIggIiLC6nuKSP1nMpmqdVxEREQatpiYGGJiYkq0Xb582aK+Fgecu+66i23btpGVlUXTpk2L23fv3g1AQECApZcq9uOPP9KqVatKz1u4cKG2lxZpRAwGQ7WOi4iISMNW1mBGUlISQUFBlfa1eIra2LFjKSwsZNmyZcVt169fJyoqiv79+9OuXTsA0tPTOXbsWIn58xcuXCh1vcTERJKSknjooYcsLUFEGglfX99qHRcREZHGy+IRnL59+zJu3Dheeuklzp8/z6233sonn3zCmTNniIqKKj7vxRdfJDo6mlOnThV/ExIcHExgYCBBQUF4eXmRlJTExx9/jK+vL3PmzKn5pxKRBm39+vX06tWrzI0GHB0dWb9+vQ2qEhERkYbAqv1Vo6OjefXVV1m5ciW//PILd955JwkJCdxzzz3F5xgMhlLTRyZMmMCGDRvYvHkz2dnZ+Pj48MQTTzB37lyLpqiJSOPSvXt3Dh8+zPDhwzlz5gwmkwmDwYCvry/r16+ne/futi5RRERE6imDqR6v1r0xz27//v1agyMiIiIi0ohZmg0sXoMjIiIiIiJS3yngiIiIiIiI3VDAERERERERu6GAIyIiIiIidkMBR0RERERE7IYCjoiIiIiI2A0FHBERERERsRsKOCIiIiIiYjcUcERERERExG4o4IiIiIiIiN1QwBEREREREbuhgCMiIiIiInZDAUdEREREROyGAo6IiIiIiNgNBRwREREREbEbCjgiIiIiImI3FHBERERERMRuKOCIiIiIiIjdsCrgXL9+ndmzZ+Pj40OTJk3o378/W7ZssahvZmYm06ZNo1WrVnh4eHDfffdx4MCBKhUtIiIiIiJSFqsCzuTJk1m4cCGPPvooixYtwsHBgaFDh/LNN99U2K+oqIhhw4YRExPDjBkzWLBgAefPnyckJISUlJRqPYCIiIiIiMgNjpaeuGfPHj777DPeeecdZs6cCcCjjz7KHXfcwaxZsyoMObGxsezcuZPY2FhGjx4NQHh4OP7+/sydO5fVq1dX8zFERERERESsGMGJjY3F0dGRadOmFbe5uLjw+9//np07d/Lzzz9X2LdNmzbF4QagZcuWhIeHs27dOvLz86tYvoiIiIiIyH9ZHHAOHDiAv78/Hh4eJdr79OkDwMGDByvsGxgYWKq9T58+ZGdnk5ycbGkZIiIiIiIi5bI44KSlpdG2bdtS7TfaUlNTa6WviIiIiIiIpSwOODk5Obi4uJRqd3V1LT5entzc3Cr3FRERERERsZTFAcfNzY3r16+Xas/NzS0+Xht9RURERERELGVxwGnbtm2ZU8nS0tIA8PHxqZW+AJGRkQwfPrzEV0xMjKWlN1iN4RmlcvociD4DAvociJk+BwKN43MQExNT6vv/yMhIi/pavE30XXfdxbZt28jKyqJp06bF7bt37wYgICCg3L4BAQFs374dk8mEwWAo0dfd3R1/f/8K771w4cIyNymwdzExMURERNi6DLExfQ5EnwEBfQ7ETJ8DgcbxOYiIiCj1jElJSQQFBVXa1+IRnLFjx1JYWMiyZcuK265fv05UVBT9+/enXbt2AKSnp3Ps2DEKCgpK9D137hxr164tbsvIyGDNmjWEhYXh5ORkaRkiIiIiIiLlsngEp2/fvowbN46XXnqJ8+fPc+utt/LJJ59w5swZoqKiis978cUXiY6O5tSpU/j6+gLmgNO/f3+mTJnC999/j7e3N0uWLMFkMvHaa6/V/FOJiIiIiEijZHHAAYiOjubVV19l5cqV/PLLL9x5550kJCRwzz33FJ9jMBhKTEMDMBqNJCYm8sILL7Bo0SJycnLo27cv0dHR+Pn51cyTiIiIiIhIo2dVwHFxcWHBggUsWLCg3HOioqJKjOjc0KxZMz766CM++ugji+93Y/voo0ePWlOm3bh8+TJJSUm2LkNsTJ8D0WdAQJ8DMdPnQKDxfg5uZILKXjFjMJlMprooqCpWr17NxIkTbV2GiIiIiIjUE6tWreKRRx4p93i9DjgZGRls2rSJTp066V05IiIiIiKNWE5ODqdOnWLw4MG0bNmy3PPqdcARERERERGxhsXbRIuIiIiIiNR3CjgiIiIiImI3FHBERERERMRuKOCIiIiIiIjdUMBpgB5//HGMRiNhYWG2LiT6K9QAAAgBSURBVEXq0NatW5k6dSr+/v64u7tz66238vjjj5Oenm7r0qQWXL9+ndmzZ+Pj40OTJk3o378/W7ZssXVZUof27t3L9OnT6dGjBx4eHnTs2JHx48dz4sQJW5cmNjR//nyMRiM9e/a0dSlSx5KSkhg+fDje3t64u7vTs2dPPvjgA1uXVS9pF7UGZt++fQQHB+Po6MgDDzzA+vXrbV2S1JHevXuTmZnJuHHj8PPz44cffmDx4sU0adKEgwcPcsstt9i6RKlBERERfP7550RGRuLn50dUVBR79+7lq6++YsCAAbYuT+rA2LFj2blzJ+PGjaNXr16kpaWxePFirl69yq5du+jRo4etS5Q69tNPP9GtWzeMRiOdO3fm8OHDti5J6sjmzZsJ+//t3V1Ik30fB/CvYFj5ErO09GhLszdnWQmVVhokSUgq1UqiyEoo0bQOAjUhtQiKWB4ZRS9o0cE06iAKVtq0rAiUAlkq0TadFQ1J5/vmdR/c6JPp8zzdN+76z+37AQ/8XdfBV7iuyXfXyz81FevXr4dGo0FAQAA6OjogSRIuXbokOp7bYcGZRSRJQnx8PFavXg29Xg+1Ws2C40UaGxuRkJAwadbQ0IBt27ahqKgIZWVlgpLRTHv37h02btyIK1eu4PTp0wD+vqITHR2N0NBQvHr1SnBCkkNTUxPi4uLg6+s7Mevo6IBarcaePXtQVVUlMB2JsH//fthsNjgcDvz48QMfP34UHYlk0Nvbi6ioKCQkJECn04mOMyvwFrVZpKqqCq2trSgvLwd7qff5vdwAwJYtWxAcHAyj0SggEbmKTqeDr68vsrOzJ2Z+fn44evQompqa0NXVJTAdyWXTpk2Tyg0AREZGYtWqVTznvZDBYEBNTQ20Wi0kSYKPj4/oSCST+/fv4/v377hw4QIAoL+/H2NjY4JTuTcWnFmir68PZ8+eRWFhIW9Fogl2ux19fX3/czVfmn2am5sRFRWFgICASfO4uDgAQEtLi4hY5AYkScK3b994znsZp9OJ3NxcHD9+nLcmeiG9Xo+goCBYLBYsX74cgYGBWLBgAU6ePInh4WHR8dwSC84sUVpaCn9/fxQUFIiOQm5Eq9VidHQUGo1GdBSaQd3d3QgLC5syH59ZrVa5I5GbuHfvHqxWK895L1NZWQmz2cxbkb1Ue3s7HA4H0tLSkJKSgtraWmRlZaGyshJHjhwRHc8t+f7/XWgmSZL0x2177ty5AIC2tjZUVFTgwYMHmDNnjivjkUz+zXHwO4PBgPPnz0Oj0SAxMXEG05Fog4OD8PPzmzIfPxYGBwfljkRuwGg0IicnB5s3b8bhw4dFxyGZ2Gw2lJSUoKSkBAsXLhQdhwSw2+0YGBjAiRMnoNVqAQBpaWkYGRnB9evXUVpaisjISMEp3Quv4Mjs5cuXmD9//h/9tLW1AQBOnTqF+Ph4pKenC05PM+XfHAe/MhqNSE9PR0xMDG7evCngLyBXmjdv3rQFeGhoaGI7eZevX79i165dUCgU0Ol0fP7CixQXF2PRokXIzc0VHYUEGf/MP3DgwKT5+O9v3ryRPZO74xUcma1cuRJ37tz5o32XLFmCFy9e4NmzZ6itrcWXL18mtjkcDgwMDMBkMiE4OBiBgYGuCUwu8U+Pg19ZLBYkJydDoVDgyZMn8Pf3d0FCEiksLGza29C6u7sBAOHh4XJHIoF+/vyJlJQU9Pb2oqGhYcpnAnmu9vZ23LhxA1qtFp2dnRPzoaEhjIyMwGQyISgoCAqFQmBKcrXw8HC0trZOeQY7NDQUANDT0yMilltjwZHZ4sWLcejQoT/e32w2AwAyMjKmbLNarVCpVNBqtcjLy5uxjOR6//Q4GGez2ZCcnIzR0VHU1dXxhRMeKjY2FvX19ejr65v05cXbt28BAGvXrhUVjWQ2NDSE1NRUdHR0QK/XY8WKFaIjkYy6urowNjaGvLy8af/Pq1Qq5Ofn4+rVqwLSkVw2bNgAvV6Pzs5OLFu2bGI+/kVYSEiIqGhui+vguDmLxYLm5uZJM0mSkJ2dDaVSiaKiIkRHR2Pp0qWCEpJc+vv7sX37dnz69Al1dXWIjY0VHYlcZHwdnMuXL+PMmTMA/rMOTkhICF6/fi04IcnB6XQiIyMDT58+xaNHj7Bz507RkUhmNpsNjY2Nk25JlCQJxcXFsNvtuHbtGiIiIvhmNQ/X0tKCdevWITMzE9XV1RPzzMxM1NTUwGQy8crub1hwZimlUomYmBgu9OlF0tLS8PjxY2RlZU15qUBgYCB2794tJhi5hEajwcOHD1FQUICIiAjcvXsX79+/x/Pnz6ddE4k8T35+PioqKpCamoq9e/dO2X7w4EEBqcgdJCYmwmazcaFPL3Ls2DHcunUL+/btw9atW1FfXw+dTofCwkKUl5eLjud2WHBmKZVKBbVazYLjRVQqFcxm87SLvCqVSnz+/FlAKnKV4eFhnDt3DtXV1ejp6cGaNWtQVlaGHTt2iI5GMklKSoLBYJj2nPfx8YHT6RSQitxBUlISbDYbPnz4IDoKycThcODixYu4ffs2rFYrlEolcnJy+IjCf8GCQ0REREREHoOviSYiIiIiIo/BgkNERERERB6DBYeIiIiIiDwGCw4REREREXkMFhwiIiIiIvIYLDhEREREROQxWHCIiIiIiMhjsOAQEREREZHHYMEhIiIiIiKPwYJDREREREQegwWHiIiIiIg8BgsOERERERF5DBYcIiIiIiLyGH8BL/qByLehI1oAAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x000000002F310630>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(-0.1,3.6)"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function average_normals!(elements, time, normal_field=\"normal\")\n",
    "    d = Dict()\n",
    "    # calculate sum of normals connecting to node k\n",
    "    for el in elements\n",
    "        c = get_connectivity(el)\n",
    "        n = interpolate(el[normal_field], time).values\n",
    "        for (ci, ni) in zip(c, n)\n",
    "            d[ci] = haskey(d, ci) ? d[ci] + ni : ni\n",
    "        end\n",
    "    end\n",
    "    # norm\n",
    "    for (ci, ni) in d\n",
    "        d[ci] /= norm(d[ci])\n",
    "    end\n",
    "    # update back to elements\n",
    "    for el in elements\n",
    "        c = get_connectivity(el)\n",
    "        #new_normals = Field(time, [d[ci] for ci in c])\n",
    "        #set_field(el, normal_field, new_normals)\n",
    "        #el[normal_field][end].values = new_normals\n",
    "        #push!(el[normal_field], new_normals)\n",
    "        el[normal_field][end].values = [d[ci] for ci in c]\n",
    "    end\n",
    "end\n",
    "\n",
    "average_normals!(Γ₁, 0.0)\n",
    "average_normals!(Γ₂, 0.0)\n",
    "\n",
    "figure(figsize=(10, 3))\n",
    "for el in Γ₁\n",
    "    plot_element(el; plot_with_normal=true)\n",
    "end\n",
    "for el in Γ₂\n",
    "    plot_element(el; plot_with_normal=false)\n",
    "end\n",
    "axis(\"equal\")\n",
    "ylim(-0.1, 3.6)\n",
    "#axis(\"off\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "calc_projection (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function newton(R, dR, x0=0.0; max_iterations=10, tol=1.0e-9)\n",
    "    x = x0\n",
    "    for i=1:max_iterations\n",
    "        dx = -R(x)/dR(x)\n",
    "        x += dx\n",
    "        if abs(dx) < tol\n",
    "            break\n",
    "        end\n",
    "    end\n",
    "    x\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "Find projection from slave nodes to master element.\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "sel :: Element\n",
    "    slave element\n",
    "mel :: Element\n",
    "    master element\n",
    "sxi :: Vector\n",
    "    projection point in slave side (typically [-1.0] or [1.0])\n",
    "\n",
    "Returns\n",
    "-------\n",
    "mxi :: Vector\n",
    "    point in master element corresponding to xi in slave\n",
    "\n",
    "\"\"\"\n",
    "function calc_projection_from_slave_to_master(sel, mel, sxi; solver_options...)\n",
    "    # slave side point\n",
    "    X1 = interpolate(sel, :Geometry, sxi)\n",
    "    N1 = interpolate(sel, :Normals, sxi)\n",
    "    # master element function and their derivatives\n",
    "    X2(xi) = interpolate(mel, :Geometry, xi)\n",
    "    dX2(xi) = dinterpolate(mel, :Geometry, xi)\n",
    "    # residual & solution\n",
    "    R(xi) = det([X2(xi)-X1 N1])\n",
    "    dR(xi) = det([dX2(xi) N1])\n",
    "    mxi = newton(R, dR; solver_options...)\n",
    "    return [mxi]\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "Find projection from master to slave element.\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "sel :: Element\n",
    "    slave element\n",
    "mel :: Element\n",
    "    master element\n",
    "mxi :: Vector\n",
    "    projection point in master side (typically [-1.0] or [1.0])\n",
    "\n",
    "Returns\n",
    "-------\n",
    "sxi :: Vector\n",
    "    point in slave element corresponding to xi in master\n",
    "\"\"\"\n",
    "function calc_projection_from_master_to_slave(sel, mel, mxi; solver_options...)\n",
    "    # slave element functions and their derivatives\n",
    "    X1(xi) = interpolate(sel, :Geometry, xi)\n",
    "    dX1(xi) = dinterpolate(sel, :Geometry, xi)\n",
    "    N1(xi) = interpolate(sel, :Normals, xi)\n",
    "    dN1(xi) = dinterpolate(sel, :Normals, xi)\n",
    "    # master side point\n",
    "    X2 = interpolate(mel, :Geometry, mxi)\n",
    "    # residual & solution\n",
    "    R(xi) = det([X1(xi)-X2 N1(xi)])\n",
    "    dR(xi) = det([dX1(xi) N1(xi)]) + det([X1(xi)-X2 dN1(xi)])\n",
    "    sxi = newton(R, dR; solver_options...)\n",
    "    return [sxi]\n",
    "end\n",
    "\n",
    "function has_projection(xi1, xi2)\n",
    "    l = abs(xi1[2]-xi1[1])[1]\n",
    "    return l > 1.0e-9\n",
    "end\n",
    "\n",
    "function calc_projection(sel, mel; clamp=true)\n",
    "    xi1a = calc_projection_from_master_to_slave(sel, mel, [-1.0])\n",
    "    xi1b = calc_projection_from_master_to_slave(sel, mel, [ 1.0])\n",
    "    xi2a = calc_projection_from_slave_to_master(sel, mel, [-1.0])\n",
    "    xi2b = calc_projection_from_slave_to_master(sel, mel, [ 1.0])\n",
    "    xi1 = Vector[xi1a, xi1b]\n",
    "    xi2 = Vector[xi2a, xi2b]\n",
    "    if clamp\n",
    "        clamp!(xi1, -1, 1)\n",
    "        clamp!(xi2, -1, 1)\n",
    "    end\n",
    "    return xi1, xi2\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "integrate_segment (generic function with 1 method)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function integrate_segment(sel, mel, xi1, xi2)\n",
    "    nsel = get_number_of_basis_functions(sel)\n",
    "    nmel = get_number_of_basis_functions(mel)\n",
    "    De = zeros(nsel, nsel)\n",
    "    Me = zeros(nsel, nmel)\n",
    "    return De, Me\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzgAAAEeCAYAAABG5rCAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl4VPX99vF7kpkJ2diCyCabCQEFRVEQiaJAA/4gVFxRqiiVAlYfRagookVNJbiLVbAUFLHaYpFqUmRfg7K4pEgChIBhlS0BsyeTmXn+GBMJWckyZ2byfl1XLug5M+fcwUJy53vO55icTqdTAAAAAOAD/IwOAAAAAAD1hYIDAAAAwGdQcAAAAAD4DAoOAAAAAJ9BwQEAAADgMyg4AAAAAHwGBQcAAACAz6DgAAAAAPAZFBwAAAAAPoOCAwAAAMBnUHAAAAAA+AwKDgAAAACfQcEBAAAA4DMoOAAAAAB8BgUHAADAg+3evVsREREKCAiQ1WpVQECAIiIitHv3bqOjAR7J5HQ6nUaHAAAAQHnJycnq3bu3iouLy+0zm83auXOnevToYUAywHNRcAAAADxURESE0tLSKt0fHh6uffv2uTER4PkoOAAAAB4qICBARUVFle63Wq0qLCx0YyLA83EPDgAAgIeq7ufQ/JwaKI+CAwAA4KFMJlOd9gONEQUHAADAAxUUSP7+Hat8TceOVe8HGiMKDgAAgIfJz5dGjpQcji/k72+u9HUPPfSQG1MB3oGCAwAA4EHy8qSYGGnLFunLL3vohx92Kjw8XFarVRaLRVarVZdeeqluu+02TZ8+XZ988onRkQGPUvmPBAAAAOBWubmucrN9u7R8uTRwoCT1qHAUtMPh0Lhx4/S73/1O/v7+uuuuu9yeF/BEFBwAAAAPkJMjDR8uffedtGKFFBVV9ev9/Py0YMEC2e123XvvvfL399ftt9/unrCAB+M5OAAAAAbLzpZuuUXaudNVbq6/vubvtdvtuu+++/Tpp59qyZIlGjVqVMMFBbwABQcAAMBAP//sKjfJydLKldJ11134MYqLizVmzBh99tlnWrp0qUaOHFn/QQEvQcEBAAAwyNmz0rBh0t690qpV0rXX1v5YNptN99xzj7744gstW7ZMw4cPr7+ggBeh4AAAABjgzBkpOlrav19avVrq06fux7TZbLrrrru0fPlyff755xo2bFjdDwp4GQoOAACAm2VkSL/5jXTokLRmjdS7d/0du6ioSHfccYdWrVqlL774QtHR0fV3cMALUHAAAADc6PRpacgQ6ehRae1a6Yor6v8chYWFuv3227V27VolJCRo8ODB9X8SwENRcAAAANzk1Clp8GDp+HFp3TqpZ8+GO1dBQYFGjRqljRs3asOGDerbt2/DnQzwIH5GBwAAAGgMTpyQbr5ZOnlS2rChYcuNJDVp0kTLli3T5MmTdfnllzfsyQAPwgoOAABAAzt+XBo0yDU1bd06qXt3oxMBvstsdAAAAABfduyYq9xkZ7tWbrp1MzoR4NsoOAAAAA3k6FHXZWn5+dLGjVJ4uNGJAN9HwQEAAGgAhw+7yo3N5io3XbsanQhoHCg4AAAA9ezgQVe5cTpdl6V16WJ0IqDxYIoaAABAPfrxR2ngQMlkotwARqDgAAAA1JMDB6SbbpLMZle56dTJ6ERA40PBAQAAqAdpaa6Vm4AA1z03l1xidCKgcaLgAAAA1FFqqqvcBAe7yk379kYnAhovCg4AAEAd7NnjKjfNmrkuS2vb1uhEQONGwQEAAKillBTXPTdhYdL69VKbNkYnAkDBAQAAqIVdu1yjoFu3dpWbiy82OhEAiYIDAABwwXbudJWbdu2kdeukiy4yOhGAEhQcAACAC5CU5Co3HTtKa9dKrVoZnQjAuSg4AAAANfTdd9KgQVLXrtKaNVLLlkYnAnA+Cg4AAEAN7NghDR4sRURIq1dLLVoYnQhARSg4AAAA1di2TfrNb6QePaRVq6TmzY1OBKAyFBwAAIAqfP21q9z07CmtXOl63g0Az0XBAQAAqMSWLVJ0tHTVVdKKFVJoqNGJAFSHggMAAFCBTZukoUOla6+Vli+XQkKMTgSgJig4AAAA51m/XrrlFum666SEBCk42OhEAGqKggMAAHCOtWul4cOlAQOk+HgpKMjoRAAuBAUHAADgF6tWSSNGSAMHSl98IQUGGp0IwIWi4AAAAEj68ktp5EjXs27+8x+pSROjEwGoDQoOAABo9BISpFtvdU1MW7pUCggwOhGA2qLgAACARu3zz6XbbnPdd/Pvf1NuAG9HwQEAAI3WsmXSHXe4Lk37178kq9XoRADqioIDAAAapX//W7rrLtfqzSefSBaL0YkA1AcKDgAAaHT+9S9p9Gjpzjulf/yDcgP4EgoOAABoVD7+WLr3XtfH4sWS2Wx0IgD1iYIDAAAajQ8/lO67z/Xx/vuSv7/RiQDUNwoOAABoFD74QHrgAdfHwoWUG8BXUXAAAIDPW7BAGjdOGj9emj9f8uM7IMBn8dcbAAD4tL/9TXroIWniRGnuXMoN4Ov4Kw4AAHzW3LnShAnSo49K77xDuQEaA/6aAwAAn/TXv0oPPyw99pj01luSyWR0IgDuQMEBAAA+5623XKs2U6ZIb7xBuQEaEwoOAADwKa+9Jj3+uDRtmvTKK5QboLGh4AAAAJ/x8svS1KnSM89Is2ZRboDGiIIDAAB8wksvuVZtnntOevFFyg3QWFFwAACA13vxRdeqzfPPuz4oN0DjRcEBAABeISEhQRaLRSaTqfTDYrFo9OgEPfecFBvrWr0B0LiZnE6n0+gQAAAAVfn888916623Vrr//vvjtWjRCDcmAuCpKDgAAMDjWSwWFRcXV7rfbDbLZrO5MREAT0XBAQAAHs9Ug5tq+JYGgMQ9OAAAAAB8CAUHAAAAgM+g4AAAAI/1c8HPemrNU1KLql9nNpvdEwiAx6PgAAAAj1PsKNbcHXMV/na43t7+tu557p4qX79s2TI3JQPg6Sg4AADAYzidTi3ft1xXzL1Cf1z+Rw2PGK7UR1L18eMfKz4+vtxKTcnwgbVr1xoRF4AHYooaAADwCDtP7NTUVVO1+sBq3dT5Jr0W/Zqubnt1le9xOBzq3LmzDh8+rDVr1mjw4MFuSgvAU1FwAACAoY7nHNez657Vgu8XKLxluF6NflUx3WJqNBpako4cOaKuXbvKbDbr2LFjat68eQMnBuDJuEQNAAAYIs+Wp9hNsQqfE66lu5fqzWFvatfDuzQycmSNy40kdejQQYsXL1Z+fr4GDBjQgIkBeANWcAAAgFs5nA59/MPHmr52uo7nHNejfR/VjBtnqEVgNaPSqjFmzBh9/PHHevTRRzVnzpx6SgvA21BwAACA22w+uFlTVk3RjmM7dFuP2zR7yGyFtwyvl2M7HA517dpVBw8e1IoVKzR06NB6OS4A70LBAQAADW5/5n5NWzNNS3cvVZ+2ffTG0Dd0Q6cb6v08x44dU+fOneXv76+jR4+qZcuW9X4OAJ6Ne3AAAECDOZN/RlNWTlGPd3po29Ft+vDWD7V9/PYGKTeS1K5dO3388ccqKCjgfhygkWIFBwAA1Dub3aa538zV8xufV2FxoZ6KekpP9H9CQZYgt5z//vvv1+LFi/Xwww/rnXfeccs5AXgGCg4AAKg3TqdT8anx+tPqPyktM03jeo/Ti4NeVJuQNm7N4XA4dOmllyo9PV3Lly/XLbfc4tbzAzAOBQcAANSLpONJemLlE1qfvl6DuwzWa9Gv6co2VxqW5/jx4+rUqZP8/Px0+PBhtWrVyrAsANyHe3AAAECdHMs+pnGfj9PV712tn3J+UsI9CVp932pDy40ktWnThvtxgEaIFRwAAFAruUW5evWrV/XyVy8ryBKk5296XuOvHi+Lv8XoaGU8+OCD2rhxo1JSUtSkSROj4wBoYBQcAABwQRxOhz7834d6Zt0zOp13Wo/1e0zTb5iu5k2aGx2tQg6HQ5Lk58eFK0BjQMEBAAA1tiF9g55Y+YS+P/697rr8LsUNjlOXFl2MjgUApcxGBwAAAJ4vNSNVT65+Up/v/Vz92vfTlnFbdP0l1xsdCwDKoeAAAIBKZeRl6IWNL+jdb95Vu9B2+uT2T3T35XfLZDIZHQ0AKsQlagAANEKLFi3SwIED1blz53L70tPTtXb9WmV1z9ILm16Q3WHX9Bum67F+jynQEuj+sABwASg4AAA0Qunp6Ro3bpwWLlxYpuT8+OOPGnnPSGUPy9Zhv8N66KqH9MLNL+jikIuNCwsAF4BxIgAANEKdO3fWwoULNW7cOKWnp0uS4rfF66r/u0q7rtul7uHd9b+J/9N7Me9RbgB4FVZwAABoxNLT03XH3XeoqGORfvjmB0WMi9Db97ytoeFDjY4GALXCkAEAABqpnKIcLTiwQEmOJNn/bdfvY3+veU/Pk9mPbw8AeC9WcAAAaGTsDrs+SPpAM9bPUOaxTLVZ00Y5R3JksVi0devWCgcPAIC34B4cAAAakTUH1ujqv12th+IfUt/gvuqzo482fr5RL730kk6cOKG77rqr9J4cAPBGFBwAABqB3ad2a8THI/Sbxb9RiDVEn0V/puwl2fr4w4/VuXNnjRs3TuHh4QoNDS0zeAAAvA0FBwAAH3Yq95QeWf6Ies3tpZRTKfr0zk+V+GCislKzyoyItlgsio2N1bp16zRp0iRt3LjR2OAAUEvcgwMAgA8qLC7UnG1z9JfNf5FTTs24YYYe7feompibVPoeh8Oha665RkFBQdq8ebNMJpMbEwNA/aDgAADgQ5xOp/6d8m9NWzNNh34+pInXTNSfB/5ZFwVfVKP3r1y5UsOGDVNCQoKGDx/ewGkBoP5RcAAA8BHbjmzTE6ue0FeHv9LwiOF6NfpVdW/V/YKO4XQ6NWjQIGVkZCgpKUl+flzNDsC78K8W4IF2796tiIgIBQQEyGq1KiAgQBEREdq9e7fR0QB4oINnD+repffqugXXKacoR6vvW62EexMuuNxIkslk0qxZs/TDDz/ok08+aYC0ANCwWMEBPExycrJ69+6t4uLicvvMZrN27typHj16GJAMgKfJKszSrM2z9MbWN9QisIX+MugvGnvlWPn7+df52Lfeeqt27typPXv2yGq11kNaAHAPCg7gYSIiIpSWllbp/vDwcO3bt8+NiQB4mmJHsRZ+v1DPrn9W2YXZmtJ/ip4c8KRCA0Lr7RzJycnq1auX3n77bf3xj3+st+MCQEOj4AAeJiAgQEVFRZXut1qtKiwsdGMiAJ5kZdpKTV09VbtO7tJ9V9ynlwa/pA5NOzTIuR544AGtWLFC+/fvV3BwcIOcAwDqG/fgAB6mup858DMJoHHadXKXhn00TMP+MUwtmrTQjvE79OGoDxus3EjSzJkzdebMGb311lsNdg4AqG8UHMDDVPfcCafNqe+ivtPBuIPK2ZVD4QF83Mnck5qYMFFXzrtS+8/s12d3faaND2zUNe2uafBzd+7cWZMmTdLs2bOVkZHR4OcDgPpAwQE8TMeOHavc3+GiDrJeZNXBFw/qm17faGuXrUp9JFUZKzJkL7C7KSWA+hAZGamJEyfKYrHIZDKVflgsFo3/w3hd1PEihc8J17+S/6VXf/Oqkh9O1qgeo9z6AM7p06fL4XBo9uzZbjsnANQF9+AAHmb37t3q1auX7PbyZeXcKWr2ArvObjirjIQMZf43UwXpBfIL8lOLIS0UNiJMYcPDFNAuwIDPAEBN/eEPf9D8+fMrf8E10mMvPqZnb3xWYUFh7gt2npkzZ2r27Nnat2+fOnRouEviAKA+UHAAD/TGG2/oiSeekMVikeS6bK1jx4764osvKhwR7XQ6lZeSp4z/ZigjPkM/f/Wz5JBC+oS4ys6IMIVeHSqTn/t+6gugehaLpcKR8CX8zf4qtlW+312ysrJ06aWXatSoUfrb3/5mdBwAqBIFB/BAY8eO1XfffacffvihVu+3ZdqU+WWma3VnRaaKzxbLcrFFYcNdZafFkBYyh5rrOTWAC1WTS8085cv0m2++qalTpyo5OVmRkZFGxwGASlFwAA9jt9vVpk0bPfTQQ5o1a1adj+ewOZT1VZYyEjKUkZChvD15MllNan5T89LVncAugfWQHEC1MjOlTZukDRukDRtk/d//ZKvmLZ7yZbqgoECRkZHq16+flixZYnQcAKgUBQfwMFu3blX//v2VmJioAQMG1Pvx8/fnu8rOfzN0dsNZOW1OBV0WVFp2mvZvKj8z80eAenFeodHOnZLTKXXpIg0cKP8PPpCjmkN40pfp999/X+PGjdOOHTt0zTUNP8UNAGqDggN4mBkzZmju3Lk6efKk/P39G/RcxdnFOrPqjOvenf9myHbSJnMLs1re0lJhI8LUcmhLWVpaGjQD4FMqKzSdO0s33eT6GDjQ9b9V/T04ZrNZNlt1azzuU1xcrCuuuEIdOnTQqlWrjI4DABWi4AAepnfv3urVq5cWL17s1vM6HU5lf5OtjHjXpWw5STmSv9Ts+mYKi3Gt7gR1D3LreFrA41VXaAYOlG6+WerUqcK3T5w4Ue+9916lh58wYYLmzZvXEMlrbdmyZbrtttu0Zs0aDR482Og4AFAOBQfwIIcPH1bHjh31z3/+U3fffbehWQqPFpau7JxZfUaOfIeadG1Seilb8xubyy+AS9nQyFRzyZluvtn1ayWF5nyRkZG6+eabtWDBgnIrOQ888IC++uor7d27t/4/jzpwOp267rrrJLkuqeWHHgA8DQUH8CDz5s3TI488otOnT6t58+ZGxyllz7fr7LqzrsKTkKHCw4XyD/FXi9+4nrnT8v9aKqANz9yBD7rAS87qQ9++fbVjxw6Fh4dr37599Xbc+rR+/XoNGjRIn332mUaNGmV0HAAog4IDeJARI0YoNzdX69evNzpKpZxOp3J/yC2dypa1NUtySqHXhrpWd2LCFNI7hJ/qwjvV8ZKz+lBQUKDAQNdkQ0/+Ej106FAdPnxYO3fulNnM2HkAnoOCA3iIvLw8hYWFKTY2VlOmTDE6To0VnSr69Zk7KzNlz7LL2s766zN3BreQf3DDDksAzjV+/HiNHTtWUVFR5fYlJiZq0aJFmj9/vmtDPV9yVl+aNGmiwsJCvfHGG3r88cfdeu6a+vbbb3XNNddo4cKFevDBB42OAwClKDiAh0hISFBMTIz27NnjtQ/RcxQ59HPiz6WrO/n78mUKMKnFINelbGHDw9SkUxOjY8LHJSYmKiYmRvHx8WVKTmJiomJGjFD8k08q6uRJt11yVhuLFy/W/fffL4vFoqKiIkOzVOXuu+/W119/rdTUVDVpwt9tAJ6BggN4iIkTJ2rt2rVKTU31mcu78lLzSp+58/Omn+Usdiq4V/Cvz9zp11Qmf9/4XOFZSkrO3FdeUdiRIwpMTlbMsmWKt9sVJbn1krPaKvl3IDc3V0FBQQanqVhqaqouu+wyvfLKK5o8ebLRcQBAEgUH8AhOp1OXXHKJ7rzzTr3xxhtGx2kQxT8XK3PVL5eyLc+U7bRN5jCzwv7vl0vZolvI0pxn7qCOzrnkLDEhQTfu369ASVY/P8VHRyvqnnsMueSsNq666iolJSWpR48eSklJMTpOpSZMmKClS5fqwIEDatq0qdFxAICCA3iCpKQkXXXVVY3muRJOu1NZ27NKn7mT+0OuTGaTmkU1K13dCYr0zJ9Yw8NUNRTg5pvVY8UK7fnpJ7388sv605/+ZHDYC+MtwwaOHj2q8PBwPfnkk3r++eeNjgMAFBzAE8TGxurll1/W6dOnZbVajY7jdgWHCkqfuXN27Vk5ChwKjAgsLTvNoprJz8ozd6ALmnKWmJiooUOHKi8vT0FBQVq5cmWFgwc8WUBAgIqKijRv3jxNmDDB6DiVevLJJ/Xuu+/qwIEDat26tdFxADRyFBzAA3z22Wfas2ePpk+fbnQUw9nz7Dqz5kzpM3eKjhXJP9RfLYe2dD1z55aWsrYuWwIXLVqkgQMHqnMFN4anp6dr48aNGjt2rJs+A9SrWk45K7kHZ8mSJYqOjtZ1112nPXv2lBs84On+/ve/a/z48bJarSosLDQ6TqUyMzPVtWtXjR07Vm+99ZbRcQA0chQcAB7L6XQqJymndCpb9o5sSVLTfk1Ln7kT3CtYBw8e1Lhx47Rw4cIyJSc9Pb3C7fBglRWaTp1+LTM33VTllLPzp6iFhYXJbDZr6dKlFU5X83TeMGxAkl566SXNnDlTqamp/H0DYCgKDgCvUXSiSBnLXZeynVl5RvYcuwIuCVDY8DDlXpurJz58Qu9/8L46d+6s0aNHa8eOHTp48KDsdnvpMcxms2JjY5WWlvbrs1AayAU9j6WxqsklZ9UUmvOd/+c+ZMgQrV27VoWFhdq+fbvX/bn37NlTycnJ6tWrl3bu3Gl0nErl5ubq0ksv1dChQ7Vo0SKj4wBoxCg4ALySo9Chs5vOugYV/DdDBQcKdCLghKb4TdHvon+nJd8v0Y+Hfqz0/XFxcZo2bVqDZqzyeSwevJLQoMXMgAdr/vWvf9Wjjz6qjz76SGPGjKm347pLXl6egoODJXn2sAFJevfdd/XII49o586d6tmzp9FxADRSFBwAXs/pdCpvb55Of3Fal02/THn2PIUpTBnKqPQ9ZrNZNputwbOdX2Y8vdxI1RezqKgoRUVFacaMGSouLi7dX+HqWDVTzjRwYIM/WDMrK0vNmjXTrbfeqmXLljXYeRpSybCB999/Xw888IDRcSpVVFSkHj16qFevXvrPf/5jdBwAjRQFB4DP2L59u/r166eZ02dq5kszq319r+3bZTGZZPXzk8VkKv299fzfn7Ot5Pcl77Ge97qSbee+bs/WrZo+erSihg7VV6tW6a+ffqr+UVGVns/fZDL8Ya9VFbNNmzbpmWeeqfS9cb/7naaFhVU75cydmjdvrqCgIB07dsyt560vc+bM0WOPPabWrVvrxIkTRsep0scff6wxY8boq6++Uv/+/Y2OA6ARouAA8Anp6ekaNWqUkpKStHfvXkVGRlb7nkdTU2VzOlXkcJT+WuR0lttmczpVVMk22y/vKXI4ZK/qZI895vpm//77pQcfrDbbueWpoiJUrpBVsP9CStj557OYTNq7bZueGT1a/YYM0fa1a0uL2eUtW6q4ioleZkm2BrzkrDYGDhyozZs3q6ioSGaz2dAsteXv7y+HwyG73S4/P88dm+5wOHTVVVepefPm2rBhg+FlHUDj453/ygPAOUqmpf32t79VSkqKunbtWqP3zYmIqNccjl/Kke28ovTVli26b/du2SSFfP65Xhk9Wpf171/j8lRyzMIKjn3+67LtdtmKiy+orFVazIKCpJ49tW7ZMunWW/WA1Spt3y7Zq6xyKpakAwfq9c+2rkaOHKlNmzbp888/1+233250nFoZPHiwVq9erfvvv18fffSR0XEq5efnp1mzZmn48OFauXKlhg0bZnQkAI0MKzgAvF7Jc3Di4uL09ddf63//+58sFkuZ+0PO5+57cHr27KnExERt3rzZ4+7BcTidKv6lDBWeU3y2bNmicSNHKj8vTwGBgXpz6VJd1r+/BrZoUe0xPe1Ly+nTp3XRRRfpjjvu0Keffmp0nFopGTbgrv/v1oXT6dTAgQOVnZ2tb7/91qNXnAD4Hv7FAeD1xo4dq86dOyslJUWXXXaZJCk2NrbK91S3vz6ce9+K2WyWn5+foqKiFB8fr5iYGCUmJjZ4hprw++XStGB/f7W0WHSx1aqD33yjSXfcoXlz50qS7rjtNj19773y27XL4LS106pVK4WGhmrr1q1GR6m1oKAgtWjRQsXFxR7z/53KmEwmzZo1S0lJSVqyZInRcQA0MhQcAD7B6XQqOTm5tOCkpaUpLi6u3PX/ZrNZcXFxSktLa/BMixYtKl2pyc7Olr+/vySVlhxPfVbIucVs9OjRklx/viXFrLqfxnvqPS49e/bU0aNH5XA4jI5Sa2+//bYk6d577zU4SfUGDBigESNGaMaMGR6/4gTAt1BwAPiEkydPKjMzU5dffrkkaf78+Zo2bZpiYmIkSS+88IKcTqdsNpumTZvmlgc9zp8/v/QytNzcXFksltJ9UVFRHvuwyXOLmdVqlZ+fn3766afSYtanT58q3++O1bHaGD58uJxOp1asWGF0lFobM2aM/Pz8dPjwYa8oan/5y1904MABLViwwOgoABoRCg4An5CSkiJJpSs4JZKSkiRJkyZNcnumc+Xm5spqtRqaoabOLWaSZLVaderUKUmuYnbllVcqLi6u3EqNO1fHauPBX6bX/eMf/zA4Sd3ceOONklwPZPV0V1xxhcaMGaPnn39eeXl5RscB0EgwZACAT3jnnXc0efJk5eXllfnGOzAwUAUFBYbf9N66dWv5+/vrp59+MjRHbYSFhSkwMFBHjhypcF9mZqbhf741FRISolatWik9Pd3oKLVW8uBSi8WioqIio+NU68CBA+revbteeOEFPfXUU0bHAdAIsIIDwCckJycrMjKy3KpCYWGhR0xwKiwsVFBQkNExaqVZs2bKzs6ucF+bNm0kuf78vUGPHj285vKuyjRt2lTNmzeXzWbTjh07jI5Tra5du2rChAmaPXu2zpw5Y3QcAI2A8V/1AaAenDtB7VxOp9MjikVRUZFH5KiNli1bKj8/v8J93bp1kyR99tln7oxUa7fccoscDoc2bNhgdJQ6ee211yRJd911l8FJaqZk0MDs2bONjgKgEaDgAPAJycnJpQMGzt0mSe3btzciUhnFxcUKCQkxOkatXHzxxZVOwSq5V+frr792Z6Ra+/3vfy9JWrx4scFJ6mbcuHEymUxKT0/3itWoiy++WJMnT9acOXN07Ngxo+MA8HEUHABe79SpUzp9+nS5FZx58+ZJkvr3729ErDLsdrtCQ0ONjlErJQXx7Nmz5fbdfffdkuSxgwXO16lTJwUGBmrTpk1GR6mz66+/XpL0xz/+0eAkNTN16lQFBgbqxRdfNDoKAB9HwQHg9SqboLZx40ZJ0kMPPeT2TOdzOp1q1qyZ0TFqpVOnTpIqvs+mQ4cOklQ6Zc0bREZG6uDBg0bHqLOEhARJ0sKFCw1OUjPNmjXT9OnTNX/+fO3bt8/oOAB8GAUHgNe0b8NdAAAXaklEQVRLTk6W2WxWREREme0l38QOGDDAiFilsrKyJEktWrQwNEdthYeHS5JSU1MrfU1ubq674tRZdHS07Ha7tmzZYnSUOmnevLmaNm2qoqIi7dy50+g4NfLwww+rbdu2eu6554yOAsCHUXAAeL2UlBR169atzIM0Jdc33SaTyaBUvyoZDe2tBScyMlKSa9xvRfz8/FRcXOzOSHUybtw4Sd5/H44kxcXFSZJGjRplcJKaCQwM1MyZM/XPf/5T33//vdFxAPgoCg4Ar1fZBDW73e4RD9csKTitWrUyOEntdO/eXZJ06NChCvdbrVaveQ6O5CpsAQEBWr9+vdFR6mzSpEkymUw6cOCAVwwbkKSxY8cqMjJS06dPNzoKAB9FwQHg9SqaoJaTkyPJ9SBKo504cUKS62Gf3qhJkyYymUyVTr/yxulwERER+vHHH42OUS/69u0rSZo8ebLBSWrGbDYrNjZWK1as8Ppx3QA8EwUHgFc7ffq0Tp48WW4FZ8GCBZKknj17GhGrjJIb8C+66CKDk9Se1WrVyZMnK9xX8rDP7777zp2R6mTIkCGy2WxelbkyJcMG3nvvPYOT1Nztt9+uPn366Omnn/aq1T8A3oGCA8Cr7d69W1L5CWr/+c9/JLm+kTLa6dOnJf1aBLxRUFBQpU+hL7mEreTP3Bs8+OCDkqQPPvjA2CD1oFWrVgoJCVFhYWHpREFPZzKZFBcXp61bt+qLL74wOg4AH0PBAeDVkpOT5e/vr27dupXZXrJqcv/99xsRq4zMzExJUrt27QxOUntNmzYtnQZ3vptuukmStG3bNjcmqpsrrrhCVqtVa9euNTpKvSh5toy3DBuQXKtogwcP1vTp02W3242OA8CHUHAAeLWUlBRFRESUGyZgt9sVGBioJk2aGJTsVyUFx1vvwZFc9zLl5+dXuK9klWz//v3ujFRnXbt29ZoHlFbn8ccfl8lkqnKUtyeaNWuWUlJS9NFHHxkdBYAPoeAA8GqVTVA7duyYxxSKs2fPSnKNU/ZWrVu3ls1mq3BfyaV33jQqWpJuvvlmFRUVVfgAU2909dVXS5KmTZtmcJKau/baa3X77bfrueeeU2FhodFxAPgI7/1qCwCqeIKaJGVnZ6tr164GJCovOzvbq8uNJLVv315Op7N0Ot35goODFRAQ4OZUdeNL9+FIvw4bmDNnjsFJLkxsbKyOHDniVUMSAHg27/6KC6BRy8zM1PHjx8ut4Ozbt09Op1NXXnmlQcnKys7OltlsNjpGnXTs2FHSr0Mdzte0aVNlZGS4M1KdjB8/XuvWrZMkvfrqqzKZTDKZTLJYLJo9e7bGjx9vcMIL16ZNGwUHB6ugoMCrLlXr3r27HnzwQcXGxio7O9voOAB8AAUHgNcq+Wb7/BWckm9cBwwY4PZMFcnJyZHFYjE6Rp1ceumlkqS9e/dWuL9169aVru54oi5duuipp54qt724uFhPPfWUwsPDDUhVd88++6wk6dZbbzU4yYX585//rKysLL3xxhtGRwHgAyg4ALxWcnKy/Pz8yk1Q27p1qyRp0KBBRsQqJz8/3+su3zpfZGSkJFV6U36HDh286h6KP//5z1XunzFjhpuS1J/IyMjSh5fu3r27zKrUxIkTS/8beqJLLrlEjzzyiF599dXSCYgAUFvefc0EgEYtJSVF4eHh5cpDcnKyzGazWrZsaVCysvLz8z1imltdlFwGePDgwQr3l9zvdOjQodLL2TxJfr508KCUnu76KC6u+uGSxcXFuu222xQcHKzg4GAFBQWV+bWmv/f393fL5ydJAwcOrPA+luLiYr333nuaMGGC27LUxtNPP6358+dr1qxZev31142OA8CLUXAAeK3KBgwcOnRILVq0MCBRxYqKijwqT20EBQXJZDLp2LFjFe7v0aOHJOn77783pOCcX2DO/zhx4tfXujqHqQbHzFdGRoZyc3NLP/Ly8pSbm6uCgoIa5bJarVWWpNqUpnOPc+69Xe+//36VWRYsWKB58+bVKLcRwsLCNHXqVMXGxurxxx/3yKIMwDtQcAB4rZSUlNJJWOfKzMxUr169DEhUsaKiIgUFBRkdo84sFotOnjxZZtv48eM1duxY9e7dW5L0ww8/6Le//a0kKTExUYsWLdL8+fPrfO4LLTCXXCJ17iz16CHdcovr9yUf7dtLFkv1I62//PLLSvfZ7Xbl5eWVFp5zy0/J7yvad+623NzceilQJYWnujHd3jDGe/LkyXr77bc1c+ZMLVy40Og4ALwUBQeAVzp79qyOHTtWboJaVlaWbDZb6YqCJyguLlZoaKjRMWpl/PjxCg8P14wZM1RcXKykpCSZTCaZzWbFxsbq+PHjiomJ0ZIlSySpdHpXYmKiYmJiFB8fX6Pz1LbAdO8uDRvm+n2XLr8WmOqG1pnN5iq/4a9u6p2/v79CQ0Mb7L+r3W5Xfn5+uUJUVaF64YUXqj3unj0PKjS0r5o27avg4F7y87NW+x53CgkJ0bPPPqvHH39cU6dOrfAZVwBQHQoOAK+UkpIiqfwEtQ0bNkiS+vXr5+5IlXI4HF5bcKqbNhYXF6dp06YpJiZGJpNJBw8eLFNuoqKiJNV9BaZTJ1eB6dJFateu+gJTndjY2Ao/r3P3G8nf318hISEKCQmp8XtqUnBycnbqxImP5HQWy2QKUGjoVaWFJzS0rwIDw2UyVX/5XkP6wx/+oNdff10zZszQZ599ZmgWAN7J5HQ6q77TEgA80N///ndNmDBBOTk5CgwMLN3+5JNP6pVXXtGuXbsqvD/H3RwOh/z9/TV69Gh98sknRse5YBaLpdqVDpvNprVrEzVkyA1q2rSVbLZijRoVL7s9qtoC06WLq7x07vzrrx061L3AVOf8lalzP5/Y2FilpaXVy6V17lTT/1Z2e75ycpKUnb1dWVnblZ29Xfn5ab+8pvk5hedahYb2VUBAG3d9CqUWL16s+++/X1u3bvWoH1YA8A4UHABe6YknnlB8fLz27dtXZnt0dLTWrFmj4uJi+fkZPwk/MzNTYWFhmjRpkt59912j41ywmvw0/+KLnb8UmOaSfpbJ9KY6dnys9JKx8z9qcgkZLtzEiRMrnKJWYsKECZUOGbDZMpWdvaO08GRlbZfN5rrfKiDgkjKrPKGhfWQ2N+yKpN1uV+/evdWqVSutW7fO8FUlAN6FLzEAvFJlE9T279+v4OBgjyg3kkqnjnn7FLWqTJokFRYm6s037br33slaunSmPvqoT+nlaXCP9evXa8KECVqwYEG5Vanf//73Wr9+faXvtVhaqmXLoWrZcqgkyel0qrDwcJnCk57+ghyOXEkmBQX1+KXw9FPTptfW+/08/v7+eumllzRy5EitXr1a0dHR9XZsAL6PggPAK6WkpOi+++4rt/3EiRNq166dAYkq9tNPP0mSLrroIoOT1IzNbtOWw1uUkJqg+NR4yV+Sver3DB7suudm1aovFRUVpQceuK3cPThoeHv37pWkehkFbTKZ1KRJRzVp0lGtW98hSXI67crNTSmz0tOQ9/OMGDFC119/vaZPn64hQ4Z4zA8tAHg+Cg4Ar5OVlaUjR46Um7DkcDiUm5ur8PBwg5KVVzJW2ZMLTmZ+plakrVB8arxWpK3Q2YKzahPSRiMiRijNkiaH3VHpe/38/MqVmaioKMXHx1NyfIzJ5K+QkF4KCemltm3HSdIv9/N8X1p4MjOX6+jROZIks7mFQkOvLS08TZv2ldV68QWcz6S4uDjdeOONWrp0qe68884G+bwA+B4KDgCvU9kEtaSkJEnSVVdd5fZMlTl16pQkqXXr1gYn+ZXT6dSe03uUkJqghH0J2nJoi+xOu65ue7X+X9//pxHdRqhPuz7yM/kpfGZ4ldPG+vTpo9dff71ciSkpOYsWLaLg+DB//0A1a3a9mjW7vnTb+ffzHDv2N9lsrql0F3o/z4EDB9S3b1+NGTNGY8aMkeQqPh07dtTcuXN19OhRjR07tmE/SQBeh4IDwOukpKTIZDIpMjKyzPaSewxuvPFGI2JVKCMjQ5LUtm1bQ3MU2Yu0+eDm0kvP9p/Zr0BzoAZ3Hax3h7+r4RHD1b5p+3LvS0tLU1xcXJXTxiorMFFRUZSbRqjy+3m2VXI/z2VlVnlc9/NYJLn+3mzfvr3cOdLS0hQdHa1Vq1a581MD4CUoOAC8TkpKirp06aKgoKAy23fs2CFJuuGGG4yIVaGSgtOmjftH7Z7OO60v932p+NR4rdy/UlmFWWof2l4juo3Qm93e1KAugxRkCaryGCWjkqdNm+aOyPBBZe/ncV1m9uv9PNuVlbVD2dnbdPz4h5Ls8vNropCQ3goN7asJE/5Z6XGdTqcmTZpUbpIiAFBwAHidyiao7d69W1artVzxMdKZM2ckSS1btmzwczmdTqWcSlF8arwSUhP09ZGv5XA6dG27azWl/xTFdItR7za9GbkLw5W9n+f3kiS7Pe+X+3l2KDt7hzIy/qsjR05WeZxDhw65Iy4AL0PBAeB1UlJSdM8995TbfuTIEYWFhRmQqHI///yzTCZTg02AKiwu1MaDGxW/N14J+xKUfjZdwZZgDek6RH8b8TcN7zZcbULcv3oEXCh//yA1azZAzZoNKN1mMlkl2Sp9D4/yA1ARCg4Ar5Kdna1Dhw5VuILz888/q2/fvgakqlxWVla9l5uTuSe1fN9yxafGa9X+VcopylHHZh01ImKERnQboZu73Kwm5ib1ek7ACNWtNrIaCaAiFBwAXmX37t2SVG5E9PHjx2W329WzZ08jYlUqOztbZnPd/ql1Op364eQPit8br/jUeG0/6rrpul+HfnpqwFOKiYxRr9a9+GYPPqdjx45KS0urcj8AnI+CA8CrlIyI7t69e5nta9eulST179/f7ZmqkpeXJ6v1wp/wXlBcoPU/ri+9n+Zw1mGFWEM09NKhmnTNJN0ScYtaB3vO6GmgIcydO1fR0dEVXopmMpk0d+5cA1IB8HQUHABeJTk5WV26dFFwcHCZ7V999ZUkafDgwUbEqtSFFJyfsn8qvfRs9YHVyrPlqUvzLrq1+62K6RajGzvdqABzQAMnBjzH0aNHtWrVKk2aNEmHDh2S0+ks9xwcADifyckdegC8yPDhw2UymZSQkFBm+4033qgtW7bIbrcblOxX48ePV3h4eJXPjpk/f76cTqe+P/596bNpvjn2jfxMfurfob9GdBuhmG4xuuyiy7j0DACAC0DBAeBVunTpojvvvFMvv/xyme2XXHKJcnJySscyG+mll17SM888U+n+sVPHKuDGACXsS9Cx7GNqGtBUw8KHaUTECN0ScYtaBbVyY1oAAHwLBQeA18jJyVFoaKg++OADjR07tsy+Jk2aKDw8XLt27TIo3a8sFkuZlZtyTFL4W+EaETFCMZExuqHjDbL4W9wXEAAAH8Y9OAC8xp49eySVn6BWVFSkwsJCRUZGGhGrnCrLjSQ5pdRHUrn0DACABtAwT54DgAZQMkGtR48eZbYnJiZKkq655hq3ZyrH4VBNagvlBgCAhkHBAeA1kpOT1alTJ4WEhJTZvmnTJknSTTfdZECqX9jt0j//KV1xhbjYDAAA41BwAHiNlJSUcpenSdK3334rSbr22mvdHUkqLpY++kjq2VO65x6pQwdZ/P2rfEtdH/wJAAAqR8EB4DUqKzj79u1TYGCge4uDzSa9/77Uo4d0331SRIS0bZu0YoWe/ctfqnxrbGysm0ICAND4UHAAeIW8vDz9+OOPuvzyy8vt++mnn9S6dWv3BCkqkubPl7p1k8aNk3r1kr79VvriC6lvX0lSWlqa4uLiyhUus9msuLg4paWluScrAACNENdJAPAKe/bskdPprHAFJzs7W3369GnYAAUF0sKFUlycdOSIdMcdrlLTq1e5l86fP1+SNG3atIbNBAAAymEFB4BXqGyC2t69e+V0OnXllVc2zInz8qS33pK6dpUefVS64QZp1y5pyZIKyw0AADAWKzgAvEJycrIuueQSNW3atMz2devWSZIGDBhQvyfMyZHmzZNeeUXKyJB+9ztp+nTXpWkAAMBjUXAAeIXKBgxs27ZNkjRo0KD6OVF2tvTOO9Jrr0lnz0pjx7qKTdeu9XN8AADQoCg4ALxCSkqKYmJiym1PTk6W2WxWy5Yt63aCs2elOXOkN9+UcnNdAwSeekrq1KluxwUAAG5FwQHg8fLz87V///4KJ6gdOnRILVq0qP3BMzNdpeatt1wT0saPl558UurQoQ6JAQCAUSg4ADxeySCBii5RO3PmjK644ooLP+ipU9Lrr0t//atkt0uTJklTp0pt29ZDYgAAYBQKDgCPVzJB7fyCk5WVJZvNVmHxqdSJE9Krr0rvviv5+UkPPyxNmSK56zk6AACgQVFwAHi85ORktW/fXs2aNSuzff369ZKkvr88YLNKR4+6JqK9955ktUqTJ7s+wsIaIjIAADAIBQeAx6tsglpiYqIkafDgwZW/+fBhafZsaf58KShImjZNeuwxqS737QAAAI9FwQHg8ZKTkzV8+PBy25OSkmQymRQZGVn+TT/+KMXFSe+/LzVtKj33nPTII9J5q0AAAMC3+BkdAACqUlBQoP3795eu4IwfP7505Wb//v0KCQmRn5/rn7LExESNv/tu14jniAhp2TIpNlZKT5eeeYZyAwBAI8AKDgCPlpqaKofDUVpwxo4dq5iYGMXHx+vEiRNq3769JCnxH/9QzIMPKt5mk9q0cd1vM2GC67I0AADQaFBwAHi05ORkSb9OUIuKilJ8fLxiYmKUl5eniDZtlDh4sGLWrVN8q1aKeu456aGHpMBAI2MDAACDcIkaAI+WkpKitm3blnmYZ1RUlF5+9FFJUvHmzYpZv17xU6Yo6sgR6dFHKTcAADRirOAA8DiLFi1S+/btNWnSJO3fv1+SFBAQoI4dO2ru44/r6N//rtZJSbJIWiXpzVdeUdSUKYZmBgAAnoEVHAAep23btoqOjlZaWpqcTqecTqeKioqUlpam6EceUftTpxT29NMKbt5cb775pmbGxpYOHgAAAI2byel0Oo0OAQDnioiIUFpaWqX727drp9y8PMXHxysqKkqJiYmlgweioqLcmBQAAHgaCg4AjxMQEKCioqIqX7N58+YyZYaSAwAAJAoOAA9ktVpls9kq3e/v76/i4uJy2xMTE7Vo0SLNnz+/IeMBAAAPRsEB4HGqW8GxWq0qLCx0YyIAAOAtGDIAwON07NixTvsBAEDjRcEB4HHmzp0rk8lU4T6TyaS5c+e6OREAAPAWFBwAHufo0aNatWqVwsPDZbVaZbFYZLVaFR4erlWrVuno0aNGRwQAAB6Ke3AAAAAA+AxWcAAAAAD4DAoOAAAAAJ9BwQEAAADgMyg4AAAAAHwGBQcAAACAz6DgAAAAAPAZFBwAAAAAPoOCAwAAAMBnUHAAAAAA+AwKDgAAAACfQcEBAAAA4DMoOAAAAAB8BgUHAAAAgM+g4AAAAADwGRQcAAAAAD6DggMAAADAZ1BwAAAAAPgMCg4AAAAAn0HBAQAAAOAzKDgAAAAAfAYFBwAAAIDPoOAAAAAA8BkUHAAAAAA+g4IDAAAAwGdQcAAAAAD4DAoOAAAAAJ9BwQEAAADgMyg4AAAAAHwGBQcAAACAz6DgAAAAAPAZFBwAAAAAPoOCAwAAAMBnUHAAAAAA+AwKDgAAAACfQcEBAAAA4DMoOAAAAAB8BgUHAAAAgM+g4AAAAADwGRQcAAAAAD6DggMAAADAZ1BwAAAAAPgMCg4AAAAAn0HBAQAAAOAzKDgAAAAAfAYFBwAAAIDPoOAAAAAA8BkUHAAAAAA+g4IDAAAAwGdQcAAAAAD4DAoOAAAAAJ9BwQEAAADgMyg4AAAAAHwGBQcAAACAz6DgAAAAAPAZFBwAAAAAPoOCAwAAAMBnUHAAAAAA+AwKDgAAAACfQcEBAAAA4DMoOAAAAAB8BgUHAAAAgM+g4AAAAADwGRQcAAAAAD6DggMAAADAZ1BwAAAAAPgMCg4AAAAAn/H/AZbnZlLoDlpFAAAAAElFTkSuQmCC",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7faed5beb5d0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(-3.0,5.0,-0.1,3.6)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function plot_all()\n",
    "    figure(figsize=(10, 3))\n",
    "    for el in Γ₁\n",
    "        plot_element(el; plot_with_normal=true)\n",
    "    end\n",
    "    for el in Γ₂\n",
    "        plot_element(el; plot_with_normal=false)\n",
    "    end\n",
    "    for sel in Γ₁\n",
    "        for mel in Γ₂\n",
    "            xi1, xi2 = calc_projection(sel, mel)\n",
    "            if has_projection(xi1, xi2)\n",
    "                X1 = interpolate(sel, :Geometry, xi1)\n",
    "                X2 = interpolate(mel, :Geometry, xi2)\n",
    "                for s=1:2\n",
    "                    x = [X1[s][1], X2[s][1]]\n",
    "                    y = [X1[s][2], X2[s][2]]\n",
    "                    plot(x, y, \"-kx\")\n",
    "                end\n",
    "                De, Me = integrate_segment(sel, mel, xi1, xi2)\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    axis(\"equal\")\n",
    "    ylim(-0.1, 3.6)\n",
    "    axis(\"off\")\n",
    "end\n",
    "\n",
    "plot_all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A little experience, raise the degree of elements, smooth geometry, and test projection for 2nd order:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "LoadError: BoundsError: attempt to access (2,)\n  at index [2]\nwhile loading In[9], in expression starting on line 13",
     "output_type": "error",
     "traceback": [
      "LoadError: BoundsError: attempt to access (2,)\n  at index [2]\nwhile loading In[9], in expression starting on line 13",
      "",
      " in getindex at tuple.jl:8",
      " in get_detJ at /home/jukka/.julia/v0.5/JuliaFEM/src/equations.jl:50",
      " in fit_derivative_field! at /home/jukka/.julia/v0.5/JuliaFEM/src/elements.jl:434",
      " [inlined code] from In[9]:14",
      " in anonymous at no file:0"
     ]
    }
   ],
   "source": [
    "# to raise degree of elements and smooth geometry\n",
    "function tangent(el, xi)\n",
    "    normal = interpolate(el, :Normals, xi)\n",
    "    [0 -1; 1 0]'*normal\n",
    "end\n",
    "for el in [Γ₁; Γ₂]\n",
    "    set_degree(el, 2)\n",
    "    for field in (:Geometry, :Normals)\n",
    "        push!(el.fields[field], [0.0, 0.0])\n",
    "        el.fields[field] = el.fields[field][1:get_number_of_basis_functions(el)]\n",
    "    end\n",
    "end\n",
    "for el in [Γ₁; Γ₂]\n",
    "    fit_derivative_field!(el, :Geometry, tangent, Int[1, 2])\n",
    "end\n",
    "plot_all()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.0",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}