cutfem_uc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2b3b8667",
   "metadata": {},
   "source": [
    "# Unique continuation for an elliptic interface problem \n",
    "\n",
    "In this notebook we revisit the elliptic interface problem treated in [cutfem.ipynb](cutfem.ipynb) in a slightly more challenging setting. Let us first recall the basic PDE problem: \n",
    " \n",
    " $$\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "          \\mathcal{L} u = & \\, f \n",
    "          & & \\text{in}~~ \\Omega_{\\pm}, \n",
    "          \\\\\n",
    "          [\\![u]\\!] = & \\, 0 \n",
    "          &  & \\text{on}~~ \\Gamma, \n",
    "          \\\\\n",
    "          [\\![-\\mu \\nabla u \\cdot \\mathbf{n}]\\!]   = & \\, 0 \n",
    "          &  & \\text{on}~~ \\Gamma,\n",
    "        \\end{aligned} \\right.\n",
    "$$\n",
    " \n",
    "for $\\mathcal{L}u_{\\pm} := -\\nabla \\cdot (\\mu_{\\pm} \\nabla u_{\\pm})$ defined in subdomains $\\Omega_{\\pm}$. If you are eagle eyed, you may have noticed that we omitted the boundary conditions on $\\partial \\Omega$ which in [cutfem.ipynb](cutfem.ipynb) were given by $u =  u_D$ for some given function $u_D$ on $\\partial \\Omega$. This is not a mistake. We will indeed assume here that boundary data is unknown which renders this problem ill-posed. \n",
    "\n",
    "<div>\n",
    "<img src=\"graphics/uc-overview-notebook.png\" width=\"600\"/>\n",
    "</div>\n",
    "\n",
    "To recover uniqueness and a limited form of stability we instead assume that measurements $q=u|_{\\omega}$ in a subset $\\omega \\subset \\Omega_{-}$ are available. Our objective is then to continue / extend the solution to a larger subset $\\omega \\subset B $ across the interface by solving a PDE constrained optimization problem. To this end, the basic CutFEM discretization from [cutfem.ipynb](cutfem.ipynb) has to be complemented by suitable regularization terms which are required to deal with ill-posed problems of this form. While we focus here on implementational aspects, we refer to the preprint [arXiv:2307.05210](https://arxiv.org/pdf/2307.05210.pdf) for mathematical details and error analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78c714ff",
   "metadata": {},
   "source": [
    "## The usual imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4efe3ca4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.geom2d import CSG2d, Rectangle\n",
    "from ngsolve import *\n",
    "from xfem import *\n",
    "from xfem.lsetcurv import LevelSetMeshAdaptation\n",
    "from math import pi\n",
    "from netgen.csg import *\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50a043c1",
   "metadata": {},
   "source": [
    "## Mesh \n",
    "For simplicity we create a mesh in which the subdomains $\\omega= [-0.5,0.5]^2$ and $B=[-1.25,1.25]^2$ are fitted. This is easy to achieve using constructive solid geometry definitions, see [CSG in 2D](https://docu.ngsolve.org/latest/i-tutorials/unit-4.1.2-csg2d/csg2d.html). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09aa23d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = CSG2d()\n",
    "omega_dom = Rectangle( pmin=(-0.5,-0.5), pmax=(0.5,0.5), mat=\"omega\", bc=\"bc_omega\")\n",
    "B_dom = Rectangle( pmin=(-1.25,-1.25), pmax=(1.25,1.25), mat=\"B\", bc=\"bc_B\")\n",
    "full_dom = Rectangle( pmin=(-1.5,-1.5), pmax=(1.5,1.5), mat=\"full\", bc=\"bc_Omega\")\n",
    "only_B = B_dom - omega_dom\n",
    "only_B.Mat(\"only_B\")\n",
    "rest = full_dom - B_dom\n",
    "geo.Add(omega_dom)\n",
    "geo.Add(only_B)\n",
    "geo.Add(rest)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.125))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50b833b1",
   "metadata": {},
   "source": [
    "Here, `mesh.Materials(\"omega\")` correspond to $\\omega$ and `mesh.Materials(\"only_B|omega\")` to the target domain $B$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6f703d9",
   "metadata": {},
   "source": [
    "## Levelset \n",
    "The unit ball in the $4$-norm will be used as the levelset function. Note that the subdomains $\\omega$ and $B$ are clearly visible in the webgui below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69efaf31",
   "metadata": {},
   "outputs": [],
   "source": [
    "r44 = x**4 + y**4\n",
    "r41 = sqrt(sqrt(r44))\n",
    "levelset = r41 - 1.0\n",
    "DrawDC(levelset, -3.5, 2.5, mesh,\"levelset\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47049c44",
   "metadata": {},
   "source": [
    "## Sample solution\n",
    "We chose the diffusivities as $\\mu_{-} = 2, \\mu_{+} = 20$ and calculate the right hand side so that `solution` solves the PDE problem (it is easy to check that the interface conditions are fulfilled as well)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e72da98",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = (2,20)\n",
    "solution = [(1/sqrt(2.0))*(1+pi*mu[0]/mu[1]) - cos(pi / 4 * r44), (mu[0]/mu[1])*(pi/sqrt(2.0)) * r41]\n",
    "coef_f = [-mu[i] * (solution[i].Diff(x).Diff(x) + solution[i].Diff(y).Diff(y) ) for i in range(2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "024ae33f",
   "metadata": {},
   "source": [
    "## Variational formulation and discretization space\n",
    "\n",
    "The variational formulation to approximate the solution $u$ is given as follows:\n",
    "Find $(u_h,z_h) \\in V_h^\\Gamma \\times V_h^0$ such that \n",
    "$$\n",
    "A_h(v_h,z_h) + s_h(u_h,v_h) + (u_h-q,v_h)_{ \\omega } + A_h(u_h,w_h) - s_h^{\\ast}(z_h,w_h) = \\ell(v_h,w_h)  \\qquad (*)\n",
    "$$ \n",
    "for all $(v_h,w_h) \\in V_h^\\Gamma \\times V_h^0$. Here, \n",
    "\n",
    "* $u_h$ is a discretization of $u$ sought in the CutFEM space $V_h^\\Gamma$ introduced in [cutfem.ipynb](cutfem.ipynb),\n",
    "* $z_h$ is a Lagrange multiplier which lives in a standard finite element space $V_h^0$ with homogeneous boundary conditions on $\\partial \\Omega$. \n",
    "\n",
    "Before we go into more detail on the terms in the variational formulation, let us first set up our discretization space $W = V_h^\\Gamma \\times V_h^0$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc45a9a2",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddbdd002",
   "metadata": {},
   "outputs": [],
   "source": [
    "order = 2 \n",
    "lsetadap = LevelSetMeshAdaptation(mesh, order=order, levelset=levelset)\n",
    "lsetp1 = lsetadap.lset_p1\n",
    "ci = CutInfo(mesh, lsetp1)\n",
    "Vh = H1(mesh, order=order, dirichlet=[],dgjumps=True)\n",
    "Vh0 = H1(mesh, order=order, dirichlet=\"bc_Omega\",dgjumps=False)\n",
    "hasneg = ci.GetElementsOfType(HASNEG)\n",
    "haspos = ci.GetElementsOfType(HASPOS)\n",
    "hasif = ci.GetElementsOfType(IF)\n",
    "Vh_Gamma = Compress(Vh, GetDofsOfElements(Vh, hasneg)) \\\n",
    "              * Compress(Vh, GetDofsOfElements(Vh, haspos)) \\\n",
    "              * Vh0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11b7fd19",
   "metadata": {},
   "source": [
    "We obtain test and trial functions in the usual way. Note that the dual variable is in fact univalued on $\\Omega$ but we prefer to define `z = [z0,z0]` for implementational / notational ease. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38a88eaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "u1,u2,z0 =  Vh_Gamma.TrialFunction()\n",
    "v1,v2,w0 =  Vh_Gamma.TestFunction()\n",
    "u = [u1,u2]\n",
    "z = [z0,z0]\n",
    "v = [v1,v2]\n",
    "w = [w0,w0]\n",
    "gradu, gradz, gradv, gradw = [[grad(fun[i]) for i in [0, 1]] for fun in [u, z, v, w]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5f1827c",
   "metadata": {},
   "source": [
    "## Unfitted Nitsche formulation\n",
    "\n",
    "Let us now elaborate on the terms in the variational formulation (*). The terms $A_h(u_h,w_h)$ and $A_h(v_h,z_h) $  simply correspond to the standard unfitted Nitsche formulation. Since the dual variable is continuous only the adjoint consistency term remains on the interface, i.e. we have\n",
    "\n",
    "$$\n",
    "A_h(u,w) = \\sum\\limits_{\\pm} (\\mu_{\\pm} \\nabla u_{\\pm}, \\nabla w_{\\pm})_{ \\Omega_{{\\pm},h} } + \\int\\limits_{\\Gamma_h}  \\{\\!\\!\\{ -\\mu \\nabla w  \\} \\!\\!\\}  \\cdot \\mathbf{n} [\\![u]\\!].\n",
    "$$\n",
    "\n",
    "We also define the corresponding right hand side. The implementation given below is well-known from [cutfem.ipynb](cutfem.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2890dcdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "dC = tuple([dCut(lsetp1, dt, deformation=lsetadap.deform,\n",
    "                     definedonelements=ci.GetElementsOfType(HAS(dt)))\n",
    "                for dt in [NEG, POS]])\n",
    "ds = dCut(lsetp1, IF, deformation=lsetadap.deform)\n",
    "\n",
    "n = 1.0 / grad(lsetp1).Norm() * grad(lsetp1)\n",
    "kappa = ( mu[1]/sum(mu) , mu[0]/sum(mu) )\n",
    "average_flux_z = sum([- kappa[i] * mu[i] * gradz[i] * n for i in [0, 1]])\n",
    "average_flux_w = sum([- kappa[i] * mu[i] * gradw[i] * n for i in [0, 1]])\n",
    "\n",
    "# a(uh,wh)\n",
    "a = BilinearForm(Vh_Gamma, symmetric=True)\n",
    "a += sum(mu[i] * gradu[i] * gradw[i] * dC[i] for i in [0, 1])\n",
    "a +=  average_flux_w * (u[0] - u[1]) * ds\n",
    "\n",
    "# a(vh,zh)\n",
    "a += sum(mu[i] * gradv[i] * gradz[i] * dC[i] for i in [0, 1])\n",
    "a +=  average_flux_z * (v[0] - v[1]) * ds\n",
    "\n",
    "f = LinearForm(Vh_Gamma)\n",
    "f += sum(coef_f[i] * w[i] * dC[i] for i in [0, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4397a3e2",
   "metadata": {},
   "source": [
    "## Data constraint \n",
    "\n",
    "Next we add the data constraint  $(u_h-q,v_h)_{ \\omega }$ to the bilinear form and the right hand side. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03e1d122",
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma_data = 1e5\n",
    "a += sum( gamma_data * u[i] * v[i] * dCut(lsetp1, dt,definedon=mesh.Materials(\"omega\"), deformation=lsetadap.deform)\n",
    "            for i,dt in zip([0, 1],[NEG,POS]) )\n",
    "f += sum(gamma_data * solution[i] * v[i] * dCut(lsetp1, dt,definedon=mesh.Materials(\"omega\"), deformation=lsetadap.deform)\n",
    "            for i,dt in zip([0, 1],[NEG,POS]) )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10ec31b0",
   "metadata": {},
   "source": [
    "## Primal stabilization terms\n",
    "The main task is to definine the terms of the primal stabilization $s_h(u_h,v_h)$. These terms can be interpreted as regularizers which are used to incorporate a priori knowledge of the solution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f048e3f",
   "metadata": {},
   "source": [
    "## Continuous interior penalty \n",
    "The first term is the well-known continuous interior penalty (CIP)\n",
    "\n",
    "$$\n",
    "J_{ \\mathrm{CIP} }(u,v) := \\sum\\limits_{ \\pm } \\sum\\limits_{ F \\in  \\mathcal{F}^{\\pm} } h \\!\\! \\int\\limits_{ F } \\!\\! \\mu_{\\pm} [\\![ \\nabla u_{\\pm} \\cdot \\mathbf{n} ]\\!] [\\![ \\nabla v_{\\pm} \\cdot \\mathbf{n} ]\\!]    \n",
    "$$\n",
    "* The facet sets  $\\mathcal{F}^{\\pm}$ contain all facets lying between elements of the active mesh corresponding to subdomain $\\Omega_{\\pm}$. In the sketch given below the facets in  $\\mathcal{F}^{-}$ are indicated by red dashed lines. We obtain them below using `GetFacetsWithNeighborTypes(mesh, a=hasneg, b=hasneg)`. \n",
    "* Obviously, there is a close connection between the 'derivative-jump' version of the ghost-penalty stabilization often employed in unfitted methods and $J_{ \\mathrm{CIP} }$. Whereas the former is usually limited to a band of facets around the interface, we employ the latter here globally. On the other hand, it is not necessary to include higher order jumps terms into $J_{ \\mathrm{CIP} }$ as is usually done for the derivative-jump' version of the ghost-penalty when `order` is larger one. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ec96e15",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"graphics/Stab-sketch-notebook.png\" width=\"600\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e343900",
   "metadata": {},
   "outputs": [],
   "source": [
    "ba_facets = { NEG: GetFacetsWithNeighborTypes(mesh, a=hasneg, b=hasneg),\n",
    "              POS: GetFacetsWithNeighborTypes(mesh, a=haspos, b=haspos) }\n",
    "dk = tuple([ dCut(lsetp1, dt, skeleton=True, definedonelements=ba_facets[dt], deformation=lsetadap.deform)\n",
    "               for dt in [NEG,POS] ])\n",
    "nF = specialcf.normal(mesh.dim)\n",
    "h = specialcf.mesh_size\n",
    "\n",
    "gamma_CIP = 5e-2\n",
    "a += sum( [ gamma_CIP * h * mu[i] * InnerProduct( (gradu[i] - gradu[i].Other()) * nF, \n",
    "             (gradv[i] - gradv[i].Other()) * nF ) * dk[i] for i in [0,1] ]  )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c47e1c66",
   "metadata": {},
   "source": [
    "## Galerkin-least squares \n",
    "Additionally, we add a Galerkin-least-squares stabilization which incorporates the PDE constraint element-wise: \n",
    "\n",
    "$$ \n",
    "J_{ \\mathrm{GLS} }(u,v) := \\sum\\limits_{\\pm} \\sum\\limits_{ T \\in \\mathcal{T}^{\\pm} } h^2 ( \\mathcal{L} u_{\\pm} , \\mathcal{L} v_{\\pm} )_{   T_{\\pm}  }\n",
    "$$\n",
    "\n",
    "Here, $T^{-} := T \\cap \\Omega_{-,h} $ is the part of the element lying in domain $\\Omega_{-}$ (as indicated in the sketch above in cyan). The differential operator $\\mathcal{L}$ is implemented below in the function `calL`.\n",
    "Since the PDE is inhomogeneous, we also have to add a corresponding term \n",
    "\n",
    "$$\n",
    "\\sum\\limits_{\\pm} \\sum\\limits_{ T \\in \\mathcal{T}^{\\pm} } h^2 ( f_{\\pm} , \\mathcal{L} v_{\\pm} )_{   T_{\\pm}  }\n",
    "$$\n",
    "\n",
    "to the right hand side to preserve consistency. Note that this stabilization is only relevant for `order` larger one since piecewise linear functions are in the kernel of $\\mathcal{L}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e823c76",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calL(fun):\n",
    "    hesse = [fun[i].Operator(\"hesse\") for i in [0,1]]\n",
    "    return (-mu[0]*hesse[0][0,0]-mu[0]*hesse[0][1,1],-mu[1]*hesse[1][0,0]-mu[1]*hesse[1][1,1])\n",
    "gamma_GLS = 5e-2\n",
    "a += sum( [ gamma_GLS * h**2 * calL(u)[i] * calL(v)[i] * dC[i] for i in [0, 1] ] )\n",
    "f += sum( [ gamma_GLS * h**2 * coef_f[i] * calL(v)[i] * dC[i] for i in [0, 1] ] )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a63e95db",
   "metadata": {},
   "source": [
    "## Stabilization on  the interface \n",
    "We use the following stabilization terms on the interfae:\n",
    "\n",
    "$$\n",
    "J_{ \\mathrm{H} }^{ \\Gamma_h }(u,v) := \\frac{ \\bar{\\mu} }{h}  \\int\\limits_{ \\Gamma_h } [\\![u]\\!] [\\![v]\\!], \\qquad \n",
    "J_{ \\mathrm{N} }^{ \\Gamma_h }(u,v) := h \\int\\limits_{ \\Gamma_h } [\\![\\mu \\nabla u \\cdot \\mathbf{n}]\\!] [\\![\\mu \\nabla v \\cdot \\mathbf{n}]\\!],   \\qquad \n",
    "J_{ \\mathrm{T} }^{ \\Gamma_h }(u,v) := h \\bar{\\mu} \\int\\limits_{ \\Gamma_h } [\\![ \\nabla_{\\Gamma_h } u ]\\!] [\\![ \\nabla_{\\Gamma_h } v  ]\\!]\n",
    "$$\n",
    "Here,\n",
    "* $J_{ \\mathrm{H} }^{ \\Gamma_h }$ is simply the stability term in the usual unfitted Nitsche formulation. \n",
    "* $J_{ \\mathrm{N} }^{ \\Gamma_h }(u,v)$ penalizes the normal jumps across the interface \n",
    "* and $J_{ \\mathrm{T} }^{ \\Gamma_h }$ the tangential jump on the interface. Here, we defined the projector $P_{\\Gamma_h} := I - \\mathbf{n}_{\\Gamma_h} (\\mathbf{n}_{\\Gamma_h})^T$ and the surface gradient $\\nabla_{\\Gamma_h} := P_{\\Gamma_h} \\nabla$. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24f1892d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def P(fun):\n",
    "    return fun - (fun * n) * n\n",
    "\n",
    "jump_flux_u =  (mu[0] * gradu[0] - mu[1] * gradu[1]) * n\n",
    "jump_flux_v =  (mu[0] * gradv[0] - mu[1] * gradv[1]) * n\n",
    "jump_tangential_u =  P(gradu[0]) - P(gradu[1])\n",
    "jump_tangential_v =  P(gradv[0]) - P(gradv[1])\n",
    "\n",
    "gamma_IF = 1e-3\n",
    "mubar = (mu[0]+mu[1])/2\n",
    "a += gamma_IF * (mubar/h) * (u[0] - u[1]) * (v[0] - v[1])  * ds\n",
    "a += gamma_IF * h * jump_flux_u * jump_flux_v * ds\n",
    "a += gamma_IF * mubar * h * jump_tangential_u * jump_tangential_v * ds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e6c52dd",
   "metadata": {},
   "source": [
    "## (Weak) Tikhonov stabilization \n",
    "Whereas the previous stabilization terms are consistent for sufficiently smooth solutions of the PDE, we will now add a weakly inconsistent penalty. By this we mean that the resulting perturbation scales with $h$ and the polynomial order `k` in a way that does not spoil the error estimates (see the paper for details). \n",
    "\n",
    "$$\n",
    "J_{ \\alpha }(u,v) := h^{2k} \\sum\\limits_{\\pm}  \\alpha_1 ( u_{\\pm}, v_{\\pm})_{ \\Omega_{\\pm,h}^{\\dagger} } + \\alpha_2 (\\nabla u_{\\pm}, \\nabla v_{\\pm} )_{ \\Omega_{\\pm,h}^{\\dagger} }.\n",
    "$$\n",
    "\n",
    "Here, $\\Omega_{i,h}^{\\dagger}$ encompasses all elements of the active mesh in the interior subdomain obtained by `ci.GetElementsOfType(HASNEG)`  (the union of the cyan and yellow region in the sketch)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c3e08e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = [1e-5,1e-2]\n",
    "dGeom = tuple([ dx(definedonelements=ci.GetElementsOfType(dt))\n",
    "                     for dt in [HASNEG,HASPOS]])\n",
    "a += sum( alpha[0] * h**(2*order) * u[i] * v[i] * dGeom[i] for i in [0, 1])\n",
    "a += sum( alpha[1] * h**(2*order) * grad( u[i]) * grad(v[i]) * dGeom[i] for i in [0, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1712890e",
   "metadata": {},
   "source": [
    "## Dual stabilization \n",
    "For the dual stabilization we make the simple choice \n",
    "$$\n",
    "s^{\\ast}_h(z_h,w_h) := (\\mu \\nabla z_h, \\nabla w_h)_{ \\Omega }.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e77de325",
   "metadata": {},
   "outputs": [],
   "source": [
    "a += sum(-mu[i] * gradz[i] * gradw[i] * dC[i] for i in [0, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8f2b4f9",
   "metadata": {},
   "source": [
    "## Solve linear system and measure the error \n",
    "\n",
    "Let us solve the linear system, measure the $L^2$-error in the target domain and plot the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8eb67295",
   "metadata": {},
   "outputs": [],
   "source": [
    "a.Assemble()\n",
    "f.Assemble()\n",
    "gfu = GridFunction(Vh_Gamma)\n",
    "gfu.vec.data = a.mat.Inverse(Vh_Gamma.FreeDofs(),inverse= \"sparsecholesky\"  )* f.vec\n",
    "err_sqr = sum( [ (gfu.components[i]   - solution[i])**2 * dCut(lsetp1, dt,definedon=mesh.Materials(\"only_B|omega\"),\n",
    "                deformation=lsetadap.deform, order=2*order)   for i,dt in zip([0, 1],[NEG,POS] ) ] )\n",
    "l2_err = sqrt(Integrate(err_sqr, mesh))\n",
    "print(\"L2-error = \", l2_err)\n",
    "\n",
    "uh = IfPos(lsetp1, gfu.components[1], gfu.components[0])\n",
    "deform_graph = CoefficientFunction((lsetadap.deform[0], lsetadap.deform[1], 4*uh))\n",
    "DrawDC(lsetp1, gfu.components[0], gfu.components[1], mesh, \"graph_of_u\", deformation=deform_graph, min=0, max=0.25)\n",
    "DrawDC(lsetp1, (gfu.components[0]-solution[0])**2, (gfu.components[1]-solution[1])**2, mesh, \"error\",min=1e-8, max=1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39d7b009",
   "metadata": {},
   "source": [
    "Notice that the error is concentrated near the boundary $\\partial \\Omega$ where data is unknown. Actually, we have to require that $B \\setminus \\omega$ does not touch $\\partial \\Omega$ to obtain practically useful error estimates (see our paper for more details)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2ccf888",
   "metadata": {},
   "source": [
    "## Now it's your turn. \n",
    "Here are some suggestions for further experiments: \n",
    "* Play around with the contrast `mu`. \n",
    "* In practice one usually only has access to noisy data. Try to add some noise to the data `f` and `q` and check how the method performs.\n",
    "* Is it necessary that the discretization spaces for primal and dual variable have the same order? Maybe it also suffices to have `Vh0 = H1(mesh, order=1, dirichlet=\"bc_Omega\",dgjumps=False)`?\n",
    "* What happens if the mesh deformation is disabled?\n"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}