unfmixed.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "80685a63",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Unfitted Mixed FEM (using `HDiv` FE spaces)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4ec015b",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Our goal is to consider the mixed Poisson problem with Dirichlet boundary conditions in an unfitted setting: <br> \n",
    "Find $u: \\Omega \\rightarrow \\mathbb{R}^d$, $p : \\Omega \\rightarrow \\mathbb{R}$ such that\n",
    "\n",
    "\\begin{alignat}{2}\n",
    "u - \\nabla p &= \\phantom{-} 0 \\quad &&\\text{in } \\Omega, \\tag{C1}\\\\\n",
    "\\operatorname{div} u &= -f \\quad &&\\text{in } \\Omega, \\tag{C2} \\\\\n",
    "p &= \\phantom{-} p_D \\quad &&\\text{on } \\partial \\Omega. \\tag{C3}\n",
    "\\end{alignat}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbbee49b",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Here $\\Omega \\subset \\mathbb{R}^d$ is a bounded domain that is not parametrized by the computational mesh, but contained in a background domain $\\widetilde{\\Omega}$ such that $\\overline \\Omega \\subseteq \\widetilde{\\Omega}$ which is triangulated by a simplicial, shape regular and quasi-uniform triangulation $\\widetilde{\\mathcal{T}_h}$. <br>\n",
    "We set \n",
    "\\begin{align}\n",
    "\\mathcal{T}_h &= \\big\\{T \\in \\widetilde{\\mathcal{T}_h} : \\text{meas}(T \\cap \\Omega) > 0 \\big\\}, \\tag{T1} \\\\\n",
    "\\mathcal{T}_h^{\\text{cut}} &= \\big\\{T \\in \\widetilde{\\mathcal{T}_h} : \\text{meas}(T \\cap \\partial \\Omega) > 0 \\big\\}, \\tag{T2}\\\\\n",
    "\\mathcal{T}_h^{\\text{int}} &= \\mathcal{T}_h \\setminus \\mathcal{T}_h^{\\text{cut}}. \\tag{T3}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "447e7abe",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The Method"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30824896",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To obtain inf-sup stable unfitted discretizations, one usually add stabilization terms like the ghost penalty stabilization. In the case of mixed problems, however, adding such terms on the $p$-$q$-coupling block would pollute the mass balance. Instead, one can either use a ghost penalty on the $\\operatorname{div}$- and $\\nabla$-block as done e.g.\n",
    "in [this preprint(arXiv:2205.12023)](https://arxiv.org/abs/2205.12023) or change the domain on which \n",
    "the $\\operatorname{div}$- and $\\nabla$-blocks act as we have done in \n",
    "[this preprint(arXiv:2306.12722)](https://arxiv.org/abs/2306.12722). We showed that one can consider the following weak formulation: <br>\n",
    "Find $(u_h,\\bar p_h) \\in \\mathbb{RT}^k(\\mathcal{T}_h) \\times \\mathbb{P}^k(\\mathcal{T}_h)$ such that \n",
    "\\begin{alignat}{2}\n",
    "a_h(u_h,v_h) + b_h(v_h,\\bar p_h) &= (v_h \\cdot n, p_D)_{\\partial \\Omega} \\qquad &&\\forall v_h \\in \\mathbb{RT}^k(\\mathcal{T}_h), \\tag{M1} \\\\\n",
    "b_h(u_h,q_h) &= - (f_h,q_h)_{\\mathcal{T}_h} \\qquad &&\\forall q_h \\in \\mathbb{P}^k(\\mathcal{T}_h), \\tag{M2}\n",
    "\\end{alignat}\n",
    "where \n",
    "\\begin{align}\n",
    "a_h(u_h,v_h) &:= (u_h,v_h)_{\\Omega} + \\gamma_u j_h (u_h,v_h), \\tag{A}\\\\\n",
    "b_h(v_h,p_h) &:= (\\operatorname{div} v_h,p_h)_{\\mathcal{T}_h}. \\tag{B}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1052d5cb",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The following code imports the basic functionality of netgen, ngsolve and ngsxfem. Afterwards, we set some general parameters for this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99a55492",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#basic imports\n",
    "from netgen.occ import *\n",
    "from ngsolve import *\n",
    "from ngsolve.internal import *\n",
    "from xfem import *\n",
    "from xfem.lsetcurv import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d4d9182",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We consider a ring geometry that is contained inside the square $[-1,1]^2$ with inner radius $r_1 = 0.25$ and outer radius $r_2 = 0.75$. Then, the geometry is described by the following signed distance function <br> <br>\n",
    "\\begin{equation}\n",
    "\\Phi(x) = \\Big\\vert \\sqrt{x_1^2 + x_2^2} - \\frac{1}{2}(r_1+r_2) \\Big\\vert - \\frac{r_2-r_1}{2}.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90b5384f",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "square = OCCGeometry(unit_square_shape.Scale((0,0,0),2).Move((-1,-1,0)), dim=2)\n",
    "mesh = Mesh(square.GenerateMesh(maxh=0.2))\n",
    "r1, r2 = 1/4, 3/4  # inner/ outer radius\n",
    "rr, rc = (r2 - r1) / 2.0 , (r1 + r2) / 2.0\n",
    "r = sqrt(x**2 + y**2)\n",
    "levelset = IfPos(r - rc, r - rc - rr, rc - r - rr)\n",
    "DrawDC(levelset,-1.0,2.0,mesh,\"x\") "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdc4292a",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Discretization parameters\n",
    "\n",
    "We choose discretization order and ghost penalty stabilization parameter (0 for no stabilization):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c54ce818",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "order = 2\n",
    "gamma_stab = 0 # 0 if no ghost penalty is to be used"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64c98663",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Geometry approximation\n",
    "Next, we construct the mesh deformation for higher order accuracy (see the `intlset`-tutorial for further details):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06652db0",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order+1, threshold=0.1,discontinuous_qn=True)\n",
    "deformation = lsetmeshadap.CalcDeformation(levelset)\n",
    "lsetp1 = lsetmeshadap.lset_p1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc63831f",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We will need several different element and facet markers for the setup of the method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1448fa2",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "ci = CutInfo(mesh, lsetp1)\n",
    "hasneg = ci.GetElementsOfType(HASNEG) # elements with negative level set (root elements)\n",
    "pos = ci.GetElementsOfType(POS)       # outside elements\n",
    "hasif = ci.GetElementsOfType(IF)      # cut elements\n",
    "hasany = ci.GetElementsOfType(ANY)    # all elements (trivial array)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fd4c882",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We define patches (for minimizing GP stabilization (if used) and postprocessings). Here, we mark all cut elements as \"bad\", i.e. they need to be supported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a82f9425",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "EA = ElementAggregation(mesh)\n",
    "EA.Update(ci.GetElementsOfType(NEG), hasif)\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)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b63cbf9",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Solving the problem "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f43ebe07",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We consider the exact solution $p = \\sin(x)$ and derive $u = \\nabla p$ and $f = - \\Delta p$. \n",
    "For simplicity, we assume here that we are given a proper extension of $f$ from $\\Omega$ to the active mesh by simply evaluating the same closed form expression. This can be generalized to deal with functions only defined on $\\Omega$. This is discussed in a separate section below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fb8b2a2",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "p_exact = sin(x)\n",
    "coeff_f = - (p_exact.Diff(x).Diff(x) + p_exact.Diff(y).Diff(y)).Compile()\n",
    "u_exact = CF((p_exact.Diff(x),p_exact.Diff(y)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7d149ed",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Finite element spaces\n",
    "\n",
    "We use a Raviart-Thomas finite element space for the flux and a standard `L2` space for the scalar variable, both are restricted to the active mesh. Note, that if we use ghost penalty stabilization we have to pay with additional couplings on cut elements, i.e. we need to set `dgjumps=True` in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8c367a9",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#FESpaces\n",
    "Shsbase = HDiv(mesh, order=order, dirichlet=[], dgjumps=(gamma_stab > 0), RT=True)\n",
    "Vhbase = L2(mesh, order=order, dirichlet=[])\n",
    "Vh = Restrict(Vhbase, hasneg)\n",
    "Shs = Restrict(Shsbase, hasneg)\n",
    "Wh = Shs*Vh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66b76065",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### GridFunction, trial- and test functions, and some standard coefs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a62396a6",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfw = GridFunction(Wh)\n",
    "gfu, gfpT = gfw.components[0:2]\n",
    "(uh,pT), (vh,qT) = Wh.TnT()\n",
    "\n",
    "h = specialcf.mesh_size     # mesh size\n",
    "n = Normalize(grad(lsetp1)) # normal to level set"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90857f99",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Differential symbols for different integration domains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12a7317d",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# active mesh integrals:\n",
    "dxbar = dx(definedonelements=ci.GetElementsOfType(HASNEG), deformation=deformation)\n",
    "# uncut mesh integrals:\n",
    "dxinner = dx(definedonelements=ci.GetElementsOfType(NEG), deformation=deformation)\n",
    "# cut element integrals (full elements, no cut integral):\n",
    "dxcut = dx(definedonelements=ci.GetElementsOfType(IF), deformation=deformation)\n",
    "# integral on zero level set:\n",
    "ds = dCut(lsetp1, IF, definedonelements=hasif, deformation=deformation)\n",
    "# integral on Omega (physical domain, this is a cut integral):\n",
    "dX = dCut(lsetp1, NEG, definedonelements=hasneg, deformation=deformation)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d96764ea",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If ghost penalties are used, we further need a facet patch integral on some facet patches. Here, we only stabilize on facets that lie within a stabilizing patch, cf. [aggregation tutorial](aggregation.ipynb) for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2c62cb4",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "if gamma_stab > 0:\n",
    "    #Integration domain for ghost penalty\n",
    "    dP = dFacetPatch(definedonelements=EA.patch_interior_facets, deformation=deformation)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e8e1411",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Next, we create the `BilinearForm` and again distinguish the ghost penalty case from the unstabilized case, because the GP case require a larger sparsity pattern, however only on the stabilizing patches:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb97da32",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "if gamma_stab > 0:\n",
    "    a = RestrictedBilinearForm(Wh, element_restriction=hasneg, \n",
    "                               facet_restriction=EA.patch_interior_facets, symmetric=True)\n",
    "else:\n",
    "    a = BilinearForm(Wh, symmetric=True,condense=False)\n",
    "f = LinearForm(Wh)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1511082a",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Definition of variational formulation:\n",
    "Now, we implement the terms in $(M1)$ and $(M2)$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63175448",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a += uh*vh * dX\n",
    "a += (div(uh) * qT + div(vh) * pT) * dxbar\n",
    "if gamma_stab > 0: # ghost penalty\n",
    "    a += gamma_stab * (uh - uh.Other()) * (vh - vh.Other()) * dP\n",
    "f += -coeff_f * qT * dxbar\n",
    "f += p_exact * vh * n * ds\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0011755",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Assembly and solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce6d5818",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a.Assemble()\n",
    "f.Assemble()\n",
    "gfw.vec.data = a.mat.Inverse(Wh.FreeDofs()) * f.vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3a488a6",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "Draw(BitArrayCF(hasneg)*gfw.components[0],mesh,\"u\",deformation=deformation, min=0, max=1.2) # flux variable on active mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3615dad9",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "DrawDC(lsetp1,gfw.components[1],-1,mesh,\"p\",deformation=deformation, min=-1, max=1)        # p on physical domain"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9594cc79",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can now measure the error and considered different measures: \n",
    "\n",
    " * `p_l2error = ` $\\Vert p - p_h \\Vert_{\\Omega}$ (the error on the physical domain)\n",
    "\n",
    " * `p_inner_l2error = ` $\\Vert p - p_h \\Vert_{\\mathcal{T}_h^{\\text{int}}}$ (the error only on uncut elements)\n",
    "\n",
    " * `u_l2error = ` $\\Vert u - u_h \\Vert_{\\Omega}$ (the error on the physical domain)\n",
    "\n",
    " * `u_l2error_bar = ` $\\Vert u - u_h \\Vert_{\\mathcal{T}_h}$ (the error on the whole active mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c854ca76",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "p_l2error = sqrt(Integrate((gfpT - p_exact)**2*dX.order(2*order+3), mesh))\n",
    "p_inner_l2error = sqrt(Integrate((gfpT - p_exact)**2*dxinner, mesh))\n",
    "u_l2error = sqrt(Integrate((gfu - u_exact)**2*dX.order(2*order+3), mesh))\n",
    "u_l2error_bar = sqrt(Integrate((gfu - u_exact)**2*dxbar, mesh))\n",
    "print(\"p_l2error = \", p_l2error)\n",
    "print(\"p_inner_l2error = \", p_inner_l2error)\n",
    "print(\"u_l2error = \", u_l2error)\n",
    "print(\"u_l2error_bar = \", u_l2error_bar)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73c87b71",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We observe several (not unexpected) effects:\n",
    "\n",
    "1. The error for $u$ is good on the physical domain, but much worse on the active mesh\n",
    "\n",
    "2. The error for $p$ is not even good on the physical domain, but only on the uncut elements\n",
    "\n",
    "The first effect can be cured when applying a ghost penalty stabilization which ensures a smooth extension of the discrete solution onto the whole active mesh. \n",
    "\n",
    "The second effect comes from the inconsistency that we introduced when changing the integration domain from $\\Omega$ to the active mesh. The latter problem can be cured by reconstructing a better, consistent approximation of $p$ based on good interior values and a good approximation of $u$ by postprocessing, which we discuss next.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b169ca24",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Patchwise Post-processing\n",
    "\n",
    "In the following post-processing we exploit that $u$ is (optimally) accurate on $\\Omega$ (independent of GP choice) and that $p$ is accurate in the interior. We then apply standard post-processing strategies on uncut elements and use an adapted patch-wise version for cut elements."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91f4e457",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "On each patch $\\omega$ with subtriangulation $\\mathcal{T}_\\omega$, we want to find $p_h^\\ast \\in \\mathbb{P}^{k+1}(\\mathcal{T}_\\omega)$ so that \n",
    "\\begin{align}\n",
    "(\\nabla p_h^\\ast, \\nabla q_h^\\ast)_{\\Omega \\cap \\omega} + j_h^\\omega(p_h^\\ast,q_h^\\ast) &= (u_h,\\nabla q_h^\\ast)_{\\Omega \\cap \\omega} \\qquad \\forall q_h^\\ast \\in \\mathbb{P}^{k+1}(\\mathcal{T}_\\omega) \\setminus \\mathbb{R} \\tag{P1-local} \\\\\n",
    "(p_h^\\ast,1)_{\\omega \\cap \\Omega^{\\text{int}}} &= (\\bar p_h,1)_{\\omega \\cap \\Omega^{\\text{int}}}. \\tag{P2-local}\n",
    "\\end{align}\n",
    "Here $j_h^{\\omega}(\\cdot,\\cdot)$ is a patch-local ghost penalty with an $H^1$ (instead of $L^2$) scaling.\n",
    "\n",
    "We can reformulate this as a global problem on the whole mesh, especially the integral constraint, with the help of a Lagrange multiplier $\\lambda \\in \\mathbb{P}^0(\\mathcal{T}_h^{\\text{int}})$: Find $p_h^s \\in \\mathbb{P}^{k+1}(\\mathcal{T}_h)$ and $\\lambda \\in \\mathbb{P}^0(\\mathcal{T}_h^{\\text{int}})$ so that\n",
    "\\begin{align}\n",
    "(\\nabla p_h^\\ast, \\nabla q_h^\\ast)_{\\Omega} + j_h^\\omega(p_h^\\ast,q_h^\\ast) + (q_h^\\ast,\\lambda)_{\\mathcal{T}_h^{\\text{int}}} &= (u_h,\\nabla q_h^\\ast)_{\\Omega} \\qquad \\forall q_h^\\ast \\in \\mathbb{P}^{k+1}(\\mathcal{T}_h) \\tag{P1} \\\\\n",
    "(p_h^\\ast,\\mu)_{\\mathcal{T}_h^{\\text{int}}} &= (p_h^\\ast,\\mu)_{\\mathcal{T}_h^{\\text{int}}} \\qquad \\forall \\mu \\in  \\mathbb{P}^0(\\mathcal{T}_h^{\\text{int}}). \\tag{P2}\n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5c70d63",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's start with the corresponding definition of the spaces (and corresponding trial/test functions and a solution GridFunction):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6f7b35f",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#FE Spaces for el-wise post-processing\n",
    "Vhpbase = L2(mesh, order=order+1, dgjumps=True, dirichlet=[]) # for p\n",
    "Lhbase = L2(mesh, order=0, dgjumps=True, dirichlet=[])        # for integral constraints\n",
    "Vhp = Restrict(Vhpbase, hasneg)\n",
    "Lh = Restrict(Lhbase, ci.GetElementsOfType(NEG)) # one constraint per patch / root element\n",
    "Zh = Vhp * Lh\n",
    "(ps,lam),(vs,mu) = Zh.TnT()\n",
    "\n",
    "gfz = GridFunction(Zh)\n",
    "gfps, gflam = gfz.components"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ab6205d",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Then, we setup the formulation $(P1)$-$(P2)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bafd9d5",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#Integration domain for ghost penalty\n",
    "dP = dFacetPatch(definedonelements=EA.patch_interior_facets, deformation=deformation)\n",
    "\n",
    "lhs_integrals = [grad(ps) * grad(vs) * dX,\n",
    "                 (lam * vs + mu * ps) * dxinner,\n",
    "                 1/h**2 * (ps - ps.Other()) * (vs - vs.Other()) * dP]\n",
    "rhs_integrals = [gfu * grad(vs) * dX,\n",
    "                 gfpT * mu * dxinner]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1cddd91",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now, we could put this into a global Matrix and vector, as usual, or we exploit the a-priori knowledge that the corresponding linear system decouples into the patches. For this we can use the `PatchwiseSolve`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0475899b",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "help(PatchwiseSolve)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca2e8f6f",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "gfz.vec.data = PatchwiseSolve(EA,Zh,sum(lhs_integrals),sum(rhs_integrals))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa4e72fc",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's take a look at the postprocessed solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cda0b04",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "DrawDC(lsetp1,gfps,-1,mesh,\"p\",deformation=deformation, min=-1, max=1)        # p on physical domain"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "895c1042",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$p_h^\\ast$ looks much better now. Let's also take a look at the error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6554b8d",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "ps_l2error = sqrt(Integrate((gfps - p_exact)**2*dX.order(2*order+3), mesh))\n",
    "print(\"ps_l2error = \", ps_l2error)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62a6166f",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "And also the error is much improved by the post-processing. This concludes the main part of this tutorial.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9eb95b83",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What's left?\n",
    "There are several follow-ups. Below, in the appendix we discuss some variants:\n",
    "\n",
    "1. How to deal with the r.h.s. if it is not provided on the active mesh, but only on $\\Omega$\n",
    "\n",
    "2. An alternative post-processing that work element-by-element in case that a ghost penalty stabilization has been used before\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f2e923f",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Several other tasks are left for **you** to play around with:\n",
    "\n",
    " * What is the convergence order for the considered fields in the discretization?\n",
    "\n",
    " * How can you implement Neumann boundary conditions (with or without polluting the mass balance)\n",
    "\n",
    " * If no GP stabilization is applied the coupling structure is the same as for body-fitted mixed methods. Hence, you can also hybridize the system yielding a symmetric positive definite system for the hybrid variables after static condensation. Try it out! \n",
    "\n",
    " Further explanations on the method, the above special cases and a demo file are provided through [our preprint(arXiv:2306.12722)](https://arxiv.org/abs/2306.12722). Have a look!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5e7ed25",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Appendix: extensions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "306436f8",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Extension 1: Patchwise approximation of $f$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "165bcf11",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Instead of assuming $f$ to be given on the whole active mesh, we can also use the ghost penalty stabilization to construct an approximation $f_h$ of $f$ by a discrete extension. For $\\gamma_f > 0$, we find $f_h \\in \\mathbb{P}^{k_f}(\\mathcal{T}_h)$, $k_f \\in \\mathbb{N}_0$, by solving \n",
    "\\begin{equation}\n",
    "(f_h,q_h)_{\\Omega} + \\gamma_f j_h (f_h,q_h) = (f,q_h)_{\\Omega} \\qquad \\forall q_h \\in \\mathbb{P}^{k_f}(\\mathcal{T}_h).\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59deceb9",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This problem can be solved patch-by-patch if $j_h(\\cdot,\\cdot)$ decouples across patches. The following code implements this approximation using patch-wise solves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7aea8d87",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "kf = order\n",
    "gfh = GridFunction(L2(mesh,order=kf))\n",
    "fh, wh = gfh.space.TnT()\n",
    "lhs_integrals = fh * wh * dX + 0.01*(fh-fh.Other())*(wh-wh.Other()) * dFacetPatch(deformation=deformation)\n",
    "rhs_integrals = coeff_f * wh * dX\n",
    "gfh.vec.data = PatchwiseSolve(EA,gfh.space,lhs_integrals,rhs_integrals)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ea9edab",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Extension 2: Element-wise Post-processing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f619fc9a",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If $u$ is well approximated not only on $\\Omega$, but on the whole active mesh, i.e. **if GP stabilization has been used**, a simpler element-by-element post-processing scheme is possible.\n",
    "\n",
    "By making use of the relation $u=\\nabla p$, we can repair the inconsistencies of the approximation of $p$ and achieve higher order convergence. <br>\n",
    "For each element $T \\in \\mathcal{T}_h$, we want to find $p_h^\\ast \\in \\mathcal{P}^{k+1}(T)$ such that \n",
    "\\begin{alignat}{2}\n",
    "(\\nabla p_h^\\ast, \\nabla q_h^\\ast)_T &= (u_h, \\nabla q_h^\\ast)_T \\qquad &&\\forall q_h^\\ast \\in \\mathcal{P}^{k+1}(T)\\\\\n",
    "(p_h^\\ast,1)_T &= (\\bar p_h,1)_T \\qquad &&\\text{if } T \\in \\mathcal{T}_h^{\\text{int}}\\\\\n",
    "(p_h^\\ast,1)_{T \\cap \\partial \\Omega} &= (p_D,1)_{T \\cap \\partial \\Omega} \\qquad &&\\text{if } T \\in \\mathcal{T}_h^{\\text{cut}}\\\\\n",
    "\\end{alignat}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1658b38c",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "<div align='center'>\n",
    "    <img src=\"graphics/elpp.png\" width=300px height=300px style=\"border:none\">\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "779f45fd",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#FE Spaces for el-wise post-processing\n",
    "Vhpbase = L2(mesh, order=order+1, dgjumps=False, dirichlet=[]) # order+1\n",
    "Lhbase = L2(mesh, order=0, dgjumps=False, dirichlet=[])\n",
    "\n",
    "Vhp = Restrict(Vhpbase, hasneg)\n",
    "Lh = Restrict(Lhbase, hasneg)\n",
    "Zh = Vhp * Lh\n",
    "\n",
    "gflh = GridFunction(Lh)\n",
    "\n",
    "gfz = GridFunction(Zh)\n",
    "gfps, gflam = gfz.components\n",
    "\n",
    "#Test- & Trialfunction\n",
    "(ps,lam),(vs,mu) = Zh.TnT()\n",
    "\n",
    "#Bilinear Form\n",
    "p = RestrictedBilinearForm(Zh, symmetric=False)\n",
    "p += grad(ps) * grad(vs) * dxbar\n",
    "p += (lam * vs + mu * ps) * dxinner\n",
    "p += (lam * vs + mu * ps) * ds\n",
    "\n",
    "# R.h.s. term:\n",
    "pf = LinearForm(Zh)\n",
    "pf += gfu * grad(vs) * dxbar\n",
    "\n",
    "pf += p_exact * mu * ds\n",
    "pf += gfpT * mu * dxinner\n",
    "\n",
    "#Assembly\n",
    "p.Assemble()\n",
    "pf.Assemble()\n",
    "\n",
    "#Solving the system\n",
    "gfz.vec.data = p.mat.Inverse(Zh.FreeDofs(),inverse=\"umfpack\") * pf.vec\n",
    "\n",
    "#For visualization:\n",
    "gflh.Set(p_exact, definedonelements=ci.GetElementsOfType(HASNEG))\n",
    "gflam.Set(gfps, definedonelements=ci.GetElementsOfType(HASNEG))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "729ef544",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "DrawDC(lsetp1,gfps,-1,mesh,\"p\",deformation=deformation, min=-1, max=1)        # p on physical domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eceef521",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ps_l2error = sqrt(Integrate((gfps - p_exact)**2*dX.order(2*order+3), mesh))\n",
    "print(\"ps_l2error = \", ps_l2error)"
   ]
  }
 ],
 "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"
  },
  "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": 5
}