aggregation.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Cell and basis aggregation in `ngsxfem`"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This unit gives a basic introduction to basis aggregation in `ngsxfem`. Especially, we treat:\n",
    "\n",
    "* definition of aggregation patches\n",
    "\n",
    "* extension of basis functions on aggregation patches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from xfem import *\n",
    "from ngsolve.webgui import *\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Introduction: Cut cells' best friends: their neighbors\n",
    "\n",
    "Let's start with a simple example (structured mesh, planar domain boundary):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from ngsolve.meshes import MakeStructured2DMesh\n",
    "mesh = MakeStructured2DMesh(quads=False, nx=3, ny=2)\n",
    "levelset = x - 0.8\n",
    "gflset = GridFunction(H1(mesh))\n",
    "InterpolateToP1(levelset, gflset)\n",
    "DrawDC(gflset, -1.0, 2.0, mesh, \"x\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In this simple configuration we have the right block of elements cut by the zero level set.\n",
    "\n",
    "The cut configuration on cut elements can become arbitrarily bad so that basis functions on the cut elements may have degenerated support \n",
    "\n",
    "* $\\leadsto$ can lead to stability or (at least) conditioning issues.\n",
    "\n",
    "To stabilize those cases we often seek for help in the **neighbor**hood. This is typically done either\n",
    "\n",
    "* through ghost penalty stabilization where essentially functions are glued together weakly in a neighborhood or\n",
    "\n",
    "* by gluing basis functions on cut elements to basis functions on uncut elements in the virtue of an extrapolation\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Patches\n",
    "To define proper neighborhood regions where support is given from \"good\" to \"bad\" elements, we will define **patches**.\n",
    "First, we will need to define:\n",
    " * which elements are \"good\" and can serve as \"root\" elements\n",
    " * and which elements are \"bad\" and need support.\n",
    "Here, we will make the simple distinction:\n",
    "$$\n",
    " \\text{uncut element} = \\text{ \"good\" element }, \\qquad\n",
    " \\text{cut element} = \\text{ \"bad\" element }\n",
    "$$  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ci = CutInfo(mesh, gflset)\n",
    "roots = ci.GetElementsOfType(NEG)\n",
    "bads = ci.GetElementsOfType(IF)\n",
    "print(\"bad element array: \", bads)\n",
    "print(\"bad elements: \", [i for i in range(len(bads)) if bads[i]])\n",
    "Draw(BitArrayCF(bads), mesh, \"bad_elements\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### `ElementAggregation`\n",
    "Next, we introduce a new class `ElementAggregation` that organizes the 'neighborhood support':\n",
    "\n",
    "Let's have a brief look at its functionality:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "help(ElementAggregation)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's create an `ElementAggregation` object and obtain a set of patches:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "EA = ElementAggregation(mesh, roots, bads)\n",
    "patch_number_field = GridFunction(L2(mesh))\n",
    "patch_number_field.vec.FV()[:] = EA.element_to_patch\n",
    "Draw(patch_number_field, mesh, \"patch_number_field\", deformation=CF((0,0,0.02*patch_number_field)))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Obviously not all elements are in patches. Indeed the root elements that are not involved in supporting a bad element form \"trivial patches\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(\"trivial root elements:\", EA.els_in_trivial_patch)\n",
    "print(\"nontrivial elements:\", EA.els_in_nontrivial_patch)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Similarly we can also relate facets that lie within a patch to the corresponding patch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(\"EA.patch_interior_facets: \", EA.patch_interior_facets)\n",
    "print(\"EA.facet_to_patch: \", EA.facet_to_patch)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "So, now we have some nice patches. These can be used to tune ghost penalty terms (reduce the number of facets where GP is applied) or to change the FEM basis which we will discuss next."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Basis aggregation"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "On each of the patches we now want to get rid of the basis functions that are not supported on \"root\" elements.\n",
    "Let \n",
    "$$\n",
    " u(x) = \\sum_{i} c_i \\varphi_i(x)\n",
    "$$\n",
    "be your usual development in your basis $\\{\\varphi\\}_i$. Now, we distinguish $\\mathcal{I}_{\\text{root}}$ and $\\mathcal{I}_{\\text{bad}}$, the index sets of basis functions with/without support on root elements. \n",
    "Then, for $i \\in \\mathcal{I}_{\\text{bad}}$ we want to express $c_i$ as a (linear) function of $(c_j)_{j \\in \\mathcal{I}_{\\text{root}}}$:\n",
    "\n",
    "We start with\n",
    "$$\n",
    " u(x) = \\sum_{i \\in \\mathcal{I}_{\\text{root}}} c_i \\varphi_i(x) + \\sum_{i \\in \\mathcal{I}_{\\text{bad}}} c_i \\varphi_i(x)\n",
    "$$\n",
    "and obtain\n",
    "$$\n",
    " u(x) = \\sum_{i \\in \\mathcal{I}_{\\text{root}}} c_i \\varphi_i(x) + \\sum_{i \\in \\mathcal{I}_{\\text{bad}}} \\sum_{c_i \\in \\mathcal{I}_{\\text{root}}} X_{ji} c_j \\varphi_i(x) = \\sum_{c_i \\in \\mathcal{I}_{\\text{root}}} c_i \\underbrace{\\sum_{j} E_{ij} \\varphi_j(x)}_{\\psi_i(x)}\n",
    "$$\n",
    "Here $\\{\\psi_i\\}_{i \\in \\mathcal{I}_{\\text{root}}}$ is the new basis obtained after aggregation. We will however characterize it always through the embedding matrix $E$.\n",
    "\n",
    "This can now be achieved by the method `AggEmbedding` which yields the embedding matrix $E$ (we explain how $E$ is obtained below). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "help(AggEmbedding)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### `L2` low order\n",
    "Let us try it out for the space of piecewise constants first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fes = L2(mesh, order=0)\n",
    "E = AggEmbedding(EA, fes)\n",
    "from helper import ShowPattern\n",
    "ShowPattern(E)\n",
    "print(E)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let us take a look at the basis functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def ExtendedBasisFunctionsAsMultiDim(EA, fes):\n",
    "    E = AggEmbedding(EA, fes)\n",
    "    gfshow = GridFunction(fes, multidim=0)\n",
    "    gf = GridFunction(fes)\n",
    "    coefvec = E.CreateRowVector()\n",
    "    for i in range(E.width):\n",
    "        coefvec[:] = 0\n",
    "        coefvec[i] = 1\n",
    "        gf.vec.data = E * coefvec\n",
    "        gfshow.AddMultiDimComponent(gf.vec)\n",
    "    return gfshow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "gfshow = ExtendedBasisFunctionsAsMultiDim(EA,fes)\n",
    "Draw (gfshow, mesh, interpolate_multidim=False, animate=False, autoscale=True);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the multidim-slider in the webgui to change the (interior) basis function and see the extended basis functions."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### `L2` higher order\n",
    "Now, higher order:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fes = L2(mesh, order=2)\n",
    "ShowPattern(AggEmbedding(EA, fes))\n",
    "gfshow = ExtendedBasisFunctionsAsMultiDim(EA,fes)\n",
    "Draw (gfshow, mesh, interpolate_multidim=False, animate=False, autoscale=True);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### `H1`\n",
    "When considering spaces with (partial) continuity, some dofs may appear in several patches but not on a root element. In this case the extensions are average:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=1)\n",
    "ShowPattern(AggEmbedding(EA, fes))\n",
    "gfshow = ExtendedBasisFunctionsAsMultiDim(EA,fes)\n",
    "Draw (gfshow, mesh, interpolate_multidim=False, animate=False, autoscale=True);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Vector-valued spaces\n",
    "And now, `HDiv`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fes = HDiv(mesh, order=1)\n",
    "ShowPattern(AggEmbedding(EA, fes))\n",
    "gfshow = ExtendedBasisFunctionsAsMultiDim(EA,fes)\n",
    "Draw (gfshow, mesh, interpolate_multidim=False, animate=False, autoscale=True, vectors = { \"grid_size\" : 20});"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, if you take a look at the basis functions which have non-trivial aggregation extension (e.g. with `multidim=9` in the GUI slider), you may wonder how this is a reasonable extension. There are however two effects that you have to keep in mind:\n",
    "\n",
    "* for `dof`s that are shared by several patches the embedding applies an average, i.e. the non-smoothness of extended basis functions on the patch can result from that. \n",
    "* Next, there is a finite element blending away from patches due to divergence-conformity (normal-continuity). \n",
    "\n",
    "To illustrate this, it may be helpful to replace `HDiv(...)` with `Discontinuous(HDiv(...))` to investigate the same (local) basis functions without normal-continuity."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### How it works\n",
    "\n",
    "The mechanism behind the patchwise extension is very generic and obviously works for `L2`, `H1`, `HDiv`, ... in a generic way. But how?\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "#### A formulation as a harmonic extension\n",
    "Let $w_h(x) = \\sum_{i \\in \\mathcal{I}_{\\text{bad}}} c_i \\varphi_i(x)$ be a finite element function only supported on bad elements and $W_h$ the corresponding subspace of the underlying FE space $V_h$ and $X_h = W_h \\setminus V_h$.\n",
    "\n",
    "\n",
    "The idea is now that for every basis function $\\varphi_i$ with $i \\in \\mathcal{I}_{\\text{root}}$ we define its extension into $W_h$ through\n",
    "$$\n",
    "\\operatorname{argmin}_{w_h \\in W_h} \n",
    "\\vert \\varphi_i + w_h \\vert_\\ast\n",
    "$$\n",
    "where $\\vert \\cdot \\vert_\\ast$ is a smoothness-measuring semi-norm. $w_h$ is like a harmonic extension from the \"root\" into the \"bad\" elements."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "#### A formulation as a patch-wise harmonic extension\n",
    "To localize this problem we proceed patch-by-patch:\n",
    "\n",
    "For every patch $\\omega$ we define $W_h^{\\omega}$ through functions that are supported on the patch, but not on a root element. Note that some basis function may appear in different patches now. \n",
    "\n",
    "Then, we solve the local patch problems:\n",
    "For every basis function $\\varphi_i$ with $i \\in \\mathcal{I}_{\\text{root}}$ and support on the patch $\\omega$ we define its extension into $W_h^\\omega$ through\n",
    "$$\n",
    "\\operatorname{argmin}_{w_h \\in W_h^\\omega} \n",
    "\\vert \\varphi_i + w_h \\vert_{\\ast,\\omega}\n",
    "$$\n",
    "where $\\vert \\cdot \\vert_{\\ast,\\omega}$ is a smoothness-measuring semi-norm. \n",
    "\n",
    "The thusly obtained extended functions are afterwards averaged.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "#### A generic semi-norm $\\vert \\varphi_i + w_h \\vert_{\\ast,\\omega}$\n",
    "\n",
    "It now remains to define a proper semi-norm. \n",
    "\n",
    "When using `AggEmbedding` as above the chosen semi-norm is defined through a ghost penalty term inside all patches (Note that patch-wise polynomials form the kernel of that operator):\n",
    "```\n",
    "def AggEmbedding(EA, fes, deformation=None, heapsize=1000000):\n",
    "    \"\"\"...\"\"\"\n",
    "    u,v = fes.TnT()\n",
    "    ghost_penalty = (u - u.Other()) * (v - v.Other()) * dFacetPatch(deformation=deformation)\n",
    "    return ExtensionEmbedding(EA, fes, ghost_penalty, heapsize=heapsize)\n",
    "```\n",
    "\n",
    "If you want to use a different semi-norm for the extension, you can use `ExtensionEmbedding` accordingly. \n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Applications of patch aggregations\n",
    "\n",
    "We now have learned about how to form patches and extend basis functions. This can be exploited in several ways:\n",
    "\n",
    " 1. To solve CutFEM problems where the resulting system is reduced to the root unknowns by the embedding, cf. the `fictdom_aggfem.py` demo of `ngsxfem`.\n",
    " \n",
    " 2. Patches can be used to minimize the facets that are used in ghost penalty stabilizations, cf. e.g. [unfmixed example](unfmixed.ipynb).\n",
    " \n",
    " 3. Patches can be used to solve local patch-wise problems, e.g. for postprocessings in unfitted mixed methods, cf. [unfmixed example](unfmixed.ipynb)."
   ]
  }
 ],
 "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.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
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}