Update documentation

_sources/preamble.ipynb

@@ -20,7 +20,7 @@
     "All the details on vectorial finite elements are found in Sabine Zaglmayr's Phd thesis from [here](http://asc.tuwien.ac.at/~schoeberl/wiki/index.php/Theses). A introduction for Hybrid Discontinuous Galerkin methods is Christoph Lehrenfeld's Master thesis, also from [here](http://asc.tuwien.ac.at/~schoeberl/wiki/index.php/Theses).\n",
     "\n",
     "\n",
-    "A recommended textbook for mixed finite element methods is [Mixed Finite  Element Methods and Applications](http://www.springer.com/us/book/9783642365188) by F. Brezzi, D. Boffi and M. Fortin.\n",
+    "A recommended textbook for mixed finite element methods is [Mixed Finite  Element Methods and Applications](http://www.springer.com/us/book/9783642365188) by D. Boffi, F. Brezzi, and M. Fortin.\n",
     "\n",
     "\n",
     "Joachim Schöberl, September 2017"
@@ -29,16 +29,14 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
+   "metadata": {},
    "outputs": [],
    "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -52,9 +50,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.7"
+   "version": "3.12.2"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }

_sources/primal/Untitled.ipynb

@@ -99,7 +99,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.2"
+   "version": "3.12.2"
   }
  },
  "nbformat": 4,

_sources/primal/boundary_conditions.ipynb

@@ -37,7 +37,7 @@
     "In the derivation of the weak form we had\n",
     "\n",
     "$$\n",
-    "\\int_\\Omega \\nabla u \\nabla v - \\int_{\\Gamma_D \\cup \\Gamma_N \\cup \\Gamma_R} \\frac{\\partial u }{\\partial n} v = \\int_\\Omega f v \\qquad \\forall \\, v \\in H^1(\\Omega)\n",
+    "\\int_\\Omega \\nabla u \\nabla v \\, dx - \\int_{\\Gamma_D \\cup \\Gamma_N \\cup \\Gamma_R} \\frac{\\partial u }{\\partial n} v \\, ds = \\int_\\Omega f v \\, dx \\qquad \\forall \\, v \\in H^1(\\Omega)\n",
     "$$"
    ]
   },
@@ -56,7 +56,7 @@
     "- \\int_{\\Gamma_R} \\alpha (u_0-u) v = \\int f v \\qquad \\forall \\, v \\in H^1(\\Omega)\n",
     "$$\n",
     "\n",
-    "We see that $\\int g v$ and $\\int \\alpha u_0 v$ are linear terms in the test-function, and thus belong to the linear-form $f$:\n",
+    "We see that $\\int g v$ and $\\int \\alpha u_0 v$ are linear terms in the test-function, and thus belong to the linear-form $f$ on the right hand side:\n",
     "\n",
     "$$\n",
     "\\text{find } u \\in H^1 : \\quad \\int_{\\Omega} \\nabla u \\nabla v + \\int_{\\Gamma_R} \\alpha u v = \\int_\\Omega f v + \\int_{\\Gamma_N} g v + \\int_{\\Gamma_R} \\alpha u_0 v\n",
@@ -89,25 +89,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from netgen.occ import unit_square\n",
     "mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))\n",
     "\n",
     "fes = H1(mesh, order=3)\n",
-    "gfu = GridFunction(fes)\n",
-    "a = BilinearForm(fes)\n",
-    "f = LinearForm (fes)\n",
-    "\n",
-    "u = fes.TrialFunction()\n",
-    "v = fes.TestFunction()\n",
+    "u, v = fes.TnT()\n",
     "\n",
     "alpha = 10\n",
     "u0 = 293\n",
-    "g = 1\n",
+    "g = 10\n",
     "\n",
-    "a += grad(u)*grad(v)*dx\n",
-    "a += alpha*u*v*ds(\"bottom\")\n",
-    "f += g * v*ds(\"right\")\n",
-    "f += alpha*u0*v*ds(\"bottom\")"
+    "a = BilinearForm(grad(u)*grad(v)*dx + alpha*u*v*ds(\"bottom\"))\n",
+    "f = LinearForm (g * v*ds(\"right\") + alpha*u0*v*ds(\"bottom\"))"
    ]
   },
   {
@@ -119,6 +111,7 @@
     "# assemble and solve\n",
     "a.Assemble()\n",
     "f.Assemble()\n",
+    "gfu = GridFunction(fes)\n",
     "gfu.vec.data = a.mat.Inverse() * f.vec"
    ]
   },
@@ -128,8 +121,32 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Draw (gfu, mesh, \"temperature\")\n",
-    "Draw (grad(gfu), mesh, \"gradient\");"
+    "The temperature:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Draw (gfu, mesh); "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The heat flux:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Draw (-grad(gfu), mesh, vectors=True);"
    ]
   },
   {
@@ -139,9 +156,10 @@
     "Essential boundary conditions\n",
     "----\n",
     "Essential boundary conditions are less natural: We have to set the solution field to the given Dirichlet values, and restrict the test-functions to 0 on the Dirichlet boundary:\n",
+    "$\\text{find } u \\in H^1, u = u_D \\text{ on } \\Gamma_D \\; \\text{such that}$\n",
     "\n",
     "$$\n",
-    "\\text{find } u \\in H^1, u = u_D \\text{ on } \\Gamma_D \\text{ s.t. } A(u,v) = f(v) \\quad \\forall \\, v \\in H^1, v = 0 \\text{ on } \\Gamma_D\n",
+    "A(u,v) = f(v) \\quad \\forall \\, v \\in H^1, v = 0 \\text{ on } \\Gamma_D\n",
     "$$\n",
     "\n",
     "We split the solution vector into two parts: The given coefficients on the Dirichlet boundary, and all other including internal coefficients and coefficients on the natural boundaries:\n",
@@ -155,7 +173,7 @@
     "$$\n",
     "A = \\left( \\begin{array}{cc} A_{DD} & A_{Df} \\\\ A_{fD} & A_{ff} \\end{array} \\right)\n",
     "\\qquad\n",
-    "f = \\left( \\begin{array}{c} f_D \\\\ f_f \\end{array} \\right)\n",
+    "f = \\left( \\begin{array}{c} f_D \\\\ f_f \\end{array} \\right).\n",
     "$$\n",
     "\n",
     "The test functions are reduced to the free nodes. Thus, the equations to solve are\n",
@@ -231,25 +249,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from netgen.occ import unit_square\n",
     "mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))\n",
     "\n",
     "fes = H1(mesh, order=3, dirichlet=\"bottom|left|right\")\n",
+    "u,v = fes.TnT()\n",
     "gfu = GridFunction(fes)\n",
-    "a = BilinearForm(fes)\n",
-    "f = LinearForm (fes)\n",
-    "\n",
-    "u = fes.TrialFunction()\n",
-    "v = fes.TestFunction()\n",
-    "\n",
-    "uD = sin(2*pi*x)\n",
-    "g = 1\n",
-    "\n",
-    "a += grad(u)*grad(v)*dx\n",
-    "\n",
-    "# assemble and solve\n",
-    "a.Assemble()\n",
-    "f.Assemble();"
+    "a = BilinearForm(grad(u)*grad(v)*dx).Assemble()\n",
+    "f = LinearForm (fes)  # nothing here"
    ]
   },
   {
@@ -265,8 +271,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "uD = sin(2*pi*x)\n",
     "gfu.Set (uD, definedon=mesh.Boundaries(\"bottom\"))\n",
-    "scene = Draw(gfu)"
+    "Draw(gfu);"
    ]
   },
   {
@@ -282,10 +289,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "r = f.vec - a.mat * gfu.vec\n",
-    "gfu.vec.data += a.mat.Inverse(fes.FreeDofs()) * r\n",
+    "inv = a.mat.Inverse(fes.FreeDofs())\n",
+    "\n",
+    "residuum = f.vec - a.mat * gfu.vec\n",
+    "gfu.vec.data += inv * residuum\n",
     "\n",
-    "scene.Redraw()"
+    "Draw (gfu);"
    ]
   },
   {