{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Element volume locking\n",
    "\n",
    "***\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## General formulation of the problem\n",
    "\n",
    "---\n",
    "\n",
    "Let us remind the relationships between $(\\lambda,\\mu)$ and $(E, \\nu)$\n",
    "\\begin{equation}\n",
    "\\lambda=\\frac{E\\nu}{(1+\\nu)(1-2\\nu)},\\,\\,\\,\\,\\mu=\\frac{E}{2(1+\\nu)}.\n",
    "\\end{equation}\n",
    "\n",
    "The basic formulation introduced in the 2D problem is based on stiffness,\n",
    "but for $\\nu\\rightarrow 0.5$ the Young modulus $E\\rightarrow\\infty$\n",
    "and hence the calculation collapses near $\\nu=0.5$.\n",
    "This singularity introduces a so called volume locking,\n",
    "which can be seen from the following simple numerical experiment.\n",
    "Let us consider a rectangular domain, fixed on the bottom and loaded on the top edge.\n",
    "Moreover let us suppose an almost critical value of Poisson ratio $\\nu=0.4999$.\n",
    "\n",
    "The problem is implemented as in *2D_Elasticity* example.\n",
    "The domain response should be symmetrical; however, the result is..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAD4CAYAAADLqNJwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2da4xkR3XH/+fe7vGbkMQTYmE2i8QjQUQBMjFS7JCABLHBIhFSpCBBEgVpv5CIpwgOKIEofEgUMChCilZgAgoBoRiUiNgYKxiQJV67xLxsYsCsiR3bu9i78+jpnn7ckw/TPdPTfR+n+lZ1Vd0+P6nl3dna6tpx39/UrfuvU8TMUBRFWZTE9wAURYkblYiiKLVQiSiKUguViKIotVCJKIpSi5aLTq+88ko+fvy4i64VRfHE6dOnf8rM67NfdyKR48eP49SpUy66VhTFE0T0YN7X9XZGUZRaqEQURamF6HaGiM4A2AYwAjBk5g2Xg1IUJR5M1kRezMw/dTYSRVGiRG9nFEWphVQiDODzRHSaiE7kNSCiE0R0iohOnTt3zt4IFUUJGqlErmPmFwC4AcDriehFsw2Y+SQzbzDzxvr63KNkRVEaikgizPzw+L9nAXwGwDUuB6UoSjxULqwS0WUAEmbeHv/6ZQD+xvnIFEVZmOuf91eVbT53j53LWPJ05ikAPkNEk/b/ysyfs/LuiqKIeOlv/m1lm3S3L+9wMKoxmqNUSoSZHwDwa9beUVEUAMCLX/Z3lW1aOwP7b2xRIICjvTOKsopc96p/ELVr79i9iAEg2RPKxrJAAJWIouTywj96n7hteydzMoa0NxS18ykQQCWirAC/+tabAQCtXfnfWXM0lrQnE04sAgFUIkpk/PK7bj7y+/aOp4GMaXUZnBJoVH1qQhMFAqhEFI884+9vzv26yYzBBa2e7BiVVld+3EpTBQKoRBQLHP/geyvbtDp+t2mle7J2E4FkLSApuZ5VIIeoRJQjHP9o9WPHA3bsf3ykFzsApD27fcY+AzlyW7UkgQAqkUbzgtveie3OxeL2g922k3GkeyRsZ9DnAgLhFKCCa8u2QKTy2G+7nBkIOTrtUiUSCb/7pTce/Ppc53I3b5IwkMkuePRSUbNQBFJGyALhFoGGfCAQGpb3TdIZyMjeY2mViAf++Ouvm/vao90rPIzkkEHP4KPgSSB58ihau4hFIMmgun0ijLP7EAigErHCO779qtyvP9T72SWP5Ci7vTWkaYbRqHxRM1aB1O3ThUBSgRQAeZydhLc6vgQCqERy+dD9v1X4Zz/u+a2Vstm7RNRutyeLS4UokNm1C18CWduWX3DtHdnFDgBpV3bBxyAQYEUkcuePf6X0zx8cXGn9PR8fXCZu+8TepaJ2E4GspSP0R8UXdMwCmWtnSSCTW572jv3FRfETGKE8gHgEAkQskQceuqr0z380eLL193xs8DOidi4FUsUqCCSZuRNoOwinSW91JgLJWgmSkkXPpgoECEwi2aPPKvyzM8NtJ+95bvQkUTvbApmVx1oyQj/Lv1i9CUQoD6BaINxi0JDmBEIF14uLOHsqXQMxFEhluwYLBPAskTJp1OHcSJaNCEUgZcQkEBqUi0QaZ7f9CHe/T3OBZC1CMsz/eyqQQ4KaiZRxRZJgO6v+pqhAipEKhLblH4tUGGe3vQbiWiDl/alApolCIo9nHVG70ARySdpHd3R4gfsQyGiUYNCRJ1HJIMqeeFpErSOQolxJ1AKZTqIuUR4TvEok+cX7K29pYhXILDYEMv1ExiTOLqaXHqxdVBGjQIrwJZCJPGhYLpGqlOpCWIzABz0TCVEgl6Z97I5nF3UE0hvNf+sf35U/1ZHiYg0kNIFkKZDMXLc+BULDaokkHQMbSpHOQizvoQlWIiEKZJpZgUzftkzzf7uy993ek88sjNdAUgZGJRe+UCBSeQDNmYG0t+QhslZHFk9PekMgTYFRtWzEsxBPAgEClciyBXJp0sdutnZEILsFUgDkcXbpLYxTgVQRqECm1y7k/ZVfIOk4X+KiULJ4J66wHRCHQIAAJeJKILtZ+cUnjbPbXgNZVCBle2JiF4ikv9nF0TWfSdQVFgggP4t3KZgKpMNrha+UMjw4uBK72VqlQHwtoq7aDITHbzkrkGSQ/2p39mWR95pGuq6x/95SMZgLhNPiy0kqEBpmSxPIDc98m+zvVxDETGRaHp2SLMiZoWx9AQAeHcpi7yEJpGhPTCwCKSsnOKElTKK2LO/EBdwKpAwTgYgJYAYywbtEfuIgzh6jQIrwIRBuMVqbMtG0ZZNHAPJbmJAEkrUJyWD+fVQgh3iXiJTHRhfjMuqjw+UXVVMEIpUHIBDIMAHtGMTYd+Up1KqCxgdtIxRIcTsVyDTBS+Qx4QIqEKZALk6HB5kQGwKZXkwd7Ng/Ysl2BgTwJ5AiKeTtiQlGIEWSkErBRB6ZHdF4l8jxqx/BmYJt/bELZJo8geStf7hJospmFhOBcMqgklxJzAKp01YskG3ZP4YG8qc11rEkECAAiRQRskAme2JmBZKXQgXkSdSOyS1MV7oGYiaQKkIVyJFciWWBtLZlITJAnkQVC8RkZiFta1EgQKASCUUgRSlUwH4S1bVAyvbExC6Qo+9bLoV0vEja3rY/C0h60oLKzREIYCARIkoBnALwMDPfaH0kY6QCkcoDOBRIWQoV8JdE1RlIPpNaq0UCmX1qsrZjf6Na0pOlW1dVIIDZTOQNAO4DIA9rGLKIQMqCZCZFlX09hZEKJFceaQbkpVY9CcSkHuq0QIqqmwFyMRidSteV9akCkSGSCBFdDeAVAN4D4M0uBjItkB4X1794oP8L4j4vb/WwM6y+mKMUSBGeBSIKnAmzJdaPdBDKA6ghkJSA0fx4mioQQD4TeT+AtwEoPGGJiE4AOAEAx44dMxrEg0PZwU2PDmXrGgBwTthnaAJJ0gzZ9GNcTwLhlLG2KRONk8BZjAIpwLpAPDzGLaPyU0dENwI4y8yny9ox80lm3mDmjfV1+2ezSAVybnhFtAKZxZZAaEigIaHVSdDqJEiGVPmSCsToSIdIBcKt+e+tCuQQyaf0WgCvJKKXA7gYwJOI6F+Y+TVuh3bItECuSHrYzvIvVqk8gOUKZHpPTG2BzKyBmCRRpZjewuQVBSpqW0VoAslv1wCB8BLP4mXmmwDcBABE9DsA3upLIGWEKpBp8gQyGs7/lBt2wk+ixj4DmT0nxrpA+oJjMl3UQ7UoBylB5ER++/j9+NKZ+VqrsQmk7FS6rR2ZaIbijXTyKg6LCIRbxU9NYhfILFUCSQYjJDsG/2hBeUQjOAMoqKodRzCSCDN/EcAXnYxkhhAEUpRAPejPchI1FIGUtpMWEzJ53GtJIFlKSEZcKpBk6sDt1o48iSpGKhDxRjrL7RwQxExkFtsCsZ1CBewvoi5DIFmLkeSkVmMVSN7BUtIkqnQfDACQ8BZmFQUCBCgRU4Hsji4qbSdNoQJHd9yWEaNAighFIGlOzY4JbXHgTHYRq0DsEpREpgXSyYrl8ODelaL+zvUvF7+3r8e4dQWStyfGt0CqAmcm8XTpGogKxB/BSORHwiTqY8ISiROBzJ5Cl0esAsnDpkBMQmQmB3CXnXE7jQrEQrslEIxEJJgKREKQApneE2MphZqMD9s2EYMUH4uoh+1UIL6JRiK2BWJyrGVoM5DZ25fWjv3Hf6a3MFVlEk3KGYYoEG4lR6uQqUAOiEIiswK5NN3LXVANWSBpK8NomMwLpODsGGkSNe15PBPGcgYE8CcQ8ewDCFcgS6yrOk3wEoltBpKXQJ0gTqJ2hYWSVSClVAlkMruQCIT6I0AcYzeYLYgj7/ZP7bNFMBL502fdjVvuv+7I10IQiI0UKmBwCxOQQDgBKOd6iF0g1D/8s2THwcHaUky28ieC/9/SdpYJRiKzuBKIrRQqcHiLUkWMAqnbblkCocG85VrCQslJV7C/ZYJ0FpKNABL8/7FdC2RJO3bzCFIipgKpeoRrux4q4GAR1bJAiuSRtwAaikCKHvu2tuUXe7orrMjuSiASXArEw2wkOIlMBLKXFVc3A4CfdH9O1N+F/iWidhOBFB1lOU2sAsnv0267Vo9zT4ybxSRwxmkCEixSrrxAPBGURH7S/3lRu7N92Z4ZU4FICFEg06lV2wJpG5xyKq+HanCkw65wDcS2QEzOhFlhgQCBSUTCIgK5tDXA7jB/ZhO7QI60qxDI5DZGWuPUBHmIbF8gokJGKpD67ZZAVBJZpRnI5BS6IoHMBs5cJFFb1muB6AykvL/4BAJEJJFZgVzW2kNnOB84sy0QkzNhFhVI0ZGV0iSqkzNhFhBI2Z4YFUhVf3EKBIhEIr5mIAsLpCCFCthPooYikNJ2KpCK/kraTfrIOYaiuD95UxsEL5EQBGIlhQp4ewpjtkHu6O+LyiRGL5BpaQyEopHOAgLc3+KSoAo3vuu5/37k964E0h+lpa+tnYsxGiYHr1JS4QcmQoEUt7MrkHR3tLBAaDCaeyVbXSRb3X1RFL0mSAXipKhyeLcmixDsTMRUIMuuhwqE+xTmsD9ZO8CvQCZMV1+fa7dlcN+WJEAmeH8XAnFx/EPgBCmRiUC6o/LA2aNdmWg2e7KZivUnMECQApkclD3BhUAkgbP2jsE6hJBkV/gNUoFYIziJPChMoj6xJ5tZLCKQsj0xsQtkliqBJANGuyP/0K9tyW5LpBvpACDdHQCtBCiZqQAqEF8EJxEJLgVSRsgCmeyJKRPI9J6Z9o79D3N7125BZWAsEAEqEH9EJ5FFBVK0JyZWgSQ514I0cOZirWQRgXBKoJJHlyqQOIhKIk2egcwmUCdICyC7OdJB1s72DEQqD0AFEgLRSCRPIHnnxIQgkKIEKmBWD9W2GJYhEG4RKCe12iiBmEih4QIBIpGI9xlISQIVkKdQgfxzYvKIUSCF/alAiv+etD9JoSNPBC8RVwKpCpHZrocKGKyBBCaQvLWLoASSF1mXygMwiKevVhJVSlCJVQD4yG985ODXJgKpSqHOJlHLMA2RcVr9IYxVIHm4FAj1h4WvZLODZLNTnD5dBBfrGismm2BnIhOB9LPyW4VzHVmN1V1PtUCAMAUyveNWKg9ALpC1TdlMIN02eHMp0lmICsQKQUrk0a6sJuoF4S2MCqQYiUCSAaO9Lc92rG0JH8325Lcc1BWeC6MCWTqVVw0RXQzgywAuGrf/N2b+a9cDq8K5QKaPspwmYIFMqoWVB84OL4q1bftPDqS3MFKBHJFHKy0/OMqXQFZUHhMkP3r3ALyEmXeIqA3gbiK6nZm/6nhshYQ+AzE60qGGQPLKC0qTqCa3MK2e3TWQhQRShQrEG5VXDzMzgEnkqT1+eXv4XUcgrVaG4cyiaigCKTrH1qTsoe2iQ/tt3QqE2ylokHOejAokGkRXEBGlAE4DeAaADzLz13LanABwAgCOHTtmc4wHxDQDSQYlgTMDMUgKGu+/t6y/kARShAokLkRXBzOPmPl5AK4GcA0RPTenzUlm3mDmjfX1ddvjLBTIRa2jP8JdCYRGVPpqbydIBnTwsoGPRdTDtosLJGvPf6xUIM3F6OkMM18gorsAXA/gu26GNI+TGUhFChWwXw8VAJK9/FPo5vuU9We/bqr8AlmZGchogcO0pX9HmkRNgot0HVA5MiJaJ6Inj399CYCXAvi+y0Hd/qIPHPzaRCCjUVL5GnbWRAJxsQaSeDsXV9quWCDcPvrvdCaQfr/8tbe3/8qy6tdwuH8xV72q2imlSGYiVwH46HhdJAHwKWb+rNth7XNhnEStYrsjK7486LpbA8naXHgbI5XHfp922/mcgaRbXVE72hIuEhlVZLe8F4Y56P0rPpE8nfk2gOcvYSxHOCtMona682fP5OFSIGWELJCsTUgGLBJI0me0t+S3G9IkKok30nkWiFJIkIlVKa4FMjmFbpbYBZL0pwNn9qfriXAznQqkGUQrEZ2BjN8npyCyNIlqcguTdoVnvahAVo4oJdJUgRQGzgzqoYorsgsFIpUHUEMgrdb+4uYsKpAoiE4itQUyuyfGoUDKHuNKyx4C5WfcThOlQIpQgURDVBLxPQOpCpFJk6ip7KHF+L1l7UITCKcpaObxqAqkmUQjkTKBpGmG0Xh2YSKQslqoE9rbMtGIMyBjgUji7LEKJA/rAjFJjapAnBKFRIxmIB6SqKYCkb23rN0yBZK1EySDw36CE8js16UHcIsx6C+xnCkxkeaS063BSuT0De/Br9/+DnS6Fx3MMsoY7Ahrovb8LKLGLpBZxAK5IFz86QszKDpbCI5gJQIAWzuyyPtAupnOsUCK9sQ0QSBJf7+fdFv+j6FN4SKRC4EIZyHMGYjC3ZcSA0FLREIoAinsL0CBTIQwoS0sZ5juyhOr1BF+gzwLRKlP1BJZhkDy9sTYFohZRfajF0iacw22t+W3JfKK7CoQJZ9oJdKUGcisQIryIBdtyS8kbiegQfVFogJRbBClROoKZPYUOlcCqXqEa5JElWK7GlkUAjF4CqMCsU90EvE5A6kqJATYP4AbANpdRtYGkoqlCxVIVZcqEBdEJZFSgbQyYFKE2UAgknNxbSdRTQUiIUiBzO6JUYE0kmgkYjIDEYlhRyoa2dsuIpCqPTFRC2QWFYg/HNeEjUIig14LEETUbSdRXQqkimAFkqbAaCSXR1f2jza60C0nURkGNVWESVTKhNkTFxGVJReSDjpl84M/eCcGnbZIIOiGIZCio4Ol/bW7HK5ABiNgMAKd366uhTp52cZ6lF2pSxQzkUoCEUjd/qTyAOoLZPYRcLoprIfa2RW1AwDsySQinoWYCERroi6N+CWyBIHkxdlDEkheJkSaQgWAVHhrEoVAdG/N0olbIg2cgZQ9xjVJonKLQIJCRioQpS7xSiQSgaQ5NVCnkdZDNaG1K5ONCkSxQZwSsSCQ6dSqiUAk5+LKCyWb3cJI4uwqEGXZxCcRyzMQacLUdhLVVCCidioQxQNxSaRCIJPZhTRIJhaD5TWQaYGUxdmNTqXzJRChPAAVSFOJRyLdNJokahm2iwkBHgTSSoHhSAXim0C+B8FL5Myf/AWe/o/vFbX1vYhaRSMEMq6DynvCxOpIKg+P0XSTmqTS4vKJrKHXZKslCQUvESkhCWTUprmnMqEJhAY5K8Sb26L+xAJRwsXiLKYREglJIHm4FkhS8MRGmkIFAOzI1kCMBBLDLKRpVKV0HdwCRS+RZQlk9pwYlwKZrYE6TXvLYB2ilYCGggtUBbIaOFpDiVoiscxAqooJrW0ZnPgmJJHuslWBrAYOF2ErJUJETwPwMQBPwf7pPSeZ+QPORiTEhkCm98SYCKQqhQoYBM52hTtsx0c6zB4glYcKRDmC46c4kpnIEMBbmPmbRHQFgNNEdCcz3+t0ZCXYnoFI8yK2k6imApHgTSBSeQAqkGWyhMfAlRJh5kcAPDL+9TYR3QfgqQC8SKRKIJPAme0kaks4U1GBVKACWR5LypEYrYkQ0XEAzwfwtZw/OwHgBAAcO3bMwtDm8ZVEdS2QojKJKhBlYZYYRBNHWIjocgC3AngjM2/N/jkzn2TmDWbeWF9ftzlG/PjP3+ItiaozkBJUIFbh0aj2C6Ns/5Wx7GUB0UyEiNrYF8jHmfnTVt7ZASEJJG9PTHACGRwdIPcMFlClUhD+RGQHkiFpEtXkvcV9ypoZ1XcNFMnTGQLwYQD3MfP73A9pMUISSG47xwLJTaACwIW5SWMxRLKLflVnFeJK9Nl+QWtJfw0o4SiZiVwL4LUAvkNE94y/9pfMfJu7YZmxLIHMrl24FEjSL/47RknU2bNfCmDpZjoTgXichVjHRCA2+4sAydOZuwEEq8sYZiCJIFfS3pLdSohvXwCDNRAVSCkqkFIiT6wK25WIYTrOLhXI2pbsw7K2JZuBSDfSTQQiirOrQOygAqkkWonYnoFIg2RSgbS6dquRRTEDMbhAVCD+uf3cP1npJ0qJVAlkEme3nkSVHiq1gEC4lYKG+TMXFYgHGi4Qm0QnEW9JVIcCKUMFsmRMLnYVCIDAj9Gc5gc3vclfElUFUtKhCsRan5ES3UykipAEkrVp7smMCiRQbAlk0s9IGCKzHV4TlmW0SaMkEpJActu5FkjBAdpGSVTpkxXpReIRJylY6z16wuLB6I2RyLIEMju7cCUQKgmbAQAubIr6M8KnQDIGEsElKm0npSGpUSMsCgRoiERimIFUSgFAuik868XgUCnxLMSFQEzEIO3PJiuwXjGH7e8hGiARGwKZjrNLBSINkrUvCJOoXeHFPhFIuz23gW4WbwIx+aCqQJaHA4EAkUvE9gzEehK1U1FcdYyxQASoQEpQgVglWolUCWQSZ7eeRN0VtltAILyWFt/2qEDsoAKxTnQSSbtHj24ow3oS1aFASlGB2KGmQMT1SULCsUCAiMJmAHDve94kFoj1RdQFBcLt/G+xCsRinxJ0BuKMqCQiJRSBFKECsdinBBWIU6K7naliGQLJK6ysArHUdtkCyVj8o9QkvGa9NKO0v9HyhdkoiTRuBlIRUzc5VEp8AUiLL7O0P1kzo7aZdAK9gjMQDzRGItaTqA4EQnuC1Op5+0lUbwIxgDMG2UyiriJZJp+xWKQREvE1A5EGydILHVE7o1sY4SwkFoEoNfG4wTF6idgUyEWb8krrrU3ZjljxZjqhQKK4hTFABWIBzzuko5ZIlUAmC6AXbdr9Jqe2d+NOC2StDfTzb5FUIMocAZRYiFYi0iCZVCDSWxinAilBBaLMEYBAgAhzIt/6wJu8JVEXFQi3C1ytAlEWJRCBABFKREooAilEBaIsSkACARoqkWUIJC/OrgIp645VIDYITCBAxGsiRegMpASPArHZlmzXEZWG1wx+5NouzUjS6L6HspWNkkhdgWStBMnUyXIuBCKRg4uaqNIPl8/ZwsrNVBpSmrExEvE1A0nPy86x4PMXZO1cCER4cVq/iDkDSLwxRd62CTjYFMhZ5qVcQSMkUiaQ0Roh7Y9LHwoEIk2hAkAiTaLuyNo1TiAu2sZA1ezCkUB8Eb1EJDOQtS15ElWK+BZGBWKvbRNomECASJ/OfOOWNwPYD5Klfa58SUk7ltdApgRC7XZhMxXIitBAgQCRSgRwkER1KJAyVCArQkMFAggkQkS3ENFZIvruMgZkE+cCKUqiqkDK26lA6ncZiEAA2UzknwFc73gc1tEZSEl3PgWyalgWCGdZUAIBBBJh5i8DeGIJY7FGHYFk7XTua94EwqwCiRkHAgkRa09niOgEgBMAcOzYMVvdGtOYGYjJB1AQJBPLQ3ixS/szSpeK39vuUp54jC7KMkpjNOKw4PJFY00izHwSwEkA2NjY8BI9FCdRbQqk3QI/ek7UnzjKbjJbWMWf8KuG7UPMAauzpOhzIhNszkASYQoVAPgJYRLVo0B8zUJMaEyNVZPkrc3Dzk2wfJsV7SPeu29968GvRUnU8z20z/eQ9IeVLzHSWxgVyNL79ILtWWEEAgFkj3g/AeArAJ5NRA8R0eusj6IGaxeGSPpZ5UuK6RoIrRWHyAAVSFV/KpACXJy/4+gQr8rbGWZ+tZN3tsDaBdmsobUrPNbB9iKqCmRp/XklBoE4JNrbGSkLCaQoRAaoQCygAikhMoEADZeIzkCW2J8QFUgJRd+bRRecl3QGcWMlogJZYn9CVCAlRDgDmdBIieQJhJeRRFWBLK0/r6hAjtCYnMgEHzMQ7vZkfRl8+EKXghTb9VUB+zVW5QlYWWpUPD5pOVSTsqkmleQsEfVM5Av/9fYjv7cpEH78/P5rr1/5UpQg8CAQoEEzEYlAkse3rL+vfC9Mg2Yh0lTmqtVNNcBbQtfBNonoJdK+IDzR24Qd4bEOqyoQm+1WEG/rQ47+n6zMj4lkW7huMRYItcr9qgKx0G4FaZpAgBWRiKlAqlCBWGi3goQkkDt6H7fWfeMlckQgJcWSVSD131cFUkxIArFNoyWiMxALqEBq02SBAA2QyB2n3537dRWIBVQgtWm6QIAGSCQPFYgFVCC1WQWBAA14xDuLTYFkHZlkjAhdCrb7M3lfYabE9vfGegJWmDCVthNDfqTVKIlIBMKPn7f/xhqqUmJCE6vzuEiiZg5uYRozC/GJJ2FLE6bB14p18FlYqR+fvLUDVITIABVIkHg8OU/6/yT4ncqaWK0Hb8kquKtAAsTjv0MFUk0jJHL7mZtL/9y6QExQgdRDBVIfx9/DRkikDCcCsfxkQgVSgAqkPppYrYcKJGJUIPXRxGo9VCARowKpjyZW66ECiRgVSH00sVoPiUB4NHKwiCqNKTbkYneBp6i97YSpvGarbWxHYGU0RiLZucflbbtdTZgqq4tlCTfmSrqj8zH7ndr+yaizEGWCr8+WJlbrk3WFNVkj+J+nREqDBAKsmERUIIp3GiYQYIUkogJRvNNAgQArIhEViOKdhgoEEEqEiK4nov8hoh8S0dur/0Y4qEAU7zRYIIBAIkSUAvgggBsAPAfAq4noOa4HtgizT2hUIIp3Gi4QQDYTuQbAD5n5AWbuA/gkgN9zO6z6qEAU76yAQABZ2OypAP536vcPAXjhbCMiOgHgBAAcO3bMyuAWwUleRFGUQqwtrDLzSWbeYOaN9fV1W90qihI4Eok8DOBpU7+/evw1RVEUkUS+AeCZRPR0IloD8IcA/sPtsBRFiYXKNRFmHhLRnwG4A0AK4BZm/p7zkSmKEgWiXbzMfBuA2xyPRVGUCFmJxKqiKO5QiSiKUguViKIotVCJKIpSC2K2X3SWiM4BeNB6xzKuBPBTT+8tJfQx6vjq0dTx/RIzzyVJnUjEJ0R0ipk3fI+jjNDHqOOrx6qNT29nFEWphUpEUZRaNFEiJ30PQEDoY9Tx1WOlxte4NRFFUZZLE2ciiqIsEZWIoii1aJREQi8oTUS3ENFZIvqu77HMQkRPI6K7iOheIvoeEb3B95hmIaKLiejrRPSt8Rjf7XtMeRBRSkT/TUSf9T2WWYjoDBF9h4juIaJTVvpsyprIuKD0/QBeiv0Sjt8A8GpmvtfrwKYgohcB2AHwMWZ+ru/xTENEVwG4ipm/SURXADgN4PcD+/4RgMuYeYeI2gDuBvAGZv6q56EdgYjeDGADwJOY+Ubf45mGiMUSXzsAAAG9SURBVM4A2GBma2G4Js1Egi8ozcxfBvCE73HkwcyPMPM3x7/eBnAf9uvrBgPvszP+bXv8CuqnIBFdDeAVAD7keyzLokkSySsoHdRFEAtEdBzA8wF8ze9I5hnfKtwD4CyAO5k5tDG+H8DbAIRa9p8BfJ6ITo+Lq9emSRJRLEBElwO4FcAbmXnL93hmYeYRMz8P+7V+ryGiYG4LiehGAGeZ+bTvsZRwHTO/APvnSL1+fItdiyZJRAtK12S8znArgI8z86d9j6cMZr4A4C4A1/seyxTXAnjleN3hkwBeQkT/4ndIR2Hmh8f/PQvgM9hfBqhFkySiBaVrMF60/DCA+5j5fb7HkwcRrRPRk8e/vgT7i+jf9zuqQ5j5Jma+mpmPY//z9wVmfo3nYR1ARJeNF81BRJcBeBmA2k8KGyMRZh4CmBSUvg/Ap0IrKE1EnwDwFQDPJqKHiOh1vsc0xbUAXov9n573jF8v9z2oGa4CcBcRfRv7PzTuZObgHqMGzFMA3E1E3wLwdQD/ycyfq9tpYx7xKorih8bMRBRF8YNKRFGUWqhEFEWphUpEUZRaqEQURamFSkRRlFqoRBRFqcX/A9osdsUCoCqAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import fenics as fe\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "# --------------------\n",
    "# Functions and classes\n",
    "# --------------------\n",
    "def bottom(x, on_boundary):\n",
    "    return (on_boundary and fe.near(x[1], 0.0))\n",
    "\n",
    "\n",
    "# Strain function\n",
    "def epsilon(u):\n",
    "    return 0.5*(fe.nabla_grad(u) + fe.nabla_grad(u).T)\n",
    "\n",
    "\n",
    "# Stress function\n",
    "def sigma(u):\n",
    "    return lmbda*fe.div(u)*fe.Identity(2) + 2*mu*epsilon(u)\n",
    "\n",
    "\n",
    "# --------------------\n",
    "# Parameters\n",
    "# --------------------\n",
    "E = 70.0e6  # Youngs modulus\n",
    "nu = 0.4999  # Poissons ratio\n",
    "lmbda, mu = E*nu/(1 + nu)/(1 - 2*nu), E/2/(1 + nu)  # Lame's constant\n",
    "\n",
    "l_x, l_y = 5.0, 5.0  # Domain dimensions\n",
    "n_x, n_y = 20, 20  # Number of elements\n",
    "\n",
    "# Load\n",
    "g_int = -10000000.0\n",
    "\n",
    "# --------------------\n",
    "# Geometry\n",
    "# --------------------\n",
    "mesh = fe.RectangleMesh(fe.Point(0.0, 0.0), fe.Point(l_x, l_y), n_x, n_y)\n",
    "\n",
    "# Definition of Neumann condition domain\n",
    "boundaries = fe.MeshFunction(\"size_t\", mesh, mesh.topology().dim() - 1)\n",
    "boundaries.set_all(0)\n",
    "\n",
    "top = fe.AutoSubDomain(lambda x: fe.near(x[1], 5.0))\n",
    "\n",
    "top.mark(boundaries, 1)\n",
    "ds = fe.ds(subdomain_data=boundaries)\n",
    "\n",
    "# --------------------\n",
    "# Function spaces\n",
    "# --------------------\n",
    "V = fe.VectorFunctionSpace(mesh, \"CG\", 1)\n",
    "u_tr = fe.TrialFunction(V)\n",
    "u_test = fe.TestFunction(V)\n",
    "g = fe.Constant((0.0, g_int))\n",
    "\n",
    "# --------------------\n",
    "# Boundary conditions\n",
    "# --------------------\n",
    "bc = fe.DirichletBC(V, fe.Constant((0.0, 0.0)), bottom)\n",
    "\n",
    "# --------------------\n",
    "# Weak form\n",
    "# --------------------\n",
    "a = fe.inner(sigma(u_tr), epsilon(u_test))*fe.dx\n",
    "l = fe.inner(g, u_test)*ds(1)\n",
    "\n",
    "# --------------------\n",
    "# Solver\n",
    "# --------------------\n",
    "u = fe.Function(V)\n",
    "fe.solve(a == l, u, bc)\n",
    "\n",
    "# --------------------\n",
    "# Post-process\n",
    "# --------------------\n",
    "fe.plot(u, mode=\"displacement\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To recover the correct physical behavior,\n",
    "we switch to mixed variational problem formulation.\n",
    "In the previous example we supposed the stress tensor is\n",
    "\\begin{equation}\n",
    "    \\boldsymbol{\\sigma}=\\lambda\\mathrm{div}(\\mathbf{u})\\mathbf{I} + 2\\mu\\boldsymbol{\\varepsilon}(\\mathbf{u})\n",
    "\\end{equation}\n",
    "and $\\mathbf{u}$ was the unknown field.\n",
    "Here we define $p=-\\lambda\\mathrm{div}\\mathbf{u}$, the hydrostatic pressure, as an additional unknown.\n",
    "Hence the previous strong form\n",
    "\\begin{align*}\n",
    "\t\\text{div} \\boldsymbol{\\sigma} + \\mathbf{b} &= 0 && \\text{in } \\Omega, \\\\\n",
    "\t\\boldsymbol{\\sigma} \\cdot \\mathbf{n} &= \\mathbf{g} && \\text{on } \\Gamma_{\\text{N}}, \\\\\n",
    "\t\\mathbf{u} &= \\bar{\\mathbf{u}} && \\text{on } \\Gamma_{\\text{D}},\n",
    "\\end{align*}\n",
    "is replaced by more equations\n",
    "\\begin{align*}\n",
    "\t- \\nabla p + \\mathrm{div}(2 \\mu \\boldsymbol{\\varepsilon}(\\mathbf{u})) + \\mathbf{b} &= 0 && \\text{in } \\Omega, \\\\\n",
    "\t\\mathrm{div}\\mathbf{u} + \\frac{p}{\\lambda} &= 0 && \\text{in } \\Omega, \\\\\n",
    "\tp \\mathbf{n} + 2 \\mu \\boldsymbol{\\varepsilon}(\\mathbf{u}) \\cdot \\mathbf{n} &= \\mathbf{g} && \\text{on } \\Gamma_{\\text{N}}, \\\\\n",
    "\t\\mathbf{u} &= \\bar{\\mathbf{u}} && \\text{on } \\Gamma_{\\text{D}},\n",
    "\\end{align*}\n",
    "Its weak formulation being: find $(\\mathbf{u},p)$ s.t. $\\mathbf{u} - \\bar{\\mathbf{u}} \\in \\mathbf{V}$, $p \\in P$ and\n",
    "\\begin{align*}\n",
    "    -\\int_\\Omega p \\text{ div}\\, \\delta\\mathbf{u} \\,\\mathrm{d}V\n",
    "    +\n",
    "    \\int_\\Omega 2 \\mu \\boldsymbol{\\varepsilon} : \\delta \\boldsymbol{\\varepsilon} \\,\\mathrm{d}V\n",
    "    &= \\int_\\Omega \\mathbf{b} \\cdot \\delta\\mathbf{u} \\,\\mathrm{d}V\n",
    "    + \\int_{\\partial \\Omega} \\mathbf{g} \\cdot \\delta\\mathbf{u} \\,\\mathrm{d}S,\n",
    "    \\quad \\forall \\delta\\mathbf{u} \\in \\mathbf{V}, \\\\\n",
    "    \\int_\\Omega \\left(\\mathrm{div}\\mathbf{u}+\\frac{p}{\\lambda} \\right) q \\,\\mathrm{d}V\n",
    "    &=\n",
    "    0,\n",
    "    \\quad \\forall q \\in P.\n",
    "\\end{align*}\n",
    "It is an example of Hellinger-Reissner mixed formulation\n",
    "and the existence is assured e.g. by the inf-sup or Ladyzhenskaya-Babuška-Brezzi condition\n",
    "(https://en.wikipedia.org/wiki/Ladyzhenskaya%E2%80%93Babu%C5%A1ka%E2%80%93Brezzi_condition)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since FEniCS requires only one form, we rewrite the system above to a one variational equality\n",
    "\\begin{equation}\n",
    "    -\\int_\\Omega p \\text{ div}\\, \\delta\\mathbf{u} \\,\\mathrm{d}V\n",
    "    +\n",
    "    \\int_\\Omega 2 \\mu \\boldsymbol{\\varepsilon} : \\delta \\boldsymbol{\\varepsilon} \\,\\mathrm{d}V\n",
    "    - \\int_\\Omega \\left(\\mathrm{div}\\mathbf{u}+\\frac{p}{\\lambda} \\right) q \\,\\mathrm{d}V\n",
    "    = \\int_\\Omega \\mathbf{b} \\cdot \\delta\\mathbf{u} \\,\\mathrm{d}V\n",
    "    + \\int_{\\partial \\Omega} \\mathbf{g} \\cdot \\delta\\mathbf{u} \\,\\mathrm{d}S,\n",
    "    \\quad \\forall \\delta\\mathbf{u} \\in \\mathbf{V},\\, \\forall q \\in P.\n",
    "\\end{equation}\n",
    "We have to also choose the discrete subspaces carefully in order to keep the Babuška-Brezi condition valid,\n",
    "even uniformly with respect to the discretization parameter;\n",
    "see the list of suitable pairs of function spaces in [Hughes, T. J. R. - The Finite Element Method (1987), Chapter 4, pp. 201]\n",
    "\n",
    "Now we are ready to proceed with the implementation itself.\n",
    "First of all it is necessary redefine stress function using new variable $p$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fenics as fe\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# --------------------\n",
    "# Functions and classes\n",
    "# --------------------\n",
    "def bottom(x, on_boundary):\n",
    "    return (on_boundary and fe.near(x[1], 0.0))\n",
    "\n",
    "\n",
    "# Strain function\n",
    "def epsilon(u):\n",
    "    return fe.sym(fe.grad(u))\n",
    "\n",
    "\n",
    "# Stress function\n",
    "def sigma(u, p):\n",
    "    return p*fe.Identity(2) + 2*mu*epsilon(u)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Definitions of material constants and mesh remain the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Parameters\n",
    "# --------------------\n",
    "E = 70.0e6  # Youngs modulus\n",
    "nu = 0.4999  # Poissons ratio\n",
    "lmbda, mu = E*nu/(1 + nu)/(1 - 2*nu), E/2/(1 + nu)  # Lame's constant\n",
    "\n",
    "l_x, l_y = 5.0, 5.0  # Domain dimensions\n",
    "n_x, n_y = 20, 20  # Number of elements\n",
    "\n",
    "# Load\n",
    "g_int = -10000000.0\n",
    "g = fe.Constant((0.0, g_int))\n",
    "\n",
    "# --------------------\n",
    "# Geometry\n",
    "# --------------------\n",
    "mesh = fe.RectangleMesh(fe.Point(0.0, 0.0), fe.Point(l_x, l_y), n_x, n_y)\n",
    "\n",
    "# Definition of Neumann condition domain\n",
    "boundaries = fe.MeshFunction(\"size_t\", mesh, mesh.topology().dim() - 1)\n",
    "boundaries.set_all(0)\n",
    "\n",
    "top = fe.AutoSubDomain(lambda x: fe.near(x[1], l_y))\n",
    "\n",
    "top.mark(boundaries, 1)\n",
    "ds = fe.ds(subdomain_data=boundaries)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The crucial aspect of the mixed formulation is the definition of function spaces.\n",
    "The space for displacements is a vector-valued,\n",
    "while the pressure variable is a scalar.\n",
    "In addition, it is necessary to select different degrees of freedom for each variable to keep Babuška–Brezzi inf-sup condition fulfiled.\n",
    "For implementation of such spaces in FEniCS it is necessary to define individual element objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "P1 = fe.VectorElement('P', fe.triangle, 2)\n",
    "P2 = fe.FiniteElement('P', fe.triangle, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where the first stands for a second-order polynomial vector-valued element\n",
    "and the second for a first-order polynomial scalar element.\n",
    "The next step is a definition of mixed element by calling the *MixedElement* constructor with an array of an individual elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "element = fe.MixedElement([P1, P2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, mixed function space is created by the *FunctionSpace* object with a mixed element as an argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = fe.FunctionSpace(mesh, element)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dirichlet boundary condition is defined as standard.\n",
    "In this case the space *V* repsesents the whole mixed function space,\n",
    "*V.sub(0)* represents the vector-valued displacement subspace and *V.sub(1)* represents the scalar space for pressure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Boundary conditions\n",
    "# --------------------\n",
    "bc = fe.DirichletBC(V.sub(0), fe.Constant((0.0, 0.0)), bottom)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to access the individual subspaces of subspace.\n",
    "For example *V.sub(0).sub(1)* is a space for displacements in the $y$-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "bc_2 = fe.DirichletBC(V.sub(0).sub(1), 0.0, bottom)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*fe.TestFunction(V)* in the mixed formulation returns a test function on defined on the whole mixed space.\n",
    "It is possible to extract the test functions from individual spaces by the function *fe.split()*,\n",
    "but a simpler solution is calling *fe.TestFunctions(V)* instead of *fe.TestFunction(V)*;\n",
    "it returns splitted functions of individual spaces. The same holds for trial functions.\n",
    "\n",
    "The weak form below is valid for nearly incompressible continuum.\n",
    "For a full incompressibility, it is neccessary to omit the $p/\\lambda$ term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Weak form\n",
    "# --------------------\n",
    "u_test, p_test = fe.TestFunctions(V)\n",
    "u_tr, p_tr = fe.TrialFunctions(V)\n",
    "\n",
    "a = fe.inner(epsilon(u_test), sigma(u_tr, p_tr))*fe.dx\n",
    "a += -p_test*(fe.div(u_tr) + p_tr/lmbda)*fe.dx\n",
    "L = fe.inner(g, u_test)*ds(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution is obtained in the standard way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# Solver\n",
    "# --------------------\n",
    "sol = fe.Function(V)\n",
    "fe.solve(a == L, sol, bc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To postprocess the solution, it is sometimes useful to access to subspaces by calling *.sub(n)*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASYAAAD4CAYAAABBh0sxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2dXahk2XXf/+tU1b39OR94Wspgadx5SAzCD5ZpnAeZYAQOiiycPEZg58UwhNggSxG29WRMHhK9yA7EEAZLGGNjEyKZBMdKLMgYMRBL6pmMFUljByWWsCaadM903+77VVXnY+WhqmburVt1zv6fW/uctavWD5qenrvOrn1P1fnV2vusvY+oKhzHcSyR9d0Bx3GcZVxMjuOYw8XkOI45XEyO45jDxeQ4jjmGMRp95pln9Pbt2zGadhxnS3j55ZffVNVbq34WRUy3b9/G3bt3YzTtOM6WICLfXfczH8o5jmMOF5PjOOZwMTmOYw4Xk+M45nAxOY5jDheT4zjmcDE5jmMOF5PjOOaIUmDJUL3xd6n4x9UpFX9QFVT8o4o7JQfVVSoeAB5XV6j4k2qfih/riIuvuPgy8vfZABV/jHDH7EtOxV/PJlT8zWxMxT9Bxt/MplT8Uxm379rTGfeZA4Crz/41fcw6PGNyHMccLqaOOaiuUfGHZEZ2SGZjhyUZT7Z/RLbPxp9U+3xGSWaIEzYD1T0ynm2fi++CXMuNttf7UI4h9WGcS2mz8bGFBKQvpVwXuUe4OEpww75NSwlIREwupKZ4W0IC0pcSKyTAspTCsCCkBebF5FJqirclpdhCAuxJKbaQ2hyTspQA42JyKTXF75aUWCEB6UupjcRSlxJgWEwupaZ4l1IdPnRrhhXS7DXiSwkwKCYXUlO8LSEBuyclVkizY9KWUldCWmBKTC6lpnhbUvL5pJB4n09qgxkxuZSa4ttL6c3iZm0sezFXKlT8vfwJKv7N/AYVDwBvCnfMk1dPzElpV+eTVtG7mF4vD6n4B+Wiy2Fdjy0YADgoSSlFzmQmLep1Uua03MNoGP7FNtER/nryLuo1DgZx32M2MzymJXlMxZ+QS3AAgPv6qad3MVmClRIrpEV86FqwR/P4LHAd2FExuxhGWdi33KNiJu39LDz7PC72cXUQvk7rqNzHtUH4urSTcoQbg/CL4rTkLlA2S5rMM+1rRJ9yHVCvwcaPdYSbCBdxSWZVU7I/AFCCy6KbSEZMBxV3smJnSm2lFMojMn4hpeD2C+78HBfcN/pRyd/aPyk5aXQlpVBYwbQ5hh3uxZbSpoW0IAkxuZTqSV1KsYUE2JNSG4ntipSABMTkUqrHpVRPm/oka1KKLSTAlpQA42LaJSnFFhLgUmqMJ4UEpC8lC/NJqzAppl0SErB7WRKQvpR8Piku5sTkUqondSn5fFIzbeqZtklKgDExuZTqcSnV4/NJYViXEkCISUQGAO4CeF1VP7LpjqQupe9OngmOfVxyfecvaC5+Sl6cxaS+knyZt6ho4KG0KHLNyc/DgItnJXw0uFyB5btHjxqP2Zb5pFUwn8iPAXgNmy3wxN+Q3+pvldzSgwdkPJ31FFw8K6XH+RUMAwsmAeCQjD/Kr2BvEF5geVLsYY9of1wOcYVovw3jcoQbw/ACyGk1xD7Rp7waAMQ1PamGVJHoKv5f/mTtz9llRLEfaFFB8D7qiHqClCsi7wHw0wB+e4OvTWNuKNaBlBgOyfgjMv6kIJdBlFwmNi6HLY6Jmx3mZCbfZmgYe3jIxrNPwakiZFWhZ/E3AfwygLU5vIg8D+B5AHjuuecu37MlXEr1bIOUuPi4QgK6kVLM+DbHWJASEJAxichHANxT1Zfr4lT1BVW9o6p3bt26tbEOAi6lJlxK9bTJklKXEpuJlcjMSAkIy5g+AOBnROTDAK4AeEJEfk9VfzZar85gSUoupHpYwbQ5JvWhG9CNlBgsCWlB41lV1U8B+BQAiMhPAvhkF1KyJCTApdRE7CxpdkzaUupiKLYNUgKM1TEtcCnVs2tS8vmkOMdYlRJAiklV/wzAn0XpyRyXUj0upXpiZ0mAPSnFFhLQrZQAYxlTylJyIW3+GGtS8vmk7jAhpk0Iad0e2g+L63R/2GNoyZBLP2JLYzrmzn+bZRAMjyb8zgmPM7IyO+fO6RH5nrHLadj4EzKef8BBvfD+1vCAao+ldzG9XjxNxd8vNlp4fgF2fVlsKR3l+8gkfAP5o3wPwyxsK14AOM73grfiBYBxMcJoEB4/LYbYG4ZXWU+LIa4Mw7fincwlfC0L3+43rwZUNfq0GuAKwuOLihN3EVn0MXijeOrcvzf9ZdW7mBjYRa+PyOGVRSlx8dy36DEZPy7IoVVBZm5k/ITMDNvMJ03JY2JLiX8qSov5JPY1Iog1GTHtkpRSFxKQvpRYIQHpS4kVEhBvWJ+EmGJKyZKQgPSlFFtIgD0psUIC0pdS7HlG02LapSwJcCk1wQoJsCelNvNJuyYlwLCYXEpN8S6lOnw+KQyLUgKMismlVBcbV0jA7knJ55MCX6PDu4fmxGRJSpaENIvfrSwJsCclH7p1gxkxWRISsDtSOpnuoSjJDyo92cv9rsfCZ3kn0/NSffLquPEYa1JihQRsp5QAA2L63vQHqPgHZFX2Wzlf+f2I3D/68ZRdLkJWBZPx05IfmqRMXgywv3e+qPTR6fr3pCgHON0j91Ens8nYy2/Y5TFFB4uIN0l6JaeR4dekxc18WCmNyeHSeEoOx3JyaFWQF0Tk+IKUdptN49iMks3EKuXWr7WZT+qb3jMmBnb49ph83BEjJRdSPaww2hzThZQY+CGuPSmVZPuxSEJMloQEuJSasJYlAfak1KYoc1ekBCQgJpdSPS6lelghAfakxAppdky6UgKMi8mlVI9LqZ7YWRJgT0rt6pNsSQkwLKaUpeRC2vwx1qTk80lxMSmmmFLyLKme1LMkwJ6UfD6Jx5SYUs6SZvEupT7jfT4pDOtSAgyJadek9OAkvNJ9MuWK74qC+6COycpvrbgPdq6k8MYtdhEgJZzncSW5KHJ96spp8DG7Op+0it7F9P3pk1T8Qc4tXXk45eIB4ICs5D4k9qg+JbMSVkp5PoAQW/GW0wFkEB5f5QPIIHzr3i7QUgDmdygFA+JtqEgRn+VgHPaFy2ZusZfHFC1uBGyS3sVkicdklnRCD93iS4mhnHLxFdm+zjO3UJHpPHOTYbj4tCSzDDaelFIbiSmdKcWNt0AyYjoqOWnQ691cSrW0lVJwPDucJAUD2JMSKyRgN6QEJCIml1I9LqV6WCEB9qTUblI8TSkBCYhpl6SUupCA9KXUZihmTUopC2mBaTFZkpJnSfXEFtLsGFtS8vmkeJgVU0wpWcqSAJdSc7zPJ8WIt4w5MVnKkgCXUhPWpLSLQ7e2x1jGlJhcSutJXUiASymo/R3Oks5iRkyWpdS0FORkTG59S1Ypl5O40pCCvNjI4jsFGS8t1sdNwwssAWBKnqOKjKeXrpDSfnr/dGulBBgQ0/3pDSr+gKzkfhhYebvgJN+jyvbHxoZimmcA8XmVPAOISnHkGVVlLYVAh2T8iOhPOf+bOE1SCsC8bWyWdIlK8VDePOX2sq+YDwX63763dzExHNHr3bhMhh26mZQSgZDxYNsnMzE2/m0phbbfYhI9tpTaSKzNRDpD31ICEhKTS2k90YUE7J6U2ggjtsQ6GLpZkBKQiJhSlpJnSZuNB2BOSq2yHpdSLY1iEpErAL6M2VMLhwD+g6r+WrQeLWFJSpayJGAHpUQKCXAphWBJSAtCMqYJgA+q6pGIjAC8JCJfVNU/j9w3SkqWsiTApbTpeJ9PioNFKQEBYlJVBXA0/+do/oe7N0tiKUsCXEqN7acuJZ9PCozv7sGZQa8kIgMReRXAPQBfUtWvrIh5XkTuisjd+/fvt+7QLkmpnA7oSW5GSpJnnJTyjJKSFEJJho0HYE5KWolLqQOCJr9VtQTwoyLyFIA/EpEfUdVvLMW8AOAFALhz506rjCplKaWeJQ2OiP7QGQYZP2mRwUxnf1dXid01fegWEN/P48Wpu3KqeiAiLwL4EIBvNMUzxJJS6NamZzmdcFnS9JSLr8hKbplk1KXNZiVZm/IBQ0gp0HnRZ3Ya9rtoTl6g5Dmdsg8hYBcRk4KpE96Ta/Yl70tKQNhduVsA8rmUrgL4KQCf3lQHHky4CtaHE140DBN2U/sJF09nPlMuPssFzGc2mwpZKf6OBEL7UxGV3FkuqPaISnFy6CbzhIpK6dn8n82Uos7YNvOoxZf3MpteHhPyqX8WwIsi8nUAX8NsjumPN9WB//gT/zY4ln48Ejk/ZElKMs0oKWW5ICOygGwqMymF9icXCNM+2Z/FMQxtpRSMYuultAlirNkLuSv3dQDv3/grk1iSksUsiYonhASAEhLQoj9kPNCRlBjarI9zKa0lkcrveFKylCUBLqUmogsJ8CwpgNg7G5gXEyOllIdugEupCZeSDbrYbsWsmCwN3QBbUkpdSG2OMTd0A1xKETEpppSl5FnSZuPbLC3x+aQ4dLkxnTkxNUlplJXnHqfsk9w18TsmJR+6xaGPnTJNicknuc/8/MxFnJGV0ENWSAUVDiEvzsGYjCf7DwBKnqOSLbBcI8naSnOXUmvMiGmVlPI1e0s/HnNV4vSe3CdcFoYxV8mdjYn5JPIiZS/qLAdXkJkDym/JHRUpAGXeYlIYdZnYukrziszeSrJP4y3e7xswIiammjv2fFLBbvyfE8tF2KFbB1LqIr4KfAuyeeZWEZJhs71NSmlTr7ENmdWmMSGmUCxJiR267aqUguPZ4SQrJCC+lNoIxqW0kmTEFFNKbbIkCpdSfXxsKbW4+KNLyYVUSxJiijnJnbKUrAmp1WsYk5IP3WxgWkyWhm5AXCl5ltSMOSn50C0aZsWUtJR86FYf7/NJTgMmxWRJSpaGboBLqRGfT9oKzInJJ7nXw0jJmpAAe1Ly+SS7mBLTtkxyZw3bsA5OSSFNuK5kbDwpGSEfEMDWY7K/LwBUsTM3cnlMGbi1bsnsUb5DmBDT0WSfio9eyX06oPbYzk7DL71synWFzmTI9gdTsvK74ONDCywBICsBZrvstyVJvAYr1laZVSCDeeU4+yxPNn4s5DXQMybExDBmM6UWTyMJve7oRbgGpUS1T2YZdDz9qKa48UBcKW0TmWx2zJqUmGJKKfrOAC6l+niXUpJsWkgLkhETI6U2WRJDTCm5kJqJLSUXUhixpAQkIqZUpeRZUkN8F1mPSykKMaUEGBdTF/NJDC6lDcb70C1ZYksJMCwmS1Ly+aQNxxuTkgspjC6EtMCkmHySezNtpy4kwKVkhS6lBBgU065IyWKWNByHxbYZJjFoiwLLLAcK8knXLqUwupYSYExMlie56yp/ByfkntynVDgGbHxkiVlDSgBD7rzSciUlJvTWtzYrxfuQEmBETMdkJXc+Jrs9tlPJTWcyZPYwIDMxuj+5QjNmzZ6iGhHxhaLcC4/vYujGXpsxr+XWleKkJyc9CWmBCTExFORQT4nMSiZx55NiSqkLIcWMB2ZSYogtpTbXZs/X82oSfG5BUmLaFSlZzJKixpNCAuxJyaSQgCSlBCQkJkZKjJCAuFLahqFb1HhjWRKwJVJKVEgLyL09+sGltKZ9l1J9vEspWUxnTDGHbgAnJUvzSQAnJWtCAuxJyeeTbGFWTC6lNW17ltR8jM8nJY9JMfkk95q2XUr18T502xrMiSlVKXUxdAutFh+MuauHrxInDcD2Z8Jf/YOcuzrZrXKLqy6lLmm0gIi8F8DvAng3Zlurv6Cq/yZGZ3qb5F76wLF7crOV3MNjMj5wqciukhVARewc22aoNzrhjuEhDSOrP8/bsod4iAkKAP9CVV8RkZsAXhaRL6nqtzbVCbaSuxpzUsrGxHzShHw8EjvUIyXDZDNspjGcZ22hqycW7Wvg6V8M9UIrvxfx5T5TKQ4UV4lK8Q4mxft8MsqiMnwZtvK76Dk9bLxiVfX7qvrK/L8PAbwG4Adjd2wdVR5xktuQlAbTbqQUq/3481UtnnQSW0oKk49ropftGRgeUpMuInIbwPsBfGXFz54Xkbsicvf+/fub6d0SqUppMOalxLCLUmKQqiMpGUMlTSkBhJhE5AaAzwP4JVV9vPxzVX1BVe+o6p1bt25tso8AOCnJJDMlJQaXUkN8B1nStkiJQmBGSkDgXTkRGWEmpd9X1S/E7dJFWCkxMFKyNJ8EcNJIXUiA0aGbQVLNks4ScldOAHwWwGuq+pn4XTqPS2lFrGdJjeyilGghASalBIQN5T4A4OcAfFBEXp3/+XDkfqHKB9GklE3EjJR8krshPvJ80s4O3QCzUgICMiZVfQkd/wqpTnIDPsm90XhrQzdgO6RkWEgLzFV+pyoldujGFOwNT0ghnbKV39wVTc8PkcWJgyn5PhVsFXeGvWPyd2iRkhTX6ENotlFKgDExxZ5PkirsXWH32GYruWNWEQ+mcaU0mCiUOPVZrqiGxLC5VJREvBQKEFv3tkIVba7oIfE+s4JR4Q6YXkVUKX3nn/7qRtszIaZqwnWD2ZMbWF8NuzKWHf4YmuTuQkoM9NCtJBf5souC26zW0Mhjty4yGOFWPwBAJdw1tmlMiImCfqRS+DvvUqqJ3zUpnRFSfi3SfoodSSn6MRF+j7TElKiUYk9yu5Qa4i8hpWhYlJKhO3vpiImQEiMkYHektGtCAvqRUuOk9zZIKfLvkIaYjEjJUiW3Z0kN8Ts8n3TZY6qrDTv0dfB72BeTS2lF2y6l2vjEhm75psoKEh66LWNbTDsgJZ9PqmfbpdSGlUPFmtcorq/4HQ1LCbAspkhS2pX5pFl8PClZExJgUEqJDN2ivMYlsSemLbjz1rSJ/oisOh4dk4I55eKHY27XfyH3/B6SwqAr9EtytUBFSpL2V/hnuLjOtn2GLZUSYE1Mm5bS/IKg9+QmK7PZ9mPCZlWslLJpRVZ+V9BB+AFSVsCQWAHASoaMjw27agAA/ZhazS5eJ+W1gPPQ4/IVG2IiKrMBIBvHO2P0fBKZhQ2JJ4bwTzuJLyUqPufipSTjE5dSG1rtIrCCwUlNQ2sedNAl/fcAwHf++SejtT0ghgUupfW4lPpnU1KqpcVrxOiXjYwpkCznPGpFSoyQZn2xIyUXkg0sSilmn5IRk0tpTTwhJc+S0sOikID4/UpCTFakFHPoNutLvCJLl1J67KqUACNzTHWskpKueerqYCK9SKlaofc6Ka2Kr5NSNTj/OzVKaSnepZQeFqXU6nFQLTGdMTGZUp2QquHFbVrrpKRDQM7E+yR3TTwhJRdSM11d+Jbmk1ZhVkybktLKeCPDN6tSyib1m21LzrXLVlgLUc9WXt+bHeNSCsPo0G0Zk2KKJaVNTnJnS9dm7Gru4TEng+FJ+E7+2biEdLFmLAKD4ym/fIVeukLGs6ypG8ovUxW+8nX4Q/qQEmBQTJuQ0qo1UzH35WYzH/phAcQSkwGbJbHLUaYFwOw3nZfAkJjKLCoqXtg1eGz21qOwR+s+s+wwrOZ0rts7qi8hLTAjpgG5J/Hw1EY5wM5JiYEd8pFr8LZZSl2xavmVGqj8NiMmBmaRZ6pSYhfiJi0lF5KzRHJiiiUlfueB3ZCSZ0lOHyQlphSlFHPoNuuLS2ltvEspWZIRk0tpVV+INW+Whm6AS2mL+PavfHzjbSYhJpfSqr70KCXV83fmIkrJhbSbmBeTS2m5Hz7JvTbepXRpNl471RLTYrIgpW0qB6ir2M7G5IblOSEw8pwI0zaAQdV8Dqvr++dfw4iUpje4W/P5Da59Nt4KZsVkSUrL6+xWsXfIZQMjtpL7KKfisxMunoIUhwWy43fe+GFEKeU3hhgdtXmoXWj7/dcYdYE5MbEb0Z+VktRcL2uraNfFH8X78A7Gke+8Nax1Ox9LCoyV0jQHhsTHrCoBBO75zQqGjbe2dAXA6IiLD+lSYTCrMiOm//2JT+CH/+VvBMcPyTeIIWqN0q5JKZQq7sLgtlIqb6zZY+dC+1zzljh7LfW9FGWBGTHFZGhkTimmlBghzeJdSmsxmCntGo0DVhH5nIjcE5FvdNGhJgaT2Z91m8Uts5BSFRB/VkqrNnM71+6pnpNSNaz/qtkZKU1zl5JzaUJm0n4HwIci9yMIem+kSJkSfedtl6TEYElKFTgpKVxKEWkcyqnql0Xkdvyu1LNtUqoGcuGx2THnk2bxhDhiDt2AuFKqiV0uG5j9T64rLqT4bOzeo4g8LyJ3ReTu/fv3N9UsAPKpJBP7Ulod61JaierGpLS6L1y4S6kbNiYmVX1BVe+o6p1bt25tqllaSlTbOyClbJLbkVJV8lJicCltDabvyrmU3kGHGaSoNl8OUMxff0yewJwQ0pSsKi/IZTRLDznQVcO1s7iQzGNWTH1LaTAF9g65CyR+NXf4BS6npAwYplNue92OkeP1H4iMlBK9dIVdr1eSl+CZ5qc3+68C/8tf3/zOAkCAmETkDwD8JIBnROR7AH5NVT8bpTdz2kgpC7jG9x7H++obsrsDEA8LAAAZh8fLmJzzYTIaMvvRPIcwld9lFf51aWS929vQi4i55qdPnK+Ib1oGNb2RYe8wsO2bXF9iE3JX7qNddGTBYAJI4HW1R1R/D0/Ip5ichIvGpbQaZYZ7wExKwY2nL6X8BiFs8tdd9UCOOvYOgfwmgkUWG1NDOWaZyZBYuJuqlBghzeITlRL5IExTUoqcJbWZ32KlBCyGrHaG5/0PUs/wzU+HjVddSqvil0RQNwc0nbqUNsFSX4obe7XhFqUkqma2gDmLKTGFcFZKjctGlqTUFL8sJR2sv7hNS6kO9g6ZdSntN9yBi4WxTEmqdlKyiqmhXBPWMyXNZOWjqndBSlHnkwDTmVIdVrMk6yQjJutSWodLaQWpTnIbzJJYUpASkMhQzqW0Kt6lFBu9dmbOyKXUKeYzphSlFCqkxdDPpbSqcUMXUUwpRRq6TZ54J+dISUgLTIupDyktLycZHZEb45N7bQtVzU2WuE/C45WIBUCtqVN2kW9JPJQz4EEE5+PjtT1kn+pCPDWmuDHClQerz/nkqfWXcaiUpk/wpQJf/0ycqm/AuJhCufIg/A3ee0xut0FAS4l4MolMIt5NY8XBSKkoIIPAPbwBaFWFV9NE3Z+JbJuczK+uBe50GMD+wcX3IyvCzzkASJVh8pSdOiazYhqOw5am7D+K90SK4XH4BehSuogWXLapTIaSsJT4LI/NxNoN3fYPwo+LLTFzk9/f/PTHsf8oTEqjY23c0vYse0cVkIXFL6QUsjm7S+kiLqWLSFVRUpJKo0upTf0TAFx5EHfeymzGtGn2iGd9mc2UMqm/cC4hJZEMqjXnKJKUKCEBkTeNI+K3IEtqdWcv3gDlHMmKaXQc/kZshZSaSDBT2iYpVdevrA23JqW2culKSkCiYnIpLeFSCou9smb5ipFMaVNCmtbdpVvRnfHT9TM6XQppQXJiciktEUtK5Ba75qW0isjzSQx1UsqfuLg4eBuzpLMkJSaX0hLbLqWEJ7kZLA7dmo756u98gm+UwNxduXXUSWn5zlyjlM7cmetdSvN3wKW0HGxESmVlRkpSqAkpdYHJjOmVf/dx/Ng/+423/91HpjQ8On/hM/ttg8iSAHCSIZaBsBXXyvSDqM6u2CevMPM4zBwOu6MB8TtmzPKcvT0MiWVF2enq2OkPXFt7TKpCWmBSTGcJkdLouML+A04Gy+LpjcnExMb+9Lq3WGgFgKtajtKNsoRgc9XZMdh76+TC/yuv7+HqG+FfBNMn93DtjZmRTt9lRwd2erKChZRGx5tV+fDRBKjZBG6Z7HAMDAJHvSenQBZ4YbHr0ybT4AJRJlviNnfjlvRoWQISeO7q6qguCzHU07KEXF1/+/8cb5+PQIkVJVC/0eWKtsOQaQFcD238IlfvNQutK3mZFdO1e+Fvyv7D8Atr+IiTQXZIrCQ+OZ39XZXNcmojpUBMSSk42I6UgiHPB/W8vHnbeuNqULhMySczT9utGb3+f7vJrM2KKZS3pSTS+AE8J6VSG7OmVlIKwaW0FBxJSuQEejQpkQ/wbJUpEbSVUjafeH/xT3+l1fEMSYvJXKYUgktpKXgLpFQ37IsoJVZIQDspZS0XBV8Gs+UCL33hk7U/71xK6+7odC2lNbeXXUpn2+WGbiYypbK8tJSKp9ffpcum5QUpjW/VDxOzQnuREmBYTHUwUmKJlimxeKbUjljzScBmpHRthQwMDt36EtKCpIdyF1gxz8RkS1GlxGRLiUmJvsBrpCTDS9yiT3iSOwQLQupifglIUEwbHcKdmQDvVUqq79Qy7bCULoUFKW14PunsHbnYUuo7Q1rG9FBueZ4p1rySZ0pniCUlrVxKZ9uNlCmtmktaxWJ+qc95pDqSyZjqpLT3xuN3/nF0sRq2FuqpIOzSikgXObG0Ita8D5cpxdtnnemHUA9QIBZfM19Wx+GfTzk8Cm93NMLwcXj8tccXJ8rHP/R0+OtFxryYbr72FhWvR8eQG9eDYquHB8iur7+TcaHtyQQS+EhqPTkJjq1Ox5C9sLkVnUyC52G0yIOrrrUsIaFV5ezwLRahFeUAtdRFixwyCHyfpzlk3T5PF2KnkGHYJaeTCWTU7eV55bsPa3/+xb/61x31xPhQjkXfuB+t7eroOLwfJ+HfitVp+DCSecSSFsRQL+acUiSU2SFAK07moc0yw2kmMyfe5+r4BBWRhenxCUB8Gb/N40P+mEtgXkxffO1fRWm3engw+zvgTXUp8bExYaUUHDo/byGZrhUpMSgZD2AmpI6lBCQwlAuFyZYWUgqKdSnRsbGIJSSgfaaUNQzjzkpJrtVnKmff56xhOiK6lJZk9MXv/xZ3/CUxnzEBzVnTspS0RiYupXmsS2kWWuTB502nuclMKXvqyfXtHp9cWkp9kHzG5JnSPNaltAgOD01sPgngMqVVQpJ3PbP+gDVC6jpbAgIzJhH5kIj8lYh8W0R+NXanVhFrrgkg32xGFkws88EnLqgLx7qUZqGJSYmZ5KazpJ7mkepoFJOIDAD8FoB/COB9AD4qIu+L3bEQ6rKl5eHczmZLLYsaXUrzWCNSCm53w8O2PrIlICxj+nEA347S8QgAAAO8SURBVFbV/6OqUwB/COAfxe3Was5mTT6Em8ducAi3EIFLaR67ZVI6N4wzmCWdJWSO6QcB/M2Zf38PwN9bDhKR5wE8DwDPPffcRjq3iphDOsdxbLCxu3Kq+oKq3lHVO7du3dpUs47j7CAhYnodwHvP/Ps98//nOI4ThRAxfQ3A3xGRvy0iewD+CYD/FLdbjuPsMo1zTKpaiMgvAvivmK2C/JyqfjN6zxzH2VmCCixV9U8A/EnkvjiO4wBIZEmK4zi7hYvJcRxzuJgcxzGHi8lxHHOIkk8sDWpU5D6A72684Xd4BsCbEdvfJCn1FfD+xiSlvgLx+/tDqrqyGjuKmGIjIndV9U7f/Qghpb4C3t+YpNRXoN/++lDOcRxzuJgcxzFHqmJ6oe8OEKTUV8D7G5OU+gr02N8k55gcx9luUs2YHMfZYlxMjuOYIykxWXgoQigi8jkRuSci3+i7LyGIyHtF5EUR+ZaIfFNEPtZ3n9YhIldE5Ksi8hfzvv56330KQUQGIvI/ROSP++5LEyLyHRH5nyLyqojc7fz1U5ljmj8U4X8B+CnMtvf9GoCPquq3eu3YGkTk7wM4AvC7qvojffenCRF5FsCzqvqKiNwE8DKAf2zx/IqIALiuqkciMgLwEoCPqeqf99y1WkTkEwDuAHhCVT/Sd3/qEJHvALijqr0UhKaUMZl5KEIIqvplAA/67kcoqvp9VX1l/t+HAF7DbL93c+iMo/k/R/M/pr9hReQ9AH4awG/33ZcUSElMqx6KYPLCSR0RuQ3g/QC+0m9P1jMfFr0K4B6AL6mq2b7O+U0Avwyg3fO0ukcB/KmIvDx/0EinpCQmpwNE5AaAzwP4JVV93Hd/1qGqpar+KGZ70P+4iJgdLovIRwDcU9WX++4LwU+o6o9h9jzJX5hPTXRGSmLyhyJEZj5f83kAv6+qX+i7PyGo6gGAFwF8qO++1PABAD8zn7f5QwAfFJHf67dL9ajq6/O/7wH4I8ymUjojJTH5QxEiMp9Q/iyA11T1M333pw4RuSUiT83/+ypmN0T+st9erUdVP6Wq71HV25h9bv+bqv5sz91ai4hcn98AgYhcB/APAHR6dzkZMalqAWDxUITXAPx7yw9FEJE/APDfAfywiHxPRH6+7z418AEAP4fZt/mr8z8f7rtTa3gWwIsi8nXMvrC+pKrmb8EnxLsBvCQifwHgqwD+s6r+ly47kEy5gOM4u0MyGZPjOLuDi8lxHHO4mBzHMYeLyXEcc7iYHMcxh4vJcRxzuJgcxzHH/wchqVA+w6CW+AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAEDCAYAAAC/Cyi3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAedElEQVR4nO3dfbBcdZ3n8ffn3oSQG54hIE8KK/jAUIBOREd8GFEsFllxp1YHa7C0UFNajgOujoMzu2u5tVqMuo5OzT5lgVksAR+AuK4okEVWhl0JQgyPQaAQNYDkcgkmAZLLvf3dP/o0dC7dfU7f7tPnoT+vqlT64XSf3+0+59O/h3POTxGBmVmdTBRdADOzYXOwmVntONjMrHYcbGZWOw42M6sdB5uZ1Y6DzazGJF0iaYukuzMs+3eSNib/7pf01CjKmAf5ODaz+pL0FmAH8M2IOL6P130SeE1EnJtb4XLkGptZjUXETcCT7Y9JermkayXdLumfJL2qw0vfD1wxkkLmYEnRBTCzkVsDfCwiHpD0euA/A6e2npT0MuBo4CcFlW9gDjazMSJpL+CNwPcktR5etmCxs4ErI2J+lGUbJgeb2XiZAJ6KiJN6LHM28IkRlScX7mMzGyMRsQ34laT3AqjpxNbzSX/b/sDPCiriUGQKNkkPS7orGQa+Le9CmdlwSLqCZki9UtJmSR8G/gz4sKQ7gHuAs9pecjbw7Rji4RKSDpC0TtIDyf/7d1jmZZI2JBlzj6SPJY9PSbpG0n3J4xdmWmeW8kt6GFgVEU/0+TeZ2ZiT9GXgyYi4UNIFwP4R8VcLltmDZh7tSvoB76bZF/gU8PqIuDFZ5gbgSxHx417rdFPUzPJ2FnBpcvtS4D0LF4iI2YjYldxdRpJNEfFMRNzYWgbYAByRtsKsgwcBXC8pgP8WEWsWLiBpNbAaYPmU/vDol4/HuMQSpS9jozE3Rsea33vXc09ExMpB3uO0ty2PmScbmZb9xZ2z9wA72x5a0ykHujgkIh5Lbv8OOKTTQpKOBK4BjgH+MiIeXfD8fsC/AL6RtsKs6fOmiHhE0sHAOkn3JQf+PS/5I9cA/MEJe8TlP+xY9tpZOVnZEfHamZ6fLLoII3PSyzb/etD3mHmywf+99rBMy04d9vDOiFjV7XlJ/xt4SYen/qb9TkREUkF6kYj4LXCCpMOA70u6MiIeT95/Cc0Dhv8+Ih5KK2+mYIuIR5L/t0haC5wM3NT7VeNhen7S4VYC4xRqZRQR7+j2nKTHJR0aEY9JOhTYkvJejybntr4ZuDJ5eA3wQER8PUt5UvvYJK2QtHfrNvBOmh17ZqXgUCu9HwAfTG5/EPifCxeQdISk5cnt/YE3Ab9M7v8HYF/g/KwrzDJ4cAhwczI0fCtwTURcm3UFZjb2LgROk/QA8I7kPpJWSbooWebVwPokZ34KfDUi7pJ0BM3m7HFA63CQj6StMLUpmrRnT0xbbpy5OVoc19bKLyJmgLd3ePw24CPJ7XXACR2W2Qz0PUTnwz2GxDvY6Pkzt27G45iMEXHNzapknmBrY1f6ghXkGtuQuRYxGv6crRcHWw680+XLn6+lcbDlxDtfPvy5WhYOthx5Jxwuf56WlYMtZ9Pzk94hh8CfofXDo6Ij4hHTxXGg5Wcu6vv5usY2QnXdiPLiz8sWy8E2Yt5Z07n5boNysBXAO253/lxsGNzHVqD2nXjc+98caDZMDraSGLfBBQeZ5cnBViILd/a6BZ3DrFzmmGCmMVV0MXLhYCuxOgSdw8yK4GCrkE4hUaawc4hZWTjYKi4tTPIIPgeY9UPSAcB3gKOAh4H3RcTWDsu9FLgIOJLmzHhnRMTDbc//PXBuROyVtk4HW805hKwELgBuaJsw+QLgrzos903gixGxLpk0+fm5ASWtAl40g3w3Po7NzPKWOmGypOOAJcklwomIHRHxTPLcJPAV4LNZV+hgM7MsDpJ0W9u/1X28NsuEya8AnpJ0taRfSPpKEmgAfw78oO09Urkpajam5mKSmfnU7qqWJ3KeMHkJzXlEXwP8hmaf3Ick/Rh4L/DHWQvaejMzs4EMYcLkzcDG1izvkr4PvIFmDe8Y4EFJAFOSHoyIY3qVx01RM8tb6oTJwM+B/SStTO6fCtwbEddExEsi4qiIOAp4Ji3UwMFmZvlLnTA5IuaBzwA3SLqL5lyi/32xK3RT1MxylWXC5OR+x0mTF7wmU6ega2xmVjuusZmNqbmYYHpun6KLkQvX2MysdhxsZlY7DjYzqx0Hm5nVjoPNzGrHo6JmY2ouJtk6t6LoYuTCNTYzqx0Hm5nVTuZgkzSZXCfph3kWyMxsUP3U2M4DNuVVEDOzYckUbJKOAN5Fc6IFM7NSyzoq+nWa1xvfu9sCyaWCVwMcergnEDEru3kmxndUVNKZwJaIuL3XchGxJiJWRcSq/Q/wmISZFSdLAp0CvFvSw8C3gVMlfSvXUpmZDSA12CLicxFxRHJZ3rOBn0TEObmXzMxqQdIBktZJeiD5v+P8oJL+VtLdyb8/bXtckr4o6X5JmyT9Rdo63WY0s7y1Jkw+Frghub8bSe8CXgucBLwe+Iyk1sXiPkRzdvhXRcSrabYce+or2CLi/0TEmf28xszGXuqEycBxwE0RMRcRTwN3Aqcnz30c+PcR0QCIiE6zXO3GNTYzyyLvCZPvAE6XNCXpIOBtNGtpAC8H/jRZ748lHZu2Qp8Ebzam5hoTzMxOZV081wmTI+J6Sa8D/h8wDfwMmE+eXgbsjIhVkv4EuITm5MpdOdjMbGBDmDCZiPgi8MXkNZcD9ydPbQauTm6vBf4xrTxuippZ3lInTE7ORT8wuX0CzWn4rk+e/j7NpinAW3kh8Lpyjc3M8nYh8F1JHwZ+DbwPmhMmAx+LiI8AS4F/kgSwDTgnIubaXn+ZpE8BO2ibi7QbB5tZhaycnE9fqGSyTJgcETtpjox2ev1TNM9VzyyXYFui3l/A9LzPJTVLU8UQK4tCamytL8wBZ9ZURIjNxwTbnls+8vWOQqFN0fYv0yFn48o1s+ErTR+bQ87GjQMtP6U83MNfuNWdt/F8lTLYwF+81Ze37fyVNtjAG4DVy8rJeW/TI1KaPrZuVk7Ou8/NKq+MgTYfYtvssqKLkYvSB5tZlZUx0MZBqZuiLa7CWxV5my1OJYLNzKwflQo2/wJaVXhbLValgg28wZhZukoOHnik1MqsKj++jZhg+649iy5GLipXY2upysZjZqNX2WAzKyP/4JZDpYPNG5GVibfHziS9V9I9khrJVXM7LXOkpBsl3Zsse17bcydJukXSxmSmqpPT1lnpYDOzSrgb+BPgph7LzAGfjojjgDcAn5DUuqLul4EvRMRJwL9L7vdUycEDs7Jxba27iNgEkMxn0G2Zx4DHktvbJW0CDgfuBQJozQq/L/Bo2jodbGZjqhHimdmlWRc/SNJtbffXRMSaHIqFpKOA1wDrk4fOB66T9FWarcw3pr1H5YPNh35Y0caktrboCZMj4kXT7fV4n72Aq4DzI2Jb8vDHgU9FxFWS3gdcDHSdxxRqEGxmVrxeEyZnJWkpzVC7LCKubnvqg0BrMOF7wEVp71WLwYMx+cU0qy01O+AuBjZFxNcWPP0ozYmSAU4FHkh7v1oEm5mVl6R/KWkz8EfANZKuSx4/TNKPksVOAT4AnJoc1rFR0hnJcx8F/qOkO4AvAavT1ummqJnlKiLWAms7PP4ocEZy+2ag47Bp8twf9rPO2tTY3Bw1sxbX2MwGUOUf1EZD7Mx+uEel1KbGZmbWkhpskvaUdKukO5JzuL4wioKZmS1WlqboLuDUiNiRHGdys6QfR8QtOZfNrNSq3Aytu9Rgi4gAdiR3lyb/Is9CmZkNIlMfm6RJSRuBLcC6iFjfYZnVySVFbpuZaQy7nGZmmWUaFY2IeeAkSfsBayUdHxF3L1hmDbAG4MQT9yikRufzRs2yixDP7azngRF9jYpGxFPAjcDp+RTHzGxwWUZFVyY1NSQtB04D7su7YGZmi5WlHnoocKmkSZpB+N2I+GG+xTIzW7wso6J30rzom5lZJfjMA7NF8DFs5VbPIREzS9eAmK3nUQSusZlZ7TjYzKx2HGxmlqssEyYny+0n6UpJ90naJOmPkscPkLRO0gPJ//unrdPBZmZ5yzJhMsA3gGsj4lXAicCm5PELgBsi4ljghuR+T7ULNo9WmZVLRGyKiF/2WkbSvsBbaE7oQkTMJmc6AZwFXJrcvhR4T9o6PSpqNq5CsCtz3SbvCZOPBqaBf5R0InA7cF5EPA0ckswUD/A74JC0N6tdsPkkeLNc5D1h8hLgtcAnI2K9pG/QbHL+2/aFIiIkpV5ko3bBZmajN4QJkzcDm9suiXYlL/SlPS7p0Ih4TNKhNC+f1lPt+tjMrHoi4nfAbyW9Mnno7cC9ye0f0JwNnuT/1Bqgg83McpVxwmSATwKXSboTOInm5MgAFwKnSXoAeEdyvyc3Rc0sV1kmTE7ubwRe1I8XETM0a3CZOdjMxlXAxGzHydcrz01RM6sdB5uZ1Y6Dzcxqx8FmZrXjYDOz2nGwmVnt+HAPs0WYnp+s/JVk1ICJnfWs29TzrzKzseYaW43NNKaG/p4HTjwz9Pc0GzYHW43kEWRp63DQWRk52CpuFGGWdf0OOSuLWgXbOFxksugg68UhZ2VRq2CrszIHWiet8jrgSixgclfRhciHg63kqhZoCzngrAg+3KOkZhpTlQ+1dq2/p05/U1H2n1hWdBH60se8oudJujtZ9vy2x7+SzDV6p6S1kvZLW2cuwTZJPa/xNArjsPOPw99ou0mdV1TS8cBHgZNpzil6pqRjkqfXAcdHxAnA/cDn0lboGltJjOPOXvW/uYjBqqrV1iDbvKLAq4H1EfFMRMwBP6UZhkTE9cljALcAR6StM7dgq+IXUISq79zD4M/AaNbq3izpQElTNC8ZfmSH5c4Ffpz2Zh48KJB35t3NNKY8yNDDsCsLCpiYzbx4zwmTB51XNCI2Sfpb4HrgaWAjsNvJuJL+BpgDLkt7v1yDbf+JZWxt1HQ8eQAOtO48ilpaPSdMHsK8okTExcDFAJK+RHOuUZL7HwLOBN4eEakTJtemj60qB+c61LKpyudUle2uCiQdnPz/Upr9a5cn908HPgu8OyIy/eKlBpukIyXdKOneZBj2vH4K6762F1RlZy0Lf14vqPJ+1Me8oldJuhf4X8AnIuKp5PF/APYG1knaKOm/pq0zS1N0Dvh0RGyQtDdwu6R1EXFv2gutyTvo4rnfrfr6mFf0zV1ef0ynx3tJrbFFxGMRsSG5vR3YBBze74rGlUNtcOM+alrl2lpR+ho8kHQU8BpgfYfnVgOrAY483P0O47wj5sW1tyFrwOTOoguRj8yDB5L2Aq4Czo+IbQufj4g1EbEqIlYddODuwTZuvzgOtfz4s7UsMgWbpKU0Q+2yiLg63yL1r0wjU97x8ufP2NJkGRUVzWNLNkXE1xa7onGotXmHG50y9buV6YfVmrLU2E4BPgCcmgy1bpR0RtqLxk1ZdrJxU/fPfRwqBHlIHTyIiJvBl+vope47V9l5UMEWqs2ZB0VxqJWDvwdrV/mT4Ivs3/DOVC5F1tzymEA572aoAiaznwRfKSOtsdWpv8ChVk7+XgzcFF0U7zzl5u/HHGx98k5TDVX/nurUuilCpYNt1P1rVd9Zxo2/r/FV6WAbJe8k1eTvbTyNfFS0ilfV9c5RbaMaLc1jZDRXjfpOmOwaWwqHmln1VDbYRtG/5lCrj1F9lz5v9MWyTHic5Urdkj4tKSQdlLbOygabWb/8Q1WYLBMet67UfRzwBuATko5rPSnpSOCdwG+yrNDB1oV3gnry9zp6WSY8znCl7r+jOaFL6gxV4GDryBt/vfn7XZSDJN3W9m/1It8ndcLjhVfqlnQW8EhE3JF1JZU8VzTPfgxv9OPBVwRpniu6ZGemChCkzCuaZcLkLBMeL7xSdzIr/F/TbIZmVslgy4tDzYZh0MM+qnjWQdqEyVkmPO5ype6XA0cDdzSvecsRwAZJJ0fE77qtz8GWcKiNH9faRqNtwuO3dpvwuNuVuiPiLuDgtuUeBlZFxBO91lm5YMujGVq3UJuZ3yvX9z9wckeu7z9KDreR+AdgGc0JjwFuiYiPSToMuCgizuCFK3XfJWlj8rq/jogfdXzHFJULtmGreqjlHWJZ11nlsMsj3Cp3FkKOuk143D5hctYrdUfEUVnWWalgG/eDH4sIsayqHnZlqrltbeyqZD9bmVQq2IatCrW1ModZmvayVyHkyhRuo6AGTO7KPCpaKWMbbGUOtSqHWTdVCblhhpubo8WpTLANsxlaxlCrY5h1U5WQs+qqTLDV1TgFWidlDLlxa5LW0didUlWW2trM/F5jH2oLlekzGdZ2Mu4DXkUZqxpbGUKtLDtumbU+o6JrcK65VVclgq0Ov3oOtP6VJeDqqs9zRStlbJqiRdXWytS8qqoiP8NhbDd1+GGumkrU2AZVRKg5zIavqBqcm6TVU/pgG/TXbtSh5kDLXxEB53Crllo3RUcZam5yjt6oP/NBtic3R0ertsE2qlBzoBWvKuFmo1PqYFvsr9woNj4HWrmM8vtY7PbVz/Zctbl3y2bkfWx5f2GjCjUrJx8i0odGMPlsPc9lLW2NrYx9Eq6lVUfe39Uoam22eKnBJukSSVsk3T2KAg0iz9qaA62ayhhu4ybLhMnJcg9LukvSRkm3LXjuk8l73CPpy2nrzFJj+x/A6Zn+giFZzK9aXhuZa2nV5++vcFkmTG55W0Sc1D4jlqS3AWcBJ0bEHwBfTVtharBFxE3Ak2nL1ZF3iPrI6wfKtbZ0WSZMTvFx4MKI2JW835a0Fwytj03S6tZkqk/MdO6QzGvgII+Ny6FWT2UIt4r2s41iwuQArpd0+4L3fwXwZknrJf1U0uvSVjK0UdGIWAOsAXjticsWfWZtv1/6sEPNgVZ/HjltUiOY3Jl5VHQUEya/KSIekXQwzRmt7ktajEuAA4A3AK8Dvivpn3WbnxQqcEpVL8MMNQfa+BlmwPV7ylWWy4ZXaVKXYUyYHBGPJP9vkbQWOBm4CdgMXJ287lZJDeAgYLrb+kp1uEc/tbVhhZoHB2xY28CYNEn71jZh8rt7TJi8QtLerdvAO4HWkRjfB96WPPcKYA9gsAmTJV0B/DHNNvZm4PMRcXGWP6js6hpo03P75L6OlUu25b6OUZuZ32vg2tuwT5avUq2thywTJh8CrE2eXwJcHhHXJq+/BLgkOeRsFvhgr2Zo6w16ioj3L/av6ceoa2tVDrVRBNcgZahy6I063MZhJquMEyY/BJzYZblZ4Jx+1lm5PrZBQ61qgVaGEOtXpzJXKeyGFW6AL3VUkJEFW69DPbLW1uoealUMsayqFnajHDktrNYWMLlzLn25CqpMjW2QUCtroNU5yLJY+PeXMegGDThfoLIYhQdb3iNDZQu1cQ+zXsocdIM0T7OE2zj0tY1S4cGWxWJqa2UJNAfZ4pUt6Hxgb3WMJNhGfdG8okPNYZaP9s+1yJBbTO3NTdLRKrTGltYMrVJNzWE2WkWHXB7h1q05WpNj2UaqtE3RfkOtiEBzmJVDUSFX9XBTo8HEM7NDfc+yyD3YujVDu9XWyl5Lc5iV26hDbjH9bm6W5q+QGtuwQm0UgeYgq65Rhly/tbde4eYR0sHlGmxZBw3KVktzmNVP6zvNM+Dat8ksIedwy8/Ia2wLa2tlqqU50OpvVIeQZG2i9tMs9SBCdiMNtjKGmsNsvOXdXM3SRO0Wbq61LV5uwZbWDM0aannV0BxotlBezdVBwi1XjUDPelR0IH1f8tuBZgXJoxY3jCuGWHa5BNs8u18Drp8mqAPNymSYIZcWbj4MZHhKc2nwPC7RPT23z/P/zAY1jO1p0G181KcnDkMfEyZ/KpkQ+W5JV0jaM3n87ZI2JBMp3yyp44Ur2+UebGm1tWEGWvuG5zCzPA2ynfXa3mcaU7vtI9Pzk3WYGyF1wmRJhwN/AayKiOOBSeDs5On/AvxZRJwEXA78m7QVlvaUqiwcXlYGi2mu9tssbYXbysn5ytXaIuL6tru3AP+qy6JLgOWSngOmgEdbbwG0PuR92x7vqtBg67em5iDb3da5FUUX4Xn7L3m66CKUQj/HyXXb/luBV7LLix8k6ba2+2uSuYT7dS7wnYUPJvOJfhX4DfAscH1bIH4E+JGkZ4FtNOcX7WmkwbbYq+COU6CVKaz60U+5xykEp+f26XvQYWFtLrdBhUbAszuzLp3rhMmS9gfOAo4GngK+J+mciPgW8CngjIhYL+kvga/RDLuucg22Xn0DWWtrdQq1qobWsKV9DnULvizHxy1cpmoXtRzChMnvAH4VEdPJ8lcDb5R0HXBiRKxPlvsOcG2H1++m1H1sVQw1h9fg6hp8C2tvnbbvhcu0am+tWlsVz0ZomzD5rd0mTKbZBH2DpCmaTdG3A7cBW4F9Jb0iIu4HTgM2pa1zZMHWTzO07IHm8CpWlYMvy7Y9ihP2Ryx1wuSkmXklsIFmc/UXNPvx5iR9FLhKUoNm0J2btsLcgm2xzdAiQ82BVQ/dvseyBl57ectaxkFkmTA5uf954PMdllsLrO1nnaVpimYJtF7B02mDcFBZu2FuD4MEUK9ybJ1b0fW925uj1ttIgi2tGZoWalk2SIeYjVLW7a0VUmnLtwda63arvy2380wbDeLpUhxKMnS5BNtc25hHpzMNWgatpZmVXZZAa79dx6ZoEfIJNia4/7mDd3usPcR6fYEOMhsHnbbzrKG22ONBx0kuwbazsfT5IFv4i9QuazXdrE4G2d5nGlP87OljaZ5yad3kEmy7YikP7XyhxjYzu/svzIF7NNv1aZ2oUM9RIhtP/fbLLTQzvxf37zyUDb9/6TCLVUv5NEUbE8zMTrHtueXPP7Ztdhn77LH7ybvtX2C3L939DlYHi6mlLTxQd/2OY/jV0wfywMxBwyxaLeUSbLONSTY/3bzk0vZde+72XHu4Zf2yHW5WZcPoalm/4xju2HoYj2/fm21PDOcyX9Fo0PCoaHaNmODx7Xvv9tjUHs8N9J6DdLaajdqw+o2n5/bhoZ0Hc8fWw/j1zAHMblvGkumlQ3nvOsunKTo/wc7ZpTy3cwlL95zrulyrr61bH1yauh+xbdUy7EGw6bl9+Pm2o9j89H7NUJtezpJtk+w5PdTV1FKmYEtOYv0GzataXhQRF/ZaPubF7Lbm/IfPQddwWxho7Y9nDbeWqp1GY9U0qhH8rXMrePCZlWx68iU8uWOK2enlLJueZOk2mJrudHEMa5cabJImgf9E86z6zcDPJf0gIu7t+qIQ7JqAZY0XPdXe57bP0me7vsViwq2TboeYmGVRxKFIW+dWsOH3L+WR7fvy+BP7ENuXsmx6kj2fgGW/D6Z+N1i3zjjIUmM7GXgwIh4CkPRtmheE6xFsMDErGkwQNGttLe19ba1R014BN2xuvlqaIo+r3Dq34vmRz+3blqOZPdjj9xPs+QRMbWmwbOscez66vbDyVUWWYDsc+G3b/c3A6xcuJGk1sDq5u+tX533m7sGLNxIHAU8UXYg+VKm8VSorVKu8rxz0DbbHk9ete+6KrMeOVOVzAYY4eJBc/3wNgKTbel1GuEyqVFaoVnmrVFaoVnkXzD+wKBFx+jDKUkZZpt97BDiy7f4RyWNmZqWUJdh+Dhwr6WhJe9Cc6+8H+RbLzGzxUpuiyaV5/xy4jubhHpdExD0pL1vMtFxFqVJZoVrlrVJZoVrlrVJZR06dJ4wxM6uuLE1RM7NKcbCZWe0MNdgknS7pl5IelHTBMN972CRdImmLpNIfbyfpSEk3SrpX0j2Sziu6TL1I2lPSrZLuSMr7haLLlEbSpKRfSPph0WVJI+lhSXdJ2jiMwz7qaGh9bMmpV60JTTfTHE19f89Trwok6S3ADuCbEXF80eXpRdKhwKERsUHS3sDtwHtK/NkKWBEROyQtBW4GzouIWwouWleS/jWwCtgnIs4sujy9SHoYWBURlTpodpSGWWN7/tSriJgFWqdelVJE3AQ8WXQ5soiIxyJiQ3J7O82ZsA8vtlTdRVNrWqWlyb/SjlJJOgJ4F3BR0WWx4RhmsHU69aq0O19VSToKeA2wvtiS9JY07TYCW4B1EVHm8n4d+Czw4qs2lFMA10u6PTmV0Rbw4EGFSNoLuAo4PyK2pS1fpIiYj4iTaJ6pcrKkUjb3JZ0JbImI24suSx/eFBGvBf458ImkW8XaDDPYfOpVjpK+qquAyyLi6qLLk1VEPAXcCJT1vMRTgHcn/VbfBk6V9K1ii9RbRDyS/L8FWEuzG8jaDDPYfOpVTpLO+IuBTRHxtaLLk0bSSkn7JbeX0xxQuq/YUnUWEZ+LiCMi4iia2+xPIuKcgovVlaQVyQASklYA7wRKP7I/akMLtoiYA1qnXm0Cvpvh1KvCSLoC+BnwSkmbJX246DL1cArwAZq1iY3JvzOKLlQPhwI3SrqT5g/euogo/WEUFXEIcLOkO4BbgWsi4tqCy1Q6PqXKzGrHgwdmVjsONjOrHQebmdWOg83MasfBZma142Azs9pxsJlZ7fx/VuFb77A4LSkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# --------------------\n",
    "# Post-process\n",
    "# --------------------\n",
    "# Plot solution\n",
    "fe.plot(sol.sub(0), mode=\"displacement\")\n",
    "plt.show()\n",
    "plot = fe.plot(sol.sub(1))\n",
    "plt.colorbar(plot)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subspaces"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, the solution is saved in *sol* variable as a mixed function.\n",
    "To acces individual subfunctions, FEniCS provides the function *split*.\n",
    "There exist two different versions of this function,\n",
    "the first is called via FEniCS package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "u, p = fe.split(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which returns some kind of link on subspace.\n",
    "Usually it is sufficient e.g. for plotting for example,\n",
    "but cannot be used for some operations with subfunctions.\n",
    "FEniCS therefore provides another split function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "u, p = sol.split()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which returns two subfunctions of a type *Function*. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complete code\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASYAAAD4CAYAAABBh0sxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2dXahk2XXf/+tU1b39OR94Wspgadx5SAzCD5ZpnAeZYAQOiiycPEZg58UwhNggSxG29WRMHhK9yA7EEAZLGGNjEyKZBMdKLMgYMRBL6pmMFUljByWWsCaadM903+77VVXnY+WhqmburVt1zv6fW/uctavWD5qenrvOrn1P1fnV2vusvY+oKhzHcSyR9d0Bx3GcZVxMjuOYw8XkOI45XEyO45jDxeQ4jjmGMRp95pln9Pbt2zGadhxnS3j55ZffVNVbq34WRUy3b9/G3bt3YzTtOM6WICLfXfczH8o5jmMOF5PjOOZwMTmOYw4Xk+M45nAxOY5jDheT4zjmcDE5jmMOF5PjOOaIUmDJUL3xd6n4x9UpFX9QFVT8o4o7JQfVVSoeAB5XV6j4k2qfih/riIuvuPgy8vfZABV/jHDH7EtOxV/PJlT8zWxMxT9Bxt/MplT8Uxm379rTGfeZA4Crz/41fcw6PGNyHMccLqaOOaiuUfGHZEZ2SGZjhyUZT7Z/RLbPxp9U+3xGSWaIEzYD1T0ynm2fi++CXMuNttf7UI4h9WGcS2mz8bGFBKQvpVwXuUe4OEpww75NSwlIREwupKZ4W0IC0pcSKyTAspTCsCCkBebF5FJqirclpdhCAuxJKbaQ2hyTspQA42JyKTXF75aUWCEB6UupjcRSlxJgWEwupaZ4l1IdPnRrhhXS7DXiSwkwKCYXUlO8LSEBuyclVkizY9KWUldCWmBKTC6lpnhbUvL5pJB4n09qgxkxuZSa4ttL6c3iZm0sezFXKlT8vfwJKv7N/AYVDwBvCnfMk1dPzElpV+eTVtG7mF4vD6n4B+Wiy2Fdjy0YADgoSSlFzmQmLep1Uua03MNoGP7FNtER/nryLuo1DgZx32M2MzymJXlMxZ+QS3AAgPv6qad3MVmClRIrpEV86FqwR/P4LHAd2FExuxhGWdi33KNiJu39LDz7PC72cXUQvk7rqNzHtUH4urSTcoQbg/CL4rTkLlA2S5rMM+1rRJ9yHVCvwcaPdYSbCBdxSWZVU7I/AFCCy6KbSEZMBxV3smJnSm2lFMojMn4hpeD2C+78HBfcN/pRyd/aPyk5aXQlpVBYwbQ5hh3uxZbSpoW0IAkxuZTqSV1KsYUE2JNSG4ntipSABMTkUqrHpVRPm/oka1KKLSTAlpQA42LaJSnFFhLgUmqMJ4UEpC8lC/NJqzAppl0SErB7WRKQvpR8Piku5sTkUqondSn5fFIzbeqZtklKgDExuZTqcSnV4/NJYViXEkCISUQGAO4CeF1VP7LpjqQupe9OngmOfVxyfecvaC5+Sl6cxaS+knyZt6ho4KG0KHLNyc/DgItnJXw0uFyB5btHjxqP2Zb5pFUwn8iPAXgNmy3wxN+Q3+pvldzSgwdkPJ31FFw8K6XH+RUMAwsmAeCQjD/Kr2BvEF5geVLsYY9of1wOcYVovw3jcoQbw/ACyGk1xD7Rp7waAMQ1PamGVJHoKv5f/mTtz9llRLEfaFFB8D7qiHqClCsi7wHw0wB+e4OvTWNuKNaBlBgOyfgjMv6kIJdBlFwmNi6HLY6Jmx3mZCbfZmgYe3jIxrNPwakiZFWhZ/E3AfwygLU5vIg8D+B5AHjuuecu37MlXEr1bIOUuPi4QgK6kVLM+DbHWJASEJAxichHANxT1Zfr4lT1BVW9o6p3bt26tbEOAi6lJlxK9bTJklKXEpuJlcjMSAkIy5g+AOBnROTDAK4AeEJEfk9VfzZar85gSUoupHpYwbQ5JvWhG9CNlBgsCWlB41lV1U8B+BQAiMhPAvhkF1KyJCTApdRE7CxpdkzaUupiKLYNUgKM1TEtcCnVs2tS8vmkOMdYlRJAiklV/wzAn0XpyRyXUj0upXpiZ0mAPSnFFhLQrZQAYxlTylJyIW3+GGtS8vmk7jAhpk0Iad0e2g+L63R/2GNoyZBLP2JLYzrmzn+bZRAMjyb8zgmPM7IyO+fO6RH5nrHLadj4EzKef8BBvfD+1vCAao+ldzG9XjxNxd8vNlp4fgF2fVlsKR3l+8gkfAP5o3wPwyxsK14AOM73grfiBYBxMcJoEB4/LYbYG4ZXWU+LIa4Mw7fincwlfC0L3+43rwZUNfq0GuAKwuOLihN3EVn0MXijeOrcvzf9ZdW7mBjYRa+PyOGVRSlx8dy36DEZPy7IoVVBZm5k/ITMDNvMJ03JY2JLiX8qSov5JPY1Iog1GTHtkpRSFxKQvpRYIQHpS4kVEhBvWJ+EmGJKyZKQgPSlFFtIgD0psUIC0pdS7HlG02LapSwJcCk1wQoJsCelNvNJuyYlwLCYXEpN8S6lOnw+KQyLUgKMismlVBcbV0jA7knJ55MCX6PDu4fmxGRJSpaENIvfrSwJsCclH7p1gxkxWRISsDtSOpnuoSjJDyo92cv9rsfCZ3kn0/NSffLquPEYa1JihQRsp5QAA2L63vQHqPgHZFX2Wzlf+f2I3D/68ZRdLkJWBZPx05IfmqRMXgywv3e+qPTR6fr3pCgHON0j91Ens8nYy2/Y5TFFB4uIN0l6JaeR4dekxc18WCmNyeHSeEoOx3JyaFWQF0Tk+IKUdptN49iMks3EKuXWr7WZT+qb3jMmBnb49ph83BEjJRdSPaww2hzThZQY+CGuPSmVZPuxSEJMloQEuJSasJYlAfak1KYoc1ekBCQgJpdSPS6lelghAfakxAppdky6UgKMi8mlVI9LqZ7YWRJgT0rt6pNsSQkwLKaUpeRC2vwx1qTk80lxMSmmmFLyLKme1LMkwJ6UfD6Jx5SYUs6SZvEupT7jfT4pDOtSAgyJadek9OAkvNJ9MuWK74qC+6COycpvrbgPdq6k8MYtdhEgJZzncSW5KHJ96spp8DG7Op+0it7F9P3pk1T8Qc4tXXk45eIB4ICs5D4k9qg+JbMSVkp5PoAQW/GW0wFkEB5f5QPIIHzr3i7QUgDmdygFA+JtqEgRn+VgHPaFy2ZusZfHFC1uBGyS3sVkicdklnRCD93iS4mhnHLxFdm+zjO3UJHpPHOTYbj4tCSzDDaelFIbiSmdKcWNt0AyYjoqOWnQ691cSrW0lVJwPDucJAUD2JMSKyRgN6QEJCIml1I9LqV6WCEB9qTUblI8TSkBCYhpl6SUupCA9KXUZihmTUopC2mBaTFZkpJnSfXEFtLsGFtS8vmkeJgVU0wpWcqSAJdSc7zPJ8WIt4w5MVnKkgCXUhPWpLSLQ7e2x1jGlJhcSutJXUiASymo/R3Oks5iRkyWpdS0FORkTG59S1Ypl5O40pCCvNjI4jsFGS8t1sdNwwssAWBKnqOKjKeXrpDSfnr/dGulBBgQ0/3pDSr+gKzkfhhYebvgJN+jyvbHxoZimmcA8XmVPAOISnHkGVVlLYVAh2T8iOhPOf+bOE1SCsC8bWyWdIlK8VDePOX2sq+YDwX63763dzExHNHr3bhMhh26mZQSgZDxYNsnMzE2/m0phbbfYhI9tpTaSKzNRDpD31ICEhKTS2k90YUE7J6U2ggjtsQ6GLpZkBKQiJhSlpJnSZuNB2BOSq2yHpdSLY1iEpErAL6M2VMLhwD+g6r+WrQeLWFJSpayJGAHpUQKCXAphWBJSAtCMqYJgA+q6pGIjAC8JCJfVNU/j9w3SkqWsiTApbTpeJ9PioNFKQEBYlJVBXA0/+do/oe7N0tiKUsCXEqN7acuJZ9PCozv7sGZQa8kIgMReRXAPQBfUtWvrIh5XkTuisjd+/fvt+7QLkmpnA7oSW5GSpJnnJTyjJKSFEJJho0HYE5KWolLqQOCJr9VtQTwoyLyFIA/EpEfUdVvLMW8AOAFALhz506rjCplKaWeJQ2OiP7QGQYZP2mRwUxnf1dXid01fegWEN/P48Wpu3KqeiAiLwL4EIBvNMUzxJJS6NamZzmdcFnS9JSLr8hKbplk1KXNZiVZm/IBQ0gp0HnRZ3Ya9rtoTl6g5Dmdsg8hYBcRk4KpE96Ta/Yl70tKQNhduVsA8rmUrgL4KQCf3lQHHky4CtaHE140DBN2U/sJF09nPlMuPssFzGc2mwpZKf6OBEL7UxGV3FkuqPaISnFy6CbzhIpK6dn8n82Uos7YNvOoxZf3MpteHhPyqX8WwIsi8nUAX8NsjumPN9WB//gT/zY4ln48Ejk/ZElKMs0oKWW5ICOygGwqMymF9icXCNM+2Z/FMQxtpRSMYuultAlirNkLuSv3dQDv3/grk1iSksUsiYonhASAEhLQoj9kPNCRlBjarI9zKa0lkcrveFKylCUBLqUmogsJ8CwpgNg7G5gXEyOllIdugEupCZeSDbrYbsWsmCwN3QBbUkpdSG2OMTd0A1xKETEpppSl5FnSZuPbLC3x+aQ4dLkxnTkxNUlplJXnHqfsk9w18TsmJR+6xaGPnTJNicknuc/8/MxFnJGV0ENWSAUVDiEvzsGYjCf7DwBKnqOSLbBcI8naSnOXUmvMiGmVlPI1e0s/HnNV4vSe3CdcFoYxV8mdjYn5JPIiZS/qLAdXkJkDym/JHRUpAGXeYlIYdZnYukrziszeSrJP4y3e7xswIiammjv2fFLBbvyfE8tF2KFbB1LqIr4KfAuyeeZWEZJhs71NSmlTr7ENmdWmMSGmUCxJiR267aqUguPZ4SQrJCC+lNoIxqW0kmTEFFNKbbIkCpdSfXxsKbW4+KNLyYVUSxJiijnJnbKUrAmp1WsYk5IP3WxgWkyWhm5AXCl5ltSMOSn50C0aZsWUtJR86FYf7/NJTgMmxWRJSpaGboBLqRGfT9oKzInJJ7nXw0jJmpAAe1Ly+SS7mBLTtkxyZw3bsA5OSSFNuK5kbDwpGSEfEMDWY7K/LwBUsTM3cnlMGbi1bsnsUb5DmBDT0WSfio9eyX06oPbYzk7DL71synWFzmTI9gdTsvK74ONDCywBICsBZrvstyVJvAYr1laZVSCDeeU4+yxPNn4s5DXQMybExDBmM6UWTyMJve7oRbgGpUS1T2YZdDz9qKa48UBcKW0TmWx2zJqUmGJKKfrOAC6l+niXUpJsWkgLkhETI6U2WRJDTCm5kJqJLSUXUhixpAQkIqZUpeRZUkN8F1mPSykKMaUEGBdTF/NJDC6lDcb70C1ZYksJMCwmS1Ly+aQNxxuTkgspjC6EtMCkmHySezNtpy4kwKVkhS6lBBgU065IyWKWNByHxbYZJjFoiwLLLAcK8knXLqUwupYSYExMlie56yp/ByfkntynVDgGbHxkiVlDSgBD7rzSciUlJvTWtzYrxfuQEmBETMdkJXc+Jrs9tlPJTWcyZPYwIDMxuj+5QjNmzZ6iGhHxhaLcC4/vYujGXpsxr+XWleKkJyc9CWmBCTExFORQT4nMSiZx55NiSqkLIcWMB2ZSYogtpTbXZs/X82oSfG5BUmLaFSlZzJKixpNCAuxJyaSQgCSlBCQkJkZKjJCAuFLahqFb1HhjWRKwJVJKVEgLyL09+sGltKZ9l1J9vEspWUxnTDGHbgAnJUvzSQAnJWtCAuxJyeeTbGFWTC6lNW17ltR8jM8nJY9JMfkk95q2XUr18T502xrMiSlVKXUxdAutFh+MuauHrxInDcD2Z8Jf/YOcuzrZrXKLqy6lLmm0gIi8F8DvAng3Zlurv6Cq/yZGZ3qb5F76wLF7crOV3MNjMj5wqciukhVARewc22aoNzrhjuEhDSOrP8/bsod4iAkKAP9CVV8RkZsAXhaRL6nqtzbVCbaSuxpzUsrGxHzShHw8EjvUIyXDZDNspjGcZ22hqycW7Wvg6V8M9UIrvxfx5T5TKQ4UV4lK8Q4mxft8MsqiMnwZtvK76Dk9bLxiVfX7qvrK/L8PAbwG4Adjd2wdVR5xktuQlAbTbqQUq/3481UtnnQSW0oKk49ropftGRgeUpMuInIbwPsBfGXFz54Xkbsicvf+/fub6d0SqUppMOalxLCLUmKQqiMpGUMlTSkBhJhE5AaAzwP4JVV9vPxzVX1BVe+o6p1bt25tso8AOCnJJDMlJQaXUkN8B1nStkiJQmBGSkDgXTkRGWEmpd9X1S/E7dJFWCkxMFKyNJ8EcNJIXUiA0aGbQVLNks4ScldOAHwWwGuq+pn4XTqPS2lFrGdJjeyilGghASalBIQN5T4A4OcAfFBEXp3/+XDkfqHKB9GklE3EjJR8krshPvJ80s4O3QCzUgICMiZVfQkd/wqpTnIDPsm90XhrQzdgO6RkWEgLzFV+pyoldujGFOwNT0ghnbKV39wVTc8PkcWJgyn5PhVsFXeGvWPyd2iRkhTX6ENotlFKgDExxZ5PkirsXWH32GYruWNWEQ+mcaU0mCiUOPVZrqiGxLC5VJREvBQKEFv3tkIVba7oIfE+s4JR4Q6YXkVUKX3nn/7qRtszIaZqwnWD2ZMbWF8NuzKWHf4YmuTuQkoM9NCtJBf5souC26zW0Mhjty4yGOFWPwBAJdw1tmlMiImCfqRS+DvvUqqJ3zUpnRFSfi3SfoodSSn6MRF+j7TElKiUYk9yu5Qa4i8hpWhYlJKhO3vpiImQEiMkYHektGtCAvqRUuOk9zZIKfLvkIaYjEjJUiW3Z0kN8Ts8n3TZY6qrDTv0dfB72BeTS2lF2y6l2vjEhm75psoKEh66LWNbTDsgJZ9PqmfbpdSGlUPFmtcorq/4HQ1LCbAspkhS2pX5pFl8PClZExJgUEqJDN2ivMYlsSemLbjz1rSJ/oisOh4dk4I55eKHY27XfyH3/B6SwqAr9EtytUBFSpL2V/hnuLjOtn2GLZUSYE1Mm5bS/IKg9+QmK7PZ9mPCZlWslLJpRVZ+V9BB+AFSVsCQWAHASoaMjw27agAA/ZhazS5eJ+W1gPPQ4/IVG2IiKrMBIBvHO2P0fBKZhQ2JJ4bwTzuJLyUqPufipSTjE5dSG1rtIrCCwUlNQ2sedNAl/fcAwHf++SejtT0ghgUupfW4lPpnU1KqpcVrxOiXjYwpkCznPGpFSoyQZn2xIyUXkg0sSilmn5IRk0tpTTwhJc+S0sOikID4/UpCTFakFHPoNutLvCJLl1J67KqUACNzTHWskpKueerqYCK9SKlaofc6Ka2Kr5NSNTj/OzVKaSnepZQeFqXU6nFQLTGdMTGZUp2QquHFbVrrpKRDQM7E+yR3TTwhJRdSM11d+Jbmk1ZhVkybktLKeCPDN6tSyib1m21LzrXLVlgLUc9WXt+bHeNSCsPo0G0Zk2KKJaVNTnJnS9dm7Gru4TEng+FJ+E7+2biEdLFmLAKD4ym/fIVeukLGs6ypG8ovUxW+8nX4Q/qQEmBQTJuQ0qo1UzH35WYzH/phAcQSkwGbJbHLUaYFwOw3nZfAkJjKLCoqXtg1eGz21qOwR+s+s+wwrOZ0rts7qi8hLTAjpgG5J/Hw1EY5wM5JiYEd8pFr8LZZSl2xavmVGqj8NiMmBmaRZ6pSYhfiJi0lF5KzRHJiiiUlfueB3ZCSZ0lOHyQlphSlFHPoNuuLS2ltvEspWZIRk0tpVV+INW+Whm6AS2mL+PavfHzjbSYhJpfSqr70KCXV83fmIkrJhbSbmBeTS2m5Hz7JvTbepXRpNl471RLTYrIgpW0qB6ir2M7G5IblOSEw8pwI0zaAQdV8Dqvr++dfw4iUpje4W/P5Da59Nt4KZsVkSUrL6+xWsXfIZQMjtpL7KKfisxMunoIUhwWy43fe+GFEKeU3hhgdtXmoXWj7/dcYdYE5MbEb0Z+VktRcL2uraNfFH8X78A7Gke+8Nax1Ox9LCoyV0jQHhsTHrCoBBO75zQqGjbe2dAXA6IiLD+lSYTCrMiOm//2JT+CH/+VvBMcPyTeIIWqN0q5JKZQq7sLgtlIqb6zZY+dC+1zzljh7LfW9FGWBGTHFZGhkTimmlBghzeJdSmsxmCntGo0DVhH5nIjcE5FvdNGhJgaT2Z91m8Uts5BSFRB/VkqrNnM71+6pnpNSNaz/qtkZKU1zl5JzaUJm0n4HwIci9yMIem+kSJkSfedtl6TEYElKFTgpKVxKEWkcyqnql0Xkdvyu1LNtUqoGcuGx2THnk2bxhDhiDt2AuFKqiV0uG5j9T64rLqT4bOzeo4g8LyJ3ReTu/fv3N9UsAPKpJBP7Ulod61JaierGpLS6L1y4S6kbNiYmVX1BVe+o6p1bt25tqllaSlTbOyClbJLbkVJV8lJicCltDabvyrmU3kGHGaSoNl8OUMxff0yewJwQ0pSsKi/IZTRLDznQVcO1s7iQzGNWTH1LaTAF9g65CyR+NXf4BS6npAwYplNue92OkeP1H4iMlBK9dIVdr1eSl+CZ5qc3+68C/8tf3/zOAkCAmETkDwD8JIBnROR7AH5NVT8bpTdz2kgpC7jG9x7H++obsrsDEA8LAAAZh8fLmJzzYTIaMvvRPIcwld9lFf51aWS929vQi4i55qdPnK+Ib1oGNb2RYe8wsO2bXF9iE3JX7qNddGTBYAJI4HW1R1R/D0/Ip5ichIvGpbQaZYZ7wExKwY2nL6X8BiFs8tdd9UCOOvYOgfwmgkUWG1NDOWaZyZBYuJuqlBghzeITlRL5IExTUoqcJbWZ32KlBCyGrHaG5/0PUs/wzU+HjVddSqvil0RQNwc0nbqUNsFSX4obe7XhFqUkqma2gDmLKTGFcFZKjctGlqTUFL8sJR2sv7hNS6kO9g6ZdSntN9yBi4WxTEmqdlKyiqmhXBPWMyXNZOWjqndBSlHnkwDTmVIdVrMk6yQjJutSWodLaQWpTnIbzJJYUpASkMhQzqW0Kt6lFBu9dmbOyKXUKeYzphSlFCqkxdDPpbSqcUMXUUwpRRq6TZ54J+dISUgLTIupDyktLycZHZEb45N7bQtVzU2WuE/C45WIBUCtqVN2kW9JPJQz4EEE5+PjtT1kn+pCPDWmuDHClQerz/nkqfWXcaiUpk/wpQJf/0ycqm/AuJhCufIg/A3ee0xut0FAS4l4MolMIt5NY8XBSKkoIIPAPbwBaFWFV9NE3Z+JbJuczK+uBe50GMD+wcX3IyvCzzkASJVh8pSdOiazYhqOw5am7D+K90SK4XH4BehSuogWXLapTIaSsJT4LI/NxNoN3fYPwo+LLTFzk9/f/PTHsf8oTEqjY23c0vYse0cVkIXFL6QUsjm7S+kiLqWLSFVRUpJKo0upTf0TAFx5EHfeymzGtGn2iGd9mc2UMqm/cC4hJZEMqjXnKJKUKCEBkTeNI+K3IEtqdWcv3gDlHMmKaXQc/kZshZSaSDBT2iYpVdevrA23JqW2culKSkCiYnIpLeFSCou9smb5ipFMaVNCmtbdpVvRnfHT9TM6XQppQXJiciktEUtK5Ba75qW0isjzSQx1UsqfuLg4eBuzpLMkJSaX0hLbLqWEJ7kZLA7dmo756u98gm+UwNxduXXUSWn5zlyjlM7cmetdSvN3wKW0HGxESmVlRkpSqAkpdYHJjOmVf/dx/Ng/+423/91HpjQ8On/hM/ttg8iSAHCSIZaBsBXXyvSDqM6u2CevMPM4zBwOu6MB8TtmzPKcvT0MiWVF2enq2OkPXFt7TKpCWmBSTGcJkdLouML+A04Gy+LpjcnExMb+9Lq3WGgFgKtajtKNsoRgc9XZMdh76+TC/yuv7+HqG+FfBNMn93DtjZmRTt9lRwd2erKChZRGx5tV+fDRBKjZBG6Z7HAMDAJHvSenQBZ4YbHr0ybT4AJRJlviNnfjlvRoWQISeO7q6qguCzHU07KEXF1/+/8cb5+PQIkVJVC/0eWKtsOQaQFcD238IlfvNQutK3mZFdO1e+Fvyv7D8Atr+IiTQXZIrCQ+OZ39XZXNcmojpUBMSSk42I6UgiHPB/W8vHnbeuNqULhMySczT9utGb3+f7vJrM2KKZS3pSTS+AE8J6VSG7OmVlIKwaW0FBxJSuQEejQpkQ/wbJUpEbSVUjafeH/xT3+l1fEMSYvJXKYUgktpKXgLpFQ37IsoJVZIQDspZS0XBV8Gs+UCL33hk7U/71xK6+7odC2lNbeXXUpn2+WGbiYypbK8tJSKp9ffpcum5QUpjW/VDxOzQnuREmBYTHUwUmKJlimxeKbUjljzScBmpHRthQwMDt36EtKCpIdyF1gxz8RkS1GlxGRLiUmJvsBrpCTDS9yiT3iSOwQLQupifglIUEwbHcKdmQDvVUqq79Qy7bCULoUFKW14PunsHbnYUuo7Q1rG9FBueZ4p1rySZ0pniCUlrVxKZ9uNlCmtmktaxWJ+qc95pDqSyZjqpLT3xuN3/nF0sRq2FuqpIOzSikgXObG0Ita8D5cpxdtnnemHUA9QIBZfM19Wx+GfTzk8Cm93NMLwcXj8tccXJ8rHP/R0+OtFxryYbr72FhWvR8eQG9eDYquHB8iur7+TcaHtyQQS+EhqPTkJjq1Ox5C9sLkVnUyC52G0yIOrrrUsIaFV5ezwLRahFeUAtdRFixwyCHyfpzlk3T5PF2KnkGHYJaeTCWTU7eV55bsPa3/+xb/61x31xPhQjkXfuB+t7eroOLwfJ+HfitVp+DCSecSSFsRQL+acUiSU2SFAK07moc0yw2kmMyfe5+r4BBWRhenxCUB8Gb/N40P+mEtgXkxffO1fRWm3engw+zvgTXUp8bExYaUUHDo/byGZrhUpMSgZD2AmpI6lBCQwlAuFyZYWUgqKdSnRsbGIJSSgfaaUNQzjzkpJrtVnKmff56xhOiK6lJZk9MXv/xZ3/CUxnzEBzVnTspS0RiYupXmsS2kWWuTB502nuclMKXvqyfXtHp9cWkp9kHzG5JnSPNaltAgOD01sPgngMqVVQpJ3PbP+gDVC6jpbAgIzJhH5kIj8lYh8W0R+NXanVhFrrgkg32xGFkws88EnLqgLx7qUZqGJSYmZ5KazpJ7mkepoFJOIDAD8FoB/COB9AD4qIu+L3bEQ6rKl5eHczmZLLYsaXUrzWCNSCm53w8O2PrIlICxj+nEA347S8QgAAAO8SURBVFbV/6OqUwB/COAfxe3Was5mTT6Em8ducAi3EIFLaR67ZVI6N4wzmCWdJWSO6QcB/M2Zf38PwN9bDhKR5wE8DwDPPffcRjq3iphDOsdxbLCxu3Kq+oKq3lHVO7du3dpUs47j7CAhYnodwHvP/Ps98//nOI4ThRAxfQ3A3xGRvy0iewD+CYD/FLdbjuPsMo1zTKpaiMgvAvivmK2C/JyqfjN6zxzH2VmCCixV9U8A/EnkvjiO4wBIZEmK4zi7hYvJcRxzuJgcxzGHi8lxHHOIkk8sDWpU5D6A72684Xd4BsCbEdvfJCn1FfD+xiSlvgLx+/tDqrqyGjuKmGIjIndV9U7f/Qghpb4C3t+YpNRXoN/++lDOcRxzuJgcxzFHqmJ6oe8OEKTUV8D7G5OU+gr02N8k55gcx9luUs2YHMfZYlxMjuOYIykxWXgoQigi8jkRuSci3+i7LyGIyHtF5EUR+ZaIfFNEPtZ3n9YhIldE5Ksi8hfzvv56330KQUQGIvI/ROSP++5LEyLyHRH5nyLyqojc7fz1U5ljmj8U4X8B+CnMtvf9GoCPquq3eu3YGkTk7wM4AvC7qvojffenCRF5FsCzqvqKiNwE8DKAf2zx/IqIALiuqkciMgLwEoCPqeqf99y1WkTkEwDuAHhCVT/Sd3/qEJHvALijqr0UhKaUMZl5KEIIqvplAA/67kcoqvp9VX1l/t+HAF7DbL93c+iMo/k/R/M/pr9hReQ9AH4awG/33ZcUSElMqx6KYPLCSR0RuQ3g/QC+0m9P1jMfFr0K4B6AL6mq2b7O+U0Avwyg3fO0ukcB/KmIvDx/0EinpCQmpwNE5AaAzwP4JVV93Hd/1qGqpar+KGZ70P+4iJgdLovIRwDcU9WX++4LwU+o6o9h9jzJX5hPTXRGSmLyhyJEZj5f83kAv6+qX+i7PyGo6gGAFwF8qO++1PABAD8zn7f5QwAfFJHf67dL9ajq6/O/7wH4I8ymUjojJTH5QxEiMp9Q/iyA11T1M333pw4RuSUiT83/+ypmN0T+st9erUdVP6Wq71HV25h9bv+bqv5sz91ai4hcn98AgYhcB/APAHR6dzkZMalqAWDxUITXAPx7yw9FEJE/APDfAfywiHxPRH6+7z418AEAP4fZt/mr8z8f7rtTa3gWwIsi8nXMvrC+pKrmb8EnxLsBvCQifwHgqwD+s6r+ly47kEy5gOM4u0MyGZPjOLuDi8lxHHO4mBzHMYeLyXEcc7iYHMcxh4vJcRxzuJgcxzHH/wchqVA+w6CW+AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAEDCAYAAAC/Cyi3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAedElEQVR4nO3dfbBcdZ3n8ffn3oSQG54hIE8KK/jAUIBOREd8GFEsFllxp1YHa7C0UFNajgOujoMzu2u5tVqMuo5OzT5lgVksAR+AuK4okEVWhl0JQgyPQaAQNYDkcgkmAZLLvf3dP/o0dC7dfU7f7tPnoT+vqlT64XSf3+0+59O/h3POTxGBmVmdTBRdADOzYXOwmVntONjMrHYcbGZWOw42M6sdB5uZ1Y6DzazGJF0iaYukuzMs+3eSNib/7pf01CjKmAf5ODaz+pL0FmAH8M2IOL6P130SeE1EnJtb4XLkGptZjUXETcCT7Y9JermkayXdLumfJL2qw0vfD1wxkkLmYEnRBTCzkVsDfCwiHpD0euA/A6e2npT0MuBo4CcFlW9gDjazMSJpL+CNwPcktR5etmCxs4ErI2J+lGUbJgeb2XiZAJ6KiJN6LHM28IkRlScX7mMzGyMRsQ34laT3AqjpxNbzSX/b/sDPCiriUGQKNkkPS7orGQa+Le9CmdlwSLqCZki9UtJmSR8G/gz4sKQ7gHuAs9pecjbw7Rji4RKSDpC0TtIDyf/7d1jmZZI2JBlzj6SPJY9PSbpG0n3J4xdmWmeW8kt6GFgVEU/0+TeZ2ZiT9GXgyYi4UNIFwP4R8VcLltmDZh7tSvoB76bZF/gU8PqIuDFZ5gbgSxHx417rdFPUzPJ2FnBpcvtS4D0LF4iI2YjYldxdRpJNEfFMRNzYWgbYAByRtsKsgwcBXC8pgP8WEWsWLiBpNbAaYPmU/vDol4/HuMQSpS9jozE3Rsea33vXc09ExMpB3uO0ty2PmScbmZb9xZ2z9wA72x5a0ykHujgkIh5Lbv8OOKTTQpKOBK4BjgH+MiIeXfD8fsC/AL6RtsKs6fOmiHhE0sHAOkn3JQf+PS/5I9cA/MEJe8TlP+xY9tpZOVnZEfHamZ6fLLoII3PSyzb/etD3mHmywf+99rBMy04d9vDOiFjV7XlJ/xt4SYen/qb9TkREUkF6kYj4LXCCpMOA70u6MiIeT95/Cc0Dhv8+Ih5KK2+mYIuIR5L/t0haC5wM3NT7VeNhen7S4VYC4xRqZRQR7+j2nKTHJR0aEY9JOhTYkvJejybntr4ZuDJ5eA3wQER8PUt5UvvYJK2QtHfrNvBOmh17ZqXgUCu9HwAfTG5/EPifCxeQdISk5cnt/YE3Ab9M7v8HYF/g/KwrzDJ4cAhwczI0fCtwTURcm3UFZjb2LgROk/QA8I7kPpJWSbooWebVwPokZ34KfDUi7pJ0BM3m7HFA63CQj6StMLUpmrRnT0xbbpy5OVoc19bKLyJmgLd3ePw24CPJ7XXACR2W2Qz0PUTnwz2GxDvY6Pkzt27G45iMEXHNzapknmBrY1f6ghXkGtuQuRYxGv6crRcHWw680+XLn6+lcbDlxDtfPvy5WhYOthx5Jxwuf56WlYMtZ9Pzk94hh8CfofXDo6Ij4hHTxXGg5Wcu6vv5usY2QnXdiPLiz8sWy8E2Yt5Z07n5boNysBXAO253/lxsGNzHVqD2nXjc+98caDZMDraSGLfBBQeZ5cnBViILd/a6BZ3DrFzmmGCmMVV0MXLhYCuxOgSdw8yK4GCrkE4hUaawc4hZWTjYKi4tTPIIPgeY9UPSAcB3gKOAh4H3RcTWDsu9FLgIOJLmzHhnRMTDbc//PXBuROyVtk4HW805hKwELgBuaJsw+QLgrzos903gixGxLpk0+fm5ASWtAl40g3w3Po7NzPKWOmGypOOAJcklwomIHRHxTPLcJPAV4LNZV+hgM7MsDpJ0W9u/1X28NsuEya8AnpJ0taRfSPpKEmgAfw78oO09Urkpajam5mKSmfnU7qqWJ3KeMHkJzXlEXwP8hmaf3Ick/Rh4L/DHWQvaejMzs4EMYcLkzcDG1izvkr4PvIFmDe8Y4EFJAFOSHoyIY3qVx01RM8tb6oTJwM+B/SStTO6fCtwbEddExEsi4qiIOAp4Ji3UwMFmZvlLnTA5IuaBzwA3SLqL5lyi/32xK3RT1MxylWXC5OR+x0mTF7wmU6ega2xmVjuusZmNqbmYYHpun6KLkQvX2MysdhxsZlY7DjYzqx0Hm5nVjoPNzGrHo6JmY2ouJtk6t6LoYuTCNTYzqx0Hm5nVTuZgkzSZXCfph3kWyMxsUP3U2M4DNuVVEDOzYckUbJKOAN5Fc6IFM7NSyzoq+nWa1xvfu9sCyaWCVwMcergnEDEru3kmxndUVNKZwJaIuL3XchGxJiJWRcSq/Q/wmISZFSdLAp0CvFvSw8C3gVMlfSvXUpmZDSA12CLicxFxRHJZ3rOBn0TEObmXzMxqQdIBktZJeiD5v+P8oJL+VtLdyb8/bXtckr4o6X5JmyT9Rdo63WY0s7y1Jkw+Frghub8bSe8CXgucBLwe+Iyk1sXiPkRzdvhXRcSrabYce+or2CLi/0TEmf28xszGXuqEycBxwE0RMRcRTwN3Aqcnz30c+PcR0QCIiE6zXO3GNTYzyyLvCZPvAE6XNCXpIOBtNGtpAC8H/jRZ748lHZu2Qp8Ebzam5hoTzMxOZV081wmTI+J6Sa8D/h8wDfwMmE+eXgbsjIhVkv4EuITm5MpdOdjMbGBDmDCZiPgi8MXkNZcD9ydPbQauTm6vBf4xrTxuippZ3lInTE7ORT8wuX0CzWn4rk+e/j7NpinAW3kh8Lpyjc3M8nYh8F1JHwZ+DbwPmhMmAx+LiI8AS4F/kgSwDTgnIubaXn+ZpE8BO2ibi7QbB5tZhaycnE9fqGSyTJgcETtpjox2ev1TNM9VzyyXYFui3l/A9LzPJTVLU8UQK4tCamytL8wBZ9ZURIjNxwTbnls+8vWOQqFN0fYv0yFn48o1s+ErTR+bQ87GjQMtP6U83MNfuNWdt/F8lTLYwF+81Ze37fyVNtjAG4DVy8rJeW/TI1KaPrZuVk7Ou8/NKq+MgTYfYtvssqKLkYvSB5tZlZUx0MZBqZuiLa7CWxV5my1OJYLNzKwflQo2/wJaVXhbLValgg28wZhZukoOHnik1MqsKj++jZhg+649iy5GLipXY2upysZjZqNX2WAzKyP/4JZDpYPNG5GVibfHziS9V9I9khrJVXM7LXOkpBsl3Zsse17bcydJukXSxmSmqpPT1lnpYDOzSrgb+BPgph7LzAGfjojjgDcAn5DUuqLul4EvRMRJwL9L7vdUycEDs7Jxba27iNgEkMxn0G2Zx4DHktvbJW0CDgfuBQJozQq/L/Bo2jodbGZjqhHimdmlWRc/SNJtbffXRMSaHIqFpKOA1wDrk4fOB66T9FWarcw3pr1H5YPNh35Y0caktrboCZMj4kXT7fV4n72Aq4DzI2Jb8vDHgU9FxFWS3gdcDHSdxxRqEGxmVrxeEyZnJWkpzVC7LCKubnvqg0BrMOF7wEVp71WLwYMx+cU0qy01O+AuBjZFxNcWPP0ozYmSAU4FHkh7v1oEm5mVl6R/KWkz8EfANZKuSx4/TNKPksVOAT4AnJoc1rFR0hnJcx8F/qOkO4AvAavT1ummqJnlKiLWAms7PP4ocEZy+2ag47Bp8twf9rPO2tTY3Bw1sxbX2MwGUOUf1EZD7Mx+uEel1KbGZmbWkhpskvaUdKukO5JzuL4wioKZmS1WlqboLuDUiNiRHGdys6QfR8QtOZfNrNSq3Aytu9Rgi4gAdiR3lyb/Is9CmZkNIlMfm6RJSRuBLcC6iFjfYZnVySVFbpuZaQy7nGZmmWUaFY2IeeAkSfsBayUdHxF3L1hmDbAG4MQT9yikRufzRs2yixDP7azngRF9jYpGxFPAjcDp+RTHzGxwWUZFVyY1NSQtB04D7su7YGZmi5WlHnoocKmkSZpB+N2I+GG+xTIzW7wso6J30rzom5lZJfjMA7NF8DFs5VbPIREzS9eAmK3nUQSusZlZ7TjYzKx2HGxmlqssEyYny+0n6UpJ90naJOmPkscPkLRO0gPJ//unrdPBZmZ5yzJhMsA3gGsj4lXAicCm5PELgBsi4ljghuR+T7ULNo9WmZVLRGyKiF/2WkbSvsBbaE7oQkTMJmc6AZwFXJrcvhR4T9o6PSpqNq5CsCtz3SbvCZOPBqaBf5R0InA7cF5EPA0ckswUD/A74JC0N6tdsPkkeLNc5D1h8hLgtcAnI2K9pG/QbHL+2/aFIiIkpV5ko3bBZmajN4QJkzcDm9suiXYlL/SlPS7p0Ih4TNKhNC+f1lPt+tjMrHoi4nfAbyW9Mnno7cC9ye0f0JwNnuT/1Bqgg83McpVxwmSATwKXSboTOInm5MgAFwKnSXoAeEdyvyc3Rc0sV1kmTE7ubwRe1I8XETM0a3CZOdjMxlXAxGzHydcrz01RM6sdB5uZ1Y6Dzcxqx8FmZrXjYDOz2nGwmVnt+HAPs0WYnp+s/JVk1ICJnfWs29TzrzKzseYaW43NNKaG/p4HTjwz9Pc0GzYHW43kEWRp63DQWRk52CpuFGGWdf0OOSuLWgXbOFxksugg68UhZ2VRq2CrszIHWiet8jrgSixgclfRhciHg63kqhZoCzngrAg+3KOkZhpTlQ+1dq2/p05/U1H2n1hWdBH60se8oudJujtZ9vy2x7+SzDV6p6S1kvZLW2cuwTZJPa/xNArjsPOPw99ou0mdV1TS8cBHgZNpzil6pqRjkqfXAcdHxAnA/cDn0lboGltJjOPOXvW/uYjBqqrV1iDbvKLAq4H1EfFMRMwBP6UZhkTE9cljALcAR6StM7dgq+IXUISq79zD4M/AaNbq3izpQElTNC8ZfmSH5c4Ffpz2Zh48KJB35t3NNKY8yNDDsCsLCpiYzbx4zwmTB51XNCI2Sfpb4HrgaWAjsNvJuJL+BpgDLkt7v1yDbf+JZWxt1HQ8eQAOtO48ilpaPSdMHsK8okTExcDFAJK+RHOuUZL7HwLOBN4eEakTJtemj60qB+c61LKpyudUle2uCiQdnPz/Upr9a5cn908HPgu8OyIy/eKlBpukIyXdKOneZBj2vH4K6762F1RlZy0Lf14vqPJ+1Me8oldJuhf4X8AnIuKp5PF/APYG1knaKOm/pq0zS1N0Dvh0RGyQtDdwu6R1EXFv2gutyTvo4rnfrfr6mFf0zV1ef0ynx3tJrbFFxGMRsSG5vR3YBBze74rGlUNtcOM+alrl2lpR+ho8kHQU8BpgfYfnVgOrAY483P0O47wj5sW1tyFrwOTOoguRj8yDB5L2Aq4Czo+IbQufj4g1EbEqIlYddODuwTZuvzgOtfz4s7UsMgWbpKU0Q+2yiLg63yL1r0wjU97x8ufP2NJkGRUVzWNLNkXE1xa7onGotXmHG50y9buV6YfVmrLU2E4BPgCcmgy1bpR0RtqLxk1ZdrJxU/fPfRwqBHlIHTyIiJvBl+vope47V9l5UMEWqs2ZB0VxqJWDvwdrV/mT4Ivs3/DOVC5F1tzymEA572aoAiaznwRfKSOtsdWpv8ChVk7+XgzcFF0U7zzl5u/HHGx98k5TDVX/nurUuilCpYNt1P1rVd9Zxo2/r/FV6WAbJe8k1eTvbTyNfFS0ilfV9c5RbaMaLc1jZDRXjfpOmOwaWwqHmln1VDbYRtG/5lCrj1F9lz5v9MWyTHic5Urdkj4tKSQdlLbOygabWb/8Q1WYLBMet67UfRzwBuATko5rPSnpSOCdwG+yrNDB1oV3gnry9zp6WSY8znCl7r+jOaFL6gxV4GDryBt/vfn7XZSDJN3W9m/1It8ndcLjhVfqlnQW8EhE3JF1JZU8VzTPfgxv9OPBVwRpniu6ZGemChCkzCuaZcLkLBMeL7xSdzIr/F/TbIZmVslgy4tDzYZh0MM+qnjWQdqEyVkmPO5ype6XA0cDdzSvecsRwAZJJ0fE77qtz8GWcKiNH9faRqNtwuO3dpvwuNuVuiPiLuDgtuUeBlZFxBO91lm5YMujGVq3UJuZ3yvX9z9wckeu7z9KDreR+AdgGc0JjwFuiYiPSToMuCgizuCFK3XfJWlj8rq/jogfdXzHFJULtmGreqjlHWJZ11nlsMsj3Cp3FkKOuk143D5hctYrdUfEUVnWWalgG/eDH4sIsayqHnZlqrltbeyqZD9bmVQq2IatCrW1ModZmvayVyHkyhRuo6AGTO7KPCpaKWMbbGUOtSqHWTdVCblhhpubo8WpTLANsxlaxlCrY5h1U5WQs+qqTLDV1TgFWidlDLlxa5LW0didUlWW2trM/F5jH2oLlekzGdZ2Mu4DXkUZqxpbGUKtLDtumbU+o6JrcK65VVclgq0Ov3oOtP6VJeDqqs9zRStlbJqiRdXWytS8qqoiP8NhbDd1+GGumkrU2AZVRKg5zIavqBqcm6TVU/pgG/TXbtSh5kDLXxEB53Crllo3RUcZam5yjt6oP/NBtic3R0ertsE2qlBzoBWvKuFmo1PqYFvsr9woNj4HWrmM8vtY7PbVz/Zctbl3y2bkfWx5f2GjCjUrJx8i0odGMPlsPc9lLW2NrYx9Eq6lVUfe39Uoam22eKnBJukSSVsk3T2KAg0iz9qaA62ayhhu4ybLhMnJcg9LukvSRkm3LXjuk8l73CPpy2nrzFJj+x/A6Zn+giFZzK9aXhuZa2nV5++vcFkmTG55W0Sc1D4jlqS3AWcBJ0bEHwBfTVtharBFxE3Ak2nL1ZF3iPrI6wfKtbZ0WSZMTvFx4MKI2JW835a0Fwytj03S6tZkqk/MdO6QzGvgII+Ny6FWT2UIt4r2s41iwuQArpd0+4L3fwXwZknrJf1U0uvSVjK0UdGIWAOsAXjticsWfWZtv1/6sEPNgVZ/HjltUiOY3Jl5VHQUEya/KSIekXQwzRmt7ktajEuAA4A3AK8Dvivpn3WbnxQqcEpVL8MMNQfa+BlmwPV7ylWWy4ZXaVKXYUyYHBGPJP9vkbQWOBm4CdgMXJ287lZJDeAgYLrb+kp1uEc/tbVhhZoHB2xY28CYNEn71jZh8rt7TJi8QtLerdvAO4HWkRjfB96WPPcKYA9gsAmTJV0B/DHNNvZm4PMRcXGWP6js6hpo03P75L6OlUu25b6OUZuZ32vg2tuwT5avUq2thywTJh8CrE2eXwJcHhHXJq+/BLgkOeRsFvhgr2Zo6w16ioj3L/av6ceoa2tVDrVRBNcgZahy6I063MZhJquMEyY/BJzYZblZ4Jx+1lm5PrZBQ61qgVaGEOtXpzJXKeyGFW6AL3VUkJEFW69DPbLW1uoealUMsayqFnajHDktrNYWMLlzLn25CqpMjW2QUCtroNU5yLJY+PeXMegGDThfoLIYhQdb3iNDZQu1cQ+zXsocdIM0T7OE2zj0tY1S4cGWxWJqa2UJNAfZ4pUt6Hxgb3WMJNhGfdG8okPNYZaP9s+1yJBbTO3NTdLRKrTGltYMrVJNzWE2WkWHXB7h1q05WpNj2UaqtE3RfkOtiEBzmJVDUSFX9XBTo8HEM7NDfc+yyD3YujVDu9XWyl5Lc5iV26hDbjH9bm6W5q+QGtuwQm0UgeYgq65Rhly/tbde4eYR0sHlGmxZBw3KVktzmNVP6zvNM+Dat8ksIedwy8/Ia2wLa2tlqqU50OpvVIeQZG2i9tMs9SBCdiMNtjKGmsNsvOXdXM3SRO0Wbq61LV5uwZbWDM0aannV0BxotlBezdVBwi1XjUDPelR0IH1f8tuBZgXJoxY3jCuGWHa5BNs8u18Drp8mqAPNymSYIZcWbj4MZHhKc2nwPC7RPT23z/P/zAY1jO1p0G181KcnDkMfEyZ/KpkQ+W5JV0jaM3n87ZI2JBMp3yyp44Ur2+UebGm1tWEGWvuG5zCzPA2ynfXa3mcaU7vtI9Pzk3WYGyF1wmRJhwN/AayKiOOBSeDs5On/AvxZRJwEXA78m7QVlvaUqiwcXlYGi2mu9tssbYXbysn5ytXaIuL6tru3AP+qy6JLgOWSngOmgEdbbwG0PuR92x7vqtBg67em5iDb3da5FUUX4Xn7L3m66CKUQj/HyXXb/luBV7LLix8k6ba2+2uSuYT7dS7wnYUPJvOJfhX4DfAscH1bIH4E+JGkZ4FtNOcX7WmkwbbYq+COU6CVKaz60U+5xykEp+f26XvQYWFtLrdBhUbAszuzLp3rhMmS9gfOAo4GngK+J+mciPgW8CngjIhYL+kvga/RDLuucg22Xn0DWWtrdQq1qobWsKV9DnULvizHxy1cpmoXtRzChMnvAH4VEdPJ8lcDb5R0HXBiRKxPlvsOcG2H1++m1H1sVQw1h9fg6hp8C2tvnbbvhcu0am+tWlsVz0ZomzD5rd0mTKbZBH2DpCmaTdG3A7cBW4F9Jb0iIu4HTgM2pa1zZMHWTzO07IHm8CpWlYMvy7Y9ihP2Ryx1wuSkmXklsIFmc/UXNPvx5iR9FLhKUoNm0J2btsLcgm2xzdAiQ82BVQ/dvseyBl57ectaxkFkmTA5uf954PMdllsLrO1nnaVpimYJtF7B02mDcFBZu2FuD4MEUK9ybJ1b0fW925uj1ttIgi2tGZoWalk2SIeYjVLW7a0VUmnLtwda63arvy2380wbDeLpUhxKMnS5BNtc25hHpzMNWgatpZmVXZZAa79dx6ZoEfIJNia4/7mDd3usPcR6fYEOMhsHnbbzrKG22ONBx0kuwbazsfT5IFv4i9QuazXdrE4G2d5nGlP87OljaZ5yad3kEmy7YikP7XyhxjYzu/svzIF7NNv1aZ2oUM9RIhtP/fbLLTQzvxf37zyUDb9/6TCLVUv5NEUbE8zMTrHtueXPP7Ztdhn77LH7ybvtX2C3L939DlYHi6mlLTxQd/2OY/jV0wfywMxBwyxaLeUSbLONSTY/3bzk0vZde+72XHu4Zf2yHW5WZcPoalm/4xju2HoYj2/fm21PDOcyX9Fo0PCoaHaNmODx7Xvv9tjUHs8N9J6DdLaajdqw+o2n5/bhoZ0Hc8fWw/j1zAHMblvGkumlQ3nvOsunKTo/wc7ZpTy3cwlL95zrulyrr61bH1yauh+xbdUy7EGw6bl9+Pm2o9j89H7NUJtezpJtk+w5PdTV1FKmYEtOYv0GzataXhQRF/ZaPubF7Lbm/IfPQddwWxho7Y9nDbeWqp1GY9U0qhH8rXMrePCZlWx68iU8uWOK2enlLJueZOk2mJrudHEMa5cabJImgf9E86z6zcDPJf0gIu7t+qIQ7JqAZY0XPdXe57bP0me7vsViwq2TboeYmGVRxKFIW+dWsOH3L+WR7fvy+BP7ENuXsmx6kj2fgGW/D6Z+N1i3zjjIUmM7GXgwIh4CkPRtmheE6xFsMDErGkwQNGttLe19ba1R014BN2xuvlqaIo+r3Dq34vmRz+3blqOZPdjj9xPs+QRMbWmwbOscez66vbDyVUWWYDsc+G3b/c3A6xcuJGk1sDq5u+tX533m7sGLNxIHAU8UXYg+VKm8VSorVKu8rxz0DbbHk9ete+6KrMeOVOVzAYY4eJBc/3wNgKTbel1GuEyqVFaoVnmrVFaoVnkXzD+wKBFx+jDKUkZZpt97BDiy7f4RyWNmZqWUJdh+Dhwr6WhJe9Cc6+8H+RbLzGzxUpuiyaV5/xy4jubhHpdExD0pL1vMtFxFqVJZoVrlrVJZoVrlrVJZR06dJ4wxM6uuLE1RM7NKcbCZWe0MNdgknS7pl5IelHTBMN972CRdImmLpNIfbyfpSEk3SrpX0j2Sziu6TL1I2lPSrZLuSMr7haLLlEbSpKRfSPph0WVJI+lhSXdJ2jiMwz7qaGh9bMmpV60JTTfTHE19f89Trwok6S3ADuCbEXF80eXpRdKhwKERsUHS3sDtwHtK/NkKWBEROyQtBW4GzouIWwouWleS/jWwCtgnIs4sujy9SHoYWBURlTpodpSGWWN7/tSriJgFWqdelVJE3AQ8WXQ5soiIxyJiQ3J7O82ZsA8vtlTdRVNrWqWlyb/SjlJJOgJ4F3BR0WWx4RhmsHU69aq0O19VSToKeA2wvtiS9JY07TYCW4B1EVHm8n4d+Czw4qs2lFMA10u6PTmV0Rbw4EGFSNoLuAo4PyK2pS1fpIiYj4iTaJ6pcrKkUjb3JZ0JbImI24suSx/eFBGvBf458ImkW8XaDDPYfOpVjpK+qquAyyLi6qLLk1VEPAXcCJT1vMRTgHcn/VbfBk6V9K1ii9RbRDyS/L8FWEuzG8jaDDPYfOpVTpLO+IuBTRHxtaLLk0bSSkn7JbeX0xxQuq/YUnUWEZ+LiCMi4iia2+xPIuKcgovVlaQVyQASklYA7wRKP7I/akMLtoiYA1qnXm0Cvpvh1KvCSLoC+BnwSkmbJX246DL1cArwAZq1iY3JvzOKLlQPhwI3SrqT5g/euogo/WEUFXEIcLOkO4BbgWsi4tqCy1Q6PqXKzGrHgwdmVjsONjOrHQebmdWOg83MasfBZma142Azs9pxsJlZ7fx/VuFb77A4LSkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import fenics as fe\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# --------------------\n",
    "# Functions and classes\n",
    "# --------------------\n",
    "def bottom(x, on_boundary):\n",
    "    return (on_boundary and fe.near(x[1], 0.0))\n",
    "\n",
    "\n",
    "# Strain function\n",
    "def epsilon(u):\n",
    "    return 0.5*(fe.nabla_grad(u) + fe.nabla_grad(u).T)\n",
    "\n",
    "\n",
    "# Stress function\n",
    "def sigma(u, p):\n",
    "    return p*fe.Identity(2) + 2*mu*epsilon(u)\n",
    "\n",
    "\n",
    "# --------------------\n",
    "# Parameters\n",
    "# --------------------\n",
    "E = 70.0e6  # Youngs modulus\n",
    "nu = 0.4999  # Poissons ratio\n",
    "\n",
    "l_x, l_y = 5.0, 5.0  # Domain dimensions\n",
    "n_x, n_y = 20, 20  # Number of elements\n",
    "\n",
    "# Load\n",
    "g_int = -10000000.0\n",
    "\n",
    "# --------------------\n",
    "# Geometry\n",
    "# --------------------\n",
    "mesh = fe.RectangleMesh(fe.Point(0.0, 0.0), fe.Point(l_x, l_y), n_x, n_y)\n",
    "\n",
    "# Definition of Neumann condition domain\n",
    "boundaries = fe.MeshFunction(\"size_t\", mesh, mesh.topology().dim() - 1)\n",
    "boundaries.set_all(0)\n",
    "\n",
    "top = fe.AutoSubDomain(lambda x: fe.near(x[1], l_y))\n",
    "\n",
    "top.mark(boundaries, 1)\n",
    "ds = fe.ds(subdomain_data=boundaries)\n",
    "\n",
    "# --------------------\n",
    "# Function spaces\n",
    "# --------------------\n",
    "P1 = fe.VectorElement('P', fe.triangle, 2)\n",
    "P2 = fe.FiniteElement('P', fe.triangle, 1)\n",
    "element = fe.MixedElement([P1, P2])\n",
    "V = fe.FunctionSpace(mesh, element)\n",
    "g = fe.Constant((0.0, g_int))\n",
    "\n",
    "# --------------------\n",
    "# Boundary conditions\n",
    "# --------------------\n",
    "bc = fe.DirichletBC(V.sub(0), fe.Constant((0.0, 0.0)), bottom)\n",
    "\n",
    "# --------------------\n",
    "# Weak form\n",
    "# --------------------\n",
    "u_test, p_test = fe.TestFunctions(V)\n",
    "u_tr, p_tr = fe.TrialFunctions(V)\n",
    "\n",
    "a = fe.inner(epsilon(u_test), sigma(u_tr, p_tr))*fe.dx\n",
    "a += -p_test*(fe.div(u_tr) + p_tr/lmbda)*fe.dx\n",
    "L = fe.inner(g, u_test)*ds(1)\n",
    "\n",
    "\n",
    "# --------------------\n",
    "# Solver\n",
    "# --------------------\n",
    "sol = fe.Function(V)\n",
    "fe.solve(a == L, sol, bc)\n",
    "\n",
    "# --------------------\n",
    "# Post-process\n",
    "# --------------------\n",
    "# Plot solution\n",
    "fe.plot(sol.sub(0), mode=\"displacement\")\n",
    "plt.show()\n",
    "plot = fe.plot(sol.sub(1))\n",
    "plt.colorbar(plot)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Uncommented other example - Mindlin beam\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEKCAYAAAAiizNaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xl4VNX5wPHvm52wBAgBgQABEpawY0RRUREV0AquFaxrVVygbrWt/GxrtVq1LlQtaFW07oDUBTcWBRUXloDskBDCFtawhTUhy/v7Yy7tGJPMBDJzZ5L38zzzcOfcc859Z5jkzb33zDmiqhhjjDHBEuF2AMYYY+oWSzzGGGOCyhKPMcaYoLLEY4wxJqgs8RhjjAkqSzzGGGOCyhKPMcaYoLLEY4wxJqgs8RhjjAmqKLcDCEXNmjXTlJQUt8MwxpiwsmjRol2qmuSrniWeCqSkpJCZmel2GMYYE1ZEZKM/9exSmzHGmKCyxGOMMSaoLPEYY4wJqoAmHhEZIiJZIpIjIvdXsD9WRCY7++eLSIrXvrFOeZaIDK5Gn8+LyEF/jmGMMSb4ApZ4RCQSGA8MBdKBkSKSXq7aTcBeVU0FxgFPOG3TgRFAN2AIMEFEIn31KSIZQGN/jmGMMcYdgTzj6QfkqGquqh4FJgHDy9UZDrzubE8FBomIOOWTVLVIVdcDOU5/lfbpJKUngd/7eQxjjDEuCGTiaQ1s9nqe55RVWEdVS4ACILGKtlX1OQaYpqrb/DzGT4jIKBHJFJHM/Px8P1+iMcaY6grk93gqOqsov852ZXUqK68oUaqItAKuBM45zjhQ1ZeAlwAyMjJsPXDjqsNHS9i4+zA7DxSRf6CIfYePcrS0jJJSz0ezQWwUDeKiaBIfQ3KTerRpGk+DWPtangkPgfyk5gFtvJ4nA1srqZMnIlFAArDHR9uKyvsAqUCOcxUtXkRynPs6lR3DmJBQVFLKsrwC5ufuZsnmArJ3HGDz3sNoNf/8SWoYS4/WCfRonUDvto05tX1T4mMsGZnQE8hP5UIgTUTaA1vwDBa4ulydacD1wA/AFcBsVVURmQa8IyLPAK2ANGABnrOXn/WpqiuBk451KiIHnaRT6TEC8YKN8deug0V8sWoH01du54d1uykqKQOgY1J9eiQncHnfZDo2r0+LRnE0bxhLk/oxxERGEBUhKHC4qJQDRcXsPniUzXsPs3nPEdbuOMDyLQXMydqJKsRERpCR0oRzOidxUc9WtG5cz90XbYwjYIlHVUtEZAwwA4gEXlXVlSLyMJCpqtOAicCbIpKD5yxkhNN2pYhMAVYBJcBoVS0FqKhPH6FUeAxjgq2opJTpK7YzJXMzP6zbTZlC26bxXH1qW07rkEi/lKY0qR/jV18J8REkxEeT3CSeXm1+OpDzUFEJizft5ZvsfL7J3sXfPlvD3z5bQ7+Upgzr3YrhvVvRMC46EC/RGL+I/fH/cxkZGWpztZmasmXfEf793XqmLspj7+FikpvU47I+rRnSvSVdWzYk0IMsN+4+xMdLt/LRkq2s3XmQ+jGRXH5yMtf1b0dq84YBPbapW0Rkkapm+KxniefnLPGYmrB2xwFe/DqXj5ZsAeCCbi0Y2a8tZ3RsRkRE8Ef0qyrL8gp444eNfLxsK0dLyhjUpTl3nZdGz+TyX38zpvos8ZwASzzmRGzafZgnZ2bx8dKt1IuOZES/Ntw8oENI3WPZfbCIt+dv4tXv1rPvcDEDOydx7/md6ZGc4HZoJoxZ4jkBlnjM8dhz6CjPz17LW/M2Ehkh/PqM9tw8oANN/bxv44YDhcW88cNGXp6by77DxVzWtzW/H9yFkxLi3A7NhCFLPCfAEo+pjrIyZdLCzTz++WoOFpVw1SltuPu8TrRoFD6/vA8UFjPhq3VMnLueyAjh9nM6MuqsDsRFR7odmgkjlnhOgCUe46/V2/bzwAfLWbxpH6d1aMpfh3cnrUX43rDftPswj09fzWfLt9MxqT6PX96TU1Kauh2WCROWeE6AJR7jS3FpGc9/uZbxX60joV40f7yoK5f2aR3wEWrB8nV2Pv/3/nK27DvCtae14/dDOtsQbOOTv4nHvtZsTDXl7DzIPZOXsHxLAZf1bc2ff5FO4/jQvY9zPM7ulMTMe87iqZlZ/Pv7DczJ2smzI3pzcjs7+zEnzhaCM8ZPqsrr32/goufmkrf3MC9e05dnftm71iWdY+rHRvHgxd2Yelt/RODKF3/g2S/WUlJa5nZoJsxZ4jHGD/sLi7n9rcU8OG0lp3dMZMY9ZzGke0u3wwqKk9s15bM7BzCsVyvGfZHNiJfmsXXfEbfDMmHMEo8xPqzZvp/h//yOWat38MCFXXn1hlNo3jB8RqzVhIZx0fxjRB/GXdWL1dv2c/Hz3/LDut1uh2XClCUeY6rwwY95XDL+Ow4VlfDuLadxy1kdas0AguNxaZ9kPhpzBgnx0VwzcT6vzM3FBiiZ6rLEY0wFysqUxz9fwz2Tl9IruTGf3Hkm/drbjXWA1OYN+Wj0GQzq0pxHPl3N3ZOXUFhc6nZYJoxY4jGmnMNHS7jj7cW8+PU6rj61LW/dfGqdu7TmS8O4aF685mR+N7gzHy3ZytUvz2P3wSK3wzJhwhKPMV527C/kqn/NY8aq7fzpF+k8ekl3oiPtx6QiERHC6IGpTPhVX1Zu3c9lL3xPbv5Bt8MyYcB+ooxx5Ow8yKXjvyM3/yCvXJfBTWe2r9P3c/x1YY+WvDvqNA4WlnDZC9+zYL0t8GuqZonHGGDp5n1c+eL3HC1VptzWn0FdW7gdUljp27YJH9xxBk3rx3DNxPnMWrXD7ZBMCAto4hGRISKSJSI5InJ/BftjRWSys3++iKR47RvrlGeJyGBffYrIRBFZKiLLRGSqiDRwym8QkXwRWeI8bg7kazbh57ucXVz98jwaxEXxn9v7062VLQ1wPNomxvOf206na8tG3PbWov+uQ2RMeQFLPCISCYwHhgLpwEgRSS9X7SZgr6qmAuOAJ5y26XiWqO4GDAEmiEikjz7vUdVeqtoT2ASM8TrOZFXt7TxeCcTrNeFp+opt3PjaQpKbxDP1ttNpl1jf7ZDCWpP6Mbx986mcktKEuycv4c15G90OyYSgQJ7x9ANyVDVXVY8Ck4Dh5eoMB153tqcCg8RzUX04MElVi1R1PZDj9Fdpn6q6H8BpXw+wLxeYKn28dCt3vL2Y7q0bMeXW/mG1jEEoaxAbxb9v7MegLs3504creOGrdW6HZEJMIBNPa2Cz1/M8p6zCOqpaAhQAiVW0rbJPEXkN2A50AZ73qne51yW4NhUFKyKjRCRTRDLz8/P9fpEmPE1bupW7Jv1IRkpT3rzpVBLibeblmhQXHckL15zMsF6teGL6GsbPyXE7JBNCApl4KhoOVP4spLI61S33bKjeCLQCVgNXOcUfAynOJbgv+N8Z1k87UX1JVTNUNSMpKamiKqaWmLZ0K3c7See1G06hfqxN0h4I0ZERjLuqN5f0bsWTM7L419d25mM8Apl48gDvs4tkYGtldUQkCkgA9lTR1mefqloKTAYud57vVtVj32x7GTj5uF+RCXuWdIIrMkJ46speXNyrFY99voZX5ua6HZIJAYFMPAuBNBFpLyIxeAYLTCtXZxpwvbN9BTBbPRM/TQNGOKPe2gNpwILK+hSPVPjvPZ6LgTXOc+8phIfhORsyddD0Fdv+m3T+faMlnWCJioxg3C97cVGPljzy6Womfrve7ZCMywL2k6eqJSIyBpgBRAKvqupKEXkYyFTVacBE4E0RycFzpjPCabtSRKYAq4ASYLRzJkMlfUYAr4tIIzyX45YCtzuh3Ckiw5x+9gA3BOo1m9D17dpd3PnuEnq3acxrN5xCfIwlnWCKiozgHyN6U6bKXz9ZRf2YSEb0a+t2WMYltvR1BWzp69pl8aa9XPPKfNo2jWfyqP42kMBFR0vKuOWNTOauzWfCr/rWmTWN6gp/l762mQtMrZa1/QA3vraQpIaxvPHrfpZ0XBYTFcEL1/Sld5vG3PnuEr5ft8vtkIwLLPGYWmvT7sNcO3E+cdERvHXTqTS37+mEhPiYKF694RRSmsVzy+uZLM8rcDskE2SWeEyttPfQUa5/bQFHS8t486ZTadM03u2QjJfG8TG88etTaRwfw/WvLbBZresYSzym1iksLuWWNzLZsu8Ir1yXQacWDd0OyVTgpIQ43rr5VAS44bWFtp5PHWKJx9QqZWXKb99bSubGvYz7ZW8yUmzV0FDWvll9Xr4+gx37Cxn15iJbybSOsMRjapUnZqzh02XbGDu0Cxf1tBFT4aBv2yaMu6o3izbu5b73llJWZiNtaztLPKbWeGveRv71dS7XnNaWUWd1cDscUw0X9mjJ2KFd+GTZNp6eleV2OCbA7Ft0plb4JjufP3+0gnO7NOcvF3ezlUPD0KizOrBh92HGz1lHu6b1+eUpFc7na2oBSzwm7K3fdYgx7yymU4uGPD+yD1GRdiIfjkSEvw7vxpZ9R/i/D5bTLjGeUzskuh2WCQD7CTVh7UBhMbe8kUlkhPDydRk2/1qYi4qM4J9X96FtYjx3vL2YLfuOuB2SCQBLPCZslZYpd09awvpdhxj/q772XZ1aolFcNC9fl8HRkjJGvZHJkaM20q22scRjwtbTM7P4cs1OHrw4ndM7NnM7HFODOiY14LmRfVi1bT+//88ybE7J2sUSjwlLHy3ZwoSv1jGyXxuuPa2d2+GYABjYpTm/G9yZj5du5V/f2Do+tYklHhN2Vm/bzx/+s4xTUprw0LDuNoKtFrv97I78omdLnpi+hjlZO90Ox9QQSzwmrOwvLOb2txbRKC6a8b/qS0yUfYRrMxHh71f0pOtJjbh70hI27znsdkimBthPrQkbqsrv31vG5r1H+OfVfWne0GabrgviY6J44Zq+lKky+p3FFJXYYINwF9DEIyJDRCRLRHJE5P4K9seKyGRn/3wRSfHaN9YpzxKRwb76FJGJIrJURJaJyFQRaeDrGCa8vDJ3PdNXbuf+IV3o197mYKtL2iXW5+kre7Esr4BHPrHV68NdwBKPiEQC44GhQDowUkTSy1W7CdirqqnAOOAJp206nmWwuwFDgAkiEumjz3tUtZeq9gQ2AWOqOoYJLwvW7+Hx6WsY0u0kbh7Q3u1wjAsu6HYSt57VgTfnbeSjJVvcDsecgECe8fQDclQ1V1WPApOA4eXqDAded7anAoPEc6d4ODBJVYtUdT2Q4/RXaZ+quh/AaV8PUB/HMGFi54FCxryzmDZN6vH3K3vaYII67L7BnTklpQlj31/O2h0H3A7HHKdAJp7WwGav53lOWYV1VLUEKAASq2hbZZ8i8hqwHegCPO/jGD8hIqNEJFNEMvPz86vzOk0AlZYpd727hP2Fxbxwzck0irOlq+uy6MgI/nl1X+JjIrn97cUcKipxOyRzHAKZeCr6s7T8t8Aqq1Pdcs+G6o1AK2A1cFU14kBVX1LVDFXNSEpKqqCJccP4OTn8kLubh4d3p2vLRm6HY0JAi0ZxPDuiD7n5B/njhyvcDscch0AmnjzAe3rZZGBrZXVEJApIAPZU0dZnn6paCkwGLvdxDBPiFm7Ywz++yOaS3q248uRkt8MxIeSM1GbcNagTH/y4hfcX57kdjqmmQCaehUCaiLQXkRg8gwWmlaszDbje2b4CmK2euTGmASOcEWntgTRgQWV9ikcq/Pcez8XAGh/HMCGs4HAxd737I22axvPXS+xLoubnxpybSr/2TfnThyvYsOuQ2+GYaghY4nHup4wBZuC59DVFVVeKyMMiMsypNhFIFJEc4F7gfqftSmAKsAqYDoxW1dLK+sRzOe11EVkOLAdaAg9XdQwTulSV+99fxs4DRTw3og8N7b6OqUBkhPCPq3oTFRnBnZN+5GhJmdshGT+J/fH/cxkZGZqZmel2GHXW2/M38sAHKxg7tAu3nt3R7XBMiJu+Yhu3vbWYW8/qwNgLu7odTp0mIotUNcNXPZu5wISU7B0HePjjVQxIa8YtA2z5auPbkO4tufrUtvzrm1y+ybYRqeHAEo8JGYXFpfzmnR9pGBfF07/sRUSE3dcx/vnTRemkNW/AvVOWsutgkdvhGB8s8ZiQ8fjna8jacYCnruxl87CZaqkXE8nzV/dhf2Ex97231NbvCXGWeExI+HbtLv79/QZuOD2Fczo3dzscE4a6nNSI/xvaha+y8nlnwSa3wzFVsMRjXFdwpJjfTV1Kh6T6/GFIF7fDMWHsuv4pnJnajEc+WW1DrEOYJR7jugc/WsHOA0WM+2Vv6sVEuh2OCWMREcKTV/YkKlK4d8oSSsvsklsossRjXPXpsm18uGQrvzk3lV5tGrsdjqkFWibU46/Du7N40z5e/Hqd2+GYCljiMa7Zub+QBz5cTq/kBEYPTHU7HFOLDO/diot6tOQfX2SzcmuB2+GYcizxGFeoKr//zzIKi0t55qreREfaR9HUHBHhkUu60zg+hnsnL6Ww2FYtDSX2025c8fb8TXyVlc/YoV3pmNTA7XBMLdSkfgx/v6InWTsO8MysbLfDMV4s8Zig27znMH/7bDUD0ppx7Wnt3A7H1GIDOzfn6lPb8vLcXBast0npQ4UlHhNUxyYAjRDh8ct72uwEJuAeuLArrRvX4w//WcaRo3bJLRREVbVTRJ7zo4/9qvrHGorH1HKTFm7mu5zdPHppd1o3rud2OKYOqB8bxROX9+RXr8xn3BfZ/J9NJOo6X2c8w4FFPh6XV9raGC9b9x3h0U9X079DIiNPaet2OKYOOSO1GSP7teWVubn8uGmv2+HUeVWe8QDjVPX1qiqISJMajMfUUqrK2PeXU1qmPGGX2IwLxl7Yha+ydvK7qcv49M4ziY2yLyu7pcozHlX9h68O/KljzH8Wb+Hr7Hz+MKQzbRPj3Q7H1EGN4qJ57LIe5Ow8yHNfrnU7nDrNr8EFzlLTz4jI+yIy7djDj3ZDRCRLRHJE5GcrfzpLW0929s8XkRSvfWOd8iwRGeyrTxF52ylfISKviki0U36OiBSIyBLn8Wd/XrOpOTv2F/Lwxys5JaUJ1/VPcTscU4ed07k5V5yczItf57Jii32x1C3+jmr7ENgAPA887fWolIhEAuOBoUA6MFJE0stVuwnYq6qpwDjgCadtOjAC6AYMASaISKSPPt8GugA9gHrAzV7HmauqvZ3Hw5igUVUe+GAFRSVl/P0KW2PHuO9PF6WTWD+G+95bastlu8TfxFOoqs+p6hxV/frYw0ebfkCOquaq6lFgEp7BCt6GA8fuIU0FBomIOOWTVLVIVdcDOU5/lfapqp+pA1gAJPv52kwATVu6lS9W7+C+CzrTvll9t8MxhoT4aB69tAdrth9gwlc5bodTJ/mbeJ4VkQdFpL+I9D328NGmNbDZ63meU1ZhHVUtAQqAxCra+uzTucR2LTDdq7i/iCwVkc9FpFtFwYrIKBHJFJHM/HxbPrcm7Dl0lIc+XkWvNo359Znt3Q7HmP86P70Fw3u3YvycHNbuOOB2OHWOv4mnB3AL8Dj/u8z2lI82FV1TKT9HeWV1qlvubQLwjarOdZ4vBtqpai88lwo/rChYVX1JVTNUNSMpKamiKqaa/vbZavYfKeaJy3sQaZfYTIj58y/SqR8bxdj3l1NmyycElb+J51Kgg6qeraoDnce5PtrkAW28nicDWyurIyJRQAKwp4q2VfYpIg8CScC9x8pUdb+qHnS2PwOiRaSZj9jNCfp+3S6mLsrjlrM60OWkRm6HY8zPJDaI5YELu5K5cS+TFm723cDUGH8Tz1KguoulLATSnBFxMXgGC5QfCTcNuN7ZvgKY7dyjmQaMcEa9tQfS8Ny3qbRPEbkZGAyMVNX/3jEUkZOc+0aISD/nNe+u5msx1VBYXMoDH6ygXWI8dw1KczscYyp1xcnJ9O+QyGOfr2bn/kK3w6kz/E08LYA1IjLD3+HUzj2bMcAMYDUwRVVXisjDIjLMqTYRSBSRHDxnKfc7bVcCU4BVeO7VjFbV0sr6dPp60Ynzh3LDpq8AVojIUuA5YIST3EyAjJ+Tw/pdh3j0kh7ERduX9EzoEhEevbQ7RSVlPPTJKrfDqTPEn9/BInJ2ReV+jGwLSxkZGZqZmel2GGFp7Y4DXPjcXH7RsxXjrurtdjjG+OX5L9fy9KxsXr0hg3O7tHA7nLAlIotUNcNXPV9T5gC1N8GYmlVW5pkWp35sFH+8yCZiNOHj1rM7Mm3pVv704UpOvSeR+rF+/Wo0x6nKS20i8omvDvypY+qGSQs3k7lxLw9c2JXEBrFuh2OM32KiInjssh5s2XeEcbZoXMD5Sutn+riXI3hmEDB13M79hTz2uWfm6StOtu/umvCTkdKUq09ty6vfreeSPq3p3jrB7ZBqLV+Jp/xMAxU5WhOBmPD20CerKCop49FLu+MMIjQm7PxhSBdmrdrB/e8v46PRZ9r3zwKkysRj93aMP77JzufTZdu49/xOdEhq4HY4xhy3hHrRPHhxOmPe+ZG35m3k+tNT3A6pVrKlr80JKSop5cFpK2nfrD63nt3B7XCMOWEX9WjJmanNeGpmFvkHitwOp1ayxGNOyMvf5LJ+1yEeGtbNFtYytYKI8NDwbhQWl/LYZ6vdDqdWssRjjtvmPYf555wcLuxxEmd1svntTO3RMakBo87qwPs/bmF+rk10UtP8XQjuDBGZJSLZIpIrIutFJDfQwZnQ9vAnq4gQ4Y8X2cBGU/uMGZhG68b1+NNHKygutXV7apK/ZzwTgWeAM4FTgAznX1NHzV6zg1mrdnDnoDRaNa7ndjjG1Lh6MZE8eHE62TsO8vr3G9wOp1bxN/EUqOrnqrpTVXcfewQ0MhOyCos9AwpSmzfg12fYOjum9jo/vQXndmnOuFnZbC+wSURrir+JZ46IPFnNheBMLfXCV+vYvOcIDw/vRkyU3SY0tZeI8JeLu1FSpjzyqU0iWlP8nZDoVOdf78nfFPC1Jo+pZTbuPsQLX69jWK9WnN7RljUytV/bxHjuOCeVcV9kM7LfLs5Itc/9ifLrz1Wvxd8GVmMhOFPLqCoPTltJTGQED9gkoKYOufXsDrRLjOdPH62gqKTU7XDCnr+j2hJE5BkRyXQeT4uITWRUx8xctYOvsvK5+7w0WjSKczscY4ImLjqSh4Z1Izf/EBO/Xe92OGHP3wv0rwIHgF86j/3Aa4EKyoSewuJS/vrJKjq3aMgNNo2IqYPO6dycC9Jb8M/ZOTbQ4AT5m3g6quqDqprrPB4CfM6PIiJDRCRLRHJE5P4K9seKyGRn/3wRSfHaN9YpzxKRwb76FJG3nfIVIvKqiEQ75SIizzn1l9mgiOPz8je55O09wl+GdSMq0gYUmLrpjxelU1KmPDF9jduhhDV/f4McEZEzjz0RkTOAI1U1EJFIYDwwFM/SCSNFpPw3DW8C9qpqKjAOeMJpmw6MALoBQ4AJIhLpo8+3gS5AD6AecLNTPhRIcx6jgBf8fM3Gsa3gCBO+WseFPU6if8dEt8MxxjVtE+O5ZUB7PvhxC4s27nE7nLDlb+K5HRgvIhtEZCPwT+A2H236ATnOGdJRYBI/X2ZhOPC6sz0VGCSeOfWHA5NUtUhV1wM5Tn+V9qmqn6kDWAAkex3jDWfXPKCxiLT083Ub4LHP1lCmytihNqDAmDvOSeWkRnH8ZdoqysrU7XDCkr+j2paoai+gJ9BDVfuo6lIfzVoDm72e5zllFdZR1RKgAEisoq3PPp1LbNcC06sRByIy6tjgifz8fB8vre5YuGEP05Zu5dazOtCmabzb4RjjuvqxUYy9sAvLtxTw3qLNvhuYn6nyezwico2qviUi95YrB0BVn6mqeQVl5f88qKxOZeUVJcryfU4AvlHVudWIA1V9CXgJICMjw/6MAUrLlL9MW0nLhDhuO6ej2+EYEzKG9WrFmz9s5O/TsxjSvSUJ9aLdDims+Drjqe/827CCh68Vv/KANl7Pk4GtldURkSggAdhTRdsq+xSRB4EkwDtR+hOHqcB7mZtZuXU/Yy/sSnyMv981Nqb2ExH+Mqwbew4f5bkv17odTtjxtQLpv5zNL1T1O+99zgCDqiwE0kSkPbAFz2CBq8vVmQZcD/wAXAHMVlUVkWnAOyLyDNAKz8CABXjOXirsU0RuBgYDg1S1rNwxxojIJDwzMBSo6jYfsdd5BUeKeXJGFqekNOHinnZLzJjyurdOYMQpbXj9+w2M7NeG1OYN3Q4pbPg7uOB5P8v+y7lnMwaYAawGpqjqShF5WESGOdUmAokikoPnLOV+p+1KYAqwCs+9mtGqWlpZn05fLwItgB9EZImI/Nkp/wzIxTNA4WXgDj9fc5323Jdr2XP4KA9e3O2/l1aNMT913wWdqRcTyUMfr8Izrsn4w9c9nv7A6UBSufs8jQCfy02q6md4fvF7l/3Za7sQuLKSto8Cj/rTp1Ne4WtxRrmN9hWr+Z+cnQd4/fsNjDilDd1b2wQVxlQmsUEs95zXiYc/WcUXq3dyfnoLt0MKC77OeGLw3MuJ4qf3d/bjuTRmahlV5eFPVlMvJpL7LujsdjjGhLxr+7cjrXkDHvl0lc3j5idf93i+Br4WkX+r6sYgxWRcNHvNTr7JzuePF3UlsUGs2+EYE/KiIyP488XpXDtxAa9+u4HbbQSoT/7e43lFRBofeyIiTURkRoBiMi45WlLGI5+upmNSfa63+diM8duAtCTO69qc8XNy2HWwyO1wQp6/iaeZqu479kRV9wLNAxOSccub8zayftch/nhROtE2H5sx1TL2wq4UFpcybla226GEPH9/u5SJSNtjT0SkHRV8CdOEr33O9xEGpDXjnM5JbodjTNjpmNSAa05rx7sLNpG1/YDb4YQ0fxPPA8C3IvKmiLwJfAOMDVxYJtie/XItBwqLeeCirjZ82pjjdNegNBrGRfPoZ6vdDiWk+TtX23SgLzAZz/drTlZVu8dTS+TmH+TNHzZy1Slt6HJSI7fDMSZsNakfw52D0vgmO5+vsna6HU7I8ncFUsGzPEFfVf0YiBeRfgGNzATN45+vITYqgnvO7+R2KMaEvWtPa0dKYjyPfrqaktIy3w3qIH8vtU0A+gMjnecH8KyLY8LcvNzdzFy1gzsGptK8oS0MnP6iAAAcAklEQVRnbcyJiomKYOyFXVm78yDvLrTZqyvib+I5VVVHA4Xw31FtMQGLygRFWZnyyKeraJUQx01ntnc7HGNqjQvSW3Bq+6aMm5XN/sJit8MJOf4mnmJn9U8FEJEkwM4hw9wHP25hxZb9/H5IF+Kifc6AZIzxk4jwp1+ks/fwUcbPyXE7nJDjb+J5DvgAaC4ijwLfAn8LWFQm4I4cLeXJGVn0Sk5gWK9WbodjTK3TvXUCl/dN5rVvN7B5z2G3wwkp/o5qexv4PfAYsA24RFXfC2RgJrBenpvL9v2F/PEX6URE2PBpYwLhvgs6ExkhPD59jduhhJQqE4+IND32AHYC7wLvADucMhOGduwv5IWv1jG0+0mckmL/jcYEykkJcdx6dgc+XbaNzA173A4nZPg641kEZDr/HtvO9No2YejpmVmUlJVx/9AubodiTK036qwOnNQojr9+soqyMpvwBXwnnmtVtQPQVVXbq2oH59HeKa+SiAwRkSwRyRGR+yvYHysik53980UkxWvfWKc8S0QG++pTRMY4ZSoizbzKzxGRAmdxOO8F4uqklVsLeG9RHjecnkK7xPq+GxhjTkh8TBT3De7M0rwCPl1uix+D78TzrPPv99Xt2BkFNx4YCqQDI0UkvVy1m4C9qpoKjAOecNqm41nWuhueL65OEJFIH31+B5wHVLR8w1xV7e08Hq7ua6ktVJW/fbaaxvWiGXNumtvhGFNnXNqnNV1bNuLvM9bYmj34TjzFIvIakCwiz5V/+GjbD8hR1VxVPQpMAoaXqzMceN3ZngoMcmZJGA5MUtUiVV2PZ9nqflX1qao/quoGv151HfV1dj7f5ezmzkFpJNSLdjscY+qMyAhh7NAubN5zhLfnbXI7HNf5Sjy/AGYAR/jffR7vR1VaA95f281zyiqso6olQAGQWEVbf/qsSH8RWSoin4tINz/q1zqlZcrjn6+hbdN4fnVqO7fDMabOOatTEgPSmvH87LV1/kulvlYg3QVMEpHVqrq0mn1XNEa3/J21yupUVl5RovR1t24x0E5VD4rIhcCHwM+uM4nIKGAUQNu2bcvvDnsf/LiFNdsP8PzIPsRE2Vo7xrjhD0O68Ivnv+XFr9bx+yF1d3CPv7+BjojIlyKyAkBEeorIH320yQPaeD1PBrZWVkdEooAEYE8Vbf3p8ydUdb+qHnS2PwOivQcfeNV7SVUzVDUjKal2rUdTWFzKMzOz6JmcwEU9WrodjjF1VvfWCVzapzUTv13PtoIjbofjGn8Tz8t41t8pBlDVZXhu/ldlIZAmIu1FJMapP61cnWnA9c72FcBsVVWnfIQz6q09njOUBX72+RMicpJz3whnRu0IYLcfr7nWeP37DWwtKOT+oV3sy6LGuOy3F3RCFZ6ZWXdXKvU38cSr6oJyZSVVNXDu2YzBc49oNTBFVVeKyMMiMsypNhFIFJEc4F7gfqftSjzr/qwCpgOjVbW0sj4BROROEcnDcxa0TERecY5xBbBCRJbimfpnhJPc6oR9zlxR53RO4vSOPzvRM8YEWXKTeG44I4Wpi/NYs32/2+G4Qvz5HSwin+P5hf+eqvYVkSuAm1R1aKADdENGRoZmZtaO78f+7bPVvDw3l8/uHEDXlrbImzGhoOBwMWc9OYc+bRvz7xtrz9JmIrJIVTN81fP3jGc08C+gi4hsAe4Gbj+B+EwQbNl3hH9/v4HL+yZb0jEmhCTERzNmYCpfZeXzXc4ut8MJOn8nCc1V1fOAJKCLqp5p35kJfU/PzEKAe21lUWNCzrX929G6cT0e+3x1nZtKp8rh1CJybyXlAKjqMwGIydSAVVv388GPWxh1VgdaNa7ndjjGmHLioiP53eDO3D15CR8v28rw3v58JbF28HXG09B5ZOC5tHbsS5y34ZmyxoSox6evoVFcNHecnep2KMaYSgzr1YpurRrx9+lZdWoqnSoTj6o+pKoPAc2Avqr6W1X9LXAyntFjJgR9l7OLb7LzGTMwlYR4mxrHmFAVESH834Vd2bLvCG/+UNE0k7WTv4ML2gJHvZ4fBVJqPBpzwsrKlMc+X03rxvW4tr9NjWNMqDsjtRlnd0ri+dk5FByuG1Pp+Jt43gQWiMhfRORBYD7/m9zThJCPl21lxZb93De4E3HRkW6HY4zxw/1Du7C/sJgJX+e4HUpQ+Duq7VHgRmAvsA+4UVUfC2RgpvqKSkp5amYW6S0bMbxX3blRaUy469qyEZf2bs2/v9vA9oJCt8MJOL9ni1TVxar6rPP4MZBBmePz9rxNbN5zxKbGMSYM3XN+J8pUefbLtW6HEnA2TXEtcaCwmOdnr+XM1Gac1al2TXJqTF3QxlmyZErmZnLzD7odTkBZ4qklXp67nr2Hi/lDHZ5q3ZhwN3pgKrFRETxdyycQtcRTC+w6WMTEublc1KMlPZIT3A7HGHOckhrGcvOADny6fBvL8va5HU7AWOKpBcbPyaGwpIx7L7CpcYwJd7cMaE+T+GienJHldigBY4knzOXtPczb8zZxRd9kOiY1cDscY8wJahgXzeiBqcxdu6vWTiBqiSfM/eOLtSBw13k/W83bGBOmrjmtHa0S4vj79DXUxuXDLPGEsbU7DvD+4jyuO62dTQRqTC0SFx3JPed3YmleAdNXbHc7nBoX0MQjIkNEJEtEckTk/gr2x4rIZGf/fBFJ8do31inPEpHBvvoUkTFOmYpIM69yEZHnnH3LRKRv4F5xcD09M5v4mCjuGGgTgRpT21zWN5m05g14cmYWJaVlbodTowKWeEQkEhgPDMUzk/VIESk/o/VNwF5VTQXGAU84bdOBEUA3YAgwQUQiffT5HXAeUH6mvaFAmvMYBbxQk6/TLUs272P6yu3cMqADTevHuB2OMaaGRUYI9w3uTG7+If6zOM/tcGpUIM94+gE5ziJyR4FJwPBydYbzvznfpgKDxLPYz3BgkqoWqep6IMfpr9I+VfXHShanGw68oR7zgMYi0rJGX6kLnpyxhsT6Mdw0oL3boRhjAuSC9Bb0aduYcbPWUlhce5ZNCGTiaQ1s9nqe55RVWEdVS4ACILGKtv70eTxxhJVv1+7iu5zdjB6YSoPYKtfyM8aEMRHhD0O6sH1/IW/8sMHtcGpMIBNPRZOFlR+eUVmd6pafaByIyCgRyRSRzPz8fB9dukdVeXLGGlo3rsevTmvrdjjGmAA7rUMiZ3dKYvycdRQcqR3LJgQy8eQBbbyeJwNbK6sjIlFAArCnirb+9Hk8caCqL6lqhqpmJCWF7lxnM1ZuZ2leAXedl0ZslC17YExd8LvBnSk4UszL3+S6HUqNCGTiWQikiUh7EYnBM1hgWrk604Drne0rgNnqGbQ+DRjhjHprj2dgwAI/+yxvGnCdM7rtNKBAVbfVxAsMtpLSMp6ckUVq8wZc1iesrxYaY6qhe+sELu7Vionfrmfn/vBfNiFgice5ZzMGmAGsBqao6koReVhEhjnVJgKJIpID3Avc77RdCUwBVgHTgdGqWlpZnwAicqeI5OE5o1kmIq84x/gMyMUzQOFl4I5AveZAe//HLazLP8R9F3QiKtK+gmVMXfLb8ztRXFrG87PDf7E4qY3fij1RGRkZmpmZ6XYYP1FYXMq5T31FUsNYPhx9Bp7Bf8aYuuSPHy5n0oLNfPnbs2mXWN/tcH5GRBapaoavevZnc5h4e/4mthYU8vshXSzpGFNH/ebcNCIjhGe/CO/F4izxhIGDRSWMn5PDGamJnJHazHcDY0yt1KJRHNefnsIHS7aQveOA2+EcN0s8YeCVubnsOXSU3w22Rd6MqetuO7sj9WOieCaMF4uzxBPidh8s4pW56xnS7SR6t2nsdjjGGJc1rR/DTWe2Z/rK7SzPK3A7nONiiSfETfhqHYePlnDfYFvkzRjjcfOA9jSOj+apmeG5WJwlnhC2reAIb87byGV9k0lt3tDtcIwxIaJhXDS3n92Rr7PzWbB+j9vhVJslnhD2z9k5qCp3DbJF3owxP3Vd/xSSGsby1IyssFsszhJPiNq85zCTF27mqlPa0KZpvNvhGGNCTL2YSH5zbioLNuzhm7XhtUS2JZ4Q9eyXa4mIEMYMtLMdY0zFRpzSltaN64XdWY8lnhCUs/Mg7y/O49rT2nFSQpzb4RhjQlRMVAR3n5fG8i0FzFgZPktkW+IJQf/4Ipu46EhuP6ej26EYY0LcpX1a0yGpPk/PzKa0LDzOeizxhJjV2/bzybJt3HhGCs0axLodjjEmxEVFRnDv+Z1Yu/Mg05ZucTscv1jiCTHPzMqmYVwUowbY2Y4xxj8Xdm9JestGjJu1luLSMrfD8ckSTwhZunkfs1bt4JYBHUiIj3Y7HGNMmIiIEO4b3IlNew4zJXOz2+H4ZIknhDw1M4sm8dHceEaK26EYY8LMwM7N6du2Mc9/mUNhcanb4VTJEk+ImJ+7m7lrd3H7OR1pGGdnO8aY6hER7hvcme37C3lr3ka3w6lSQBOPiAwRkSwRyRGR+yvYHysik53980UkxWvfWKc8S0QG++rTWQ57voisdfqMccpvEJF8EVniPG4O5Gs+HqrK0zOzSWoYy7WnpbgdjjEmTJ3esRlnpjZjwlfrOFhU4nY4lQpY4hGRSGA8MBRIB0aKSHq5ajcBe1U1FRgHPOG0TQdGAN2AIcAEEYn00ecTwDhVTQP2On0fM1lVezuPVwgx3+bsYsGGPYwZmEq9mEi3wzHGhLH7Bndmz6GjvPbterdDqVQgz3j6ATmqmquqR4FJwPBydYYDrzvbU4FB4lleczgwSVWLVHU9kOP0V2GfTptznT5w+rwkgK+txqgqT83MpnXjeozo18btcIwxYa53m8ac17UFL83NpeBwsdvhVCiQiac14D28Is8pq7COqpYABUBiFW0rK08E9jl9VHSsy0VkmYhMFZGQ+u3+xeqdLN28jzsHpRIbZWc7xpgT99sLOnGwqIR/fbPO7VAqFMjEIxWUlf9abWV1aqoc4GMgRVV7Al/wvzOsnwYiMkpEMkUkMz8/v6IqNa6sTHl6ZhYpifFc1jc5KMc0xtR+XVs24uKerXjtuw3sPFDodjg/E8jEkwd4n10kA1srqyMiUUACsKeKtpWV7wIaO3385FiqultVi5zyl4GTKwpWVV9S1QxVzUhKSqrGyzx+n63YxprtB7j7vE5ER9oAQ2NMzbnn/E4cLS1jwpzQO+sJ5G+7hUCaM9osBs9ggWnl6kwDrne2rwBmq2eK1WnACGfUW3sgDVhQWZ9OmzlOHzh9fgQgIi29jjcMWF3Dr/O4lJSW8cysbNKaN+DiXq3cDscYU8u0b1afK/om8878TWzZd8TtcH4iYInHud8yBpiB55f9FFVdKSIPi8gwp9pEIFFEcoB7gfudtiuBKcAqYDowWlVLK+vT6esPwL1OX4lO3wB3ishKEVkK3AncEKjXXB0fLtlKbv4h7j2/E5ERFV0pNMaYE3PneZ5lVf45e63LkfyUhNMaDsGSkZGhmZmZAev/aEkZg575ikZx0XzymzPxDMozxpia9+BHK3hr/iZm//Zs2iXWD+ixRGSRqmb4qmc3Flzw3qLNbN5zhPsu6GxJxxgTUKMHphIVITz7Zeic9VjiCbLC4lKe/zKHvm0bc07n4AxiMMbUXc0bxXH96Sl8+OMWcnYecDscwBJP0L0zfxPb9xfa2Y4xJmhuPasD9aIjGfdFaJz1WOIJosNHS5jwVQ79OyRyemozt8MxxtQRiQ1iufGM9ny6bBurtu53OxxLPMH07+83sOvgUe4b3MntUIwxdcwtAzrQMC6KZ2Zlux2KJZ5g2V9YzL++zmVg5yRObtfU7XCMMXVMQnw0owZ04IvVO1i6eZ+rsVjiCZKJc9dTcKSYe8/v7HYoxpg66sYz29MkPpqnXT7rscQTBHsPHWXit+sZ0u0keiQnuB2OMaaOahAbxW1nd+Sb7HwWrN/jWhyWeILgxW/WcehoCfdeYPd2jDHuuq5/Cs0axPL0zCzcmkDAEk+A7TxQyOvfb2B4r1Z0atHQ7XCMMXVcvZhIxgzsyPz1e/h+3W5XYrDEE2AT5qyjuFS56zw72zHGhIaRp7alVUIcT7l01mOJJ4C27jvCO/M3cUXfZNo3C+wcScYY46/YqEjGnJvGj5v2MSdrZ9CPb4kngJ6fvRZF+c2gVLdDMcaYn7gyI5m2TeN5emY2ZWXBPeuxxBMgG3YdYkpmHlf3a0tyk3i3wzHGmJ+IjozgrkFprNy6nxkrtwf12JZ4AuS5L9cSFSGMHmhnO8aY0HRJn9Z0TKrPuC+yKQ3iWY8lngBYu+MAHyzZwvWnp9C8UZzb4RhjTIUiI4S7z+tE9o6DfLJsa9COG9DEIyJDRCRLRHJE5P4K9seKyGRn/3wRSfHaN9YpzxKRwb76dJbDni8ia50+Y3wdI1D+8cVa4qMjufWsDoE+lDHGnJCLerSky0kNGTcrm5LSsqAcM2CJR0QigfHAUCAdGCki6eWq3QTsVdVUYBzwhNM2HRgBdAOGABNEJNJHn08A41Q1Ddjr9F3pMQJlxZYCPl2+jV+f2Z7EBrGBPJQxxpywiAjh3vM7sWH3Yd5fvCU4xwxg3/2AHFXNVdWjwCRgeLk6w4HXne2pwCDxLFIzHJikqkWquh7IcfqrsE+nzblOHzh9XuLjGAExblY2jeKiuHmAne0YY8LD+ekt6JmcwLNfruVoSeDPegKZeFoDm72e5zllFdZR1RKgAEisom1l5YnAPqeP8seq7Bg1bvGmvXy5Zie3nt2RhHrRgTiEMcbUOBHhtxd0Zsu+I0zO3Oy7wQkKZOKp6Kyi/LCJyurUVLm/cSAio0QkU0Qy8/PzK2jinwFpzbjh9JTjbm+MMW44K60Zw3q1okl84P9ojgpg33lAG6/nyUD5YRPH6uSJSBSQAOzx0bai8l1AYxGJcs5qvOtXdoyfUNWXgJcAMjIyjmtcYd+2TXjzplOPp6kxxrhKRHhuZJ+gHCuQZzwLgTRntFkMnsEC08rVmQZc72xfAcxWz8RB04ARzoi09kAasKCyPp02c5w+cPr8yMcxjDHGuCBgZzyqWiIiY4AZQCTwqqquFJGHgUxVnQZMBN4UkRw8ZyEjnLYrRWQKsAooAUarailARX06h/wDMElEHgF+dPqmsmMYY4xxh9gf/z+XkZGhmZmZbodhjDFhRUQWqWqGr3o2c4ExxpigssRjjDEmqCzxGGOMCSpLPMYYY4LKEo8xxpigslFtFRCRfGDjcTZvhucLraEmVOOC0I3N4qoei6t6amNc7VQ1yVclSzw1TEQy/RlOGGyhGheEbmwWV/VYXNVTl+OyS23GGGOCyhKPMcaYoLLEU/NecjuASoRqXBC6sVlc1WNxVU+djcvu8RhjjAkqO+MxxhgTVJZ4qkFEhohIlojkiMj9FeyPFZHJzv75IpLitW+sU54lIoODHNe9IrJKRJaJyJci0s5rX6mILHEe5ZetCHRcN4hIvtfxb/bad72IrHUe15dvG+C4xnnFlC0i+7z2BfL9elVEdorIikr2i4g858S9TET6eu0L5PvlK65fOfEsE5HvRaSX174NIrLceb9qdOZdP+I6R0QKvP6//uy1r8rPQIDj+p1XTCucz1RTZ19A3i8RaSMic0RktYisFJG7KqgTvM+XqtrDjweeZRjWAR2AGGApkF6uzh3Ai872CGCys53u1I8F2jv9RAYxroFAvLN9+7G4nOcHXXy/bgD+WUHbpkCu828TZ7tJsOIqV/83eJbfCOj75fR9FtAXWFHJ/guBz/GsqnsaMD/Q75efcZ1+7HjA0GNxOc83AM1cer/OAT450c9ATcdVru7FeNYIC+j7BbQE+jrbDYHsCn4eg/b5sjMe//UDclQ1V1WPApOA4eXqDAded7anAoNERJzySapapKrrgRynv6DEpapzVPWw83QenhVaA82f96syg4FZqrpHVfcCs4AhLsU1Eni3ho5dJVX9hgpWx/UyHHhDPebhWXW3JYF9v3zGparfO8eF4H2+/Hm/KnMin82ajisony9V3aaqi53tA8BqoHW5akH7fFni8V9rYLPX8zx+/h/33zrqWYK7AEj0s20g4/J2E56/ao6JE5FMEZknIpfUUEzViety57R+qogcW9Y8JN4v55Jke2C2V3Gg3i9/VBZ7IN+v6ir/+VJgpogsEpFRLsTTX0SWisjnItLNKQuJ90tE4vH8Av+PV3HA3y/x3ALoA8wvtyton6+ArUBaC0kFZeWHBFZWx5+2x8vvvkXkGiADONuruK2qbhWRDsBsEVmuquuCFNfHwLuqWiQit+E5WzzXz7aBjOuYEcBUdVa/dQTq/fKHG58vv4nIQDyJ50yv4jOc96s5MEtE1jhnBMGwGM8ULgdF5ELgQyCNEHm/8Fxm+05Vvc+OAvp+iUgDPInublXdX353BU0C8vmyMx7/5QFtvJ4nA1srqyMiUUACnlNuf9oGMi5E5DzgAWCYqhYdK1fVrc6/ucBXeP4SCkpcqrrbK5aXgZP9bRvIuLyMoNxlkAC+X/6oLPZAvl9+EZGewCvAcFXdfazc6/3aCXxAzV1i9klV96vqQWf7MyBaRJoRAu+Xo6rPV42/XyISjSfpvK2q71dQJXifr5q+iVVbH3jODnPxXHo5dkOyW7k6o/np4IIpznY3fjq4IJeaG1zgT1x98NxMTStX3gSIdbabAWupoZusfsbV0mv7UmCes90UWO/E18TZbhqsuJx6nfHc6JVgvF9ex0ih8pvlF/HTm78LAv1++RlXWzz3LU8vV14faOi1/T0wJIhxnXTs/w/PL/BNznvn12cgUHE5+4/9UVo/GO+X87rfAP5RRZ2gfb5q7I2uCw88oz6y8fwSf8ApexjPWQRAHPCe80O4AOjg1fYBp10WMDTIcX0B7ACWOI9pTvnpwHLnB285cFOQ43oMWOkcfw7Qxavtr533MQe4MZhxOc//Ajxerl2g3693gW1AMZ6/Mm8CbgNuc/YLMN6JezmQEaT3y1dcrwB7vT5fmU55B+e9Wur8Pz8Q5LjGeH2+5uGVGCv6DAQrLqfODXgGHHm3C9j7hefypwLLvP6fLnTr82UzFxhjjAkqu8djjDEmqCzxGGOMCSpLPMYYY4LKEo8xxpigssRjjDEmqCzxGGOMCSpLPMaEMBFJEZEjIrKkmu2ucqa3/yRQsRlzvCzxGBP61qlq7+o0UNXJwM0+KxrjAks8xrhERE5xZuaOE5H6zgJd3X20SRGRNSLyirOI2Nsicp6IfOcs0hW0udCMOV42O7UxLlHVheJZxfQRoB7wlqpWuGplOanAlcAoYCFwNZ4pUYYB/wcEe7kGY6rFEo8x7noYT/IoBO70s816VV0OICIrgS9VVUVkOZ7JKY0JaXapzRh3NQUa4FmOOM7PNkVe22Vez8uwPyZNGLDEY4y7XgL+BLwNPOFyLMYEhf11ZIxLROQ6oERV3xGRSOB7ETlXVWf7amtMOLNlEYwJYSKSAnyiqlWOdquk7TnAfar6ixoOy5gTYpfajAltpUDC8XyBFJiAZ4E2Y0KKnfEYY4wJKjvjMcYYE1SWeIwxxgSVJR5jjDFBZYnHGGNMUFniMcYYE1T/D+5xfgZhcIHvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f10ee073128>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import fenics as fe\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# --------------------\n",
    "# Parameters\n",
    "# --------------------\n",
    "E = 50.0e9\n",
    "G = 25.0e9\n",
    "b = 0.1\n",
    "h = 0.01\n",
    "Ar = b*h\n",
    "I = (1.0/12.0)*b*h**3\n",
    "\n",
    "k = 1\n",
    "f_lin = 1.0\n",
    "n = 80000\n",
    "l = 2.0\n",
    "F = 1.0\n",
    "\n",
    "# --------------------\n",
    "# Define geometry\n",
    "# --------------------\n",
    "mesh = fe.IntervalMesh(n, 0, l)\n",
    "\n",
    "# --------------------\n",
    "# Define spaces\n",
    "# --------------------\n",
    "P1 = fe.FiniteElement('P', fe.interval, 2)\n",
    "element = fe.MixedElement([P1, P1, P1])\n",
    "V = fe.FunctionSpace(mesh, element)\n",
    "W = fe.FunctionSpace(mesh, 'P', 2)\n",
    "\n",
    "d_u, d_w, d_phi = fe.TestFunctions(V)\n",
    "u, w, phi = fe.TrialFunctions(V)\n",
    "u_ = fe.Function(V)\n",
    "\n",
    "# --------------------\n",
    "# Boundary conditions\n",
    "# --------------------\n",
    "bc_u = fe.Constant(0.0)\n",
    "\n",
    "tol = 1e-14\n",
    "def left_end(x):\n",
    "    return fe.near(x[0], 0.0)\n",
    "\n",
    "def right_end(x):\n",
    "    return fe.near(x[0], 2.0)\n",
    "\n",
    "bc1 = fe.DirichletBC(V.sub(0), fe.Constant(0.0), left_end)\n",
    "bc2 = fe.DirichletBC(V.sub(1), fe.Constant(0.0), left_end)\n",
    "bc3 = fe.DirichletBC(V.sub(1), fe.Constant(0.0), right_end)\n",
    "bc4 = fe.DirichletBC(V.sub(0), fe.Constant(0.0), right_end)\n",
    "bc = [bc1, bc2, bc3, bc4]\n",
    "\n",
    "# --------------------\n",
    "# Initialization\n",
    "# --------------------\n",
    "X = fe.Function(V)\n",
    "f = fe.Function(V)\n",
    "f.interpolate(fe.Constant((0.0, f_lin, 0.0)))\n",
    "\n",
    "# --------------------\n",
    "# Solution\n",
    "# --------------------\n",
    "dx_shear = fe.dx(metadata={\"quadrature_degree\": 2})\n",
    "A = d_u.dx(0)*E*Ar*u.dx(0)*fe.dx + \\\n",
    "    d_w.dx(0)*G*Ar*w.dx(0)*dx_shear + d_w.dx(0)*G*Ar*phi*dx_shear + \\\n",
    "    d_phi.dx(0)*E*I*phi.dx(0)*fe.dx + d_phi*G*Ar*(phi + w.dx(0))*dx_shear + \\\n",
    "    fe.Constant(0.0)*d_w*fe.dx\n",
    "\n",
    "A_ass, b_ass = fe.assemble_system(fe.lhs(A), fe.rhs(A), bc)\n",
    "\n",
    "pointForce = fe.PointSource(V.sub(1), fe.Point(1.0), 1)\n",
    "pointForce.apply(b_ass)\n",
    "\n",
    "fe.solve(A_ass, X.vector(), b_ass)\n",
    "\n",
    "\n",
    "u_, w_, phi_ = X.split(deepcopy=True)\n",
    "\n",
    "#plot(w_)\n",
    "fe.plot(X.sub(1))\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"deflection [m]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}