tracefem_vector.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Vector SurfacePDEs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Vector-Laplace problem\n",
    "We want to solve a problem of the form\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "          -2 \\, \\textbf{P} \\, \\textbf{div}_\\Gamma (E_s(\\textbf{u})) + \\textbf{u} = & \\, \\textbf{f} \n",
    "          & & \\text{on}~~ \\Gamma,          \n",
    "        \\end{aligned}\n",
    "$$\n",
    "\n",
    "with $\\textbf{P} \\textbf{u} = \\textbf{u}$ and $\\Gamma := \\{ \\phi = 0 \\}$. Here, $\\textbf{P} = \\textbf{I} - \\textbf{n}\\textbf{n}^T$ is the projection onto the tangential space and \n",
    "\n",
    "\\begin{equation*}\n",
    "\\nabla_{\\Gamma} \\textbf{u} = \\textbf{P} \\nabla \\textbf{u} \\textbf{P}, \\qquad E_s(\\textbf{u}) = \\frac{1}{2} (\\nabla_{\\Gamma} \\textbf{u} + \\nabla_{\\Gamma}^T \\textbf{u})\n",
    "\\end{equation*} \n",
    "are the surface gradient for vector-valued functions and the surface rate-of-strain tensor. Moreover, the surface divergence operators are defined as follows:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\operatorname{div}_{\\Gamma}\\textbf{u} = \\operatorname{tr}(\\nabla_{\\Gamma}\\textbf{u}), \\qquad \\textbf{div}_{\\Gamma}\\textbf{A} = (\\operatorname{div}_{\\Gamma} (\\textbf{e}_1^T \\textbf{A}), \\ldots, \\operatorname{div}_{\\Gamma}(\\textbf{e}_n^T \\textbf{A}))^T.\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# ngsolve stuff\n",
    "from ngsolve import *\n",
    "# basic xfem functionality\n",
    "from xfem import *\n",
    "# basic geometry features (for the background mesh)\n",
    "from netgen.occ import *\n",
    "# for isoparametric mapping\n",
    "from xfem.lsetcurv import *\n",
    "# pi\n",
    "from math import pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Geometry approximation\n",
    "First, 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": [
    "cube = OCCGeometry(unit_cube.shape.Scale((0,0,0),3).Move((-1.5,-1.5,-1.5)), dim=3)\n",
    "mesh = Mesh (cube.GenerateMesh(maxh=0.4, quad_dominated=False))\n",
    "Draw(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "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)\n",
    "\n",
    "To compute the corresponding mapping we use the `LevelSetMeshAdaptation` class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "phi = Norm(CoefficientFunction((x,y,z)))-1\n",
    "lsetmeshadap = LevelSetMeshAdaptation(mesh, order=2)\n",
    "deformation = lsetmeshadap.CalcDeformation(phi)\n",
    "lset_approx = lsetmeshadap.lset_p1\n",
    "mesh.SetDeformation(deformation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Weak formulation with penalty term for tangential flow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Isoparametric TraceFE / CutFE space\n",
    "For the discretization we use standard background FESpaces restricted to the tetrahedrons cut by $\\Gamma_h$ and transformed to the deformed mesh:\n",
    "\n",
    "\\begin{equation*}\n",
    "U_{h,\\Theta}^k = \\left\\{ \\textbf{u}_h \\circ \\Psi_h^{-1} \\mid \\textbf{u}_h \\in (V_h^k)^3|_{ \\Omega^{\\Gamma}} \\right\\} \\end{equation*}\n",
    "\n",
    "In NGSolve we will simply take the space $(V_h^k)^3$, which is of order $k$, but mark the irrelevant dofs using the CutInfo-class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# FESpace \n",
    "Uhbase = VectorH1(mesh, order=2)\n",
    "ci = CutInfo(mesh, lset_approx)\n",
    "\n",
    "Uh = Compress(Uhbase,GetDofsOfElements(Uhbase,ci.GetElementsOfType(IF)))\n",
    "\n",
    "# expressions of test and trial functions:\n",
    "u,v = Uh.TnT()\n",
    "# overwrite freedofs of VhG to mark only dofs that are involved in the cut problem\n",
    "freedofs = Uh.FreeDofs()\n",
    "# declare a grid function to store the solution\n",
    "gfu = GridFunction(Uh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Isoparametric TraceFEM / CutFEM discretization\n",
    "\n",
    "Find $\\textbf{u}_h \\in U_{h,\\Theta}^k$ s.t.\n",
    "    \\begin{align*}\n",
    "      \\int_{\\Gamma_h} \\mathrm{tr}(\\tilde{E}_s(\\textbf{u}_h) \\tilde{E}_s(\\textbf{v}_h)) &+ \\textbf{u}_h \\cdot \\textbf{v}_h \\,ds_h + \\eta \\int_{\\Gamma_h} (\\textbf{u}_h \\cdot \\textbf{n}_h)  (\\textbf{v}_h \\cdot \\textbf{n}_h) \\,ds_h  + \\rho\\int_{\\Omega_{\\Theta}^\\Gamma}(\\nabla \\textbf{u}_h \\textbf{n}_h)\\cdot(\\nabla\\textbf{v}_h \\textbf{n}_h)\\,dx \\\\\n",
    "      &= (\\mathbf{f}, \\textbf{v}_h)_{0,\\Gamma} \\qquad \\text{ for all }\\textbf{v}_h\\in U_{h,\\Theta}^k.\n",
    "    \\end{align*}\n",
    "with $\\tilde{E}_s(\\textbf{u}) = \\frac{1}{2} (\\nabla_{\\Gamma_h} \\textbf{u} + \\nabla_{\\Gamma_h}^T \\textbf{u})$, $\\eta \\sim h^{-2}$ and $h\\lesssim \\rho \\lesssim h^{-1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# declare the integration domains\n",
    "ds = dCut(levelset=lset_approx, domain_type=IF, definedonelements=ci.GetElementsOfType(IF),deformation=deformation)\n",
    "dX = dx(definedonelements=ci.GetElementsOfType(IF),deformation=deformation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### helper coefficients / parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# calculate normal vector n, mesh size h and parameters eta and rho\n",
    "n = 1.0/Norm(grad(lset_approx)) * grad(lset_approx)\n",
    "h = specialcf.mesh_size\n",
    "eta = 100.0/(h*h)\n",
    "rho = 1.0/h\n",
    "\n",
    "# define solution and right-hand side\n",
    "functions = {\n",
    "           \"extvsol1\" : ((-y-z)*x + y*y + z*z),\n",
    "           \"extvsol2\" : ((-x-z)*y + x*x + z*z),\n",
    "           \"extvsol3\" : ((-x-y)*z + x*x + y*y),\n",
    "           \"rhs1\" : -((y + z) * x - y * y - z * z) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),\n",
    "           \"rhs2\" : ((-x - z) * y + x * x + z * z) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),\n",
    "           \"rhs3\" : ((-x - y) * z + x * x + y * y) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),\n",
    "}\n",
    "uSol = CoefficientFunction((functions[\"extvsol1\"],functions[\"extvsol2\"],functions[\"extvsol3\"]))\n",
    "rhs = CoefficientFunction((functions[\"rhs1\"],functions[\"rhs2\"],functions[\"rhs3\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The bi- and linear forms:\n",
    "    \\begin{align*}\n",
    "      \\int_{\\Gamma_h} \\mathrm{tr}(\\tilde{E}_s(\\textbf{u}_h) \\tilde{E}_s(\\textbf{v}_h)) & + \\textbf{u}_h \\cdot \\textbf{v}_h \\,ds_h + \\eta \\int_{\\Gamma_h} (\\textbf{u}_h \\cdot \\textbf{n}_h)  (\\textbf{v}_h \\cdot \\textbf{n}_h) \\,ds_h  + \\rho\\int_{\\Omega_{\\Theta}^\\Gamma}(\\nabla \\textbf{u}_h \\textbf{n}_h)\\cdot(\\nabla\\textbf{v}_h \\textbf{n}_h)\\,dx \\\\ &= (\\mathbf{f}, \\textbf{v}_h)_{0,\\Gamma} \\qquad \\text{for all }\\textbf{v}_h\\in U_{h,\\Theta}^k.\n",
    "    \\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# define projection matrix\n",
    "Pmat = Id(3) - OuterProduct(n,n)\n",
    "\n",
    "# bilinear forms:\n",
    "a = BilinearForm(Uh, symmetric = True)\n",
    "a += (InnerProduct(Pmat*Sym(grad(u))*Pmat,Pmat*Sym(grad(v))*Pmat))*ds\n",
    "a += u * v * ds\n",
    "a += (eta*((u * n) * (v * n))) * ds\n",
    "a += rho * ((grad(u)*n) * (grad(v)*n)) * dX\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(Uh)\n",
    "f += rhs * v * ds\n",
    "f.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can now solve the problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu.vec[:] = 0.0\n",
    "c = a.mat.CreateSmoother(freedofs)\n",
    "inv = CGSolver(a.mat, c, maxsteps=300)\n",
    "gfu.vec.data = inv * f.vec\n",
    "print(\"CG with Jac. preconditioner, steps:\", inv.GetSteps())\n",
    "    \n",
    "l2error = sqrt(Integrate(InnerProduct(gfu-uSol, gfu-uSol)*ds, mesh))\n",
    "print (\"l2error : \", l2error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Visualization\n",
    "Create vtk-output file to be read by Paraview."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "vtk = VTKOutput(ma=mesh,\n",
    "                coefs=[lset_approx,deformation,gfu,uSol,rhs],\n",
    "                names=[\"P1-levelset\",\"deform\",\"uh\",\"uex\",\"rhs\"],\n",
    "                filename=\"vectorlaplacepenalty3d\",subdivision=2)\n",
    "vtk.Do()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<center>\n",
    "![title](graphics/tracefem_vector1.png)\n",
    "</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "ls -al vectorlaplacepenalty3d.vtu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "#%%bash\n",
    "#paraview vectorlaplacepenalty3d.vtk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Weak formulation with Lagrange-multiplier term for tangential flow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Isoparametric TraceFE / CutFE space\n",
    "For the discretization we use standard background FESpaces restricted to the tetrahedrons cut by $\\Gamma_h$ and transformed to the deformed mesh:\n",
    "\n",
    "\\begin{equation*}\n",
    "U_{h,\\Theta}^k = \\left\\{ \\textbf{u}_h \\circ \\Psi_h^{-1} \\mid \\textbf{u}_h \\in (V_h^k)^3|_{ \\Omega^{\\Gamma}} \\right\\}, \\qquad M_{h,\\Theta}^k = \\left\\{ \\lambda_h \\circ \\Psi_h^{-1} \\mid \\lambda_h \\in V_h^k|_{ \\Omega^{\\Gamma}} \\right\\}\n",
    "\\end{equation*}\n",
    "\n",
    "In NGSolve we will simply take the space $(V_h^k)^3$ for the solution and $V_h^k$ for the lagrange multiplier but mark the irrelevant dofs using the CutInfo-class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# FESpace \n",
    "Vhbase = VectorH1(mesh, order=2)\n",
    "Mhbase = H1(mesh, order=1)\n",
    "#ci = CutInfo(mesh, lset_approx)\n",
    "\n",
    "VhG = Compress(Vhbase,GetDofsOfElements(Vhbase,ci.GetElementsOfType(IF)))\n",
    "Mh = Compress(Mhbase,GetDofsOfElements(Mhbase,ci.GetElementsOfType(IF)))\n",
    "\n",
    "Qh = VhG * Mh\n",
    "(u,p), (v,q) = Qh.TnT()\n",
    "\n",
    "# overwrite freedofs of VhG to mark only dofs that are involved in the cut problem\n",
    "freedofs = Qh.FreeDofs()\n",
    "# declare a grid function to store the solution\n",
    "gfu = GridFunction(Qh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Isoparametric TraceFEM / CutFEM discretization \n",
    " Find $(\\textbf{u}_h,p_h) \\in U_{h,\\Theta}^k \\times M_{h,\\Theta}^l$ s.t.\n",
    "    \\begin{align*}\n",
    "      A_h(\\textbf{u}_h,\\textbf{v}_h) + B_h(\\textbf{v}_h,p_h) &= (\\mathbf{f}, \\textbf{v}_h)_{0,\\Gamma} &&\\text{for all }\\textbf{v}_h\\in U_{h,\\Theta}^k \\\\[1ex]\n",
    "      B_h(\\textbf{u}_h,q_h) &= 0 &&\\text{for all } q_h\\in M_{h,\\Theta}^l\n",
    "    \\end{align*}\n",
    "    with bilinear forms using stabilization\n",
    "    \\begin{align*}\n",
    "      A_h(\\textbf{u}_h,\\textbf{v}_h) &:= \\int_{\\Gamma_h} \\mathrm{tr}(\\tilde{E}_s(\\textbf{u}_h) \\tilde{E}_s(\\textbf{v}_h)) + \\textbf{u}_h \\cdot \\textbf{v}_h \\,ds_h + \\rho\\int_{\\Omega_{\\Theta}^\\Gamma}(\\nabla \\textbf{u}_h \\textbf{n})\\cdot(\\nabla\\textbf{v}_h \\textbf{n})\\,dx,\\\\[1ex]\n",
    "      B_h(\\textbf{u}_h,q_h) &:= (\\textbf{u}_h\\cdot\\textbf{n},q_h)_{0,\\Gamma} + \n",
    "        \\rho \\int_{\\Omega_{\\Theta}^\\Gamma}(\\textbf{n}^T\\nabla \\textbf{u}_h\\textbf{n})(\\textbf{n}\\cdot\\nabla q_h)\\,dx, \n",
    "    \\end{align*}\n",
    "    with $\\tilde{E}_s(\\textbf{u}) = \\frac{1}{2} (\\nabla_{\\Gamma_h} \\textbf{u} + \\nabla_{\\Gamma_h}^T \\textbf{u})$ and $h\\lesssim \\rho \\lesssim h^{-1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# declare the integration domains\n",
    "lset_if = {\"levelset\": lset_approx, \"domain_type\": IF, \"subdivlvl\": 0}\n",
    "def SurfaceIntegral(form):\n",
    "    return SymbolicBFI(levelset_domain = lset_if , form = form.Compile(), definedonelements=ci.GetElementsOfType(IF))\n",
    "def ActiveMeshIntegral(form):\n",
    "    return SymbolicBFI(form = form.Compile(), definedonelements=ci.GetElementsOfType(IF))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### helper coefficients / parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# calculate normal vector n, mesh size h and parameters eta and rho\n",
    "n = 1.0/sqrt(InnerProduct(grad(lset_approx),grad(lset_approx))) * grad(lset_approx)\n",
    "h = specialcf.mesh_size\n",
    "rho = 1.0*h\n",
    "\n",
    "# define solution and right-hand side\n",
    "functions = {\n",
    "           \"extvsol1\" : ((-y-z)*x + y*y + z*z),\n",
    "           \"extvsol2\" : ((-x-z)*y + x*x + z*z),\n",
    "           \"extvsol3\" : ((-x-y)*z + x*x + y*y),\n",
    "           \"rhs1\" : -((y + z) * x - y * y - z * z) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),\n",
    "           \"rhs2\" : ((-x - z) * y + x * x + z * z) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),\n",
    "           \"rhs3\" : ((-x - y) * z + x * x + y * y) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),\n",
    "}\n",
    "uSol = CoefficientFunction((functions[\"extvsol1\"],functions[\"extvsol2\"],functions[\"extvsol3\"]))\n",
    "rhs = CoefficientFunction((functions[\"rhs1\"],functions[\"rhs2\"],functions[\"rhs3\"]))\n",
    "\n",
    "# define projection matrix\n",
    "Pmat = Id(3) - OuterProduct(n,n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "bilinear form\n",
    "    \\begin{align*}\n",
    "      A_h(\\textbf{u}_h,\\textbf{v}_h) &:= \\int_{\\Gamma_h} \\mathrm{tr}(E_s(\\textbf{u}_h) E_s(\\textbf{v}_h)) + \\textbf{u}_h \\cdot \\textbf{v}_h \\,ds_h + \\rho\\int_{\\Omega_{\\Theta}^\\Gamma}(\\nabla \\textbf{u}_h \\textbf{n})\\cdot(\\nabla\\textbf{v}_h \\textbf{n})\\,dx,\n",
    "    \\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a = BilinearForm(Qh, symmetric = False)\n",
    "a += (InnerProduct(Pmat*Sym(grad(u))*Pmat,Pmat*Sym(grad(v))*Pmat))*ds\n",
    "a += u * v * ds\n",
    "a += rho* ((grad(u)*n) * (grad(v)*n)) * dX"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "bilinear form\n",
    "    \\begin{align*}\n",
    "      B_h(\\textbf{u}_h,q_h) &:= (\\textbf{u}_h\\cdot\\textbf{n},q_h)_{0,\\Gamma} + \n",
    "        \\rho \\int_{\\Omega_{\\Theta}^\\Gamma}(\\textbf{n}^T\\nabla \\textbf{u}_h\\textbf{n})(\\textbf{n}\\cdot\\nabla q_h)\\,dx, \n",
    "    \\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a += (p * (v * n) + q * (u * n)) * ds\n",
    "a += (rho* ((n*grad(p)) * (n*(grad(v)*n))+ (n*grad(q)) * (n*(grad(u)*n)))) * dX\n",
    "a.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "linear form"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = LinearForm(Qh)\n",
    "f += rhs * v * ds\n",
    "f.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can now solve the problem (recall that freedofs only marks relevant dofs):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu.vec[:] = 0.0\n",
    "gfu.vec.data = a.mat.Inverse(freedofs) * f.vec\n",
    "    \n",
    "uNum = gfu.components[0]\n",
    "l2error = sqrt(Integrate(InnerProduct(Pmat*(uNum-uSol), Pmat*(uNum-uSol))*ds, mesh))\n",
    "print (\"l2error : \", l2error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Visualization with VTK\n",
    "Create vtk-output file to be read by Paraview."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mesh.UnsetDeformation()\n",
    "vtk = VTKOutput(ma=mesh,\n",
    "                coefs=[lset_approx,deformation,uNum,uSol,rhs],\n",
    "                names=[\"P1-levelset\",\"deform\",\"u\",\"uSol\",\"rhs\"],\n",
    "                filename=\"vectorlaplacelagrange3d\",subdivision=2)\n",
    "vtk.Do()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<center>\n",
    "![title](graphics/tracefem_vector2.png)\n",
    "</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "ls -al vectorlaplacelagrange3d.vtu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "#%%bash\n",
    "#paraview vectorlaplacelagrange3d.vtk"
   ]
  },
  {
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}