NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(advanced:Hdiv)=\n",
    "# Mixed formulation for second order equations\n",
    "\n",
    "We have used $H(div)$ elements in the adaptivity tutorial. Here is used to introduce a mixed formulation for second order equations.\n",
    "\n",
    "Why mixed formulations?:\n",
    "* exact flux conservation\n",
    "* useful for a posteriori error estimates\n",
    "* model problem for $4^{th}$ order problems, Stokes, ...\n",
    "\n",
    "\n",
    "We consider the diffusion equation\n",
    "$$\n",
    "\n",
    "\\begin{array}{rcll}\n",
    "-\\text{div} \\lambda \\nabla u & = & f & \\text{ in } \\Omega \\\\\n",
    "u & = & u_D & \\text{ on } \\Gamma_D \\\\\n",
    "\\lambda \\frac{\\partial u}{\\partial n} & = & g & \\text{ on } \\Gamma_N\n",
    "\\end{array}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "***Primal variational formulation***\n",
    "\n",
    "Find $u \\in H^1,  u = u_D$ on $\\Gamma_D$ such that\n",
    "\n",
    "$$\n",
    "\\int_\\Omega \\lambda \\nabla u \\nabla v = \\int_\\Omega f v + \\int_{\\Gamma_N} g v\n",
    "\\quad \\forall v \\in H^1, v = 0 \\text{ on } \\Gamma_D\n",
    "$$\n",
    "\n",
    "We can recast the fisrt equation as a set of two equations introducing the flux $\\sigma = \\lambda \\nabla u$:\n",
    "\n",
    "Find scalar $u$ and the flux $\\sigma$ such that\n",
    "\n",
    "$$\n",
    "\\lambda^{-1} \\sigma = \\nabla u, \\quad \\text{div}  \\sigma = -f\n",
    "$$\n",
    "\n",
    "with boundary conditions\n",
    "\n",
    "$$\n",
    "\\sigma \\cdot n = g \\text{ on } \\Gamma_N, \\quad \\text{ and } \\quad\n",
    "u = u_D \\text{ on } \\Gamma_D \n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "***Mixed variational formulation***\n",
    "\n",
    "Find $\\sigma \\in H(\\text{div} )$ and $u \\in L_2$ such that $\\sigma_n = g$ on $\\Gamma_N$ and\n",
    "\n",
    "\\begin{array}{rlll}\n",
    "\\int_\\Omega \\lambda^{-1} \\sigma \\tau &+ \\int_\\Omega \\text{div}  \\tau u & = & 0 \\\\[0.4 cm]\n",
    "\\int_\\Omega \\text{div}  \\sigma v   & & = &  -\\int_\\Omega f v + \\int_{\\Gamma_D} u_D \\tau_n\n",
    "\\end{array}\n",
    "\n",
    "for all test-functions $\\tau \\in H(\\text{div} )$ and $v \\in L_2$ with $\\tau_n = 0$.\n",
    "\n",
    "Here $\\sigma_n$ is the normal trace operator $\\sigma \\cdot n |_{\\Gamma_N}$, which is meaningful in $H(\\text{div} )$.\n",
    "\n",
    "\n",
    "***The big-Bilinear Form***\n",
    "\n",
    "A Compact notation (and meaningful) notation in the context of FEM is the single-liner\n",
    "\n",
    "Find $(\\sigma, u) \\in H(\\text{div} ) \\times L_2$ such that $\\sigma_n = g$ on $\\Gamma_N$ and\n",
    "\n",
    "$$\n",
    "\\int_\\Omega \\lambda^{-1} \\sigma \\tau + \\text{div}  \\sigma v + \\text{div}  \\tau u = \n",
    "-\\int_\\Omega f v + \\int_{\\Gamma_D} u_D \\tau_n\n",
    "$$\n",
    "\n",
    "for all test-functions $(\\tau, v) \\in H(\\text{div} ) \\times L_2$ with $\\tau_n = 0$.\n",
    "\n",
    "the above can be seen as \n",
    "\n",
    "\\begin{array}{rlll}\n",
    "a((\\sigma, u), (\\tau, v)) & = & L((\\tau, v))\n",
    "\\end{array}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import netgen.occ\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MakeGeometryOCC():\n",
    "    # create a rectangle with boundary conditions\n",
    "    square = Circle((0, 0), 2).Face()\n",
    "    square.edges.name = \"outer\"\n",
    "    square.faces.name = \"air\"\n",
    "\n",
    "    # create a rectangle with boundary conditions\n",
    "    el1 = MoveTo(-0.4, 0.2).Rectangle(0.8, 0.1).Face()\n",
    "    el1.edges.name = \"el1\"\n",
    "    el1.vertices.name = \"el1\"\n",
    "\n",
    "    # create a rectangle with boundary conditions\n",
    "    el2 = MoveTo(-0.4, -0.2).Rectangle(0.8, 0.1).Face()\n",
    "    el2.edges.name = \"el2\"\n",
    "    el2.vertices.name = \"el2\"\n",
    "\n",
    "    dielec = MoveTo(-0.9, -0.1).Rectangle(1.8, 0.3).Face()\n",
    "    dielec.faces.name = \"dielec\"\n",
    "\n",
    "    air = square  # subtract the rectangles from the air rectangle\n",
    "    shape = Glue([air - dielec, dielec])\n",
    "    shape = shape - el1 - el2\n",
    "\n",
    "    # adding extra specifications on the shape\n",
    "    # predefined mesh size for the shape\n",
    "    Draw(shape)\n",
    "    return OCCGeometry(shape, dim=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "mesh = Mesh(MakeGeometryOCC().GenerateMesh(maxh=0.1))\n",
    "mesh.Curve(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lam = CoefficientFunction([1, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setup and solve primal problem:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving the mixed problem\n",
    "Define spaces for mixed problem. Note that essential boundary conditions  are set to the $H(\\text{div})$-component on the opposite boundary. Creating a space from a list of spaces generates a product space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order_flux=1\n",
    "V = HDiv(mesh, order=order_flux, dirichlet = \"outer\")\n",
    "Q = L2(mesh, order=order_flux-1)\n",
    "fesm = V*Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The space provides now a list of trial and test-functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma, u = fesm.TrialFunction()\n",
    "tau, v = fesm.TestFunction()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The normal vector is provided as a *special* coefficient function (which may be used only at the boundary). The orientation depends on the definition of the geometry. In 2D, it is the tangential vector rotated to the right, and is the outer vector in our case. Since every CoefficientFunction must know its vector-dimension, we have to provide the spatial dimension:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "normal = specialcf.normal(mesh.dim)\n",
    "print (normal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the forms on the product space. For the boundary term, we have to use the Trace operator, which provides the projection to normal direction. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the boundary conditions\n",
    "electrode = mesh.BoundaryCF({\"el1\": 1, \"el2\": -1}, default=0)\n",
    "ud = electrode\n",
    "\n",
    "am = BilinearForm((1/lam*sigma*tau + div(sigma)*v + div(tau)*u)*dx).Assemble()\n",
    "fm = LinearForm( ud*(tau.Trace()*normal)*ds(\"el.*\")).Assemble()\n",
    "\n",
    "gfm = GridFunction(fesm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The proxy-functions used for the forms know to which component of the product space they belong to. To access components of the solution, we have to unpack its components. They don't have own coefficient vectors, but refer to ranges of the big coefficient vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfsigma, gfu = gfm.components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just to try something:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now set the essential boundary conditions for the flux part:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = fm.vec.data - am.mat * gfm.vec\n",
    "gfm.vec.data += am.mat.Inverse(freedofs=fesm.FreeDofs(), inverse=\"pardiso\") * res\n",
    "# solvers.BVP(bf=am, lf=fm, gf=gfm)\n",
    "Draw (gfsigma, mesh, \"flux-mixed\")\n",
    "Draw (gfu, mesh, \"u-mixed\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the difference:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Post-processing for the scalar variable\n",
    "\n",
    "The scalar variable is approximated one order lower than the vector variable, what is its gradient. Knowing the gradient of a function more accurate, and knowing its mean value, one can recover the function itself. For this post-processing trick we refer to [Arnold+Brezzi 85]\n",
    "\n",
    "\n",
    "find $\\widehat u \\in P^{k+1, dc}$ and $\\widehat \\lambda \\in P^{0, dc}$ such that\n",
    "\n",
    "$$\n",
    "\\begin{array}{ccccl}\n",
    "\\int \\lambda \\nabla \\widehat u \\nabla \\widehat v & + & \\int \\widehat \\lambda \\widehat v & = & \\int \\sigma \\nabla \\widehat v & \\forall \\, \\widehat v\\\\\n",
    "\\int \\widehat u \\widehat \\mu & & & = & \\int u \\widehat \\mu & \\forall \\, \\widehat \\mu\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fespost_u = L2(mesh, order=order_flux+1)\n",
    "fespost_lam = L2(mesh, order=0)\n",
    "fes_post = fespost_u*fespost_lam\n",
    "\n",
    "(u,la),(v,mu) = fes_post.TnT()\n",
    "\n",
    "a = BilinearForm( (lam*grad(u)*grad(v)+la*v+mu*u)*dx).Assemble()\n",
    "f = LinearForm((gfsigma*grad(v)+gfu*mu)*dx).Assemble()\n",
    "\n",
    "gfpost = GridFunction(fes_post)\n",
    "gfpost.vec.data = a.mat.Inverse() * f.vec\n",
    "\n",
    "\n",
    "# if we zoom in the deformation we can observe theat the elements are not conforming\n",
    "Draw (gfpost.components[0], mesh, \"upost\", deformation=True)\n"
   ]
  },
  {
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}