{
"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$:
\n",
"[myIntegrator.hpp](myIntegrator.hpp)\n",
"[myIntegrator.cpp](myIntegrator.cpp)\n",
"\n",
"\n",
"* put together element matrices to the global system matrix:
\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
}