{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 8.3 Unfitted FEM discretizations\n",
    "\n",
    "We want to solve a geometrically unfitted model problem for a *stationary* domain.\n",
    "\n",
    "We use a level set description (cf. [basics.ipynb](basics.ipynb)): \n",
    "\n",
    "$$\n",
    "  \\Omega_{-} := \\{ \\phi < 0 \\}, \\quad\n",
    "  \\Omega_{+} := \\{ \\phi > 0 \\}, \\quad\n",
    "  \\Gamma := \\{ \\phi = 0 \\}.\n",
    "$$\n",
    "\n",
    "and use a piecewise linear approximation as *a starting point* in the discretization (cf. [intlset.ipynb](intlset.ipynb) for a discussion of geometry approximations). \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We first import the related packages: \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# the constant pi\n",
    "from math import pi\n",
    "# ngsolve stuff\n",
    "from ngsolve import *\n",
    "# basic xfem functionality\n",
    "from xfem import *\n",
    "# basic geometry features (for the background mesh)\n",
    "from netgen.geom2d import SplineGeometry\n",
    "# visualization stuff\n",
    "from ngsolve.webgui import * "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "![bubble](graphics/bubble_coarse.png)\n",
    "\n",
    "## Interface problem \n",
    "\n",
    "We want to solve a problem of the form: \n",
    "\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "          - \\nabla \\cdot (\\alpha_{\\pm} \\nabla u) = & \\, f \n",
    "          & & \\text{in}~~ \\Omega_{\\pm}, \n",
    "          \\\\\n",
    "          [\\![u]\\!] = & \\, 0 \n",
    "          &  & \\text{on}~~ \\Gamma, \n",
    "          \\\\\n",
    "          [\\![-\\alpha \\nabla u \\cdot \\mathbf{n}]\\!]   = & \\, 0 \n",
    "          &  & \\text{on}~~ \\Gamma,\n",
    "          \\\\\n",
    "          u = & \\, u_D\n",
    "          &  & \\text{on}~~ \\partial \\Omega.\n",
    "        \\end{aligned} \\right.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "square = SplineGeometry()\n",
    "square.AddRectangle([-1.5,-1.5],[1.5,1.5],bc=1)\n",
    "mesh = Mesh (square.GenerateMesh(maxh=0.4, quad_dominated=False))\n",
    "\n",
    "levelset = (sqrt(x*x+y*y) - 1.0)\n",
    "DrawDC(levelset, -3.5, 2.5, mesh,\"levelset\")\n",
    "\n",
    "lsetp1 = GridFunction(H1(mesh,order=1))\n",
    "InterpolateToP1(levelset,lsetp1)\n",
    "DrawDC(lsetp1, -3.5, 2.5, mesh, \"lsetp1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Cut FE spaces\n",
    "\n",
    "For the discretization we use standard background FESpaces restricted to the subdomains: \n",
    "\n",
    "$$\n",
    "V_h^\\Gamma \\quad = \\qquad V_h |_{\\Omega_-^{lin}} \\quad \\oplus \\quad V_h |_{\\Omega_+^{lin}}\n",
    "$$\n",
    "\n",
    "| composed | inner | outer |\n",
    "|-|-|-|\n",
    "| ![spaceboth](graphics/SpaceBoth.png) | ![spaceinner](graphics/SpaceInner.png) | ![spaceouter](graphics/SpaceOuter.png) |\n",
    "\n",
    "In NGSolve we simply take the product space $V_h \\times V_h$ but mark the irrelevant dofs using the CutInfo-class:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "Vh = H1(mesh, order=2, dirichlet=\".*\")\n",
    "VhG = FESpace([Vh,Vh])\n",
    "\n",
    "ci = CutInfo(mesh, lsetp1)\n",
    "freedofs = VhG.FreeDofs()\n",
    "freedofs &= CompoundBitArray([GetDofsOfElements(Vh,ci.GetElementsOfType(HASNEG)),\n",
    "                              GetDofsOfElements(Vh,ci.GetElementsOfType(HASPOS))])\n",
    "\n",
    "gfu = GridFunction(VhG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Let us visualize active dofs:\n",
    "* active dofs of first space are set to -1\n",
    "* active dofs of second space are set to 1\n",
    "* inactive dofs are 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu.components[0].Set(CoefficientFunction(-1))\n",
    "gfu.components[1].Set(CoefficientFunction(1))\n",
    "for i, val in enumerate(freedofs):\n",
    "    if not val:\n",
    "        gfu.vec[i] = 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "Draw(gfu.components[0], mesh, \"background_uneg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "Draw(gfu.components[1], mesh, \"background_upos\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Only the parts which are in the corresponding subdomain are relevant. The solution $u$ is: \n",
    "\n",
    "$$\n",
    " u = \\left\\{ \\begin{array}{cc} u_- & \\text{ if } {\\phi}_h^{lin} < 0, \\\\ u_+ & \\text{ if } {\\phi}_h^{lin} \\geq 0. \\end{array} \\right.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "DrawDC(lsetp1, gfu.components[0], gfu.components[1], mesh, \"u\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Improvement: use `Compress` to reduce unused dofs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Vh = H1(mesh, order=2, dirichlet=[1,2,3,4])\n",
    "ci = CutInfo(mesh, lsetp1)\n",
    "VhG = FESpace([Compress(Vh,GetDofsOfElements(Vh,ci.GetElementsOfType(cdt))) for cdt in [HASNEG,HASPOS]])\n",
    "\n",
    "freedofs = VhG.FreeDofs()\n",
    "gfu = GridFunction(VhG)\n",
    "gfu.components[0].Set(1)\n",
    "gfu.components[1].Set(-1)\n",
    "DrawDC(lsetp1, gfu.components[0], gfu.components[1], mesh, \"u\")\n",
    "print(Vh.ndof, VhG.components[0].ndof, VhG.components[1].ndof)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Nitsche discretization\n",
    "\n",
    "For the discretization of the interface problem we consider the Nitsche formulation:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "         \\sum_{i \\in \\{+,-\\}} & \\left( \\alpha_i \\nabla u \\nabla v \\right)_{\\Omega_i} + \\left( \\{\\!\\!\\{ - \\alpha \\nabla u \\cdot n \\}\\!\\!\\}, [\\![v]\\!] \\right)_\\Gamma + \\left( [\\![u]\\!],\\{\\!\\!\\{ - \\alpha \\nabla v \\cdot n \\}\\!\\!\\} \\right)_\\Gamma + \\left( \\frac{\\bar{\\alpha} \\lambda}{h}  [\\![u]\\!] , [\\![v]\\!] \\right)_{\\Gamma} \\\\\n",
    "         & = \\left( f,v \\right)_{\\Omega}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "\n",
    "for all $v \\in V_h^\\Gamma$.\n",
    "\n",
    "For this formulation we require:\n",
    "\n",
    "* a suitably defined average operator $\\{ \\cdot \\} = \\kappa_+ (\\cdot)|_{\\Omega_{+}} + \\kappa_- (\\cdot)|_{\\Omega_{-}}$\n",
    "* a suitable definition of the normal direction\n",
    "* numerical integration on $\\Omega_{+}^{lin}$, $\\Omega_{-}^{lin}$ and $\\Gamma^{lin}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Cut ratio field\n",
    "\n",
    "For the average we use the \"Hansbo\"-choice: \n",
    "\n",
    "$$\n",
    "\\kappa_- = \\frac{|T \\cap \\Omega_-^{lin}|}{|T|} \\qquad \n",
    "\\kappa_+ = 1 - \\kappa_- \n",
    "$$\n",
    "\n",
    "This \"cut ratio\" field is provided by the CutInfo class: \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "kappaminus = CutRatioGF(ci)\n",
    "kappa = (kappaminus, 1-kappaminus)\n",
    "Draw(kappaminus, mesh, \"kappa\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Normal direction\n",
    "\n",
    "The normal direction is obtained from the (interpolated) level set function:\n",
    "\n",
    "$$\n",
    "  n^{lin} = \\frac{\\nabla \\phi_h^{lin}}{\\Vert \\nabla \\phi_h^{lin} \\Vert}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "n = Normalize(grad(lsetp1))\n",
    "Draw(n, mesh, \"normal\", vectors={'grid_size': 20})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Averages and jumps\n",
    "\n",
    "Based on the background finite elements we can now define the averages and jumps: \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "h = specialcf.mesh_size\n",
    "\n",
    "alpha = [1.0,20.0]\n",
    "\n",
    "# Nitsche stabilization parameter:\n",
    "stab = 20*(alpha[1]+alpha[0])/h\n",
    "\n",
    "# expressions of test and trial functions (u and v are tuples):\n",
    "u,v = VhG.TnT()\n",
    "\n",
    "average_flux_u = sum([- kappa[i] * alpha[i] * grad(u[i]) * n for i in [0,1]])\n",
    "average_flux_v = sum([- kappa[i] * alpha[i] * grad(v[i]) * n for i in [0,1]])\n",
    "\n",
    "jump_u = u[0] - u[1]\n",
    "jump_v = v[0] - v[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Integrals\n",
    "\n",
    "To integrate only on the subdomains or the interface with a symbolic expression, you have to use the `dCut` differentail symbol, cf. [intlset.ipynb](intlset.ipynb). (Only) to speed up assembly we can mark the integrals as undefined where they would be zero anyway: \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "dx_neg = dCut(levelset=lsetp1, domain_type=NEG, definedonelements=ci.GetElementsOfType(HASNEG))\n",
    "dx_pos = dCut(levelset=lsetp1, domain_type=POS, definedonelements=ci.GetElementsOfType(HASPOS))\n",
    "ds = dCut(levelset=lsetp1, domain_type=IF, definedonelements=ci.GetElementsOfType(IF))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We first integrate over the subdomains: \n",
    "\n",
    "$$\n",
    "\\int_{\\Omega_-} \\alpha_- \\nabla u \\nabla v \\, d\\omega \n",
    "$$\n",
    "\n",
    "$$\n",
    "\\int_{\\Omega_+} \\alpha_+ \\nabla u \\nabla v \\, d\\omega \n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a = BilinearForm(VhG, symmetric = True)\n",
    "a += alpha[0] * grad(u[0]) * grad(v[0]) * dx_neg\n",
    "a += alpha[1] * grad(u[1]) * grad(v[1]) * dx_pos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We then integrate over the interface: \n",
    "\n",
    "$$\n",
    "        \\int_{\\Gamma} \\{\\!\\!\\{ - \\alpha \\nabla u \\cdot \\mathbf{n} \\}\\!\\!\\} [\\![v]\\!] \\, d\\gamma \n",
    "         + \\int_{\\Gamma} \\{\\!\\!\\{ - \\alpha \\nabla v \\cdot \\mathbf{n} \\}\\!\\!\\} [\\![u]\\!] \\, d\\gamma \n",
    "         + \\int_{\\Gamma} \\frac{\\lambda}{h} \\bar{\\alpha} [\\![u]\\!] [\\![v]\\!] \\, d\\gamma\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a += (average_flux_u * jump_v + average_flux_v * jump_u + stab * jump_u * jump_v) * ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Finally, we integrate over the subdomains to get the linear form: \n",
    "\n",
    "$$\n",
    "f(v) = \\int_{\\Omega_-} f_- v \\, d\\omega + \n",
    "       \\int_{\\Omega_+} f_+ v \\, d\\omega\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "coef_f = [1,0]    \n",
    "\n",
    "f = LinearForm(VhG)\n",
    "f += coef_f[0] * v[0] * dx_neg\n",
    "f += coef_f[1] * v[1] * dx_pos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Assembly\n",
    "\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): \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# homogenization of boundary data and solution of linear system\n",
    "def SolveLinearSystem():\n",
    "    gfu.vec[:] = 0\n",
    "    f.vec.data -= a.mat * gfu.vec\n",
    "    gfu.vec.data += a.mat.Inverse(freedofs) * f.vec\n",
    "SolveLinearSystem()\n",
    "\n",
    "DrawDC(lsetp1, gfu.components[0], gfu.components[1], mesh, \"u\", min=0, max=0.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Higher order accuracy\n",
    "In the previous example we used a second order FESpace but only used a second order accurate geometry representation (due to $\\phi_h^{lin}$). \n",
    "\n",
    "We can improve this by applying a mesh transformation technique, cf. [intlset.ipynb](intlset.ipynb):\n",
    "\n",
    "![lsetcurv](graphics/lsetcurv.jpg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# for isoparametric mapping\n",
    "from xfem.lsetcurv import *\n",
    "lsetmeshadap = LevelSetMeshAdaptation(mesh, order=2)\n",
    "deformation = lsetmeshadap.CalcDeformation(levelset)\n",
    "Draw(deformation, mesh, \"deformation\")\n",
    "\n",
    "# alternatively to passing the deformation to dCut me can do the mesh deformation by hand\n",
    "mesh.deformation = deformation\n",
    "a.Assemble()\n",
    "f.Assemble()\n",
    "mesh.deformation = None\n",
    "\n",
    "SolveLinearSystem()\n",
    "\n",
    "\n",
    "DrawDC(lsetp1, gfu.components[0], gfu.components[1], mesh, \"u\", deformation=deformation, min=0, max=0.25)\n",
    "\n",
    "uh = IfPos(lsetp1, gfu.components[1], gfu.components[0])\n",
    "deform_graph = CoefficientFunction((deformation[0], deformation[1], 4*uh))\n",
    "DrawDC(lsetp1, gfu.components[0], gfu.components[1], mesh, \"graph_of_u\", deformation=deform_graph, min=0, max=0.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## XFEM spaces\n",
    "\n",
    "Instead of the CutFEM space \n",
    "$$\n",
    "V_h^\\Gamma = V_h |_{\\Omega_-^{lin}} \\oplus V_h |_{\\Omega_+^{lin}}\n",
    "$$\n",
    "we can use the (same) space with an XFEM characterization:\n",
    "$$\n",
    "V_h^\\Gamma = V_h \\oplus V_h^x\n",
    "$$\n",
    "with the space $V_h^x$ which adds the necessary enrichments. \n",
    "\n",
    "In `ngsxfem` we can also work with this XFEM spaces: \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Vh = H1(mesh, order=2, dirichlet=[1,2,3,4])\n",
    "Vhx = XFESpace(Vh,ci)\n",
    "VhG = FESpace([Vh,Vhx])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "| original | after cut |\n",
    "|-|-|\n",
    "| ![xfem1](graphics/xfem1.png) | ![xfem2](graphics/xfem2.png) |\n",
    "\n",
    "\n",
    "* The space `Vhx` copies all shape functions from `Vh` on cut (`IF`) elements (and stores a sign (`NEG`/`POS`))\n",
    "* The sign determines on which domain the shape function should be supported (and where not)\n",
    "* Advantage: every dof is an active dof (i.e. no dummy dofs)\n",
    "* Need to express $u_+$ and $u_-$ in terms of $u_h^{std}$ and $u_h^x$:\n",
    "\n",
    "  * $u_- = u_h^{std} +$ `neg`($u_h^x$) and $u_+ = u_h^{std} +$ `pos`($u_h^x$)\n",
    "\n",
    "* `neg` and `pos` filter out the right shape functions of `Vhx`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "### express xfem shape functions as cutFEM shape functions:\n",
    "(u_std,u_x), (v_std, v_x) = VhG.TnT()\n",
    "\n",
    "u = [u_std + op(u_x) for op in [neg,pos]]\n",
    "v = [v_std + op(v_x) for op in [neg,pos]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Similar examples and extensions (python demo files)\n",
    "In the source directory (or on http://www.github.com/ngsxfem/ngsxfem ) you can find in the `demos` directory a number of more advanced examples (with fewer explanatory comments). These are:\n",
    "\n",
    "* `unf_interf_prob.py` : Similar to this notebook. The file implements low and high order geometry approximation and gives the choice between a CutFEM and XFEM discretisation. \n",
    "\n",
    "* `fictdom.py` : Fictitious domain/CutFEM diffusion problem (one domain only) with ghost penalty stabilization.\n",
    "\n",
    "* `fictdom_dg.py` : Fictitious domain/CutFEM diffusion problem with a discontiunous Galerkin discretisation and ghost-penalty stabilization.\n",
    "\n",
    "* `ficdom_mlset.py` : Fictitious domain/CutFEM diffusion problem with a geometry described by multiple level set funcions.\n",
    "\n",
    "* `stokesxfem.py` : Stokes interface problem with using XFEM and a Nitsche formulation.\n",
    "\n",
    "* `tracefem.py` or [tracefem_scalar.ipynb](tracefem_scalar.ipynb) : A scalar Laplace-Beltrami problem on a level set surface in 3d using a trace finite element discretization (PDE on the interface only).\n",
    "\n",
    "* `lsetgeoms.py` : Shows a number of pre-implemented 3d level set geometries.\n",
    "\n",
    "* `moving_domain.py` : A scalar convection-diffusion problem posed on a moving domain discretised with an Eulerian time-stepping scheme.\n",
    "\n",
    "* `aggregates/fictdom_aggfem.py` : Fictitious domain/CutFEM diffusion problem (one domain only) with aggregation of FE spaces\n",
    "\n",
    "* `aggregates/fictdom_dg_aggfem.py` : Fictitious domain/CutFEM diffusion with DG discretization and cell aggregation\n",
    "\n",
    "* `spacetime/spacetimeCG_unfitted.py` : A scalar unfitted PDE problem on a moving domain, discretized with CG-in-time space-time finite elements.\n",
    "\n",
    "* `spacetime/spacetimeDG_unfitted.py` : As `spacetimeCG_unfitted.py` but with a DG-in-time space-time discretisation.\n",
    "\n",
    "* `spacetime/spacetimeDG_fitted.py` : A fitted FEM heat equation solved using DG-in-time space-time finite elements.\n",
    "\n",
    "* `spacetime/spacetime_vtk.py` : Demonstration to generate space-time VTK outputs.\n",
    "\n",
    "* `mpi/mpi_nxfem.py` : Same problem as `unf_interf_prob.py` (XFEM + higher order) campatible with MPI parallisation."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.9.1-final"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}