spacetime_unfitted.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unfitted space-time finite elements\n",
    "\n",
    "In this tutorial, we want illustrate the unfitted space-time functionality of ngsxfem. To this end, we consider a convection-diffusion problem on a moving domain with a manufactured solution (and according right-hand side source term). To illustrate the setting, the simulation end result from ngsxfem is illustrated in the following video:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%html\n",
    "<iframe width=100% height=600px src=\"https://www.youtube.com/embed/16emZemDZak?rel=0\" frameborder=\"0\" allow=\"autoplay; encrypted-media\" allowfullscreen></iframe>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "In this example we consider a moving domain problem with homogeneous Neumann boundary conditions:\n",
    "\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "\\partial_t u + \\mathbf{w} \\cdot \\nabla u - \\alpha \\Delta u &= f \\quad \\text{ in } \\Omega(t),  & \\\\\n",
    "~ \\partial_{\\mathbf{n}} u &=  0  \\quad \\text{ on } \\partial \\Omega(t), & \\\\\n",
    "u &= u_0  \\quad \\text{at } t=0, & \\\\\n",
    "\\end{aligned}\\right.\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\operatorname{div}(\\mathbf{w}) = 0  \\quad \\text{ in } \\Omega(t),  \\quad \\mathbf{w} \\cdot n = \\mathcal{V}_n \\text{ on } \\partial \\Omega(t).\n",
    "$$\n",
    "\n",
    "We consider the Discountinuous Galerkin space-time discretization as discussed in [1, 2], which is of high order in space and time.\n",
    "\n",
    "**Literature:**\n",
    "\n",
    "[1]: F. Heimann, C. Lehrenfeld, J. Preuß. Geometrically Higher Order Unfitted Space-Time Methods for PDEs on Moving Domains. *SIAM Journal on Scientific Computing*, 45(2):B139--B165, 2023. doi: [10.1137/22M1476034](http://dx.doi.org/10.1137/22M1476034)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from netgen.occ import *\n",
    "from xfem import *\n",
    "from math import pi\n",
    "from xfem.lset_spacetime import *\n",
    "ngsglobals.msg_level = 1\n",
    "\n",
    "#import netgen.gui \n",
    "DrawDC = MakeDiscontinuousDraw(Draw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly, we pick some discretisation parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DISCRETIZATION PARAMETERS:\n",
    "\n",
    "# Parameter for refinement study:\n",
    "i = 3\n",
    "n_steps = 2**i\n",
    "space_refs = i\n",
    "\n",
    "# Polynomial order in time\n",
    "k_t = 2\n",
    "# Polynomial order in space\n",
    "k_s = k_t\n",
    "# Polynomial order in time for level set approximation\n",
    "lset_order_time = k_t\n",
    "# Integration order in time\n",
    "time_order = 2 * k_t\n",
    "# Time stepping parameters\n",
    "tstart = 0\n",
    "tend = 0.5\n",
    "delta_t = (tend - tstart) / n_steps\n",
    "maxh = 0.5\n",
    "# Ghost-penalty parameter\n",
    "gamma = 0.05"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Background geometry and mesh:\n",
    "\n",
    "* We consider a simple square as background domain an use a simple mesh for that.\n",
    "* The space-time method uses tensor-product elements. Hence, we do not need space-time meshes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Outer domain:\n",
    "rect = OCCGeometry(Rectangle(1.2,2).Face().Move((-0.6,-1,0)),dim=2)\n",
    "ngmesh = rect.GenerateMesh(maxh=maxh, quad_dominated=False)\n",
    "for j in range(space_refs):\n",
    "    ngmesh.Refine()\n",
    "mesh = Mesh(ngmesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Handling of the time variable\n",
    "\n",
    "For the handling of the space-time integration we use the following rules:\n",
    " * every time step is formulated with respect to the reference interval $[0,1)$ in time\n",
    " \n",
    " * Example: $t_{n-1} = 0.4$, $t=0.55$, $\\Delta t = 0.2$ $\\quad \\longrightarrow \\quad$ $\\hat{t} = 0.75$.\n",
    " \n",
    " * $\\hat{t}$ is the `ReferenceTimeVariable` (`tref`)\n",
    " \n",
    " * We define $t_{old}(=t_{n-1})$ and $\\Delta t$ as a `Parameter`, s.t. we can change the time interval later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Map from reference time to physical time\n",
    "told = Parameter(tstart)\n",
    "t = told + delta_t * tref\n",
    "t.MakeVariable()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Data functions (depending on $t$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Level set geometry\n",
    "# Radius of disk (the geometry)\n",
    "R = 0.5\n",
    "# Position shift of the geometry in time\n",
    "rho = (1 / (pi)) * sin(2 * pi * t)\n",
    "# Convection velocity:\n",
    "w = CoefficientFunction((0, rho.Diff(t)))\n",
    "# Level set\n",
    "r = sqrt(x**2 + (y - rho)**2)\n",
    "levelset = r - R\n",
    "\n",
    "# Diffusion coefficient\n",
    "alpha = 1\n",
    "# Solution\n",
    "u_exact = cos(pi * r / R) * sin(pi * t)\n",
    "# R.h.s.\n",
    "coeff_f = (u_exact.Diff(t)\n",
    "           - alpha * (u_exact.Diff(x).Diff(x) + u_exact.Diff(y).Diff(y))\n",
    "           + w[0] * u_exact.Diff(x) + w[1] * u_exact.Diff(y)).Compile()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### A View on the time-dependent level set function on $[0,$`tend`$]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_one_timestep = sqrt(x**2 + (y - (1 / (pi)) * sin(2 * pi * tref * tend))**2)\n",
    "levelset_one_timestep = r_one_timestep - R\n",
    "#TimeSlider_Draw(levelset_one_timestep,mesh,autoscale=False,min=-0.02,max=0.02,deformation=True)\n",
    "\n",
    "from helper import ProjectOnMultiDimGF\n",
    "lset_to_show = ProjectOnMultiDimGF(levelset_one_timestep,mesh,order=3,sampling=8)\n",
    "Draw(lset_to_show,mesh,\"u_exact\",autoscale=False,min=-0.5,max=1, interpolate_multidim=True, animate=True, deformation=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "u_exact_one_timestep = cos(pi * r_one_timestep / R) * sin(pi * tref * tend)\n",
    "\n",
    "from helper import ProjectPairOnMultiDimGF\n",
    "lset_cf_to_show = ProjectPairOnMultiDimGF(levelset_one_timestep, u_exact_one_timestep,mesh,order=3,sampling=8)\n",
    "\n",
    "Draw(lset_cf_to_show,mesh,\"u_exact\",eval_function=\"value.x>0.0?value.z:value.y\",autoscale=False,\n",
    "     min=-0.5,max=1, interpolate_multidim=True, animate=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Space-Time finite elements\n",
    "\n",
    "* For the construction of a space-time `FESpace` we can combine any spatial `FESpace` with a scalar `FiniteElement` in a tensor-product fashion.\n",
    "* Here, we use a nodal `FiniteElement` to simplify the extraction of spatial functions at fixed times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Spatial FESpace for solution\n",
    "fes1 = H1(mesh, order=k_s, dgjumps=True)\n",
    "# Time finite element (nodal!)\n",
    "tfe = ScalarTimeFE(k_t)\n",
    "# (Tensor product) space-time finite element space\n",
    "st_fes = tfe * fes1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Levelset mesh adaptation class\n",
    "\n",
    "In order to achieve higher order accuracy in space and time, a variant of the isoparametric mapping for stationary domains is applied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Space time version of Levelset Mesh Adapation object. Also offers integrator\n",
    "# helper functions that involve the correct mesh deformation\n",
    "lsetadap = LevelSetMeshAdaptation_Spacetime(mesh, order_space=k_s,\n",
    "                                            order_time=lset_order_time,\n",
    "                                            threshold=0.5,\n",
    "                                            discontinuous_qn=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Space-Time version of the `CutInfo` class\n",
    "The `CutInfo` class also works for space-time geometries. Its initialization is trivial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "ci = CutInfo(mesh, time_order=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the argument `time_order=0` which makes the `CutInfo` looking for space-time cut information afterwards.\n",
    "\n",
    "In addition, we define a Bitarray for the facets of the mesh for later use in the definition of the Ghost-penalty stabilisation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ba_facets = BitArray(mesh.nfacet)\n",
    "active_dofs = BitArray(st_fes.ndof)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "To Update the slab geometry later on (for updated values of `told`) we do the following:\n",
    " * update of the isoparametric mapping\n",
    " * update of the cut information\n",
    " * update of facet markers\n",
    " * update of dof markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def UpdateTimeSlabGeometry():\n",
    "    lsetadap.CalcDeformation(levelset)\n",
    "\n",
    "    # Update markers in (space-time) mesh\n",
    "    ci.Update(lsetadap.levelsetp1[INTERVAL], time_order=0)\n",
    "\n",
    "    # re-compute the facets for stabilization:\n",
    "    ba_facets[:] = GetFacetsWithNeighborTypes(mesh,\n",
    "                                              a=ci.GetElementsOfType(HASNEG),\n",
    "                                              b=ci.GetElementsOfType(IF))\n",
    "    active_dofs[:] = GetDofsOfElements(st_fes, ci.GetElementsOfType(HASNEG))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that here the call of CalcDeformation of lsetadap entails the calculation of the P1 projection of the levelset function internally. The space-time P1-in-space level set approximation of `lsetadap` can be accessed by `lsetadap.levelsetp1[timetype]` where `timetype` is either `INTERVAL` which yields the space-time function or `TOP` or `BOTTOM` which yields the spatial P1 function that is obtained by restriction to either `tref=1` or `tref=0`. Similarly the access to the deformation is organized as `lsetadap.deformation[timetype]`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solution GridFunction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(st_fes)\n",
    "u_last = CreateTimeRestrictedGF(gfu, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Variational formulation\n",
    "\n",
    "Now we would like to derive a suitable variational formulation on the time slabs $Q^{n}$. \n",
    "\n",
    "We start by multiplying the equation  \n",
    "\\begin{equation*}\n",
    "\\partial_{t} u- \\alpha \\Delta{u} + w \\cdot \\nabla{u} = f \\quad  in \\quad \\Omega(t),   \\qquad  t \\in [0,T] \n",
    "\\end{equation*}\n",
    "by a test function $v$ and perform integration by parts. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Due to homogeneous Neumann boundary conditions this leads to: \n",
    "\\begin{equation*}\n",
    "(\\partial_{t} u, v)_{Q^n} + \\alpha (\\nabla{u},\\nabla{v})_{Q^n}   + (w \\cdot \\nabla{u},v)_{Q^n} = (f,v)_{Q^n}.\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Upwind DG in time\n",
    "In order to impose weak continuity in time, an upwind stabilisation is added, yielding\n",
    "\n",
    "\\begin{equation*}\n",
    "(\\partial_{t} u, v)_{Q^n} + \\alpha (\\nabla{u},\\nabla{v})_{Q^n}  + (w \\cdot \\nabla{u},v)_{Q^n} + (u_{n−1}^+,v_{n−1}^+)_{\\Omega^{n−1}} = (f,v)_{Q^n} + (u_{n−1}^-,v_{n−1}^+)_{\\Omega^{n−1}}.\n",
    "\\end{equation*}\n",
    "\n",
    " ![limits-time-slab](graphics/limits-time-slab.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Ghost penalty stabilization\n",
    "To gain sufficient control on all space-time d.o.f.s we add a so-called *Ghost-Penalty* stabilization \n",
    "as in [1]. Adding the stabilization, the variational formulation on the time slabs becomes:\n",
    " \n",
    "\\begin{aligned}\n",
    " &(\\partial_t u, v)_{Q^{n}} + \\alpha (\\nabla{u},\\nabla{v})_{Q^{n}} + (v, \\nabla{u} \\cdot w)_{Q^{n}}  + (u_{n−1}^+,v_{n−1}^+)_{\\Omega^{n−1}} + s_h(u,v) \\\\\n",
    " &= (f,v)_{Q^{n}}  +  (u_{n−1}^-,v_{n−1}^+)_{\\Omega^{n−1}}          \\\\\n",
    "\\end{aligned}  \n",
    "\n",
    "$$\n",
    "\\text{with} \\qquad\\qquad\n",
    "s_h(u,v) =   \\sum\\limits_{F \\in F_{h}}{ \\gamma_{j} \\int\\limits_{t_{n-1}}^{t_{n}}{   \\int\\limits_{\\omega_F}{  h^{-2} [\\![ u]\\!] \\, [\\![ v]\\!]         \\, d\\mathbf{x} \\, dt.  } }\t\t}                 \\\\\n",
    "$$\n",
    "where $[\\![u]\\!]$ is the difference of $u|_{T_1}$ and $u|_{T_2}$ (interpreted as polynomials $\\in \\mathbb{R}^d$) and $F_h$ is the set of facets on which the stabilization shall be applied.\n",
    "![macro-element](graphics/macro-element.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Implementation of space-time integrals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "u,v = st_fes.TnT()\n",
    "h = specialcf.mesh_size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Transformation from reference interval to $(t_{n-1},t_n)$:\n",
    "The transformation\n",
    "$$\n",
    "(x,\\hat{t}) \\to (x,t_{n-1} + \\hat{t} \\Delta t), \\qquad v(x,t) = \\hat{v}(x,\\hat{t}), \\quad u(x,t) = \\hat{u}(x,\\hat{t}).\n",
    "$$\n",
    "implies the following for the time derivative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dt(u):\n",
    "    return 1.0 / delta_t * dtref(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we define integration region symbols, which are the numerical counterparts of the regions introduced above. Note that the levelset deformation is included to achieve higher order in space. The `definedonelements` information is not necessary (for the first three symbols), but helpful for efficiency reasons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dQ = delta_t * dCut(lsetadap.levelsetp1[INTERVAL], NEG, time_order=time_order,\n",
    "                    deformation=lsetadap.deformation[INTERVAL],\n",
    "                    definedonelements=ci.GetElementsOfType(HASNEG))\n",
    "dOmold = dCut(lsetadap.levelsetp1[BOTTOM], NEG,\n",
    "              deformation=lsetadap.deformation[BOTTOM],\n",
    "              definedonelements=ci.GetElementsOfType(HASNEG), tref=0)\n",
    "dOmnew = dCut(lsetadap.levelsetp1[TOP], NEG,\n",
    "              deformation=lsetadap.deformation[TOP],\n",
    "              definedonelements=ci.GetElementsOfType(HASNEG), tref=1)\n",
    "dw = delta_t * dFacetPatch(definedonelements=ba_facets, time_order=time_order,\n",
    "                           deformation=lsetadap.deformation[INTERVAL])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we setup the bilinear form and linear form corresponding to the previously described discrete variational formulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a = RestrictedBilinearForm(st_fes, \"a\", check_unused=False,\n",
    "                           element_restriction=ci.GetElementsOfType(HASNEG),\n",
    "                           facet_restriction=ba_facets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "First integral:  $(\\partial_t u, v)_{Q^{n}}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a += v * (dt(u) - dt(lsetadap.deform) * grad(u)) * dQ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that here, due to the time-dependent mesh deformation, the partial derivative in physical coordinates that we used in the formulation before corresponds to the partial derivative w.r.t. the reference configuration minus an additional mesh velocity contribution. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Second integral: $\\alpha (\\nabla{u},\\nabla{v})_{Q^{n}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a += alpha * InnerProduct(grad(u), grad(v)) * dQ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Third integral: $(v, \\nabla{u} \\cdot w)_{Q^{n}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a += v * w * grad(u) * dQ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Fourth integral: $(u_{n−1}^+,v_{n−1}^+)_{\\Omega^{n−1}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a += u * v * dOmold"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Fifth integral:\n",
    "$$ s_h(u,v) =   \\sum\\limits_{F \\in F_{h}}{ \\gamma_{j} \\int\\limits_{t_{n-1}}^{t_{n}}{   \\int\\limits_{\\omega_F}{  h^{-2} [\\![ u]\\!] \\, [\\![ v]\\!]         \\, d\\mathbf{x} \\, dt.  }}} =   \\sum\\limits_{F \\in F_{h}}{ \\Delta t \\ \\gamma_{j} \\int\\limits_{t_{n-1}}^{t_{n}}{   \\int\\limits_{\\omega_F}{  h^{-2} [\\![ \\hat u]\\!] \\, [\\![ \\hat v]\\!]         \\, d\\mathbf{x} \\, dt.  }}}  $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a += h**(-2) * (1 + delta_t / h) * gamma * \\\n",
    "    (u - u.Other()) * (v - v.Other()) * dw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Sixth integral: $(f,v)_{Q^{n}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = LinearForm(st_fes)\n",
    "f += coeff_f * v * dQ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Seventh integral: $(u_{n−1}^-,v_{n−1}^+)_{\\Omega^{n−1}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f += u_last * v * dOmold"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solution of linear systems in every time step\n",
    "* setup the new linear system\n",
    "* solve the linear system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def SolveForTimeSlab():\n",
    "    a.Assemble(reallocate=True)\n",
    "    f.Assemble()\n",
    "    inv = a.mat.Inverse(active_dofs)\n",
    "    gfu.vec.data = inv * f.vec.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### At the end of every time step, we\n",
    "* store the solution at $t_n$ into a (purely) spatial `GridFunction` (to be used in next time step)\n",
    "* compute the error\n",
    "* update visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def FinalizeStep(scene=None):\n",
    "    # Evaluate upper trace of solution for\n",
    "    #  * for error evaluation\n",
    "    #  * upwind-coupling to next time slab\n",
    "    RestrictGFInTime(spacetime_gf=gfu, reference_time=1.0, space_gf=u_last)\n",
    "\n",
    "    # Compute error at final time\n",
    "    l2error = sqrt(Integrate((u_exact - u_last)**2 * dOmnew, mesh))\n",
    "    print(\"\\rt = {0:12.9f}, L2 error = {1:12.9e}\".format(told.Get(), l2error), end=\"\")\n",
    "    if scene:\n",
    "        scene.Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The final time loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scene = DrawDC(lsetadap.levelsetp1[TOP],u_last,-2, mesh,\"u\", autoscale=False, min = -2, max = 1, deformation = lsetadap.deformation[TOP])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to the mesh deformation changing (discontinuously) between time slabs, we have to project solutions from one time slab to the other. This is automatically done for all `GridFunctions` registered in `lsetadap` by `ProjectOnUpdate`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Set initial values\n",
    "u_last.Set(fix_tref(u_exact, 0))\n",
    "# Project u_last at the beginning of each time step\n",
    "lsetadap.ProjectOnUpdate(u_last)\n",
    "told.Set(0)\n",
    "gf_lset_to_show = GridFunction(H1(mesh, order=1),multidim=0)\n",
    "gf_u_to_show = GridFunction(u_last.space,multidim=0)\n",
    "timestep = 0\n",
    "md_comps = 0\n",
    "with TaskManager():\n",
    "    while tend - told.Get() > delta_t/2:\n",
    "        UpdateTimeSlabGeometry()\n",
    "        if timestep % 1 == 0:\n",
    "            gf_lset_to_show.AddMultiDimComponent(lsetadap.levelsetp1[BOTTOM].vec)\n",
    "            gf_u_to_show.AddMultiDimComponent(u_last.vec)\n",
    "            md_comps += 1\n",
    "        SolveForTimeSlab()       \n",
    "        timestep += 1\n",
    "        told.Set(told.Get() + delta_t)\n",
    "        FinalizeStep(scene)\n",
    "gf_lset_to_show.AddMultiDimComponent(lsetadap.levelsetp1[TOP].vec)\n",
    "gf_u_to_show.AddMultiDimComponent(u_last.vec)\n",
    "print(\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helper import MultiDimPairToGF\n",
    "gf_pair_to_show = MultiDimPairToGF(gf_lset_to_show,gf_u_to_show,mesh,sampling=md_comps)\n",
    "Draw(gf_pair_to_show, mesh,\"uh\",eval_function=\"value.x>0.0?value.z:value.y\",autoscale=False,\n",
    "     min=-0.5,max=1, interpolate_multidim=True, animate=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Play around suggestions:\n",
    "\n",
    "* use other orders in space or time (and a coarse grid)\n",
    "* use different level set evolutions\n",
    "\n",
    "Further material:\n",
    "\n",
    "Take a look at the demos in `demos/spacetime` of the `ngsxfem` repository."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.11.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}