Analytic FEM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analytic Finite Element Method\n",
    "\n",
    "This is a FEM implementation in a CAS and is intended to be useful\n",
    "for understanding the implementation of the method.\n",
    "\n",
    "The code solves the Navier-Cauchy equations for elasticity in 2D."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IPython console for SymPy 1.2 (Python 3.6.8-64-bit) (ground types: gmpy)\n",
      "\n",
      "These commands were executed:\n",
      ">>> from __future__ import division\n",
      ">>> from sympy import *\n",
      ">>> x, y, z, t = symbols('x y z t')\n",
      ">>> k, m, n = symbols('k m n', integer=True)\n",
      ">>> f, g, h = symbols('f g h', cls=Function)\n",
      ">>> init_printing()\n",
      "\n",
      "Documentation can be found at http://docs.sympy.org/1.2/\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from sympy import *\n",
    "init_session()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def const_model(E, nu, const=\"plane_stress\"):\n",
    "    if const == \"plane_stress\":\n",
    "        fac = E/(1 - nu**2)\n",
    "        C = fac*Matrix([\n",
    "            [1, nu, 0],\n",
    "            [nu, 1, 0],\n",
    "            [0, 0, (1 - nu)/2]])\n",
    "    elif const == \"plane_strain\":\n",
    "        fac = E*(1 - nu)/((1 - 2*nu)*(1 + nu))\n",
    "        C = fac*Matrix([\n",
    "            [1, nu/(1 - nu), 0],\n",
    "            [nu/(1 - nu), 1, 0],\n",
    "            [0, 0, (1 - 2*nu)/(2*(1 - nu))]])\n",
    "    return C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "r, s = symbols(\"r s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolation functions and matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIMAAAByCAMAAABHjZLoAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRCKJ3btm7818bGRdMYwAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAWFSURBVGgF7ZvrkqMgEIXxOrsxiTO7vv+zLjQ0t+4WcMymakp/JAT7NEdE8QuJUtm2ZJ/3PnZh50GZTdBtsPUu3T2k3St1X2bvdMOYQ7KHbVqpbhtGvU02W191PN3QPyB+cOHHZKtpuN+Mh6hHu/3jwc7SPWA9zLbmoAwOY8099JEf7Obo/enLzoP6go44KOM92KPzLenC2q/riF1NPYyrCaYype7+/CpJBu3k/bDQZOttu6mHG3s02WROBiNTw011GI7v/hQqkIEFlXu4hWgboJZ5jOpCEc8FKKhMqedzmmeXRJAJHmCITU/cTJJ+hNj+S2+f5gU+owe4Rm+MrOs3M9j3ZLwH33t2N7x+RsOUOSDTPCOb9Fh17nWf+HRo3bqG6vxczHQ8LPrS8VuU7NNWruZiZmQP7XxA94LMZcivzQeeQt9scq4x2TI8NnuUeujpjcpWfTX5m6gkg0byflB4u/QW0gImC7Vwyz4q4z1Epy60E0rEwwS3h/iMh+BQEmUmhPSDwqQhQVzKJ5N5cHsPyoyaelBmQFdvox8+B2XUg51K/9erO1KmH6r74KzAkzyQMdfi7xwP0+aHRUvjLvYUD/P9/R7G+e0euun9Hu7q7R7W5e0eZv1E8+5+mAa9bYObug5cmtx8cSDN8qLrojRz4SOSfqhen1vvKK14AEHmQ8V71A5uZaS5E+mbMYVMFvZJHpBpQiSWctKUI1EB77ks2il4YPARnt5B6R6NLWkykVF6/WjtPyUyX6sLggcGH8PUiI/nQJpMZJw/ejRMZHGM4ME+4deQJmUBnb6WNK0T3oPDxwrS5ECznjT3PFikqCHNBD6wg6tJc9eDG0ll0nTIlgJqNWnueUB8LJMmRmIXmPd60gwe5il/EnP4WEGaDGiqetK0HhZ9j02+j7LVFh+Tk43XZkaaFDRVPWmGfmA8MNSJHqzOvAJpMpEhQpcEWRzDX5uclCRzdEfq4/RMIgqFkgcamc+jSJo0MjEhyaIgyUOZOj1pNoGm8rJggvXwct4M7ZsS6yENefmnkzwUBub+YZzj4eJN6OUfwpsf28f+oCns/S5v/sm/Iy20x+y+eBM65afzJnPmpaoW3mRyFO+T+dyb5YgQtj4yy1HyIANtjrD1kZmF0rwpAi1B2PrI3ELBAwVaCWHrI4mFggcKtGGSThG2PrLVAwBtDfniGmsb7Do7ZkzOzDMe7LZAW0G+uMZavayadMaiv9lj+MLGAOPUkK9yNNQIu1E/yB7sCCyTr3Lo2wi7NR4c0JbJ162xtsJujQcLtBXk69ZYW2G3xoNdOa0hXxvZCrtVHijQhvuDS2DJl1ljlSJRF94L8wVJRCoQN8kOUoGRoXFXKnggunxyRPKlS7NiZKuHIvqG21sJfUNkbmKvH16Bvnn75vOeBy7+FXUneSADsMWr8fDx63eLhIn9HvP+/bUzZzGtsVXXGit0y7XGCt1wMS90w09n3nzeY+8MtjJm3pIsAlTMKN6rZXxEqX3PsHNHlkWGNJIHER+D1JRy7JRleWSUR/BA8THS6KKAnYxMiIzTCR4oPsaieFUiwU5GFmbUJDJOJ3hgVk5rsBNlbdjJe+BWTiuwE2WN2Ml7SJDCdlsNdqKsETsFD3YkpSunZezE38g2YifvgVs5VWXsdLJW7OQ9cCunFdjpZK3YyXtgfqKLXzHY0YFXHLvg2oqdggcKmvEFHd8fsF5acEW3GOeWZsNHmS+oNFbR1VsHhVRGagg+ioxDIxMP+eyI2EllUmSUTjgX9gcEUVyh6GGyRJ3XGqvck+K5kCXcHjLyuCCp7hwP3+PNUzxcvAmn+OJN6IaLN6EbzuJN+BYwkIB0MxHrD/+m1/93cjZ/YBzDn4jEpqQdMW9KMXw9/HdyHNU/Y6hpxZEkDcUAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}\\frac{\\left(- r + 1\\right) \\left(- s + 1\\right)}{4}\\\\\\frac{\\left(r + 1\\right) \\left(- s + 1\\right)}{4}\\\\\\frac{\\left(r + 1\\right) \\left(s + 1\\right)}{4}\\\\\\frac{\\left(- r + 1\\right) \\left(s + 1\\right)}{4}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡(-r + 1)⋅(-s + 1)⎤\n",
       "⎢─────────────────⎥\n",
       "⎢        4        ⎥\n",
       "⎢                 ⎥\n",
       "⎢(r + 1)⋅(-s + 1) ⎥\n",
       "⎢──────────────── ⎥\n",
       "⎢       4         ⎥\n",
       "⎢                 ⎥\n",
       "⎢ (r + 1)⋅(s + 1) ⎥\n",
       "⎢ ─────────────── ⎥\n",
       "⎢        4        ⎥\n",
       "⎢                 ⎥\n",
       "⎢(-r + 1)⋅(s + 1) ⎥\n",
       "⎢──────────────── ⎥\n",
       "⎣       4         ⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = S(1)/4 *Matrix([\n",
    "[(1 - r)*(1 - s)],\n",
    "[(1 + r)*(1 - s)],\n",
    "[(1 + r)*(1 + s)],\n",
    "[(1 - r)*(1 + s)]])\n",
    "\n",
    "display(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interpolation matrix is a matrix with the interpolation\n",
    "functions reordered for each degree of freedom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAADkCAMAAABXLJq9AAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRCKJ3btm7818bGRdMYwAAAAJcEhZcwAADsQAAA7EAZUrDhsAAA3pSURBVHgB7Z3tmqsqEoXtTifnTH9l75nc/7WOiAhFAQJrkdh5Oj86qFVF1YoYzdvoNMnXh1wsL71smzvdNv8DNj4rc3r5awzP76t5p1tlZ/cye7ktr9elv9e6T/fl9Pq12J+sfafbvUrc6efLKjBNL7fT2/w6G/uXnU/X6mUsz1aKy7Km083EOcLraup/vRkp/Ih/9c1kjt/b2lWK6a/ZLTrdtmBHaFwjKexHHWZ2fb1e37b9X0vxdp2ttdu88nPd0+Zmxi3s5+HtSIoPXdP1/fY+fbmjo67pPI+QhNs0nd6nF2fu3rdxNRm3Y70iKd59zmueH5e3cJ1vuwFiXLTb7P39fb5c1ihpt2NLsXxW52/3MqW8vtmUX//Orz/mz7LCSWG+VN/tRyz9Xl5v5iBUcGNKcTl9Xj93DnTl/qK9YtujA68/YQeJj3eWIeV2ng+nTsXEsWJVL+gHa/6dvwAvX3VnAumeIiku+ljxMR9Y/SuQ4o9de52/fhNu09es4MmpmHbzYeHWdcnmcznt6wwWSTF9udG9xZPHAVfTx+nrZj/z+eg4f4Mot+k6f/G4g+22V0RuWydw49uO7JvOozp0LIU90cq7Oym8xfJBdLr5KGjrtkjxflvOE/uCxVJsn18mnJLibE4rwhOHpGPOLWncsfJyOxmv99t6iO8IMSkp1tJyseLj0mVJYT5dsIq0uuXsW9d/3JYLhhf71upt7ZUUkzny17/e3ODsdKvvqWj5YfcKohT2Au1ef4vFtW0cMkDaUjiMtT1snqmHza7i1FGxKwri9L2kcGV+mXalc0Yy6OpROdlTrBPzFEv1UbHi8vl4KaYvc+ItrhEqMhcm+htEbK5aeLscQIrL/OvI3g9I5WoIUrycjyBFucyarQQpPqdfKazU149fKawSl/ms/3evWLQ4n+bX7VS+BKkZqg+3IRwr5l95D/ANgivJkOL6fXsFfifIFdFyWeh+LZtjdbpN+iI9l9jd1+9gOpdPhG6nTrc5HGOvcFlR3zcIVYwao9up0810clQpdiCshQ2mAIFu95Bvzs0EOqwUO+fQ/lLY4ZgF3e6x25zboaXQEEKwW13Tgm6T7LYW3R50gCQgrGC3WoqFwSbc6tHtQaWQ8MXsv5LdaikWD+3WgG5xKXBYuYxU+cdiRMlgV+pYYrApt2p0i0uBw0qpwrKUgrDh7zKJvcJ8OyTc6tEtLAUBViakSEBYwW4DKQJ0m2K39ejWSHE5O5iRSGpvFQFWprrQEFYcB5wUMYPVbvXo9mO+eAj/FyuVVnEdAVam4u9AWCeFd7U/8Ha6LWHQAcJgMb6eoKWLDTZqsutAZaebiY1KwYCVYY1b29W2rRCN+PLTodsddpt1M8FxKRZ+DMFKUeS20ARhN3Tbhny9G0GKMQPkXtTW/RPvoj+6V0wEWLntCI9twFIQYCVNgfIxc68bWAoCrNzLsXY7iG5hKQiwsrbUHTsU3eJS4LByp8TazSi6xaWozXS0HYxun0cKGN0aKf65/TP6IxsfH0e3/43mg4zPeUwPBHT7LAOEgG6fRQqzr4Ho9omkQNFtVor4ejYe4gxgG8d87HJOigKF5QHbx5Ye956RIk9hmcA2Tuaxy2kpEvA2R14TpmFFObfQ5iDttBQJeOsvgHnA9iAarGmkpbDwdjiw/QFSrBR2OLA9nhQX8WvnnJ+lL+OB7bGk+Jj/7VIhITcHdDSwZUqBU+zkscJR2NHAlikFTrGTUqzwdjiwJUpBoNhJKdYZtMOBLVEKAsVOS5GgsP68whWAA1sXifBOoNhpKTSe1Wsc1NQahZWprc4tNMLbDEiXkUInHF+pMoAtLoGLwKDYGSn2Kaw/F+nkvK4Kyjt/yq1Nazi8pVQvggwcIKKfn7BAoNi5AdJWvjo4trkTrAkUmyIFyG0JSkwEim2kQJEQym0ZUhAoNgMJodyWIgVOsQkDBOa2FCnwIAQpYG6LV0GJgEuBc1tKIXgQWAoCt8WroESApSBwW0oheBBYCpMCyG3xKigRGFKg3LZYSHxFHBt7drtnKTy927aaIcUWbECjnt0WLEVeMfLdNh5cinp2m7fcajUNhXz91mNLkQCyGQibsPRFzq2MW2hjpFBIKDR4aLue3SYsReb+2lki38AoiYSC7Y9t1rNbPUV3zrx2sq0t8tADpJ7dpuba1k+2/QFS1LNbgWzcntx4n2R8r8BhpUtdvdez29VSztGtnmzL2itwWKkkcCvq2a2zdJ7mvX6yrfWC9woCrAzzF+16dpuYott8n2RYCgKsFOWLBTuDVhwI3LdiNNlWz7Vtvk8yLAUBVorqxUI9u01YikhOQb9S3YQRlYLBYnx6qqUrUGtWpqnWy1hqs0ahqBQMWCmTFks64fj607FbbSkC5dwCI1yKUVNubZK7QHZjt7uWQdXTtLn5tagUIwfIndmtkQJCQgRY6T+Yh7ZgJESAlTQB1LGxKTI6QBiwsinhgjGIbmEpCLCyUF3LJhTd4lLgsLKl3oItim5xKQrJ3XUTjG6fRwoY3T6NFDi6fRYpCOj2WaQgoNtnkcIcoUF0uydFfEUXfyl49lhvGccgLaPo1khRQEIFEBmxx3pLUuX0MGUklAeRMXust6SXwApYHCAJEJlhj/WWrMT5cYpSJECkv/gT7LHekl8CK2JRigVEitmmwbQQJ4V4tGsbpWQVwYlTksKCSDHbNCFF+GjXn/Y0V6FhSYqFP8jZpgkpjJUjFY2UUmTy8IWiFPYY6R5KWrg/sIObjZTy4dWLBEpSrCAynG2a2itmwVbLVkopMgEXcIpdksKCSDHbNJQiuD/wiizrbwkM1p1wxyl2SQo73dQdB2z37ss0iSzrbwmcqAVbRaDYRSkSINJJ4RNf2GO9pXektggUuyhFMBxc3kqKLLLMWbpI3HcCxTZS5JGQBpHx9WceWeYsuRKs0RiQbgcJ7YLIjT3WW47QgkGxywNkRNZDYo6acrsmO4LeDtFhDsoYIE+yVzBuHM2RQn1bjPr0s3EJFJsiBchts/U1bCBNuVW3/WlIwZii3Laxu7T571NuN11wik0YIDC33cp5bIMgBcxtH6vA1jsuBc5tt2Qe2zBS5K9BKnIjcNuKXu5hsnMNsp8Cgdvud3IXC3yAzGmC3PYuhe53wpAC5bb7Wc5ql408xW6wjEIypIhCjljMw+mIYheechtbxnn+DCmycDqm2Pmn3CrLWAn4CVMq4IgVGk5nKHbiKbc5S53nj9grNJz2l8IO3eaecpuz/KFSaIytCxRPue2i2D9hr0hgbC1F+JTbPor9E6RIYGwtxWK04qs+io1LgcNKPWzlmhVOV980eeqj2EaKwr+lyZySSzisTIYNViYwdmKvMF8V1rKTYpf/LS3IJ9ckwMpcaLc+gbEDKQKKvT7ltpNiwwOEACtdydn3ZT6twNhOiohiW+DdSbFhKQiwMiuB26DhtJPCWUxT7qbJOUvv6VqoFAwW43LJv6t61ApHd9UGtcJZqt5QKRiwUiWlVqj04wtVR7H1U26zlqoTXIqxU27XhPfg9Eaxd++v7C1jLVAp7jBA7sVuUSkYsDL+eB60DEtBgJW00tUhsikyLAUBVjYlXDAG0S0sxe+UW//h4LDSx4JajCm3EBKCsmc6w+gWRkLMaqBYMLrFjxVQ/jxnHN0+ixQEdPssUhDQ7bNIYUYaiG6zUsRXdMVRzUCWxQ5qNqLoNidFHlLKrGIQmfeLLWWcAyxlpMhCSpmyApFZP2UpAx1hKS2FhpQy1xyI1H45SxnvEEtpKTSklMn6S0AespQ9PGApLUXipr9i5q2WAkeWDyhedmmkUEgoddNfMfNWS4EjS5nXA5aSSEggB5uUnHmrpcCR5QOKl10mB0jypr+jkaXMq3kJR7dJKVZIKdMJZ94m9goYWcrempdwdJuUInXTXzHzNpCChiybixcOBHSblGJ9tKvoSxw/nBRMZCl6a14goNu0FBpSytycFH4tjix9rJ4WAd2mpUjMtRX5KSkcyVMb1ApnKeLBCwwylZFiJ+H4spWBLCE5GOg2I8UuepSJexBZDzdlBHBp8DxTMLu7uo8aICN4rYg5QCXC3ZJzA2RAtkNDEtCtkeKff/8Dpqm+J8B47e4EdPu/f+d/8bsFv022ZzE/ieN26XGj+hxjnukhptzi6JZxrEC5LXXv6A9GkALmtv3ZUz0JUsDcllpQfzBcCpzb9mdP9YSlIHBbakH9wWApCNy2P3uqJyyFyQbkttSC+oMxpEC5bTH7+PeAorE/Vdxz85ZbQIYUW7ABjTyOlp1FcLrgFln6MAeXIoujfQWmFcPpvFtsGcQ5thQaRwepz80MnE64ZSzDcMeWohNjJ9z8tbME3oEWx5aiE2M7t7aZt4eWohNjO7fGmbeHlkJgKLsr12Bs59Y48xaXAue2wXiVzU6M3Xn/aFwKnNvK+oOlToy9urXOvIWlIHDboHjZ7MTYnfePhqUgcFtZf7iUeHatOxAsZu4rMsLY1q115i0sBYHbhsXLdifGTrg50Xz8n/+UW1/L3FIFrqxXrc9aBuHQvYLBbYN04mYnxtZu8ZXqBrx9j7gUY6fc7uFoX4lpbRh7122z9AFQKRjc1mcjW4KzjliQ3aFS/E659XoSuK0PBrb00bIlILxXELhtS74lWxDdwlL8Trn1nw7ObX0sqIWiW3yvgNInOsPo9nmkgNHt00iBo9tnkYKAbq0Uy6mc/32cOILvFQpDt1/2XNbMElpe7/dKe1Q//ej2ahUYldjd4xLQ7f8BlYfp0atZHuEAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}\\frac{\\left(- r + 1\\right) \\left(- s + 1\\right)}{4} & 0\\\\0 & \\frac{\\left(- r + 1\\right) \\left(- s + 1\\right)}{4}\\\\\\frac{\\left(r + 1\\right) \\left(- s + 1\\right)}{4} & 0\\\\0 & \\frac{\\left(r + 1\\right) \\left(- s + 1\\right)}{4}\\\\\\frac{\\left(r + 1\\right) \\left(s + 1\\right)}{4} & 0\\\\0 & \\frac{\\left(r + 1\\right) \\left(s + 1\\right)}{4}\\\\\\frac{\\left(- r + 1\\right) \\left(s + 1\\right)}{4} & 0\\\\0 & \\frac{\\left(- r + 1\\right) \\left(s + 1\\right)}{4}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡(-r + 1)⋅(-s + 1)                   ⎤\n",
       "⎢─────────────────          0        ⎥\n",
       "⎢        4                           ⎥\n",
       "⎢                                    ⎥\n",
       "⎢                   (-r + 1)⋅(-s + 1)⎥\n",
       "⎢        0          ─────────────────⎥\n",
       "⎢                           4        ⎥\n",
       "⎢                                    ⎥\n",
       "⎢(r + 1)⋅(-s + 1)                    ⎥\n",
       "⎢────────────────           0        ⎥\n",
       "⎢       4                            ⎥\n",
       "⎢                                    ⎥\n",
       "⎢                   (r + 1)⋅(-s + 1) ⎥\n",
       "⎢        0          ──────────────── ⎥\n",
       "⎢                          4         ⎥\n",
       "⎢                                    ⎥\n",
       "⎢ (r + 1)⋅(s + 1)                    ⎥\n",
       "⎢ ───────────────           0        ⎥\n",
       "⎢        4                           ⎥\n",
       "⎢                                    ⎥\n",
       "⎢                    (r + 1)⋅(s + 1) ⎥\n",
       "⎢        0           ─────────────── ⎥\n",
       "⎢                           4        ⎥\n",
       "⎢                                    ⎥\n",
       "⎢(-r + 1)⋅(s + 1)                    ⎥\n",
       "⎢────────────────           0        ⎥\n",
       "⎢       4                            ⎥\n",
       "⎢                                    ⎥\n",
       "⎢                   (-r + 1)⋅(s + 1) ⎥\n",
       "⎢        0          ──────────────── ⎥\n",
       "⎣                          4         ⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "H = zeros(2, 8)\n",
    "for cont in range(4):\n",
    "    H[0, 2*cont] = N[cont]\n",
    "    H[1, 2*cont + 1] = N[cont]\n",
    "\n",
    "display(H.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The local derivatives matrix is formed with the derivatives of the interpolation functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAToAAAA0CAMAAAAt4QSqAAAANlBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABHL6OuAAAAEXRSTlMA74lUMhB2mUS7qyLdZs3V90DuHGwAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAS8SURBVGgF7Zttc5swEIRlA7bjt5b//2eL4mKvdAu3IhdPOxO+RIi9R6sdHAdfnFLz0d+bS7SCKHAUx7j+Irg/nq+GGTERBY7imD0Z8G6cjr3RrUx03xNdSlHgKI7JYAbvc2a7tOv6vh+MamViBqxItl2KAkdxzC5m8DBF1k3R9U/F/nDonifLgxmwrKivhIE1ULvB1A7uMbrjKfUf9a7JebOzMLAIajaYNoCL6D4+ukF55TY7CwOLoGaDaQO4iK4/j+Pr1Utut79Tzc4YuDvYBVywCHI5ZmkGNqLyfQyj6y4pncmGDKLVGQVviU4FtRpMFGz2vRzddbrjjtJdtyPUlSkKptE5YBU0vfm1HRRMEAjGu+52ux1ORF9NXY7XUbo5n3UUbKPzwRrI5zydzQMKni8+f5ZgjO4pecvARrdx2TBQ4/o/0TUG9pKvRTfcr8/jrvwO/MT6Vft7Pq67zx/nl5fm0UaQ79BXTF4huvxY9nmsb2FWjesyuMorzOuMy4CTuKIGcRVy6rFYUckguhr4zef1jjcvFwZqdFBHpzyHTUuIMjBjKviOjQwQj6FRUJBRGU49IVa8ZFV03ag8iCVRBu5sBd2xlQHjc2gVDGRVNac+FytAVkY37KXoRBm4IxVsx0QGkDwkCgIiqopTn4oVKCujOwxSdKIM7JEKsuNEZADJQ6IgIKKqOPWpWIGyIrq+k6ITZeCOVQzTE3N1MFkpYQoLYqqSU5+JFYWsiG6fpOhEGdgTK3yZr8irairwp1YUYIzudpGiE2XgTKzwZb4ir6qpwJ9aUYIhumH6vEm460QZOBMrfJmvyKtqKvCnVlRgiK47Tsd4vCGUjEUZVIoVvsxX5FU1FfhTKyowRJdhF+kdtv8lydBdGFgCvcVgGd3tYzy7LbH+eP8tyDC4FAZWQG8yWEZXbHf5pPnz62VUeSUKHMUp3U1nCP6JzsSzNrEU3e08fcJ+tn+pGhYCzEUyEQZWQW8xiHfd7TSe0lXoTjQ7iwKrDt9iEKK7DIfXJyrk3nlNNToLA8ugtxiE6MQm7JRgo7NAMOvEkcf/GIMOuIgO/nHndYuRUbOzMDADOTsk/u0U4yYHjNFd1HYDNnKtDTsTBqYgtsPGFjbl0ugAjNGdtF91ZSPXBmVnwsAUZKKLMuiAMTq75/9ixuwwyrUD/oluOegvRCc1cs3KfpWvMFA+8W0tbAkMd13VoeVu00IfeUmd58PAHFTfHFy1waADhujW4P/ytXqHYV4dcB2d9iYb0cLmO/TXNwq6Q6Piy8GsrXDAVXTQoQWqGYoyqBMrfJlVsB1aFXihQ1LhgMvosENLF3hMijIgiBW+jCjIDokKzLAhq3DAZXTYoWUL/J0TZUAQK3wZUZAdEhWYYUNW4YCL6IoOLVvgMSfKACBW+DKm+LYWtgMuois6tLDzaijKoEqs8GW+Iq+qqcCfWlGAMbqyQ4vkYizKoEas8GW+Iq+qqcCfWlGCIbqqQ4toHIsyKBErfJmvyKtqKvCnVlRgiK7q0CIbx6IMSsQKX+Yr8qqaCvypFRU4Rwdf6pQaxGqnG92FgTWQptpucP5SZ37oe3yVWGkQT+uJMnAmVvgyX5FX1VTgT62YwY+vEv8BMl89VZQJB5gAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}\\frac{s}{4} - \\frac{1}{4} & - \\frac{s}{4} + \\frac{1}{4} & \\frac{s}{4} + \\frac{1}{4} & - \\frac{s}{4} - \\frac{1}{4}\\\\\\frac{r}{4} - \\frac{1}{4} & - \\frac{r}{4} - \\frac{1}{4} & \\frac{r}{4} + \\frac{1}{4} & - \\frac{r}{4} + \\frac{1}{4}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡s   1    s   1  s   1    s   1⎤\n",
       "⎢─ - ─  - ─ + ─  ─ + ─  - ─ - ─⎥\n",
       "⎢4   4    4   4  4   4    4   4⎥\n",
       "⎢                              ⎥\n",
       "⎢r   1    r   1  r   1    r   1⎥\n",
       "⎢─ - ─  - ─ - ─  ─ + ─  - ─ + ─⎥\n",
       "⎣4   4    4   4  4   4    4   4⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dHdr = zeros(2, 4)\n",
    "for cont in range(4):\n",
    "    dHdr[0, cont] = diff(N[cont], r)\n",
    "    dHdr[1, cont] = diff(N[cont], s)\n",
    "\n",
    "display(dHdr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gauss integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss_int2d(f, x, y):\n",
    "    acu = 0\n",
    "    pts = [-1/sqrt(3), 1/sqrt(3)]\n",
    "    w = [1, 1]\n",
    "    for i in range(2):\n",
    "        for j in range(2):\n",
    "            acu += f.subs({x: pts[i], y: pts[j]})*w[i]*w[j]\n",
    "    \n",
    "    return acu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Local matrices generation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Jacobian of the transformation is given by the product of the derivative\n",
    "matrix and the coordinates of each element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def jaco(dHdr, coord_el):\n",
    "    return simplify(dHdr * coord_el)\n",
    "\n",
    "\n",
    "def jaco_inv(dHdr, coord_el):\n",
    "    jac = jaco(dHdr, coord_el)\n",
    "    return Matrix([[jac[1, 1], -jac[0, 1]], [-jac[1, 0], jac[0, 0]]])/jac.det()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can re-arrange the derivatives of the interpolation function as a matrix that\n",
    "_translates_ from displacements to strains, commonly denoted by $B$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def B_mat(dHdr, coord_el):\n",
    "    dHdx = jaco_inv(dHdr, coord_el) * dHdr\n",
    "    B = zeros(3, 8)\n",
    "    for cont in range(4):\n",
    "        B[0, 2*cont] = dHdx[0, cont]\n",
    "        B[1, 2*cont + 1] = dHdx[1, cont]\n",
    "        B[2, 2*cont] = dHdx[1, cont]\n",
    "        B[2, 2*cont + 1] = dHdx[0, cont]\n",
    "\n",
    "    return simplify(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With these elements we can form the local stiffness and mass matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def local_mass(H, coord_el, rho):\n",
    "    det = jaco(dHdr, coord_el).det()\n",
    "    integrand = rho * det * expand(H.T * H) \n",
    "    mass_mat = zeros(8, 8)\n",
    "    for row in range(8):\n",
    "        for col in range(row, 8):\n",
    "            mass_mat[row, col] = gauss_int2d(integrand[row, col], r, s)\n",
    "            mass_mat[col, row] = mass_mat[row, col]\n",
    "\n",
    "    return mass_mat\n",
    "\n",
    "\n",
    "def local_stiff(dHdr, coord_el, C):\n",
    "    det = jaco(dHdr, coord_el).det()\n",
    "    B = B_mat(dHdr, coord_el)\n",
    "    integrand = det * expand(B.T * C * B)\n",
    "    stiff_mat = zeros(8, 8)\n",
    "    for row in range(8):\n",
    "        for col in range(row, 8):\n",
    "            stiff_mat[row, col] = gauss_int2d(integrand[row, col], r, s)\n",
    "            stiff_mat[col, row] = stiff_mat[row, col]\n",
    "\n",
    "    return stiff_mat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assembly process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def assembler(coords, elems, mat_props, const=\"plane_stress\"):\n",
    "    ncoords = coords.shape[0]\n",
    "    stiff_glob = zeros(2*ncoords, 2*ncoords)\n",
    "    mass_glob = zeros(2*ncoords, 2*ncoords)\n",
    "    for el_cont, elem in enumerate(elems):\n",
    "        E, nu, rho = mat_props[el_cont]\n",
    "        C = const_model(E, nu, const=const)\n",
    "        coord_el = coords[elem, :]\n",
    "        stiff_loc = local_stiff(dHdr, coord_el, C)\n",
    "        mass_loc = local_mass(H, coord_el, rho)\n",
    "        for row in range(4):\n",
    "            for col in range(4):\n",
    "                row_glob = elem[row]\n",
    "                col_glob = elem[col]\n",
    "                # Stiffness matrix\n",
    "                stiff_glob[2*row_glob, 2*col_glob] += stiff_loc[2*row, 2*col]\n",
    "                stiff_glob[2*row_glob, 2*col_glob + 1] += stiff_loc[2*row, 2*col + 1]\n",
    "                stiff_glob[2*row_glob + 1, 2*col_glob] += stiff_loc[2*row + 1, 2*col]\n",
    "                stiff_glob[2*row_glob + 1, 2*col_glob + 1] += stiff_loc[2*row + 1, 2*col + 1]\n",
    "                # Mass matrix\n",
    "                mass_glob[2*row_glob, 2*col_glob] += mass_loc[2*row, 2*col]\n",
    "                mass_glob[2*row_glob, 2*col_glob + 1] += mass_loc[2*row, 2*col + 1]\n",
    "                mass_glob[2*row_glob + 1, 2*col_glob] += mass_loc[2*row + 1, 2*col]\n",
    "                mass_glob[2*row_glob + 1, 2*col_glob + 1] += mass_loc[2*row + 1, 2*col + 1]\n",
    "\n",
    "    return stiff_glob, mass_glob"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Single element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = Matrix([\n",
    "    [-1, -1],\n",
    "    [1, -1],\n",
    "    [1, 1],\n",
    "    [-1, 1]])\n",
    "\n",
    "elems = [[0, 1, 2, 3]]\n",
    "\n",
    "mat_props = [[S(8)/3, S(1)/3, 1]] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "stiff, mass = assembler(coords, elems, mat_props, const=\"plane_strain\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAADKCAMAAACBkYe9AAAAdVBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA268pkAAAAJnRSTlMAMquZdlQQQO0wRCJmu++Jzd3T8a+538Vww9X34cfRpdfBXCB8bJjBZ8oAAAAJcEhZcwAADsQAAA7EAZUrDhsAABVrSURBVHgB7V17f9s4DpRjJ91N0tppd7vP2+edv/9HPEoy+IAAcYYKvf7d2X9UJEVgMEAsm+JUHobla/ewe3l9XY6vjBwP+4enlfPuqePjU3i5p4sTu1PRxTt9CTWHNSCWz+fz+QGnGmbuQ+n2HymTy+SHgHU+I6a7/cNbC0Kw6UmoPSzDcjdlo0j+8+MO/FuV5OyPw/B6lh5zfHgOr48Y3GNzMboSag5rEMu3uQbDsDvvD+H1mGXw+TnrwM22d8aI+wjiSexwSDKxL6HmsGIxXscKPIQ/5t15JzHL8TmUZ4/9sYrJMDx9WvhJJ1dbT+FtBb2aWfcl1BxWLMbEfryyGMV4CX+pr8hlPEvh4XTIelQT/uRvZt2XUHNYUDHGVO7OL1RGWz/AAxJc9nbWXQm1h1VYOu+M8U/16QxeyGPJXs75504crjaO++qUy4QidtRonNeXUHNY0DvjU3hTvJyJD43p8+LlzC1NJJsn9COjjF3MkWNfQn2LcQx1wP9ex2y8hfIdPhHly1L4Ca7h46fMjGn2JdQc1lBYOpep4XgcFw7Ea3c8HB/YD5mL/zfweviyfzs/NH5L6EioPSxl6RWDKMN96ntl4F6M98rkO/i5ejGeTm/xdWpdJFaI8xi8RSUE63QVxC7GfKMk/Gv5NMd4C3GDW+IzxXc88qa4BT4zhnNpaEu7GNrq3r9KBrxiNO1OwAtpRe14PMBfblsxwvdDHETiA8GakmVgOMVo2p14ZFaJEks4np6HF3T90IpBgUhsIFhTsiwMrxgNuxNPx7ZiHMIu1hO4emjFGBgQyRMK1r6VE4jnSXOKMcbD7k4cntqK8Ym4odWKMTAgUgwGjE2WieEXg92d2D02FuP8CO+dNGMMBIjkiQFjk2VjuMWgdyeO4TZvy62pl3O4Z/uMbae2YoSbnjiIJIoAo5NlY7jFYC9Try+NxXia7rtDt+ubMUJo47UQApE8cWBtlymFsVIMandi/ARue2fMKUJuvm/AwEGkFiQYlSwPwy4GvzvxuA+v8/glj36dRiPkj3YDxgCDSPg4GJ8sD8MuxvC2P+w/f2E+AoJK7MP5R4Fhjs8fd6dHUJy2+7npcyl8KP162P/rNyYsnNDby+7UtpUT9Hvfnv+OUTnF+O77w/HrD3EW0BhVYucH4ltq9Ln76esXTImy259+acMY9r/vj3/8GTGBBk7or98+f2nbypn1e58lGqcY/DqGl71JCEE1hX2XCgb4zOR9avUl1BqW0u85xRgJcF8Q2lRic8pwLvjM2XPxbz9CrWGNF5JMv+cXg1zHNKnELrnCueAzizJMnY6ENoSV6/fcYrDrmBaVmCQM54LPFN/x2JPQhrBy/Z5bDPYyNZLmZW9zqnAu+MxYhNTgLlOjHUqoPaxCv7dSDG4d06YSmzOFc8FnphrEVj9C7WEVeii7GA3rGFolFnPEfEdqZd2ZUGtYYS8nF0TZxWiQpPEqsVSNQsmVho0WPlMZ8xo7hlBzWEOh33OK0SBJo1Viki6l5JJh44jPXBj3JLQhrKHQ7znFWJC5D1whA/diXCHJKMTVi1FVcqGRr8zjMXiLFXjvVBVkLMbTo749q9VVnvs0zluILW6JzxTf8cib4hb4zBjOpaEtX8LtEeO/kWmze/8aGVi5TIEKLonyvZRc4s86bsIYSEKU7I31LewKaZ1fDFDBJV7fTcklDo3jJoyBJETJ3ljfQq7U77nFKNRVYrty5DcMkjMUawtGKRdL2G6Lkb2hBDSYwnCLwSi4BIO/EzdbMlitGAMDMobFyN5Y35IwheEVg1FwiWtyw0DMBgarFYMCmSIjZG8Mgch6bCgMrxiEgkvcsxsGYjcQWM0YDMgUGCN7IwhE0mNDYzjFUOqqwoXfabuEcFhtGAMHEigSsjfat2RQY9jFIBVc4pzbMLhYkVhNGJPKndTYTUquzto6hWEXA1dwSR0aNgwupjhWO8aAgwghXPbW4FtAlLTOLsY4eXpAAv6UtLRhgDxcbATIHo4GPowhYRCrvxgOCBICmxR1z2Edh2rr5mRlhEZ+a69L/ArDLcbrx1EuRjwl7fKf8o2HizlRpYejzVjOtGw4bkrgq78UDgwij3t73R8xbd34AKIxWYlQFrPZjPGXGGMxPpw/mCZhUKmsvGnFOLoHuUH2xqz+0HASB95itsUJOfH/OxRj5UahUlmlgFdaKJctsrcAj36tQsNJjHiL2ZYjZMTvXqZibLnKKg6uNVAuW2Rv4asn+tg3NJxEibeYbSlCVvz1YuQqqxTwSgvlskX2NuCrPzScxIi3mG0ZQmb81WIUKqsU8EqL4YKqxAw4421uzGJ0QGLOEBAbOeKEjPirxShUVgK5ekS5bJG9hQDQ1R8aTuLEW8y2HCEj/moxCpVVCnilhXJpl71Rqz80nMSIt5htYUJe/NViFCqrFPBKC1V0MSoxBZdWf+qE0UXDSaa8xWyLE3LirxajUFmlgN0Woehqlr0NcfXnhiEniHAuJryFYOGPr3PirxYjIt0b3TNwL0b3FOMAeTGqIqul2waTpZPKyDUwKiEYp3tENRbDvzelVVZGUGoIt8BnKoiwWymvxRk9IBMDS/DFW4hj3NKbWbk3JUD34zUykF+mTLxCZWXOWA6Cii5iU2KJgUrSmkBAAouomGQZGLVilCqrBbo5ACq64k1900ltsCcI6HsRIpMsC6NSDKWyWsBbA6iiy7mpb7lcjvUEQX3rqJhkmRiVYiiVlYY3+4yiy7hbZvpcDPYEYXzngTHJMjEqxVAqqxzaazOKLuumvue3GO8JwvgugiKSZWOsF0OrrApsp0Mousyb+o7bcrgnCOG7CIpJlo2xXgytsirA7Q6n6Gq8TPUE4XxnSSCS5WCsF+M2H2DGStKMnYMsh6pJaupyayVJy0+VbQ+jUgylsip9Wj1c0eXd1Le8qrGeILhvFVT4D/ZbHynnFoNXcrFqsXhTH16X5SoxUJIWQaAf7rzkl9TvDbYkbVGsYmDCyAmFs2Mxng5PxbzQidqvUmWlp2X9aHFRdGWnnKbc1McXf5lKDJak0do6kaTh+r1IAE6WYGSExiS9hKdumropfu+Rt5irhC/+cJWYUX82PFy/hxPQYSlC7mWKjb1FhpFiw75VcSqx5H1qsYQ4/R5GQIU0KEK3UQxw8UepxDRxthjBHtfvgQR0TIrQTRQDXfwxKjHNu+VZk7B+DyWgg1KEbqIYuHB2/IBr/P05/p3B6PfaLlNjcTJCN1IMbF3GqcTUnyFfDEa/hxFQIekfOf3ni0Es/mCVmCY99vligPo9goCOSxHyi4H+8E4EaNV+pXVZdOU0cJWY4YAPD9Xv4QR0WIqQUwxeycVbSGSy+JP+yrFd9tYSHqrfIwhobiUhpxja6N6/RgbuxbhGlkGMqxejh/hLc+UxeAuNCfSrIGMxliI2T2XlI/IW4gu3xGeK73jkTXELfGYM59LQlncRm87QP9hfuUwZKqu1QOFdCcsJiLUJA5W9xfg2StKin7VGgeEXw1JZrbiNN/VX5rinQKxNGPcnsbnpL06Yiq5ixtxp3zgI9iiI4G6WpImjlaPCcN8Zpspqxe94qvVuGYPVinF/EluleJfTtqLLtm3cOAg3R9lfO90sSbMJFKMKw3tn2CqrwpPutN7UZx6S1ozBgEzEtkvSdH6WfY3hFMNRWS39FSNtlxAOqw3j/iS2olBux1N0OQZtGwckyIi9WZLmEMiHFYb9zuCVXO039XGsdoz/jSexKZVVXlXVjjf1iYVZs+wtfB6jvz/6sHt5nX9xFpS9tUvS8MfWSfzGk9gMEZujslIFyLtyUx9fmLXL3pLGLo/AbCeZGCp7iwRoSRoue4vUSwxXxDZTUyork285yCzM+J3QiIWa9iUQw5kauOzN2QS2PzMiiFJZxfH1BvqNB82ogYaa9iVQBsbI3sz4a8XYHw7oUxNjZPDCzIwoulltoKZKJrbqM56ECUQLacCyNzP+SjGUykowV4/4wsyMaNV3PIma9iUQw5EGLHsz468UYwTJVFaCWTne0GWqLwGdB1z21lKMNtkYujAzI9IE7T5q2peAjg2XvZnxV94ZSmWlwZd9amFmRrT0aY2gpn0J6MhA2VswM+OvFEOprDS40Y+LP+OcHuJ1ZdEDatqXQAzn0kBlb6EYlkawUgz84WISlyz+pO8fW3RlF2+EaSkT86OJZ3AC0SQ2UNmbE3+tGBHn3uifgXsx+ucYRsiLURVZLb02mCydVEaugVEJwTjdI6qxGEsRm4BrlZWM+0fcAp+5QMNN8ZkCwlvwlh7GXcQmubyBY36ZMsMpVFbmjOXgLUnSiO2VRAQkkAwuLSZZBkatGMzDxSS2W5Kkxd0JCQ45ggQWrphkWRiVYiiV1QLeGkDVYszOxwKnJwjqWwfFJMvEqBSDebiYxHZzkjT0vmULAbEZj0yyzCRViqFUVjm012bUYs0bBz1BGN9FEohk2RjrxdAqqwLb6RDyN3znQ2P1BCF8F2ExybIx1otBPFxM4rpBSRp3meIICO1wJJLlYKwXA1dySVCkWgzd+RD387EnCOk7D0xJ0vJTZdvDqBQDfriYoN2YJI3aXpk44ASEcjzCyfIwKsV4/vWw/9dvEQ5oBNnbh/OPwMTh7WV3OnxaPHcMMQ1ysW/PfyMz3/aH/ecvDAhOQOM/f9yd0B8g3f18XgY1FsMUsc1I+9/3xz/+1Khr/VE1Nv62Zv3112+fvzy0PZVllot9rmMM331/OH79AZgYp+AEosmlsfvp6xdMSrPbn34xklQRsfELM0Y1Zu49aopWH5eL9SWgY8MJmTMrl6kRjfsywqjGzIg0QavPyMV6EtCx4YTMmfVikAszRjVmRqQJOn1YLob/8OiMxBDQseGEzJnVYrALM0Y1ZkakCTp9WC7Wk4CODSdkzqwWg32Xj/GhsjczIk3Q7uNysZ4EdGw4IXMmUAxuYcaoxsyINEG7j8vF8B8enZEYAjo2nJA5c70Y/JppYFRjZkSaoN0H5WKdCejYcELmzPVihIXZMHALM0Y1Ziq5NEG7j8rF+hLQseGEzJmVYjQoumDVmKPk0gTtPioX60lAR4YTcmZWiqHx7v2eGbgXo2d2Sd95MRp0WQ0mZIDh1tnpLb5OO9q8k0GPqMZi3EVsWcE8gVk2xWnilt7Mu4jNSe0/MZxfphx8Q23lzJyHNyq5Vn3Hkx1BmmRvEhiYLBujXgxLbSXQ1nGrksvyuRjrCNIke5MAwWQ5GNVimGorwTaOm5Vchs/FUE8QfgckhYcmy8GoFsNUWyX4RWuzkmvh0RjoDcJt4aQAmWQZGLVi2GqrBL9obVZyLTwaA51ByC2cGCCTLAujVgxbbRXhF43tSq6Fy+VAZxB2ByQGSCTLxKgUw1FbRfhFY7uSa+FyOdAdxLiELKNYjHDJMjDWi+GprRZxpIHNSq7kym/1BuG2cC5xkskyMNaL4amt/DRt/03NFd/x1Ga5WPS0aDTsgFx84MnyMNxixEedoQ8wuzwcTT1cbEG2GLB+U7OYoDqyVuJB8KekpR2QmAIVhe5mj6sDkxUxhNDsciyGIWKLjwrDfxx0//A2uSwfLqYDL/rzw9HSw9KKk0YnrZVoEOIpafQPkCYC6NPeZJclEZrYuiI2c1/QSFAa4i1mW1z25qyVUgh+C5e9JR8oIZxA8j23FCH3MoVGkvzzFrMtI3trkXpMKJzsbQ4MJcQRSAmbW9m3qlsoBvO0N2utpOk5fVz2dnEAF4MhoIPLCd1AMRjZm7lW0vycPix7E3u0GAwB8S3HgtANFGOMC5W9NV+mAgR4d1vS5DwTKp0uWjiBwix0busyxanGjLWSpmf3Gdnb7AF9Z3AEdHQZoRt4Z8CyN2+tpOnZfVD2lhmjxYAJZL6npibkF8N6VJj2VvRNXVYxw+7gsre4VrIdrY+isrfkBSWEE0i+55Yi5BTDUVlpZ1mft4jGsOxN1krRkmmgsjfxSRCCCYhvOSpCTjFk9v14zQzci3HNbFewrl6MHuIvzZHH4C00JtCvgozF+PDNt8qVp7JS07IubyHGuCU+U3zHI2+KW+AzYziXhrb8zzfjiutmRJM63P+vvneZKm+0gzmh17gXvx0VaSlyBkSsQEJNyTIwnGKoG+1iuH4EFVwLJx0VaQmLARErkFBTsiwMrxjH8EM/4eOEeaEKLu2zpyItYjEgYoQSUrsSYg4dCwynGKOj7A4W5JdRcOUOeyvSJiwGRIJjCLHJMjH8YuQ32sV07cgouAo/nRVpMxYBIsExhNhk2RhuMYob7WK7diQUXIWbzoq0GYsBkegIQnSybAy3GOxlilNwSTDh2F2RNmIRIBIaR6jtMqUwVoqR3WiXAP0jqeDKHfVWpE1YMIhERhKikuVh2MXQN9rF2j/iCq6Fj46KtIQFg4gJTohPlodhFyP7z/iokisATAouXC3W/rudQyYbE2LucSbAyN4i5W2SNDeicM/D/pFTpxhyoz2J2VZcX07NCi5cLRbXSrQibUiysVpYkQAMEi1w/d5F9hYJ1YLK4i9kb04xkjt071EscLXYhrUSIxtjCXBCBKE9HnFCTvzvXgxOLdb2JWRgZGPXK8ZUEOh2lhP/uxcjRISrxVrXSszT0q5aDJCQE3+PYsBqsea1EiMbu2YxUEJO/B2KwajFGi9T4wUBlY1dsxjMStmIv0MxGLVY01ppGBjZ2HWLgRFy4u9QDFAt1r5WCr9UEZ5JBi4C7B/uHN9a7osv3+SKIOTEXy8GLWZD1WJKweXmxjjByMZQKVqC4S1mW5yQE3+lGISSK5JB1WKysIyGRAOWjfEEeAuJmyBkx18phuDcj9fIwL0Y18gyiJEXoyqyWvpsMFk6qYxcA6MSgnG6R1R5MZaQWmW1nKFHcAt8psYYcFN8poDwFrylh7FeDMG5H6+SAaAY0K2vLFhGLcb6jjAdQd5LkhZjNRo2Rr0YoJIrIjJqMdb3NUDwXYkYTWqAhByMajEKlVVCdVuMWoz1HUF7guC7EjGc2EAJORjVYjBKrjEoRi3G+o6ke4O03r9kCBkYtWIwSq4pV4RajPYdi9EZBNyViOFIgyFkYdSKQSi5ppAYtRjrWziHe4T7YXienxsTB51GAwi6K7FAJLBMjEoxlMpqAb8YINRitO8I1h3EuIREcL/BETIw1otBKrnGOGG1WIPvmIfeINiuRAxnbpCEDIz1YuBKrhgYrBZr8H0NEGJXIoYzN3BCHsZcjGl9/qCcX7roJs5lOqMWgzeIdGQ9QfBdCR1V6IPJWmK8zbdIgprjML2eDefhv8x8NH730px5GYTVYg2+I25HEGJXIoYjDTRZS4zXuQbDfwFgvbKV7OQF5AAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}\\frac{5}{3} & \\frac{3}{4} & - \\frac{7}{6} & \\frac{1}{4} & - \\frac{5}{6} & - \\frac{3}{4} & \\frac{1}{3} & - \\frac{1}{4}\\\\\\frac{3}{4} & \\frac{5}{3} & - \\frac{1}{4} & \\frac{1}{3} & - \\frac{3}{4} & - \\frac{5}{6} & \\frac{1}{4} & - \\frac{7}{6}\\\\- \\frac{7}{6} & - \\frac{1}{4} & \\frac{5}{3} & - \\frac{3}{4} & \\frac{1}{3} & \\frac{1}{4} & - \\frac{5}{6} & \\frac{3}{4}\\\\\\frac{1}{4} & \\frac{1}{3} & - \\frac{3}{4} & \\frac{5}{3} & - \\frac{1}{4} & - \\frac{7}{6} & \\frac{3}{4} & - \\frac{5}{6}\\\\- \\frac{5}{6} & - \\frac{3}{4} & \\frac{1}{3} & - \\frac{1}{4} & \\frac{5}{3} & \\frac{3}{4} & - \\frac{7}{6} & \\frac{1}{4}\\\\- \\frac{3}{4} & - \\frac{5}{6} & \\frac{1}{4} & - \\frac{7}{6} & \\frac{3}{4} & \\frac{5}{3} & - \\frac{1}{4} & \\frac{1}{3}\\\\\\frac{1}{3} & \\frac{1}{4} & - \\frac{5}{6} & \\frac{3}{4} & - \\frac{7}{6} & - \\frac{1}{4} & \\frac{5}{3} & - \\frac{3}{4}\\\\- \\frac{1}{4} & - \\frac{7}{6} & \\frac{3}{4} & - \\frac{5}{6} & \\frac{1}{4} & \\frac{1}{3} & - \\frac{3}{4} & \\frac{5}{3}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡5/3   3/4   -7/6  1/4   -5/6  -3/4  1/3   -1/4⎤\n",
       "⎢                                              ⎥\n",
       "⎢3/4   5/3   -1/4  1/3   -3/4  -5/6  1/4   -7/6⎥\n",
       "⎢                                              ⎥\n",
       "⎢-7/6  -1/4  5/3   -3/4  1/3   1/4   -5/6  3/4 ⎥\n",
       "⎢                                              ⎥\n",
       "⎢1/4   1/3   -3/4  5/3   -1/4  -7/6  3/4   -5/6⎥\n",
       "⎢                                              ⎥\n",
       "⎢-5/6  -3/4  1/3   -1/4  5/3   3/4   -7/6  1/4 ⎥\n",
       "⎢                                              ⎥\n",
       "⎢-3/4  -5/6  1/4   -7/6  3/4   5/3   -1/4  1/3 ⎥\n",
       "⎢                                              ⎥\n",
       "⎢1/3   1/4   -5/6  3/4   -7/6  -1/4  5/3   -3/4⎥\n",
       "⎢                                              ⎥\n",
       "⎣-1/4  -7/6  3/4   -5/6  1/4   1/3   -3/4  5/3 ⎦"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAADKCAMAAACFSJXFAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRLuJZiLvzd18bLCbjx4AAAAJcEhZcwAADsQAAA7EAZUrDhsAAApDSURBVHgB7Z3rmqo4EEXjBc8Zb93O8P7POgnKpTYhVdWJDalPf7QQZZssElRWB52D2wXW1avXw+VwUm8FG+y+oUC/ml+Pps1sx/Xu3D6zJbvD/qFvPN3iJ/XYtd1t/0w6XXNZPDyLXev/ZN2abBaKejyeCEK9D0d/a55VP56yWZydu7f+T9atAAt5PW4BwL4NLHZjrXdNNosQlj3QXD4LbT1uyOLq0ixOh+vtOmE3UiRLjyNZhRVJCMdCkuGcph7I4nZnWHz7oXQKIzF5ux6SD0tCOBaSDKeqB7A4+d2Z7Be3r9DIK/MucUujEIUwLEQZunoAi+bgb+3htrhXL927DXM0OF+dOyfGkSiEYSHJUNYDWAQG99T7SNuxOLevN50os9130zSXxIcUSYhruh4YfYFQKMjQ1mPO4nZp94stPbVd7z+3qUPjV/dmvdgKPwj5kPvh0e4TLyLJ0NZjzmK5Df6Re+u7f3gX7u6ST11+sETIOzLULLpdmssiP+T+7FtZFcEMJQtJ11zuEK9HSoS8IwNYPD+W+79LDXoesprUsZPNeB33UiGrZACLJQRD+aX7Qn9LvdUMz11cKBHyhowZi+ZwTR0Xnx9xDunPWkyGk4SskIEsjv7zw7E7tC3s1Uf4DP6V+CDlt+cyHB+yRgayCF+1kx+2Tlf/5TaJwrEZjg9ZIwNY3Luv73nnHqrNABanJ4vFz50LA4cUV5sBLNyXx7Brl7+bkVYvrNSagSwa/555zBsjrtYMZOF2x2P2edtKM5BF+KZ9zhQktWYgi4v/nnxJv2V6WmlYggzW4ggyHGeTBBmkHsjifDweuJOZ3EluPoO3OHwGb5P4DFoPZLHw1kCKrdqkn7CwapN+wMKsTQIWEgHD2SQnCdFZHDJEhxXmTLm6HsBCIGBYm+QEIV6xpL4MO1EGx0JbD8pCIGBYm/Q6O8HoJJ3FGboCWWBYCBrj3LQelIVAwLA2yQlCnNLiEATDCsNCXQ/KQiBgQk2SJzjeYXGG5pOF0jaJsBCeXE7aJJEJ0locguC1Ut4mERYlBEwRnVSiIvoMYJEvcTyL/JB1MggL4RiJ9dixrETIOhmBxanpnbhlEzTurWEJjNS9O6fXf0kvIWBciZBVMsgYEUkct4LFGfbkdKF4PSgLgcSp1gRNMfbL1EgBC17iVGuC+vZP76mRAhbTJ8aXqzVBkeZAW9QsqjVBERbQFjULwzZJz6JWExTpF2C19Czs2iQ9i1pNUKxf0LbMWBB7EtteYGD8ZjXaJGRB7UmMBW9g/FbcP0mzRqpEPbQ2KbD40/4ZGq2YizNsM1+o0yb9C/NHHvK5OHMEQ0mdNgnHSGgO08GHFi8uVGqTYiySc3EkIqhSmxRhkbY4AgFTq02as5jak/koEAiYam3SjAVjcQQCplqbhCy4uTiWbRKyYCyO8AR1nTYJWcyPEKREL2DI5s+VEiHvyFCzyBdBm7VJShbCMRLpDJOiEiHvyAgsJt9HwJ5MGvBatGyT8PvIvPW0pITE2apNmo0RxsBIpgTVapOQBbUntE90a/yUoGptErKg9iTCwrBNAhZgTyIs+KJqM4AF2BO+4ZFnVJsBLAyboMhegyJkYdcEQcMjq8jCrgmKNB6KkAW1J/Bk4WqtGchCYILM2qTA4nQMe/J5401QCYvjX4s52f479aA26e7nl5FrjfVQFu/t2iQcI4sIhgfs2iQ9iwCF6eADt8WFTdokYCExQb6BSZskmc+zSZsELAQmyKNI2yTBnKBt2iTKQmCCPIq0TRLMK9qoTaIsBCbIT/HOvsLcRm0SZSExQZxNkswr8p1rg3OTCAvRyWXGJonmFflxlrrSnSijRD3o3CTCooSA+cwr8p2+v60zJ6h/9fFeXw/SL0RjZHy1+FKJkHUyAovRFVk2QZE9B2aMuiLLJijCAorIGPnMK5p+ZzdsgqATdKvUjNF+Ibg43Gde0Qi1WhM0NmFYgrZAvxietrhQrQmKtAjaomZh2CbpWdi1SXoWdm2SnkWtJihyvHC0LXMWpa/sFquDwEj5zX57blJgMXVF1J7E2vE7Fse/MnOyvXw95q6IuSJPDA+U1WqT5mOkAItK5ya9g0XoJkwHh54UWV3BJqlZGLZJahaGbZKWhWWbpGVh2SZFWCR/J8iyTUIW1J7Mju+iE9QlLM4aNglZzFpPC0zbJDWLz7yivneIxkj/5KX7EiHvyAj9YnRFDuzJvDWWbRJ1RfO2Y4llmzQ7XnzmFQ27n9qToXhcMGyTsF985hX116z0/0gSFvN+i6XaDOgXYE/GoaFYqjYDWBg2QfzeRBZ2TZCehV0TpGdB7Qm/fewZtWbgGJFYHKs2KbCYuiLewNi1SXNXFOv1tKyAQfE/D5X9K2r+n4jyfsbVz/0g9cAxQlsdXyvAYpM2aR0WgfH2bBKwEJkgrl+IQjY4NwlYiEwQx0IUssG5SZSFzAQxLEQhW5ybRFmITJBjWEhCNjk3ibKQmCB/1EvaJMm8om3OTSIsRCeXGZv0a3OCytskwsK0CYp/VpqWAgvDJmja6vgyYSEaI/GcsbREyDoZgcXoiiyboHFvDUtgxqgrsmyCBgKLC2SMfOYVfeYVvXoK7RefeUWqa8RUa4IihwloC/SLyAZQVK0JgnaEVWiLmoVhm6RnYdcm6VnYtUl6FrWaoMjxgptXVOQKc3XapNAvpq6oxBXmarVJ6IqoPYn1K0kZcxZQEFGmHjqbhMeLMleYK8BiBZuELMIuy7Y43NlhQb8oUg/l3KQYi+QV5kQiiOsXopBft0kRFmmLIxJBHAtRSLoeb7jS3ZxF2uKIRBA3RkQh6Xq840p3MxaMxZGIIH/EeSQPCZIQph6Oz1Bf6Q5ZcBbHsk1CFsycINEJ6lptErJI9u1wRUF/kcbwPzXdHfPkpYdLhLwjQ83CsE1SshCNkaX+0JeXCHlHRmDx5+8/fTXBnvTF471lm/TfX931wS3bpNkY+cwrGgbBZ17RMJeGv8Lc5/eKhn4jWAADI9hi/pRVMuB4AfZkXklBSbUZwMKwCeL3IrKwa4L0LOyaID0LuyZIz0Iwr8isTcLjBT+vyK5NQhZ8TypjcZizgHw1YE6QYIPoU6Y26QcsilicAiyK1ON4Gudp6VkEvDZtErCQSZx8i8OdKZf85pHfJ0mrJckgv5sELEQSJ//3ilhrIDBBHkVhm0RZiCRO/u8VhVGWNiiiipS2SZQFL2B8M/ItDs9CUpH8eoBNoiwkJoizSZJ5Rb5fbHBuEmEhOrnM2CTRvKJt2iTCooSAKaKTSlREnwEsDJsgf5BiboSFaIwwgaIxss0MwuJ12GvaZrmyrE0ShGw0g7KwbIKWd2//CGXx/Ihz+O4fjd4zNkk0OWmTGZSFM3yFudh+pWYMWBg2QTEW9Ip7wCK2AS1bxeLQKnRrb6iHmkW1JijCE9qiZmHYJulZ2LVJehZ2bZKehV2b9GTRfSjeR44ukSKBTYpsRYs2l/F4fi0I84q6mz/HLrnxNolP2VzG7YnA/Q+jYRPD3NYiagAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}\\frac{4}{9} & 0 & \\frac{2}{9} & 0 & \\frac{1}{9} & 0 & \\frac{2}{9} & 0\\\\0 & \\frac{4}{9} & 0 & \\frac{2}{9} & 0 & \\frac{1}{9} & 0 & \\frac{2}{9}\\\\\\frac{2}{9} & 0 & \\frac{4}{9} & 0 & \\frac{2}{9} & 0 & \\frac{1}{9} & 0\\\\0 & \\frac{2}{9} & 0 & \\frac{4}{9} & 0 & \\frac{2}{9} & 0 & \\frac{1}{9}\\\\\\frac{1}{9} & 0 & \\frac{2}{9} & 0 & \\frac{4}{9} & 0 & \\frac{2}{9} & 0\\\\0 & \\frac{1}{9} & 0 & \\frac{2}{9} & 0 & \\frac{4}{9} & 0 & \\frac{2}{9}\\\\\\frac{2}{9} & 0 & \\frac{1}{9} & 0 & \\frac{2}{9} & 0 & \\frac{4}{9} & 0\\\\0 & \\frac{2}{9} & 0 & \\frac{1}{9} & 0 & \\frac{2}{9} & 0 & \\frac{4}{9}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡4/9   0   2/9   0   1/9   0   2/9   0 ⎤\n",
       "⎢                                      ⎥\n",
       "⎢ 0   4/9   0   2/9   0   1/9   0   2/9⎥\n",
       "⎢                                      ⎥\n",
       "⎢2/9   0   4/9   0   2/9   0   1/9   0 ⎥\n",
       "⎢                                      ⎥\n",
       "⎢ 0   2/9   0   4/9   0   2/9   0   1/9⎥\n",
       "⎢                                      ⎥\n",
       "⎢1/9   0   2/9   0   4/9   0   2/9   0 ⎥\n",
       "⎢                                      ⎥\n",
       "⎢ 0   1/9   0   2/9   0   4/9   0   2/9⎥\n",
       "⎢                                      ⎥\n",
       "⎢2/9   0   1/9   0   2/9   0   4/9   0 ⎥\n",
       "⎢                                      ⎥\n",
       "⎣ 0   2/9   0   1/9   0   2/9   0   4/9⎦"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: $2 \\times 2$ mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = Matrix([\n",
    "        [-1, -1],\n",
    "        [0, -1],\n",
    "        [1, -1],\n",
    "        [-1,  0],\n",
    "        [0,  0],\n",
    "        [1,  0],\n",
    "        [-1,  1],\n",
    "        [0,  1],\n",
    "        [1,  1]])\n",
    "\n",
    "elems = [[0, 1, 4, 3],\n",
    "        [1, 2, 5, 4],\n",
    "        [3, 4, 7, 6],\n",
    "        [4, 5, 8, 7]]\n",
    "\n",
    "mat_props = [[16, S(1)/3, 1]]*4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "stiff, _ = assembler(coords, elems, mat_props)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "stiff_exact = Matrix([\n",
    "    [8, 3, -5, 0, 0, 0, 1, 0, -4, -3, 0, 0, 0, 0, 0, 0, 0, 0], \n",
    "    [3, 8, 0, 1, 0, 0, 0, -5, -3, -4, 0, 0, 0, 0, 0, 0, 0, 0], \n",
    "    [-5, 0, 16, 0, -5, 0, -4, 3, 2, 0, -4, -3, 0, 0, 0, 0, 0, 0], \n",
    "    [0, 1, 0, 16, 0, 1, 3, -4, 0, -10, -3, -4, 0, 0, 0, 0, 0, 0], \n",
    "    [0, 0, -5, 0, 8, -3, 0, 0, -4, 3, 1, 0, 0, 0, 0, 0, 0, 0], \n",
    "    [0, 0, 0, 1, -3, 8, 0, 0, 3, -4, 0, -5, 0, 0, 0, 0, 0, 0], \n",
    "    [1, 0, -4, 3, 0, 0, 16, 0, -10, 0, 0, 0, 1, 0, -4, -3, 0, 0], \n",
    "    [0, -5, 3, -4, 0, 0, 0, 16, 0, 2, 0, 0, 0, -5, -3, -4, 0, 0], \n",
    "    [-4, -3, 2, 0, -4, 3, -10, 0, 32, 0, -10, 0, -4, 3, 2, 0, -4, -3], \n",
    "    [-3, -4, 0, -10, 3, -4, 0, 2, 0, 32, 0, 2, 3, -4, 0, -10, -3, -4], \n",
    "    [0, 0, -4, -3, 1, 0, 0, 0, -10, 0, 16, 0, 0, 0, -4, 3, 1, 0], \n",
    "    [0, 0, -3, -4, 0, -5, 0, 0, 0, 2, 0, 16, 0, 0, 3, -4, 0, -5], \n",
    "    [0, 0, 0, 0, 0, 0, 1, 0, -4, 3, 0, 0, 8, -3, -5, 0, 0, 0], \n",
    "    [0, 0, 0, 0, 0, 0, 0, -5, 3, -4, 0, 0, -3, 8, 0, 1, 0, 0], \n",
    "    [0, 0, 0, 0, 0, 0, -4, -3, 2, 0, -4, 3, -5, 0, 16, 0, -5, 0], \n",
    "    [0, 0, 0, 0, 0, 0, -3, -4, 0, -10, 3, -4, 0, 1, 0, 16, 0, 1], \n",
    "    [0, 0, 0, 0, 0, 0, 0, 0, -4, -3, 1, 0, 0, 0, -5, 0, 8, 3], \n",
    "    [0, 0, 0, 0, 0, 0, 0, 0, -3, -4, 0, -5, 0, 0, 0, 1, 3, 8]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAHBCAMAAAC2QPTkAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRIlmzd0i77t8bBwggJIAAAAJcEhZcwAADsQAAA7EAZUrDhsAABVsSURBVHgB7Z3rehs3EgWp2FaylmLHWb7/s65mAFDxuA/Rmf560TYrP0IRtwEKhwOaVpmXp+v+328X/oPAvyDwpeXmcnm6fvj49t+nf9GXphC4vG6p+e26BegJHBA4R+CVAJ0DR69GQAXo+cPL64u+Ma2svqy8+Mprl1y4CtCfb++Inr98Va+zldWXlRdfee2SCxcBev1ri87LnyJAK6svKy++8to1Fy4C9G3/M/2n67OdoJXVl5UXX3ntmgsXAbruAfp8FX+0X1l9WXnxldeuuXA7QM/XD9ut5/P1o3kHWll9WXnxldcuunA7QF+vL1tyntrDDyFaWX1ZefGV1y66cBWg/Q6kA7Su+o3juouvvHbRhdsBWnmrnly76J18v0tPph6rLrpwO0D9/dqn+2+i11RXnlt7j53EpebCRYC+fdteUq/yj/ELqy+F55Y6tZoLFwFqH5h9uP9B4prq/nnamouv5FJz4SJAly/bX2X8Jf8ybGV15bmlcim5cBWg55e3X/WQ+bmsrF56cRZ++FBHBejQjKcQsAkQIJsLpU4CBMgJimY2AQJkc6HUSYAAOUHRzCZAgGwulDoJECAnKJrZBAiQzYVSJwEC5ARFM5sAAbK5UOokQICcoGhmE1ABWinQTa5d0q/rdCdTj1WXXLgK0EqBbnLtkn5dD9Bk6rHqkgsXAVop0E2uXdOvawGaTD1WXXPhIkArBbrJtWv6dS1Ak6nHqmsuXARopUA3uXZNv64FaDL1WHXNhdsBivkDqb2Lygl7gB5x4XaAVgp0k2sX9ev2AE2mHqsuunAVoLLyXlG/rgcoEVvRhW8Bev50/Ec4Uu/FscE5wmr9gwVf3+wL499IbG/3kgy52OA1/bp/volOwlZz4fYRVtNha5tUeW6IhX2PajpsfXLt4zjEwrFX43ENF3EHKumwDVKp9l7hwUtuigrQSoFucm3EwvFKOjxOuOVUqwAd5sZTCNgECJDNhVInAQLkBEUzmwABsrlQ6iRAgJygaGYTIEA2F0qdBAiQExTNbAJbgH6//m5XUgqBGYG/+b6wGSLq7xHgCLtHh7opAQI0RUSDewRUgGIKXGrvkn5dZ/x4C1cBiilwqb1L+nU9QI+3cBGgmAKX2rumX9cC9IALFwGKKXCpvWv6dS1AD7hwEaCYApfau6Zf1wL0gAu3AxQTJ1J7Y2XUsjLsAMUUuNTeRf26/Q70iAtXAUo05GJf+1fUr+sBSsRWdOFbgJ4/Iha29zCVz8eic/v69r3MiIU9Pm8PMekxtXfNudlHWGV5r/LcEAvHS3GNpNauPrl2ZelxMvVYdc2FiztQSYdtpLuw+5c6tZKbogKUY6H1CMQGRywcr6TDYwzryd4qQIe58RQCNgECZHOh1EmAADlB0cwmQIBsLpQ6CRAgJyia2QQIkM2FUicBAuQERTObwBYgxEKbDaUOAoiFDkg00QQ4wjQbahwECJADEk00ARWgVEMuNjhiodjOGNaTvVWAUg252OCIhSJAMawne4sAtd9ceflTTHVlNWJhqU0RAUo15GKDIxYef4O9ByqG9WxvEaBUQy42OGLh2/fjWP/FsJ7tbQcoVQ2MDV5UTth3NLaySe+iC7cDlGrIxQZHLLy+WPefVVxUgBINOcRClYC70BELx+tmcqueVBe9kz/sEYZYOHLdHtt7yaTvHIwNjlg4tiqm3yEWqj/Gf9sAv17/v9X2e6CaDlsPYEzP+3l719wUEaCSDtu4g6Xae4UHL7kpKkAnNbO+x6m9EQvHK+nwmEpdDa4CdJgbTyFgEyBANhdKnQQIkBMUzWwCBMjmQqmTAAFygqKZTYAA2VwodRIgQE5QNLMJbAFCLLTZUOoggFjogEQTTYAjTLOhxkGAADkg0UQTUAE6qZn1C6X2RiwU25lKXQ2uAnRSM+tLS+2NWCgClEpdDS4CtNIcnFwbsVDkZ8ItqVoE6Kxm1taW2huxUP3K4W8b/U/yNxJzqkWAzmpmLUCpvRELEQt3feXzuW/ew8o4x20iu5yutu9AMfcvtfcqgW6/t6aubDJ40YWrAN113GJqYKx3Ub+u5ysRW9GFbwHiGwvbW7e3/5++k+8jpPYuOjfEwlt49h9i7l9qb8TCsVWIhYPE948TLjWNSvs9UE2HreNun4h9uP+vp/2C1TU3RQSopMM2Xq+F3b/UqZXcFBUg5ZH1TVxZjVg4XkmHxyWbogJ0mBtPIWATIEA2F0qdBAiQExTNbAIEyOZCqZMAAXKCoplNgADZXCh1EiBATlA0swlsAUIstNlQ6iCAWOiARBNNgCNMs6HGQYAAOSDRRBNQAVIeWR9pZTViodjOJZuiAqQ8sj73ldWIhSJASzZFBCjJQmsrjw2OWCjyE8N6trcIUKoaGBscsRCxcHckP1/PCXKIhee4Jdme9h0oVS+IDV5UTtiPldjKJr2LLtwO0MRxW1ld1K/bA5TKpejCVYASDTnEQr6x8Bqx2ye36kl10Tv5wx5hiIXf/6E4VQ2MDY5YOLZqItBNqmv6dW1tk6nHqmsu3H4PVNNh6wFsH3n9gubgZGE1N0UEqKTDNu5gqfZe4cFLbooK0BJJrSdkcm3EwvFKOjxOuOVUqwAd5sZTCNgECJDNhVInAQLkBEUzmwABsrlQ6iRAgJygaGYTIEA2F0qdBAiQExTNbAJbgBALbTaUOgggFjog0UQT4AjTbKhxEFABWuIY9flOro0XJvZ1wi2nWgVoiWPUwUyujRcmAjThllMtAnTWEmpLS+2NFybyk0pdDi4CFFO3UnvjheGF4YWZN5GJuVVTiLPvQJNfbF9ZzS/V84Vz111reWoPP7wWJ3JVUT1qX8Zk6rHqogu370AxdSu19xvHss5a6tSKLtwO0MozanJtjrB6R5jxd2Exgym1d009qp3ED7hw8XdhMYMptXdNPaoF6AEXbh9hNRWk/n67faaFF9Zx3B7WcBEBKqkgDVSF1a3UqZXcFBWgHImoRyA2OF7YeCUdHmNYT/ZWATrMjacQsAkQIJsLpU4CBMgJimY2AQJkc6HUSYAAOUHRzCZAgGwulDoJECAnKJrZBLYAPX8Uv+Nmd6EUAu8ExD+y+d6AnyBwjwBH2D061E0JEKApIhrcI6AClGOh9ZnEBkcsFBsaw3qytwpQjoXWVx4bHLFQBCiG9WRvESDpkbW5r6xGLBT5WbMpIkCpamBscMRC8aFLDOvZ3iJAE8dtZXVNv67dFVK51Fy4HaCJGbGyGiujnpXxdH06HKsxBS61d1G/bgf4iAu370Cphlxs8KJ+XQ9QovNYdOF2gFaeUZNrc4TVO8IQC99P8FQ1MDZ4TaMSsfA9PNtPqWpgbPCac7OPMMTC72N1e7ZG3uuXX3lxeW0RoJIO29jGVHuv8OAlN0UF6KRm1vc4tTdi4XglHR5TqavBVYAOc+MpBGwCBMjmQqmTAAFygqKZTYAA2VwodRIgQE5QNLMJECCbC6VOAgTICYpmNoEtQIiFNhtKHQQQCx2QaKIJcIRpNtQ4CBAgBySaaAIqQCc1s36h1N6IhWI7U6mrwVWATmpmfWmpvRELRYBSqavBRYDWSGoNzOTaiIUiPxNuSdUiQGc1s7a21N6IhYiFfGOheROZaImIhZ3aRLuYVGNl1LMyEAtvN4RUNTA2eFGj0n4PFHP/UnsX9ev2CD7iwu0ATY6RldUcYfWOMMTC2xFWU97r04t5iUm9EQvfw7P9FHP/UnvXnJt9hCEWfh+r27P2adyab0usuSkiQCUdtrGNhd2/1KmV3BQVIOWR9U1cWY1YOF5Jh8clm6ICdJgbTyFgEyBANhdKnQQIkBMUzWwCBMjmQqmTAAFygqKZTYAA2VwodRIgQE5QNLMJbAFCLLTZUOoggFjogEQTTYAjTLOhxkGAADkg0UQTUAFSHlkfaWU1YqHYziWbogKkPLI+95XViIUiQEs2RQQoyUJrK48Njlgo8hPDera3CFCqGhgbHLEQsRCx0LyJIBZ2LKnSBlZGPSsDsfB2Q4i5f6m9EQvHLsX0O8TC68sg+d1jDOvp3vab6NRDKDY4R1i9Iwyx8P2lnKTftQvEBq8pPSIWvodn+ylVDYwNXnNu9hFW02HrW90+8loj9628ds1NEQEq6bCNe0WqvVd48JKbogK0RFLrCZlcG7FwvJIOjxNuOdUqQIe58RQCNgECZHOh1EmAADlB0cwmQIBsLpQ6CRAgJyia2QQIkM2FUicBAuQERTObwBYgxEKbDaUOAoiFDkg00QQ4wjQbahwECJADEk00ARWgJZJan+bk2oiFYjsn3HKqVYCWSGodzOTaiIUiQBNuOdUiQGc1s7a01N6IhSI/qdTl4CJAMfcvtTdiIWIhYqF5E0Es7Fhi2sWkN1ZGPSsDsfB2Q0hVA2ODIxaOXTrtsLUBYt1/3t5FjUr7TfTkGFlZzRFW7whDLBw3x0tNea9PL+YlJvVGLHwPz/ZTzP1L7V1zbvYRVtNh61vdPtNCLPw++ZdFeyYCVNJhG8QKu3+pUyu5KSpAORZaj0BscMTC8Uo6PMawnuytAnSYG08hYBMgQDYXSp0ECJATFM1sAgTI5kKpkwABcoKimU2AANlcKHUSIEBOUDSzCWwBQiy02VDqIIBY6IBEE02AI0yzocZBgAA5INFEE1AByrHQ+jxigyMWiu2MYT3ZWwUox0LrK48NjlgoAhTDerK3CJD0yNrcV1YjFor8rNkUEaBUNTA2OGIhYiFioXkTQSzsWFKlDayMelYGYuHthhBz/1J7IxaOXYq5fUX9un1xsZVNehdduP0mOvUQig3OEVbvCEMsHDdHxMLrp3cW//xJaomIhf/EhFh4VZ8RfNs4vf5YbR9hiyS1tpcTcbDy3CZTj1XXXLgIUEmHbdwrUu29woOX3BQVoJOaWd/j1N6IheOVdHhMpa4GVwE6zI2nELAJECCbC6VOAgTICYpmNgECZHOh1EmAADlB0cwmQIBsLpQ6CRAgJyia2QS2ABl/F2Y3phQCRwLi78KOzXgOAZsAR5jNhVInAQLkBEUzm4AK0EnNrF8ktTdiob2Va7ioAJ3UzPrSUnsjFooApVJXg4sArZHUGpjJtRELRX4m3JKqRYBi7l9qb8RC9UuDu2736cdfGmyBS9oUEaCJ47ayun8fyuf7v737C1bXXLgdoJg4kdobK6OelYFYeHtXkaoGxgZHLBy7NBHoJtVF/bp9cZOpx6qLLnw7wn78RzZTD6HY4BxhtY4w8Y9sSo+s3UZWVlf+UsFULjUXbr+JrvnleP0MTP1awMKD19wUEaD2qVPJbwWs6de1cKdiq7lwEaCSDlu/A1WeW6qWWHLhKkDKI+ubuLIasXC8kg6PSzZFBegwN55CwCZAgGwulDoJECAnKJrZBAiQzYVSJwEC5ARFM5sAAbK5UOokQICcoGhmE9gChFhos6HUQQCx0AGJJpoAR5hmQ42DAAFyQKKJJqAClKoGxgZfI9B1hLGpx3qXXLgKkPLIOseV1YiF4n6wZFNEgJIstLby2OCIhSI/Maxne4sAJVlobemxwRELEQv5xkLzJjLRNRELOzWsjI9mfiZciuoo9hEWU+BSexf16/ZIPOLCVYA+bESeri/mayVmyMV6F/XreoASsRVd+BYgxMLby2RyjKysLnqEIRbewrP/kKoGxgZHLBxbFZP3avp1bW2xlU1611y4/R6opsPWA9g+8iopPaZOreamiACVdNjGHSzV3is8eMlNUQFaIqn1hEyujVg4XkmHxwm3nGoVoMPceAoBmwABsrlQ6iRAgJygaGYTIEA2F0qdBAiQExTNbAIEyOZCqZMAAXKCoplNYAsQYqHNhlIHAcRCBySaaAIcYZoNNQ4CBMgBiSaagApQTIFL7V3Sr+uEH2/hKkBLJLW+C5NrIxaK+8GEW061CNBZzawtLbU3YqHITyp1ObgIUMz9S+2NWIhYiFho3kQQCzuWVHehqJywL/0RF24fYamGXGxwxEIp6+0Sn3b5cqpVgBINOcTCU7omYuE492M3eo4wvrHwV/Tr2qsjtrJJb8TCcQuaCHST6pp+XVvbZOqx6poLt98D1XTYegDbZ1qIheP1OB7XcBEBKumwDVKF3b/UqZXcFBWgHAutRyA2OGLheCUdHmNYT/ZWATrMjacQsAkQIJsLpU4CBMgJimY2AQJkc6HUSYAAOUHRzCZAgGwulDoJECAnKJrZBLYAIRbabCh1EEAsdECiiSbAEabZUOMgQIAckGiiCagApRpyscERC8V2xrCe7K0ClGOh9ZXHBkcsFAGKYT3ZWwRIemRt7iurEQtFftZsighQqhoYGxyxELEQsdC8iSAWdiwx7WLSGyujlpVhH2Ex9y+1N2IhYuE1Yi0W9ev2m29MmZz0Lrrw7Q7ENxbe3pJMzs+V1UXPbr6x8Bae/YeJ3LeyGrFwbNWv6Ne1tcVWNumNWDgCFDPgKkuPsZVNetdcuP2nsEtJh20EMNXeKzx4yU1RATqpmfU9Tu2NWDheSYfHVOpqcBWgw9x4CgGbAAGyuVDqJECAnKBoZhMgQDYXSp0ECJATFM1sAgTI5kKpkwABcoKimU1gCxBioc2GUgcBxEIHJJpoAhxhmg01DgIEyAGJJpqACtBJzaxfKLU3YqHYzlTqanAVoJOaWV9aam/EQhGgVOpqcBGgNZJaAzO5NmKhyM+EW1K1CFDM/UvtjViIWIhYaN5EEAs7llR3oaicsC/9ERduH2GpamBscMRCxELEQvME+3nFwt//+M9hRan34tjgHGG13Pj//nG5PF2fDgGq6bD1Sa6U+1Zeu+am2O+BajpsPUAT/+6Xra65KSJA7VOnkt8KWNOva+FOxVZz4SJAJR22cc4Wdv9Sp1ZyU1SAlEfWN3FlNWLheCUdHpdsigrQYW48hYBNgADZXCh1EiBATlA0swkQIJsLpU4CBMgJimY2AQJkc6HUSYAAOUHRzCZAgGwulDoJECAnKJrZBAiQzYVSJwEVIKUB9WFXVuOFib1dsikqQEoD6nNfWY0XJgK0ZFNEgJIkorby2OB4YSI/Maxne4sApZpdscHxwvDC8MLMmwheWMcS+635SW9+qb7WL9XbR1hM3UrtjReGF4YXZp5gP68X9qPWMzlGVlZzhP0MR1hNBam/ble6WSuvXXNT7PdANRWkHqBfVvyaLKzmpogAtY+V8MKOb1ZWcsELu+1GzJ4qqUf1tcVWNuldcuHiDlTZvao8t1Q1q+TCVYBudwt+gMA9AgToHh3qpgQI0BQRDe4RIED36FA3JUCApohocI8AAbpHh7opgRag6/bf/ls60w40gEAn8GWPzdsXzj1/3P/7DBkI/BsCry03l/8B2ros3UR10BIAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{array}{cccccccccccccccccc}0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\end{array}\\right]$$"
      ],
      "text/plain": [
       "⎡0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎤\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎢0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎥\n",
       "⎢                                                    ⎥\n",
       "⎣0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0⎦"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stiff_exact - stiff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Example: $2 \\times 1$ cantilever beam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = Matrix([\n",
    "        [0, 0],\n",
    "        [1, 0],\n",
    "        [2, 0],\n",
    "        [0, 1],\n",
    "        [1, 1],\n",
    "        [2, 1]])\n",
    "\n",
    "elems = [[0, 1, 4, 3],\n",
    "        [1, 2, 5, 4]]\n",
    "\n",
    "mat_props = [[1, S(3)/10, 1]]*4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "stiff, mass = assembler(coords, elems, mat_props, const=\"plane_stress\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAC0AAAErCAMAAAB96voBAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRIlmzd0i77t8bBwggJIAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAKzSURBVGgF7VvreoMgDMVrO6914/2fdVxGuEgEv7pV1/hDwzyNmIQDJJMVXB0l2zomDWKs4FUtjmYLzAYJKblEF5tAe3NA0G3VD/1KCYYeRb/aabZqlYSgh4e824956EUZqOGtD0d0c4XueGCnOLrllVTa8TpH98x7CSv0xf4irnvWujPR+3rC9Fs2eW/JlkV2dsi0oPZOlekdNknPP8JAidtEIHsRnyFYdO2IiLUecSXS7VpDy+e3CbFPOEhexT43fgsDCGGfz/goJvZ5ei4OHXDGUdzWwQQd77T46yxm2+fXJ3H1xzAbsQ+xjxtfWFTR2ieME9dqroxZ0MVYeT+a2MdaT0qIBWnntZqLt3ZetPbJiiofBC0kBuG+L1wZTeyT5UtiH2IfL1CuPOZft4MpnNSctGBklgIrF1U5QYMhO3QLYI2DTnuH0I7pfsRGJbd1I2XBuZp4CbnnFNp/1n40zFLtOMExqhW8LjKIs37GwTv0vloqSA6k+t2LIkIJEZ5Cy5JDwU3dIYnuRPaBi5M6UmgJsmWHHPS0x5e9qj1k9mSwYGxtrzSpUyfqFJ3ZmaX6XYxN0yzGPSn0Q0WKeZREb7GPwelrmn1cfKonLhbd7/ggaGG6Ke9josaY6risM7CPUc2wnRfCPkjWGfElZZ1Xa00wfCAgFgxQpnkuNLGP8Yu+It55+50XsU9WnPggaCFRBfd94X3QxD5Znif2WdVGae3zuvyJG7NX5ipiH9eT6B6Q2IfYxwuUK4/587DmudmHMjm/l8lZxyC2mpEj7Xb/8EafaCCZnK979P/CKJPzLzM5YRVrHVU2bPwqVpq/qS5lbWeknTY5bRUL4mRdxTKv+nNNxElQ80qg9+kWX7dQzSswmWz+ac1LOTj6dRXUvODrqlZ+OFXXpjTodR1qXurrqrpm39+yZ4lNWRgFAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}0\\\\0\\\\0\\\\0\\\\0\\\\- \\frac{1}{2}\\\\0\\\\0\\\\0\\\\0\\\\0\\\\- \\frac{1}{2}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡ 0  ⎤\n",
       "⎢    ⎥\n",
       "⎢ 0  ⎥\n",
       "⎢    ⎥\n",
       "⎢ 0  ⎥\n",
       "⎢    ⎥\n",
       "⎢ 0  ⎥\n",
       "⎢    ⎥\n",
       "⎢ 0  ⎥\n",
       "⎢    ⎥\n",
       "⎢-1/2⎥\n",
       "⎢    ⎥\n",
       "⎢ 0  ⎥\n",
       "⎢    ⎥\n",
       "⎢ 0  ⎥\n",
       "⎢    ⎥\n",
       "⎢ 0  ⎥\n",
       "⎢    ⎥\n",
       "⎢ 0  ⎥\n",
       "⎢    ⎥\n",
       "⎢ 0  ⎥\n",
       "⎢    ⎥\n",
       "⎣-1/2⎦"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "load = zeros(12, 1)\n",
    "load[5] = -S(1)/2\n",
    "load[11] = -S(1)/2\n",
    "load"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "stiff2 = stiff.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "stiff2[0, :] = eye(12)[0, :]\n",
    "stiff2[:, 0] = eye(12)[:, 0]\n",
    "stiff2[1, :] = eye(12)[1, :]\n",
    "stiff2[:, 1] = eye(12)[:, 1]\n",
    "stiff2[6, :] = eye(12)[6, :]\n",
    "stiff2[:, 6] = eye(12)[:, 6]\n",
    "stiff2[7, :] = eye(12)[7, :]\n",
    "stiff2[:, 7] = eye(12)[:, 7]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxkAAAAyBAMAAADCTgukAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEImZRO/dMlQiu6vNZnZmcXX2AAAACXBIWXMAAA7EAAAOxAGVKw4bAAALJ0lEQVR4Ae1bb2xcRxGfd+ez73w++5TK/pKoPjkUpAaSS9wKVFfKSQiJD5XiSEmhIagmFRGESjk1El99AlUI8cGnRlQKAvmg/FP4EIeoVCi0uQKiiNLGIOUDEm0PoggqUB2XQvPHrZmd2X3/bt7uu0i174P3w9vZmd/OzL65fbtv3xwAlm1VdcXyIa62rht1B7wlbek7DWNy4OE6kxkTFSPZqj/oO/CQNlA6UtbUggnCrg/a9pb++B0otTXn6S8wUbilGfmDmuipumvql0786Qd+gpipv+s5KOKfmvqpjzklIgQmAVOjBQXg7dvVlPhhXmH6LHqe2f2RMFPThal7G0Y2ZH7WAi7CIiCjf2D03GaiuKoZ4zVN9FJ5n4dv43CeLFs64dNxTw121LLtZFChCgvLBnMiGYcSssUGCWhHW1UBjC8XHgGw/xC3Q+Yd8H4LDwq6TkP+PSMbbAkAn0U22BABGT3e0ID9GHAsY22qAJ7RdU9VcRaKVe/CJVs0imXIzaL64Vay6tEbMDanMaNHknFAtrRBAlrRFkUsehbgAmyfMb9IEX8e4DAMV+GsIF1pwP+M7FpLABgW2dCGCEgXyM9qwHyHiLElbuetHuk+XRXGMnMdYMUWjdwcFN8duNnVN8wYOASTbY353h/Cki6abNGFgA50V/coYw1gD8CodewPA1yqH2hGO+rWp2reGrCs8NeWCNFMskEXAhr0MS2erBIxucTtQW5qYdpqsgWj7zqiMbKK0RiedalcWNaYStpoENCBtlv1cNFcqTuicbmM0fhRkiJ8UrEsM9JKwih+EA0CGvRCkzvpMJhovKHZNo3dshzOjXcc0cBeg6tjL+2+r7t7iOMdA8aUmimjQUAXOmRBIo/i3Kg5ooH9TpbXnrhakxTARAVY9pe00SCgQZtJEIvGCdGWi4mPqZH33NE4UJl8BXINi7bCh/cCY+6GlNEgoAttMalE+Ji6vOyMxsDb3toyfFbS9dRjZZZ5nZTRIKCPLumHZDQaQ/gLv5PyKnzthjsaX4HJVcgetBq4p8aYStpoENCFtpoEGOnk9zSc0Rhc8tbLcH9ZUpY5xLIMpIwGAX209z4rjUYj63yuS64AlN56wrVuAGSrasdUsgc8d5Aw+XrKaBDQiZa9Driv/2OPa90AmAbcOcEbtaBbiDpcJtk30kaDgAH6j6wqGo3cUshAT6R6YbHuqQCu4srRhhI+0pJKoQnDNwizDVJGg4BOdJLBgI+7c/ueCkarAI9jNJaDToY6A3CuoWR/qqSMhqeAdGEdK02qo9GYF0wx3HUdWXJFI1+FM/hiYpsbY6sYDcJ86+LF/75gs0mRxwsBnWibJpa9rPc7ycivQ6GJLybS3MAH2LmGkn364sVLn+sk6+CIY9jzCvjPAL3INz4aDR0iizpRVHoE5jG49rmBZxdvDuG6MSuqIGauBYM3Dcb29qdtsUEC2tHJJklypTmEL0L2uYFnGKPNeVw3BF1o/KSR4RgsJdjhQq6FOLqoGn/MWKLRSNxPEzbxMvqi9woKrdEoHJuamYOfw3gjUQ2UlmGhYzC3knHGFhskoB1tVQXw7/IO/G3ao/H01M5fQLHiSXuq1yD/vpGNtWzGQtEgoEEPz1GvaDQwyHdUdk7XAGYefRnHlFRG1tfX5yAz/VwSQPE/Of1DPH0jzL3rz1uQZIsNEtCOtigiUWYa34Kyz67ZLF5eX8d1eudMXVA2MP1ATcsy525bbgLZYEME9NHqKANLJBref4jn7ftzl02JR+A+vtz1se4DV4nXB0PQq2kkGkO87R1peF+Oeyjx4ph+a/8eJppxnyReHLMJbX3n9XHhYkW5kL9BjlwB+CgRoYvEC4n7kSzNQqYdc0zixSCb0tRPJb2M8L6thK9wWH4GsFCOOSXxYpB+aw4vQdeRtMTrC7/5O9MILx/HyaVRPvBGyYFGzEeJF4P0WxOn/QBP9sAziRdIN5HiPSEeBWPJcBi4UrPmQCfqmcSLIvqvhR8LBvT5j++cxPOFm0ngMbIql9RCd09NkZDF4w2AwtsA8xVFBUXiBdL+pHAtHMKhRIrEiwA2q4Efs1QpHgLYcV7T9Ngq4NzoiobA4z79e12sCtEQeH0xgi9qL57swMReposUDempJPFgx3FVXgLAFztdeBUKje8OMawupIdIQZkxjHUcLT2VJB6MBUriOrzfqSE+XrdhwjLSFL8JkpLAJKNNNL5aAbi/Tm7wkwpwxV5oxPySeDFIvzVxxdZb9sAziRdIN5HibRRk8EkFWd7a6sX8LJ5GlmOeSbwYpN+amJ1ivqn5rkk8X7iZhF7F8TMLFr3D5aDgm96Dcc8kXhzTZ+3SLGTbMZ8kXgyyOU0+IdR7wJWackKflgw2vM9gKo3i+EXi+cIwEesXFvl0GgyBUwMJjWk0kfIlGK+n4YU7xfFhmaGdGCfAaApqPiEc5kUEd35Y8vyV1Nv9rzoM6pxQ3UHiaVG0ivWLCnUrDYagqYGs+Hj0+brtredw1rt5ER9j+IjMNJwYJ8BoMrXHW/GxNjFwr4FFn10RB77JVeQq8XzAq+VTDWwkYzhdVeGTMehEFUb+NvQJBbMDUXwKQFsldL5JVeSSlmeMSnitkLOEsSFiwr6IgMCtyBCJre985AwXeL5wt0rQ26ckHglVRuyj6xTgRAxQuirBkzH47bwFufV1jbABUdMJ3yo28EiBrtFLKp7KjNVGJTxr5CxhRYuYsC8iIPArMkRi681fNBr4pd0Ur2OooJZ4SsoZsb9+Qf00kzAoWlHpqqpYMADXWjCy603C2YH4tQ7XPrbK8Lu5ilzT8CgzVhuV8KyQs4QVLWEivkiAkFfXQkNktt7NRqOxP+hSCkifkngsXMG5sURkMgYoXVWBLBhQqalFVmoHIkZl37JV7rBs+oXqVDz1gVQblfCsjrKEiZQwEV8kQOBTZIjMLl6nOhqNrteMQIWdCqJhxfEZpRWiUlP9aFiRAJVYNBzwZHEoGskglCwk3+YefBGGmGuR3Wg05hvE7P2iovH86x93dZyouBCgUlOLr12tO4E4c1SubgqrblUUDbdRzBJOKr34IgxxsUOKo9G44+w2FY3vw2ItyVnmY7qqHaCWFIzGcLnAe207mrJv3VbtSkiqouE2qrKEk0oPvkhDFLPbsteTrDn4Kho4oDkHjI5hrBiTmvpjK4qFOvvWadWtSkUDi9Oo/vIgKOzBF2mI2nJ0bgyk+UkKvuhUquxNSRbmHXZNDpOaur8Z7ibSeZ2r67Yqdg8zdTScRnMHw71CdC++CEOUs6LhsZCFXkicG8U2/YnD0uuMSle1yFHkUSLrr+j/FHYk/sNd5eq6rbrUKLmKhtMoZQknaOvBF08Yov41RPOphLPbBOtxNkYj19bfDuMyv03pqn5LJCg/tXMUYL9rEgGn4LqtinZiTHU3nEYpSzjW0TQpEzidL9IQzXodfVJBrmr091ZjNLKUgm7rhu9qJ90PoFwL0InQe6hF5ZEUVi3dfZGKhtMoOjZoeRT34Atqilpb0LdFRwO/wlDpSnlhtvOK0cCDgytaaQKc0lUTZAF7rIVHKJhnnabcSmE1jR4VDafREmUJJ6rrwZeuIf5Ga+XjQvWY4WL4upmyoozY07v22uGcrmrHgEpNLey7UHbASKyyb51WUyiizFi3UcoSTlLXgy9dQ1THZFTmK1QVzRRcrDN/67qBd2C8po2t8N0v4HdvKvk5TWxVG3cHntGm1F+lqax0NDGj661qw+5AqaVNTZgnVum2flIXOxvmxZYhvgMP6RuRP+qvltuXNO++rZu0sXfAO6/tfXdZEf8HfBbj9EAuTMIAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left\\{\\left ( 0, \\quad 0, \\quad - \\frac{91}{15}, \\quad - \\frac{26}{3}, \\quad - \\frac{364}{45}, \\quad - \\frac{1144}{45}, \\quad 0, \\quad 0, \\quad \\frac{91}{15}, \\quad - \\frac{26}{3}, \\quad \\frac{364}{45}, \\quad - \\frac{1144}{45}\\right )\\right\\}$$"
      ],
      "text/plain": [
       "⎧⎛      -91          -364   -1144         91         364  -1144 ⎞⎫\n",
       "⎨⎜0, 0, ────, -26/3, ─────, ──────, 0, 0, ──, -26/3, ───, ──────⎟⎬\n",
       "⎩⎝       15            45     45          15          45    45  ⎠⎭"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol = linsolve((stiff2, load))\n",
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "/* Based on Lorena Barba template available at: https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css*/\n",
       "@font-face {\n",
       "font-family: \"Computer Modern\";\n",
       "src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "div.cell{\n",
       "width:800px;\n",
       "margin-left:16% !important;\n",
       "margin-right:auto;\n",
       "}\n",
       "h1 {\n",
       "font-family: 'Alegreya Sans', sans-serif;\n",
       "}\n",
       "h2 {\n",
       "font-family: 'Fenix', serif;\n",
       "}\n",
       "h3{\n",
       "font-family: 'Fenix', serif;\n",
       "margin-top:12px;\n",
       "margin-bottom: 3px;\n",
       "}\n",
       "h4{\n",
       "font-family: 'Fenix', serif;\n",
       "}\n",
       "h5 {\n",
       "font-family: 'Alegreya Sans', sans-serif;\n",
       "}\t\n",
       "div.text_cell_render{\n",
       "font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "line-height: 135%;\n",
       "font-size: 120%;\n",
       "width:600px;\n",
       "margin-left:auto;\n",
       "margin-right:auto;\n",
       "}\n",
       ".CodeMirror{\n",
       "font-family: \"Source Code Pro\";\n",
       "font-size: 90%;\n",
       "}\n",
       "/* .prompt{\n",
       "display: None;\n",
       "}*/\n",
       ".text_cell_render h1 {\n",
       "font-weight: 200;\n",
       "font-size: 50pt;\n",
       "line-height: 100%;\n",
       "color:#CD2305;\n",
       "margin-bottom: 0.5em;\n",
       "margin-top: 0.5em;\n",
       "display: block;\n",
       "}\t\n",
       ".text_cell_render h5 {\n",
       "font-weight: 300;\n",
       "font-size: 16pt;\n",
       "color: #CD2305;\n",
       "font-style: italic;\n",
       "margin-bottom: .5em;\n",
       "margin-top: 0.5em;\n",
       "display: block;\n",
       "}\n",
       ".warning{\n",
       "color: rgb( 240, 20, 20 )\n",
       "}\n",
       "</style>\n",
       "<script>\n",
       "MathJax.Hub.Config({\n",
       "TeX: {\n",
       "extensions: [\"AMSmath.js\"]\n",
       "},\n",
       "tex2jax: {\n",
       "inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "},\n",
       "displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "\"HTML-CSS\": {\n",
       "styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "}\n",
       "});\n",
       "</script>\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open('./styles/custom_barba.css', 'r').read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "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"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}