{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 8.6 Integration over domains described by multiple level sets\n",
    "\n",
    "The `ngsxfem` add-on library has the ability to integrate over domains which are defined by muliple level set functions.\n",
    "\n",
    "We recall that previously, our domain of interest was described by a single level set function, i.e., for a continuous level set function $\\phi$, we considered the domain $\\Omega_{-} := \\{x : \\phi(x) < 0 \\}$.\n",
    "\n",
    "Now, let us consider the set $\\{\\phi_i \\}$ of level set functions. With this we define the domain of interest by the intersection $\\Omega := \\bigcap_i \\{\\phi_i < 0\\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The libraries are imported as usual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Import geometry features, NGSolve and xfem\n",
    "from netgen.geom2d import SplineGeometry\n",
    "from ngsolve import *\n",
    "from xfem import *\n",
    "\n",
    "# Visualisation\n",
    "from ngsolve.webgui import *\n",
    "DrawDC = MakeDiscontinuousDraw(Draw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For a better interactive experience set the following parameter to 2 (1 leads to coarse grids). It is set to 1 for the documentation to keep file sizes at a managable level:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interactive = 1"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generate a background mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = SplineGeometry()\n",
    "geo.AddRectangle([-1, -0.5], [1, 1.5], bc=1)\n",
    "mesh = Mesh (geo.GenerateMesh(maxh=0.3, quad_dominated=False))\n",
    "Draw(mesh, 0, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Integrating over regions of co-dimension 0\n",
    "\n",
    "As an example, we shall consider an isosceles triangle, described by the level set functions $\\phi_1 = y - 1$, $\\phi_2 = 2x - y$ and $\\phi_3 = -2x - y$. \n",
    "\n",
    "As in the single level set setting, `ngsxfem` can only handle piecewise linear level set functions. We thefore interpolate each smooth level set into an P1 GridFunction $\\phi^\\ast_i \\in \\mathcal{P}^1(T), T \\subset \\Omega$. In order to use these more easily later on, we store them in a `tuple.`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "level_sets = (y-1, 2*x-y, -2*x-y)\n",
    "nr_ls = len(level_sets)\n",
    "level_sets_p1 = tuple(GridFunction(H1(mesh, order=1)) for i in range(nr_ls))\n",
    "\n",
    "for i, lset_p1 in enumerate(level_sets_p1):\n",
    "    InterpolateToP1(level_sets[i], lset_p1)\n",
    "    DrawDC(lset_p1, -3.5, 2.5, mesh, \"lset_p1_{}\".format(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The domain described where all three level-sets are negavive is then"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw(CoefficientFunction((level_sets_p1[0], level_sets_p1[1], level_sets_p1[2], 0)),\n",
    "     mesh, \"NEG-NEG-NEG\", min=-1, max=1,\n",
    "     eval_function=\"value.x<0.0?(value.y<0.0?(value.z<0.0?1.0:-1.0):-1.0):-1.0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To integrate, we still use the function `Integrate`. In the multiple level set context, the input to the `dCut(levelset, domain_type, order)` is very similar to the single level set setting. Here, \n",
    "* the `levelset` argument  must be a `tuple` of the P1-GridFunction level set functions and \n",
    "* the `domain_type` can be a \n",
    "  * `tuple({NEG, POS, IF})` of the same length as the levelset tuple which describes the domain whith respect each level set function\n",
    "  * `list(tuple)` of such tuples when the domain is described via (the union of) multiple regions\n",
    "  * `DomainTypeArray` (see `xfem.mlset` convininece layer below).\n",
    "* the remaining (keyword) arguments remain as in the single level set case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "area = Integrate(CoefficientFunction(1) * dCut(level_sets_p1, (NEG, NEG, NEG), order=0), mesh=mesh)\n",
    "error = abs(area - 0.5)\n",
    "print(\"Result of the integration: {}\".format(area))\n",
    "print(\"Error of the integration: {:5.3e}\".format(error))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Since the smooth level set functions are linear, the integration is exact up to mashine precision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Integrating over regions of  co-dimension > 0\n",
    "\n",
    "We can also integrate over subdomains of higher codimensions. For example, the lenghths of the sides of the triangle can be computed as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len_top = Integrate(CoefficientFunction(1) * dCut(level_sets_p1, (IF, NEG, NEG), order=0), mesh)\n",
    "print(\"Top side-length of the triangle\");\n",
    "print(\"Result of the integration: {}\".format(len_top))\n",
    "print(\"Error of the integration: {:5.3e}\".format(abs(len_top - 1)))\n",
    "\n",
    "len_right = Integrate(CoefficientFunction(1) * dCut(level_sets_p1, (NEG, IF, NEG), order=0), mesh)\n",
    "print(\"\\nRight side-length of the triangle\")\n",
    "print(\"Result of the integration: {}\".format(len_right))\n",
    "print(\"Error of the integration: {:5.3e}\".format(abs(len_right - sqrt(5)/2)))\n",
    "\n",
    "len_left = Integrate(CoefficientFunction(1) * dCut(level_sets_p1, (NEG, NEG, IF), order=0), mesh)\n",
    "print(\"\\nLeft side-length of the triangle\")\n",
    "print(\"Result of the integration: {}\".format(len_left))\n",
    "print(\"Error of the integration: {:5.3e}\".format(abs(len_left - sqrt(5)/2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To compute the perimeter of the triangle, we can pass a list containing the tuples describing each side to `Integrate` via the `domain_type` entry in the `levelset_domain` dictionary: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "domain_list = [(IF, NEG, NEG), (NEG, IF, NEG), (NEG, NEG, IF)]\n",
    "len_together =  Integrate(CoefficientFunction(1) * dCut(level_sets_p1, domain_list, order=0), mesh)\n",
    "\n",
    "\n",
    "print(\"Perimeter of the triangle\")\n",
    "print(\"Result of the integration: {}\".format(len_together))\n",
    "print(\"Error of the integration: {:5.3e}\".format(abs(len_together - 1 - sqrt(5))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Similarly, we can also perform point evaluations at (codim 2) intersections of level sets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "point_val = Integrate((x**2 + y**2) * dCut(level_sets_p1, (IF, IF, NEG), order=0), mesh)\n",
    "\n",
    "print(\"Result of the integration: {}\".format(point_val))\n",
    "print(\"Error of the integration: {:5.3e}\".format(abs(point_val - 1.25)))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convenience Layer: `xfem.mlset`\n",
    "\n",
    "In order to work with more complex domains, described by multiple level sets, we provide a convenience layer in the module `xfem.mlset`. At its core, this module provides a container class `DomainTypeArray` which allows us to perform operations and manipulate domain regions. A full description of this can be cound in the docstring:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xfem.mlset import *\n",
    "\n",
    "DomainTypeArray?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Working with `DomainTypeArray`s\n",
    "A `DomainTypeArray` is initialised either with a `tuple({NEG,POS,IF,ANY})` or with a list of such tuples. While integration does not work with `COMBINED_DOMAIN_TYPES`, the `DomainTypeArray` class accepts `ANY` as an entry, and expands this internaly to a list of tuples containing only `{IF, POS, NEG}`, such that the result can then be used for integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "domain = DomainTypeArray((ANY, NEG, POS))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The list of regions contained in this can then be acessed via the `as_list` attribute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "domain.as_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "and the codimension of the resulting domain can be accesed via the `codim` attribute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "domain.codim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Geometrical features\n",
    "\n",
    "The main feature of the `DomainTypeArray` class is that we can treat it as a geometrical set and thus can perform logical operations with them. Available operators are:\n",
    "\n",
    "Operator | Operation | Remark\n",
    ":-- |:-- | :--\n",
    "`~` | Returns a `DomainTypeArray` describing the region **NOT** part of the original region. | The result keeps the codimension.\n",
    "` \\| ` | Returns a `DomainTypeArray` describing the region defined by the **UNION** of two regions. | Only possible for domains of the same codimension.\n",
    "`&` | Returns a `DomainTypeArray` describing the region defined by the **INTERSECTION** of two regions. | Intersections can increase the codimension.\n",
    "\n",
    "\n",
    "Note: \n",
    " - For `|` and `&`, the tuples in each `DomainTypeArray` must have the same length.\n",
    " - The operators `&=` and `|=` are also available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "triangle = DomainTypeArray((NEG,NEG,NEG))\n",
    "outside_triangle = ~triangle\n",
    "outside_triangle.as_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dta1 = DomainTypeArray((POS, NEG)) | DomainTypeArray([(POS,POS),(NEG,POS)])\n",
    "dta1.as_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "dta2 = outside_triangle & DomainTypeArray([(POS, NEG, POS), (POS, POS, POS)])\n",
    "dta2.as_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dta3 = DomainTypeArray((NEG, IF, POS)) & DomainTypeArray((NEG, POS, IF))\n",
    "dta3.as_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Generating the boundary**\n",
    "\n",
    "In order to obtain the boundary of the region described by a `DomainTypeArray` we can use the class `.Boundary()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary = triangle.Boundary()\n",
    "boundary.as_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Since each section in a `DomainTypeArray` is considered as a separate domain, we can even go one step further and construct the corners:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corners = boundary.Boundary()\n",
    "corners.as_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Visualisation Tools\n",
    "\n",
    "To help visualise the regions descibed by a DomainTypeArray for debugging purposes, we can use the `.Indicator(lsets)` method for domains of codim=0 and the `IndicatorSmoothed(lsets, eps=0.03)` method for domains of codim>0. \n",
    "\n",
    "These functions return a `CoefficientFunction` which has the value 1 inside the region of interest and is 0 elsewhere. For regions of co-dimension>0, the result of `IndicatorSmoothed(lsets, eps)` is a `CoefficientFunction` with the value 1 in the `eps` region around the domain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For illustration purposes, we shall consider fine meshes here, in order to render the indicator functions accurately within the NGSolve webgui. Alternaltievly, it is possible to visualise the indicator functions on coarse meshes in the netgen gui or in paraview (using vtk outputs). In these cases, a larger number of subdivisions is then necessary to get a resonalbly sharp representation on elements with multiple cuts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if interactive > 1:\n",
    "    nref = 4\n",
    "else:\n",
    "    nref = 2\n",
    "for i in range(nref):\n",
    "    mesh.ngmesh.Refine()\n",
    "mesh=Mesh(mesh.ngmesh)    \n",
    "\n",
    "for i, lset_p1 in enumerate(level_sets_p1):\n",
    "    lset_p1.Update()\n",
    "    InterpolateToP1(level_sets[i], lset_p1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DrawDC(outside_triangle.Indicator(level_sets_p1), -3.5, 2.5, mesh, \"outside_indicator\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Similarly, `DomainTypeArray.IndicatorSmoothed(level_sets, eps)` returns a CoefficientFunction which has the value 1 in the `eps`-strip around the interface described in `DomainTypeArray`. This assumes that each level set function is approximately a distance function around the interface:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DrawDC(boundary.IndicatorSmoothed(level_sets_p1, 0.05/interactive), -3.5, 2.5, mesh, \n",
    "       \"boundary_indicator\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DrawDC(corners.IndicatorSmoothed(level_sets_p1, 0.06/interactive), -3.5, 2.5, mesh, \"boundary_indicator\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Removing unnecesary regions\n",
    "\n",
    "These operations take place on a pure logical `DOMAIN_TYPE` level, without taking ang geometrical information into account. \n",
    "\n",
    "Let us again consider the above case of a Zalesak Disk consisting of the 5 regions, marked with a $+$:\n",
    "\"<center><img src=\"./graphics/z_disc.png\" alt=\"\" width=\"500\"/></center>\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "geo = SplineGeometry()\n",
    "geo.AddRectangle((-1.1, -1.1), (1.1, 1.1), bc=1)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.2/interactive))\n",
    "\n",
    "level_sets = [x * x + y * y - 1, -x - 1 / 3, x - 1 / 3, y - 0.5]\n",
    "level_sets_p1 = tuple(GridFunction(H1(mesh, order=1)) \n",
    "                      for i in range(len(level_sets)))\n",
    "\n",
    "for i, lsetp1 in enumerate(level_sets_p1):\n",
    "    InterpolateToP1(level_sets[i], lsetp1)\n",
    "    DrawDC(lsetp1, -3.5, 2.5, mesh, \"lset_p1_\" + str(i))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The easiest way to generate this is to use the `DomainTypeArray` functionality:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_disc =    DomainTypeArray((NEG,ANY,ANY,ANY)) \\\n",
    "         & ~DomainTypeArray((NEG,NEG,NEG,NEG))\n",
    "z_disc.as_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DrawDC(z_disc.Indicator(level_sets_p1), -3.5, 2.5, mesh, \"zdisc\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "However, we note that there are 7 regions, rather than the expected 5. However, on closer inspection, we see that `(NEG, POS, POS, POS)` and `(NEG, POS, POS, NEG)` are regions that do not exist with respect to our level sets."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We have two options to remove such irrelevant regions from a `DomainTypeArray`:\n",
    "\n",
    "First: \n",
    "\n",
    "If the object has already been created, we can use the `.Compress(lsets, persistent=False)` method. Given a tuple of level set functions `lsets`, the `.Compress` method removes any domain tuples from the instance, which have zero measure with respect to these level set functions. If `persistent=True` is passed to `.Compress`, then any `DomainTypeArray`s derived from the compressed one will also be compressed with respect to the here given set of level set fucntions. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "So that we are able to integrate internally, `lsets` needs to be of the type `tuple(ngsolve.GridFunction)`, such that we have acess to the mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_disc.as_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "z_disc.Compress(level_sets_p1)\n",
    "z_disc.as_list"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This removes the unnecessary regions `(NEG, POS, POS, POS)` and `(NEG, POS, POS, NEG)`."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Second:\n",
    "\n",
    "On initialisation of a `DomainTypeArray`, we can also pass a tuple `lsets` of level sets for compression at initilisation, and a flag `persistent_compress` which indicates, whether the result when using operators schould again be compressed with respect to the same set of level sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp1 = DomainTypeArray((NEG,ANY,ANY,ANY), \n",
    "                       lsets=level_sets_p1, \n",
    "                       persistent_compress=True)\n",
    "tmp2 = DomainTypeArray((NEG,NEG,NEG,NEG), \n",
    "                       lsets=level_sets_p1, \n",
    "                       persistent_compress=True)\n",
    "z_disc2 = tmp1 & ~tmp2\n",
    "z_disc2.as_list"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The difference between `persistent_compress=True/False` becomes more apparent when computing the boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_disc = DomainTypeArray((NEG,ANY,ANY,ANY)) & ~DomainTypeArray((NEG,NEG,NEG,NEG))\n",
    "bnd1 = z_disc.Boundary()\n",
    "DrawDC(bnd1.IndicatorSmoothed(level_sets_p1,0.08/interactive), -3.5, 2.5, mesh, \"bnd1\")\n",
    "bnd1.as_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "bnd2 = z_disc2.Boundary()\n",
    "DrawDC(bnd2.IndicatorSmoothed(level_sets_p1,0.08/interactive), -3.5, 2.5, mesh, \"bnd2\")\n",
    "bnd2.as_list"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Generating outward pointing unit normal vectors\n",
    "\n",
    "To implement boundary conditions using Nitsche's method, for example in a CutFEM problem, the outward pointing unit normal vectors are needed. Since in the multiple level set setting the boundary can consist of many sections (8 in the case of the Zalesak disk), we have a member function to compute the outward pointing unit normal vector on each section of the boundary of a codim=0 domain. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "normals = z_disc2.GetOuterNormals(level_sets_p1)\n",
    "normals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ngsolve.internal.visoptions.showsurfacesolution=1\n",
    "Draw(normals[(IF,NEG,NEG,POS)], mesh, \"normal_1\", vectors={'grid_size': 14})\n",
    "DrawDC(DomainTypeArray((IF,NEG,NEG,POS)).IndicatorSmoothed(level_sets_p1, eps=0.08/interactive), -3.5, 2.5, mesh, \"boundary_1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Draw(normals[(NEG,NEG,NEG,IF)], mesh, \"normal_2\", vectors={'grid_size': 14})\n",
    "DrawDC(DomainTypeArray((NEG,NEG,NEG,IF)).IndicatorSmoothed(level_sets_p1, eps=0.08/interactive), -3.5, 2.5,  mesh, \"boundary_2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Draw(normals[(IF,NEG,POS,NEG)], mesh, \"normal_3\", vectors={'grid_size': 14})\n",
    "DrawDC(DomainTypeArray((IF,NEG,POS,NEG)).IndicatorSmoothed(level_sets_p1, eps=0.08/interactive), -3.5, 2.5,  mesh, \"boundary_3\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Draw(normals[(NEG,NEG,IF,NEG)], mesh, \"normal_4\", vectors={'grid_size': 14})\n",
    "DrawDC(DomainTypeArray((NEG,NEG,IF,NEG)).IndicatorSmoothed(level_sets_p1, eps=0.08/interactive), -3.5, 2.5,  mesh, \"boundary_4\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Draw(normals[(IF,POS,NEG,POS)], mesh, \"normal_5\", vectors={'grid_size': 14})\n",
    "DrawDC(DomainTypeArray((IF,POS,NEG,POS)).IndicatorSmoothed(level_sets_p1, eps=0.08/interactive), -3.5, 2.5,  mesh, \"boundary_5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Draw(normals[(IF,POS,NEG,NEG)], mesh, \"normal_6\", vectors={'grid_size': 14})\n",
    "DrawDC(DomainTypeArray((IF,POS,NEG,NEG)).IndicatorSmoothed(level_sets_p1, eps=0.08/interactive), -3.5, 2.5,  mesh, \"boundary_6\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Draw(normals[(NEG,IF,NEG,NEG)], mesh, \"normal_7\", vectors={'grid_size': 14})\n",
    "DrawDC(DomainTypeArray((NEG,IF,NEG,NEG)).IndicatorSmoothed(level_sets_p1, eps=0.08/interactive), -3.5, 2.5,  mesh, \"boundary_7\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Draw(normals[(IF,NEG,POS,POS)], mesh, \"normal_8\", vectors={'grid_size': 14})\n",
    "DrawDC(DomainTypeArray((IF,NEG,POS,POS)).IndicatorSmoothed(level_sets_p1, eps=0.08/interactive), -3.5, 2.5,  mesh, \"boundary_8\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Integrating using `DomainTypeArray`s\n",
    "\n",
    "In order to use DomainTypeArrays for integration purposes we can simply pass the `DomainTypeArray` to the `domain_type` argumet of the SymbolicDifferentialSymbol `dCut`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from math import pi\n",
    "area_zd = pi - 2 / 3 * (0.5 + sqrt(2) * 2 / 3)\n",
    "area_zd += - (2 * asin(1 / 3) - sin(2 * asin(1 / 3))) / 2\n",
    "\n",
    "area = Integrate(CoefficientFunction(1) * dCut(level_sets_p1, z_disc2, order=0), mesh)\n",
    "\n",
    "print(\"Area 2 = {:12.10f}\".format(area))\n",
    "print(\"Error: {:4.2e}\".format(abs(area - area_zd)))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Since the smooth level sets are not all piecewise linear, integration is only accurate \n",
    "up to $\\mathcal{O}(h^2)$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Constructing new domais from `DomainTypeArrays` defined with respect to different level sets**\n",
    "\n",
    "Some complicated gemetries consist of the union or the intersection of geometrical objects which are easily descriped be a small number of level sets. For this purpose `xfem.mlset` provides the two functions `TensorUnion` and `TensorIntersection`, which generate the `DomainTypeArray` description from the union or intersection of different `DomainTypeArray`s, each defined with respect to different level sets.\n",
    "\n",
    "For example, consider the follwoing two triangular domains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = SplineGeometry()\n",
    "geo.AddRectangle((-2.5, -1), (1, 2.5), bc=1)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.3/interactive))\n",
    "\n",
    "level_sets1 = [-2*x + y - 1, 2*x + y - 1, -y + 2]\n",
    "dta1 = DomainTypeArray((POS, POS, POS))\n",
    "DrawDC(dta1.Indicator(level_sets1), -3.5, 2.5, mesh, \"dta1\")\n",
    "\n",
    "level_sets2 = [-2*x + y - 4, 2*x + y + 2, -y]\n",
    "dta2 = DomainTypeArray((NEG, NEG, NEG))\n",
    "DrawDC(dta2.Indicator(level_sets2), -3.5, 2.5, mesh, \"dta2\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can construct the Union of these two domains as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta_union = TensorUnion(dta1, dta2)\n",
    "DrawDC(dta_union.Indicator(level_sets1 + level_sets2), -3.5, 2.5, mesh, \"union\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note, that a the same result would have been achieved \"by hand\" with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta3 = DomainTypeArray([(POS, POS, POS, ANY, ANY, ANY),\n",
    "                        (ANY, ANY, ANY, NEG, NEG, NEG)])\n",
    "DrawDC(dta3.Indicator(level_sets1 + level_sets2), -3.5, 2.5, mesh, \"union\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unsing `TensorIntersection` we can also generate the domain description of \"outside both triangles\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta_outside = TensorIntersection(~dta1, ~dta2)\n",
    "DrawDC(dta_outside.Indicator(level_sets1 + level_sets2), -3.5, 2.5, mesh, \"outside\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notehowever, that the description can be very long and contain non existent regions. This is because the `DomainTypeArray` does not have any information on the actual level sets used to define the domains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(dta_outside)\n",
    "print(\"Number of regions in DomainTypeArray:\", len(dta_outside))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In such cases, using the `.Compress()` memeber of `DomainTypeArray` can make sense, to remove regions which do not contibute to the domain for the given level sets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P1 = H1(mesh, order=1)\n",
    "levelsets_p1 = tuple([GridFunction(P1) for _l in level_sets1 + level_sets2])\n",
    "for lsetp1, lset in zip(levelsets_p1, level_sets1 + level_sets2):\n",
    "    InterpolateToP1(lset, lsetp1)\n",
    "\n",
    "dta_outside.Compress(levelsets_p1)\n",
    "print(\"Number of regions in DomainTypeArray after compression:\", len(dta_outside))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DrawDC(dta_outside.Indicator(levelsets_p1), -3.5, 2.5, mesh, \"outside\")"
   ]
  }
 ],
 "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.10.10"
  },
  "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
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}