{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.6 Scalar PDE on surfaces (with HDG)\n",
    "We are solving the scalar linear transport problem\n",
    "\n",
    "$$\n",
    "  \\partial_t u - \\varepsilon \\Delta_\\Gamma u + \\operatorname{div}_\\Gamma ( \\mathbf{w} u ) = 0 \\quad \\text{ on } \\Gamma,\n",
    "$$\n",
    "\n",
    "where $\\Gamma$ is a closed oriented surface in $\\mathbb{R}^3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Mesh of (only) the surface:\n",
    "We consider the unit sphere (only surface mesh!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "sphere = Sphere((0,0,0),1).faces[0]\n",
    "sphere.name = \"sphere\"\n",
    "geo = OCCGeometry(sphere)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.5))\n",
    "Draw(mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "order = 4\n",
    "if order > 0:\n",
    "    mesh.Curve(order)\n",
    "Draw(mesh)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Tangential velocity field:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "t = Parameter(0.)\n",
    "b = CoefficientFunction((y,-x,0))\n",
    "Draw (b, mesh, \"b\")\n",
    "eps = 5e-5"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### some local quantities\n",
    "normal, tangential, co-normal and mesh size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "n = specialcf.normal(mesh.dim)\n",
    "h = specialcf.mesh_size\n",
    "local_tang = specialcf.tangential(mesh.dim)\n",
    "con = Cross(n,local_tang) #co-normal pointing outside of a surface element\n",
    "bn = InnerProduct(b,con)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Discretization"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Surface `FESpaces`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fesl2 = SurfaceL2(mesh, order=order)\n",
    "fesedge = FacetSurface(mesh, order=order)\n",
    "\n",
    "fes = fesl2*fesedge\n",
    "u,uE = fes.TrialFunction()\n",
    "v,vE = fes.TestFunction()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### surface gradients, co-normal derivatives and jumps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gradu = u.Trace().Deriv()\n",
    "gradv = v.Trace().Deriv()\n",
    "jumpu = u - uE\n",
    "jumpv = v - vE\n",
    "dudn = InnerProduct(gradu,con)\n",
    "dvdn = InnerProduct(gradv,con)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Diffusion formulation (Hybrid Interior Penalty on Surface):\n",
    "\n",
    "$$\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 u (u-\\widehat u)\n",
    "+ \\frac{\\alpha p^2}{h} \\sum_F \\int_F (u-\\widehat u)(v-\\widehat v)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a = BilinearForm(fes)\n",
    "# diffusion\n",
    "\n",
    "dr = ds(element_boundary=True)\n",
    "\n",
    "alpha = 10*(order+1)**2\n",
    "a += eps*gradu * gradv * ds\n",
    "a += - eps*dudn * jumpv * dr\n",
    "a += - eps*dvdn * jumpu * dr\n",
    "a += eps*alpha/h * jumpu * jumpv * dr\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Convection formulation (Hybrid Upwinding):\n",
    "\n",
    "$$\n",
    "  \\sum_T \\int_T - u \\mathbf{w} \\cdot \\nabla v\n",
    "+ \\sum_T \\int_{\\partial T} (\\mathbf{w} \\cdot \\mathbf{n}) u^{hyb-upw} v \n",
    "+ \\sum_T \\int_{\\partial T_{out}} (\\mathbf{w} \\cdot \\mathbf{n}) (\\hat{u} - u) \\hat{v}\n",
    "$$\n",
    "with\n",
    "$$\n",
    "    u^{hyb-upw} = \\left\\{ \\begin{array}{cll} u & \\text{if } \\mathbf{w} \\cdot \\mathbf{n} & (outflow),  \\\\ \\hat{u} & \\text{else}  & (inflow) . \\end{array} \\right.\n",
    "$$\n",
    "\n",
    "|  outflow   | inflow     |\n",
    "|------|------|\n",
    "|   ![elementupwind](elementupwind.png)  | ![facetupwind](facetupwind.png) |\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a = BilinearForm(fes)\n",
    "# convection (surface version of Egger+Schöberl formulation):\n",
    "a += - u * InnerProduct(b,gradv) * ds\n",
    "a += IfPos(bn, bn*u, bn*uE) * v * dr + IfPos(bn, bn, 0) * (uE - u) * vE * dr\n",
    "a.Assemble()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Mass matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "m = BilinearForm(fes)\n",
    "m += u*v*ds\n",
    "m.Assemble()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### $M^\\ast$-matrix for implicit Euler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "dt = 0.02\n",
    "\n",
    "mstar = m.mat.CreateMatrix()\n",
    "mstar.AsVector().data = m.mat.AsVector() + dt * a.mat.AsVector()\n",
    "mstarinv = mstar.Inverse()\n",
    "\n",
    "from math import pi\n",
    "T = 8*pi"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Time stepping"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Visualization of solution (Deformation in normal direction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "gfu.components[0].Set(1.5*exp(-20*(x*x+(y-1)**2+(z)**2)),definedon=mesh.Boundaries(\"sphere\"))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The time loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "sceneun = Draw (gfu.components[0]*n, mesh, \"un\", deformation = True, autoscale=False)\n",
    "t.Set(0)\n",
    "for i in range(int(T/dt)):\n",
    "    gfu.vec.data = mstarinv @ m.mat * gfu.vec\n",
    "    sceneun.Redraw()\n",
    "    t.Set(t.Get()+dt)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Remark:\n",
    "* We can apply static condensation as in the volume"
   ]
  }
 ],
 "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}