{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8.2 Integration on level set domains\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 \n",
    " * Cut differential symbols to conveniently define integrals on cut domain \n",
    " * for higher order level set geometry approximations "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    " * simple example: calculate area of the unit circle using level set representations\n",
    " * background domain: $[-1.5, 1.5]^2 \\subset \\mathbb{R}^2$.\n",
    "\n",
    "Let $\\phi$ be a continuous level set function. We recall:\n",
    "$$\n",
    "  \\Omega_{-} := \\{ \\phi < 0 \\}, \\quad\n",
    "  \\Omega_{+} := \\{ \\phi > 0 \\}, \\quad\n",
    "  \\Gamma := \\{ \\phi = 0 \\}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# basic geometry features (for the background mesh)\n",
    "from netgen.geom2d import SplineGeometry\n",
    "# ngsolve\n",
    "from ngsolve import *\n",
    "# basic xfem functionality\n",
    "from xfem import *\n",
    "# for isoparametric mapping\n",
    "from xfem.lsetcurv import *\n",
    "# visualisation\n",
    "from ngsolve.webgui import *\n",
    "# the contant pi\n",
    "from math import pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We generate the background mesh of the domain and use a simplicial triangulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "square = SplineGeometry()\n",
    "square.AddRectangle([-1.5,-1.5], [1.5,1.5], bc=1)\n",
    "mesh = Mesh(square.GenerateMesh(maxh=0.8, quad_dominated=False))\n",
    "Draw(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "On the background mesh we define the level set function. In this example we choose a levelset function which measures the distance to the origin and substracts a constant $r$, resulting in an circle of radius $r$ as $\\Omega_{-}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "r = 1\n",
    "levelset = sqrt(x**2+y**2) - r\n",
    "DrawDC(levelset, -3.5, 2.5, mesh, 'levelset')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* Area of the unit circle:\n",
    "$$\n",
    "A = \\pi = \\int_{\\Omega_-(\\phi)} 1 \\, \\mathrm{d}x\n",
    "$$\n",
    "\n",
    "## Piecewise linear geometry approximation\n",
    "\n",
    "* The ngsxfem integration on implicitly defined geometries is only able to handle functions $\\phi$ which are in $\\mathcal{P}^1(T)$ for all $T$ of the triangulation.\n",
    "\n",
    "* Therefore we replace the levelset function $\\phi$ with an approximation $\\phi^\\ast \\in \\mathcal{P}^1(T), T \\subset \\Omega$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "V = H1(mesh,order=1,autoupdate=True)\n",
    "lset_approx = GridFunction(V,autoupdate=True)\n",
    "InterpolateToP1(levelset, lset_approx)\n",
    "DrawDC(lset_approx, -3.5, 2.5, mesh, 'levelset_p1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Notice that this replacement $\\phi \\mapsto \\phi^\\ast$ will introduce a geometry error of order $h^2$. We will investigate this numerically later on in this notebook.\n",
    "\n",
    "As a next step we need to define the integrand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = CoefficientFunction(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As in NGSolve, we integrate by using the `Integrate` function. For this to be aware of the fact that we are integration over a level set domain we use the `dCut` differential symbol, to which we need to pass: \n",
    "\n",
    " * `levelset`: level set function which describes the geometry (best choice: P1-GridFunction)\n",
    " * `domain_type`: decision on which part of the geometry the integration should take place (`NEG`/`POS`/`IF`)\n",
    " * `order`: The order of integration rule to be used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "order = 2\n",
    "integral = Integrate(f * dCut(levelset=lset_approx, domain_type=NEG, order=order), mesh=mesh)\n",
    "error = abs(integral - pi)\n",
    "print(\"Result of the integration: \", integral)\n",
    "print(\"Error of the integration: \", error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Convergence study\n",
    "\n",
    "To get more accurate results\n",
    "\n",
    " * save the previous error and execute a refinement \n",
    " * repeat the steps from above\n",
    " * Repeatedly execute the lower block to add more refinement steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "errors = []\n",
    "\n",
    "def AppendResultAndPostProcess(integral):\n",
    "    error = abs(integral - pi)\n",
    "    print(\"Result of the integration: \", integral)\n",
    "    print(\"Error of the integration: \", error)\n",
    "\n",
    "    errors.append(error)\n",
    "    eoc = [log(errors[i+1]/errors[i])/log(0.5) for i in range(len(errors)-1)]\n",
    "\n",
    "    print(\"Collected errors:\", errors)\n",
    "    print(\"experimental order of convergence (L2):\", eoc)    \n",
    "\n",
    "AppendResultAndPostProcess(integral)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "scene = Draw(lset_approx, mesh, \"lsetp1\", min=0, max=0.0)\n",
    "def RefineAndIntegrate():\n",
    "    # refine cut elements only:\n",
    "    RefineAtLevelSet(gf=lset_approx)\n",
    "    mesh.Refine()\n",
    "\n",
    "    InterpolateToP1(levelset,lset_approx)\n",
    "    scene.Redraw()\n",
    "    integral = Integrate(f * dCut(levelset=lset_approx, domain_type=NEG, order=order), mesh=mesh)\n",
    "    AppendResultAndPostProcess(integral)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RefineAndIntegrate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We observe only a second order convergence which results from the piecewise linear geometry approximation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Higher order geometry approximation with an isoparametric mapping\n",
    "In order to get higher order convergence we can use the isoparametric mapping functionality of xfem.\n",
    "\n",
    "We apply a mesh transformation technique in the spirit of isoparametric finite elements:\n",
    "![title](graphics/lsetcurv.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Video of the mesh transformation**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%html\n",
    "<iframe width=100% height=600px src=\"https://www.youtube.com/embed/Mst_LvfgPCg?rel=0\" frameborder=\"0\" allow=\"autoplay; encrypted-media\" allowfullscreen></iframe>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To compute the corresponding mapping we use the LevelSetMeshAdaptation class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mesh = Mesh(square.GenerateMesh(maxh=0.4, quad_dominated=False))\n",
    "lsetmeshadap = LevelSetMeshAdaptation(mesh, order=2, threshold=1000, discontinuous_qn=True)\n",
    "deformation = lsetmeshadap.CalcDeformation(levelset)\n",
    "Draw(deformation,mesh,\"deformation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can observe the geometrical improvement in the following sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "scene1 = Draw(deformation, mesh, \"deformation\", autoscale=False, min=0.0, max=0.028)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scene2 = DrawDC(lsetmeshadap.lset_p1, -3.5, 2.5, mesh, \"lsetp1_ho\", deformation=deformation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=100\n",
    "\n",
    "deformation.vec.data = 1.0/N*deformation.vec\n",
    "for i in range (1, N+1):\n",
    "    deformation.vec.data = (i+1)/i * deformation.vec\n",
    "    scene1.Redraw()\n",
    "    scene2.Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Convergence study with isoparametric mapping applied\n",
    "We now apply the mesh deformation (2nd order mapping) to improve the accuracy of the integration problem.\n",
    "\n",
    "**error tables**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "order = 2\n",
    "mesh = Mesh (square.GenerateMesh(maxh=0.7, quad_dominated=False))\n",
    "\n",
    "levelset = sqrt(x*x+y*y)-1\n",
    "referencevals = { POS : 9-pi, NEG : pi, IF : 2*pi }\n",
    "\n",
    "lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=0.2, discontinuous_qn=True)\n",
    "lsetp1 = lsetmeshadap.lset_p1\n",
    "errors_uncurved = dict()\n",
    "errors_curved = dict()\n",
    "eoc_uncurved = dict()\n",
    "eoc_curved = dict()\n",
    "\n",
    "for key in [NEG,POS,IF]:\n",
    "    errors_curved[key] = []\n",
    "    errors_uncurved[key] = []\n",
    "    eoc_curved[key] = []\n",
    "    eoc_uncurved[key] = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**refinements**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = CoefficientFunction (1.0)\n",
    "refinements = 5\n",
    "for reflevel in range(refinements):\n",
    "    if(reflevel > 0):\n",
    "        mesh.Refine()\n",
    "\n",
    "    for key in [NEG,POS,IF]:\n",
    "        # Applying the mesh deformation\n",
    "        deformation = lsetmeshadap.CalcDeformation(levelset)\n",
    "\n",
    "        integrals_uncurved = Integrate(f * dCut(levelset=lsetp1, domain_type=key, order=order), mesh=mesh)\n",
    "        # dCut deals with the mesh deformation internally\n",
    "        integrals_curved = Integrate(f * dCut(levelset=lsetp1, domain_type=key, order=order, deformation=deformation), mesh=mesh)\n",
    "        \n",
    "        errors_curved[key].append(abs(integrals_curved - referencevals[key]))\n",
    "        errors_uncurved[key].append(abs(integrals_uncurved - referencevals[key]))\n",
    "    # refine cut elements:\n",
    "    RefineAtLevelSet(gf=lsetmeshadap.lset_p1)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Convergence rates**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "for key in [NEG,POS,IF]:\n",
    "    eoc_curved[key] = [log(a/b)/log(2) for (a,b) in zip (errors_curved[key][0:-1],errors_curved[key][1:]) ]\n",
    "    eoc_uncurved[key] = [log(a/b)/log(2) for (a,b) in zip (errors_uncurved[key][0:-1],errors_uncurved[key][1:]) ]\n",
    "\n",
    "print(\"errors (  curved):  \\n{}\\n\".format(  errors_curved))\n",
    "print(\"   eoc (  curved):  \\n{}\\n\".format(     eoc_curved))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note: The integration on level set domains (and the mapping) can be extended to bilinear and linear form integrators with the same syntax."
   ]
  }
 ],
 "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.9.1-final"
  },
  "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
}