From f1c7d7b64473345e56debb55bb3ebdb12637f3bb Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Tue, 15 Sep 2015 22:16:00 +0300 Subject: [PATCH] 2d contact segmentation --- notebooks/2015-09-15-2d-segmentation.ipynb | 323 +++++++++++++++++++++ src/elements.jl | 74 +++++ 2 files changed, 397 insertions(+) create mode 100644 notebooks/2015-09-15-2d-segmentation.ipynb diff --git a/notebooks/2015-09-15-2d-segmentation.ipynb b/notebooks/2015-09-15-2d-segmentation.ipynb new file mode 100644 index 0000000..3366f6d --- /dev/null +++ b/notebooks/2015-09-15-2d-segmentation.ipynb @@ -0,0 +1,323 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Contact segmentation in 2d\n", + "\n", + "Author(s): Jukka Aho\n", + "\n", + "**Abstract**: Contact segmentation in 2d" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING: could not import Base.help into PyCall\n" + ] + } + ], + "source": [ + "using JuliaFEM: PSeg, get_field, set_field, interpolate, new_field!, push_field!\n", + "using JuliaFEM: calculate_normals!, average_normals!, fit_derivative_field!\n", + "using JuliaFEM: set_degree, get_number_of_basis_functions, dinterpolate\n", + "using PyPlot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create some test boundaries" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "srand(42)\n", + "\n", + "nsl = 5\n", + "nm = 5\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", + "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 = PSeg(con)\n", + " set_field(el, :Geometry, Vector[X[:, i], X[:, i+1]])\n", + " push!(Γ, el)\n", + " end\n", + " return Γ\n", + "end\n", + "Γ₁ = create_elements([x1 y1]', 0)\n", + "Γ₂ = create_elements([x2 y2]', nsl);\n", + "\n", + "# calculate and average normals like se did in last notebook\n", + "for el in [Γ₁; Γ₂]\n", + " calculate_normals!(el)\n", + "end\n", + "average_normals!(Γ₁)\n", + "average_normals!(Γ₂)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "calc_projection (generic function with 1 method)" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\"\"\" Find projection from slave nodes to master element. \"\"\"\n", + "function calc_projection_slave_nodes_to_master_element(sel, mel)\n", + " X1 = get_field(sel, :Geometry)\n", + " N1 = get_field(sel, :Normals)\n", + " X2(xi) = interpolate(mel, :Geometry, xi)\n", + " dX2(xi) = dinterpolate(mel, :Geometry, xi)\n", + " R(xi, k) = det([X2(xi) - X1[k] N1[k]]')\n", + " dR(xi, k) = det([dX2(xi) N1[k]]')\n", + " xi2 = Vector[[0.0], [0.0]]\n", + " for k=1:2\n", + " xi = xi2[k]\n", + " for i=1:3\n", + " dxi = -R(xi, k)/dR(xi, k)\n", + " xi += dxi\n", + " if abs(dxi) < 1.0e-9\n", + " break\n", + " end\n", + " end\n", + " xi2[k] = xi\n", + " end\n", + " clamp!(xi2, -1, 1)\n", + " return xi2\n", + "end\n", + "\n", + "\"\"\" Find projection from master nodes to slave element. \"\"\"\n", + "function calc_projection_master_nodes_to_slave_element(sel, mel)\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", + " X2 = get_field(mel, :Geometry)\n", + " R(xi, k) = det([X1(xi) - X2[k] N1(xi)]')\n", + " dR(xi, k) = det([dX1(xi) N1(xi)]') + det([X1(xi) - X2[k] dN1(xi)]')\n", + " xi1 = Vector[[0.0], [0.0]]\n", + " for k=1:2\n", + " xi = xi1[k]\n", + " for i=1:3\n", + " dxi = -R(xi, k)/dR(xi, k)\n", + " xi += dxi\n", + " if abs(dxi) < 1.0e-9\n", + " break\n", + " end\n", + " end\n", + " xi1[k] = xi\n", + " end\n", + " clamp!(xi1, -1, 1)\n", + " return xi1\n", + "end\n", + "\n", + "function has_projection(sel, mel)\n", + " xi1 = calc_projection_master_nodes_to_slave_element(sel, mel)\n", + " l = abs(xi1[2]-xi1[1])[1]\n", + " return l > 1.0e-9\n", + "end\n", + "\n", + "function calc_projection(sel, mel)\n", + " xi1 = calc_projection_master_nodes_to_slave_element(sel, mel)\n", + " xi2 = calc_projection_slave_nodes_to_master_element(sel, mel)\n", + " return xi1, xi2\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualize results" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": [ + "iVBORw0KGgoAAAANSUhEUgAAA1IAAAEbCAYAAADQ7rWKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XdYVGf6//H3DDNDLyoCihJQLCgxKvZuTDTGAqSamMhGY5obU7bE3WR3s4mb6GaN5ptsjGviTzS9gWI3iQ1BNBo1NhQVu1jpbYaZ3x+HQYkFEJgzB+7XdXGplJkPKnDu8zzPfetsNpsNIYQQQgghhBDVplc7gBBCCCGEEEJojRRSQgghhBBCCFFDUkgJIYQQQgghRA1JISWEEEIIIYQQNSSFlBBCCCGEEELUkBRSQgghhBBCCFFDUkgJIYQQQgghRA1JISWEEEIIIYQQNSSFlBBCCCGEEELUkBRSQgghhBBCCFFDUkgJIYQQQgghRA1JISWEEEIIIYQQNSSFlBBCCCGEEELUkBRSQgghhBBCCFFDUkgJIYQQQgghRA1JISWEEEIIIYQQNSSFlBBCCCGEEELUkBRSQgghRB2aPHkyM2fOxGg0otPpKl6MRiMzZ85k8uTJakcUQghRB6SQEkIIIepQWFgY06ZNw2KxVHq9xWJh2rRphIeHq5RMCCFEXdLZbDab2iGEEEKIhsJoNF5TRF3NYDBgNpsdmEgIIUR9MKgdQAghhGgIbDYbFkv2TYsooMq3CyGE0AYppIQQQoibsNmsmM0XKC09Q0nJGUpLr7z89s9Wa7HacYUQQjiIbO0TQgjRKFmtFszmrCqLo9LSs9hslVeRDIZmuLq2wGRSXlxdW1b8PjDwoSqfe/To0cTExDB27FiaN29eX5+iEEKIeiSFlBBCiAbFai2htPTsVQXR6esWS2bzOeDqH4F6TKaA6xZHyktQ+euC0OtNN3z+qs5I6XQ6BgwYQHJyMjqdjv79+xMbG0tMTAxhYWF19xchhBCiXkkhJYQQQhPKygoqFUQlJaevu5JksVyq9HE6nRGTKeiqAqlFpQLpSsEUgE7nUuucM2fOZNq0aTd8+4wZM3jllVc4d+4cSUlJJCYmsnbtWkpKSrjjjjuIiYkhJiaGO+64A51OV+s8Qggh6ocUUkIIIVSjNGjIqbI4Ki09Q1lZXqWP1evdMJlaVlEctcBobIpO57hpH5MnTyY8PJzXXnut0sqUwWBg+vTpZGRkMH/+/Eofk5+fz6pVq0hMTGTZsmXk5OQQGhpKTEwMsbGx9O/fHxeX2hd5Qggh6o4UUkIIIeqc0qDh4nUKpLPXFEu/bdDg4uJzg+KoRaXCycXFx+lXbL766ivGjRvHxx9/zKRJk6r1MaWlpaxfv57ExESWLFnC6dOn8ff3Z8yYMcTGxnLXXXfh7u5ez8mFEEJURQopIYTQuMmTJxMXF8eAAQOueVtycjLx8fHXrIDcKqVBw7kbdLC7ukCqukHD9YojkykIFxfPOsnqDM6ePUuLFi145plnmDt3bo0/3mq18vPPP5OQkEBCQgLp6el4enpyzz33EBMTw6hRo2jSpEk9JBdCCFEVKaSEEELjkpOTGTNmDElJSZWKqRu9/nqubdBwbXGkNGg4D1iv+sjKDRpuvJIUhF7vWj9/AU7OxcWF/v37s3Hjxlo/1oEDB0hMTCQhIYGtW7diMBgYMmQIMTExREdH06pVqzpILIQQojqkkBJCiAbgt0WT/c+JiV/Ts2fYdc4dna6iQYOhGsVRC4zGAPR6GUl4M02aNMHb25vjx4/X6eOeOnWKpUuXkpCQwLp167BYLPTs2bPiXFXHjh2dfuujEEJomRRSQghxi26lqUBd+G2DBntBtHz5RqZMSWLYMH/S0i7x1luuREYWVfpYpUHD9QqklpWKJaOxmUMbNDRkERERHDt2jMLCwnp7juzsbJYvX05iYiIrV66koKCA9u3bV7RV79WrF3q9/HsKIURdkkJKCCFu0VtvvcWrr756w7fb21xX17UNGm40KPb0dRo0eGMwBNKvXwYAr746hKeeuveaWUgGg6+sUjjY6NGjWb58OWVlZQ4pZoqLi/nhhx9ITExk6dKlnD9/nhYtWhAdHU1MTAxDhw7FZLrxHCwhhBDVI4WUEELcoqoGrxoMBsxmcxUNGq4UR6WlWdhs5t88RtObrhzZB8UaDF4kJyczaNAgfHyUbnbVORsl6t+rr77KW2+9xS+//ELXrl0d+txlZWWkpqZWNKs4evQoPj4+jBo1ipiYGEaOHIm3t7dDMwkhREMhhZQQQtyiqlZ2DAZYvz7wOg0adBiNATc8d3R1geTi4latLPYzUc2bN+fo0aOsW7eu2o0mRP1asWIFo0aNYs6cObzwwguq5bDZbPz6668VzSp27tyJyWTirrvuIiYmhrFjxxIYGAg4thOkEEJolRRSQghRDZZ8C0XpRRQeKFRe0guJ/Cbyph+j18Phw69fMyi2rhs0XN1o4rPPPuOjjz4iMzOTEydOSDHlBPLz8/H29mbChAnEx8erHadCZmYmS5YsISEhgU2bNmGz2ejXrx8xMTG0atWKZ599tladIIUQoqGTQkoIIcrZbDZKT5deKZaueik5WVLxfqaWJjw6eNBtXbdqPWZ9u3r1YMmSJcTExDB79mxefPFFWT1wEkajkW7durF161a1o1zXhQsXSEpKIjExkTVr1lBcXExYWBhnzpzhgw8+YMKECaSlpUkRJYQQV5FCSgihGXW13aisuIyiQ0UUplculorSiyjLLwNAZ9ThHu6OR4QHHh098Ojgofy+gwcGH2U1qbpnpBypuLgYd3d3oqOjSUxMdOhzixtr3rw5BoOBM2fOqB2lSvn5+axevbpiC2BBQQEArq6uvPPOOzz77LMYDNLyXgghpJASQmhGTQbP2mw2zBfMlVeWygun4qPFFUeWDE0NFQXS1b+6hbqhN9y8w9rMmTOZNm3aDd9e0659dcXLy4uAgACOHDni8OcW19e1a1f2799PSUlJ1e/sRMxmM23atOHkyZP4+vqSk5ND06ZNGTNmDDExMQwfPhwPDw+1YwohhCrklpIQQjMGDBhAUlJSpaJp4/qNjI0ey6fTPiUkJYQDCw5UFE2WS+WrRXpwb+OOewd3msc2x72De0XBZPK/9TbQGRkZzJgx46ZzpNQQGhqq2nOL62vXrh27du2itLRUU63H09LSOHfuHKDcnJg3bx7Hjx8nISGB+Ph43N3dGTFiBLGxsYwePZqmTZuqnFgI57B//37Gjh3L8ePHsdls6HQ6QkJCWLp0KREREWrHE3VEVqSEEJpgybFUrCxtWLuBRz97lLamtpwpPcNbvMXt3I6LlwvuHdzxjPBUtuOVrzC5h7ujd63f+T0mk4lu3bqRlpZWr89THRMmTGDx4sWcP38ef39/teMI4O233+avf/0rGzZsYNCgQWrHqRb7Su+kSZOYNWsWb775JrNmzaq4iXHw4MGK7X9btmzBxcWFQYMGERsbS3R0NCEhIWp/CkKoYu/evXTt2vW6W78NBgO7d++WYqqBkEJKCOE0bFYbxceLK3fHK38pPVta8X7mYDPDTw0H4PUHXmfqs1Px6OCBqaVJtWGz9ovIdevWqfL8V1u0aBFxcXF8/PHHTJo0Se04AqUoGThwINOnT7/pEGdncfV2WQ8PD6KionjttdcYMWLEdbfXnjlzpqID4Lp16zCbzURFRRETE0NsbCydOnWSQdCi0WjXrt1NdwWEh4dz6NAhByYS9UUKKSGEw5UVlFF4sPDaZg8Hi7AWKYeX9G563Nu7X3N2yaO9By/95SXef/99JkyYwNKlS52ii5her2fkyJEsX75c1RwA2dnZNGnShHHjxvHFF1+oHUcApaWluLq6cv/99/Ptt9+qHadKVzd2sVgsGI3GiuxVNXbJyclhxYoVJCYmsmLFCvLz82nXrh0xMTHExMTQp08f9Pr6XSEWQk2urq6Ulpbe8O0mk0lz5yXF9UkhJYSoFzabjdKzN2glfvzKDxBjoPHaYqmDB24hbuhcrr2DnZyczJAhQygrK8NsNrNlyxanaMms0+l48MEH+frrr1XLcDV3d3dCQ0PZv3+/2lFEOTc3Nzp06MCuXbvUjlJjJpOJLl268PPPP9fo40pKSvjxxx9JTExkyZIlnDt3jsDAQMaOHUtsbCx33nknrq6u9ZRaiPqVmwv798O+fcrL3r3Kr8eOmYAbd2w1Go03LbSEdkghJYSoFWuplaKMa7fiFaYXUpZb3krcoLQSd+/grpxd6nilpbixibHaz2XfblRSUkLTpk05efJkpderVUzZVxueeOIJFixY4PDnv5527dpx+vTpitbVQn0tW7aktLSUCxcuqB2lxvz9/XF1deXUqVO3/BhlZWVs2bKl4lzV4cOH8fb2ZuTIkcTGxnLvvffi4+NTh6mFqBs5OVeKpasLphMnlLfrdBAWBp06QefOMGuWKxaLrEg1BtK1TwhRLeaL5kotxCu24x0pAqVewuCntBL37OxJ8weaK40eOrjj3tYdvbH2W3ni4+P55JNPuP/++ysd2Ld384uPj1elkLp06RKgtB13Fl26dCEjI4PCwkJpT+0kWrduzY4dO9SOcUuaN2/OsWPHavUYLi4u9O/fn/79+/Pvf/+bvXv3kpCQQEJCAo888ghGo5Fhw4YRGxvL2LFjCQoKqqP0QlTP5cvXFkv79oH9/oFOB23bKgXT+PFK0dSpE3ToAJ6eVx7nu+9CbnpGShqxNBxSSAkhKtjKbBRnFlOwv6DSkNrCA4WYL5RvU9CBW5gbHh08aDa62ZXteB09MDY31uuB8vnz5/P3v/8dgCeeeKLS2wYMGKDa1r6LFy8C4O3trcrzX8+wYcP4/vvvWbFiBQ888IDacQTQsWNHtm7dSm5uruZWXkJCQkhPT6+zx9PpdERGRhIZGcnf/vY3jh07VtGs4rnnnuOZZ56hT58+Fc0q2rVrV2fPLbQrPj6ewYMHExoaes3bMjMz2bBhA3FxcVU+zqVLlQsl++/t87L1eggPV4qk3/1O+dVeMLm7V51z6dKldOnS5YZd+5YuXVr1gwhNkK19QjRCllzLlZWlq1aYig4VYStVviXoPfUV2+8qfo1QWom7uLuolr179+7s3r2b0tJSpzmwbu/INnPmTP785z+rHQeA06dPExwczKRJk/j444/VjiOA999/n6lTp5KUlMTo0aPVjlMjzz//PB988AHHjh2r97vpFy9eZNmyZSQmJrJ69WqKioro3LlzRVHVvXt36QDYSGVmZjJx4kQWLFhQqZi60evPn7/+lrysLOXtLi7Qrt2VLXn2gql9e3Bzq11WmSPVOMiKlBANlM1mo+RkyXWbPZSevrJ32xRswqODB36D/Wj5TMuKosm1latTXqzs37+fNm3aOE0RBVe29vn6+qqc5IqWLVtiNBpr3BxA1J/BgwcDkJKSorlC6vbbbwdg69at9V5INWvWjLi4OOLi4igsLGTNmjUkJCTw4Ycf8q9//YvWrVsTHR1NTEwMgwYNwmis/jlLoW2hoaEsWLCgUtF09Ggmjz8+kSlTFrBsWWilwun8eeXjDAalOOrUCZ5++krR1L491Nd87IiICGlx3ghIISWEA8XHxxMcHMyzzz57zV2quXPncurUqWptS7haWVEZRYeKKjV5sP9qLVBaietMOjzaK+eVgp4IqrTSZPDWzreBPXv2UFxczF133aV2lEqys7MB8PPzUzlJZUFBQRw9elTtGKJcZGQkgCa79vXq1QuAX375xaFbRT08PCrappvNZjZt2kRiYiKJiYl88MEHNGnShNGjRxMTE8OIESPwvPqgimhQbDY4exYOHw5lwID5dOrUl+DgN8jM/AKLZQGbN4diNCrb7zp1gjvvvLLC1K4dSL0t6oN2rqCEaABatGjB8OHD+e2O2oyMDIYPH86aNWuu+3E2mw3zOfO1xdKBQoozi6H84YzNlVbi3lHeBI4PrCiW3EKv30pca+bNmwfAU089pXKSynJycgDnK6Q6d+7MqlWrsFgsGAzy7V5ter0ed3d3jhw5onaUGuvSpQsABw4cUC2D0Wjkzjvv5M477+S9995jx44dFUXV4sWLcXNzY/jw4cTGxjJ69Gj8/f1Vyypunc2mnFW6eiue/dfLl5X3MRjysVjOkpHxFE88sZ7Ro0Pp1Ek51yTf6oQjyX83IRxoypQp1xRRdjabjWeefobdK3ZXNHuwN3ooPFCIJbv80KoLuLdR2og3v795RaMHjw4eGJs17Ftua9aswdXVla5du6odpZLc3FwAmjZtqnKSygYPHsyqVatYt24dd999t9pxBEr3uzP2E+0aotfrMZlMZGZmqh0FUJpVREVFERUVxZtvvklGRkZFW/WJEyei0+kYNGhQxWrWbbfdpnZk8Rs2m9IN77fF0r59SrtxUM4pdeyobMW7915ldSkiAr76ahmvv65j1apVzJjxT7p3X3DdBhRC1DdpNiGEA1U17dyAgbWsBcDFx6VSkwf7793D3dGbnOd8kCMZjUY6d+7Mzp071Y5SySuvvMK///1vMjIyaNu2rdpxKqSnp9OxY0emTp3Ke++9p3YcgVLcJicnU1ZWpnaUGvP398dkMnH69Gm1o9zU2bNnWbp0KYmJifz444+UlpbSrVu3imYVkZGRTnn+s6Gy2ZR5S9crmPLylPdxd1cKJHvTh4gI5dewMKUhxNUyMzPp06cPgYGB7Nq164aNJoRwBFmREsKBqrxv4QJ3rL0Dj44emIJM8sP+KuvXr8disTBq1Ci1o1wjr/xqoFmzZionqaxDhw64uLiQlpamdhRRrlOnTmzcuJHTp0/TsmVLtePUSEBAgNOsSN1MUFAQTz31FE899RS5ubmsWrWKhIQE/vOf//CPf/yDtm3bVqxU9e3bF5ffXqmLW2K1wrFj13bI278f8vOV9/HwuFIkxcZeKZxCQ5WW41WxF01ubm4MGTIEuH4DCiEcRQopIRyoqsJI76KnydAmDkqjLfYW3k8//bTKSa6VX36V4IyzgZo3by6do5xIr169+Oijj9iwYQOPPPKI2nFqJCQkRNUzUrfCx8eHhx56iIceeoiSkhLWrVtHYmIin332GbNmzSIgIICxY8cSExPDsGHDcKttz+tGwGqFzMxrV5j274fCQuV9vLyuFEwPPnilS15ISPUKphvZsGED//73v+nZsyf9+vWreL29mNqwYYMUUsKhZGufEA7Url27m047Dw8Pl4veGwgJCeHSpUsVRYszue+++0hISKh6xVEFQ4cOZcOGDVgsFqdqGd9YHT16lDZt2mhyu+XUqVN5//33yczM1PyZI6vVypYtWyrOVWVkZODl5cXIkSOJiYlh1KhRTjXOQA1lZXDkyLVzmA4cgKIi5X18fK4UTPbVpc6doVUrqK8NFUuWLCEmJobjx4/TunXr+nkSIapJVqSEcKC5c+det2sfKKtVc+fOVSGV87NYLJw8ebLSHUhnUlBQ4LTbMPv168f69evZtm0bvXv3VjtOoxcWFoZOp2PPnj1qR6mxO+64A4C0tDTNF1J6vZ5+/frRr18/Zs6cyb59+yqKqvHjx2M0Ghk6dCgxMTFER0drbhtmTVgscPjwteeXDhyAkhLlfXx9lQIpKgomTLjSVjw4uP4KphtJTU0lODhYiijhFOT2pBAOdOrUKRITE9HpdLi4uGA0GjGZTISHh7NmzRpOnTqldkSn9P3332Oz2bjvvvvUjnJdhYWFTrvaM3bsWACWLl2qchJh5+Xlpcn5XvZC3NmavdSWTqejc+fOvPrqq/z8888cP36cd999F4vFwvPPP09wcDB9+vRh5syZpKenqx33lpnNyva7776DN96AceOgSxfw9FQ64913H7z7rtIYolcvmDkT1qxROutdvgybN8P8+fDiizB8eP2uOt1MSkqK095UE42PbO0TwsEWLlzIxIkTyczMJCQkRO04mhAdHc3SpUu5fPmy081qAujRowe//vorJfbbt07EarViMBgYNGgQ69evVzuOQNnie+bMGafcpnozVqsVFxcXYmJiSEhIUDuOQ1y6dInly5eTmJjIqlWrKCwspGPHjsTGxhITE0OPHj2c7iZKaSkcOnTtCtPBg0oxBdCsWeXtePYVpsBAdYqj6jKbzfj4+PDWW2/x0ksvqR1HCCmkhHC0oUOHotfr+fHHH9WOohkBAQGUlZVx8eJFtaNcV2RkJEeOHKHQftLayTRr1gyTyaTJ+UUN0YgRI1i7dq0mz625urrSqVMnfvnlF7WjOFxhYSFr164lMTGRpKQkLl68SHBwMNHR0cTExDBkyBCMRsfN8ispUYqj33bJO3RI2a4HEBBw5QzT1QVTQIDDYtapbdu20atXL7Zs2SJblYVTkDNSQjjQsWPHWL9+PQsXLlQ7imbk5+dz/vx5RowYoXaUGyouLnboBVRNtWvXjp9//lntGKJcZGQka9as4dChQ3To0EHtODXi4+PD2bNn1Y6hCg8PD6Kjo4mOjsZisZCcnFxxrurDDz/Ez8+PUaNGERsby4gRI/Dy8gJg8uTJhIeH89prr2GxVziAwWBg+vTpZGRkMH/+/Bs+b3ExpKdfWzBlZCgNIQCCgpQCadgwmDr1SvHk71+vfyUOl5qaiqurK926dVM7ihCAFFJCONSnn36Kh4cH999/v9pRNGPRokUAPProoyonubGSkhKnLqR69+5NWloa6enpmrtwb4j69u0LKK2ctfbvERAQwJEjR9SOoTqDwcCQIUMYMmQIs2fP5pdffmHJkiUkJCTw2Wef4ebmxt13301MTAwBAQFMmzbtmsewWCxMmzaNGTNmAEonvPT0a9uKHz6stBwHaNlSKZhGjICXX75SMDVt6sjPXj0pKSn06NEDk8mkdhQhANnaJ4TD2Gw2OnToQJ8+fSqKA1E1e/vu4uJip/3hGRgYiF6vd9qtc2vXrmX48OG8/fbb172gE4514cIFmjdvzpNPPnnTlQhndO+997Jy5UqnbPXvLA4fPkxiYiKJiYls3ry5Gn9XLWjb9jRHjoD9XVu1urIN7+oteU54RNShQkJCePjhh3nnnXfUjiIEIF37hHCYtLQ0Dh06xIQJE9SOoik7duwgKCjIaYsogNLSUlxdXdWOcUNDhw4FlBUQoT5/f39cXFzYt2+f2lFqzL6CduzYMZWTOK+2bdvyhz/8gU2bNlXz5sp5oqOVjnipqZCdrXTOW70aZs+GJ5+Efv2kiDp58iQnTpyQjn3CqcjWPiEcJD4+nlatWlVc1IqqnT17ltzcXKc+HwVKJylnLqQMBgM+Pj7s3btX7SiinI+PjyaLkS5dugDKFiutz5JyhMDAwGq8l4VZs+o9iualpqYCV7bGCuEMZEVKCAcoKSnhyy+/5LHHHsPFxUXtOJrxv//9D4Df/e536gapQllZGW5ubmrHuKmwsLBG2yTAGbVo0cJpu1DeTK9evYCGN0tKOL/U1FTCwsIICgpSO4oQFaSQEsIBkpKSyM7Olm19NbR06VL0ej333HOP2lFuymKx4O7urnaMm+rRowdms5mTJ0+qHUWgdFIsLi6u1MVNCyIiIgA0PZjW0QyGm2/+qertQiGDeIUzkkJKCAdYtGgRPXv2rLgIEdWzb98+QkNDnX7WjtVqxcPDQ+0YN2UvRhMTE1VOIuDKFjmtzWPS6/WYTCYyMzPVjqIZ06dPr9XbhTJiYseOHbKtTzgd5746EaIBOHfuHCtXriQuLk7tKJqSnp5OUVERw4YNUztKlWw2m9MXUvfeey8A69atUzmJAOjfvz8AGzduVDlJzfn6+pKVlaV2DM3IyMhgxowZ16w8GQwGZsyYQUZGhkrJtGP79u2YzWZZkRJOR9aThahnX3zxBTqdjnHjxqkdRVM++ugjAJ566imVk1TNZrNVDN90Vh4eHnh4eLB79261owhg4MCBgNKVUmsCAwPl4r8G7C3up0yZgre3tybb3qstJSUFT09Pbr/9drWjCFGJrEgJUc/i4+MZPXo0zZo1UzuKpqxatQqTyUSPHj3UjnJTxcXFAE6/IgXKDBY5I+UcPDw8MBqNmjxrFBISQnFxMVb7lFhRLfZW6E0by/TcOpSamkqvXr3kPJlwOlJICVGPfv31V3755RdpMnELMjIyKmbWOLNLly4B4O3trXKSqnXr1o3i4mKys7PVjiIAPz8/Tpw4oXaMGrN/XR49elTlJNpiL6T8/f1VTqItNpuNlJQUOR8lnJIUUsIhsrKy6N+/vyYHUNbG4sWLadasWcX5FFE9ycnJWCwWRo4cqXaUKtkLKR8fH5WTVM1+3mzp0qUqJxEArVq14vLly2rHqDH7TYMPP/xQ5STaYh8/EBAQoHISbcnMzCQrK0vORwmnJIWUcIgPPviAXbt20aJFC7WjOIzFYuHTTz/lkUcewWQyqR1HUz7++GMAnnnmGZWTVE1LhVRsbCwAP/zwg8pJBCgrO2azuWJ7qFb07NkTgE2bNqmcRFvOnTsHSCFVUykpKQD06dNH5SRCXEsKKVHvCgsL+fDDD5k0aRJNmjRRO47D/Pjjj5w5c0a69d2C9evX4+HhQVhYmNpRqmRfUfD19VU5SdWaNm2Kq6sr27dvVzuKALp37w7A5s2bVU5SM/ZW+sePH1c5ibZcuHABoFHdUKwLqampdOjQQc4ZC6ckhZSodwsXLiQ7O5sXX3xR7SgOFR8fT0REBFFRUWpH0RSr1cqJEycq5uw4u5ycHEA576IFLVu25NixY2rHEMCAAQMAZSurltgP/Nv/74vquXjxIqB8DYrqk0G8wplJISXqVVlZGbNnz+aBBx7QxOpCXcnNzSUhIYEJEyag0+nUjqMpS5YswWq1VmxDc3b2i0mtrLZ26dKFgoICzW0na4jsW+R27dqlcpKa0+l0lJSUqB1DU+zbgKXZRPXl5+eza9cuaTQhnJYUUqJeLV26lIyMDP7whz+oHcWhvvnmG0pKSnjsscfUjqI5ixcvBuDJJ59UOUn1aK2QGjJkCACrV69WN4jAYDDg5ubGoUOH1I5SYyaTCZvNpnYMTbF/r9Dr5dKrurZt24bVapUVKeG05KtZ1KtZs2YxcOBAevXqpXYUh1q0aBHDhg2jVasXU6EJAAAgAElEQVRWakfRnM2bN+Pn56eZWSu5ubkAmtm/b1/pW7VqlcpJBCjn1k6fPq12jBqzN1exWCwqJ9GO3NxcXFxc1I6hKSkpKfj6+hIREaF2FCGuSwopUW+2bNnC5s2bG91q1NGjR9m4caM0mbgFxcXFnDt3zumH8F4tLy8P0M6K1G233YbBYGDr1q1qRxEow221eNaodevWAKxbt07lJNqRl5cnA2VrKDU1lT59+sgqnnBa8j9T1JtZs2bRrl07xowZo3YUh1q8eDFeXl6aOePjTBYtWgTAuHHjVE5SffZCSgvtz+0CAwM5cuSI2jEEEBERQVlZWcX5Ga2w7zL4/vvvVU6iHQUFBRiNRrVjaIbNZiM1NVXORwmnJoWUqBdHjhzh+++/5+WXX25Ud5JsNhuLFi3igQcewNPTU+04mvP1118DMH78eJWTVF9BQQGgrXMPnTp1Ijs7G6vVqnaURs/ecGLDhg0qJ6mZ++67D0BWNmugsLAQV1dXtWNoxsGDB7l06ZKcjxJOTTs/+YWmzJkzh6ZNmzJhwgS1ozhUSkoKhw8fbnSfd13Zvn07QUFBuLm5qR2l2goKCjTXmXHQoEEAbNy4UeUkwv5vkZqaqnKSmhk6dCgAJ06cUDmJdhQVFeHu7q52DM1ISUlBp9PRu3dvtaMIcUNSSIk6d+nSJT755BOee+45PDw81I7jUIsWLSIkJITBgwerHUVzLly4QHZ2tubuPhYWFmpqNQogJiYGgKSkJJWTCPsh+t27d6ucpGbsZ33szVZE1UpKShrdz8TaSE1NJTIyUlPbpkXjo62f/kIT5s2bR1lZGVOmTFE7ikMVFRXx1Vdf8fjjj2vuwtoZzJ8/H0Bzq3nFxcWa68QVGRmJXq9ny5Ytakdp9PR6PZ6enpo8s6bT6SgtLVU7hmaYzWa8vLzUjqEZMohXaIFc7Yla279/P+3atcPV1RWj0chf//pXXF1dK6a4NxZJSUnk5ORorhBwFomJiej1es01JykqKtJkJy5/f39Nzi9qiPz9/cnKylI7Ro25urrKLKkasFgseHt7qx1DE7Kzs9m7d680mhBOTwopUSt79+6lS5cuZGRkUFpaWjFTJDc3ly5durB//36VEzpOfHw8ffr0oX379mpH0aQ9e/YQEhKiudW84uJiTRZS7du3b3Q3O5xVWFhYRfdHLZg8eTIzZ86kuLgYUFamdDodRqORmTNnMnnyZJUTOier1Yqvr6/aMTQhLS0NQFakhNPT1hWLcDoxMTE3HMhosVgYO3asgxOp4+zZs6xevVpWo27R4cOHKSwsrDjAriUlJSWabGncr18/rFYrO3fuVDtKoxcZGYnNZuP48eNqR6mWsLAwpk2bds3rLRYL06ZNIzw8XIVUzs3+c9LPz0/lJNqQkpKCv7+//F8STk8KKVErVf3g18qFQW19/vnnuLi48PDDD6sdRZM++ugjAE3eyS4tLcVkMqkdo8bsWygTExNVTiLsXcnWr1+vbpBq+sc//nHTt7/22msOSqIdZ8+eBaBp06YqJ9EG+/worXVEFY2PFFKiVqraH19WZqMxbKFftGgRY8aMkR+St2jlypUYjUZN7oc3m82aLKT69euHTqcjOTlZ7SiN3pAhQ4Ar25mc3Y12IVT37Y3RmTNnAGjWrJnKSZxfWVkZW7Zs0eTPA9H4SCElaqWqu0VlZTp69IDPPwez2UGhHGzXrl3s2rWLuLg4taNo1sGDB2nXrp3aMW6J2WzW1NwrO71ej5+fH/v27VM7SqPXqlUr9Ho9e/fuVTtK3dDrSTx/nizp6FfBXkg1b95c5STOb9++feTl5cn5KKEJUkiJWgkJCbnp21u0CMHfH8aPh7ZtYdYsyMlxUDgHWbRoEc2bN+eee+5RO4ompaWlYTabNfv3Z7FYNFlIAbRp04bz58+rHUMAXl5eZGZmqh3j5s6fhz//uer30+uJ3buXoJQU2m7ZwuP79/PhqVPszMvDYrXWf04ndO7cOQACAgJUTuL8UlJScHFxoUePHmpHEaJKUkiJWlm6dOkNO5YZDAZ+/HEpq1fDrl1w553wl79A69bwxz/CiRMODlsPLBYLn332GY8++qgmGw44A/v8qGeeeUblJLemrKwMd3d3tWPckt69e2OxWDh8+LDaURq9wMBA5y1qL12CV1+FNm1g7lwMVexEMADH+/Thy06dGNOsGemFhbyQkUG37dvxS05m2M6dvHbkCCsuXuRSQ92q8Bv2f9ugoCCVkzi/lJQUunbtiqenp9pRhKiSFFKiViIiIti9ezfh4eGYTCaMRiMmk4nw8HB2795NREQEAF26wMKFkJkJU6bAJ59AWJiyUrVjh6qfQq2sWbOGrKws6dZXCz/99BPu7u6a3dpntVrx8PBQO8YtGTlyJCANJ5xB27ZtKSwsxOpMKzbZ2fD668o36zlz4Lnn4OhRpr/99k0/bPr06bR2c+PhgADmtGvH1qgocgcMYGPXrvw9NBQfg4H5Z84w6tdfabZ5MxFbtzLpwAE+Pn2afQUFWBvgwdoLFy4A0LJlS5WTOL/U1FTZ1ic0QwopUWsREREcOnSIkpISSktLKSkp4dChQxVF1NVatoS331ZWo959F1JSICpKWa1asQKc6RqiOhYtWkRkZCTdunVTO4omWa1Wjh07RmRkpNpRbpmWC6nhw4cD2ukW15B16dIFwDnOrOXlwfTpSgE1cyZMngxHjyq/9/cnIyODGTNmXLMbwWAwMGPGDDIyMq55SHcXFwb6+fHnkBASIiM5268fGb17s7hjR+7082NHfj5PHzxI523baLZ5MyN37+bNzEx+uHSJvAbQvOLSpUsAtGjRQuUkzu38+fMcOnRIGk0IzdDeFEnRIHh5wdSpyg3O779Xzk6NGgUREfDyy/DYY+Dsx06ys7NJTEzkzTfflBatt2jFihVYrVaio6PVjlIrWt2CYjKZ8PLy4tdff1U7SqNnvwO/ceNG9W4s5OfDf/8L77yj/P7pp2HaNPjNxb99O+4rr7xyy0+l0+lo6+5OW3d3Hivf7pZnsbA1L4/UnBxScnOZffIkf7dY0AO3e3rS19eXvj4+9PXxIdzdXVPfd7OzswE02eHTkbZs2QLIIF6hHVJICVUZDPDQQ/Dgg7B5s1JQPfWUsh3/+efh2WfBWbvFfvPNN5jNZsaPH692FM2Kj48HtDk/CqCwsBBQGgVoVVhYGOnp6WrHaPQGDx4MwLZt2xz/5IWFMHeusuKUnQ1PPgl//Su0auXQGN4GA8OaNGFYkyYAWG020gsLScnNJSUnh43Z2Xx0+jQA/kYj/cqLqn6+vvTw9sbDxcWheWsiJycHvV42AVUlJSWFFi1aVNnISghnIYWUcAo6HQwYoLwcPAizZ8O//gVvvQVPPAEvvQTONuA8Pj6eu+++W/a810JycjK+vr6a7WRl367j7e2tcpJb1717d3799VfOnTun2X+HhsDPzw8XFxf279/vuCctLob//U/Zb33+vPLN9rXX4LbbHJfhJvQ6HRGenkR4ejKpfFXsstnMltxcUsuLq38dP05+WRkGnY6uXl5KYVVeXLV2dXWaVavc3NwbNmYSV9jPRznLv5sQVZHbI8LptG+v3Bw9flzZVfLNN8rr7rtPWbVyhnPIhw8fZvPmzdJkohaKi4s5e/Ys3bt3VzvKLWsIhdSIESMAaTjhDHx9fTl+/Hj9P1FpqfJNNjxcuUs1YgSkp8P8+U5TRN1IE6ORkc2a8UZYGD907Ur2gAHs7NGD98PDifDwYMXFizyyfz+3bdlCq9RUHty7l9knTrAlJ4cSFQ/hFhQUSCFVBbPZzNatW+V8lNAUKaSE02reHP7+dzh2DObNg/37lRWrvn3h22+hrEy9bIsWLcLb25uYmBj1QmjcF198AcDDDz+scpJbZy+kfH19VU5y68aMGQPADz/8oHISERwcXPF/ql6YzfDxx8qdqSlTYMgQ5RvrwoXKoD8NctHpuMPLi2eCg1kUEUFGnz5k9etHYmQkjwcGklVayl+PHqXvL7/gu2kT/Xfs4E+HD5Nw/jxnS0oclrOgoABXV1eHPZ8W7dq1i6KiIjkfJTRFbo8Ip+furjSNmjRJ6ew3a5ZypiosTLmZ+sQTSvMKR7FarSxatIgHH3xQs93anMGXX34JwOOPP65ykluXUz5dWsuFlJeXF+7u7uzatUvtKI1eeHg4v/76KxaLpW5XLywW+OwzeOMNOHJEOZi6YgV06lR3z+FEAkwmov39ifb3B6DUamVXfj6publszsnhy3Pn+E/5IMMwN7eKc1Z9fXzo4umJoR7OMhUXF2t2cLejpKamYjKZNL1LQTQ+siIlNEOvh9GjYd06+Pln6NNHKaRCQpRz0WfOOCZHcnIymZmZxMXFOeYJG6ht27YREBCg6WLU3onLz89P5SS106pVK040hAnZGte1a1cA0tLS6uYBy8rg88+hc2f43e+UgX67dsFXXzXYIup6THo9PX18mNqqFV917syJvn050acPX3fqRLS/PxlFRbyUkUHU9u34JiczdOdOXj1yhGUXLnCxjgYGl5SUaHZwt6OkpKQQFRUlK3dCU6SQEpoUFaVcHxw+rKxIffABhIYqv9+zp36fe9GiRYSGhjJgwID6faIG7NKlS1y+fFnze+HthVST8i5jWnXHHXdQVFREfn6+2lEaNfv3lM2bN9fugaxW5XBply7K1PP27WH7dkhIUF4naOXmxoMBAcwODyetfGBwcrduvB4aShODgU/OnGHMnj34b95Mx7Q0njhwgPmnT7MnP/+WBgaXlpZqdkyCo8ggXqFFUkgJTbvtNmWr34kTyvzItWvh9tvhnnvghx/qvjFFYWEhX3/9NRMmTJBWtrXwySefAPDYY4+pnKR2GsLWPoC77roLgOXLl6ucpHGzX0Tu2LHj1h7AZoPEROjWTdm+17o1bNkCSUkg26Vuys3Fhf6+vvwpJITvIyM5068fh3v35tOICIY1acKu/HyePXiQ23/+mSbJyYzYtYt/Zmay9tIlcm8wMHjy5MkkJycDYLFYKjWlSU5O1uzYh/pw+vRpjh07pvmba6LxkStB0SD4+sKf/qRs/1+8GLKy4O67leuJxYuVJlV1YcmSJeTl5Wn6XI8zSEhIQKfTab5ZR15eHgDNnHXYWTXZByKvWrVK5SSNm5ubGyaTqeZzvWw2WL4cevSA2Fjw94fkZFi1Cnr3rp+wDZxOp6ONuzvjAwP5b/v27OjRg+wBA/jpjjv4c0gIRp2O/zt5kuG7d+OXnEyXbdt4Oj2d+LNnOVhYiM1mIy4ujjFjxpCcnExZWRk+Pj6AUkSNGTNGtodf5ZtvvgGgT58+KicRomZ0NpszNJMWom7ZbPDTT8pq1cqVEBwMU6cqw35rc5zlnnvuIT8/v+Iuo7g1Xl5eNGvWjGPHjqkdpVaeeeYZ5s2bR05OTsVFkla5urrSoUMHdu/erXaURi0oKAir1cq5c+eqfmebTVmG//vfIS1NaWv65ptKNz5R76w2GwcLC5WZVuVzrfaVD+n2Nxrp4+NDywMH+HzyZPJzchg3bhxTpkxhzJgxJCUlyfbwq/To0YOdO3diucHqnhDOSlakRIOk08GwYUpjqj17lDEpf/ubstPlpZcgM7Pmj3n69GnWrl0rs6Nq6dixYxQUFDCkAVzs2VekvBzZNrKetGjRgsxb+cIQdapVq1YVZ+9u6qefYOBA5ZubTqcUVBs3ShHlQHqdjo6enjzRogXzO3Rgb69eXO7fn1VduvBcy5aUWK18ERRE/iuvAJC0bx93jRrF1IULaRUVhdzHviI9PZ1WrVqpHUOIGpNCSjR4nTvDJ58o86heeAHi45U5lOPGwbZt1X+czz//HKPRyEMPPVR/YRuBjz76CIBJkyapnKT2CgoKABrEebnIyEjy8vIorat9sOKWdOjQAbPZfOPGH5s2wdChyp2ikhJlS19KCtx1l1JQCVX5GY2MaNqUf4aFseaOO7g8YACbHn+cVp07U7B7N94TJ/KGry9haWkEp6Zy/549zDpxglSVBwarKTs7m/z8fHr06KF2FCFqTPs//YWopqAgpSHFiRPw3ntKC/VevWDwYFi6VGl0dSM2m434+Hiio6M13+pabcuXL8dgMDBo0CC1o9RaQUEBugZy8Tp48GBABvOqxd6YICoqCoBNmzZVvC05OZnJY8fC8OEwaBBcvgxLlsDWrXDvvVJAOTEXnQ4yM8k/dYo5c+ZgWbiQpVYrSyMjiQsK4oLZzN+OHqXfL7/gs2kT/Xbs4A8ZGXx3/jynHTgwWE1fffUVgObPzIrGSQop0eh4esKUKZCeDt99p8yqjI6GiAj43/+gqAji4+MrbXPauXMne/bsYcKECWRmZhIfH6/eJ6Bx6enphIeHqx2jThQVFeHi4qJ2jDpx3333AbBy5UqVkzRO9sYE9s5uKSkpACTPn8+YoUOJS0qC06fh229hxw4YO1YKKA2wN5ZISkrihRdeICkpiQn330+TAwd4u00bNnTrRs6AAWzr3p3/tG1LqJsb354/zwN79xKcmkpoaiqP7tvH+ydPsj0vD3MDXLVatmwZcOV7kBBaIs0mhABSU+E//1HGrPj7w6OPZrJ9+0QWL15AaGgoL774Il988QWbN2/mqaeeYsEC5fWiZnbs2EFUVBRTp07lvffeUztOrUVFRbF3716Ki4vVjlInDAYD3bt3Z+vWrWpHaZTsF93Z2dmMHTyYP5WVMSY5maTWrRnwzjvw4IPKZHKhCVcXUVc3lrjR6692qqSE1JycikYWO/LyKLXZ8NDr6entTT9fX/r6+NDXxwd/k8lRn1K9CA4OpqCgoHpnA4VwMlJICXGVw4dh9mz4f/8PLJZM/P0nsmDBPB5/vD9jxozh6NGjUkTVgr3L3b59+4iIiFA7Tq116tSponlGQ9CyZUuKi4u5dOmS2lEareT4eAb+7ncEAiV6PUl/+QsD/vlPaCArn43J5MmTiYuLu26xlJycTHx8PPPnz6/WYxWXlbEjP18prHJySMnN5Wz5ecZ27u708/Ghr68v/Xx86OTpqWwp1ACr1YrRaKR3794Vq7BCaIkUUkJcx8WL8NFHMGdOJhcujAV+pWPHnqxY8TVhYaHqhtOw8PBwTp06RVFRkdpR6kTbtm25cOFCxWBerRs2bBjr1q3DYrE0iAYamnLgAPzzn/DVV3QxGPjVbGbOu+/ywksvqZ1MOCGbzcax4uKKtuupubnsys+nDPBxcaF3+WpVP19fent742c0qh35upKTkxk4cCCvvfYab775ptpxhKgx+UkpxHU0awavvgonT4bywAOPAXDgwDs8/HAoX32lnKsSNWO1WsnMzKRz585qR6kzJSUlGJ30AuVW9O/fH5vNRmpqqtpRGo+MDJgwQWkvunkzyX/8Iyc8PZkzZw6vv/GGzKwT16XT6Qh1d+fRwEA+aN+e7T16kDNwIOvuuINpISGYdDo+OHWKe3bvpunmzURu3cpT6eksPHOG9PKBwc7gyy+/BOCxxx5TOYkQt8agdgAhnNmZM5lcvLiKdevW8+KL/8RoXMC4caHcdhu8+CJMmgTlZ8NFFdauXUtZWRljxoxRO0qdKS0txaTx8wlXi4mJ4c0332TZsmX0799f7TgN29GjyvDcRYsgMBDef5/kDh0Y88ADFWdnoqKiZHirqDZPFxeGNGnCkCZNAGXV6lBRUcWKVUpuLh+fOYMNaGYw0Oeq7YA9vb3xMjj+knDjxo0YjUY6dOjg8OcWoi7I1j4hbiAzM5OJEydWnImy//kPf1jAl1+G8uWXSgfAp5+GqVMhOFjtxM7tkUce4csvv+TMmTMEBQWpHadO+Pn50aRJE44ePap2lDphtVoxGAwMGDCAjRs3qh2nYTp+HP71L1iwQFn6njYNnn6a5O3bb7kxgRDVlWOxkJabW3HWaktuLrllZbgAXby8Kp21CnVzq/fxDh4eHrRs2ZKMjIx6fR4h6osUUkJcx2+LqOu93mAI5b33lJbphYXwyCPwhz/AHXeol9uZBQcHk5eXR25urtpR6oyXlxetW7dm//79akepM82bN0ev15OVlaV2lIbl9Gl46y2YP19Zxn7lFXjuOeVuDHXbmECI6iqz2dhfUFCxYpWak0N6+RnWQKORfuVFVV9fX6K8vHCrw6Yn586dIzAwkEcffZTPPvuszh5XCEeSQkqI64iPj2fw4MHX7c6XmZnJhg0biIuLAyA3Fz7+WBnye/w43HUX/PGPyuxMjTROqnelpaW4ubkxcOBANmzYoHacOuPm5kZERAS//PKL2lHqTL9+/UhLS6OsrEztKA3D2bMwcybMnQseHvCnP8Hvfy97goXTulBayhZ7YZWby9bcXAqtVow6HVHe3koTi/LiKtjV9ZafZ86cObz00kt8++233H///XX4GQjhOFJICVFHzGZlVuZ//qPMy4yMVFaoHnkEavGzpkH49NNPefzxx3n//ff5/e9/r3acOmM0GunZs2eDatv78ssvM3v2bPbs2dOgGoM43Pnz8M478MEHYDQq3wxeeAF8fdVOJkSNWKxWdhUUVJprlVk+Oy/E1bWiO2BfHx+6enlhvEHHz/j4eIKDg3n22Wc5fvw4ZrMZm81GmzZtmDdvHqdOnaq4QSmEVkghJUQds9lgwwaYNQuWLYMWLeD555WzVE2bqp1OHffeey8rV64kLy8PLy8vtePUGRcXF4YOHcoPP/ygdpQ689NPPzFs2DCmT5/Oq6++qnYc7bl0Sfnif+89ZUn6hReUIqq8AYAQDcGZkpKKc1apubn8XD4w2F2vp4e3d8WKVV8fHwLKG/KsWbOGe+6557odA3U6HWvWrOGuu+5y9KciRK1IISVEPdq/Xxnwu2iRMk9z0iSl21+bNmoncyx/f390Oh3nz59XO0qd0ul0jB49mqSkJLWj1BmLxYLRaGT48OGsXr1a7TjakZ0Nc+YoX/AWi3L35I9/BH9/tZMJUe9KrFZ25OWRWr4dcHNODmfKBwaHu7vT18eHlSNHcuHYsRs+Rnh4OIcOHXJUZCHqhBRSQjjAuXPw3/8qL5cvw333KddYvXurnaz+5ebm4uvry6hRo1i2bJnaceqUTqdj3LhxfPHFF2pHqVN+fn54e3tz4sQJtaM4v7w8+L//U/b0FhcrDSReeQUCAtROJoRqbDYbx0tKSM3JqRgavL1v35sOYTSZTJSUlDgwpRC1JwN5hXCAgAD45z+VZhT//S/s2gV9+sCAAZCYCA35XP+CBQsAGD9+vMpJ6lZhYSFAg9qqaNemTRvOnj2rdgznVlAA//43hIXBG2/A44/DkSPKtj4pokQjp9PpuM3NjXGBgfxfu3b83KMHxiq6L8l9faFFUkgJ4UAeHvDMM3DgACQkKEcoYmOhY0elqVf5tXmD8t1336HT6XjwwQfVjlKnLly4AIB3A+y+1qtXLywWC8ePH1c7ivMpKoJ331UKqNdegwcfhMOHlVWpFi3UTieE06pqJlV9z6wSoj5IISWECvR6iImBTZsgLQ26dVM6IoeEwN/+Bg1phM/OnTsJDg7GYDCoHaVOXbp0CQAfHx+Vk9S9ESNGAJCQkKByEidSXAzvv68ccPzzn2HsWDh4ULkD0qqV2umEcHohISG1ersQzkgKKSFU1qsXfP01ZGTA+PHKWfXbboMnn4R9+9ROVzsnT54kPz+fQYMGqR2lzmVnZwMNs5AaOXIkAOvXr1c3iDMoLYV586BdO6VTzIgRkJ6uDI+7zpw5IcT1zZ0794arTjqdjrlz5zo4kRC1J4WUEE4iLEzpmHziBLz+OqxYAZ07w+jRsG6d0lZda+bNmwfAxIkTVU5S9y5fvgyAbwOcC+Tm5oanpye7d+9WO4p6zGb45BNo3x6efRYGDVLubCxcCG3bqp1OCM05deoUa9asITw8HJPJhNFoxGQyER4ezpo1azh16pTaEYWosYa110aIBqBJE5g2DV5+Gb74Qjm7fued0L27Mo7mwQeV+Z5asGzZsopZSw2NfUXKz89P5ST1IyQkhCNHjqgdw/EsFvj8c6WBxOHDyhfcihXQqZPayYTQNPuwXWlxLhoSWZESwkmZTBAXp3T4W71aGUczfrxyM3zWLMjNVTth1fbt20ebNm3Q32DSvZbl5OQADbeQ6t69OyUlJRVnwRq8sjLlzkVkpPKF16WL8sX39ddSRAkhhLiuhnd1I0QDo9PB8OFKMbVrl7I69Ze/QOvWyiwqZx318/3331NaWtogzxDBlUKqSZMmKiepH3fffTcAS5YsUTlJPbNa4ZtvlMLp0UeVs1Dbt8P33yuvE0IIIW5ACikhNKRLF+WIRmamMvfzk0+UJmLjx8OOHWqnq+z1118HYNq0aeoGqSe55UuCzZo1UzlJ/YiOjgbghx9+UDlJPbHZlCFu3brBQw8pnfe2bIGkJGUfrRBCCFEFKaSE0KCWLeHtt5XVqFmzICUFoqKU1aoVK5Sb7Go7cOAAAA888IDKSepHXl4e0HALKT8/P9zc3NjhbBV6bdlssHw59OihDHFr1kyZQ7B6NfTurXY6IYQQGiKFlBAa5uUFU6fCoUPKUY78fBg1Sjnm8cknyugbNVgsFsxmM25ubuoEcID8/HwAPDw8VE5Sf4KDgznhrHtHa8pmgzVroG9fpRWmh4fSDvOnn2DAALXTCSGE0CAppIRoAAwGpblYWhps3Kh0bJ48WZlH9eabcOGCY/O89957APTo0cOxT+xABQUFAA2ykYZdly5dKCgooLCwUO0otbNuHQwcqMyA0umUgmrjRhgyRO1kQgghNKzhXgEI0QjpdMr1YmIiHDgA990Hb70FISEwZYoy9NcRPvjgAwDeffddxzyhCgoLCxt0EQVUtK1fuXKlykluUXIyDB2q7HktLlb2vaakwN13K18sQgghRC007KsAIRqx9u1h7lw4flyZS/XNN8rr7rsPNm+u3+e2bwfr2bNn/T6RihpDIRUbGxjwSiAAAAteSURBVAvA6tWrVU5SQ1u2KK0uBw6Ey5dhyRLYtg1GjpQCSgghRJ1p2FcBQgiaN4e//x2OHYN582D/fuVISN++8O23yviculRcXExZWRne3t51+8BOpqioCBcXF7Vj1KtWrVphNBrZtm2b2lGq5+eflUOCffvC6dPKf/AdO2DsWCmghBBC1DkppIRoJNzdlXNTe/cqHZ5dXZVzVe3bw/vvK40q6sIrr7wCXNkW1lA1hkIKICgoiKNHj6od4+Z27oToaOjZEw4fhs8/V4au3X8/NPBVQyGEEOqRnzBCNDJ6vdK0bP16ZbdTr17w0kvKOaq//hXOnKnd43/55ZcAvP/++7UP68RKSkowGo1qx6h3nTp1IicnB4vFonaUa+3dCw88oMyC2rsXFi2CPXvgkUegERS5Qggh1CWFlBCNWI8e8MUXyk383/1OWZkKDYUnnlCuR2/F+fPn0el0hISE1GVUp9NYCqmoqCgAFi5cqG6Qq6WnK8XS7bfD9u1Kr/8DB+Dxx5UWlkIIIYQDSCElhOC22+Ddd5UBv9Onw9q1yjXqPffADz8oI3iq4+zZs9hstgY7pPZqpaWlmEwmtWPUu759+wLw4YcfqpwEpe1kXBx06qR05Js7VymqJk6UAkoIIYTDSSElhKjg5wd/+hMcOQKLF0NWltIpuls35c+lpTf/+Oeffx640u2tITObzbi6uqodo96NHj0agIMHD6oXIjMTJk2Cjh2VKv+995Si6umnoREUs0IIIZyTzmar7r1mIURjY7PBTz/BrFmwciUEB8PUqfDUU0rRNXnyZMLDw3nttdcqnaHR6XS8/fbbZGRkMH/+fBU/g/rj5eVF69at2b9/v9pR6p29zbvVanXsE584Af/6l7J1r1kzpY//008rnVOEEEIIlcmKlBDihnQ6GDZMmWO6Z48ymudvf4PWrZUGFU2ahDFt2rRrGhHYbDamTZtGeHi4Ssnrn8Viwb2RXND7+Phgs9kc13Di9Gl4/nkID4fvvoO331aWSV98UYooIYQQTkNWpIQQNXL2rNKUYu5cuHzZCNz44tpgMGA2mx0XzoGMRiO9evVic31PN3YC/fr1IzU1lYULFxIXF1d/T5SVBTNnKv+53N2Vfaa//z008JlkQgghtElWpIQQNRIUpOy2OnECblZEAc7ZMruOWK3WRrMiNXHiRAA+/vjj+nmCCxfgz3+GsDBlG9+0aXD0KPzlL1JECSGEcFqyIiWEuGU6ne7m76CH4HeCCfQKJNAz8MqvV/++/NdmHs3Q67Rzb0en0zF27FiWLFmidpR6Z7FYMBqN+Pn5cfny5bp74EuXlAN4//d/yp9ffBFefhmaNKm75xBCCCHqifSLFULUHxtM7DaRrPwssgqyOHDhABsyN5BVkEWhubDSu7roXGju2ZwAz4Aqi67mns0x6NX79mVvuuDp6alaBkcyGAzodDpyc3Pr5gGzs2HOHJg9GywW5TzUH/8I/v518/hCCCGEA0ghJYS4ZQaD4abb9wwuBt4Y+sZ135Zfmk9WfhbnCs6RVZBVUWzZf83MziTtZBpZBVnkllS+gNehw9/Dv1orXQGeARhd6mZw7uTJk4mLi6Nr166A0rnPLjk5mfj4+AbbpdDDw4OCgoLaPUhenrL69J//QHExPPecsqUvMLBuQgohhBAOJIWUEOKWTZ8+nWnTpt307TfiZfLCq6kXbZu2rfJ5isxFnCs4x9n8s9ctuk7mnmT76e1kFWSRXZx9zcc3dW9aUVxVrHiV/znIK6hS4eVquPFsqLi4OMaMGcP//vc/ALzLz+8kJyczZswYkpKSqvxctKpt27bs3r2btWvXcvfdd9fsgwsK4IMP4J13lGLq6aeV808tWtRPWCGEEMIBpJASQtyyjIwMZsyYcc0cKYPBwPTp08nIyKiT53E3unOb323c5ndble9bYimpWOU6V3DumqIrqyCLvef2klWQxYXCC9d8vK+r741XuvwDmfHJDJ6IewJQ2oJfXUQNGDCgTj5fZ/TQQw+xe/duZs+eXf1CqqhI6cA3Y4aynW/SJHj1VWjVqn7DCiGEEA4gzSaEEI2WuczM+cLz1xZbVxVd9t9fKLyA1VY+kPb/t3f/rlHfcRzH38mdltOAaAhnB0FoUKyDdLaTQ4dOXQQpDlkylNL+AxGNKCGj4CjYv6N/gZt20EESaHFIySkq2rQh9nIdQtNGL8m9kJioj8eYb758vzd9efL5daeqfq5qfN6o3q+9OvXjqTrxxYn1qYT9RrpG9o9svznHHra8vFytVqvGxsaq0+ls989Vt26tnf/U6VRNTFRdulR1/Pi7eFUAeCeEFMAAuqvdevLnk1pcWqz7c/dr+ofpmvtlrs59d66OfXXsjRDr9rob7m81WwOt6WqPtOvQJ4f2ZHQNDQ1Vo9HYfF3cykrV7dtr++MvLFRdvFh1+XLVZ9tP3wSA942pfQADaAw31mJnpF0v5l7U498e140bN2p6erqufnu1vvzmv2l9q73VevrX0y1Huu7+fnf9b69WNx5avL+xf6Dgah9s15HWkR2PrsnJyRofH6+qqm63u/689SmcDx/WrbNnq65dq3r0qOrChaorV6pOntzR9wKA3WRECiDw+pqot10j1ev16vny8/7B1WeK4fLfyxvubw43B9oyvj3SrtHWaDWGG/E7zszM1NTU1KbXfzp8uCaePas6f34toE6fjp8BAO8bIQUwoM2i6V1tONHr9erlysuBgmvxj8VaerVxu/LhoeEaOzDWd4rh0ZGjG3Y1/P9ZXfv27dtym/tPh4Zq4d69qjNnduy3A8BeI6QABvTvOVL9YmkvniO1tLLUN7he30a+s9Tpe1bX6IHRah9s14PvH2z7LJ8SAD42QgqA9bO6+oXXza9vbnu/TwkAHxshBcCWBtnMwqcEgI/N8G6/AAB7W7O59Qav210HgA+RkAJgS9evX3+r6wDwIRJSAGxpfn6+Zmdn3xh5ajabNTs7W/Pz87v0ZgCwe6yRAgAACBmRAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgJCQAgAACAkpAACAkJACAAAICSkAAICQkAIAAAgJKQAAgNA/FInWzYXzUKMAAAAASUVORK5CYII=" + ], + "text/plain": [ + "PyPlot.Figure(PyObject )" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "(-0.5,4.0)" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function plot_element(el; plot_with_normal=false)\n", + " # create a array of vectors\n", + " xis = Vector[[xi] for xi in linspace(-1, 1)]\n", + " coords = interpolate(el, :Geometry, xis)\n", + " normals = interpolate(el, :Normals, Vector[[-1.0], [1.0]])\n", + " ncoords = interpolate(el, :Geometry, Vector[[-1.0], [1.0]])\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:2\n", + " p0 = ncoords[i]\n", + " p1 = ncoords[i]+0.1*normals[i]\n", + " plot([p0[1], p1[1]], [p0[2], p1[2]], \"-k\")\n", + " end\n", + " end\n", + "end\n", + "\n", + "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", + " if has_projection(sel, mel)\n", + " xi1, xi2 = calc_projection(sel, mel)\n", + " X1 = interpolate(sel, :Geometry, xi1)\n", + " X2 = interpolate(mel, :Geometry, xi2)\n", + " #println(\"$X1\\n\\n$X2\")\n", + " for s=1:2\n", + " plot([X1[s][1], X2[s][1]], [X1[s][2], X2[s][2]], \"-kx\")\n", + " end\n", + " end\n", + " end\n", + " end\n", + " axis(\"equal\")\n", + " axis(\"off\")\n", + " ylim(-0.0, 0.5)\n", + " xlim(-0.5, 4.0)\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": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": [ + "iVBORw0KGgoAAAANSUhEUgAAA1IAAAEbCAYAAADQ7rWKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XdYVGf6//H3DFMYiqBgV2woooi9Y0ui0QiKLWuSjW7MmjVls8nmu4nZuNkka7LpZZNf3DRXdFOtiIlGzcaCgL0iURGxKzbqMMwMM78/DhBNbAjMmQP367rONTgDcz6oMHOf53nuR+d2u90IIYQQQgghhLhperUDCCGEEEIIIYTWSCElhBBCCCGEEJUkhZQQQgghhBBCVJIUUkIIIYQQQghRSVJICSGEEEIIIUQlSSElhBBCCCGEEJUkhZQQQgghhBBCVJIUUkIIIYQQQghRSVJICSGEEEIIIUQlSSElhBBCCCGEEJUkhZQQQgghhBBCVJIUUkIIIYQQQghRSVJICSGEEEIIIUQlSSElhBBCCCGEEJUkhZQQQgghhBBCVJIUUkIIIYQQQghRSVJICSGEEEIIIUQlSSElhBBCVKPp06fz2muvYTQa0el0FYfRaOS1115j+vTpakcUQghRDaSQEkIIIapRmzZtmDlzJk6n84r7nU4nM2fOJDw8XKVkQgghqpPO7Xa71Q4hhBBC1BZGo/FXRdTlDAYDDofDg4mEEELUBIPaAYQQQojq5HaXUlpqxeWyXnZbjMtVjMtlu+wowe2243LZcbsdZR87cLuduN1OoBS3uxS32wW4Km5/TQfo0On0gP66RRRww8eFEEJogxRSQgghvIbL5cDpvIjDcbHs9hJOZ+4VR2lpHk5nPqWl+WW3hZSWFpTdFuJyFVfyrDp0OhN6vQmdzohOZyi79ak4lJnw5cWS7rKvdQNulMkd5cWWEEKIukCm9gkhhKgxLpcdu/0sdvtZHI6z2O05OBw5Zbfnyo7zZccFSksLrvo8Op0JgyG47AjCYAjCxycQH596GAyB+PgElP05AB8ff/R6f3x8/PHx8UOv90Ovt+DjY0Gv90WnM6PXK4dSQFXvNUWdTnfDz/n973/PpEmTGDZsGEajsVrPL4QQwjOkkBJCCFFpLpcTu/0MJSUnKCk5gd1+ipKSU9jtp7DbT5d9fAan8+KvvtZgqI/R2AiTqSFGY0OMxlCMxhAMhhCMxpCyjxtgNNbHYGiAwRCMj49Fhe/y1txojZROp6Nt27YcPnyY+vXrEx8fz6RJk7j99tsxmUweTCqEEKIqpJASQggNmz59OlOnTiUmJuZXjyUnJ5OQkMAnn3xS6ed1OvOw2bKx2Y5isx2lpOQYNtsxSkqOlRVPp7h8vZBOZ8Zsbo7J1BSzuRkmU1NMpiZlt40rDqOxIXp97S4WXnvtNWbOnHnNx1999VWefvppdu/ezcKFC1m4cCGHDh0iKCiIsWPHMnHiREaMGIHZbPZgaiGEEJUlhZQQQmhYcnIycXFxJCUlXVFMXev+ci5XCTbbUYqLs7DZDlNcfASbLQub7Qg2WzZOZ27F5+p0Jnx9wzCbwzCbW+Lr2xKzuSVmc4uyozkGQ4ObmtJWF0yfPp3w8HBmzZp1xciUwWBg9uzZZGZmXlHcut1u9u3bx8KFC1m0aBEZGRkEBgYSFxfHxIkTGTlyJBaLdkbkhBCirpBCSgghNO6XRVP5nxMTl9CrVzOs1oMUFx8qOzIpLs7EZjtG+YiSTmfE17cVvr5t8fVtja9vm7Lb1vj6tsJkalzWZEFU1tdff83kyZP59NNPefDBB2/qa9LT01m8eDGLFi1i7969+Pv7M3r0aCZOnMioUaMICAio4dRCCCFuhhRSQgihcQ7HRdau/Zpx4/5EeHh9jh27yBtvNCYi4gxQCoBeb8FiCcdiaV92267i1mxuUdaZTlS3M2fO0LRpU2bMmMGcOXMq/fUHDx6sKKp27NiBr68vo0aNYsKECcTGxhIUFFQDqYUQQtwMKaSEEEIj7PbzWK3pFBUph9W6n6KiDByOswDccw+cOQNPP92R6dNHYrFE4OcXgcXSHrO5mYwqqcTHx4eBAweyYcOGKj1PVlYWS5YsYdGiRWzevBmTycTw4cOZMGECY8aMISQkpJoSCyGEuBlSSAkhhJcpLbVhtaZTWLiHoqK9FBXtpbBwb0XBpNMZsFg64O/fGT+/Tvj7R7JrVzFjxsygpKSEwMBAvvvuu6uujRKeV79+fQIDAzl27Fi1Pefx48criqpNmzah1+sZNmwYEyZMYNy4cTRu3LjaziWEEOLqpJASQohbVNmmAldTUnKGwsJdFUdR0R6s1gMo65d0+Pq2JSCgC/7+XfD374y/fxQWS/srOt+Vr4l66aWXePzxx5k2bRpLliy5ZqMJ4VmRkZEcPXoUq9VaI89/+vRpli1bxqJFi1i/fj0ul4uYmBgmTJjA+PHjadmyZY2cVwgh6joppIQQ4ha98sorPPfcc9d8/NVXX+WZZ54BlM5sJSXHKCjYQWHhjrLbndjtpwHw8QkkIKAr/v5dy2674O8fhcFw/cYClzeaGDBgAD4+PgwfPpznn3/+ul37hOfExsby7bffUlpail5fs9Mrz58/z/Lly1m8eDFr1qzB4XDQp08fxo8fz4QJEwgPD6/R8wshRF0ihZQQQtyiG2286u+vZ8+eZyko2EZBwTaczgtlX9eYwMAeBAR0JyCgO4GB3fH1bXNLa5h+uY9U/fr18ff358SJE1XaR0pUn+eee45XXnmFnTt30q1bN4+dNy8vjxUrVrB48WJWrVpFcXEx0dHRFUVV586dpWW9EEJUgRRSQghxi270JtRggPXrmxAY2IuAgJ4EBiqH2dysxjL16tWL3bt343A4auwconK+++47Ro8ezbvvvsuf/vQnVTIUFRWxatWqiimfBQUFdOjQgXHjxjF+/Hh69+59xf/nmtroWQghahMppIQQ4gZcJS6KjxRTnFlM8aGyI7OYbmtuPLrgcrk8etX/4Ycf5t///jfZ2dm0atXKY+cV11ZYWEhgYCBTpkwhISFB7TiUlJTwww8/sGTJEpYtW8aFCxdo2bJlRVEVExNDamrqLW30LIQQdYkUUkIIAbgcLmxHbBQfKsZ6yPpzwXSoGNsxW/netegteizhFizhFros7XLD5/X0r9jExETi4+N55513eOKJJzx6bnFtRqOR7t27s2XLFrWjXMHpdJKcnMySJUtYsmQJJ0+epGHDhowdO5aIiAhefvnlX230LEWUEEIopJASQmhGVacbuUvd2I4rxVLxwWKsBy8rmI4Ul+9di97352LJ0qHstr0Fv/Z+mJqa0OmVEaYbrZEyGAwen2Jns9mwWCyMGTOGxMREj55bXFvDhg0xGAycPn1a7SjX5HK52LZtG0uWLGHx4sVkZmZisViw2+2Eh4dz5swZVqxYIUWUEEKUkUJKCKEZ17oifvn9AwcOxJHjUIqk8mKp/PZwMe4S5VeezqjD0k4pkMoPvw5+WNpbMDc3VxRL1/Paa68xc+bMaz5+edc+TwoICKBRo0ZkZWV5/Nzi6rp160ZGRgYlJSVqR7kpbreb9PR0lixZwssvv4zdbsfHx4e77rqLcePGERcXR2hoqNoxhRBCVVJICSE0pbxo+udL/2R0l9HsWLeDKa9N4YOYD+iU2wnrQSul+WVDSzrwbeWrFEkRSpFk6aCMLJlbmdEbqtaKujr2kaoJUVFRZGZmYrPZPH5ucXWTJk1i0aJFlJSUYDKZbvwFXiI5OZnbb78du92O2Wymffv2pKeno9frGTRoEOPHjyc+Pl72qhLiFzIyMhgzZgzHjh3D7Xaj0+kICwtj+fLlREZGqh1PVBMppIQQXsvldFFytATrASvWA2UjSwespO5J5ZELj2DBgg8+vBb0Gv2i+lUUSxVFU7gFH18fj2Q1mUz06NGDtLQ0j5zveqZMmcKCBQs4d+6cjBp4iX/+85/89a9/Zf369QwePFjtODel/KLFtGnTePvtt5k9ezZvvvkm8+bN4+zZsyxdupQffvgBh8NBr169GDduHOPGjZM3iaLOS09Pp1u3bled+m0wGNizZ4/8nNQSNbszoBBC3ATHJQd5qXmcnnearGez2DduH1s6bWGj30Y2h29m7+i9HHnuCHkpeRhDjQycNhCAYop58ZUXeST3EXok96DjZx1pNbMVDcc1JCAqwGNFFEBpaSkWi8Vj57ueESNGAMgaKS8yaNAgADZu3Khykptz+XTZ++67D1C6/SUlJfG73/2OTp06sXLlSs6dO8fnn39O69ateeWVV+jUqRMdO3bk2WefZcuWLbhcLpW/EyE8Lz4+/prrZ51OJ2PGjPFwIlFTZERKCOERLqfSFc96wIr1JyvFB4qx/qSMNDnO/dyQwRxmxi/CTxlVirBUfGxu8fO6pSeeeIL33nuPKVOmsHz5cq/oIqbX6xk1ahTffvutqjkA8vPzCQoKYvLkyXz55ZdqxxFQMTVuwoQJLFq0SO04N3R5Yxen04nRaKzIfq3GLsXFxaxdu5alS5eyfPlyLly4QPPmzRk7dizx8fEMGTJEU9MahbhVZrMZu91+zcdNJpNm1kuK65NCSghRrZz5zopiyXrAijWjrHA6VIzbofy60fvrKwokv44/31raW/Dxu/4oUnJyMkOHDqW0tBSHw0FaWppXtGTW6XRMmjSJb775RrUMl7NYLLRu3ZqMjAy1o4gyvr6+REREsHv3brWjVJrJZCI6Oppt27bd1Oc7nU42bdrEsmXLWLp0KUePHiUoKIjY2Fji4+MZOXIkAQEBNZxaCM86fRpWr4YHHjDhdl+7Y6vRaLxuoSW0QwopIUSlud1u7KfsWH+yUpRRpBRNPylFk/3Uzy8O5hZmZVSpvFiKvGx06RY2qS2fblRSUkKDBg04ceLEFferVUyVtxx/4IEHmDt3rsfPfzXt27fn1KlTFBUVqR1FlGnWrBl2u53z58+rHaXSQkJC8PX15eTJk5X+Wrfbze7du1m6dCnLli1jz549mM1mhg8fTnx8PHFxcTRq1KgGUgtRs6xW2LAB1qxRjr17lft1OjNut4xI1QUGtQMIIbyXy+nCdtj2c7GU8fMIU2mB0hlPZ9QpDR4i/WjyQBP8I/0rpuQZAqv3V0xCQgJz585l/PjxVyzYj4mJISkpiYSEBFUKqYsXLwIQGBjo8XNfS3R0NJmZmVitVvz8/NSOI4CWLVuyY8cOtWPckkaNGnH06NFb+lqdTke3bt3o1q0bL774IllZWSxbtozExEQeeugh3G43AwcOrJgCGB4eXs3phageTids2wY//ABr10JKCtjt0KwZjBgBzz4Lt98OAweGkZmZec3nCQsL82BqUZNkREoIQWlx6c/T8MqOoowiZTqeXfkV4VPPB79IP/wj/ZWRpY7KCJNvG98qtxGvjBdeeIEXX3yRtWvXcvvtt3vsvNezb98+unTpwqxZs/jHP/6hdhwAPvzwQx599FEWLlzIxIkT1Y4jgKlTpzJ//nzy8vKoV6+e2nEqZcSIEaxdu7bam0ecO3eOpKQkEhMTWb16NTabjU6dOhEfH098fDw9e/ZEr5e+WEIdLhfs2QM//gj/+x+sXw8FBRAYCEOHKkXT8OEQGQmXT7LIyMggOjpauvbVAVJICVGHOAucSpG0vwjrfmvFx7YjNij7TWBqYlIKpbLDv5NSOJmamG5pOl5169mzJ7t378Zut3vNG6wNGzYwZMgQXn/9df7yl7+oHQeAU6dO0bx5cx588EE+/fRTteMI4P333+fxxx8nKSmJ2NhYteNUyh//+Ec++OADjh8/TosWLWrkHEVFRXz//fckJiaSlJTEpUuXaN68OWPGjGHs2LEMGzZMmlWIGlVaqhRO69crx4YNcPEi+PrCwIFw220wbBj07g2GG0y4kH2k6gaZ2idELeS45LiiYCq/LTn+85xs39a++HXyI3RcaEWx5BfphzHYqGLyG9u/fz9t27b1miIK4NKlSwAEBQWpnORnzZo1w2g03nRzAFHzhgwZAkBKSormCqkuXboAkJaWVmMjnP7+/owfP57x48fjdDpJTk4mMTGRZcuWMWfOHOrVq8eoUaMYO3Yso0aNIjg4uEZyiLrDboft25WCaeNGSE6GvDwwm6F/f/jjH5XCqV8/5b7KiIyM5NChQzUTXHgNKaSE8KCEhASaN2/Oww8//KurVHPmzOHkyZNMnTr1pp/PccFB0f4iitIvK5jSrdjPlC1y1YOlnQW/Tn40/m1jZYSpsz9+EX74+Htuj6Xqkp6ejs1m44477lA7yhVyc3MBvO6NXZMmTThy5IjaMUSZqKgoAE127evTpw8AO3fu9MhUUYPBwNChQxk6dChvv/02e/fuJTExkcTERO69996Kx8tHq2TNibgZly5BWhps2qQUTZs3g80G/v5K4fR//weDB0OfPsoolBA3IoWUEB7UtGlTRowYwS9n1GZmZjJixAhWr1591a9zXHRQlF5WMKVbK4onx9my9qo+YAm34N/Zn6a/b4pfZ7+Kpg+e3JS2pn300UcAPPTQQyonuVJ5IVW/fn2Vk1wpKiqKlStX4nQ6MdxoHoqocXq9HovFQlZWltpRKq28CPzpp588fm6dTkd0dDTR0dH87W9/4/jx4yxfvpzExESeeuopHn/8cbp168bYsWMZO3Ys3bp184ppyKJ6JSQkMGTIEFq3bv2rx7Kzs1m/fv0VFyJdLjh4EFJTlSMlBdLTlccaNoRBg+Dll5Xbbt3A6N2TMYSXkjVSQnhQ+/btr9vJp22btuxYsKOiaCovnMpHmHSGsg55nZSRJf9O/vh19sOvgx96k/dMdaspHTt2JDs7G5vNpnaUK8yePZu//e1vbN26lV69eqkdp8Krr77Ks88+y6pVq7jzzjvVjiOAVq1akZeXV1F8a4nJZKJLly5s375d7SgV8vLyWLlyJcuXL+e7774jLy+Pli1bVoxUySbAtUd2djbTpk1j7ty5VxRT5fe/+eZczpxpzebNykjTli3KCJROB507KyNOAwcqR7t2VzaHEOJWSSElhAfdaLdzAwbWsOaKEabyoy4VTNdiNBrp3Lkzu3btUjvKFZ5++mneeOMNsrKyaNOmjdpxKhw4cICOHTvy+OOP895776kdR6Csk0pOTqa0tFTtKJUWGhqKyWTi1KlTake5KrvdzoYNGypGq44dO1axrmrMmDGMGjXK60aNReWUF03PPjsbvb4PP/xwgs8+m4bJNJcTJ1oDEBqqTM3r319Z29S7N3jR8lVRy8hcDyE86IbXLXyg1/Ze+HX0Q2+uuwXT1axbtw6n08no0aPVjvIrBQUFgLJpqTeJiIjAx8eHtLQ0taOIMp06dWLDhg2cOnWKZs2aqR2nUho1akR2drbaMa7JZDJxxx13cMcdd/Dee++xZ8+einVV9913HwaDgUGDBjF27Fji4uJo27at2pHFTcjNhZ07YccO2L69NdnZcxkxIhKIQ68/T48ecxk0qDV9+ihFU9u2MtokPEfeqQnhQTeat6/30RPQNUCKqKuYO3cuAH/4wx9UTvJr5YVUQECAykl+rWHDhtedTio8q7xpw/r161VOUnlhYWFeN632WnQ6HV27duX5559n+/btHD9+nPfffx+z2czTTz9Nu3bt6NKlC8899xxpaWnVvj+WqDy3G7KzITERXnoJxo+HNm2gfn2l7fjzzyuPDxvmB9iAhaxZ83e2bm3N22/D5MkyZU94nrxbE8KDbtRZSjpPXdu6desICAjwyr+joqIiAK9qyV6uY8eOXLp0Sd4oeomhQ4cCaHKUsEOHDrjdbo4ePap2lEpr0aIFM2bMYOXKlZw/f55FixbRo0cPPvroI/r370+zZs148MEHSUxMrPh5FjXn0iWl3fiHH8LDD0NMDAQHK4VTfDz861+Qnw8TJ8J//wv79yt/TkmBvn2XAfD1118ze/aLXj1KKmo/mdonhAfNmTPnql37QLmCOmfOHBVSeT+n08mJEycYMGCA2lGuqqioyGu7hA0YMIB169axdetW+vbtq3acOq9NmzbodDr27dundpRKi46OBmDz5s20atVK5TS3LjAwkAkTJjBhwgRKS0tJSUkhKSmJ5cuXM3fuXHx9fbn99tsZM2YMsbGxmpuC6U0uXICMDOXYv1/pmpeeDuXL7AwGiIiALl0gNla57dYNmjW7+shSdnY2L7/8Mo0aNWLSpEn06dPnqg0ohPAU77t8KkQtdvLkSZYtW4ZOp8PHxwej0YjJZCI8PJzVq1dz8uRJtSN6pSVLluB2uxk/frzaUa7KarXi4+Odbebj4+MBWL58ucpJRLmAgABNXkUvn5bobc1eqsLHx4dBgwbx+uuv89NPP3HgwAFmz55NYWEhjzzyCM2bN6dXr168+OKL7Nix48brXOsgux0OHICkJHjzTZg+XdmLqWFDpfHDoEEwYwZ8952yX9MDD8CXX8KePVBUBPv2KX+eORNGj4bmza9dRE2bNo3GjRszePBgdDodrVu3Zu7cuUybNk2TP1NC+6RrnxAeNm/evIpf+t44Tc0bjR07luXLl3Pp0iWv2/QWoGfPnhWbBXsbl8uFwWBg8ODBrFu3Tu04AmUbhNOnT1NYWKh2lEpxuVz4+Pgwbtw4lixZonacGnfx4kVWrlxJUlISK1euJD8/n+bNmxMbG0tcXBy33XYbFotF7ZgeYbPBkSNw+LByZGYqx6FDyrql8iaU/v7QoQN07KiMNEVGKkf79lXf4DYhIYH+/fvTtWtXXnnlFZ588smKx662j5QQniCFlBAeNmzYMPR6PT/88IPaUTSjUaNGlJaWcuHCBbWjXFWnTp04evSo166tCA0NxWg0cvr0abWjCODOO+9kzZo1OJ1Or1xXdz1ms5lOnTqxc+dOtaN4lMPhYOPGjSQlJZGUlMThw4exWCzccccdxMXFERsbS9OmTdWOectsNjh+HI4eVQqj8tsjR5Tj8o73vr5KU4fwcOW2QwelUIqIuPaUvOqydetW+vTpQ2pqKv369au5Ewlxk2SNlBAedPToUdatW8e8efPUjqIZVquVc+fOMXLkSLWjXFNJSQkGg/f+Og0PD2fbtm1qxxBloqKiWL16NYcOHSIiIkLtOJVSr149zpw5o3YMjzMajdx2223cdtttvP322/z0008kJSWxYsUKZsyYgcvlolevXsTGxhIbG0uPHj0q1k1Onz6d8PBwZs2ahdPprHhOg8HA7NmzyczM5JNPPqmx7EVFSiF06hScPAknTijH8eM/Hzk5P3++TqcURG3aKMfttyu3bdsqhVPTpqBW/Z+SkoLZbKZ79+7qBBDiF7z3lV+IWmjBggX4+fkxYcIEtaNoRkJCAgD33HOPykmuraSkBKPRqHaMa+rbty+bN28mIyODyMhItePUef379weUFuhaK6QaNWpEVlaW2jFUpdPpiIyMJDIykqeffpoLFy6wcuVKVqxYwTvvvMMLL7xAs2bNGD16NHFxcbRo0YKZM2f+6nmcTiczZ87k1VdfrdT53W4oLIRz534+zp5ViqGzZ+HMGeU4fVo5ynZnqBAUpKxDatlSaewwZgyEhSlHq1bK/SZTVf6Gak5qaio9e/bEbDarHUUIQKb2CeExbrebiIgI+vXrx/z589WOoxnDhg1j/fr12Gw2TF766t6oUSMMBgOnLp//4kXWrFnDiBEjeOWVV3j22WfVjlPnnT9/noYNG/L73/++RkciasKoUaNYtWqVNF24BofDQXJyMitWrCApKYlDhw7d8Gt8fPzYsKGIvDxl89ny49Il5bh4Uel+V357/rzS4OGXQkKgceOfj6ZNfz6aN1dGmVq0UNYxaVVYWBi/+c1veOONN9SOIgQghZQQHpOWlkb//v1Zu3Ytt99+u9pxNCMoKAh/f3+vLVIAgoODqV+/PkeOHFE7ylU5nU6MRiMjR45k5cqVascRKNO6+vbty6ZNm9SOUilPPvkk7777LtnZ2Zpuge4pBw8evIlRRxNQUvEnHx9lT6XgYGjQQNmQNiREORo0ULrhhYQot5cfXjwoXi1OnDhBy5YtWbx4sdd2cBV1j0ztE8JDEhISaNGiRcWGnOLGcnJyyM/P584771Q7ynU5nU58q9qSqgYZDAbq1atHenq62lFEmXr16mlyY9vyvaTS0tKkkLoJHTp0uInPsrNvnzLlLigIAgJqtmGDVqWmpgI/T40Vwhtoq12QEBpVUlLCV199xf333++1+w15o48//hiABx54QOUk1+fthRQoG8HWxSYB3qpp06Ze24Xyesr3ktqxY4fKSWqXzp2VaXeBgVJEXUtKSgpt2rTRdHdEUftIISWEByQlJZGbm8v999+vdhRNSUxMRK/Xe/2IVGlpqdfvJ9O7d28cDgcnTpxQO4pA2UvKZrNd0cVNC8qblRw4cEDlJNpxo46e3tzx05ukpqbKaJTwOlJICeEB8+fPp3fv3tIxrZLS09Np3bq11++143K58PPzUzvGdZUXo8uWLVM5iYCfp8hpbT8mvV6PyWQiOztb7SiaMXv27Co9LsBms7Fjxw4GDBigdhQhruDd706EqAVycnJYuXKl7LheSQcOHKC4uFgTjTlcLhf+Xt4K66677gLgxx9/VDmJABg4cCAAGzZsUDlJ5QUFBXH27Fm1Y2hGZmYmr7766q9GngwGA6+++iqZmZkqJdOO7du343A4ZERKeB0ZTxaihn355ZfodDomT56sdhRN+eijjwB46KGHVE5yc7y9kPLz88PPz4/du3erHUUAgwYNArS51qhx48by5r8SylvcP/roowQGBmqy7b3aUlJS8PPzqxjJFcJbyIiUEDUsISGB2NhYQkJC1I6iKatWrcJkMtGrVy+1o1yX1WoFICAgQOUkNxYWFsbJkyfVjiFQCluj0ajJtUZhYWHYbDZcLpfaUTTl9OnTADRo0EDlJNqTmppKnz59ZD2Z8DpSSAlRg/bu3cvOnTuZMmWK2lE059ChQzex/4r6zp8/D0BgYKDKSW6se/fu2Gw2cnNz1Y4iUPYfO378uNoxKq3859Jb903zVuWFVGhoqMpJtMXtdpOSkiLro4RXkkJKeMTZs2cZOHAg+/fvVzuKR82fP5+QkJCK9SlKI1LnAAAgAElEQVTi5mzatAmn08moUaPUjnJDFy9eBJR9gbxd+Xqz5cuXq5xEALRo0UKTRW35RYMPP/xQ5STaUr79QKNGjVROoi3Z2dmcPXtW1kcJrySFlPCIDz74gN27d9ep/R+cTif//e9/ueeeezCZTGrH0ZRPP/0UgBkzZqic5MbKC6mgoCCVk9zYuHHjAFi7dq3KSQQoIzt2ux2bzaZ2lErp3bs3ABs3blQ5ibbk5OQAUkhVVkpKCgD9+vVTOYkQvyaFlKhxVquVDz/8kAcffJD69eurHcdj1q5dy5kzZ6Rb3y348ccf8fPzo02bNmpHuaFLly4B2iikGjRogNlsZvv27WpHEShTLUEZgdWSkSNHAnDs2DGVk2hL+TTgunRBsTqkpKQQEREhUyKFV5JCStS4efPmkZubyxNPPKF2FI+aP38+kZGR9OzZU+0omuJyuTh+/LhmujOVT80KDg5WOcnNad68OUePHlU7huDnzn3JyckqJ6mc8gX/eXl5KifRlgsXLgDQrFkzlZNoi2zEK7yZFFKiRpWWlvLOO+8wceJETYwuVJe8vDyWLl3K1KlT0el0asfRlMTERFwuV8U0NG9X/mZSK6OtXbp0oaioSHPTyWqj8ilyWmxJr9PpKCkpUTuGppRPA5aRlZtXWFjI7t27pdGE8FpSSIkatXz5cjIzM3nqqafUjuJRixYtoqSkhPvuu0/tKJqzYMECAH7/+9+rnOTm5OfnA9oZkRo6dCgA33//vbpBBAaDAV9fXw4dOqR2lEozmUy43W61Y2hKbm4uOp0OvV7eet2srVu34nK5ZERKeC35aRY16q233mLQoEH06dNH7SgeNX/+fO644w5atGihdhTNSUlJITg4WDN7rZQXUlrJWz7St3LlSpWTCFD+35w6dUrtGJVW3qXS6XSqnEQ78vPzpYiqpJSUFOrVq0enTp3UjiLEVclPtKgxaWlpbNq0qc6NRh05coQNGzbI3lG3wGazcfbs2YopT1pQUFAAaGe6TqtWrTAYDGzdulXtKAJlc1strjVq2bIloDSGETenoKBANpStpNTUVPr16ycFqPBa8j9T1Ji33nqL9u3bExcXp3YUj1qwYAEBAQGaWePjTebPnw/Ab37zG5WT3LzCwkIA/Pz8VE5y8xo3bkxWVpbaMQQQGRlJaWlpxfoZrSifZbBkyRKVk2hHUVGRbIVRCS6Xi9TUVFkfJbyaFFKiRmRlZbFkyRL+/Oc/16krSW63m/nz5zNx4kT8/f3VjqM533zzDYCm1paVF1Ja0qlTJ3Jzc3G5XGpHqfPKR1/Xr1+vcpLKGT9+PICMbFaC1WrFbDarHUMzDh48yMWLF6WQEl6t7rzDFR717rvv0qBBgzo3vS0lJYXDhw/Xue+7umzfvp0mTZrg6+urdpSbZrVaNXexYPDgwQBs2LBB5SSi/N8iNTVV5SSVM2zYMED2kqoMm82mqd9taktNTUWn09G3b1+1owhxTdp69ReacPHiRT777DMeeeQRTU13qg4JCQmEhYUxZMgQtaNozvnz58nNzdXc1UctFlLx8fEAJCUlqZxEREZGArBnzx6Vk1RO+Vqf8mYr4sZKSkrq3GtiVaSkpBAVFVXR2EQIb6StV3+hCR999BGlpaU8+uijakfxqOLiYr755hvuv/9+zb2x9gaffvopAFOnTlU5SeUUFxfj4+OjdoxKiYqKQq/Xk5aWpnaUOk+v1+Pv76/JNWs6nQ673a52DM1wOBwEBASoHUMzZCNeoQXybk9UWUZGBu3bt8dsNmM0GvnrX/+K2Wyu2MW9rli+fDl5eXkyre8WLV26FL1eT2xsrNpRKsVms2myE1doaCgHDx5UO4ZA+bc4e/as2jEqzWw2y15SleB0OgkMDFQ7hibk5uaSnp6uuRkKou6RQkpUSXp6OtHR0WRmZmK32yv2FMnPzyc6OpqMjAyVE3rO/Pnz6devHx06dFA7iibt27ePsLAwzY3m2Ww2jEaj2jEqrUOHDprrFFdbtWnTpqKNvhZMnz6d1157DZvNBigjUzqdDqPRyGuvvcb06dNVTuidXC4XQUFBasfQhM2bNwPIiJTwetp6xyK8Tnx8/DU3ZHQ6nYwZM8bDidRx5swZvv/+e81NS/MWhw8fxmq1Vixg1xK73a7JQmrAgAG4XC527dqldpQ6LyoqCrfbrZnGDW3atGHmzJm/ut/pdDJz5kzCw8NVSOXdyl8ng4ODVU6iDSkpKYSEhNC+fXu1owhxXVJIiSq50Qu/Vt4YVNUXX3yBj48Pd999t9pRNOmjjz4C0OSVbLvdrsm9Ycr3d1u2bJnKSUR5V7J169apG+Qm/f3vf7/u47NmzfJQEu04c+YMAA0aNFA5iTaUr4/S6XRqRxHiuqSQElVyo/nxdWX+fEJCAmPGjJEXyVu0cuVKjEajJqdxOBwOTe4NM2DAAHQ6HcnJyWpHqfOGDh0K/DydydtdaxZCxeOlpXx/8SKb8/M5aLVyzm7HWcf3LDt9+jQAISEhKifxfqWlpaSlpcn6KKEJ2lshLbzKja4WORw6Fi6EiROhtl5Y2r17N3v27GH27NlqR9GsAwcOaHYKh9Pp1OTeMHq9nuDgYPbv3692lDqvRYsW6PV60tPT1Y5SPfR6Rl6lnXuQjw+hRmPF0dBkopHRSGOTicYmE00uOxoYDLVqNKK8kGrYsKHKSbxfeno6BQUFUkgJTZBCSlRJWFgYmZmZ13zczy+Mu++Gvn3h9dehbO/JWiUhIYGGDRsycuRItaNo0tatW3E4HJr9+3M6nVgsFrVj3JJ27drJGikvERAQQHZ2ttoxri8tDV599cafV1rK0X79yHU6ueR0csnh4KLTyQWHg/MOBxccDs45HBywWtlot3PW4aCwtPSKpzDrdDQzm2luNtPisiPMbKal2UwrX19CjUbNFFs5OTkANGrUSOUk3i81NRUfHx969eqldhQhbkgKKVEly5cvJzo6+qpTPQwGA9u2LefMGfjLX2DIEIiNVV6HO3dWIWwNcDqdfP7559x7772abDjgDT7++GMAZsyYoXKSW1NaWqrZQqpPnz5s27aNw4cP065dO7Xj1GmNGzfm5MmTasf4Nbcb1q6Ff/4TfvwRIiIw6PXXnapnMBgI8/UlrBKnKSot5azdzhm7nVMlJZy22zllt3OipISTJSVsLyjgeEkJtsvO66fX08rXlzblh8VCO19f2lkstLVY8Pei/d3OnTsHQJMmTVRO4v1SUlLo1q0b/v7+akcR4oakkBJVEhkZyZ49exgzZgzHjh3D7Xaj0+kICwtj+fLlREZGEhkJW7bAN9/Ac89BdDRMnQovvggtW6r9HVTN6tWrycnJkW59VfC///0Pi8Wi2al9LpcLPz8/tWPcklGjRvHhhx+ybNkynnrqKbXj1Gnt2rXj0KFDuFwu79gCwOWCxER45RXYtg169oRFiyA+ntlvvnnVrn3lbmWas7+PD23LCqBrcbvdnHc4OFZSwjGbjezLjvV5ecw7c4aiywqtpiYTHSwW2vv50cFiIcLPjwg/P9r6+mL08N/x+fPnAWjWrJlHz6tFqamp3HnnnWrHEOKmSCElqiwyMpJDhw5d93P0epg8GcaPh48/hpdegi++gMceg2efBa2uv01ISCAqKopu3bqpHUWTXC4X2dnZ9OzZU+0ot8ztdmv2yumIESMA+PHHH6WQUll0dDSrVq1i//79REVFqRfE4YAvv1SmDmRkKFMJvv8ehg+vWOiamZnJq6++yqxZs66YjWAwGJg9e/Z1p3tXhU6no6HJREOTiZ5X2djW7XaT43BwuLiYzLLjoNXK9oICvsrJqZg+aNTpaG+xEOnnRyd/fzr7+9PZz48Ofn6YaqjAKt+zrWnTpjXy/LXFuXPnOHToEC+++KLaUYS4KVJICY8ymZTiaepUePttePNN+OQTePpp+NOfICBA7YQ3Lzc3l8TERGbPnq2Zefre5rvvvsPlcjF27Fi1o9wSV9nVb62OSJlMJgICAti3b5/aUeq88oX1GzZsUKeQstngP/9RFrNmZyvzsD/9FK6y4P+TTz4B4JlnnvFwyOvT6XQVjSsG/GLjW7fbzWm7nZ+sVn6yWtlfVESG1conp09zxm4HwKDTEWGx0CUggGh/f6IDAujq709zs7nKv+Nzc3MBNLlVgielpaUBshGv0A4ppIQqAgPh73+HRx5RZo689BK8/74y9e+hh0AL3aS/+eYbHA4H9957r9pRNCshIQHQ5v5RAPn5+QAEXuXquFa0adOGAwcOqB2jzhsyZAigNF/xqIIC+Pe/lStbOTlw993KlL7oaM/mqGG6suYVzcxmbqtf/4rHLjocpBcVsa+oiL1lx3cXLpBfNoIVYjDQLSCAHoGB9Ci7DbdY0N+guJo+fTpTp04lJiaGvLy8K6ZsJicnk5CQUFGUCkVKSgpNmzalVatWakcR4qZIISVU1bAhvPMOPPGEUlg98YTyev7CC/Db34IXrRX+lfnz5zN8+HCZ814FmzZtIigoSLOdrC5cuABou5Dq0aMHe/fuJScnR7P/DrVBcHAwPj4+ZGRkeOaEFy8qV6/eew8KC2HKFHjmGdDoWsWqaGA0Mig4mEHBwRX3ud1ujpWUsKuwsOL4OieHN44fByDQx4ceAQH0Cgykd7169A4MpI2v7xUjV1OnTiUuLo6kpCTy8/MxGJS3XMnJyRX3iyulpKRU7HEnhBZ4wYpWIaBVK5g3D/buhe7d4Xe/gy5dYPFipWmUt8nMzGTTpk3SZKIK7HY7p0+fpkePHmpHuWXl6x7q1auncpJbV76oe9myZSonEUFBQRw7dqxmT3L2rFIwtWqlrIP67W/h8GFlGl8dLKKuRafT0crXl7Ghofy9dWuWRkVxtH9/zg8cyJroaP4aFkZDk4lF584xef9+2m3eTMNNm7hrzx5eys7m+4sXierbl6SkJOLi4sjJycFoNF5RRMXExKj9bXoVh8PB1q1bZVqf0BQZkRJepVMnWLIEtm6FWbOUjXx79IDZs2HkSO/Z1HfBggUEBgZqdm2PN/j8888B+M1vfqNykltXXkgF/WI9hpbExcUBsHbtWh566CGV09RtzZo1u2Hjnlt24gS88YbS7cdohEcfhSefhMaNa+Z8tVSI0cgdDRpwR4MGFfeds9vZWlDA5vx8thQU8O6JE1xyOtEBnfz8iPngA769/37Mvr5SRF3H7t27KS4ulo14haZIISW8Uu/eSqOo9euVguquu5Q1z7Nnw7Bh6mZzuVzMnz+fu+++W7NNBrzBV199BcD999+vcpJbV76AXMuFVEBAABaLhd27d6sdpc5r3749+/btw+l0VkwDq7LDh+G115Qh/4AApU3qH/8Iv1gnJG5dQ5OJu0JCuKus/azL7eZQcTGpeXmk5ueT0qEDbrMZW3ExhunTeSc4mG3HjzM4OJiuAQH4eMsVQpWlpqZiMpk0PUtB1D0ytU94tSFDYMMGWLkS7Ha47TblSElRL1NycjLZ2dlMmTJFvRC1wNatW2nUqJGmi9HyQir4srUVWtSiRQuOl639EOop30Zh8+bNVX+yjAxl3VNEhNI84uWX4ehReP55KaJqmF6nI8LPj981bcpHERHMKSkhwGhk3MMP45g3j8wtW5iZlUXP7dtpkJzM6D17eOPYMbbm51PqjXPZPSQlJYWePXti1kK3KSHKSCElvJ5Op0zr27IFli2DCxdg4MCf7/O0hIQE2rRpI1MzquDixYtcunRJ83Ph8/LyAGhw2TQfLeratSvFxcUUFhaqHaVOK/+dsmnTplt/kj17lM57nTvD//6ndO/Jzoa//EVplyo8qnxN1MrvvmPJhx+y9ttvOfbUU3wLbOjWjb+0bInd7eaF7Gz67NhBSHIyY/bu5d3jx9lbWIi7DhVWqampmn9NEHWPFFJCM3Q6GDsWdu6Eb76BY8egb19lu5Pt2z2TwWq1snDhQu6///4rWtmKyvnss88A+O1vf6tykqopL6Tqa/wK/x133AHAt99+q3KSuq18bciOHTsq/8Vbtyq/ILt2VT7+97+VaX2PPw4WSzUnFTfjao0lYmJiSEpKYmJ8PLp9+5jVujVrunblUkwMyd27838tW1JYWsrMrCyit22jaUoK9+3fz7zTpzlhs6n8HdWcU6dOcfToUVkfJTRH3gkKzdHrYdIkpcPf559DZib06gVjxtR8QbVs2TIKCgpkWl8VLV26FJ1OR3x8vNpRqqSgoACAkLK1EVpV3jRl1apVKiep23x9fTGZTJXb1yslBUaNgj594KeflLVQBw9qZ0O+WiwhIeGqjSXKi6nyffQATHo9A4OCmNW6Nf/r1o1LMTGs7dqVB5o25WBxMdMOHKBlWhqdt2zhycxMVl24QHHZPle1wcKFCwHo16+fykmEqByduy6NG4taqbQUvvxS2dT30CGIi1P2pOrZs/rPNXLkSAoLC0lOTq7+J69DAgICCAkJ4ejRo2pHqZIHHniAefPmUVJSgslkUjtOlZjNZiIiItizZ4/aUeq0Jk2a4HK5yMnJuf4nrl8P//gH/PCD0u70b39TrjB58+Z74pZdcDj44dIlvr94ke8vXuSk3Y6vXs/Q4GDuatCAu0JCaKfhkcdevXqxa9cunE6n2lGEqBQZkRKa5+OjbIWyfz/Mn69clO3VSymotm2rvvOcOnWKNWvWyN5RVXT06FGKiooYOnSo2lGqrKioCEDzRRRA06ZNyc7OVjtGndeiRYuKJia/4nYrhdOQITB0KJw/D4sWKcPzkydLEVWLhRiN3N2oEZ917Mjx/v3Z17s3s9u0weFy8dThw4Rv3kzE5s38OTOT/126hMPlUjtypRw4cIAWLVqoHUOISpNCStQaBgPcf79SUC1YoMxu6d0bRo+G6miC9fnnn2M0Gpk0aVLVn6wO++ijjwB48MEHVU5SdUVFRehqSeviqKgoCgoKsNvtakepk6ZPn05ycjIRERE4HI4rGn8kb9zI9FGjICYG7rgDioqUTnw7d8KECcp8Z1Fn6HQ6Ovv781TLlqzt1o0LAweyLCqKIcHBfJ2Tw+27dxO6aROT09P54uxZLjkcake+rtzcXAoLC+nVq5faUYSoNPntK2odg+HnEaovvoAjR6BfPxgxAjZuvLXndLvdJCQkEB8fr/lW12pbsWIFBoOBwYMHqx2lyqxWa60ppIYMGQIoG/MKz5s6dSpxcXEVjUs2btwIbjfJr79O3LBhTF21SpnH/N13SjOJMWO8Z4dyoapAg4GxoaF8HBHBif792dGzJ0+1bElmcTH3ZWTQcNMmbt+1i/dPnOCYFzas+PrrrwE0v2ZW1E2yRkrUeqWlsHixspnv3r3KrJhZs+D222/+fciOHTvo2bMn3377LXfddVfNBq7lzGYzbdu2JSMjQ+0oVda3b1927txZK0ZxDh8+THh4OI899hjvv/++2nHqpOTkZO666y4KCgqYNWkSd+7eTdzBgyRFRRHz9tvKaJQUT6ISTthsJF24QOL58/wvNxeH203PgADGN2zI+NBQOvr7qx2RuLg4VqxYQVFRkab3FRR1k4xIiVrPx0fZVmXXLli6FAoLYfhwGDAAVqxQlh3cyPz582ncuDEjRoyo+cC12I4dO7Db7bXm77G4uBifWrIupV27dvj4+JCWlqZ2lDorZuBAVjzzDABLFy4k7vBhkt59l5g9e5RfWlJEiUpq4evLw82bs6prV84NHMiXkZG0s1h45ehRIrdupfOWLfz9yBH2qLhn1Y4dOwgKCpIiSmiSFFKiztDrIT5emRWzcqVSYMXFQffuyr5U1+ok63A4+OKLL7jvvvswGAyeDV3LfPLJJwDMmDFD5STVw2az1ar/E40aNeLw4cNqx6h7XC5YsgS6d2fwrFlEWCykAy+89RYxf/qTFFCiWgQZDExu3JivO3fm3MCBJEZF0SswkPdOnKDrtm1EbtnC80eOeHQjYJfLxZkzZ4iMjPTI+YSoblJIiTpHp4ORI5X1UuvWQaNG8JvfQGQkzJ0Ln32WcEX3slWrVnHu3DmmTJlCdnb2FXt/iMpZu3Ytvr6+teZFs6SkBKPRqHaMahMZGUlubi4ujXX80qzLCigmTIDQUJL/9S/Oms28++67vPDCC7LVgqgRFh8fxoSGkhAZSc7AgXzbpQv9g4L414kTRG/bRtTWrbyUnc1Bq7VGc6SkpOByuSo2BRdCa6SQEnWWTqesl1q9GrZsgagoePBBmDVrCMOHT2P//mxA2VSxa9euBAUFMW3atIpF+aJyXC4XR44coXPnzmpHqTZ2u71WFVIDBw7E7XaTmpqqdpTa7ZcFVKNGsHEjyS+8QNzzz5OUlMSf/vQnkpKSiIuLk2JK1CiTXs9dISH8p2NHzg4cSFJUFD0CA3nj+HEitmyh57ZtvHX8OCdqoFHFl19+CcBvf/vban9uITxBCikhUNqkL1kC+/bB8OGtycqaS9eu05gxYxdJSUnExsYybdo05s6dS+vWrdWOq0lr1qyhtLSUuLg4taNUG7vdjtlsVjtGtSnvmrVixQqVk9RS1yigWLOGZJRF90lJScTExAAQExMjxZTwKLNeT2xoKAsiI8kZMIBFnTvTxteX57KyCEtL47Zdu/js9Glyq6ml+saNGzEajURERFTL8wnhadK1T4irOHoUnn8+mwUL7sTtzqRp0/4sXPhfBg5srXY0zbrnnnv46quvOH36NE2aNFE7TrWoV68ejRs35tChQ2pHqRYulwuDwUBMTAwbNmxQO07t4XLBsmXw4ouwZ4/SMvSFF5R9ocpMnz6dqVOnVhRRl0tOTiYhIaFijaEQnpbndLL03Dk+z8nhh0uXMOl0jAkN5beNGzOyQQNMt7iXmZ+fH82aNSMzM7OaEwvhGVJICXEd//nPYqZNm0hAwDpstiHccw88/bQyDVBUTosWLcjPzyc/P1/tKNXGz8+Pdu3asXfvXrWjVJuGDRui1+s5e/as2lG0z+3+uYDavfuqBZQQWnOqpIQvc3JYcOYMu4uKCDUaubdRI6Y2aUL3gICb3lsvJyeHxo0bc++99/L555/XcGohaoZM7RPiGrKzs1mw4P+xbt06evR4kWefzWbdOujSBe66C3788eZapwtlCtypU6fo3r272lGqldPpxNfXV+0Y1ap9+/acP39e7Rja5nZDYiL06AHjx0NICGzYAGvXShElNK+Z2cxTLVuyq3dvdvfqxdTGjfk6J4ee27fTdds23jl+nHM3sbfeF198AcD48eNrOrIQNUYKKSGuIjs7u2JN1JAhQ0hImEty8jR++CGbBQvgxAm47Tbo0we+/hqcTrUTe7dvvvkGt9vNpEmT1I5SrVwuV63b+6Rfv364XC7S09PVjqI9bjckJUGvXspeC8HBsH49/PADDBqkdjohql10QABvhodzon9/VnTpQoSfH89kZdEsNZXx+/bx7YULOMu6gCYkJLB27Vrat2+P2Wzmz3/+MwBPP/00a9eulY64QpOkkBLiFy4vosobS7Ru3Zq5c+fyhz9MIyYmm927YdUq5X3S5MkQHg7vvQcFBepm91blVx5/97vfqRukmtXGQio2NhaAZcuWqZxEQ9xuZXO6Pn1gzBjw91eGrH/8EQYPVjudEDXOoNczOiSEhZ07c3rAAN5q146s4mJi9+6lVVoas7KyoEEDRowYQWZmJna7vWKvqqysLEaMGEHz5s1V/i6EqDxZIyXELyQkJDBkyJCrdufLzs5m/fr1TJ06teK+nTvhrbeUkSl/f/jDH+CPf4QWLTwY2suFhoai0+k4d+6c2lGqlU6nY8KECSxatEjtKNXG6XRiNBoZMWIE33//vdpxvJvbDWvWwN//DmlpMGAAvPSSMlwtm+iKOs7tdrOjsJDPTp/m87NnyZ88GU6duubnh4eH15rGPaLukEJKiGpy/Dj861/w8cdgtSqb/D75JPTsqXYydeXn5xMUFMTo0aNrVVvt8oJj6tSpzJs3T+041So4OJjAwECOHz+udhTv9eOP8PzzkJwMffsqDSVGjJACSoirKCotJdjPD+d11k6ZTCZKSko8mEqIqpOpfUJUk5Yt4Y03lILq9ddh0yZlqcSQIUrjrtJStROqY+7cuQDcd999KiepXhcvXgQgMDBQ5STVr23btpw5c0btGN4pORmGDVNGnYqLYcUKSE2FO++UIkqIa/D38UF3g+v2cl1faJEUUkJUs3r1lJGoQ4dg0SKlgBo3Djp0UNZR1aLu3zdl8eLF6HS6Wtdo4sKFC0DtLKT69OmD0+nk2LFjakfxHlu2KMXSoEFw8aJydWTrVhg9WgooIW7Cjdqi32zbdCG8iRRSQtQQgwEmTFAuYG/eDP36wf/9n7J26okn4PBhtRN6xq5du2jevDkGg0HtKNUqNzcXUDblrW3uvPNOAJYuXapyEi+wa5fSQKJvX2W4eeFCZWHk2LFSQAlRCWFhYVV6XAhvJIWUEB7Qpw98/jlkZyuNKP77X2jfXnl/tnZt7d2P6sSJExQWFjK4FnYuu3TpEgBBQUEqJ6l+o0aNAmDdunXqBlFTejpMnAjdu8NPPyk/tHv3Kvfp5aVTiMqaM2fONUeddDodc+bM8XAiIapOXg2E8KDmzeHll5UL2598AkePwvDhEBUFc+ZAYaHaCavXxx9/DMC0adNUTlL9ygup4OBglZNUP19fX/z9/dmzZ4/aUTzv0CG47z5l5+3t22HuXNi/X7nPx0ftdEJo1smTJ1m9ejXh4eGYTCaMRiMmk4nw8HBWr17NyZMn1Y4oRKVJ1z4hVOR2K/t1vv++suQiMBAeeAAeeUQZsdK67t27s3fvXux2O/padhX/ww8/5NFHH+W7776rGMGpTTp16kRWVhY2m03tKJ6RnQ3/+AckJECTJjBrFkybBiaT2smEEEJ4qdr1zkYIjdHpYOhQWLwYjhyBhx9WZhB16AAjRyoNwbTc7W///v20bdu21hVRAHl5eQA0aNBA5SQ1o0ePHibeM88AABDqSURBVJSUlFR0J6y1Tp2Cxx5TfuhWrIA334TMTJgxQ4ooIYQQ11X73t0IoVFhYfDPfyrT/ubNgwsXIC5OGZl6/XU4f17thJWzZMkS7HZ7rWzGAMr+WFB7C6nhw4cDkJiYqHKSGnL+vNL9pV07+OILZR+ow4eVTjC+vmqnE0IIoQFSSAnhZXx9YepUpbPy5s1Kt+Xnn1e6/U2ZomxZo4UJuS+88AIAM2fOVDdIDSkvpEJCQlROUjPGjh0LwNq1a1VOUs3y8pQfqDZt4KOP4OmnleHgZ5+FgAC10wkhhNAQWSMlhAacP6+sef/3v5X3fF27KtMA771XWVfljUwmEw6Ho9Zusjh16lTmz59PaWlprZy6CGCxWGjdujUZGRlqR6m6oiJlMeLrrysb6T72GDzzDISGqp1MCCGERtXOV38hapnQUOXCeWYmrFwJrVopDSmaNVMKqt271U54JafTicPhwLcWT5EqLGuxWFuLKIDmzZtz/PhxtWNUTUmJUkC1a6eMRN1zjzKF7403pIgSQghRJbX3HYAQtZBerzShSExURqaefFL5uFs3ZcPf//wHrFa1U8J7770HQK9evVROUnOKioquuSdKbREdHU1RURFWb/hPVVlOp/ID0aGDsu5p5Eg4cAD+3/9TrkAIIYQQVSSFlBAaFRYGL72k7EW1ZAkEB8ODDyrvER97DNTcAuiDDz4A4O2331YvRA2zWq21ejQK4LbbbgNg5cqVKiepBJcLFi5UNmebNg1694Z9+5QOLm3aqJ1OCCFELVK73wUIUQcYjTBuHKxapcxYeuwxpZ16167Qt6+y8W9BgWczlU8H6927t2dP7EFWqxWfWr5Ba3x8PADff/+9yklugtut/BD07g13360UTdu2waJFEBmpdjohhBC1kBRSQtQibdrA7Nlw7BgsXaosAfnDH5RRqunTIS2t5jv+2Ww2SktLCfTWLhjVxGaz1fpCqkWLFhiNRrZu3ap2lOvbtEnZkG3UKLBYlF2uV66Enj3VTiaEEKIWk0JKiFrIaIT4ePj2W8jOhr/8BVavhv79oUsXePttOHeuZs79zDPPADBs2LCaOYGXsNlsGAwGtWPUuCZNmnDkyBG1Y1zd7t0QGwsxMUpb8xUrYONGGDxY7WRCCCHqACmkhKjlwsKUZmVZWfD999C5M8ycCc2bw4QJSrHldFbf+b766isA3n///ep7Ui9UUlKC0WhUO0aN69SpE3l5eTir8z9JVWVmKr3/u3VTGkh8+SXs2AGjR0MtbwAihBDCe0ghJUQd4eMDI0bA11/DqVPw5pvK+9HYWKXYmjkTfvqp6uc5d+4cOp2OsLCwqj+ZF7Pb7ZhMJrVj1LieZdPj5s2bp24QgNOnlb7/kZHK9L2PPoL9+2HyZKWlpRBCCOFB8sojRB0UGgqPPw67dsH2/9/e3cbmXdd7HP+0vbqt3VY6WlZuBmFbFRwyTgYqM0PUkIBREnjgXfA4bzLwAZwHHs+cyRIimTAiCjGGBEeiTE0WTUhwKonEAXGKiojTATLLzdnYDbvturu2V9vrPLjOJhwOcz9d27G9XsmedOv/+l170Pbd3////T1V35n6znfqP5/On1//+bS3t/y627ZtS61WS0dHx/Ff9AlmcHAwEydOHO9ljLr58+cnSe69997xW8SePclXvlI/C2rVquRrX6v/FuDGG+v3sQLAOBBScApraEjmzaufV7p1a31qdEdH/Zf+Z56ZfPzjyc9/fuy3/t1yyy1Jkuuvv34UV31iqFarp0RIfeQjH0mSbNiwYexf/ODB5M47k1mzkm99K/niF+v3qC5eXB8qAQDjqKFWG+0ZXsBbzdatyQ9/WD9655lnkq6u5IYbkk9/uj5W/bBFixalu7s7S5cufd0zNA0NDbnjjjvS09OTFStWjP0bGAMtLS15+9vfnnXr1o33Ukbd4fOyRkZGxuYFq9X6Ybpf/WqyfXt99OTSpfW6B4AThB0p4A3OOiv50peSv/ylfhTPJz6RfP/79Wf7585Nvv71ZPPmZObMmVmyZMkbBhHUarUsWbIk3d3d4/QORt/w8HBaTpFdkba2ttRqtdEfOFGr/f0w3Ztuqo80f/755NvfFlEAnHDsSAHHpFqtT/1buTL5yU+SwcGkVmtO8uY/XFcqlVSr1bFb5BhqamrKlVdemTVr1oz3Ukbde9/73jzxxBP53ve+l4ULF47Oi/zyl/WJJ3/4Q/08qNtvr5c7AJyg7EgBx6S5uT7h70c/SrZtS+p37B19h+KEGpl9nNVqtbS2to73MsbE5z73uSTJ/ffff/wv/vTTydVXJ1ddVR8t+dhj9QfzRBQAJzghBRRrb08+//nxXsX4qtVqmTJlyngvY0x85jOfSZKsX7/++F30hRfqZ0HNm5ds3Jg8+GDyxBPJlVcev9cAgFEkpIDR05Rcfv/lue3x2/LUlqcyUhujYQWjrL+/P0lOmZCqVCppaGhIX1/fv36x7duTW25JLrywfhbUihX1h/Guv95hugC8pVTGewHAW1elUjnq7XuNExpz7mnn5htPfCO3PnZrpk+enqtnX51ruq/JVbOuyvTJ08dwtcfP7t27kyRTp04d55WMndbW1hw4cOCfv8C+fck3v1k/CbqpKVm2rB5Up8jtkQCcfOxIAf+0ZcuWHfXvb7/19vz4oz/Ozv/amUcXPprP/ttns+7VdbnhwRvSdVdX5t03L4sfWZxfvPCLHKweHKNV/+tOxZCaPXt2kuSRRx4p+8RqNbn33qS7O7njjuQLX6ifBfXlL4soAN7ShBTwT+vp6cny5ctTqbx+c7tSqWT58uXp6elJkjQ3Nef9578/y69annVfWJctX9ySldetzEXTL8oP/vyDXP2Dq9O+vD3v++77cuujt2bNS2tyqHpoPN7SMTkVQ+pjH/tYkuTuu+8+tk84PMp8zpzk5puTa65JNmyoz84//fRRXCkAjA3jz4FxVavV8tzO57LmpTV59OVH89jLj2X3od1pbmzOu855V64474osOG9B5s+Yn47WjvFebpLkoYceynXXXZf77rsvN95443gvZ0z09/enpaUlnZ2d2bFjx9H/8eOPJ4sXJ7//fX2U+fLl9QPIAOAk4hkpYFw1NDRkzhlzMueMObn53TdnpDaSZ7Y/k8f/+/H8auOvsnLdytz56zuTJBd2Xpj5M+bn8hmX5z3nvCcXTb8olcax/zLW29ubJGlvbx/z1x4vkyZNSpLs2bPnzf/RM8/Uz4L66U+Tyy5L1qxJPvCBMVohAIwtIQWcUBobGnNx18W5uOvi3Pzum1Or1fJS70v5zabf5Ncbf53fbv5tVq5bmeHacFoqLbn07Etz2VmX5dKzL828s+blgo4L0tTYNKprPBxS06ZNG9XXOVEsWrQo3d3dSZLh4eE0/O90vUqlkmXLlqVn3bqsaG1Nvvvd5Pzzk1Wrko9+NGl09zgAJy8hBZzQGhoaMmvarMyaNiufmvupJMmBwQN5autTeXLzk3lyy5NZvWF17vndPUmSlkpL5nbNzSVdl+SSMy/J3K65eef0d6Z90r+2e7Ro0aIsXLgwCxYsyN69e5P8PaTWrl2bBx54ICvqpxSfdGbOnJklS5a84eNDQ0NZsmRJHqhU6oeL3XNPctNNyYQJ47BKABhbnpECTgq9/b3549Y/5k/b/pSntz2dP7/65zy749kMjdTHs89om1G/hbBzTt5xxjtyQccFuaDzgnRN7jqyw3I0a9euzbXXXpvVq1fnoYceyl133ZUXX3wxmzdvPvLxBQsWjPbbHBfNzc1HHXPf2dCQHb29SVvbGK4KAMaXkAJOWoPDg/nrzr9m/fb1Wb99fZ7d8Wye2/lcenb3HDkcuG1iW7pP787bTn9bZk+bndmnz87M9pk5v/38zGibkeam5iPXOxxTV1xxRVavXp2HH344n/zkJ0/qiEpyTKHpWwkApxohBZxyBoYG8sKeF/L8zuezYdeG/G3339Kzuycv7nkxr/S9klrqXxYbGxpz9tSzc27buZnRNiMz2mZk4MWB3Pcf92V4aDhT26Zm1YOr8qEPfuiYYuNEM1Ibye5Du/Pq/lezbf+2bN2/NVv3bc2WfVvyyr5Xsmnvpmzq25Qt/7nlH17LtxIATjVCCuA1BoYGsnHvxrzU+1Je7n05m/Zuysa+jdnctzmv9L2Szfs2Z//d+5MdSa5JcnlSaayks7Uzna2d6WjpSEdrR6ZNmpb2Se05beJpaZvYlraJbZk6cWqmTJiSyc2T09LckpZKSyZVJmVC04RMrExMc2NzKo2VVBoraWxoPBJntVottdQyPDKc4dpwqsPVVEeqGRweTP9QfwaGBnJo6FAOVQ/lQPVA9g/uz76Bfekb6EvfQF96+3vTO9CbPYf2ZNehXdl1cFd2HtyZXYd2HdmZO+y0iafl7KlnZ0bbjJzTdk7Oazsvt33wtn/4/+ZbCQCnGiEFUGDt2rX58Ic/nOs+cV0eXPVgFn97cTov7MyOgzuy48CO7O7fnV0Hd2VP/57s7d+b3v7e9A305dDQ2B8wPGXClEydMDXtk9qP/Olo7UhHS0c6WztzRusZmT55erqmdKVrclfOnHJmJk+Y/Ibr/KNnpCqVSqrV6mi+FQA44ZjaB3CMDj8j9bOf/SwLFizIon9f9PdBE1ce/Rmp6nA1+wf350D1QA4MHjiyg9Q/1J/B4cEMDA9kaGQo1eFqhmvDGamNvG63qLGhMY0NjWlqaEpzU3OaG5uP7GRNqkxKa3NrWiotmTxh8pFdr+M1Bn7ZsmX/79S+1/49AJxq7EgBHIPXTu177WCJN/v4yeTwOVJLly593c7UkXOkenpO2tHvAPBmhBTAMXjtOVL/18l+jhQA8EZCCgAAoFDjeC8AAADgrUZIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIX+B1EwuEtDkD9xAAAAAElFTkSuQmCC" + ], + "text/plain": [ + "PyPlot.Figure(PyObject )" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "(-0.5,4.0)" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "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.5.0-dev", + "language": "julia", + "name": "julia-0.5" + }, + "language_info": { + "name": "julia", + "version": "0.5.0" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/src/elements.jl b/src/elements.jl index 0da9e8f..28cb335 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -271,6 +271,9 @@ Examples >>> interpolate(el, :temperature, [0.0, 0.0]) 15.0 """ +function interpolate(el::Element, field, xi::Number) + interpolate(el, field, [xi]) +end function interpolate(el::Element, field, xi::Vector) field = get_field(el, field) sum(get_basis(el, xi) .* field) @@ -280,6 +283,12 @@ function interpolate(el::Element, field, xis::Array{Vector, 1}) interpolate_(xi) = sum(get_basis(el, xi) .* field) map(interpolate_, xis) end + +""" +""" +function dinterpolate(el::Element, field, xi::Number) + dinterpolate(el, field, [xi]) +end function dinterpolate(el::Element, field, xi::Vector) fld = get_field(el, field) dbasis = get_dbasisdxi(el, xi) @@ -462,3 +471,68 @@ function fit_derivative_field!(el::Element, field, f, fixed_coeffs=Int[]) return end + +""" Find projection from slave nodes to master element. """ +function calc_projection_slave_nodes_to_master_element(sel, mel) + X1 = get_field(sel, :Geometry) + N1 = get_field(sel, :Normals) + X2(xi) = interpolate(mel, :Geometry, xi) + dX2(xi) = dinterpolate(mel, :Geometry, xi) + R(xi, k) = det([X2(xi) - X1[k] N1[k]]') + dR(xi, k) = det([dX2(xi) N1[k]]') + xi2 = Vector[[0.0], [0.0]] + for k=1:2 + xi = xi2[k] + for i=1:3 + dxi = -R(xi, k)/dR(xi, k) + xi += dxi + if abs(dxi) < 1.0e-9 + break + end + end + xi2[k] = xi + end + clamp!(xi2, -1, 1) + return xi2 +end + +""" Find projection from master nodes to slave element. """ +function calc_projection_master_nodes_to_slave_element(sel, mel) + X1(xi) = interpolate(sel, :Geometry, xi) + dX1(xi) = dinterpolate(sel, :Geometry, xi) + N1(xi) = interpolate(sel, :Normals, xi) + dN1(xi) = dinterpolate(sel, :Normals, xi) + X2 = get_field(mel, :Geometry) + R(xi, k) = det([X1(xi) - X2[k] N1(xi)]') + dR(xi, k) = det([dX1(xi) N1(xi)]') + det([X1(xi) - X2[k] dN1(xi)]') + xi1 = Vector[[0.0], [0.0]] + for k=1:2 + xi = xi1[k] + for i=1:3 + dxi = -R(xi, k)/dR(xi, k) + xi += dxi + if abs(dxi) < 1.0e-9 + break + end + end + xi1[k] = xi + end + clamp!(xi1, -1, 1) + return xi1 +end + +function has_projection(sel, mel) + xi1 = calc_projection_master_nodes_to_slave_element(sel, mel) + l = abs(xi1[2]-xi1[1])[1] + return l > 1.0e-9 +end + +""" +Calculate projection between 1d boundary elements +""" +function calc_projection(sel, mel) + xi1 = calc_projection_master_nodes_to_slave_element(sel, mel) + xi2 = calc_projection_slave_nodes_to_master_element(sel, mel) + return xi1, xi2 +end +