{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 8.7 Unfitted FEM for domains described by multiple level sets\n",
    "\n",
    "In order to solve a partial differential equation on a stationary domain defined by multiple level set functions, we proceed similarly to the demo [cutfem.ipynb](cutfem.ipynb), with a few key differences in order to deal with corners.\n",
    "\n",
    "For simplicity we consider a Poisson problem on the unit square in order to show the differences between the single and multiple level set setting. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The libraries are imported as usual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Basic NGSolve things\n",
    "from netgen.geom2d import SplineGeometry\n",
    "from ngsolve import *\n",
    "from ngsolve.solvers import PreconditionedRichardson as PreRic\n",
    "\n",
    "# ngsxfem and the mlset convenience layer \n",
    "from xfem import *\n",
    "from xfem.mlset import *\n",
    "\n",
    "# Visualisation\n",
    "from ngsolve.webgui import *\n",
    "DrawDC = MakeDiscontinuousDraw(Draw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We want to solve the problem\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    -\\Delta u &= f &&\\text{in }\\Omega\\\\\n",
    "    u &= 0 &&\\text{on } \\partial\\Omega\n",
    "\\end{aligned}$$\n",
    "\n",
    "where the domain $\\Omega\\subset\\widetilde\\Omega$ is a subset of the meshed background domain $\\widetilde\\Omega$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = SplineGeometry()\n",
    "geo.AddRectangle((-0.2, -0.2), (1.2, 1.2), bcs=(\"bottom\", \"right\", \"top\", \"left\"))\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We define the domain as the set $\\Omega := \\bigcap_i \\{\\phi_i < 0\\}$ with $\\{\\phi_i\\} = \\{ -y, x - 1, y - 1, -x \\}$ and use a `DomainTypeArray` as a container for this set and to generate the boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "level_sets = [-y, x - 1, y - 1, -x]\n",
    "level_sets_p1 = tuple(GridFunction(H1(mesh, order=1)) for i in range(len(level_sets)))\n",
    "for i, lsetp1 in enumerate(level_sets_p1):\n",
    "    InterpolateToP1(level_sets[i], lsetp1)\n",
    "    DrawDC(lsetp1, -3.5, 2.5, mesh, \"lsetp1_{}\".format(i))\n",
    "\n",
    "square = DomainTypeArray((NEG, NEG, NEG, NEG))\n",
    "boundary = square.Boundary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Finite Element spaces are defined as usual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 3\n",
    "\n",
    "V = H1(mesh, order=k, dgjumps=True)\n",
    "gfu = GridFunction(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Multiple level set cut-information\n",
    "\n",
    "\n",
    "As before, we use standard FE spaces restricted to the active part of the mesh. In order to handle multiple level sets, we need the `MultiLevelsetCutInfo` class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlci = MultiLevelsetCutInfo(mesh, level_sets_p1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Marking the correct elements\n",
    "\n",
    "Marking the correct active elements needs to be handled differently in the multipe level set setting, since elements of the type `(IF,IF)` are sometimes not part of the active domain, as can be seen in the following image:\n",
    "<center><img src=\"graphics/mlset_cuts.png\" alt=\"\" width=\"500\"/></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In order to select the correct active elements, the `MultiLevelsetCutInfo` class has a member function `GetElementsWithContribution()`. This takes in a `tuple(ENUM)`, a list of such tuples or a `DomainTypeArray`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "els_hasneg = mlci.GetElementsWithContribution(square)\n",
    "els_if = mlci.GetElementsWithContribution(boundary)\n",
    "\n",
    "Draw(BitArrayCF(els_hasneg), mesh, \"els_hasneg\")\n",
    "Draw(BitArrayCF(els_if), mesh, \"els_if\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To deal with arbitraty mesh-interface cuts we use ghost-penalty stabilisation. The elements needed for this can be obtained in the same manner as before and similarly, the free degrees of freedom can also be determined as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "facets_gp = GetFacetsWithNeighborTypes(mesh, a=els_hasneg, b=els_if, \n",
    "                                       use_and=True)\n",
    "\n",
    "freedofs = GetDofsOfElements(V, els_hasneg) & V.FreeDofs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We implement the homogeneous Dirichlet boundary conditions using Nitsche's trick. As each boundaty section has a different unit-normal vector, we have to decompose the Nitsche term into the different segments. Each will be added as a separate weak form, so we also require element markers for each side seperately."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Since the easient and most flexible nomenclature for each segment is the tuple of `DOMAIN_TYPE`s, we store these element markers in a dictionary using the tuple as the key:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "els_if_single = {}\n",
    "for i, dtt in enumerate(boundary):\n",
    "    els_if_single[dtt] = mlci.GetElementsWithContribution(dtt)\n",
    "    Draw(BitArrayCF(els_if_single[dtt]), mesh, \"els_if_singe\"+str(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Defining the weak form\n",
    "\n",
    "As in any other unfitted problem, we define some parameters and get the trial and test functions and mesh-size coefficient functions as usual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma_n = 10\n",
    "gamma_s = 0.5\n",
    "\n",
    "u, v = V.TnT()\n",
    "h = specialcf.mesh_size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For the normal vectors, we use the `DomainTypeArray.GetOuterNormals()` functionality, see [mlset_basic.ipynb](mlset_basic.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "normals = square.GetOuterNormals(level_sets_p1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We then define the cut `SymbolicDifferentialSymbols` for the volume, boundary and ghost-penalty facet-patch integrators, see also [mlset_basic.ipynb](mlset_basic.ipynb). Note that we pass the `BitArrays` on which the integrators are defined to the `SymbolicDiffeentialSymbol`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = dCut(level_sets_p1, square, definedonelements=els_hasneg)\n",
    "ds = {dtt: dCut(level_sets_p1, dtt, definedonelements=els_if_single[dtt])\n",
    "      for dtt in boundary}\n",
    "dw = dFacetPatch(definedonelements=facets_gp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Integrators\n",
    "\n",
    "As we have an unfitted problem, we use a `RestrictedBilinearForm`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = RestrictedBilinearForm(V, element_restriction=els_hasneg, facet_restriction=facets_gp, check_unused=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The diffusion and ghost_penalty forms can then be added normally"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a += InnerProduct(Grad(u), Grad(v)) * dx\n",
    "a += gamma_s / (h**2) * (u - u.Other()) * (v - v.Other()) * dw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For the Nitsche terms, we loop over the normals dictionary. This gives us acess to both the boundary domain-tuples and the normal vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for bnd, n in normals.items():\n",
    "    a += -InnerProduct(Grad(u) * n, v) * ds[bnd]\n",
    "    a += -InnerProduct(Grad(v) * n, u) * ds[bnd]\n",
    "    a += (gamma_n * k * k / h) * InnerProduct(u, v) * ds[bnd]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We want the exact solution to be $u = 16 x (1 - x) y (1 - y)$. To add the appropriate forcing term to the right-hand side, we do not have anything special"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm(V)\n",
    "f += 32 * (y * (1 - y) + x * (1 - x)) * v * dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Solving the system\n",
    "\n",
    "Now the system has been set up, we can assemble the systrem and sove in the usual manner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f.Assemble()\n",
    "a.Assemble()\n",
    "inv = a.mat.Inverse(freedofs=freedofs, inverse=\"sparsecholesky\")\n",
    "\n",
    "gfu.vec.data = PreRic(a=a, rhs=f.vec, pre=inv, freedofs=freedofs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can now visualise the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_ex = 16 * x * (1 - x) * y * (1 - y)\n",
    "grad_u_ex = CoefficientFunction((16 * (1 - 2 * x) * y * (1 - y), \n",
    "                                 16 * x * (1 - x) * (1 - 2 * y)))\n",
    "DrawDC(square.Indicator(level_sets_p1), 0, gfu, mesh, \"solution\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "And compute the error. To integrate with quadrature rules of the correct order we generate new cut `SumbolicDifferentialSymbol`s using the `.order(k)` method with returns a new `DifferentialSymbol` which tells the `Integrate` to integrate with order `k`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "err_l2 = sqrt(Integrate((gfu - u_ex)**2 * dx.order(2 * k), mesh=mesh))\n",
    "err_h1 = sqrt(Integrate((Grad(gfu) - grad_u_ex)**2 * dx.order(2 * (k - 1)), mesh=mesh))\n",
    "print(\"L2 error = {:1.3e}\".format(err_l2))\n",
    "print(\"H1 error = {:1.3e}\".format(err_h1))"
   ]
  }
 ],
 "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"
  },
  "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
}