08-fe-assembly-template.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e8a88f1-21b1-4ccb-b92d-924697f00ab9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import sparse\n",
    "from scipy.sparse import linalg as sla\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import tri\n",
    "\n",
    "import plotly.offline as py\n",
    "import plotly.figure_factory as ff\n",
    "import plotly.io as pio\n",
    "pio.renderers.default = 'iframe'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2281c06-8d2f-4d78-8d91-f7ee4c42a45c",
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://raw.githubusercontent.com/lukeolson/cs555-demos/main/mesh.e\n",
    "!wget https://raw.githubusercontent.com/lukeolson/cs555-demos/main/mesh.v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cce58e0a-6a3f-40ca-8d49-aeb05052819a",
   "metadata": {},
   "outputs": [],
   "source": [
    "V = np.loadtxt('mesh.v')\n",
    "E = np.loadtxt('mesh.e', dtype=int)\n",
    "t = tri.Triangulation(V[:, 0], V[:, 1], E)\n",
    "mesh = tri.UniformTriRefiner(t)\n",
    "t = mesh.refine_triangulation(subdiv=1)\n",
    "\n",
    "X = t.x\n",
    "Y = t.y\n",
    "E = t.triangles\n",
    "#stats = tri.TriAnalyzer(triangles)\n",
    "plt.triplot(X, Y, E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70971cc7-ae22-4f3e-b5aa-8a2bcdb0a284",
   "metadata": {},
   "outputs": [],
   "source": [
    "def kappa(x, y):\n",
    "    if np.sqrt(x**2 + y**2) <= 0.25:\n",
    "        return 100.0\n",
    "    return 2.0\n",
    "\n",
    "def f(x, y):\n",
    "    if np.sqrt(x**2 + y**2) <= 0.25:\n",
    "        return 100.0\n",
    "    return 0.0\n",
    "\n",
    "def g(x, y):\n",
    "    return 0*x + 0.0\n",
    "\n",
    "def g1(x, y):\n",
    "    return 2.0 * np.cos(np.pi * y / 2)\n",
    "\n",
    "def g2(x, y):\n",
    "    return (1 + x) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ade85cc-9496-4ba6-b267-3881d0d7fbfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "tol = 1e-14\n",
    "id1 = np.where(abs(Y+1) < tol)[0]\n",
    "id2 = np.where(abs(Y-1) < tol)[0]\n",
    "id3 = np.where(abs(X+1) < tol)[0]\n",
    "id4 = np.where(abs(X-1) < tol)[0]\n",
    "bc = [{'id': id1, 'g': g2},\n",
    "      {'id': id2, 'g': g2},\n",
    "      {'id': id3, 'g': g1},\n",
    "      ]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d888e741-b07e-41a1-865e-d4954884ecd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ne = E.shape[0]\n",
    "\n",
    "# allocate sparse matrix arrays\n",
    "m = 3  # for linears\n",
    "AA = np.zeros((ne, m**2))\n",
    "IA = np.zeros((ne, m**2), dtype=int)\n",
    "JA = np.zeros((ne, m**2), dtype=int)\n",
    "bb = np.zeros((ne, m))\n",
    "ib = np.zeros((ne, m), dtype=int)\n",
    "jb = np.zeros((ne, m), dtype=int)\n",
    "\n",
    "# Assemble A and b\n",
    "for ei in range(0, ne):\n",
    "    # Step 1: set the vertices and indices\n",
    "    K = E[ei, :]\n",
    "    x0, y0 = X[K[0]], Y[K[0]]\n",
    "    x1, y1 = X[K[1]], Y[K[1]]\n",
    "    x2, y2 = X[K[2]], Y[K[2]]\n",
    "\n",
    "    # Step 2: compute the Jacobian, inv, and det\n",
    "    J = np.array([[x1 - x0, x2 - x0],\n",
    "                  [y1 - y0, y2 - y0]])\n",
    "    invJ = np.linalg.inv(J.T)\n",
    "    detJ = np.linalg.det(J)\n",
    "\n",
    "    # Step 3, define the gradient of the basis\n",
    "    dbasis = np.array([[-1, 1, 0],\n",
    "                       [-1, 0, 1]])\n",
    "\n",
    "    # Step 4\n",
    "    dphi = invJ.dot(dbasis)\n",
    "\n",
    "    # Step 5, 1-point gauss quadrature\n",
    "    Aelem = # TODO\n",
    "\n",
    "    # Step 6, 1-point gauss quadrature\n",
    "    belem = # TODO\n",
    "\n",
    "    # Step 7\n",
    "    AA[ei, :] = # TODO\n",
    "    IA[ei, :] = # TODO\n",
    "    JA[ei, :] = # TODO\n",
    "    bb[ei, :] = belem.ravel()\n",
    "    ib[ei, :] = K[np.arange(m)]\n",
    "    jb[ei, :] = 0\n",
    "\n",
    "# convert matrices\n",
    "A = sparse.coo_matrix((AA.ravel(), (IA.ravel(), JA.ravel())))\n",
    "A.sum_duplicates()\n",
    "b = sparse.coo_matrix((bb.ravel(), (ib.ravel(), jb.ravel()))).toarray().ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24a74e99-0064-497a-a5a5-ae79496b6ff3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set BC in the right hand side\n",
    "# set the lifting function (1 of 3)\n",
    "u0 = np.zeros((A.shape[0],))\n",
    "for c in bc:\n",
    "    idx = c['id']\n",
    "    g = c['g']\n",
    "    u0[idx] = g(X[idx], Y[idx])\n",
    "\n",
    "# lift (2 of 3)\n",
    "b = b - A * u0\n",
    "\n",
    "# fix the values (3 of 3)\n",
    "for c in bc:\n",
    "    idx = c['id']\n",
    "    b[idx] = u0[idx]\n",
    "\n",
    "# set BC to identity in the matrix\n",
    "# collect all BC indices (1 of 2)\n",
    "Dflag = np.full((A.shape[0],), False)\n",
    "for c in bc:\n",
    "    idx = c['id']\n",
    "    Dflag[idx] = True\n",
    "\n",
    "# write identity (2 of 2)\n",
    "for k in range(0, len(A.data)):\n",
    "    i = A.row[k]\n",
    "    j = A.col[k]\n",
    "    if Dflag[i] or Dflag[j]:\n",
    "        if i == j:\n",
    "            A.data[k] = 1.0\n",
    "        else:\n",
    "            A.data[k] = 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71270925-5b78-414c-9760-cc8509677e64",
   "metadata": {},
   "outputs": [],
   "source": [
    "u = sla.spsolve(A.tocsr(), b)  # see CS556"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e66c6bf8-4429-403a-a341-d52f64a01669",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = ff.create_trisurf(x=X, y=Y, z=u,\n",
    "                        simplices=E)\n",
    "py.iplot(fig, filename='surf')\n",
    "fig.show(renderer=\"colab\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8a80480-ec06-4565-ab45-38f569959bcb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80fd9965-608f-4989-8370-bc785612a938",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}