3D7 - Basic FEM - heat equation.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "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.5.4"
    },
    "colab": {
      "name": "3D7 - Basic FEM - heat equation.ipynb",
      "provenance": [],
      "include_colab_link": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/chrisrichardson/FEMexamples/blob/main/3D7%20-%20Basic%20FEM%20-%20heat%20equation.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q-9rHoChR-ej"
      },
      "source": [
        "# Basic FEM Example: Heat Equation\n",
        "\n",
        "First import the required libraries for basic algebra and plotting"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6k4fzuB4R-el"
      },
      "source": [
        "from numpy import *\n",
        "set_printoptions(2, suppress=True)\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YhPPo81AR-eu"
      },
      "source": [
        "### Mesh\n",
        "We define the \"mesh\" in terms of its geometry and topology. The geometry is just a set of $(x, y)$ points - a $(p \\times 2)$ floating point array, where $p$ is the number of points. The topology is the list of points which make up each triangle, hence it is a $(n\\times 3)$ integer array, where $n$ is the number of triangles. Note that numbering starts from zero. Also, the order of the points (i.e. triangle orientation) is important."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DShH-MaxR-ew"
      },
      "source": [
        "geometry = array([[0.0, 0.0],\n",
        "                  [2.0, 0.0],\n",
        "                  [0.0, 1.0],\n",
        "                  [2.0, 1.0]])\n",
        "\n",
        "topology = array([[0, 1, 2],\n",
        "                  [1, 3, 2]])\n",
        "\n",
        "mesh = (geometry, topology)"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Bz2hm-wYR-e5"
      },
      "source": [
        "Here is a simple plotting routine to display the mesh"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kVmeEPjER-e7"
      },
      "source": [
        "def plot(mesh):\n",
        "    geom, topo = mesh\n",
        "    x = geom[:,0]\n",
        "    y = geom[:,1]\n",
        "    plt.triplot(x, y, topo)\n",
        "    xmax = x.max()\n",
        "    xmin = x.min()\n",
        "    ymax = y.max()\n",
        "    ymin = y.min()\n",
        "    dx = 0.1*(xmax - xmin)\n",
        "    dy = 0.1*(ymax - ymin)\n",
        "    plt.xlim(xmin-dx, xmax+dx)\n",
        "    plt.ylim(ymin-dy, ymax+dy)\n",
        "    return\n"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lxFZG4pCR-fB",
        "outputId": "bee25de5-2686-4e06-fb1b-2bb73625f40c",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        }
      },
      "source": [
        "plot(mesh)"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3RVZdr+8e+dSu+hSAsdQ5ESOiTOSFdBERXsooBIj++MOoWfZea1jG8AEaSIOlgABcWIKMQyCV1ChyAQmjQREOkCkef3R+JMBoMcIMlOzrk+a7HWOXs/K+daz5y51nafkzvmnENERAq+IK8DiIhIzlChi4j4CRW6iIifUKGLiPgJFbqIiJ8I8eqFy5Ur5yIjI716eRGRAmnlypWHnHMR2Z3zrNAjIyNJSUnx6uVFRAokM9t1sXO65SIi4idU6CIifkKFLiLiJ1ToIiJ+QoUuIuInVOgiIn5ChS4i4idU6CIifkKFLiLiJ1ToIiJ+QoUuIuInVOgiIn5ChS4i4icuWehm9rqZfW9mGy5y3szsZTNLM7N1ZtYs52OKiMil+HKF/ibQ9TfOdwPqZP4bALx69bFERORyXXIeunMu2cwif2NJT2Cac84By8yslJlVcs7tz6GMV23kzDUUCQuma8OKXkcRET+ydvePlC0WTt+W1byOAuTMH7ioDOzO8nxP5rFfFbqZDSDjKp5q1fJmA86fd3y4ei8A7yz/Nk9eU0QCiz8Vus+cc5OByQDR0dEuL14zKMh+dez5Xo2oXb5YXry8iPix3hOXUqFEuNcx/i0nCn0vUDXL8yqZx/KNZtVKUSQshF7NKvPM3FRGfbSR4R3rMCCmJqHB+qKPiFyZZtVKUTTcs7/k+Ss50WYJwH2Z33ZpDRzNT/fPf2EGvZpVIXFkLJ0aVOAf8zfT85XFbNh71OtoIiI5wpevLU4HlgL1zGyPmT1kZo+Y2SOZS+YB24E0YArwaK6lzQERxcMZf1czJt3bnIMnztBz/GJe+Owbfjr3s9fRRESuii/fcul7ifMOGJxjifJIlwYVaV2jLP87bxOv/msb8zd8x/O3NaZljTJeRxMRuSIBfQO5ZJFQXujdmLcfasXZn89zx6Sl/HXOBk6cSfc6mojIZQvoQv9F+zrlWDAyhn7tavD28l10jk/iq83fex1LROSyqNAzFQkLYdTNUcx6pC1FwkN48I0VxM1cw5GTZ72OJiLiExX6BZpXL80nw9oz7Pe1SVi7j06jk/hk3X4yPioQEcm/VOjZCA8JJq5zPT4e2p5KJQsz+N1VDHxrJQeO/eR1NBGRi1Kh/4ZrK5Xgw0fb8mS3+iRtOUjH+CRmrvhWV+siki+p0C8hJDiIgbG1+GxEDNdWKsHjs9dzz9TlfHv4lNfRRET+iwrdRzXKFWVG/9b87ZaGrN19lC5jkpm6aAc/n9fVuojkDyr0yxAUZNzTujoLRsbQumYZnp2bSu+JS9h64LjX0UREVOhX4ppShXn9gRaM7dOEnYdOcuPLi3j5i62cTT/vdTQRCWAq9CtkZvRsUpnP42Lp0rAi8Ylb6PHKItbt+dHraCISoFToV6lssXDG9W3KlPuiOXLqLLeMX8xz8zZx+qyGfYlI3lKh55BOURVIjIvlzhZVmZS8nW5jk1m2/bDXsUQkgKjQc1CJQqE816sx7z7civMO+kxexp8/XM/xn855HU1EAoAKPRe0rV2O+SNieLh9DaZ//S2dRyfz5TcHvI4lIn5OhZ5LCocF85ebopg9qC3FC4XQ780URsxYzQ8a9iUiuUSFnsuaVivN3KEdGNGxDp+s30/H+CQS1u7T+AARyXEq9DwQFhLEiI51mTu0A1XLFGHY9NX0n5bCd0c17EtEco4KPQ/Vq1icDwa15S83XsuitEN0ik9i+tca9iUiOUOFnseCg4yHO9Rk/ogYGlYuyZMfrOeuKcvZdfik19FEpIBToXuketmivNu/Fc/1asSGvRnDvl5buF3DvkTkiqnQPWRm9G1ZjcS4WNrXLsffPtlEr1eXsPk7DfsSkcunQs8HKpYsxJT7ohnXtyl7fjjFTeMWMjpxi4Z9ichlUaHnE2bGzdddQ2JcLDc2qsTYL7Zy07iFrNmtYV8i4hsVej5TpmgYY/o05fUHojn+Uzq9Jizmb3NTNexLRC5JhZ5P/b5+BRaMjKFvy2q8tmgHXcYks2TbIa9jiUg+pkLPx4oXCuXvtzZixoDWBBncNWU5T36wjmMa9iUi2VChFwCta5bl0+ExDIypycwVu+kUn8TnqRr2JSL/zadCN7OuZrbZzNLM7Ilszlczs6/MbLWZrTOz7jkfNbAVDgvmye7XMmdwO0oXCePhaSkMnb6aQyfOeB1NRPKJSxa6mQUD44FuQBTQ18yiLlj2F+A951xToA8wIaeDSobGVUqRMKQ9j3Wqy/wN39EpPok5q/dqfICI+HSF3hJIc85td86dBWYAPS9Y44ASmY9LAvtyLqJcKCwkiKE31OGTYe2JLFeUETPX8NA/U9j342mvo4mIh3wp9MrA7izP92Qey+op4B4z2wPMA4Zm94PMbICZpZhZysGDB68grmRVp0JxZj3SllE3RbF022E6j07m7WW7OK/xASIBKac+FO0LvOmcqwJ0B94ys1/9bOfcZOdctHMuOiIiIodeOrAFBxn92tdg/ogYrqtakr/M2UDfKcvYcUjDvkQCjS+FvheomuV5lcxjWT0EvAfgnFsKFALK5URA8U21skV4+6FWvHhbY1L3H6PrmGQmJW0j/WeNDxAJFL4U+gqgjpnVMLMwMj70TLhgzbfADQBmdi0Zha57KnnMzLijRVU+j4sltm4Ez336DbdOWELqvmNeRxORPHDJQnfOpQNDgPnAJjK+zbLRzJ4xsx6Zyx4D+pvZWmA68IDT1y48U6FEISbd25zxdzVj/9HT9HhlEf+3YDNn0jU+QMSfhfiyyDk3j4wPO7MeG5XlcSrQLmejydUwM25sXIm2tcry7CepjPsyjU83fMcLtzWmefXSXscTkVyg3xT1c6WLhhF/RxPeeLAFp86k03viEp7+eCOnzqZ7HU1EcpgKPUD8rl55FsTFcm/r6ryxeCedRyezaKuGfYn4ExV6ACkWHsIzPRvy3sA2hAUHcc/U5fxx1lqOntawLxF/oEIPQC1rlGHe8A4Mur4Ws1ftpVN8EvM3fud1LBG5Sir0AFUoNJjHu9bno8HtKFcsnIFvrWTwO6s4eFzDvkQKKhV6gGtYuSQfDWnHH7rUIzH1AB3jk5i9co+GfYkUQCp0ITQ4iMG/q8284R2oXb4Yj72/lgfeWMFeDfsSKVBU6PJvtcsX4/2BbXjq5ihW7PyBzvFJTFu6U8O+RAoIFbr8l6Ag44F2GcO+mlUvzaiPNnLn5KVsO3jC62gicgkqdMlW1TJFmNavJS/dfh1bDpyg29iFTPhXGuc07Esk31Khy0WZGb2bVyExLoYb6pfnxc82c8v4xWzYe9TraCKSDRW6XFL54oV49Z7mvHp3Mw4cO0PP8Yv5x/xv+Omchn2J5CcqdPFZt0aV+DwuhlubVmb8V9vo/vJCUnb+4HUsEcmkQpfLUqpIGC/dfh3T+rXkzLnz3D5pKU8lbOTkGQ37EvGaCl2uSEzdCBaMjOH+NpH8c2nGsK/kLfqbJiJeUqHLFSsaHsJTPRrw/sA2FAoN4r7Xv+Z/3l/Lj6fOeh1NJCCp0OWqRUeW4ZNhHRjyu9p8uHovHeOT+XT9fq9jiQQcFbrkiEKhwfxPl3okDGlHhRLhDHpnFY+8tZLvj/3kdTSRgKFClxzV4JqSfDS4HY93rc+Xm7+nY3wS76fs1rAvkTygQpccFxIcxKDra/Hp8A7Uq1icP8xax32vf83uH055HU3Er6nQJdfUiijGzAFteLZnA1btOkKXMcm8uXiHhn2J5BIVuuSqoCDj3jaRLIiLpUVkGZ76OJXbJy0l7fvjXkcT8TsqdMkTlUsV5s0HWxB/x3VsO3iC7mMX8cqXWzXsSyQHqdAlz5gZvZpVIXFkLJ0aVOClBVvo8YqGfYnkFBW65LmI4uGMv6sZk+5tzqETGcO+nv9Uw75ErpYKXTzTpUFFPh8ZS+9mVZiYtI3uYxfy9Q4N+xK5Uip08VTJIqG80Lsxbz/UinPnz3PHpKX8dc4GTmjYl8hl86nQzayrmW02szQze+Iia+4ws1Qz22hm7+ZsTPF37euUY/6IGPq1q8Hby3fROT6JrzZ/73UskQLlkoVuZsHAeKAbEAX0NbOoC9bUAZ4E2jnnGgAjciGr+LkiYSGMujmK2YPaUjQ8hAffWEHczDUcOalhXyK+8OUKvSWQ5pzb7pw7C8wAel6wpj8w3jl3BMA5p0sruWLNqpVm7rD2DPt9bRLW7qNjfBJz1+3T+ACRS/Cl0CsDu7M835N5LKu6QF0zW2xmy8ysa3Y/yMwGmFmKmaUcPKjZ2XJx4SHBxHWux8dD23NNqcIMeXc1A99ayQEN+xK5qJz6UDQEqANcD/QFpphZqQsXOecmO+einXPREREROfTS4s+urVSCDx9ty5Pd6pO05SAd45OYueJbXa2LZMOXQt8LVM3yvErmsaz2AAnOuXPOuR3AFjIKXuSqhQQHMTC2Fp+NiCGqUgken72ee6Yu59vDGvYlkpUvhb4CqGNmNcwsDOgDJFywZg4ZV+eYWTkybsFsz8GcItQoV5Tp/Vvz91sbsnb3UbqMSWbqoh38rGFfIoAPhe6cSweGAPOBTcB7zrmNZvaMmfXIXDYfOGxmqcBXwB+cc4dzK7QErqAg4+5W1UmMi6FNrbI8OzeV3hOXsPWAhn2J+HQP3Tk3zzlX1zlXyzn398xjo5xzCZmPnXMuzjkX5Zxr5JybkZuhRSqVLMzU+6MZ26cJOw+dpPvLC3n5i62cTdewLwlc+k1RKbDMjJ5NKvN5XCxdG1YiPnELPV5ZxNrdP3odTcQTKnQp8MoWC2dc36ZMuS+aI6fOcuuExTw3bxOnz2rYlwQWFbr4jU5RFUiMi+XOFtWYlLydbmOTWbZdH+VI4FChi18pUSiU53o14t3+rXBAn8nL+POH6zn+0zmvo4nkOhW6+KW2tcrx2fAY+neowfSvv6Xz6GS+/OaA17FEcpUKXfxW4bBg/nxjFB882o4ShULp92YKw2es5vCJM15HE8kVKnTxe02qluLjoe0Z0bEO89bvp9PoZBLWatiX+B8VugSEsJAgRnSsy9yhHahapgjDpq+m/7QUvjuqYV/iP1ToElDqVSzOB4Pa8pcbr2VR2iE6xScx/WsN+xL/oEKXgBMcZDzcoSbzR8TQsHJJnvxgPXdNWc6uwye9jiZyVVToErCqly3Ku/1b8XyvRmzYmzHsa0rydg37kgJLhS4Bzczo07IaiXGxtK9djr/P20SvCYvZ/J2GfUnBo0IXASqWLMSU+6IZ17cpe46c5qZxCxmduEXDvqRAUaGLZDIzbr7uGhLjYrmxUSXGfrGVm8YtZI2GfUkBoUIXuUCZomGM6dOU1x+I5vhP6fSasJi/zU3VsC/J91ToIhfx+/oVWDAyhrtaVeO1RTvoMiaZJdsOeR1L5KJU6CK/oXihUP52SyNmDGhNkMFdU5bz5AfrOHpaw74k/1Ghi/igdc2yfDYihoGxNZm5YjedRyeRmKphX5K/qNBFfFQoNJgnu13LnMHtKF0kjP7TUhjy7ioOadiX5BMqdJHL1LhKKRKGtOexTnVZsPEAneKTmLN6r8YHiOdU6CJXICwkiKE31OGTYe2JLFeUETPX8NA/U9j342mvo0kAU6GLXIU6FYoz65G2jLopiqXbDtN5dDJvL9vFeY0PEA+o0EWuUnCQ0a99DRaMjKFJ1VL8Zc4G+kxZxo5DGvYleUuFLpJDqpYpwlsPteTF2xqzaf8xuo5JZmLSNtJ/1vgAyRsqdJEcZGbc0aIqn8fFEls3guc//YZbJywhdd8xr6NJAFChi+SCCiUKMene5oy/qxn7j56mxyuL+L8FmzmTrvEBkntU6CK5xMy4sXElEkfG0qPJNYz7Mo0bX17Eyl1HvI4mfsqnQjezrma22czSzOyJ31h3m5k5M4vOuYgiBVvpomHE39GENx9swemzP9N74hKe/ngjp86mex1N/MwlC93MgoHxQDcgCuhrZlHZrCsODAeW53RIEX9wfb3yzB8Zw72tq/PG4p10Hp3Moq0a9iU5x5cr9JZAmnNuu3PuLDAD6JnNumeBFwD9GXWRiygWHsIzPRvy3sA2hAUHcc/U5fxx1lqOntKwL7l6vhR6ZWB3lud7Mo/9m5k1A6o65z75rR9kZgPMLMXMUg4ePHjZYUX8RcsaZZg3vAODrq/F7FV76Tg6ic82fOd1LCngrvpDUTMLAuKBxy611jk32TkX7ZyLjoiIuNqXFinQCoUG83jX+nw0uB0RxcJ55O2VDH5nFQePa9iXXBlfCn0vUDXL8yqZx35RHGgI/MvMdgKtgQR9MCrim4aVS/LRkHb8oUs9EjcdoGN8ErNX7tGwL7lsvhT6CqCOmdUwszCgD5Dwy0nn3FHnXDnnXKRzLhJYBvRwzqXkSmIRPxQaHMTg39Vm3rAO1C5fjMfeX8sDb6xgr4Z9yWW4ZKE759KBIcB8YBPwnnNuo5k9Y2Y9cjugSCCpXb4Y7w9sw9M9GrBi5w90jk9i2tKdGvYlPgnxZZFzbh4w74Jjoy6y9vqrjyUSuIKCjPvbRvL7+uX504frGfXRRj5eu4/nb2tMrYhiXseTfEy/KSqST1UtU4Rp/Vry0u3XseXACbqNXciEf6VxTsO+5CJU6CL5mJnRu3kVEuNiuKF+eV78bDO3jF/Mhr1HvY4m+ZAKXaQAKF+8EK/e05xX727GgWNn6Dl+Mf+Y/w0/ndOwL/kPFbpIAdKtUSW+iIulV9PKjP9qG91fXkjKzh+8jiX5hApdpIApWSSUf9x+HdP6teTMufPcPmkpTyVs5OQZDfsKdCp0kQIqpm4EC0bGcH+bSP65NGPYV9IWjdQIZCp0kQKsaHgIT/VowPsD21AoNIj7X/+ax95by4+nznodTTygQhfxA9GRZfhkWAeG/K42c9bspWN8Mp+u3+91LMljKnQRP1EoNJj/6VKPhCHtqFgynEHvrOKRt1by/TFNtA4UKnQRP9PgmpLMebQdj3etz5ebv6djfBLvp+zWsK8AoEIX8UMhwUEMur4Wnw3vQP2KJfjDrHXc9/rX7P7hlNfRJBep0EX8WM2IYswY0JpnezZg1a4jdBmTzBuLd/Czhn35JRW6iJ8LCjLubRPJgrhYWkSW4emPU7lj0lLSvj/udTTJYSp0kQBRuVRh3nywBfF3XMe2gyfoPnYRr3y5VcO+/IgKXSSAmBm9mlUhcWQsnRpU4KUFW+jxioZ9+QsVukgAiigezvi7mjHp3uYcPpEx7Ov5TzXsq6BToYsEsC4NKpIYF0vvZlWYmLSN7mMX8vUODfsqqFToIgGuZOFQXujdmHcebsW58+e5Y9JS/jpnA8d/Oud1NLlMKnQRAaBd7XLMHxFDv3Y1eHv5LrqMTuarzd97HUsugwpdRP6tSFgIo26OYvagthQND+HBN1YQN3MNR05q2FdBoEIXkV9pVq00c4e1Z9gNdUhYu4+O8UnMXbdP4wPyORW6iGQrPCSYuE51+XhoeyqXLsyQd1cz8K2VHNCwr3xLhS4iv+naSiX4YFBb/tS9PklbDtIxPomZK77V1Xo+pEIXkUsKCQ5iQEwt5o+IIapSCR6fvZ67X1vOt4c17Cs/UaGLiM8iyxVlev/W/P3Whqzbc5QuY5KZukjDvvILFbqIXJagIOPuVtVJjIuhTa2yPDs3ldteXcKWAxr25TUVuohckUolCzP1/mjG9mnCrsMnufHlhbz8xVbOpmvYl1d8KnQz62pmm80szcyeyOZ8nJmlmtk6M/vCzKrnfFQRyW/MjJ5NKvN5XCzdGlYiPnELPV5ZxNrdP3odLSBdstDNLBgYD3QDooC+ZhZ1wbLVQLRzrjEwC3gxp4OKSP5Vtlg4L/dtymv3RfPjqXPcOmEx/ztvE6fPathXXvLlCr0lkOac2+6cOwvMAHpmXeCc+8o598vH3cuAKjkbU0QKgo5RFVgQF8OdLaoxOXk73cYms3TbYa9jBQxfCr0ysDvL8z2Zxy7mIeDT7E6Y2QAzSzGzlIMHD/qeUkQKjBKFQnmuVyPe7d8KB/Sdsow/fbieYxr2lety9ENRM7sHiAb+kd1559xk51y0cy46IiIiJ19aRPKZtrXK8dnwGPp3qMGMr7+lc3wyX35zwOtYfs2XQt8LVM3yvErmsf9iZh2BPwM9nHNnciaeiBRkhcOC+fONUXzwaDtKFg6l35spDJ+xmsMnVBG5wZdCXwHUMbMaZhYG9AESsi4ws6bAJDLKXPM2ReS/NKlaio+Htmdkx7rMW7+fTqOTSVirYV857ZKF7pxLB4YA84FNwHvOuY1m9oyZ9chc9g+gGPC+ma0xs4SL/DgRCVBhIUEM71iHuUM7ULVMEYZNX03/aSnsP3ra62h+I8SXRc65ecC8C46NyvK4Yw7nEhE/Va9icT4Y1JY3Fu/gpQWb6RyfzJPdr6VPi6oEBZnX8Qo0/aaoiOS54CDj4Q41mT8ihoaVS/KnD9dz12vL2HnopNfRCjQVuoh4pnrZorzbvxXP92rExr3H6Do2mSnJ2zXs6wqp0EXEU2ZGn5bVSIyLpX3tCP4+bxO9Jixm83ca9nW5VOgiki9ULFmIKfc1Z1zfpuw5cpqbxi1kdOIWDfu6DCp0Eck3zIybr7uGxLhYbmp8DWO/2MpN4xay+tsjXkcrEFToIpLvlCkaxug7m/D6A9Ec/ymdXq8u4dm5qZw6m+51tHxNhS4i+dbv61dgwcgY7m5VjamLdtB1zEKWpB3yOla+pUIXkXyteKFQ/nZLI2YMaE2QwV2vLeeJ2es4elrDvi6kQheRAqF1zbJ8NiKGgbE1eS9lN51HJ5GYqmFfWanQRaTAKBQazJPdrmXO4HaULhJG/2kpDHl3FYc07AtQoYtIAdS4Ssawr8c61WXBxgN0ik9izuq9AT/sS4UuIgVSaHAQQ2+owyfD2hNZrigjZq6h35sr2Pdj4A77UqGLSIFWp0JxZj3SllE3RbFs+w90Hp3MW8t2cT4Axweo0EWkwAsOMvq1r8GCkTE0qVqKv87ZQJ8py9gRYMO+VOgi4jeqlinCWw+15MXejflm/zG6jklmYtI20n8OjPEBKnQR8Stmxh3RVfk8Lpbr60Xw/KffcOuEJaTuO+Z1tFynQhcRv1S+RCEm3tOcCXc3Y//R0/R4ZRH/t2AzZ9J/9jparlGhi4jfMjO6N6pE4shYejS5hnFfpnHjy4tYucs/h32p0EXE75UuGkb8HU1488EWnD77M70nLuHpjzdy8ox/DftSoYtIwLi+Xnnmj4zh3tbVeWPxTrqMSWbh1oNex8oxKnQRCSjFwkN4pmdD3hvYhrDgIO6d+jV/nLWWo6cK/rAvFbqIBKSWNcowb3gHHr2+FrNX7aXj6CQ+2/Cd17GuigpdRAJWodBg/ti1Ph8NbkdEsXAeeXslg99ZxcHjBXPYlwpdRAJew8ol+WhIO/7QpR6Jmw7QMT6J2Sv3FLhhXyp0EREyhn0N/l1t5g3rQO3yxXjs/bXc/8YK9hw55XU0n6nQRUSyqF2+GO8PbMPTPRqQsvMHuoxOZtrSnQVi2JcKXUTkAkFBxv1tI1kwMobmkWUY9dFG7py8lG0HT3gd7Tep0EVELqJK6SL888EWvHT7dWw5cIJuYxcy4V9pnMunw758KnQz62pmm80szcyeyOZ8uJnNzDy/3MwiczqoiIgXzIzezauQGBdDx2vL8+Jnm7ll/GI27D3qdbRfCbnUAjMLBsYDnYA9wAozS3DOpWZZ9hBwxDlX28z6AC8Ad+ZG4Ctx4kw65352/HjqrNdRRKSACgsO4n9vbURMnQj++tEGbhq3CICoSiU8TvYflyx0oCWQ5pzbDmBmM4CeQNZC7wk8lfl4FvCKmZnLJ9/52XIg475Xk2cSPU4iIv4mdf8xnHOYmddRfCr0ysDuLM/3AK0utsY5l25mR4GywKGsi8xsADAAoFq1alcY+fJ1b1SRIyfP0blBhTx7TRHxf09/nErpIqH5oszBt0LPMc65ycBkgOjo6Dy7ep9wd/O8eikRCSAPtqvhdYT/4suHonuBqlmeV8k8lu0aMwsBSgKHcyKgiIj4xpdCXwHUMbMaZhYG9AESLliTANyf+bg38GV+uX8uIhIoLnnLJfOe+BBgPhAMvO6c22hmzwApzrkEYCrwlpmlAT+QUfoiIpKHfLqH7pybB8y74NioLI9/Am7P2WgiInI59JuiIiJ+QoUuIuInVOgiIn5ChS4i4idU6CIifkKFLiLiJ1ToIiJ+QoUuIuInVOgiIn5ChS4i4idU6CIifkKFLiLiJ8yrKbdmdhDYlYcvWY4L/oKSANqX7GhPsqd9yV5e70t151xEdic8K/S8ZmYpzrlor3PkN9qXX9OeZE/7kr38tC+65SIi4idU6CIifiKQCn2y1wHyKe3Lr2lPsqd9yV6+2ZeAuYcuIuLvAukKXUTEr6nQRUT8hN8Vupl1NbPNZpZmZk9kcz7czGZmnl9uZpF5nzJv+bAnD5jZQTNbk/nvYS9y5jUze93MvjezDRc5b2b2cua+rTOzZnmdMa/5sCfXm9nRLO+VUdmt8zdmVtXMvjKzVDPbaGbDs1nj/fvFOec3/4BgYBtQEwgD1gJRF6x5FJiY+bgPMNPr3PlgTx4AXvE6qwd7EwM0AzZc5Hx34FPAgNbAcq8z54M9uR6Y63VOD/alEtAs83FxYEs2/z/y/P3ib1foLYE059x259xZYAbQ84I1PYF/Zj6eBdxgZpaHGfOaL3sSkJxzycAPv7GkJzDNZVgGlDKzSnmTzhs+7ElAcs7td86tynx8HNgEVL5gmefvF38r9MrA7izP9/DrTf/3GudcOnAUKJsn6bzhy54A3Jb5n4mzzKxq3kTL93zdu0DTxszWmtmnZtbA6zB5LfM2bVNg+QWnPH+/+Fuhy5X5GIh0zjUGEvnPf8GIXGgVGbNErgPGAXM8zpOnzKwYMBsY4Zw75nWeC/lboe8Fsl5dVsk8lqoJhygAAAFTSURBVO0aMwsBSgKH8ySdNy65J865w865M5lPXwOa51G2/M6X91NAcc4dc86dyHw8Dwg1s3Iex8oTZhZKRpm/45z7IJslnr9f/K3QVwB1zKyGmYWR8aFnwgVrEoD7Mx/3Br50mZ9o+KlL7skF9/l6kHF/UDL26b7Mby+0Bo465/Z7HcpLZlbxl8+czKwlGR3izxdEQMY3WICpwCbnXPxFlnn+fgnJyxfLbc65dDMbAswn49sdrzvnNprZM0CKcy6BjP9R3jKzNDI+/OnjXeLc5+OeDDOzHkA6GXvygGeB85CZTSfjWxvlzGwP8P+AUADn3ERgHhnfXEgDTgEPepM07/iwJ72BQWaWDpwG+vj5BdEv2gH3AuvNbE3msT8B1SD/vF/0q/8iIn7C3265iIgELBW6iIifUKGLiPgJFbqIiJ9QoYuI+AkVuoiIn1Chi4j4if8PFiCIKyfkWSMAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PevXNxprR-fH"
      },
      "source": [
        "### Element matrix\n",
        "The \"element stiffness matrix\", $K$ depends on the conductivity tensor $D$ and the triangle geometry.\n",
        "\n",
        "$$ K_e = \\int_{\\Omega_e} B_e^T D B_e d\\Omega $$\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4l1Tuw6IR-fI"
      },
      "source": [
        "def Kmat(D, p, geometry):\n",
        "    ''' Calculate K from the D matrix and the points '''\n",
        "    # D contains the conductivity tensor\n",
        "    # p contains the indices of the three points\n",
        "    # making up the triangle\n",
        "    assert(len(p) == 3)\n",
        "    x0, y0 = geometry[p[0]]\n",
        "    x1, y1 = geometry[p[1]]\n",
        "    x2, y2 = geometry[p[2]]\n",
        "\n",
        "    # Element area Ae\n",
        "    Ae = 0.5*abs((x0 - x1)*(y2 - y1) - (y0 - y1)*(x2 - x1))\n",
        "\n",
        "    # 'B' Matrix - representing the 'gradient' operator\n",
        "    B = array([[y1 - y2, y2 - y0, y0 - y1],\n",
        "               [x2 - x1, x0 - x2, x1 - x0]])/(2*Ae)\n",
        "\n",
        "    K = Ae*matmul(B.transpose(), matmul(D, B))\n",
        "    return K"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jDC4biQcR-fN"
      },
      "source": [
        "### Matrix assembly (LHS)\n",
        "The process of combining the local element matrices into a larger global matrix is called \"assembly\". It is really just a case of relabelling the indices with their global values, and adding the local matrices together."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "W_NetHAJR-fO"
      },
      "source": [
        "def assemble_matrix(mesh, Dmat):\n",
        "    geom, topo = mesh\n",
        "    np = len(geom)\n",
        "    Kglobal = zeros((np, np))\n",
        "    \n",
        "    for tri in topo:\n",
        "        K = Kmat(Dmat, tri, geom)\n",
        "        for i, idx in enumerate(tri):\n",
        "            for j, jdx in enumerate(tri):\n",
        "                Kglobal[idx, jdx] += K[i, j]\n",
        "    return Kglobal"
      ],
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2U5aWd3IR-fV"
      },
      "source": [
        "Dmat = array([[5.0, 0.0],\n",
        "              [0.0, 5.0]])"
      ],
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vzdr63-2R-fa"
      },
      "source": [
        "Kglobal = assemble_matrix(mesh, Dmat)"
      ],
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vr2RuCqGR-fe"
      },
      "source": [
        "The fully assembled left-hand side (LHS) looks like this:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IfraziGQR-ff",
        "outputId": "5174d09f-12e0-40c3-9f2c-3593a4c22ce0",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 85
        }
      },
      "source": [
        "print(Kglobal)"
      ],
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[[ 6.25 -1.25 -5.    0.  ]\n",
            " [-1.25  6.25  0.   -5.  ]\n",
            " [-5.    0.    6.25 -1.25]\n",
            " [ 0.   -5.   -1.25  6.25]]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RAMBIlFcR-fk"
      },
      "source": [
        "### Vector assembly (RHS)\n",
        "Now we need to perform RHS assembly on the vector $f$. This is similar to the matrix assembly above.\n",
        "\n",
        "The individual vector for each element is defined in `fvec()` and they are combined together in `assemble_vector()`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NezaFbmTR-fk"
      },
      "source": [
        "def fvec(s, p, geom):\n",
        "    assert(len(p) == 3)\n",
        "    x0, y0 = geometry[p[0]]\n",
        "    x1, y1 = geometry[p[1]]\n",
        "    x2, y2 = geometry[p[2]]\n",
        "\n",
        "    # Element area Ae\n",
        "    Ae = 0.5*abs((x0 - x1)*(y2 - y1) - (y0 - y1)*(x2 - x1))\n",
        "    return s*Ae/3.0*ones(3)\n",
        "\n",
        "def assemble_vector(mesh, s):\n",
        "    geom, topo = mesh\n",
        "    fg = zeros(len(geom))\n",
        "    for tri in topo:\n",
        "        f = fvec(s, tri, geom)\n",
        "        for i, idx in enumerate(tri):\n",
        "            fg[idx] += f[i]\n",
        "    return fg\n",
        "\n",
        "fglobal = assemble_vector(mesh, 3.0)"
      ],
      "execution_count": 10,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lkXPd6oPR-fp"
      },
      "source": [
        "### Boundary flux terms\n",
        "There are also some boundary flux terms in the example. Here we will add them directly to $f$."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kgpuqc2WR-fp"
      },
      "source": [
        "fglobal += array([0.0, -10.0, 0.0, -10.0])"
      ],
      "execution_count": 11,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7C6J0d1sR-fu"
      },
      "source": [
        "### Boundary conditions \n",
        "\n",
        "Dirichlet boundary conditions can be enforced by zeroing a row of the $K$ matrix, setting the diagonal entry to 1, and putting the desired value in the RHS $f$ vector"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "y9qEEUw2R-fv"
      },
      "source": [
        "def set_bc(K, f, row, val):\n",
        "    # Set BC\n",
        "    K[row] = 0.0\n",
        "    K[row, row] = 1.0\n",
        "    f[row] = val\n",
        "\n",
        "set_bc(Kglobal, fglobal, 0, 0.0)   \n",
        "set_bc(Kglobal, fglobal, 2, 0.0)"
      ],
      "execution_count": 12,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "I8fN0CkcR-f1"
      },
      "source": [
        "Finally, we can refer the $K.u = f$ problem to a linear algebra solver which will use an LU method to get the answer"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jflTNLivR-f2",
        "outputId": "81f83765-0a33-407d-e7ff-b794b1f518a0",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "u = linalg.solve(Kglobal, fglobal)\n",
        "print(u)\n"
      ],
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[-0.   -6.76  0.   -6.84]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xYzZ76vcR-f5"
      },
      "source": [
        ""
      ],
      "execution_count": 13,
      "outputs": []
    }
  ]
}