diff --git a/notebooks/2015-12-10-3d-mortar-assembly.ipynb b/notebooks/2015-12-10-3d-mortar-assembly.ipynb new file mode 100644 index 0000000..c41949c --- /dev/null +++ b/notebooks/2015-12-10-3d-mortar-assembly.ipynb @@ -0,0 +1,234 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 3d mortar assembly\n", + "\n", + "Author: Jukka Aho\n", + "\n", + "**Abstract**: Mortar assembly for 3d problems." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "using JuliaFEM" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "1-element Array{JuliaFEM.Core.Element{E<:JuliaFEM.Core.AbstractElement},1}:\n", + " JuliaFEM.Core.Element{JuliaFEM.Core.Tri3}([4,5,6],Dict{ASCIIString,JuliaFEM.Core.Field{A<:Union{JuliaFEM.Core.Continuous,JuliaFEM.Core.Discrete},B<:Union{JuliaFEM.Core.Constant,JuliaFEM.Core.Variable},C<:Union{JuliaFEM.Core.TimeInvariant,JuliaFEM.Core.TimeVariant}}}(\"geometry\"=>JuliaFEM.Core.Field{JuliaFEM.Core.Discrete,JuliaFEM.Core.Variable,JuliaFEM.Core.TimeInvariant}([[-1.0,1.0,0.1],[2.0,-0.5,0.1],[2.0,2.0,0.1]])))" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "nodes = Vector{Float64}[\n", + " [0.0, 0.0, 0.0],\n", + " [3.0, 0.0, 0.0],\n", + " [0.0, 3.0, 0.0],\n", + " [-1.0, 1.0, 0.1],\n", + " [ 2.0, -0.5, 0.1],\n", + " [ 2.0, 2.0, 0.1]]\n", + "\n", + "mel = JuliaFEM.Core.Tri3([4, 5, 6])\n", + "mel[\"geometry\"] = Vector{Float64}[nodes[4], nodes[5], nodes[6]]\n", + "\n", + "sel = JuliaFEM.Core.Tri3([1, 2, 3])\n", + "sel[\"geometry\"] = Vector{Float64}[nodes[1], nodes[2], nodes[3]]\n", + "R = [\n", + " 0.0 1.0 0.0\n", + " 0.0 0.0 1.0\n", + " 1.0 0.0 0.0]\n", + "sel[\"nodal ntsys\"] = Matrix{Float64}[R, R, R]\n", + "sel[\"master elements\"] = JuliaFEM.Core.Element[mel]" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(\n", + "2x3 Array{Float64,2}:\n", + " -1.0 2.0 -1.0\n", + " -1.0 -1.0 2.0,\n", + "\n", + "2x3 Array{Float64,2}:\n", + " -2.0 1.0 1.0\n", + " 0.0 -1.5 1.0,\n", + "\n", + "2x6 Array{Float64,2}:\n", + " -1.0 0.0 1.0 1.0 0.25 -1.0 \n", + " -0.5 -1.0 -1.0 0.0 0.75 0.333333,\n", + "\n", + "3x3 Array{Int64,2}:\n", + " 1 0 1\n", + " 1 1 0\n", + " 0 1 1,\n", + "\n", + "[0.039743589743589755,-0.1952991452991453])" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "time = 0.0\n", + "x0, Q = JuliaFEM.Core.create_auxiliary_plane(sel, time)\n", + "S = Vector{Float64}[]\n", + "for p in sel(\"geometry\", time)\n", + " push!(S, JuliaFEM.Core.project_point_to_auxiliary_plane(p, x0, Q))\n", + "end\n", + "\n", + "M = Vector{Float64}[]\n", + "for p in mel(\"geometry\", time)\n", + " push!(M, JuliaFEM.Core.project_point_to_auxiliary_plane(p, x0, Q))\n", + "end\n", + "\n", + "#S = Matrix[x' for x in S]\n", + "#S, M\n", + "S = reshape([S...;], 2, 3)\n", + "M = reshape([M...;], 2, 3)\n", + "P, neighbours = JuliaFEM.Core.clip_polygon(S, M)\n", + "C = JuliaFEM.Core.calculate_polygon_centerpoint(P)\n", + "S, M, P, neighbours, C" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAG7CAYAAABQNZVMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlcVXX+x/HXveyC4AKIiLgBsrpg7jVpjWW5Vu65hLlki+XP0ZycsmwdSy3HbHFBzTLTNNcZHVObXLISVFYRdwUEFXBhv/f8/jiKIoKg3Hvuhc/z8fDBnO3eN8TcD99zvotOURQFIYQQoobSax1ACCGE0JIUQiGEEDWaFEIhhBA1mhRCIYQQNZoUQiGEEDWaFEIhhBA1mhRCIYQQNZqt1gGq2oULF9i6dStNmzbFyclJ6zhCCCE0kpuby8mTJ3n88cdxd3cv87xqVwi3bt3K8OHDtY4hhBDCQqxYsYJnn322zOPVrhA2bdoUUL/xoKCgMs+bNGkSc+fONVMqIT9v85OfuXnJz9u8KvLzTkhIYPjw4cV1oSzVrhDeuB0aFBREeHh4mee5ubmVe1xULfl5m5/8zM1Lft7mVZmf990ek0lnGSGEEDWaSQvhH3/8wcsvv0xISAguLi40adKEwYMHc/To0Qpdn5WVxbhx4/Dw8MDFxYVHHnmE6OhoU0YWQghRw5j01ug///lP9u3bx8CBA2nVqhWpqanMnz+f8PBwfvvtN0JCQsq81mg00qtXLw4fPszUqVOpX78+CxYsoFu3bhw4cAA/Pz9TRhdCCFFDmLQQTp48mfbt22Nre/NtBg8eTFhYGB999BHffPNNmdeuWbOGffv2sWbNGp5++mkABg0aREBAADNmzODbb7+9r2xDhw69r+tF5cjP2/zkZ25e8vM2r6r8eeu0WI+wXbt26PV6/vjjjzLPGTRoELt37yYlJaXE/hdeeIEVK1aQmZmJnZ1dqeuioqJo164dBw4ckAfXGsvIyOC996YSG/s7Ol0RimJLaGgH/vGPWXh4eGgdTwhRzVW0Hpi916iiKJw/f56wsLByz4uOjr5j8Pbt2/P111+TlJRU7q1Voa309HT69+/C8OHH6N8fdDowGiExMZ5+/X5l/fp9UgyFEBbB7L1Gv/32W1JSUhg8eHC556WmptKwYcNS+2/su72lKCzL+++/zvDhxwgOVosggF4PwcHw7LPHeO+9qdoGFEKI68xaCBMTE3nppZfo0qULo0aNKvfcvLw8HBwcSu13dHQE1KlzhOWKjf2dsuYzCApSjwshhCUw263RtLQ0evXqRd26dVmzZg26G82EMjg5OZGfn19qf15eXvHx8kyaNAk3N7cS+4YOHSoPtM1EpyuirP/Eer16XAghqsrKlStZuXJliX3Z2dkVutYshTA7O5snnniCy5cv8+uvv+Ll5XXXaxo2bHjH25+pqakAeHt7l3v93LlzpbOMhhTFFkXhjsXQaFSPCyFEVblTQ+dGZ5m7Mfmt0by8PPr06UNycjKbNm0iMDCwQte1adOGqKgobu/Uun//fpydnQkICDBFXFFFQkM7kJBw52MJCRAQ4GPeQEIIUQaTFkKDwcDgwYPZv38/q1evpmPHjnc8Ly0tjcTERIqKbt4uGzBgAOfPn2ft2rXF+y5cuMDq1avp06fPHYdOCMvxj3/MYsWKFsTFqS1AUL/GxcHixfDkk/u4dq2MSimEEGZk8gH1GzdupE+fPly4cIEVK1aUOH5juaRp06axfPlyTp48ia+vL6AWwk6dOhEREUF8fHzxzDKKovDOO++YMraoAh4eHqxfv4/33pvKihUr0OuLMBpt8ff34q23zlK79hViYnoRHr4fe3sZRiGE0I5JC+GhQ4fQ6XRs3LiRjRs3ljim0+mKC6FOpyvVeUav17NlyxamTJnCvHnzyM3NpUOHDixfvhx/f39TxhZVxMPDg88+i2Tv3v9SUHAOBwcv2rdP4ODBh7h69SB5eSeIje1P69Y/Y2PjqHVcIUQNZdJbozt37sRgMGA0Gkv9MxgMxedFRkZiMBiKW4M31KlTh4ULF5KRkcHVq1fZsWOHdICxcra2LoSGbsTeXu3sdPnyXo4cGV3qWbAQQpiLLMMkzM7R0YewsI3o9bUASE9fycmTb2sbSghRY0khFJqoXTuc4ODvAPWW+KlTM0lLK3sSdiGEMBUphEIz7u79aNFidvH2kSPPk5X1Pw0TCSFqIimEQlM+Pq/h7f0CAIpSSGzsU+TkVGzhZiGEqApSCIWmdDodfn7/om7dxwAoKrpETEwvCgsvaZxMCFFTSCEUmtPrbQkJ+YFatdRltXJzjxIb+zRGY4HGyYQQNYEUQmERbG3daNVqM3Z2ngBkZ//CkSPjZFiFEMLkpBAKi+Ho2ISwsA3o9erg+vPnl3H69IcapxJCVHdSCIVFcXXtSGDgzWEUJ05MJz19lYaJhBDVnRRCYXE8PQfQrNnNlmBCwiiys/dpmEgIUZ1JIRQWydf3dby8RgOgKPnExvYjN/eExqmEENWRFEJhkXQ6HQEBX1CnTncACgszrg+ryNI4mRCiupFCKCyWXm9PSMiPODm1BCAnJ4H4+IEYjYUaJxNCVCdSCIVFs7OrS6tWm7G1rQ9AZuZ2jh59SYZVCCGqjBRCYfGcnFoQGvoTOp09AKmpCzlzZvZdrhJCiIqRQiisQp06DxIYGFm8ffz4VDIy1mmYSAhRXUghFFajQYNhNG369vUthYSEZ7l8+U8tIwkhqgEphMKqNGnyFg0aDAfAaMwlNrYPeXlnNE4lhLBmUgiFVdHpdLRsuQg3twcBKChIIyamN0VFVzROJoSwVlIIhdXR6x0ICVmHo2MLAK5dO0x8/GCMxiKNkwkhrJEUQmGV7O3drw+rqAvApUv/5tixSRqnEkJYIymEwmrVqtWSkJC16HS2AJw7N5+zZ/+lcSohhLWRQiisWt263QgIWFi8nZz8GhcvbtYwkRDC2kghFFavYcPn8PV94/qWkbi4wVy9ekjTTEII6yGFUFQLzZq9i4fHQACMxmvExPQmPz9F41RCCGsghVBUCzqdnsDAZdSu3RGA/PyzxMT0wWC4pnEyIYSlk0Ioqg0bGyfCwtbj6NgUgKtXo4iPfxZFMWgbTAhh0aQQimrF3r4BYWGbsLFxBeDixfUcO/a6xqmEEJZMCqGodpydQwgJWQPYAHD27GxSUr7SNpQQwmJJIRTVUr16PQgIWFC8nZT0EpcubdMwkRDCUkkhFNWWt/c4fHwmX98yEBc3kGvX4jTNJISwPFIIRbXWosU/cXfvD4DBcJnDh3tRUHBe41RCCEsihVBUazqdDUFBK3BxCQcgP/8UMTH9MBhyNU4mhLAUUghFtWdj40xY2EYcHHwAuHJlP4mJo1AUo8bJhBCWQAqhqBEcHLyvD6twASAjYzUnTrypcSohhCWQQihqDBeX1gQHf8+NX/vTpz8gNTVS21BCCM1JIRQ1Sv36vfDz+7R4OylpHJmZOzVMJITQmhRCUeP4+LxCo0YvA6AoRcTFPUNOzhGNUwkhtCKFUNRILVrMpV69JwEoKsq8PqzigsaphBBakEIoaiS93pbg4O9xdm4FQF7eMeLinsJozNc4mRDC3ExeCK9du8aMGTPo2bMn9erVQ6/Xs2zZsgpdu3TpUvR6/R3/paenmzi5qO5sbWsTFrYJe3svALKzd5OY+DyKomicTAhhTramfoOMjAzeffddmjRpQps2bdi1axc6na5Sr/Huu+/SrFmzEvvc3NyqMqaooRwdGxMaupGDB/+C0ZhLevq31KrlT9OmM7SOJoQwE5MXQm9vb9LS0vD09OTAgQO0b9++0q/xxBNPEB4eboJ0QoCr6wMEBX1HXNzTgMLJk2/j5ORHgwbPah1NCGEGJr81am9vj6enJ8A933JSFIUrV65gMMgCq8I0PDz607z5rOLtxMTRZGXt1jCREMJcrKKzTPfu3XFzc8PZ2Zl+/fqRnJysdSRRDTVuPJmGDccBoCgFxMb2Jzf3mMaphBCmZvJbo/fD2dmZiIgIunfvjqurK3/++Sdz5syhS5cuREVF4ePjo3VEUY3odDr8/eeTl3eCzMz/UlR0kcOHexEevg87u7paxxNCmIhFtwgHDhzI4sWLGT58OH379mXmzJls3bqVixcv8v7772sdT1RDer0dwcE/UKtWMAC5uUeIi3sGo7FA42RCCFOx6BbhnXTt2pWOHTuyffv2cs+bNGlSqZ6lQ4cOZejQoaaMJ6oBO7s6hIVtIiqqI4WFGWRl7SQpaQItWy6qdI9nIYR5rFy5kpUrV5bYl52dXaFrra4QAvj4+JCUlFTuOXPnzpWepuKeOTk1IzR0AwcPdkNR8klLW4KTkz9NmkzTOpoQ4g7u1NCJioqiXbt2d73Wom+NluX48eN4eHhoHUNUc25unQgKWl68feLE30lPX6NhIiGEKVhMIUxLSyMxMZGioqLifRkZGaXO27JlC1FRUfTs2dOc8UQN5ek5iGbNbj6PTkwcweXL+zVMJISoama5NTp//nyysrJISUkBYMOGDZw+fRqAiRMn4urqyrRp01i+fDknT57E19cXgC5duhAeHk67du1wc3MjKiqKJUuW4OvryxtvvGGO6ELg6/t3cnKSOH9+GUZjHjExfQkP34+TU1OtowkhqoBZCuHs2bM5deoUoHZRX7duHWvXrkWn0zFy5EhcXV3R6XSlOiIMGTKEzZs3s23bNnJycvD29mb8+PHMmDFDbo0Ks9HpdLRs+TV5eSfJzv6FwsJ0YmJ6Ex6+B1tbmepPCGunU6rZDMM3Ho4eOHBAOstYiL17fSgoOIeDgw+dO5/ROs49Kyy8RFRUJ3JzjwJQt+5jhIVtRq+3yj5nQlR7Fa0HFvOMUAhLZ2dXj7Cwzdja1gMgM3MbycmvyGoVQlg5KYRCVEKtWv6Ehv6ETmcHQErKl5w9+6nGqYQQ90MKoRCVVKfOQ7RsuaR4+9ixyVy4sF7DREKI+yGFUIh74OU1nCZN3rq+pRAfP4wrV6I0zSSEuDdSCIW4R02bvo2npzqThdGYQ0xMH/LyzmqcSghRWVIIhbhH6rCKJbi6dgGgoCCFmJjeFBVd1TiZEKIypBAKcR9sbBwJDf0JR8fmAFy7doiEhKEoiiwiLYS1kEIoxH2yt/cgLGwzNjbq4PqLFzeRnDxZ41RCiIqSQihEFXB2DiQ0dC06nTq4/ty5zzh37nONUwkhKkIKoRBVpG7dRwgI+LJ4++jRiVy8uEXDREKIipBCKEQVatjweRo3fv36lpH4+MFcvXpY00xCiPJJIRSiijVv/gHu7s8AYDBcJSamN/n5qRqnEkKURQqhEFVMp9MTFLSc2rU7AJCff4bY2L4YDNc0TiaEuBMphEKYgI1NLUJD1+PgoK6teeXKnyQkjEBRjBonE0LcTgqhECbi4OB1fVhFbQAuXFjH8eN/1ziVEOJ2UgiFMCEXl1BCQlYDNgCcOTOLlJSF2oYSQpQghVAIE6tX73H8/f9VvH306ItcurRdw0RCiFtJIRTCDBo1moCPz2sAKEoRcXEDuHYtXuNUQgiQQiiE2bRo8Qn16/cBwGDIJiamNwUF6RqnEkJIIRTCTHQ6G4KCvsPFpS0AeXkniI3tj8GQp3EyIWo2KYRCmJGtrQthYRuxt/cG4PLlfSQmPifDKoTQkBRCIczMwaERYWGb0OudAcjIWMXJkzM0TiVEzSWFUAgN1K7dluDglYAOgFOn3iMtbbm2oYSooaQQCqERd/c+tGgxp3j7yJExZGX9omEiIWomKYRCaMjH51W8vV8EQFEKiY19ipycJI1TCVGzSCEUQkM6nQ4/v8+oV68nAEVFmcTE9KKw8KLGyYSoOaQQCqExvd6W4OBVODuHApCbm0xs7FMYjfkaJxOiZpBCKIQFsLV1JSxsE3Z2DQDIzv6VI0fGoiiKxsmEqP6kEAphIRwdmxAWthG93gmA8+e/4dSp9zVOJUT1J4VQCAvi6tqeoKBvirdPnnyT8+dXaphIiOpPCqEQFsbD4xmaN/9n8XZiYgTZ2Xs1TCRE9SaFUAgL1LjxFLy8ngdAUfKJje1Hbu5xjVMJUT1JIRTCAul0OgICvqBOnUcAKCy8cH1YRZbGyYSofqQQCmGh9Ho7QkLWUKtWIAA5OYnExQ3AaCzUOJkQ1YsUQiEsmJ1dXcLCNmNn5w5AVtbPJCVNkGEVQlQhKYRCWDgnp+aEhq5Hp3MAIC1tMWfOfKJxKiGqDymEQlgBN7cuBAZGFm8fP/46GRlrNUwkRPUhhVAIK9GgwVCaNp15fUshIWE4ly//oWkmIaoDKYRCWJEmTf5BgwYjADAac4mN7Ute3mmNUwlh3aQQCmFFdDodLVsuxM3tIQAKCtKIielNUdFljZMJYb1MXgivXbvGjBkz6NmzJ/Xq1UOv17Ns2bIKX5+VlcW4cePw8PDAxcWFRx55hOjoaBMmFsKy6fUOhIauw8nJD4Br12KIjx+M0VikcTIhrJPJC2FGRgbvvvsuR44coU2bNoD6V21FGI1GevXqxcqVK5k4cSKzZs0iPT2dbt26kZycbMrYQlg0O7v6hIVtxta2LgCXLv2H5ORXZViFEPfA5IXQ29ubtLQ0Tpw4wccff1ypa9esWcO+fftYtmwZb775Ji+++CK7du3CxsaGGTNmmCixENahVq0AQkPXodPZAZCSsoBz5+ZpnEoI62PyQmhvb4+npydApf9aXbNmDV5eXjz99NPF+9zd3Rk0aBDr16+nsFBm2BA1W506D9Oy5cLi7eTkSVy4sFHDRDVLRkYGUyIi6BUSQt+WLekVEsKUiAgyMjK0jiYqwaI7y0RHRxMeHl5qf/v27cnJySEpKUmDVEJYFi+vUfj6Tr++pRAfP5QrV+Q5uqmlp6czuHNnnlm6lE3x8WxISmJjfDzPLF3K4M6dpRhaEYsuhKmpqTRs2LDU/hv7UlJSzB1JCIvUrNlMPDwGA2A0XiMmpg/5+ec0TlW9ffz663xw7BidgBu9HvRAJ+D9Y8eYNXWqduFEpVh0IczLy8PBwaHUfkdHRwByc3PNHUkIi6TT6QkMjMTVtRMABQXniInpQ1HRVY2TVV/xv/9OxzKOdbx+XFgHW60DlMfJyYn8/PxS+/Py8oqPl2XSpEm4ubmV2Dd06FCGDh1atSGFsBA2Nk6Ehq4nKqojeXknuXo1moSEZwkNXYtOZ6N1vOojJwd+/BGbEycoq/+7HrApkuEs5rRy5UpWrlxZYl92dnaFrrXoQtiwYcM73v5MTU0F1B6pZZk7d+4dny8KUZ3Z23sSFraZqKguGAzZXLy4gWPHpuDnN0fraNYvOhoWLYJvv4XsbAyAAncshkbAYGvRH6/Vzp0aOlFRUbRr1+6u11r0rdE2bdoQFRVVqrfp/v37cXZ2JiAgQKNkQlguZ+dgQkLWAGor8OzZuZw794W2oaxVVhZ88QW0awfh4bBgAVxvZQQD+8u4bD8Q3KGDuVKK+2QxhTAtLY3ExESKbrmdMGDAAM6fP8/atTdn2b9w4QKrV6+mT58+2NnZaRFVCItXr95fCQi4WfyOHn2FS5e2apjIiigK/PorjBoF3t7w4osQFXXzeK1a8NxzTN24kTdatGAfaguQ61/3AdNbtGDqrFnmzy7uiVna7vPnzycrK6v4NueGDRs4fVqdKHjixIm4uroybdo0li9fzsmTJ/H19QXUQtipUyciIiKIj4+nfv36LFiwAEVReOedd8wRXQir5e09ltzco5w58zFgIC5uIG3b7sXFJVTraJbp/HlYvly9/XmnoVnt28Pzz8OQIeDmhgewqmNHZj39NO/t3o0NYPDyIrhnT1bNmoWHh4e5vwNxj8xSCGfPns2pU6cAdXq1devWsXbtWnQ6HSNHjsTV1RWdTldq6jW9Xs+WLVuYMmUK8+bNIzc3lw4dOrB8+XL8/f3NEV0Iq9a8+Ufk5iZz4cI6DIYrxMT0Ijx8Pw4OXlpHswwGA2zbpha/DRvg9g4uderAiBFqAWzdutTlHh4efDx2LOzere546y2YMMEMwUVVMkshPHHixF3PiYyMJDIystT+OnXqsHDhQhYuXHiHq4QQ5dHp9AQFreDgwYe5cuVP8vNPExvbjzZtdmJjU0vreNo5eRIiI2HJEjh7tvTx7t1hzBh46ikop3e6qB6kW5MQ1ZyNTS1CQzcQFdWR/PwzXLnyOwkJIwkJ+QGdzmK6CZhefr7a6lu0CP77X/VZ4K28vCAiAkaPBj8/bTIKTUghFKIGcHBoSFjYJqKju2IwXOXChR85cWI6zZt/qHU004uPh8WL1ed/Fy6UPKbXw5NPwtix8MQTIB3waiQphELUEC4urQgOXkVMTB/AyOnTH+Hk5E/DhqO1jlb1rl6F1ath4ULYt6/08WbN1Fufo0ZBo0bmzycsihRCIWqQ+vWfxM/vM5KTXwEgKWk8jo5NqVv3EY2TVQFFgT/+UG99fv89XLlS8ri9PTzzjFoAu3VTW4NCIIVQiBrHx+dlcnOPcu7cPBSliLi4Z2jbdh/OzoFaR7s3ly7BihVqAYyJKX08NFS99fnss1C/vvnzCYsnhVCIGsjPbw65uce4dGkzRUVZ14dV/Ia9vZWMfTMa4Zdf1OL3449qR5hbubjA0KFq6699e9CVNSuoEFIIhaiRdDobgoNXEh39ENeuHSIv7zixsU/RuvV2bGwctY5XtpQUWLpUHfZw7Fjp4507q8Vv0CC1GApRAVIIhaihbG1rExa2iaioDhQUpHL58h6OHHmeoKAVpSa30FRREWzZorb+tmxRB8Hfqn59GDlSHfQeEqJNRmHVpBAKUYM5OvoQFraR6Oi/YDTmkJ7+HU5O/jRr9rbW0dQW3+LFagvw+oozJfToobb++vWDO6xbKkRFSSEUooarXbsdwcHfERv7FKBw6tQ7ODn54eU13Pxh8vJg7Vq19bdzZ+njjRqpA94jItQhEEJUASmEQgjc3fvRosUnHDs2GYAjR57H0bEJdeo8ZJ4Ahw+rxW/FCsjMLHnM1hb69FFbf48/DjayyLCoWlIIhRAA+PhMIicnidTUr1CUAmJjnyI8/Ddq1TLRdGOXL6vj/RYtUsf/3c7fXy1+I0eq058JYSJSCIUQgLoyjL//v8jLO0Fm5jaKii5eH1axDzu7elXzJoqizvSyaBGsWgU5OSWPOzrCwIFqAXzoIRn2IMxCCqEQopheb0dIyA9ERXUlJyeO3Nwk4uKeoVWrrej19vf+whkZ8M03agFMSCh9vE0bddD7sGHq0kdCmJEUQiFECba2bteHVXSksDCdrKxdHDkyjsDAyMoNqzAaYft2tfj99BMUFpY87uqqzvYyZgyEh1ftNyFEJUghFEKU4uTUlLCwDRw82A2jMY/z55dRq1YATZq8cfeLz5y5udbf9QW5S3joIbX4DRgAtWrwmojCYsiss0KIO3J17Uhg4PLi7V27ptOsmRcODg7Y29vj4OCAv78/CQkJamtv7Vp1SaMmTWDGjJJF0MMDpkyBxET43//UDjBSBIWFkBahEKJMnp4Dyc39gB073mDsWDAYzpc4npycTKvQUA7XqUPQpUslL9bpoGdPtfXXu7e6+oMQFkgKoRCiXL6+03j77Q8xGK7c8XiR0UjfS5c4evMCdbqziAho3NhsOYW4V1IIhRDl0ul0pKXll3vOabg57OHRR2XQu7AqUgiFEGXLyoLvvsNYUFDuaYqdHfzwg5lCCVG1pBAKIUpSFPj1V3XYw+rV7MjLo+gul1jUahVCVJIUQiGE6vx5WLZMXfEhKal4d0fAHbhQzqW+vr6mTieEyUghFKImMxhg61a19bdxo7r2363q1MF5xAj+1707rQYNouj246iPA5cunWimwEJUPRlHKERNdPIkvPUWCY0a8X6vXijr1pUsgt27w7ffqivCz5tH0FNPcfjwYXx9fbG1tcXOzg47Oxu8vWHBAoiLm8KVKwc0+3aEuB/SIhSipsjPh/XrufbVV6zesYNFwB6gPjAMaOblpQ55GD0a/EqvOBEUFMSpWwbJK4pCQsJw3nrrO/7zn3ycnXvyzDNRODrKkAlhXaQQClHdxcWhLFrEgchIFmVn8x1wBegBrNLp6PfEEziMH6/OCmNb8Y8EnU5Hy5aLmTjxGFFR+5ky5QLu7j159NHfsLWtbarvRogqJ7dGhaiOrl5V5/rs0gVCQ3ni009pn53NJuA14Hjjxmz74AMGnT2Lw+bN0LdvpYrgDTY2jnTqtInZsxuTlweTJ8cTHT0Qo/Fu/UyFsBxSCIWoLhQFfv8dxo2Dhg3V2V327QNgOLDJ1paTQ4Yw8+efaXbyJPz97+Dtfd9va2/vzmOPbeOjj1xITobXX9/K0aOT7vt1hTAXuTUqhLW7dAlWrFB7fsbElD4eFsbwsWPVJY/qVdECu7dxdg5k4MANpKX9lRkzjLz//nw++CAAH59XTPJ+QlQlKYRCWCOjEXbtoujrr9n644+4FxXR8dbjLi7qIrdjxsADD5hlpfe6dbszduwizp0bzddfQ5Mmr/Laa82pX7+Xyd9biPshhVAIa5KSAkuXcuLLL1ly5gyRwDlgIurAd7p0UYvfwIFqMTSzhg0jmDo1CUX5iDZtFOLjh9C27W5cXFqbPYsQFSXPCIWwdEVFsGED+b16scrHhx7Tp9P8zBnmAX2AP93c+HTSJIiLgz171CEQVVwEjx8/zvvvv09WVtZdz23e/H1efXUg9euDwXCVmJje5OenVGkeIaqSFEIhLFVyMrzxBvj6sr1fP7y3bGGIopAPLANSH32UL374gXbnz6ObMweCgyv18qNHj6Zly5aMHj36ruf++uuvvPnmm9hUYFUJnU5PYOAyatdWb9bm558lJqYvBsO1SuUTwlzk1qgQliQvT13pfdEi2LmzeHcwMBp4vkEDAsePV1t9TZve89uMHj2aTZs2kZGRQWZmJqNHj2bJkiVlnh8dHY2fnx+1a1dsfKCNjRNhYes5cKAj+fmnuHr1AAkJwwkJWYNOJ0s0CcsihVAIS3D4sFr8VqyAzMySx2xt8e7bl4/HjIHHHquStf727NlDRkYGABkZGezZs6fc86Ojo2nTpk2l3sPevgGtWm0mKqoLBsNlLlz4iePHp9GixceKWLnWAAAgAElEQVT3nFsIU5Bbo0Jo5fJl+PprLoaHs6p1a/jXv0oWwYAAmDULzp6FH3+EJ56osgVvu3btioeHBwAeHh507dq1zHMVReHgwYO0bdu20u/j7BxCSMhqQM195swnpKR8fU+ZhTAVKYRCmJOiwN69GCMi2O7pyZDx4/GOjmYkcB7AyQlGjoT//Q8SE2HKFGjQoMpjLFmyhN69exMQEEDv3r3LvS164sQJLl++XG4hLO95Y716jxEQ8DkAO3bAl1++wKVL/73/b0KIKiK3RoUwh4wM+OYbzn7xBUuTk1kMnASCgA+BEWFheEyYAEOHQp06ZolUXvG71cGDBwHKvDVakeeN3t7jyclJYv/+OezcqeDl1Z9Ro37H2Tnk/r4JIaqAFEIhTMVohO3bYdEilHXrGFxUxI+AIzAYGOPsTOeRI9GNGQPh4RqHLVt0dDReXl54eXnd8XhFnze2aDGLDz5IYvToTbzxRg6eno/Rt28U9vZV3+IVojLMcms0Pz+f119/HW9vb2rVqkWnTp3Yvn37Xa9bunQper3+jv/S09PNkFyIe3DmDMycCc2bw+OPw+rV6IqKaAV8AaR26cKS5cvpkp6ObsECiy6CAPXr12fAgAFlHq/o80adzoY2bb5n9uwwXFxg8uQU9u7thcGQa5LcQlSUWVqEzz33HD/++COTJk3C39+fyMhInnzySXbu3FnuQ/ob3n33XZo1a1Zin5ubm6niClF5BQWwaZPa8/M//1GfBd7K05N/PPecutZfy5aaRLxXr732WrnHlyxZwujRo9mzZw9du3Yt95arjY0zDz30bz75JJxx49KZPPkAixePoHXrH9DppMuC0IbJC+Hvv//OqlWr+OSTT/i///s/AEaMGEFoaChTp069a7dtgCeeeIJwC/+rWdRQR47A4sUcXLSI2pmZtLj1mF4PPXuqU5717g12dlqlNLmKPm8EcHBoRK9eW3nvvc5MnpzH9Ok/Mm/em7Ro8b4JEwpRNpP/CbZmzRpsbW0ZN25c8T4HBweef/559u3bx7lz5+76GoqicOXKFQwGgymjClExOTmwbBnZXbrwZWAgD3z8MW0zM1lw43iTJuqt0ZMnYfNmeOqpal0E70Xt2m0YNmw1kyfr2LIFtmz5gNTUpVrHEjWUyQthdHQ0AQEBuNw292H79u2Bmz3SytO9e3fc3NxwdnamX79+JCcnmySrEOU6cADlhRfY7eHBc889R8N9+3gZ8AbW29jwzwEDYNs2OH4c3nwTGjfWOrFFc3fvzcsvf8rChRAWBklJ48jM3KV1LFEDmfzWaGpqKg0bNiy1/8a+lJSyJ+N1dnYmIiKC7t274+rqyp9//smcOXPo0qULUVFR+Pj4mCy3EABkZcF338GiReyIjuZF4AjQHHgTGOXvj/eECTBiBLi7a5vVCjVq9Ap/+UsSKSmfoyiFxMU9TXj4PmrVsq7nqMK6mbwQ5ubm4uDgUGq/o6Nj8fGyDBw4kIEDBxZv9+3bl8cff5y//OUvvP/++3zxxRdVH1gIRVEHtC9eDKtXq/N/Au5AOLDAwYFuw4ahHzsWOnUqtdZfRTuOCNDpdPj5fUpe3nEuXfo3RUWZHD7ci/Dw37C3lz8shHmYvBA6OTmRn59fan/e9Q8XJyenSr1e165d6dix412HX0yaNKlUz9KhQ4cydOjQSr2fqEHS0mDZMrUAHj1a6nCrDh34bswYGDwYXF3v+BKVncxagF5vS3Dw90RHP8i1azHk5R0jLu4pWrfejl5f+o9oIe5k5cqVrFy5ssS+7OzsCl1r8kLYsGHDO97+TE1NBcDb27vSr+nj40NSUlK558ydO1d6moq7KyqCrVsp+PprojZvptPtHbLq1lVvez7/PLRqddeXq+xk1pYsKyuLgoICPD09Tf5etrauhIVtIiqqIwUFaWRn7+bIkTEEBi5Hd1uLW4g7uVNDJyoqinbt2t31WpN3lmnbti1JSUlcuXKlxP79+/cDZU/bVJ7jx48XD+AV4p6cOAFvvklCo0b8rXdvGm3YwMMGA8V/Pz7yiPpsMCUFPvusQkUQKjeZtaVbuXIljRo1uuMdHVNwdPQlNHQjer16l+js2RWcPDnTLO8tajaTF8IBAwZgMBj4+uubM87n5+cTGRlJp06daNSoEQBpaWkkJiZSVFRUfN6Nv6xvtWXLFqKioujZs6epo4vqJj8fVq3i2iOPsLR5cx587z2C09NZCowAotzdcXvjDXVB3J9/Vuf9vP4su6IqM5m1pYuOjiYkJOSOz/hNxdX1AYKCvqWgAF59FWbPfpvz578z2/uLmsnkt0Y7dOjAwIED+fvf/056ejotWrRg2bJlnD59msjIyOLzpk2bxvLlyzl58iS+vr4AdOnShfDwcNq1a4ebmxtRUVEsWbIEX19f3njjDVNHF9VFXBwsWoRx+XJeunSJb4ErQA9glU5HvyefxGH8eHWZI9v7/7+ENRe/W93LGoRVwcPjKQIDZ9Gq1VQWLIBGjUbx4otNcHOz3ta1sGxmmWJt+fLlvPnmm3zzzTdkZmbSunVrNm3axIMPPlh8jk6nK/UsYMiQIWzevJlt27aRk5ODt7c348ePZ8aMGXJrVJTv6lVYtUqd8uy33wD19kcO8BoQ4etLsxdegFGj4B6eU1d3RUVFxMTEMHz4cE3ev3HjvzFjRhLnzi1i5swivLx6M2TInzg5tbj7xUJUllLNHDhwQAGUAwcOaB1FXLdnTyNl506UvXt9TPtGRqOi/PaboowZoyguLoqiDoS4+c/BQVGGDVOUHTsUxWAwbRYrFxMTowDKrl27NMtgMBQo+/Z1VwICUNzdUTZsaK4UFFzSLE+Zli27+Tu2YIHWacQtKloPZJZbYf0uXoTPPiM1KIiLnTqprcCrV28eDwuDefPUji/ffgvdu6vzgIoy3W0NQnPQ6+1o124tc+f6odPB5MnH+eOPfhiNBZplEtWTfBoI62Q0ws8/UzR4MJu8vOj/2ms0PnKEz28cr10bxo2D33+HQ4fglVegXj0tE1uV6OhomjdvrvkqL3Z2dejefRsff1yXs2dhypRfSUgYj3L76h5C3AdZmFdYl3PnYOlSjn/5JUvOniUSSAHaAv8ChnboABMmwMCB4OysbVYrdujQIU1bg7dycmpG//6bSUv7CwsXFnH06FJq1w7E1/d1raOJakIKobB8hYWwZQssWsTuzZt5W1H4GXADngWer1OH8OefVwe9BwVpHLZ62LhxI1lZWVrHKObm1pmIiBU88MAQbGzg+PFpODq2wNOz7AWDhagoKYTCciUnq9OdLV2qTn8GXAUKgOXAM3/9K7XGj4e+fcHeXsOg1Y+Tk1Olpz80NU/Pwfj5JXPixD8ASEwcgaOjL66uHTROJqydFEJhWXJzYe1atcPLrl2lDvds3Jieo0dDRIS67p+oUXx93yAnJ4nz55djNOYRE9OXdu324+govwvi3kkhFJbh0CGUhQs5+803NL58ueQxW1vo109d6b1HD7Cx0Saj0JxOp6Nly6/JyztJdvb/KCw8f321ij3Y2mrbsUdYL+k1KkwmISEBf39/unU7x1//Cg8/fBZ/f38SEhLUEy5fhq++4kKbNnzapg1hn39O28uXKZ7ZsmVL+PhjOHsW1qyBnj2lCAr0egdCQ9fi5OQPQE5OHHFxgzEai+5ypRB3Ji1CYRJxcXG0adOmxNyxBgMkJyfTKiyMg717k/qf/7AoP591gAI8Bcy1t8duyBAYOxa6di211p8QAHZ29QkL20xUVCeKii6xd+9WHBxepmXLL2S1ClFp0iIUJtG/f/8SRfBWRQYDbdavp0d+PjHAR8C5Vq1Y9cUX9EhPR79sGTz4oBRBUcro0aNp2bIlo0ePplYtf0JD15GWZssrr8C7737F2bOfaR1RWCFpEQqTOH36dLnHjcBeZ2c6jRqFbswYaNvWPMGE1Spr0eNu3RYzbtwovvgCGjWaxNSpzXF376t1XGFFpEUoqt7p0yhltAZvsLGxoXNGBrrPP7/nInhr60BUjSVLljBs2DCtY9xRWYsee3mN5G9/+wd9+sCcObBixSCuXInSMqqwMlIIRdUoKIAff1SXMmraFJ3RWO7pOhsbuI9xajdaB0lJSWzatEmKYRXZsWMHx48f1zrGHZW36HGzZjN5991BtG0Lb76Zz6ZNPcnLO6tVVGFlpBCK+5OYCFOmgI8PDBgA//kPKAq+d7msqKiIhQsXYrxLwSxLWa0Da2NprdqDBw/S1kJvU5e36LFOpyM0dBmffNIed3f4298y2L37CYqKrpbzikKopBCKyrt2DZYtg4ceUqc0++QTuF6UAGjalA0TJ2JbxiK3tra29O7dm3HjxtGpUyd+//33Skcor3VgLSytVZubm0tiYqLFFkJQi+GRI0fuuPixjY0jnTtvZs4cH3JzYdKkWGJihqAoBg2SCmsihVBUjKLAgQPqhNbe3vDcc7B7983j9vYweDD8979w7BhBn33G4cOH8fPzw85OHf5nZwd+fn4cPnyY9evXs3v3bgoLC+nUqRNjx47lwoULFY5TXuvAWmjVqi2rFRoTE4PBYLCYybbvhb29B489to0PP3SmRw/Izt7MsWN/0zqWsHTmWR7RfGRh3ip26ZKizJ+vKG3alF7oFhQlOFhR5s5VlIyMMl9i7doGyoABKD/95FXqWFFRkfL5558rderUUerWrassXrzYlN+NRYmIiFA8PDwUQPHw8FAiIiI0fc+vvvpKsbGxUXJyckyew9QuXdqu7Nplq+zcibJzJ8rZs5+b7s1kYV6LJQvzinunKPDLLzBihNr6e/lluL5QK6Aub/T887BvH8TGwmuvgbt7mS+XmWlkzRr16+1sbGx48cUXSUpK4plnniEzM9MU35FF0qJVW14rNDo6msDAwBKTbVvaM8yKqlv3Ufz9vyjePnr0FS5e/I+GiYQlk3GE4qa0NPXZ36JF6soPt+vYUZ3vc/BgdeHbKuTh4cHChQur9DWtgblv6Xbt2pXMzEwyMjJKPVs9ePBgiduiZY3bsxbe3mPIzT3KmTOzACPx8YNo23YPLi5hWkcTFkYKYU1XVARbt6rFb+NGdR60W9Wrp7YMn38ewuQDxNotWbKE0aNHs2fPHrp27VqisI0dOxYfH5/i7erQM7d58w/JzU3mwoW1GAxXiInpRXj47zg4eGkdTVgQKYQ11YkTsGQJREaqq77f7tFH1dZf//7g6Gj+fMJkymrV3X77s7zWo7XQ6fQEBX3DwYNnuHLlD/LzzxAb25c2bXZhY1NL63jCQsgzwpokPx9WrVKXMmreHN57r2QR9PaG6dPh2DHYvh2GDLGoIrhnzx6GDRtGSkqK1lFqhOrQMxfAxqYWoaEbcHDw5fJl+PjjPzh8+FkU5d7GsIrqR1qENUFsrLrS+/LlcOlSyWM2NtC7t9r669lTXfvPQmVmZvLzzz/TsmVL3nrrLV599VXsZWV6k7LW4nc7BwcvwsI28c03nVi7NoesrJ+YP38afn6ztI4mLIC0CKurq1fV4te5s/ps79NPSxbBFi3gww/hzBn46Se1GFpwEQTo3bs3R44cYfTo0UybNo3WrVuzfft2rWMJK+HiEsagQT8yZYqOrVvho48+JiVlkdaxhAWQQlidKAr89pu6ll/Dhmor77ffbh53cIBnn4WdO+HoUZg2TT3PxJycdISEqF/LU5Gu+nXq1OGzzz4jOjoaDw8PevTowYABA+662oUQAPXr9+TFF+czapT6d+LXX48nM/NnrWMJjUkhrA4uXlRbfGFhagtw0SK1RXhD69bwr39BaiqsWAHdupl1rb8mTeyYPx98fe3KPKey0421atWKX375hRUrVrBnzx7Cw8O5du1aVUcXJqD12MRGjV7kjTcm8uij8MEHRlat6se1awmaZBGWwbLvhYmyGY1qy27RIli7Vl394Va1a8OwYWqrsF07i1/k9l666ut0Op599ln69OnD/v37cXZ2NnVMcZ8sZWyiv/8cPvroKBER/+bvf7+Gp+dj9O59AHt7T7NnEdqTFqG1OXcO3n8f/Pzgr3+F778vWQS7dlWHRKSmwpdfwgMPWHwRhPubRNvV1ZUePXqYKlq1t3v3brM9a7WUsYk6nQ1t2vzAnDmh1KoFs2efJTa2PwZDniZ5hLakRWgNCgthyxZYuBD+/W+1NXgrd3cYNUod9B4UpE3G+1TeQG9hWp9++imXLl3ir3/9q8nfy5LGJtrauvDQQ/9m7tx2ODikc/nyPo4ciSAo6Ft0Omkj1CRSCC3Z0aPqE/1ly9Tpz26l08Fjj6m3Pvv2VVd/sHKmLH6HDh2iVatW6KygdWxu0dHR9O/f3yzvZWl/8Dg6+vDEE/8mOvohjMYc0tO/x8nJn2bNZmqaS5iXFEJLk5urrvS+aJE68fXtGjeG0aMhIgKaNDF/Pit05MgR2rZty+OPP868efPw9/fXOpLFyM7O5vjx42Zdg1Dr4ne72rXDCQ5eSWxsf0Dh1Kl3cXLyw8trpNbRhJlI+99SHDoEr7yizu4yYkTJImhre3P19xMn4O23pQhWQkBAAOvWrSMxMZHQ0FCmT58uPUyvO3ToEIBVr0FYFdzd+9Kixezi7SNHxpCV9T8NEwlzkkKopexstUNL+/bQpg3Mnw9ZWTePBwaqq7+fOwerV8Pjj6szwYhK0el09OvXj/j4eKZNm8bs2bMJDAxk9erVKIqidTxNRUdH4+DgQGBgoNZRNOfj8xre3hMAUJRCYmOfIifnqMaphDlIITQ3RVFXdn/uOXUw+4QJ8OefN487Od1c/T0+HiZPBk/r7tKdnFzIwIFw7FihpjmcnJx45513iI+Pp23btgwaNIgePXqQfKclp2qI6OhowsLCsLXwWYXMQafT4ec3j7p1HwegqOgSBw8+SWHhRY2TCVOTQmgu6elq6y4oCB56SO0Ak5t78/gDD6itw9RUdfhD165WMeyhIoqKFC5cUL9agubNm7NhwwY2bdpESkoKhtuXnroLrQeEV6WDBw+a9fmgpdPrbQkJ+QFn51B27ICIiGT27++L0Vhw94uF1ZJCaEoGg/pcb8AAaNQIpkyBI0duHq9TR139PToa/vgDxo8HNzft8tYwvXr1IjY2lpYtW1b4msrOgGPJjEYjDg4OdOjQQesoFsXW1pWwsE34+9fn7Fn429/2Ehf3fI2/jV6dyf0QUzh1Sm3VLVmiTmp9u27d1GEPTz+t3goVmtHrK/e3oKUMCK8Ker2e/fv3ax3DIjk6NqFfvy2cP/8gr79eyNtvr2DOnECaNJmudTRhAlIIq0pBAWzYoA572LZNfRZ4qwYN1CEPo0eDdN+3WpY0IFyYlqtrB0aO/I6UlIHMnQs+Pv9g+nQ/PD0Hax1NVDEphPcrIeHmWn/XWwrF9Hp48km19ffkk2BX9qTTwjK99957+Pj4MHLkSPR6vcUNCBem5ek5gNde+4izZ6cxfz54ew/nhRd8cXPrrHU0UYXkGeG9uHYNli6FBx+E4GCYPbtkEWzWTF39/fRp2LgR+vWTImiFFEXh6NGjRERE8OCDDxIdHQ2oA8KPHDkiRbAaKa8DVOPGU5kxI4JOneCdd4r46ade5OYe1yClMBWTF8L8/Hxef/11vL29qVWrFp06darwBL9ZWVmMGzcODw8PXFxceOSRR4o/jMxOUdRhDi+8oA57iIiAW58P2dvDkCGwfTskJ8P06WoHGWG1dDody5YtY+fOnVy+fJkHHniAF198kUu3LnAsrN7dOkDpdDqCgr7kk08extcXYmIyiYnpTWFhVhmvKKyNyQvhc889x9y5cxkxYgTz5s3DxsaGJ5988q6dDIxGI7169WLlypVMnDiRWbNmkZ6eTrdu3cw77iszUx3o3ratOvD9q6/gypWbx0NC1LUAU1Jg5Up49FH1lqgo5ulpw4svgoeHdU4G0K1bN6Kjo5k9ezbffvstAQEBLFy4EOPtk58Lq1SRDlB6vT0PPLCOJUta0rs35OQkEBc3AKNR27GxooooJrR//35Fp9Mps2fPLt6Xl5en+Pn5KV26dCn32lWrVik6nU758ccfi/dlZGQodevWVYYNG1bmdQcOHFAA5cCBA/ce3GhUlJ07FeXZZxXFwUFR1PbgzX/OzooyZoyi/Pabeq4o1549jZSdO1H27vXROkqlREREKAEBAUpERETxvtTUVGXkyJEKoEyYMEHDdKKqREREKB4eHgqgeHh4lPjvfbucnGRl9253ZedOlJ07URITxyjGpUtvfjYsWGDG5OJuKloPTFoIp0yZotjZ2SlXrlwpsf/DDz9UdDqdcvbs2TKvHThwoNKwYcNS+8ePH684OzsrBQUFd7yuvG88Pj5eedjPTwm1t1fC7eyUUHt75WE/PyU+Pl49ISVFUT78UFH8/EoXP1CUTp0UZdEiRbl8uRI/BWGNhfBuH467d+9WYmNjNUp3/3Jzc7WOYFHu9EdPWbKydiu7dtkrO3eiLF2K0rS+s2IPih0o9jY2it+tnylCUxUthCbtNRodHU1AQAAuLi4l9rdv3x5QZ7VoVMZztOjoaMLDw0vtb9++PV9//TVJSUmEhIRUOEtcXBzPtGnD0qIiOgI6wAjsT07mqbAw1j38MEG//KIOgr9VvXowcqS61l9oaIXfT1i3u90us/ZhE08//TT16tVjxYoVWkexCJXp+OTm1pXAwKVs2TKMsWPBYLhlAneDgeTkZFq1asXhw4cJstL1QWsakz7MSk1NpWHDhqX239iXkpJikmvv5KX+/YksKqITahEE9ZvvDCwxGJiwY0fJInhj9feUFJg7V4pgDdO1a1c8PDwAquV4wejoaJo1a6Z1DKvVoMFQ3nmnbqm/m28oKiqib9++5g0l7plJW4S5ubk4ODiU2u/o6Fh8vCx5eXn3fO2dXDx9mk5lHOsEXAR1CaQba/01b16p1xfVy/2OF7x27RpXr16lQYMGJkp479LS0khLS6vxSy/dr9TU8pfyOn36tJmSiPtl0kLo5OREfn5+qf15eXnFx01xLcCkSZNwu2XezjOFhXwPDL3DuXrA3sZGnRpNZuEX193POMHZs2cze/ZsZs6cyUsvvWRRqzscPHgQQCbbvkdGo5GEhIS7TtauyNykZrVy5UpWrlxZYl92dnaFrjXp/zsbNmx4x1uYqampAHh7e5vkWoC5c+eWeMYY5uDAkII7zyBvBAoUBc6fl7F/okq8/PLLpKWlMWnSJBYtWsT8+fN5+OGHtY4FqLdFXV1d5dZoJf30009ERkaye/fuCo0l1VWT1WOsxdChQxk6tGRTJyoqinbt2t31WpM+I2zbti1JSUlcuXXcHRRP9FverZk2bdoQFRVV6q+q/fv34+zsTEBAQKWy1Pf15bcyjv0G1DcaoWVL+OgjuENLVNy7nBwjsbHq15qiXr16LFiwgD///BMXFxe6devGsGHDKv1s2xQOHjxImzZt5IO6ks6cOcPVq1d55ZVX2L59O83v8vikcePGZkom7pspu67eGEf4ySefFO+7MY6wc+fOxftSU1OVhIQEpbCwsHjfjXGEa9asKd6XkZGh1KlTRxk6dGiZ71lWd9n4+Hilpa2tsgcUw/XhEAZQ9oDSEpT4W4dJtGihKBs3yhjBKrJ4sToMITLSU+somjAYDEpkZKTi6empuLi4KLNmzSrxu14VKtP939/fX3n11Ver9P2t2YULF5SffvpJOX78eKWui4+PV2xtbRWg1D8bG5QdO94wUWJRURYxjlBRFGXQoEGKnZ2dMnXqVOWrr75SunTpotjb2yu//vpr8TmjRo1SdDqdcurUqeJ9BoNB6dy5s1K7dm1l5syZyueff66EhIQobm5uSlJSUpnvd0/jCPftU5SXXlIUvb7kuMEnnlCUxMSq/YHUQDW9EN6QmZmpTJw4UenYsaNSVFRUZa9bmQHhly9fvv7fIrLK3r8qVaag36tTp04pK1asUMaPH68EBwcXF6/58+dX+rXi4+MVP0/PW8YR6hVvb3V84c6deiUjY6MJvgNRURZTCPPy8pQpU6YoDRs2VBwdHZWOHTsq27ZtK3HOc889p+j1+hKFUFHUD44xY8Yo7u7uirOzs9K9e/e7fkP3NbPMoUOK8vDDJYuhnZ2iTJmiKNnZlX89oSiKFMLblTUZxL0KCAgo0RoJCAgo89zCwkLljz/+UNLT06s0Q1WoTEG/F5MnT1YaN25c/HMKCgpSxo0bp3zzzTfKyZMn7/2Fly0rMbPMsWNvFM8888svzsrly9FV902ISrGIAfUADg4OzJo1i1mzZpV5TmRkJJGRkaX216lTh4ULF7Jw4UJTRrypVSvYuRNWr4bJk+HsWSgshI8/hm++gX/+E4YPl7lExX2xq+KVSCqzRqKtrS0PPPBAlb5/VTH1ose1a9dm4MCBPPTQQyXGiVa1Zs3eJTc3mYyMHzAarxET05t27X7HwaH8Dn5CO/KJfjudDgYNgsREePNNuDGWMS0NRo2Crl3hjz+0zSjELZYsWULv3r0JCAigd+/e9zXso7zliEztXiYxyMnJYceOHbzzzjulOuXdbsaMGcyePZv+/ftXeRGMBQ5d/986nZ7AwKW4uqojlwsKzhET06fkDDTCokghLIuzM8ycqS68+9RTN/f/9ht07Kgutpuerl2+akDLD11LlZ+fz4QJEzh+vHLr3VXFGol3W47I1CpS0C9evMiGDRuYMmUKnTp1ws3NjUcffZTPPvuMo0ePmjXvrWYA027ZtrFxIjT0JxwdmwJw9WoU8fHDUJTyxx4KbUghvJtmzWDtWti2DW7MG6go6qr0/v7qEkyFshRLZWn9oWupjh07xqZNmwgODmbGjBmVnkHpfpj61mRFlFXQCwoKaNWqFe7u7vTr14/vv/+e5s2bM2/ePGJiYrhw4cId5ybWkr19A8LCNmNj4wrAxYsbOHZsqsapxJ1IIayoHj3g0CF13lFX9Reby5dh0iRo3Rr++4dAY0YAACAASURBVF9t81kZS/jQtUTBwcEkJibyf//3f3z44YcEBwezfv16s8xSYsnzq9rb2/P000+zfPlyTpw4wenTp/nuu++YMGECoaGh6C30ub2zczAhIWsAdS3Os2fncO7cl9qGEqVY5m+PpbKzg9deg6NH1dUobgxITkiAxx5Tb6FW8pZWTdCihR0//ADNm9/sJGLJH7pac3Z25oMPPiA2NpaWLVvSv39/nnzySZPf+qvKZ40VUVRUxB9//MGcOXOYMGHCXc9/++23GTFiBE2bNrWqyQDq1etBQMCC4u2jR1/m0qWtGiYSt5NCeC88PWHRIti/HzrdMpX3Tz9BcLDayeaaPBi/wc5Oh4eH+vUGc37oWuuzyICAAP7973+zbt06EhISCA0NZe/evSZ9z6p41liW3Nxcdu3axcyZM+nRowd16tShQ4cOTJ8+nYSEhOJ5hC1NVfz+eHuPo3Hjv13fMhAXN4irV2OrJqC4b1II70f79rBnDyxbBl5e6r78fHjvPQgMhFWr1OeJ4o5M+aF7g7U/i9TpdPTv35/4+HhmzZpVvJbnvZg+fTo//PBDFaaruIMHD+Lm5kb37t2ZO3cujo6OvPXWW+zdu5fs7Gx27dpVvLKMJanK35/mzT/C3b0/AAbDZWJielNQcL6qoor7IIXwfun16sK9R47AlCnq7VNQxyAOGQLdu8Phw9pmrMGqy7PIWrVq8eqrr97zGERFUfj8889JTk6u4mQVExgYyKeffsrhw4e5ePEiGzduZOrUqXTu3Bl7e3tNMlVEVf7+6HQ2BAWtwMVFnQQ6P/8UMTF9MRjM1yFK3JkUwqri6gqzZkFMDPTseXP/L79A27bw8stQgRnrRdWSZ5GqkydPkp2dXWVLLymKQmJiIgsXLmTkyJFMnDix3PMdHR158cUXCQsLq1THFq1va1f174+NjTNhYRtwcPAB4MqV30lMHImi1JwJ6S2RFMKq1rIlbNkCGzdCixbqPqMRPv9cHW7x5ZeUuay1qHLm7gCipfPnz5fZuzQ6Ohr+v707D6/pWh84/j0nM2IIEUFRkpBRUXNrak2NeYrUUBQdlNatoqVFKUUvqqbqgJZrbE0drlLE+DNrQhIxa0sSaooQMuzfH1sSrgwnyTlnn+H9PE8ePWvvs/eb3eS8WXuv9S7yXvElL2lpaRw+fJjZs2fTrVs3vLy88Pf35/XXX+fkyZOULVu20HHnxhJuaxv687MEWGXgMV1cKhIc/BMODiUAuHp1HefPjzdOwKJQLGe1UFui00GHDuqUi9mz1WeGyclqj/CNN9Rk+MUX8PzzWkdqF6wx+Q0aNIi9e/fStGlTg+JPT0/nhRdeoGzZssybN4/g4ODHth8/fhwvLy+8vb0LFc/ixYsZNmwYLi4uNGzYkKFDh/L888/TuHFjSmZOJzIyS7mtbcj1L+gVKFGiNgEBq4iK6gRkcOnSNNzcfPH2HlioGEXRSI/QlFxcYOxY9flhnz7Z7X/8Ac2awcsvq88ShXhEYXpCDg4OzJo1i4SEBOrUqcM777zDpEmT0Ol06HQ6Jk+eTEJCAjqdLse6v/nNU+zatSt79uzh1q1bREREMGXKFNq2bWuyJAi2f1u7bNlQfHzmZL2OixvKjRs7NIzIfkkiNIdKlWD5cti9W31emGnlSvVW6tSpYKFDx40hISGN+fMhMTFN61CsQmF7Qm3atCEyMpKpU6eyYMECJk6cmON+Y8aMYdy4caxatYphw4YREhLClClT8jy2t7c3TZs2xSWz9q4Z2MNt7cqVh1Op0nAAFCWNkye7kZwcq3FU9kcSoTk995xasPvLLyHzmcrduzBuHAQGwqZNNjnd4saNDNatU/8V+StKT8jZ2ZnRo0eTmk/Zv6lTpxIeHs7WrVupX78+DRs2LFLMpmKOKTZaq1FjFh4eLwGQlnaTqKhQHjy4pnFU9kUSobk5OMDQoWp1muHD1degVqTp3Bnat1dXvhB2y1w9oStXrhAXF8c333xDmzZtTHIOkT+93pGAgFUULx4CQErKOU6c6EJ6uu3eJbI0kgi1UqYMzJ0Lx46pcw0zbdkCwcHqeoi3bmkXn9BUYXtCFy9ezHPtz0dVyCwCITTn6OhOcPBPODurg5lu397LqVOvmqXGrJBEqL3gYPj9d3Ux4CpV1La0NJg1C/z8YMkSdfqFEHnYuHEjTZo0oVq1akyYMEHrcEQhuLo+RXDwZvT6YgAkJv6HCxcmaRyVfZBEaAl0OujRQy3ePWECZJaaSkyEQYOgcWO1rqkQubh9+zblypVjxYoVJMo6mWa3GJhnhOO4u9fD338FoNblvXhxEgkJK4xwZJEXSYSWpFgxmDhRTYjdu2e3HzyoFvceOBDi4zULT1iufv36sWnTJl5++WXc3d2ZPn16vu9ZuHChGSKzD1uAn410LE/PLtSoMTPrdWzsIG7e3GOko4ucSCK0RNWqwbp1sG2bOpo009Kl6u3Sf/8bHjzQKjphRsnJyaxevZq1a9cW6H2jR4/ONRl++umnvPPOO7z55pvMm2d4P0brcmf2pHLlf+HtPRQARXnAiRNduHtXmzqx9kASoSV74QV1MM3nn0OpUmpbUhKMGgUhIerAGivg5qYjMFD9V+QvJSWFDRs20Lt3b8qXL0/v3r3ZsGFDgY8zevRoFEXh/fffp1q1aiiKgqIojBkzhlmzZjF27FgqVapk0LEsodyZPdHpdPj6zqNMmdYApKX9Q1RUKKmpNzSOzDZJIrR0Tk4wYoQ63WLIkOzFgE+dUot7d+4MZ89qG2M+qlZ1Yt48qFKlcCsn2IPU1FT++9//MmDAALy8vOjatSsxMTGMHz+es2fPsmJF4Z8TJScnU7x48cfadDod06ZNo2vXrgYdw1LKndkTvd6JwMC1FCsWAMC9e3GcPNmNjAy5G2RskgithacnLF6sTshv3Di7fdMmdTHgcePgzh3t4hNFsnXrVtq3b8/+/ft55513OHnyJH/88Qfvv/8+1atXL9KxnZ2dC11jNJOtlzuzVI6OpQgO/hknp/IA3Ly5k7i412VahZFJIrQ29eqpiwEvXw6ZH24PHqhl2mrVUsu2yS+J1XnxxRc5evQosbGxTJo0iYCAAKMde+bMmWzdurVIx7CHcmeWys2tGkFBG9Hr1dHk8fFLuHTpU42jsi2SCK2RTqcW8T51Si3qnblY699/q4W8mzWD48e1jVFkURSFy5cv57mPs7MzderUQaez3Oeo9lDuzFKVKtWIWrWWZb0+f/4DEhMLNoBK5E4SoTVzd4dp0+DkSXXZp0x79qg9xzfegGtSs7CwijpKMiYmhgkTJlCrVi2eeeYZ0tKMW3TclKM4b9++LbffCiAQCDHxOcqX78XTT3+S9To2tj+3bv2fic9qHyQR2gJfX3Uh4J9/Vv8b1Go0ixap0y3mz1er1QiDFXaU5Llz55g2bRq1a9cmICCAOXPm0KRJE77//nuj9vZMOYozJSWFxo0bM2rUKEmGBvoYyH/mZtFVqfI+FSoMACAjI4UTJzpz794FM5zZtkkitCUvvQQnTsCMGVBCXf2aGzfgrbegbl3YuVPT8KxJQUdJXr9+nYYNG1KjRg2mTJlCQEAAGzZsICEhgSVLltC2bVscMgusaxBfQbi6uvLmm28ya9Ys3nnnHUmGFkSn0+Hn9yWlS7cAIDU1kaioUNLSpC5xUUgitDXOzvDeexAXB/37Z7dHRanFvcPC4NIl7eKzEgUdJVmmTBnq1q3LqlWrSExMZOXKlXTu3BnXzHJ5GsdXUMOGDWPRokXMnTuXYcOGkSH1bi2GXu9MYOAPuLn5AXD3bjQnT/YkIyPvpbdE7iQR2ipvb1i2DPbtU58XZlqzRh1dOnky3LtnllDOnEmlZ084e9Z6flELOkpSp9OxcOFCwsLCnpizZwnxFcZrr73G119/zaJFi3j99dclGVoQJycPgoN/xtHRA4AbN7Zy+vRw6b0XkiRCW9e4sVqr9Ouv1bmIoCbAjz5S5x+uX2/y6RZpaQrXrqn/WpNvv/2WI0eO8OKLL9K1a1duWdiyWMYcxZnbwJtXX32VJUuW8PXXXzN48GDS09OLfC5hHMWK+RAUtAGdzhmAK1e+5K+/ZmsclXWSRGgP9Hp49VX1dunbb2cvBnzhAnTrBm3aQHS0piFaknv37vHjjz/Sq1cvypcvT58+fUhISODKlStah1ZgiqLw3HPPsXnz5lz3yW/gzSuvvML333/PsmXLmDp1qqlDFgVQuvTz1Kz5Tdbrs2dHcfVqwcvx2TtJhPakdGmYMwf++EOtY5pp2za1dunIkXDzpnbxaezXX3+lf//+eHl50b17d86cOcPEiRM5f/48+/bto1atWlqHWGAPHjxg79693LiRe41KQwbe9OnTh02bNvHWW2+ZLFZROBUq9KVq1Y8evlKIielDUtIRTWOyNpII7VFgIGzdCj/8oK50AZCeriZJPz/45htNFgPWenWDadOmcejQIUaNGkVsbCxHjx5l9OjRVMu8RhYQY0ElJycD5Pnc0tCBN6GhoZQpU8b4QYoiq1ZtIuXLhwOQkXGXqKiOpKT8qXFU1kMSob3S6dTbotHR8PHH4Oamtl+9CoMHQ8OGsH+/2cKxhNUNNm/eTHR0NB999BE1a9a0yBgL6u7du0DeiVDKpxXdACBMw/PrdDpq1vyWkiXVP2IePLhCVFRH0tKSNIzKekgitHdubvDhhxAbC716ZbcfPgxNmsArr4AZno2Zcl6coigcPnyYhISEPPcrVapUnpPerXEFBkN6hGA75dO06rEnAbfNesYnOTi4EhS0HldXtUh7cvIfREeHk5EhxTTyI4lQqKpUgdWrYccOCA7Obv/uO/V26cyZJl0M2BTz4k6cOMH48ePx9fWlfv36RVrKyFQxmpqhidAWWGOP3dicnT0fTqsoDcD16z9z9uy7Gkdl+SQRise1aAFHj8K8eZD5POjOHRg9Wk2Qv/5a4EOWL+/Am2+Cp2fulVWMdXvu9OnTTJkyhaCgIIKDg5k/fz7Nmzfnt99+Y8SIEYU6prFjNKfMRFisWDGTnmfLli3cM9O81NxYY4/dFIoXr0Vg4A/odI4A/P33XP76a57GUVk2R60DEBbI0RGGDVOr0Hz4IXz5pTrXMC5OLePWoQPMng0+PgYdzsPDgZ49wcUl7xJjRU0sY8aMYcaMGRQvXpwuXbrw6aef0qZNG5ydnYt03EdZQ/J7lDl6hNeuXaNHjx40atSIjRs3mjzp5qZp06bcuHGDq1evWk2P3VTKlGmFn9+XnDr1KgBnzryNm1t1ypZ9SePILJP0CEXuypWDhQvhyBF47rns9p9+Ukeevv++RS0G3KVLF9auXUtiYiLLly+nQ4cORk2C1qhy5coMHz7cpKM9y5Urx88//8z+/fsJDQ3ljkY/E9bYYzclb+9BVKky9uGrDKKjw7hzJ1LTmCyVyRPhzZs3GTp0KJ6enpQoUYJWrVpx7Ngxg947ceJE9Hr9E19umSMchXnUqQO7dsF//gOVKqltDx7Ap59CzZqwYoVZFgPOr3xU48aN6dGjh2Y9EksUFBTE3LlzKZFZhN1EmjVrxpYtWzh8+DDt27cnKUmb0Yq2MujHWJ5++hM8PXsAkJ5+h6ioDty/b32FIUzNpIkwIyOD0NBQVq5cyYgRI5gxYwaJiYm0aNGCM2fOGHycRYsWsXz58qyvpUuXmi5okTOdDsLD1dGlH3ygFvcGuHwZ+vZVe4xHjxr9tLdv3+b7778nNDSUsWPH5v8GoZmmTZuydetWIiMjadu2Lbdvaz2OUuh0emrV+g539wYA3L//J1FRHUlPT9Y4Msti0kS4bt069u/fz7Jly/jwww9588032blzJw4ODkyYMMHg4/To0YOXX3456yssTMsZO3auRAn45BN1/mGnTtnt+/bBs8/Ca6+pcxGL4O7du6xdu5bu3btTvnx5+vfvz+3btwl+dDSrsEiNGjVi27ZtxMTE0Lp1a27aSaWiIcAwrYPIhYODG0FBG3FxqQLAnTtHiInph6JIEfVMJk+EFSpUoFu3bllt5cqVo1evXmzcuJHUVMNWI8jIyJAVsy1NjRqwcaM6ijRz8rmiwOLF6nSLL74o8GLAkZGR9OnTh/Lly9OrVy8uXbrEJ598wqVLl9i9ezd9+/Y1wTcijK1+/fr8/vvvnDlzhl8LMcrYGrUDOmgdRB5cXCoQHPwzDg7uAFy7tp5z5+QOSyaTJsJjx45Rt27dJ9rr16/P3bt3iYuLM+g41atXp3Tp0pQsWZJ+/fqRmJho7FBFYbVrB5GR8Nln4K7+knHzJowYQUytWvhWqkSLFn/z4ovQvPlf+Pr6EhMTk+Oh7ty5Q1RUFB988AGnT5/m0KFDvPvuuzz11FNm/IaEMdStW5dTp04RHh6udShGZ21l9jKVKBFEYOBaQB29/eefM7l8+Sttg7IQJk2EV65cwdvb+4n2zLbLly/n+X4PDw+GDx/O4sWL+eGHHxg8eDCrV6/m+eef1+xhvMiBszO8+646vWLAAABOAiFnz3Lm8mVSU9VSpqmpcObMGUJCQnJMho0bNyYyMpIPPvgAHwOnZgjLVa5cOa1DMDprn7Tv4dEWX98vsl7Hxb3B9evbNIzIMhg8j1BRFO7fv2/QvpmrcqekpODi4pLr9vwm4P7vBOiuXbvSoEED+vTpw4IFCxgzZoxB8QgzqVABliyB11+nS7NmpOVSiSYtLY1OnTpx+vTpx9rzKm8mhCXIcdJ+ixbaBlVAlSq9wb17px+uXZjOyZM9qFt3H8WLB2gdmmYM7hFGRERQrFgxg74yb3m6ubnlmDxTUlKythdUeHg4FSpU4Pfffy/we4WZNGzIpXx2uXQpvz2EMSQmJnL9+nWtw7AZ1lhmLyc1asykbFl1sFt6+i2iokJ58MB+HzkZ3CP09/c3eNpChQoVAPUWaE63PzMXOK1YsaKhp39M5cqV8/3lHjlyJKVKlXqsLTw83CafWVii/AY2KampcOIEBAUV+NiDBg1i7969NG3aVOaL5WPQoEE4ODiwceNGrUOxCd9+++2TP3/ffad1WAWm0zng77+C48ebcefOMVJSLnDiRGdq196Og4N1ztNeuXIlK1eufKzt1q1bhr1ZMaGePXsqFSpUUDIyMh5rHzJkiFKiRAnlwYMHBT5mRkaG4unpqbRr1y7H7UeOHFEA5ciRI4WKWRiHs7OzAuT5tUKnUzKGD1eU69cNPu7AgQMVT09PBVA8PT2VgQMHmvC7sH4tWrRQwsPDtQ7jMcnJycqAAQOUixcvah2KcSxbpijqmGlFWbBA62gKJCXlL2Xv3krKjh0oO3agnDgRpmRkpGsdltEYmg9MOlimR48eJCQk8OOPP2a1Xbt2jbVr19KxY0ecnJyy2i9dukRsbOxj77+aw3y0hQsXcu3aNdq1a2e6wEWRValSJc/txYA+isLIL75Qp1t89ZU6oiYfUli5YO7evWtxK0/8888/7Ny5k+bNm3PhwgWtwzGKE8AfWgdRCC4ulQgO3oxer/6MXL26mvPnP8rnXbbH5ImwUaNGDBw4kMmTJ7NgwQJatGiBoihMmjTpsX379+9PQMDjD2urVq3KoEGDmDVrFgsWLODll19m+PDh1KlTh9dee82UoYsi2rRpE46OOd95d3R05PDbb7PVxYVXAa5dg6FDoUEDyCex2cozGnNJTk62uET41FNPERERgYODA82bN+fs2bNah1RkEwBrnZXn7l6HgICVZKaDS5c+IT5+mbZBmZlJE6Fer+eXX34hLCyMuXPnMnr0aMqXL8/27dvx9fV9bF+dTvfEqMG+ffty8OBBJk2axMiRIzly5Ahjxoxh165dWSNPhWXy9/cnMjISHx8fnJzAwQGcnMDHx4fIyEj858zhxTNnCH70me3Ro2qptn791NJtOZDCygVjiYkQ1DsGERERuLi40KJFiydGEGey1jl71qZcuY74+MzKen3q1BBu3NipXUDmZp47teYjzwgtT+YziH37Kue8Q0SEooSEZD9nAUUpXlxRPv1UUVJSzBuskQ0cOFDx8/PT7Fmmp6enMmXKFE3ObYjLly8rtWrVUry9vZXY2NjHtlnN8+Bly5RuoLSzwmeEj8rIyFBOnXoz63nh7t1llOTkU1qHVSQW8YxQCIM0a6Yu9bRgAXh4qG3JySSPHculmjXh55+1ja+QLGHytaX2CDN5e3uzc+dOPDw8aN68OdHR0Vnb5Hmweel0Onx8PsfDQx1/kZZ24+G0imsaR2Z6kgiFZXB0hDfeUKvTvPkm6PXMAWpdvMjHHTpwr21bdZsVMecHeU63EDMyMixysMz/8vLyYseOHdSvX/+xJbTkebD56fWOBASspnhxtcD9vXtnOHmyGxkZhhVTsVaSCIVlKVsW5s+Ho0cZ0aQJI4ApQMBvv7E+IADlvffASpb3MdcHeV49z4MHD9Lp0VVCLJSnpyebN2+mWrVqWW3yPFgbjo4lCQ7+CScnLwBu3drNqVODbXrRA0mEwjLVro37nj18umoVJ7y8CAC6pafT5rPPiK5RQ53EnGHZy8iY64M8t56nXq+nfv36eHl5meS85iAL7WrD1bXKw2kV6uT6hITlXLw4ReOoTEcSobBcOh2EheF39iw/jx/PT46OXABCrl3jnVde4V7jxnD4sNZR5skcH+RyC1GYQsmS9fH3X571+sKFj0hIWJnHO6yXJEJh+YoXh8mTCT11ihMdO/IJ6uRll4MH1bmHgweDHS/NVdCep6VOSbDUuAyxBFildRAm4OnZjerVp2e9jo0dyK1b+zSMyDQkEQrrUb06Lps2MWbLFrbXrKn+8CoKfPONWp1mzhx1rSc7ZGjP0xJGsuYkv7iSk5M1iswwJYFS+e5lnZ566j28vQcDoCj3OXGiM/fundM4KuOSRCisT5s26KKiYNYsKFlSbbt1C0aOhNq1YZusr5YbS52SkFdchw4d4umnnyYiIkKr8OyaTqfD13cBpUu/AEBq6jWiokJJTb2hcWTGI4lQWCcnJzXxxcXBoEHq80SAmBho3Rq6dYPz57WN0QJZ6vPEvOIKDAwkJCSE9u3bs337dq1CtGt6vROBgesoVqwWAHfvxnLyZA8yMmzjDowkQmHdvLzUW6MHDkDDhgDcB15av55f/Pzgo4/g7l1tY7QgljolIa+4ihUrxubNm2nWrBmhoaH89ttvGkZqv5ycShMc/DNOTuofLDdvbicu7g2bmFYhiVDYhvr1Yd8+WLaM6+XKcR8ITUujw+TJnK5RA9asUZ8n2pHjx4/z2WefPfFBZalTEvKKy83NjQ0bNtCqVSs6derEr7/+qkGEws2tOkFBG9DpXACIj/+GP/+cqXFURSeJUNgOvR7698f77Fm2jRrFDw4OnAAC4+MZGxbGnWbNIDJS6yjNZs+ePXzwwQdPFLO3Vq6urvz444+0bduWLl26sHnzZq1DskulSjWhVq2lWa/PnRvD1as/aBeQEUgiFLanZEl0M2fS7eRJYlq3ZjzwOVBzzx5WPPMMyrBhcP261lGanDWUVysoFxcX1q5dS4cOHejevTuRdvSHjSXx8upNtWqTs17HxPTl9u2DGkZUNJIIhe2qWRO3LVv4aPNmYqtUoSnQV1FYtWAB+PrCokUGLQZsrSy94HZhOTs7s2rVKhYvXkxwcLDW4QCwGJindRBmVrXqOLy8+gOQkZFCVFQnUlIuahxV4UgiFLZNp4MOHagaF8eaadPY7epKD1B7hG+8AfXqwe7dWkdpEraaCAGcnJwYMGCAxdz23QJY5xophafT6ahZczGlSjUDIDU1gaioDqSlWUct4EdJIhT2wcUFxo7luTNncOrTJ7v9jz/UZaBefhn++ivPQ1hb5ZPk5OTHVnMQwtj0eheCgn7Ezc0HgOTkE0RHh5GRkaZxZAUjiVDYl0qVYPlytRdYp052+8qVULMmTJ0KKSlPvM1SK7LkxZZ7hMJyODmVJTj4Zxwd1bVEr1//L2fOjLCqaRWSCIV9eu45OHQIvvxSXfoJ4O5djo0bxwU/P9i06bHpFpZakSUvWg6WsbbesyiaYsX8CAr6EZ3OCYDLlxfy11+faxyV4SQRCvvl4ABDh8Lp0zB8ODg4MBbw//NPJnbuzN02bSA2FrDciix5KVOmDFWrVjX7ebXuPaenp7NNyuyZXenSzalZ8+us12fP/otr16xjioskQiHKlIG5c+HYMX54/nlGAtMA/23bWBcYiPLuu3w7Z45FVmTJy1dffcXixYvNfl6te8+rV6+mdevWzJ8/36znFVChQn+qVh3/8JVCdHQ4SUnHNI3JEJIIhcgUHEyJiAimrl3LSW9vagM9MzJ4YdYsop5+mm+bNeNUTIxVJEEtad17Dg8P59133+Wtt95izpw5Zj23gGrVJuHpGQZARkYyUVEduH//b42jypskQiEepdNBjx74nDnDpgkT+MXJib+BOtevM3bgQGjcGA5a78Rhc9C6nqlOp2PmzJmMGTOGkSNH8tlnn5n8nIFAiMnPYh10Oj21ai2lZMnGADx4cJmoqI6kpd3ROLLcSSIUIifFisHEibSPiyOqa1c+BbxBTYING6orXiQkaByk5dK6nqlOp2PatGmMHz+e9957j2nTppn0fB8D0/Pdy344OLgSFLQBV9enAbhz5xgxMS+jKJZZwEISoRB5qVYN5x9/ZNS2bbwdEJDdvmSJuhjwrFnw4IF28Ylc6XQ6Jk+ezMSJE/nggw/4+OOPtQ7Jrjg7lyc4+CccHNQli//5ZzNnz47SOKqcSSIUwhAvvADHj8Pnn0Oph2uR374N776rLgYsSwNZrAkTJjBlyhRmz57N5cuXtQ7HrhQvHkBg4Dp0OkcA/vprDn//vUDjqJ4kiVAIQzk5wYgR6nSLIUOyFwOOjYW2bUnv3BnOndM2RpGjcePGERMTQ8WKFbUOxe54eLyIr2928jt9egT//PNfDSN6kiRCIQrK0xMWL1YnVYN+YAAAG1NJREFU5DdWBwT8H+C/aRM/1ayJMm4cJCdrG6N4QoUKFQr1PikOUHQVKw7hqafee/gqnejoXty5E6VpTI+SRChEYdWrB3v3wvffU8bTk6pAx7Q0QqdO5VT16rBqlWaLAV+9epXatWuzZ88eTc5vK7QuDmBLqlf/lHLlugGQnp70cFpFvMZRqSQRClEUOh307UvNs2f5bfRo1js4EAMEJyYyOjyc202bqs8WzSwpKYnIyEgeyECeItG6OIAt0en0+Pt/j7v7swDcv3+JEyc6kZ5+V+PIJBEKYRzu7uimT6dLTAzR7drxEer6dDX37+e7unXJeOMN+Ocfs4WT/PDWrBTdLhqtiwPYGgeHYgQFbcLF5SkAkpIOERPTH0XJ0DQuSYRCGJOvL26//sr4n34itlo1mgNDFYU/Fy1SFwOePx/STL9EjSRCw0VHRzNq1CgyMp78MDa0OMAAIMy0YdoMFxfvh9Mq3AG4du0Hzp37QNOYJBEKYQqhoVSJjWXV9OnEFStGVYAbN+Ctt6BuXdi5s9CHNmTwRmYilPUI8xcVFcXs2bMZOHAg6elPTvg2pDhAEmB9y9Fqp0SJEAICVpOZgv78czpXrnyjWTySCIUwFRcXGD2aKqdPQ79+2e1RUdCyJYSFwaVLBTqkoYM3pEdouLCwMFasWMGKFSvo378/aWbosQsoW7Y9vr5zs17Hxb3OjRu/axKLJEIhTK1iRfjuO3WEab162e1r1kCtWjB5Mty7Z9ChDB28cfeuOgBBEqFhevfuzcqVK1mzZg19+vQhNTVV65DsQqVKw6hU6W0AFCWNEye6k5wcY/Y4JBEKYS5NmsCBA/DVV1CunNp27x7ffPQR5319Yf36fKdbGDp4Q26NFlzPnj1Zs2YN69evJzw8XEbcmomPz7/x8AgFID39FlFRoTx4cNWsMUgiFMKcHBxg8GC1Os3bb3NXr2cy4P/333zUrRt3X3wRYnL/i9jQwRshISGMHz8evV5+xQuia9eu/PDDD2zevJlevXpx//59rUOyeTqdAwEBKylevDYAKSnnOXGiC+npKWaLQX5LhNBC6dIwZw7FIiM52aIF7wEzgFrbt7MmKAhl5Ei4dSvHtxoyeKN+/fpMnjzZNLHbuI4dO7JhwwZcXFzQZZbREybl6OhOcPBPODt7A3D79j5OnRqEYqaCFJIIhdBSYCDFt29n8g8/EF2pEvWAsIwMWs6ZQ2T16vDtt5DDsH5jkNJhuV+D9u3bs3r1apydnTWKzP64ulYmOHgzer16Oz8xcSUXLkw0y7klEQqhNZ0OunWj+unTrP/4Y7Y4OxOPuhjwwldfVdc//L//M+oppXSYca/BEGCY8UKzW+7u9QgI+A+g9sQvXvyY+PjvTX5eSYRCWAo3N/jwQ9qcPk1k9+7MBJoBHD6sFvceMACuXDHKqaR0mHGvQTugg5HisnflynWmRo3Psl6fOvUqN2/uMuk5TZoI4+PjGTt2LC1btsTd3R29Xk9ERESBjvH333/Tq1cvypQpQ6lSpejSpQvnz583UcRCWIAqVXBet45/7dhBYHBwdvuyZVCzJsycWeTFgC2pdJhWt2gt6RqIx1WuPJKKFV8HQFFSOXGiK3fvnjbZ+UyaCGNjY5kxYwZXrlwhJCQEoEAPn+/cuUPLli3ZvXs348aNY9KkSRw7dozmzZtz/fp1U4UthGVo0QKOHoV586BMGbUtKQlGj4bgYPj110If2tDRp6am5S3awl4DRVFISTHfiEZ7pNPp8PGZS5kybQBIS7tOVFQoqakm+txXTCgpKUm5ceOGoiiKsnbtWkWn0ykREREGv3/69OmKTqdTDh8+nNUWGxurODo6Kh988EGO7zly5IgCKEeOHCla8MJo9u6tpOzYgbJvX2WtQ7FeV68qyuuvK4pOpyjqbEPlDij32rdXlNOntY6u0Pz8/BQg68vPz0/rkPI1depUpWHDhlmfbcqyZVn/T5QFC7QNzsakpt5UDhwIVHbsQNmxA+Xo0eZKevp9g99vaD4waY+wRIkSlC5dutDvX7duHQ0aNKDeI9U4atasyQsvvMCaNWuMEaIQ1qFcOVi4EI4cgeeeA2ASEPDrr2z090cZOxbu3Mna/dy5cyQmJmoUrOGs8fZkmzZtiIuLo3Xr1ty4cUPrcGyao2MpgoN/wsmpPAC3bkVw6tRQo0+rsNjBMhkZGURGRvLss88+sa1+/fqcPXs2q3qGsEwxMTE0buxL375/M3Qo9OnzF40b+xKTx4RxkY86dWDXLvjPfxjk6Ykf0CUtjXbTpxNTvToxM2bg6+ODj48P3t7euLi44OtrudfcUm7RFkS9evXYvn0758+fp2nTpjz97ru4AM6Ay/DhFn29rZGbWzWCgzeh17ty8yZMmrSM5s29ePHFmrzwQiBvvz0wa9BTYTkaKVaju379Og8ePMDb2/uJbZltly9fxtfX19yhCQOcPHmSl156hjFj0vD3V2cIZGRATMwZ2rUL4b//jcTf31/rMK2TTgfh4dTq2JFfP/mEnz77jHfS0gi6ehXGjCFz1qGiKDx48IAzZ84QEhJCZKRlXnNrSH7/65lnnuGrr76iW7duj29IT7f4622NSpZsSNmy83jnncEMHgz+/lfR6a6SkQGxsdF07rybjRv3Z91dKCiDE6GiKAaXG3J1dS1UMI+697AIsYuLS67Hv2dgoWJhfoMHd2H06DQCArLb9HoIDIT33kujT59nWLw4IPcDCMN0h2oNfNn7778I3JtEbkMJ0tLSePHF2mzcGJjroVJSMrh5M++VF7y8nPIc8HbzZhopKbkXAHB11VO6dN4fO/HxeY+ILVXKATc3h1y3m+v7GDEiLtftaWlptGnzDOvXy8+4scyZc5HBg3niMyUgAPr0OcuUKaP5/PMlhTq2wYkwIiKCVq1aGbRvbGwsfn5+hQook5ubG0COyTdzxFbmPsLyJCZeeuwH9lEBAXDz5gPu3Dlu3qBs1J0y8M8USHoReHI5vSwJCal5XvM9e+DDD/M+19at4JjHp8a0abB9e+7bmzaFKVPyPkeXLpDDsoBZxo2DF1/Mfbu5vo/4+LzPER8vP+PGdOYMvPpqztv8/eHHHw8W+tgGJ0J/f3+WLl1q0L4VKlQobDxZPDw8cHFx4UoOE4gz2ypWrJjr+0eOHEmpUqUeawsPDyc8PLzIsYn8OTkp5PYHt16vfgjpdE/29kXhZZDPHRsl72seFKQwc2begxAcHHR59qT69Mmgffvc31+6NOh0eQ9NmD49I89FOKpXzzsGc30f//d/eS/VpORzvUXBODjcz/MzJSHhOp06dXqs/VYu9Xr/l8GJ0MvLi/79+xu6e5Hp9XqCg4M5dOjQE9sOHDhAjRo18lxrbfbs2dStW9eUIYo8pKbqHn4QPLktIwPS051p3lzmYhmTg4ML6em531Z0yIDm8zrCZ59B1aomiaF5c8s4hjliGDYsn+vtID/jxuTqGoiiROf6meLl5cGmTZseaz969Ohjsw5yYzGjRi9dukRsbOxjbT169ODQoUMcOXIkq+3UqVPs2LGDnj17mjtEUQDly1chOjrnbdHR6nZhXFWq5H1NqwCsW6cuBjxpksGLAYuc5Xu989kuCiYoqEGuK5TFxKjbC0unGHtCxv+Y8vCBwMmTJ1m9ejWDBg2iWrVqAIwfPz5rvxYtWrBr1y4yHqm0f+fOHerUqUNSUhKjRo3C0dGRWbNmoSgKx48fp2zZsk+cL/MvgCNHjkiPUEMxMTG0axfCe++pA2b0evWvtuhomDnTUUaNmkBMTAwhISGkpT05UMRRryeydGn8H63IVLUq/Pvf0K1bzl13kac8r7ejo4waNbKrV6/SuXNj+vQ5i79/9mdKTAysWFEjx1GjBueDws/5N4xOp1P0ev1j/2b+96NatGjxRJuiKMpff/2l9OzZUylVqpTi7u6udOrUSTl79myu55PKMpYjOjpaadTIR6le3VmpWdNJqV7dWWnUyEeJjo7WOjSbFR0drfj4+CjOzs6Kk5OT4uzsrPj4PLzmN28qyr/+pSiOjtmVUEBRWrVSlKgorUO3StHR0YpP+fKKMyhOoDg7OGRfb2F0iYmJyogRA5RWrQKUF17wU1q1ClBGjBigJCYm5ri/ofnA5D1Cc5MeoRD5iImBt99Wh05mcnCAYcNg4sTsuqbCMN99B6+8ov73ggXwxhvaxiOyGJoPLOYZoRDCTPz9YcsW2LABnn5abUtPh7lzwc8Pvvoq7/kLQtgYSYRC2COdDjp3Vh/aTpkCxdRVwbl2DYYOhQYNYN8+bWMUwkwkEQphz1xd1RnqsbHQu3d2+9Gj6qzxfv3g8mXt4hPCDCQRCiHgqadg5UqIiICHa4cCsHy5ert0+nQwsMSiENZGEqEQIluzZupSTwsWgIeH2pacDGPHQlAQ/PyztvEJYQKSCIUQj3N0VEc+xsXBm2+qE7ZALfbYoQOEhqrbhLARkgiFEDkrWxbmz1efFzZrlt3+yy9q73DMGEhK0i4+IYxEEqEQIm+1a8POnbBqFVSurLalpsKMGerzw++/V0t8CGGlJBEKIfKn00FYmDq6dPx4yFwnND4e+veH556Dw4e1jVGIQpJEKIQwXPHiMHmyOv+wS5fs9v371bmHQ4ZAYqJ28QlRCJIIhRAFV706rF+vVqipVUttUxT4+mv1dunnn6u3T4WwApIIhRCF16YNREbCrFlQsqTadusWvPMOPPMMbNumbXxCGEASoRCiaJycYORIdUrFoEHZ7dHR0Lo1dO8OFy5oFp4Q+ZFEKIQwDi8v+OYbOHAAGjbMbv/xR7XQ94QJcPeudvEJkQtJhEII48os2L10qZocAVJS4OOP1eeJa9eqzxOFsBCSCIUQxqfXq2v0xcXBqFFqtRqAP/+EXr2gVSuIitI2RiEekkQohDCdkiVh5kw16bVtm92+c6c6mOatt+D6dc3CEwIkEQohzKFWLfj1V9i0SZ16AWo1mvnz1ekWixbJYsBCM5IIhRDmodNBx45w8iRMnZq9GPA//6hFvp99Fnbv1jZGYZckEQohzMvVFd5/H06dgpdfzm4/flwt7v3yy/DXX9rFJ+yOJEIhhDYqV4YVK9Re4DPPZLevXAk1a6q9xpQU7eITdkMSoRBCW5kFuxctUpd+AnW+4bhxEBioPleU6RbChCQRCiG05+AAr72mTrd4663sxYDPnYPOneGll9RbqUKYgCRCIYTl8PCAL75Qnxe2aJHd/t//qosBv/ce3L6tWXjCNkkiFEJYnuBg2L5drUJTpYralpYGn32mTrdYtkwWAxZGI4lQCGGZdDro0QNiYtQ6pa6uantCAgwYAE2awKFDmoYobIMkQiGEZStWDCZOVBNi9+7Z7QcOqHVNBw1Sk6MQhSSJUAhhHapVg3Xr1DUOAwKy25csUW+XzpoliwGLQpFEKISwLi+8oA6m+fxzKFVKbbt9G959F0JC4LfftI1PWB1JhEII6+PkBCNGwOnTMGSI+jwRIDZWLe7dpYs69UIIA0giFEJYL09PWLxYHTTTuHF2+8aN6u3T8eMhOVm7+IRVkEQohLB+9erB3r3w/ffg7a223b8Pn3yirnyxapVUpxG5kkQohLANOh307atWoBkzRr19CmoB7/BwaN4c/vhD2xiFRZJEKISwLe7u8Omn6nJPoaHZ7bt3Q9268Oab6tJPQjwkiVAIYZt8feGnn9QvX1+1LSMDFi5UXy9YoFarEXZPEqEQwraFhkJUFEyfDiVKqG03bsCwYeqzxYgIbeMTmpNEKISwfS4uMHq0+vywX7/s9shItbh3797w55+ahSe0JYlQCGE/KlaE775TR5jWrZvdvnq1uhjwlCmyGLAdkkQohLA/TZrAwYPw1VdQrpzadu8efPihOv9wwwaZbmFHJBEKIeyTgwMMHqwuBvz22+prgPPnoWtXtUJNTIy2MQqzMGkijI+PZ+zYsbRs2RJ3d3f0ej0RBXgwPXHiRPR6/RNfbm5uJoxaCGFXypSBOXPUOYatWmW3b92q1i7917/g1i3t4hMm52jKg8fGxjJjxgz8/PwICQlh//796DJrAhbAokWLKJE52gtwyPzLTQghjCUwUF3ZYv16NfldvKhOr5g9G1asgGnT1HUQ9XIjzdaYNBE+++yzXL9+ndKlS7Nu3Tr2799fqOP06NEDDw8PI0cnhBD/Q6eDbt2gfXuYOVNNfikpkJgIr74KixbB3LnQqJHWkQojMumfNiVKlKB06dJFPk5GRga3b99GkYfXQghzcHODjz5SV7Po2TO7PbO494ABEB+vWXjCuKyij1+9enVKly5NyZIl6devH4mJiVqHJISwB1Wrwpo1sH07BAVlty9bBn5+XJ04kfe+/JJQoBMQ+vHHvDdwIFevXtUqYlEIJr01WlQeHh4MHz6cxo0b4+Liwq5du5g/fz4HDx7k8OHDuLu7ax2iEMIetGwJx46pt0Y//BBu3iQxKYnekyYxFZgB6ICM+HgOLl1K2O7drN6/H09PT40DF4YwOBEqisL9+/cN2tfV1bXQAT1qxIgRj73u2rUrDRo0oE+fPixYsIAxY8YY5TxCCJEvR0d46y21Cs348cz88kumAo8+LdQ/fP3J2bPMGD2amUuWaBOrKBCDb41GRERQrFgxg77i4uJMFnB4eDgVKlTg999/N9k5hBAiV+XKwaJFRFevTsNcdmkIRB88aM6oRBEY3CP09/dn6dKlBu1boUKFwsZjkMqVK3P9+vU89xk5ciSlSpV6rC08PJzw8HBThiaEsBMOjo7kNhlMDzjIyhZmtXLlSlauXPlY2y0D538anAi9vLzo379/wSIzAUVRuHDhAvXq1ctzv9mzZ1P30VqCQghhROmOjiiQYzLMeLhdmE9OHZ2jR4/mmyvAgkaNXrp0idjY2Mfachp5tXDhQq5du0a7du3MFZoQQjwhoEEDDuSy7cDD7cI6mPxPlilTpgBw8uRJAL777jt27doFwPjx47P269+/P7t27SIjIyOrrWrVqvTu3ZugoCBcXV3Zs2cPq1evpk6dOrz22mumDl0IIXI1esYMwnbv5pOzZ2mI2qvIQE2C42rUYPWMGdoGKAynmJhOp1P0ev1j/2b+96NatGjxRNuQIUOUwMBApWTJkoqzs7Pi5+envP/++8qdO3dyPd+RI0cUQDly5Eiecf3nP/8p/DclCkyut/nJNTe9xMREZdSAAcpLAQFK3QoVlJcCApRRAwYoiYmJWodm8wz5+TY0H5i8R/hoDy8vO3bseKJt8eLFxg4ny8qVK2XgjBnJ9TY/ueam5+npmTVFolOnTmzatEnjiOyHMX++LeYZoRBCCKEFSYRCCCHsmiRCIYQQds3mJrrcu3cPgJh8Vpa+desWR48eNUdIArneWpBrbl5yvc3LkOudmQcy80JudIpiW2sbrVixgr59+2odhhBCCAuxfPly+vTpk+t2m0uE165dY8uWLVSrVg03NzetwxFCCKGRe/fuceHCBdq2bUu5cuVy3c/mEqEQQghREDJYRgghhF2TRCiEEMKuSSIUQghh1yQRCiGEsGuSCIUQQtg1SYTA77//zqBBg/Dz86N48eLUqFGDIUOGEB8fr3VoNik+Pp6xY8fSsmVL3N3d0ev1REREaB2WTbh//z5jxoyhYsWKFCtWjEaNGrFt2zatw7JJycnJTJgwgXbt2uHh4YFer2fZsmVah2WzDh06xFtvvUVgYCAlSpSgatWqhIWFcfr06SIfWxIhMGbMGHbt2kX37t354osv6N27N2vWrKFOnTokJCRoHZ7NiY2NZcaMGVy5coWQkBAAdLqc1vkWBTVgwABmz55Nv379mDt3Lg4ODrz00kvs3btX69BsztWrV5k8eTKnTp3imWeeAeTn2JSmT5/O+vXrad26NXPnzmXo0KHs2rWLunXrZq13W2hFXxXK+u3evfuJtl27dik6nU4ZP368BhHZtqSkJOXGjRuKoijK2rVrFZ1Op0RERGgclfU7cOCAotPplH//+99ZbSkpKYqPj4/SpEkTDSOzTffv31cSEhIURVGUw4cPKzqdTlm2bJnGUdmuffv2KampqY+1nT59WnF1dVX69u1bpGNLjxB47rnnnmh7/vnn8fDwIDY2VoOIbFuJEiUoXbq01mHYnHXr1uHo6MjQoUOz2lxcXHj11VfZv38/f//9t4bR2R5nZ2fKly8PgCJ1SUyucePGODo+Xh7bx8eHgICAIn9OSyLMxZ07d0hKSsqzLI8QluTYsWP4+flRokSJx9rr168PwPHjx7UISwiTURSFhISEIn9OSyLMxZw5c0hNTSUsLEzrUIQwyJUrV/D29n6iPbPt8uXL5g5JCJNasWIFly9fLvLntM0tw6QoCvfv3zdoX1dX1xzbd+3axaRJkwgLC6NFixZGjM72GON6C+O4d+8eLi4uT7RnXvf8lqIRwprExsYybNgwmjRpwiuvvFKkY9lcjzAiIoJixYoZ9BUXF/fE+2NjY+natSshISF8/fXXGnwH1qWo11sYj5ubW45/lKSkpGRtF8IWxMfHExoaSpkyZVi3bl2RR+vaXI/Q39+fpUuXGrRvhQoVHnv9559/0qZNG8qUKcMvv/xC8eLFTRChbSnK9RbG5e3tnePtzytXrgBQsWJFc4ckhNHdunWL9u3bc/v2bXbv3m2UzxWbS4ReXl7079+/wO/7559/aNOmDampqezYsQMvLy8TRGd7Cnu9hfHVqVOHnTt3kpSUhLu7e1b7gQMHALLmuglhrVJSUujYsSNnzpxh27Zt1KpVyyjHtblbo4WRnJzMSy+9xJUrV/jll1+oUaOG1iEJUWA9evQgPT2dxYsXZ7Xdv3+fJUuW0KhRIypVqqRhdEIUTXp6OmFhYRw4cIC1a9fSsGFDox3b5nqEhdGnTx8OHTrEoEGDOHny5GNVCtzd3encubOG0dmmKVOmAGRd6++++45du3YBMH78eM3ismYNGjSgZ8+evP/++yQmJlKjRg2WLVvGpUuXWLJkidbh2aR58+Zx8+bNrFvSmzZt4tKlSwCMGDGCkiVLahmeTXn33XfZvHkzHTt25Nq1ayxfvvyx7X379i30sWWFeuDpp5/m0qVLOU6KrVatGufOndMgKtum1+vR6XQoipL1L6glqtLT0zWOznrdv3+fDz/8kOXLl3Pjxg1q167N5MmTad26tdah2aSnn36aixcvAtnl1TJ/ps+fP0+VKlW0DM+mtGzZkl27duX4OV3Uzw1JhEIIIeyaPCMUQghh1yQRCiGEsGuSCIUQQtg1SYRCCCHsmiRCIYQQdk0SoRBCCLsmiVAIIYRdk0QohBDCrkkiFEIIYdckEQohhLBrkgiFEELYNUmEQggh7Nr/A3d7c6QobPPhAAAAAElFTkSuQmCC", + "text/plain": [ + "PyPlot.Figure(PyObject )" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "(-1.6,2.1)" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using PyPlot\n", + "fig = figure(figsize=(5.0, 5.0))\n", + "plot(S[1, [1, 2, 3, 1]]', S[2, [1, 2, 3, 1]]', \"-yo\", lw=2.0)\n", + "plot(M[1, [1, 2, 3, 1]]', M[2, [1, 2, 3, 1]]', \"-ro\", lw=2.0)\n", + "#plot(P[1, [1, 2, 3, 4, 5, 6, 1]]', P[2, [1, 2, 3, 4, 5, 6, 1]]', \"-k\", lw=2.5)\n", + "npts = size(P, 2)\n", + "# loop over cells\n", + "ips = JuliaFEM.Core.get_integration_points(JuliaFEM.Core.Tri3, Val{5})\n", + "#info(\"ips = $ips\")\n", + "SM = zeros(3, 3)\n", + "SS = zeros(3, 3)\n", + "for i=1:npts\n", + " xvec = [C[1], P[1, i], P[1, mod(i, npts)+1]]\n", + " yvec = [C[2], P[2, i], P[2, mod(i, npts)+1]]\n", + " plot(xvec, yvec, \"--ko\", lw=1.0)\n", + " X = hcat(xvec, yvec)'\n", + " geom = JuliaFEM.Core.Field(Vector{Float64}[X[:,i] for i=1:size(X,2)])\n", + " #info(\"geom = $geom\")\n", + " for ip in ips\n", + " dN = JuliaFEM.Core.get_dbasis(JuliaFEM.Core.Tri3, ip.xi)\n", + " J = sum([kron(dN[:,i], geom[i]') for i=1:length(geom)])\n", + " w = ip.weight*det(J)\n", + " x = vec(JuliaFEM.Core.get_basis(JuliaFEM.Core.Tri3, ip.xi)*geom)\n", + " plot(x[1], x[2], \".k\")\n", + " theta1 = JuliaFEM.Core.project_point_from_plane_to_surface(x, x0, Q, sel, time)\n", + " N1 = sel(theta1[2:3], time)\n", + " theta2 = JuliaFEM.Core.project_point_from_plane_to_surface(x, x0, Q, mel, time)\n", + " N2 = mel(theta2[2:3], time)\n", + " SS += w*N1'*N1\n", + " SM += w*N1'*N2\n", + " end\n", + "end\n", + "plot(C[1], C[2], \"ko\")\n", + "xlim(-2.1, 2.1)\n", + "ylim(-1.6, 2.1)\n", + "#axis(\"off\")" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: SS = \n", + "[1.024777091906741 0.5824974279835513 0.4791452331961691\n", + " 0.5824974279835513 0.8193479938271759 0.47546939300412483\n", + " 0.4791452331961691 0.47546939300412483 0.4983174725651672]\n", + "INFO: SM = \n", + "[0.8084490740740896 0.7707818930041307 0.5071887860082416\n", + " 0.43219521604939215 0.7584104938271751 0.686709104938285\n", + " 0.4931520061728493 0.3671039094650286 0.5926761831275834]\n" + ] + } + ], + "source": [ + "info(\"SS = \\n$SS\")\n", + "info(\"SM = \\n$SM\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 0.4.1", + "language": "julia", + "name": "julia-0.4" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "0.4.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/src/integrate.jl b/src/integrate.jl index 102e00b..8f65d8a 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -54,24 +54,46 @@ function get_integration_points(::Type{Seg3}) return get_integration_points(Seg3, Val{3}) end -### 2d elements +### 2d triangular elements -function get_integration_points(::Type{Tri3}) +# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF + +typealias TriangularElements Union{Type{Tri3}, Type{Tri6}} + +function get_integration_points(::TriangularElements, ::Type{Val{1}}) # http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html [ IntegrationPoint([1.0/3.0, 1.0/3.0], 0.5) ] end -function get_integration_points(::Type{Tri6}) +function get_integration_points(::TriangularElements, ::Type{Val{2}}) # http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html [ - IntegrationPoint([2.0/3.0, 1.0/6.0], 1.0/6.0) - IntegrationPoint([1.0/6.0, 2.0/3.0], 1.0/6.0) + IntegrationPoint([2.0/3.0, 1.0/6.0], 1.0/6.0), + IntegrationPoint([1.0/6.0, 2.0/3.0], 1.0/6.0), IntegrationPoint([1.0/6.0, 1.0/6.0], 1.0/6.0) ] end +function get_integration_points(::TriangularElements, ::Type{Val{5}}) + # http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF + [ + IntegrationPoint([0.33333333333333, 0.33333333333333], 0.22500000000000), + IntegrationPoint([0.47014206410511, 0.47014206410511], 0.13239415278851), + IntegrationPoint([0.47014206410511, 0.05971587178977], 0.13239415278851), + IntegrationPoint([0.05971587178977, 0.47014206410511], 0.13239415278851), + IntegrationPoint([0.10128650732346, 0.10128650732346], 0.12593918054483), + IntegrationPoint([0.10128650732346, 0.79742698535309], 0.12593918054483), + IntegrationPoint([0.79742698535309, 0.10128650732346], 0.12593918054483) + ] +end + +function get_integration_points(::Type{Tri3}) + return get_integration_points(Tri3, Val{1}) +end + + function get_integration_points(::Type{Quad4}) [ IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0), @@ -102,4 +124,3 @@ function get_integration_points(::Type{Tet10}) IntegrationPoint([b, b, b], w) ] end - diff --git a/src/mortar.jl b/src/mortar.jl index 7820087..6411c79 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -124,6 +124,434 @@ function project_from_master_to_slave{S,M}(slave::Element{S}, master::Element{M} error("find projection from master to slave: did not converge") end +### Mortar projection calculation for 3d cases + +""" +Construct auxiliary plane for surface. + +Parameters +---------- +x::Array{Float64, 2} + Node coordinates +ximp::Array{Float64, 1} + Element mid-point in dimensionless mother element coordinates ξ +normals::Array{Float64, 2} + Normal directions in nodes + +Returns +------- +x0, Q + x0::Array{Float64, 1} - origo of auxiliary plane + Q::Array{Float64, 2} - orthogonal basis, first vector is normal direction + and two rest vectors create orthonormal right-handed basis. + +Examples +-------- +Calculate auxiliary plane given nodal coordinates, midpoint of mother element, +node normals and suitable function space: + +julia> xquad = [ +... -2.5 2.5 2.0 -2.0 +... -2.0 -2.0 2.3 2.0 +... 1.0 0.7 0.0 1.0] +julia> m_midpoint = [0.0, 0.0] +julia> normals = [ +... 0.05989060 0.0590504 0.225612 0.2445800 +... -0.00748633 0.1670810 0.182034 -0.0305725 +... 0.99817700 0.9841730 0.957059 0.9691470] +julia> basis(xi) = [ +... (1-xi[1])(1-xi[2])/4 +... (1+xi[1])(1-xi[2])/4 +... (1+xi[1])(1+xi[2])/4 +... (1-xi[1])(1+xi[2])/4]' +julia> x0, Q = create_auxiliary_plane(xquad, mmidpoint, normals, basis) +julia> x0 +3-element Array{Float64,1}: + 0.0 + 0.075 + 0.675 +julia> Q +3x3 Array{Float64,2}: + 0.148586 0.988899 0.0 + 0.0784519 -0.0117877 0.996848 + 0.985783 -0.148118 -0.0793325 + +Notes +----- +- Midpoint in mother element typically (0, 0) for quadrangles and (1/3, 1/3) + for triangles. +- Uses Gram-Schmidt process to find orthogonal basis +- [1](http://www.math.umn.edu/~olver/aims_/qr.pdf) +- [2](http://www.ecs.umass.edu/ece/ece313/Online_help/gram.pdf) +- [3](http://www.terathon.com/code/tangent.html) + +""" +# function create_auxiliary_plane(x, ximp, normals, basis) +function create_auxiliary_plane(element::Element{Tri3}, time::Real) + proj(u, v) = dot(v, u) / dot(u, u) * u + xi = [1.0/3.0, 1.0/3.0] + x0 = element("geometry", xi, time) + n = element("nodal ntsys", xi, time)[:, 1] + n /= norm(n) + # gram-schmidt + u1 = n + j = indmax(abs(u1)) + v2 = zeros(3) + v2[mod(j,3)+1] = 1.0 + u2 = v2 - proj(u1, v2) + u3 = cross(u1, u2) + t1 = u2/norm(u2) + t2 = u3/norm(u3) + new_basis = [n t1 t2] + return x0, new_basis +end + +""" +Project point q onto a plane given by a point p and normal n. + +Parameters +---------- +q::Array{Float64, 2} + point to project (row vector) +x0::Array{Float64, 2} + origo of plane +n::Array{Float64, 2} + normal vector of plane + +Returns +------- +y::Array{Float64, 2} + projected point + +Examples +-------- +julia> p = [-0.5 -1.0 4.0]' +julia> x0 = [0.0 0.075 0.675]' +julia> n = [0.1485860 0.0784519 0.9857830]' +julia> project_node_to_auxiliary_plane(p, x0, n) +3-element Array{Float64,1}: + 0.963455 + -1.2447 + 0.925247 + +Notes +----- +[1](http://stackoverflow.com/questions/8942950/how-do-i-find-the-orthogonal-projection-of-a-point-onto-a-plane) + +""" +function project_point_to_auxiliary_plane(p::Vector, x0::Vector, Q::Matrix) + n = Q[:,1] + ph = p - dot(p-x0, n)*n + qproj = Q'*(ph-x0) + @assert isapprox(qproj[1], 0.0) + return qproj[2:3] +end + +""" +Find edge intersections of two planar arbitrary shape polygons. + +Parameters +---------- +S::Array{Float64,2} +M::Array{Float64,2} + +Matrices with size (2, n) where n is number of vertices of each polygon. + +Returns +------- +P::Array{Float64,2} + Intersection points of polygons +n::Array{Float64,2} + Neighbour info matrix with size (ns, mn). This keeps information which + edges of polygons are intersecting. See further explanation in example + below. + +Examples +-------- +Find intersection points of two triangles: + +julia> S = [0 0; 3 0; 0 3]' +julia> M = [-1 1; 2 -1/2; 1 3/2]' +julia> P, n = get_edge_intersections(S, M) +julia> P +2x4 Array{Float64,2}: + 1.0 1.75 0.0 0.0 + 0.0 0.0 0.5 1.25 +julia> n +3x3 Array{Int64,2}: + 1 1 0 + 0 0 0 + 1 0 1) + +So intersection points are: (1.00, 0.00), (1.75, 0.00), (0.00, 0.50), (0.00, 1.25). +"Neighbour matrix" can be interpreted as following: + + 1 1 0 <--> First edge of S intersects edges 1 and 2 of M + 0 0 0 <--> Second edge of S doesn't intersect at all + 1 0 1 <--> Third edge of S intersects with edges 1 and 3 of M + +""" +function get_edge_intersections(S::Matrix, M::Matrix) + ns = size(S, 2) + nm = size(M, 2) + P = zeros(2, 0) + n = zeros(Int64, ns, nm) + k = 0 + for i=1:ns + for j=1:nm + b = M[:,j]-S[:,i] + A = [S[:,mod(i,ns)+1]-S[:,i] -M[:,mod(j,nm)+1]+M[:,j]] + if rank(A) == 2 + r = A\b + if (r[1]>=0) & (r[1]<=1) & (r[2]>=0) & (r[2]<=1) # intersection found + k += 1 + f = S[:,i]+r[1]*(S[:,mod(i,ns)+1] - S[:,i]) + f = f'' + P = hcat(P, f) + n[i, j] = 1 + end + end + end + end + return P, n +end + + +""" +Find any points laying inside or border of triangle. + +Parameters +---------- +Y::Array{Float64, 2} + Triangle coordinates in 2×3 matrix +X::Array{Float64, 2} + List of points to test in 2×n matrix + +Returns +------- +P::Array{Float64, 2} + List of points in triangle in 2×m matrix, where m is number of points inside triangle + +Examples +-------- +julia> S = [0.0 0.0; 3.0 0.0; 0.0 3.0]' # triangle corner points +julia> pts = [-1.0 1.0; 2.0 -0.5; 1.0 1.5; 0.5 1.5]' # points to tests +julia> points_in_triangle(S, pts) +2x2 Array{Float64,2}: + 1.0 0.5 + 1.5 1.5 +""" +function get_points_inside_triangle(Y::Matrix, X::Matrix) + @assert size(Y, 2) == 3 # "Point in TRIANGLE..." + P = zeros(2, 0) + v0 = Y[:,2] - Y[:,1] + v1 = Y[:,3] - Y[:,1] # find interior points of X in Y + d00 = (v0'*v0)[1] + d01 = (v0'*v1)[1] + d11 = (v1'*v1)[1] # using baricentric coordinates + id = 1/(d00*d11 - d01*d01) + for i=1:size(X, 2) + v2 = X[:,i] - Y[:,1] + d02 = (v0'*v2)[1] + d12 = (v1'*v2)[1] + u = (d11*d02-d01*d12)*id + v = (d00*d12-d01*d02)*id + if (u>=0) & (v>=0) & (u+v<=1) # also include nodes on the boundary + P = hcat(P, X[:,i]'') + end + end + return P +end + + +""" +Make polygon clipping of shapes S and M. + +Parameters +---------- +S::Array{Float64, 2} +M::Array{Float64, 2} + Shapes to clip. Needs to be triangles at the moment. + +Returns +------- +Array{Float64, 2}, Array{Float64, 2} +- Polygon vertices in 2×n matrix, sorted in counter-clockwise order. +- 3×3 "neighbouring" matrix, see example. + +Examples +-------- +julia> S = [0 0; 3 0; 0 3]' +julia> M = [-1 1; 2 -1/2; 2 2]' +julia> P, n = clip_polygon(S, M) +julia> P +2x6 Array{Float64,2}: + 0.0 1.0 2.0 2.0 1.25 0.0 + 0.5 0.0 0.0 1.0 1.75 1.33333, +julia> n +3x3 Array{Int64,2}: + 1 0 1 <- first edge of M ([-1 1; 2 -1/2]') intersects with edges 1 and 3 of S ([0 0; 3 0]' and [0 3; 0 0]') + 1 1 0 <- second edge of M ([2 -1/2; 2 2]') intersects with edges 1 and 2 of S + 0 1 1 <- third edge of M ([2 2; -1 1]') intersects with edgse 2 and 3 of S + +""" +function clip_polygon(S::Matrix, M::Matrix) + P1, neighbours = get_edge_intersections(M, S) + P2 = get_points_inside_triangle(M, S) + P3 = get_points_inside_triangle(S, M) + P = hcat(P1, P2, P3) + meanval = mean(P, 2) + tmp = P .- meanval + angles = atan2(tmp[2,:], tmp[1,:]) + angles = reshape(angles, length(angles)) + order = sortperm(angles) + P = copy(unique(P[:, order], 2)) + return P, neighbours +end + + +""" +Calculate polygon geometric center point + +Parameters +---------- +P::Array{Float64, 2} + Polygon vertices in 2×n matrix + +Returns +------- +Array{Float63, 2} + Center point + +Examples +-------- +julia> P +2x6 Array{Float64,2}: + 0.0 1.0 2.0 2.0 1.25 0.0 + 0.5 0.0 0.0 1.0 1.75 1.33333, +julia> C = get_polygon_cp(P) +2x1 Array{Float64,2}: + 1.039740 + 0.804701 + +""" +function calculate_polygon_centerpoint(P::Matrix) + n = size(P, 2) + A = 0.0 + for i=1:n + A += 1/2*(P[1,i]*P[2,mod(i,n)+1] - P[1,mod(i,n)+1]*P[2,i]) + end + Cx = 0.0 + Cy = 0.0 + for i=1:n + inext = mod(i, n)+1 + Cx += 1/(6*A)*(P[1,i] + P[1,inext])*(P[1,i]*P[2,inext] - P[1,inext]*P[2,i]) + Cy += 1/(6*A)*(P[2,i] + P[2,inext])*(P[1,i]*P[2,inext] - P[1,inext]*P[2,i]) + end + return Float64[Cx, Cy] +end + +""" +Project point from auxiliary plane to parametric surface given by (ξ₁, ξ₂) + +Parameters +---------- +p::Array{Float64,1} + point in auxiliary plane, in (n,t1,t2) coordinate system +x0::Array{Float64,1} + origo of auxiliary plane cs +Q::Array{Float64,2} + basis of auxiliary plane cs +x::Array{Float64,2} + surface node coords +basis::Array{Float64,2} + surface basis functions +dbasis::Array{Float64,2} + partial derivatives of surface basis functions + +Returns +------- +Array{Float64,2} + solution vector (d, ξ₁, ξ₂) where d is distance to surface + +Examples +-------- +Define surface with node points, basis + dbasis + +julia> xquad = [ +... -2.5 -2.0 1.0 +... 2.5 -2.0 0.7 +... 2.0 2.3 0.0 +... -2.0 2.0 1.0]' +julia> basis(xi) = [ +... (1-xi[1])(1-xi[2])/4 +... (1+xi[1])(1-xi[2])/4 +... (1+xi[1])(1+xi[2])/4 +... (1-xi[1])(1+xi[2])/4] +julia> dbasis(xi) = [ +... -(1-xi[2])/4 -(1-xi[1])/4 +... (1-xi[2])/4 -(1+xi[1])/4 +... (1+xi[2])/4 (1+xi[1])/4 +... -(1+xi[2])/4 (1-xi[1])/4] + +We aim to find point p, which we first project to auxiliary plane defined as following +julia> p = [-2.5 -2.0 1.0]' +julia> x0 = [0.0 0.075 0.675]' +julia> Q = [ +... 0.1485860 0.9888990 0.0000000 +... 0.0784519 -0.0117877 0.9968480 +... 0.9857830 -0.1481180 -0.0793325] + +Our projected point is therefore +julia> n = Q[:,1] # first component is normal direction +julia> ph = project_node_to_auxiliary_plane(p, x0, n) +julia> ph = Q'(ph-x0) +julia> ph +3x1 Array{Float64,2}: + 1.33264e-7 + -2.49593 + -2.09424 + +Our point ph is now in auxiliary plane in n,t1,t2 coordinate system. Next we +project it back to surface defined by xquad*basis + +julia> theta = project_point_from_plane_to_surface(ph, x0, Q, xquad, basis, dbasis) +julia> theta +3x1 Array{Float64,2}: + -0.213874 + -0.999999 + -1.0 + +We see that our ξ₁ = ξ₂ = -1 so we found first point of xquad +[-2.5 -2.0 1.0]' correctly. + +julia> xquad*basis(theta[2:3]) +3-element Array{Float64,1}: + -2.5 + -2.0 + 1.0 + +""" +function project_point_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix, element::Element{E}, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-9) + basis(xi) = get_basis(E, xi) + dbasis(xi) = get_dbasis(E, xi) + x = element("geometry", time) + ph = Q*[0; p] + x0 + theta = Float64[0.0, 0.0, 0.0] + n = Q[:,1] + for i=1:max_iterations + b = ph + theta[1]*n - basis(theta[2:3])*x + J = [n -dbasis(theta[2:3])*x] + dtheta = J \ -b + theta += dtheta + if norm(dtheta) < iter_tol + return theta + end + end + error("project_point_to_auxiliary_plane: did not converge in $max_iterations iterations!") +end + + ### Mortar problem """ @@ -142,7 +570,7 @@ end # Mortar assembly function assemble!(assembly::Assembly, problem::BoundaryProblem{MortarProblem}, slave_element::Element, time::Number) - + # get dimension and name of PARENT field field_dim = problem.parent_field_dim field_name = problem.parent_field_name @@ -167,7 +595,7 @@ function assemble!(assembly::Assembly, problem::BoundaryProblem{MortarProblem}, # projected integration point xi_projected = project_from_slave_to_master(slave_element, master_element, xi_gauss) - # add contribution to left hand side + # add contribution to left hand side N1 = slave_element(xi_gauss, time) N2 = master_element(xi_projected, time) S = w*N1'*N1 @@ -184,4 +612,3 @@ function assemble!(assembly::Assembly, problem::BoundaryProblem{MortarProblem}, end end end - diff --git a/test/test_mortar.jl b/test/test_mortar.jl index 4bb88ad..3e69004 100644 --- a/test/test_mortar.jl +++ b/test/test_mortar.jl @@ -5,10 +5,18 @@ module MortarTests using JuliaFEM.Test -using JuliaFEM.Core: Element, Seg2, Quad4, MortarProblem, Assembly, assemble! +using JuliaFEM.Core: Element, Seg2, Quad4, Tri3, MortarProblem, Assembly, assemble! using JuliaFEM.Core: PlaneStressElasticityProblem, DirichletProblem, DirectSolver + +# 2d stuff using JuliaFEM.Core: project_from_slave_to_master, project_from_master_to_slave +# 3d stuff +using JuliaFEM.Core: create_auxiliary_plane, project_point_to_auxiliary_plane, + get_edge_intersections, get_points_inside_triangle, + clip_polygon, calculate_polygon_centerpoint, + project_point_from_plane_to_surface + function get_test_2d_model() # this is hand calculated and given as an example in my thesis N = Vector[ @@ -118,7 +126,7 @@ function test_create_flat_2d_assembly() info("size of B = $(size(B))") info("B matrix in first slave element = \n$(B[10:11,:])") info("B matrix expected = \n$(B_expected[10:11,:])") - @test isapprox(B, B_expected) + @test isapprox(B, B_expected) fill!(B_expected, 0.0) empty!(assembly) @@ -184,7 +192,7 @@ function test_2d_mortar_multiple_bodies_multiple_dirichlet_bc() dy1 = Seg2([1, 2]) dy1["geometry"] = Vector[N[1], N[2]] dy1["displacement 2"] = 0.0 - + boundary2 = DirichletProblem("displacement", 2) push!(boundary2, dy1) @@ -282,7 +290,7 @@ function test_2d_mortar_three_bodies_shared_nodes() dy1 = Seg2([1, 2]) dy1["geometry"] = Vector[N[1], N[2]] dy1["displacement 2"] = 0.0 - + bc2 = DirichletProblem("displacement", 2) push!(bc2, dy1) @@ -343,4 +351,113 @@ function test_2d_mortar_three_bodies_shared_nodes() end #test_2d_mortar_three_bodies_shared_nodes() +function test_auxiliary_plane_transforms() + nodes = Vector{Float64}[ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0]] + e1 = Tri3([1, 2, 3]) + # local coordinate system N, T1, T2 in node + R = [0.0 1.0 0.0 + 0.0 0.0 1.0 + 1.0 0.0 0.0] + e1["geometry"] = Vector{Float64}[nodes[1], nodes[2], nodes[3]] + e1["nodal ntsys"] = Matrix{Float64}[R, R, R] + time::Real = 0.0 + x0, Q = create_auxiliary_plane(e1, time) + info("x0 = $x0") + info("Q = $Q") + @test isapprox(x0, [1.0/3.0, 1.0/3.0, 0.0]) + @test isapprox(Q, R) + p1 = Float64[1.0/3.0+0.1, 1.0/3.0+0.1, 1.0] + p2 = project_point_to_auxiliary_plane(p1, x0, Q) + info("point in auxiliary plane p2 = $p2") + @test isapprox(p2, [0.1, 0.1]) + theta = project_point_from_plane_to_surface(p2, x0, Q, e1, time) + info("theta = $theta") + @test isapprox(theta[1], 0.0) + X = e1("geometry", theta[2:3], time) + info("projected point = $X") + @test isapprox(X, Float64[1.0/3.0+0.1, 1.0/3.0+0.1, 0.0]) +end +test_auxiliary_plane_transforms() + + +function test_get_edge_intersections() + # first case, two triangles + S = [ 0.0 0.0; 3.0 0.0; 0.0 3.0]' + M = [-1.0 1.0; 2.0 -0.5; 1.0 1.5]' + P, n = get_edge_intersections(S, M) + P_expected = [ + 1.00 1.75 0.00 0.00 + 0.00 0.00 0.50 1.25] + n_expected = [ + 1 1 0 + 0 0 0 + 1 0 1] + @test isapprox(P, P_expected) + @test isapprox(n, n_expected) + + # slave 4 vertices non-convex, master triangle + S = [ 0.0 0.0; 2.5 0.0; 1.0 1.0; 0.0 2.0]' + M = [-1.0 1.0; 2.0 -0.5; 1.0 1.5]' + P, n = get_edge_intersections(S, M) + P_expected = [ + 1.0 1.75 1.375 0.60 0.00 0.00 + 0.0 0.00 0.750 1.40 0.50 1.25] + n_expected = [ + 1 1 0 + 0 1 0 + 0 0 1 + 1 0 1] + @test isapprox(P, P_expected) + @test isapprox(n, n_expected) + + # slave 3 triangle, master 4 vertices + S = [ 0.0 0.0; 3.0 0.0; 0.0 3.0]' + M = [-1.0 1.0; 2.0 -0.5; 1.0 1.5; -1.0 2.0]' + P, n = get_edge_intersections(S, M) + P_expected = [ + 1.00 1.75 0.00 0.00 + 0.00 0.00 0.50 1.75] + n_expected = [ + 1 1 0 0 + 0 0 0 0 + 1 0 1 0] + @test isapprox(P, P_expected) + @test isapprox(n, n_expected) +end +#test_get_edge_intersections() + + +function test_get_points_inside_triangle() + S = [0.0 0.0; 3.0 0.0; 0.0 3.0]' + pts = [-1.0 1.0; 2.0 -0.5; 1.0 1.5; 0.5 1.5]' + P = get_points_inside_triangle(S, pts) + @test isapprox(P, [1.0 1.5; 0.5 1.5]') +end +#test_get_points_inside_triangle() + + +function test_polygon_clipping() + S = [0 0; 3 0; 0 3]' + M = [-1 1; 2 -1/2; 2 2]' + P, n = clip_polygon(S, M) + @test isapprox(P, [0.0 0.5; 1.0 0.0; 2.0 0.0; 2.0 1.0; 1.25 1.75; 0.0 4/3]') + @test isapprox(n, [1 0 1; 1 1 0; 0 1 1]) +end +#test_polygon_clipping() + + +function test_calculate_polygon_centerpoint() + P = [ + 0.0 1.0 2.0 2.0 1.25 0.0 + 0.5 0.0 0.0 1.0 1.75 1.33333] + C = calculate_polygon_centerpoint(P) + info("Polygon centerpoint: $C") + @test isapprox(C, [1.0397440690338993, 0.8047003412233396]) +end +#test_calculate_polygon_centerpoint() + + end