basics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geometry description with the level set function\n",
    "\n",
    "This unit gives a basic introduction to a few fundamental concepts in `ngsxfem`. Especially, we treat:\n",
    "\n",
    "* level set geometries and its (P1) approximation (for higher order approximations see [intlset.ipynb](intlset.ipynb))\n",
    "* Visualization of discontinuous (across level set zeros) functions\n",
    "* handling of cut configurations inside the mesh (active elements / facets / dofs)\n",
    "* Restriction of `FESpace`s, `integrators` and `Matrices` to active parts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Most features of the `ngsxfem` add-on to NGSolve are based on a level set description.\n",
    "\n",
    "Let $\\phi$ be a continuous level set function. \n",
    "\n",
    "We define the domains and the interface\n",
    "\n",
    "$$\n",
    "\\begin{array}{rccc}\n",
    "  \\text{domain:} & \\Omega_{-} := \\{ \\phi < 0 \\}, & \\quad\n",
    "  \\Omega_{+} := \\{ \\phi > 0 \\}, & \\quad\n",
    "  \\Gamma := \\{ \\phi = 0 \\}. \\\\\n",
    "  \\texttt{ngsxfem} \\text{ keyword:} & \\texttt{NEG} & \n",
    "  \\texttt{POS} & \n",
    "  \\texttt{IF}  \n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let us import `Netgen/NGSolve` stuff first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#import netgen.gui\n",
    "# ngsolve stuff\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import *\n",
    "# basic geometry features (for the background mesh)\n",
    "from netgen.occ import *\n",
    "# visualization stuff\n",
    "from ngsolve.internal import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We add the basic functionality of the `ngsxfem` package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# basic xfem functionality\n",
    "from xfem import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The background mesh\n",
    "We generate the background mesh of the domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "geo = OCCGeometry(unit_square_shape.Scale((0,0,0),3).Move((-1.5,-1.5,0)), dim=2)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.5,quad_dominated=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "On the background mesh we define the level set function. In the visualization -- using `DrawDC` -- we distinguish the negative and the positive part of the level set values corresponding to the subdomains that are described."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "levelset = (sqrt(sqrt(x**4+y**4)) - 1.0)\n",
    "DrawDC(levelset, -1.0, 2.0, mesh, \"x\")\n",
    "#help(DrawDC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Geometry approximation through level set approximation\n",
    " * Arbitrary level set functions can be difficult to handle (cut topology)\n",
    " * To obtain something feasible, we interpolate into the space of piecewise linear functions\n",
    " * The thusly obtained approximation is straight within each (simplex) element\n",
    " * Simplifies numerical integration\n",
    " * We will discuss enhancements later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "lsetp1 = GridFunction(H1(mesh,order=1))\n",
    "InterpolateToP1(levelset,lsetp1)\n",
    "DrawDC(lsetp1,-1,1,mesh,\"lsetp1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* `InterpolateToP1` takes vertex values (in contrast to `Set`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Cut Information"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "In unfitted FEM it is important to know which elements are cut by the interface $\\{ \\phi = 0 \\}$. These informations are gathered in the cut information class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "ci = CutInfo(mesh, lsetp1)\n",
    "help(CutInfo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Marking elements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "### `DOMAIN_TYPE`\n",
    "We can ask for BitArrays corresponding to the differently marked elements \n",
    "\n",
    "  * `NEG` : True if $\\phi < 0$ everywhere on this element, else False\n",
    "  * `POS` : True if $\\phi > 0$ everywhere on this element, else False\n",
    "  * `IF` : True if $\\phi = 0$ somewhere on this element, else False  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dts = [IF,NEG,POS]\n",
    "def DrawMarkedElements():\n",
    "    for dt in dts: \n",
    "        print(\"Drawing\", dt); \n",
    "        Draw(BitArrayCF(ci.GetElementsOfType(dt)),mesh,\"elements_\"+str(dt))\n",
    "        yield\n",
    "dme = DrawMarkedElements()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "next(dme)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "next(dme)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "next(dme)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### `COMBINED_DOMAIN_TYPE`\n",
    "There are also some predefined combinations of the previous markers:\n",
    "  \n",
    "| IF|POS|NEG| combined flag | description    |\n",
    "|---|---|---|---------------|----------------|\n",
    "|  0|  0|  0|            NO | False on all elements\n",
    "|  0|  0|  1|           NEG | True if $\\phi < 0$ everywhere on this element, else False   \n",
    "|  0|  1|  0|           POS | True if $\\phi > 0$ everywhere on this element, else False\n",
    "|  0|  1|  1|         UNCUT | True if $\\phi \\neq 0$ everywhere on this element, else False   \n",
    "|  1|  0|  0|            IF | True if $\\phi = 0$ somewhere on this element, else False  \n",
    "|  1|  0|  1|        HASNEG | True if $\\phi \\leq 0$ somewhere on this element, else False\n",
    "|  1|  1|  0|        HASPOS | True if $\\phi \\geq 0$ somewhere on this element, else False\n",
    "|  1|  1|  1|           ANY | True on all elements\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "cdts = [NO,CDOM_NEG,CDOM_POS,UNCUT,CDOM_IF,HASNEG,HASPOS,ANY]\n",
    "for dt in cdts: print(int(dt),\": {0:03b} ,\".format(int(dt)),str(dt) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "for dt in cdts: print(dt);Draw(BitArrayCF(ci.GetElementsOfType(dt)),mesh,\"elements_\"+str(dt),min=0,max=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Updating the CutInformation\n",
    "* level set function can be changed with `Update`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "levelset = sqrt((x-0.5)**2+y**2)-1\n",
    "scene = Draw(levelset,mesh,\"levelset\") # overwrite visualization\n",
    "InterpolateToP1(levelset,lsetp1)\n",
    "ci.Update(lsetp1);scene.Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Marking with facets \n",
    " * `GetFacetsWithNeighborTypes`\n",
    " * `GetElementsWithNeighborFacets`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Sometimes want to mark some facets depending on their neighbor states (NEG/POS/IF)\n",
    "* Needed for some stabilizations in unfitted FEM (CutFEM/XFEM/..)\n",
    "* Facet-marking based on two element-based BitArrays. \n",
    "* There are two options to connect two BitArrays A and B:\n",
    " * use_and=True (default): the result is True for a facet if (A(left) and B(right)) or (B(left) and A(right))\n",
    " * use_and=False     (or): the result is True for a facet if (A(left) or B(right)) or (B(left) or A(right))\n",
    " \n",
    "| A(left) | B(right) | partial result (`use_and=True`) | partial result (`use_and=False`) |\n",
    "|----------|-----------|-----------------------|------------------------|\n",
    "| False    | False     | False                 | False                  |\n",
    "| False    | True      | False                 | True                   |\n",
    "| True     | False     | False                 | True                   |\n",
    "| True     | True      | True                  | True                   |\n",
    " \n",
    "* with bnd_val_a and bnd_val_b (default: True) we can prescribe how the BitArray should be interpreted if there is only one neighbor (boundary)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 1:\n",
    "\n",
    "* All facets that are on the outer boundary of the NEG-domain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#help(GetFacetsWithNeighborTypes)\n",
    "ba_facets = GetFacetsWithNeighborTypes(mesh,\n",
    "                                       a=ci.GetElementsOfType(NEG),\n",
    "                                       b=ci.GetElementsOfType(IF))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* We can also go from facets to elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#help(GetElementsWithNeighborFacets)\n",
    "ba_surround_facets = GetElementsWithNeighborFacets(mesh,ba_facets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def SelectFacetsToDraw(dt1,dt2): \n",
    "    ba_facets = GetFacetsWithNeighborTypes(mesh,a=ci.GetElementsOfType(dt1), b=ci.GetElementsOfType(dt2))\n",
    "    ba_surround_facets = GetElementsWithNeighborFacets(mesh,ba_facets)\n",
    "    Draw(BitArrayCF(ba_surround_facets), mesh, \"surrounding_facets\")    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "SelectFacetsToDraw(CDOM_IF,CDOM_NEG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 2:\n",
    "We can successively apply this to extend from some part of the domain by direct neighbors, etc .... \n",
    "\n",
    "(execute the next cell several times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "ba_facets = GetFacetsWithNeighborTypes(mesh,a=ba_surround_facets,b=ba_surround_facets,\n",
    "                                           bnd_val_a=False,bnd_val_b=False,use_and=False)\n",
    "ba_surround_facets = GetElementsWithNeighborFacets(mesh,ba_facets)\n",
    "Draw(BitArrayCF(ba_surround_facets), mesh, \"ba_surround\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now learned how to mark elements and facets. Next, we want to apply this to adapt finite element formulations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Restriction to *active* mesh/`dof`s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "* Unfitted FEM: Solve a PDE problem *restricted* to a subdomain, e.g.\n",
    "  * a fictitious domain method on $\\Omega_{-} = \\{ \\phi < 0 \\}$ or\n",
    "  * a trace finite element on $\\Gamma = \\{ \\phi = 0\\}$\n",
    "* Standard `FESpace` but only a selection of `dof`s (`active_dofs`)\n",
    "* We have helper functions for that as well"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Marking dofs (`GetDofsOfElements`)\n",
    "With `GetDofsOfElements` we can select a subset of dofs corresponding to a set of elements.\n",
    "With `Compress` we can reduce a finite element space to the corresponding dofs:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vhbase = H1(mesh, order=1, dirichlet=[])\n",
    "VhC = Compress(Vhbase,GetDofsOfElements(Vhbase,ci.GetElementsOfType(HASNEG)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, you can directly `Restrict` a finite element space to a set of elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VhR = Restrict(Vhbase,ci.GetElementsOfType(HASNEG))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can compare:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(Vhbase.ndof, VhC.ndof, VhR.ndof) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Compress` and `Restrict` yield the same amount of dofs. However, there is a subtle difference: `VhR` is not defined away from the marked elements whereas `VhC` can be seen as the background `FESpace` with the inactive dofs set to zero. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Next:\n",
    "\n",
    " * `Compress` and `Restrict` allow to *restrict* vectors and matrices to \"active\" subset of the background FE space\n",
    " * What to do if you want to set up the finite element formulation only on a selection of elements / facets / unknowns?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Marking integrators (`definedonelements`)\n",
    "We can define integrators on a set of elements only (this is a standard NGSolve feature):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "integrator = VhR.TrialFunction()*VhR.TestFunction() * dx(definedonelements=ci.GetElementsOfType(HASNEG))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ci.GetElementsOfType(HASNEG))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Note that the BitArray taken in the `definedonelements` argument is taken by reference, i.e. later changes to the `CutInfo` `ci` will automatically change the set of elements the `integrator` will act on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Restricted Bilinearforms\n",
    "Furthermore, to avoid that matrices are set up for the full standard FESpace we can restrict matrix couplings to elements (and facets for dg terms). This may be attractive if you do not work with `Restrict`ed or `Compress`ed finite element spaces (e.g. because of changing domains) but want the linear systems to be of reduced size. Further, in unfitted FEM often stabilizations of \"Ghost Penalty\" type are applied which may lead to dgjumps, i.e. couplings of dofs over facets. To avoid that these couplings are introduced on the whole domain, you may want to restrict the facet couplings explicitly. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#help(RestrictedBilinearForm)\n",
    "a_full = BilinearForm(Vhbase, check_unused=False)\n",
    "a = RestrictedBilinearForm(Vhbase,element_restriction=ci.GetElementsOfType(IF),\n",
    "                           facet_restriction=None, check_unused=False)\n",
    "a_full.Assemble()\n",
    "a.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* This does not change the dimension:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(a.mat.height,a.mat.width)\n",
    "print(a_full.mat.height,a_full.mat.width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* But the number on non-zero entries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(len(a_full.mat.AsVector()))\n",
    "print(len(a.mat.AsVector()))"
   ]
  },
  {
   "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
}