{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# 9.2 Implement our own system assembling\n",
    "\n",
    "In this tutorial we\n",
    "\n",
    "* write integrators for $\\int_T f v dx$ and $\\int_T \\nabla u \\nabla v dx$: <br>\n",
    "[myIntegrator.hpp](myIntegrator.hpp)\n",
    "[myIntegrator.cpp](myIntegrator.cpp)\n",
    "\n",
    "\n",
    "* put together element matrices to the global system matrix: <br>\n",
    "[myAssembling.cpp](myAssembling.cpp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "import pyngcore   # for timers \n",
    "\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import unit_square\n",
    "\n",
    "from ngsolve.fem import CompilePythonModule\n",
    "from pathlib import Path\n",
    "\n",
    "m = CompilePythonModule(Path('myassemblemodule.cpp'), init_function_name='mymodule')\n",
    "dir (m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))\n",
    "fes = H1(mesh, order=1, dirichlet=\".*\")\n",
    "u, v = fes.TnT()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "### use our own integrators for element matrix calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = BilinearForm(fes)\n",
    "a += m.MyLaplace(1)\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += m.MySource(x)\n",
    "f.Assemble();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec\n",
    "Draw (gfu);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "### we can inspect the element matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "mylap = m.MyLaplace(1)\n",
    "ei = ElementId(VOL,17)\n",
    "mylap.CalcElementMatrix(fes.GetFE(ei), mesh.GetTrafo(ei))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "### use our own matrix assembling function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "for l in range(4): mesh.Refine()\n",
    "fes.Update()\n",
    "gfu.Update()\n",
    "f.Assemble();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "pyngcore.ResetTimers()\n",
    "print (\"num elements =\", mesh.ne, \", ndof =\", fes.ndof)\n",
    "with TaskManager(pajetrace=10**8):\n",
    "    # using our integrator\n",
    "    mymatrix = m.MyAssembleMatrix(fes, m.MyLaplace(CF(1)), parallel=False)\n",
    "    \n",
    "    # using NGSolve built-in symbolic integrator\n",
    "    # mymatrix = m.MyAssembleMatrix(fes, (grad(u)*grad(v)*dx)[0].MakeBFI(), parallel=True)\n",
    "\n",
    "# print (\"my matrix = \", mymat)\n",
    "\n",
    "if fes.ndof < 100000:\n",
    "    gfu.vec.data = mymatrix.Inverse(fes.FreeDofs()) * f.vec\n",
    "    Draw (gfu);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "### Compare timings for \n",
    "\n",
    "  * `MyLaplace` and `grad(u)*grad(v)*dx`\n",
    "  * sequential and parallel\n",
    "  * for different polynomial orders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in pyngcore.Timers(): \n",
    "    if t[\"name\"].startswith(\"MyAssemble\"): \n",
    "        print (t[\"name\"],\"   \", t[\"time\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "### Exercise:\n",
    "\n",
    "\n",
    "   * Implement `MyAssembleVector` for building the right hand side vector for  $\\int_\\Omega f v \\, dx$\n",
    "   * Implement `MyNeumannIntegrator` for evaluating $\\int_{\\partial \\Omega} g v \\, ds$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "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.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}