NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(advanced:hdg)=\n",
    "## Hybrid Discontinuous Galerkin methods\n",
    "use additionally the *hybrid* facet variable on the skeleton:\n",
    "\n",
    "$$\n",
    "A(u,\\widehat u; v, \\widehat v) = \n",
    "  \\sum_T \\int_T \\nabla u \\nabla v\n",
    "- \\sum_T \\int_{\\partial T} n \\nabla u (v-\\widehat v)\n",
    "- \\sum_T \\int_{\\partial T} n \\nabla v (u-\\widehat u)\n",
    "+ \\frac{\\alpha p^2}{h} \\sum_F \\int_F (u-\\widehat u)(v-\\widehat v)\n",
    "$$\n",
    "\n",
    "the jump-term is now replaced by the difference $u - \\widehat u$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No additional matrix entries across elements are produced. Dirichlet boundary conditions are set as usual to the facet variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *\n",
    "\n",
    "def AtomicShape(maxh = 0.1):\n",
    "    ax2d = gp_Ax2d(gp_Pnt2d(0.5, 0.5), gp_Dir2d(1, 0))\n",
    "    outeli = Ellipse( ax2d, 0.5, 0.1)\n",
    "    ineli = Ellipse( ax2d, 0.3, 0.06)\n",
    "\n",
    "    shape = outeli.Face() - ineli.Face() \n",
    "\n",
    "\n",
    "    shape += shape.Rotate(gp_Ax1(gp_Pnt(0.5, 0.5, 0), gp_Dir(0,0,1)), 120)\n",
    "    shape += shape.Rotate(gp_Ax1(gp_Pnt(0.5, 0.5, 0), gp_Dir(0,0,1)), 240)\n",
    "    shape += Circle((0.5,0.5),0.2).Face() \n",
    "\n",
    "    Draw(shape)\n",
    "    geo = OCCGeometry(shape, dim = 2)\n",
    "    return Mesh(geo.GenerateMesh(maxh = maxh)).Curve(3)\n",
    "\n",
    "mesh = AtomicShape()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order=2\n",
    "V = L2(mesh, order=order)\n",
    "F = FacetFESpace(mesh, order=order, dirichlet=\".*\")\n",
    "fes = V*F\n",
    "u,uhat = fes.TrialFunction()\n",
    "v,vhat = fes.TestFunction()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the jump is the difference between element-term and facet-term:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jump_u = u-uhat\n",
    "jump_v = v-vhat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_centr = (x-0.5)**2+(y-0.5)**2\n",
    "\n",
    "f = LinearForm(100*exp(10*r_centr)*v*dx).Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 4\n",
    "condense = False\n",
    "h = specialcf.mesh_size\n",
    "n = specialcf.normal(mesh.dim)\n",
    "\n",
    "a = BilinearForm(fes, condense=condense)\n",
    "dS = dx(element_boundary=True)\n",
    "a += grad(u)*grad(v)*dx + \\\n",
    "    alpha*order**2/h*jump_u*jump_v*dS + \\\n",
    "    (-grad(u)*n*jump_v - grad(v)*n*jump_u)*dS\n",
    "\n",
    "b = CF( (20,1) )\n",
    "uup = IfPos(b*n, u, uhat)\n",
    "a += -b * u * grad(v)*dx + b*n*uup*jump_v *dS\n",
    "a.Assemble()\n",
    "\n",
    "#f = LinearForm(1*v*dx).Assemble()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "print (\"A non-zero elements:\", a.mat.nze)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not condense:\n",
    "    inv = a.mat.Inverse(fes.FreeDofs(), \"pardiso\")\n",
    "    gfu.vec.data = inv * f.vec\n",
    "else:\n",
    "    fmod = (f.vec + a.harmonic_extension_trans * f.vec).Evaluate()\n",
    "    \n",
    "    inv = a.mat.Inverse(fes.FreeDofs(True), \"pardiso\")\n",
    "    gfu.vec.data = inv * fmod\n",
    "    \n",
    "    gfu.vec.data += a.harmonic_extension * gfu.vec\n",
    "    gfu.vec.data += a.inner_solve * f.vec\n",
    "\n",
    "Draw (gfu.components[0], mesh, \"u-HDG\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remarks on sparsity pattern in NGSolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.sparse as sp\n",
    "import matplotlib.pylab as plt\n",
    "A = sp.csr_matrix(a.mat.CSR())\n",
    "plt.spy(A, markersize=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Remark 1: The sparsity pattern is set up a-priorily\n",
    "* The sparsity pattern of a sparse matrix in NGSolve is independent of its entries (it's set up a-priorily). \n",
    "* We can have \"nonzero\" entries that have the value 0\n",
    "\n",
    "Below we show the reserved memory for the sparse matrix and the (numerically) non-zero entries in this sparse matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes2 = L2(mesh, order=order, dgjumps=True)\n",
    "u,v=fes2.TnT()\n",
    "a3 = BilinearForm(fes2)\n",
    "a3 += u*v*dx + (u+u.Other())*v*dx(skeleton=True)\n",
    "a3.Assemble();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.sparse as sp\n",
    "import matplotlib.pylab as plt\n",
    "plt.rcParams['figure.figsize'] = (12, 12)\n",
    "A = sp.csr_matrix(a3.mat.CSR())\n",
    "fig = plt.figure(); ax1 = fig.add_subplot(121); ax2 = fig.add_subplot(122)\n",
    "ax1.set_xlabel(\"numerically non-zero\"); ax1.spy(A)\n",
    "ax2.set_xlabel(\"reserved entries (potentially non-zero)\"); ax2.spy(A,precision=-1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Remark 2: Sparsity pattern with and without `dgjumps=True` is different"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a1 = BilinearForm(L2(mesh, order=order, dgjumps=False)); a1.Assemble()\n",
    "a2 = BilinearForm(L2(mesh, order=order, dgjumps=True)); a2.Assemble()\n",
    "A1 = sp.csr_matrix(a1.mat.CSR())\n",
    "A2 = sp.csr_matrix(a2.mat.CSR())\n",
    "fig = plt.figure(); ax1 = fig.add_subplot(121); ax2 = fig.add_subplot(122)\n",
    "ax1.set_xlabel(\"dgjumps=False\"); ax1.spy(A1,precision=-1)\n",
    "ax2.set_xlabel(\"dgjumps=True\"); ax2.spy(A2,precision=-1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Remark 3: Dof numbering of higher order FESpaces "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "* In `NGSolve` `FESpace`s typically have a numbering where the first block of dofs corresponds to a low order subspace (which is convenient for iterative solvers). \n",
    "* For L2 this means that the first dofs correspond to the constants on elements. \n",
    "\n",
    "* You can turn this behavior off for some spaces, e.g. for L2 by adding the flag `all_dofs_together`.\n",
    "\n",
    "We demonstrate this in the next comparison:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (15, 15)\n",
    "fig = plt.figure()\n",
    "ax = [fig.add_subplot(131), fig.add_subplot(132), fig.add_subplot(133)]\n",
    "\n",
    "for i, (order, all_dofs_together, label) in enumerate([(0,False, \"non-zeros (p=0)\"),\n",
    "                                                    (1,False,\"non-zeros (p=1, low order + high order)\"),\n",
    "                                                    (1,True,\"non-zeros (p=1, all_dofs_together)\")]):\n",
    "    a = BilinearForm(L2(mesh,order=order,dgjumps=True,all_dofs_together=all_dofs_together))\n",
    "    a.Assemble()\n",
    "    ax[i].spy(sp.csr_matrix(a.mat.CSR()),markersize=3,precision=-1)\n",
    "    ax[i].set_xlabel(label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}