{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6.1.2 Surface PDEs\n",
    "## Surface Poisson equation\n",
    "### Homogeneous Dirichlet data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from netgen.csg import *\n",
    "from netgen.meshing import MeshingStep\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a half-sphere we want to solve the surface Poisson equation $$\\int_S\\nabla_{\\Gamma}u\\cdot\\nabla_{\\Gamma}v\\,ds = \\int_S fv\\,ds,\\qquad u=0 \\text{ on } \\Gamma_D$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo          = CSGeometry()\n",
    "sphere       = Sphere(Pnt(0,0,0), 1)\n",
    "bot          = Plane(Pnt(0,0,0), Vec(0,0,-1))\n",
    "finitesphere = sphere * bot\n",
    "\n",
    "geo.AddSurface(sphere, finitesphere.bc(\"surface\"))\n",
    "geo.NameEdge(sphere,bot, \"bottom\")\n",
    "\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.3))\n",
    "mesh.Curve(2)\n",
    "Draw(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, we define the usual $H^1$ finite element space and use the <i>dirichlet_bbnd</i> flag to indicate the BBoundary on which the Dirichlet conditions are prescribed. The test- and trial-functions are given as usual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=2, dirichlet_bbnd=\"bottom\")\n",
    "u, v = fes.TnT()\n",
    "print(fes.FreeDofs())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the (bi-)linear form we have to take care that we have to define the according integrators on the boundary and that the <i>Trace</i> operator has to be used to obtain the tangential/surface derivative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = BilinearForm(fes, symmetric=True)\n",
    "a += grad(u).Trace()*grad(v).Trace()*ds\n",
    "a.Assemble()\n",
    "\n",
    "force = sin(x)*y*exp(z)\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += force*v*ds\n",
    "f.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving is done as usual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "gfu.vec.data = a.mat.Inverse(fes.FreeDofs())*f.vec\n",
    "Draw(gfu, mesh, \"u\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inhomogeneous Dirichlet data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To solve the same problem with non-homogenous Dirichlet data, $u=u_D$ on $\\Gamma_D$ the same technique as in the volume case is used, where we have to set a function on the BBoundary instead on the boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu.Set(x, definedon=mesh.BBoundaries(\"bottom\"))\n",
    "r = f.vec.CreateVector()\n",
    "r.data = f.vec - a.mat*gfu.vec\n",
    "gfu.vec.data += a.mat.Inverse(fes.FreeDofs())*r\n",
    "Draw(gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finite element spaces for surfaces\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We differ between two types of finite element spaces for surfaces. The first class consists of spaces, where the restriction of a 3D element to the surface leads to a valid 2D element of the same type: H1, HCurl, HCurlCurl, NumberSpace.\n",
    "\n",
    "These spaces can directly be used, one has to take care using the Trace operator. Otherwise an exception is thrown during assembling. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "a += grad(u)*grad(v)*ds\n",
    "a.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The NumberSpace is an exception as it represents only a number, where no Trace operator has to be used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second class is given by\n",
    "\n",
    "| Space     |  Surface Space |\n",
    "|:----------|:---------------|\n",
    "| L2 \t       | SurfaceL2\n",
    "| HDiv         | HDivSurface\n",
    "| HDivDiv \t   | HDivDivSurface\n",
    "| FacetFESpace | FacetSurface\n",
    "\n",
    "Here, a 2D reference element is mapped directly onto the surface. To be consistent, also here the Trace operator has to be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}