tracefem_scalar.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Scalar SurfacePDEs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Laplace-Beltrami problem\n",
    "We want to solve a problem of the form\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "          - \\Delta_{\\Gamma} u + u = & \\, f \n",
    "          & & \\text{on}~~ \\Gamma,          \n",
    "        \\end{aligned}\n",
    "$$\n",
    "\n",
    "with $\\Gamma := \\{ \\phi = 0 \\}$. Here $\\textbf{P} = \\textbf{I} - \\textbf{n}\\textbf{n}^T$ is the projection onto the tangential space and\n",
    "\n",
    "\\begin{equation*}\n",
    "\\nabla_{\\Gamma} u = \\textbf{P} \\nabla u, \\qquad \\operatorname{div}_{\\Gamma}u = \\operatorname{tr}(\\nabla_{\\Gamma}u)\n",
    "\\end{equation*}\n",
    "\n",
    "are the surface gradient for scalar-valued functions and the surface divergence operator. Moreover, the Laplace-Beltrami operator is defined as follows:\n",
    "\\begin{equation*}\n",
    "\\Delta_{\\Gamma} u = \\operatorname{div}_{\\Gamma} (\\nabla_{\\Gamma} u).\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# ngsolve stuff\n",
    "from ngsolve import *\n",
    "# visualization stuff\n",
    "from ngsolve.webgui import *\n",
    "# basic xfem functionality\n",
    "from xfem import *\n",
    "# basic geometry features (for the background mesh)\n",
    "from netgen.occ import *\n",
    "# pi\n",
    "from math import pi\n",
    "\n",
    "# material parameter:\n",
    "reac_coeff=1\n",
    "diff_coeff=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Geometry approximation\n",
    "We use a piecewise linear levelset approximation. \n",
    "\n",
    "First, we generate the background mesh of the domain and use a simplicial triangulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "cube = OCCGeometry(unit_cube.shape.Scale((0,0,0),3).Move((-1.5,-1.5,-1.5)), dim=3)\n",
    "mesh = Mesh (cube.GenerateMesh(maxh=0.2, quad_dominated=False))\n",
    "Draw(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "On the background mesh we define the level set function and do a P1-interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "phi = sqrt(x**2+y**2+z**2)-1\n",
    "lsetp1 = GridFunction(H1(mesh,order=1))\n",
    "InterpolateToP1(phi,lsetp1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## TraceFE / CutFE spaces\n",
    "For the discretization we use standard background FESpaces restricted to the tetrahedrons cut by $\\Gamma_h$:\n",
    "$$\n",
    "V_h^\\Gamma = V_h |_{\\Omega^{\\Gamma}_h},\n",
    "$$\n",
    "with $\\Omega^{\\Gamma}_h = \\cup \\{ T \\in \\mathcal{T}_h \\mid T \\cap \\Gamma_h \\neq \\emptyset \\}$. In NGSolve we will simply take the space $V_h$ but mark the irrelevant dofs using the CutInfo-class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# FESpace \n",
    "Vh = H1(mesh, order=1)\n",
    "\n",
    "ci = CutInfo(mesh, lsetp1)\n",
    "VhG = Compress(Vh,GetDofsOfElements(Vh,ci.GetElementsOfType(IF)))\n",
    "#VhG = Restrict(Vh,ci.GetElementsOfType(IF))\n",
    "freedofs = VhG.FreeDofs()\n",
    "\n",
    "u,v = VhG.TnT()\n",
    "\n",
    "gfu = GridFunction(VhG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Some helper coefficient functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# normal vector\n",
    "n = Normalize(grad(lsetp1))\n",
    "\n",
    "# mesh size parameter\n",
    "h = specialcf.mesh_size\n",
    "\n",
    "# define tangential projection\n",
    "P = lambda u: u - (u*n)*n\n",
    "\n",
    "# integration domains (and integration parameter \"subdivlvl\" and \"force_intorder\")\n",
    "ds = dCut(levelset=lsetp1, domain_type=IF, definedonelements=ci.GetElementsOfType(IF))\n",
    "dX = dx(definedonelements=ci.GetElementsOfType(IF))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## TraceFEM / CutFEM discretization\n",
    "Find $u_h \\in  V_h^{\\Gamma}$ such that\n",
    "\\begin{equation*}\n",
    "\\int_{\\Gamma_h} \\nabla_{\\Gamma_h} u_h \\cdot \\nabla_{\\Gamma_h} v_h \\, ds_h + \\int_{\\Gamma_h} u_h v_h \\, ds_h + \\rho \\int_{\\Omega_{h}^{\\Gamma}} (n_h \\cdot \\nabla u_h) (n_h \\cdot \\nabla v_h)  \\, dx = \\int_{\\Gamma_h} f_h v_h \\, ds_h ~~~ \\text{for all}~ v_h \\in V_h^{\\Gamma}\n",
    "\\end{equation*}\n",
    "with \n",
    "$$\n",
    "  n_h = \\nabla \\phi_h^{lin}/\\Vert \\nabla \\phi_h^{lin} \\Vert.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# expressions of test and trial functions:\n",
    "a = BilinearForm(VhG, symmetric = True)\n",
    "a += (diff_coeff * P(grad(u)) * P(grad(v)) + reac_coeff * u * v) * ds\n",
    "a += (diff_coeff/h+reac_coeff*h)*(grad(u)*n) * (grad(v)*n) * dX\n",
    "\n",
    "f_coeff = (sin(pi*z)*(diff_coeff*pi*pi*(1-z*z)+reac_coeff)+diff_coeff*cos(pi*z)*2*pi*z)\n",
    "f = LinearForm(VhG)\n",
    "f += f_coeff * v * ds\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a.Assemble()\n",
    "f.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can now solve the problem (recall that freedofs only marks relevant dofs):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu.vec[:] = 0.0\n",
    "gfu.vec.data = a.mat.Inverse(freedofs) * f.vec\n",
    "\n",
    "exact = sin(pi*z)\n",
    "l2error = sqrt( Integrate( (gfu-exact)**2 * ds.order(2), mesh=mesh))\n",
    "print (\"l2error : \", l2error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Visualization on the active mesh\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Draw(gfu, mesh, \"u\", clipping={\"x\": 0, \"y\": -0.6, \"z\": -0.5})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "well... not so nice."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Visualization with VTK\n",
    "\n",
    "Create vtk-output file to be read by Paraview or `pyvista`.\n",
    "To install pyvista, use \n",
    "```\n",
    "!PIP_BREAK_SYSTEM_PACKAGES=1 pip3 install pyvista vtk ipywidgets trame trame-vtk trame-vuetify --break-system-packages --upgrade\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "vtk = VTKOutput(ma=mesh,coefs=[lsetp1,gfu,exact],names=[\"P1-levelset\",\"u\",\"uSol\"],filename=\"laplacebeltrami3d\",subdivision=0)\n",
    "vtk.Do()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    import pyvista\n",
    "    pyvista.set_jupyter_backend('html')\n",
    "    visobj = pyvista.read('laplacebeltrami3d.vtu')\n",
    "    plot = pyvista.Plotter()\n",
    "    contour = visobj.contour([0], scalars=\"P1-levelset\", rng=[-1,1])\n",
    "    plot.add_mesh(contour, scalars=\"u\", cmap=\"jet\")\n",
    "    plot.show()\n",
    "except:\n",
    "    print(\"pyvista plotting not successful. No pyvista installed?\")\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "ls -l *.vtu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#%%bash\n",
    "#paraview laplacebeltrami3d.vtk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![alt](graphics/tracefem_scalar.png)"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}