{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.7 Facet spaces and hybrid methods\n",
    "\n",
    "Mixed methods for second order problems lead to saddle point problems, and indefinite matrices. By hybridization one obtains a positive definite system again. It's structure is similar to the non-conforming $P^1$ method, but hybridization works for any order. See text-book by Brezzi and Fortin.\n",
    "\n",
    "One skips the normal-continuity of the $H(div)$ variable, and reinforces it by a Lagrange parameter. This leads to the following discrete system:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find $\\sigma, u, \\widehat u \\in \\Sigma_h \\times V_h \\times F_h$:\n",
    "\n",
    "$$\n",
    "\\DeclareMathOperator{\\Div}{div}\n",
    "\\begin{array}{ccccccll}\n",
    "\\int \\lambda^{-1} \\sigma \\tau & + & \\sum_T \\int_T \\Div \\tau \\, u & + & \\sum_F \\int_F [\\tau_n] \\widehat u & = & 0 & \\forall \\, \\tau \\in \\Sigma \\\\\n",
    "\\int \\Div \\sigma \\, v &&&&& = & \\int f v & \\forall \\, v \\in V_h \\\\\n",
    "\\int [ \\sigma_n ] \\, \\widehat v &&&&& = & \\int_{\\Gamma_n} g \\widehat v & \\forall \\, \\widehat v \\in F_h\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "where $\\Sigma_h$ is an discontinuous $H(div)$ finite element space, $V_h$ a sub-space of $L_2$, and $F_h$ consists of polynomials on every edge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "same example as in 'mixed':"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source = sin(pi*x)\n",
    "ud = CF(5)\n",
    "g = mesh.BoundaryCF( {\"left\" : y*(1-y)}, default=0)\n",
    "lam = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "define spaces: \n",
    "\n",
    "* The *Discontinuous* FESpace-wrapper generates an element-wise discontinuous space\n",
    "* FacetFESpace lives only on facets (i.e. faces in 3D, edges in 2D, points in 1D)\n",
    "\n",
    "Boundary conditions are now posed for the facet-space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order = 3\n",
    "V = Discontinuous (HDiv(mesh, order=order))\n",
    "Q = L2(mesh, order=order-1)\n",
    "F = FacetFESpace(mesh, order=order, dirichlet=\"bottom\")\n",
    "X = V*Q*F\n",
    "print (\"sigmadofs:\", X.Range(0))\n",
    "print (\"udofs:    \", X.Range(1))\n",
    "print (\"uhatdofs: \", X.Range(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assemble forms. The jump-term is rewritten as\n",
    "$$\n",
    "\\sum_F \\int_F [\\sigma_n] v = \\sum_T \\int_{\\partial T} \\sigma_n v\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma,u,uhat = X.TrialFunction()\n",
    "tau,v,vhat = X.TestFunction()\n",
    "\n",
    "a = BilinearForm(X, condense=True)\n",
    "a += (1/lam * sigma*tau + div(sigma)*v + div(tau)*u) * dx\n",
    "n = specialcf.normal(mesh.dim)\n",
    "a += (-sigma*n*vhat-tau*n*uhat) * dx(element_boundary=True)\n",
    "\n",
    "c = Preconditioner(a, \"bddc\")\n",
    "\n",
    "f = LinearForm(X)\n",
    "f += -source*v * dx - g*vhat.Trace() * ds\n",
    "\n",
    "a.Assemble()\n",
    "print (\"A non-zero\", a.mat.nze)\n",
    "\n",
    "gfu = GridFunction(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve system. Either we leave everything to the sparse direct solver, or use CG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f.Assemble()\n",
    "gfu.components[2].Set(ud, BND)\n",
    "\n",
    "if a.condense:\n",
    "\n",
    "    fmod = (f.vec + a.harmonic_extension_trans * f.vec).Evaluate()\n",
    "    solvers.CG(mat=a.mat, pre=c.mat, rhs=fmod, sol=gfu.vec, initialize=False)\n",
    "    \n",
    "    gfu.vec.data += a.harmonic_extension * gfu.vec\n",
    "    gfu.vec.data += a.inner_solve * f.vec\n",
    "\n",
    "else:\n",
    "    \n",
    "    r = f.vec - a.mat * gfu.vec\n",
    "    inv = a.mat.Inverse(freedofs=X.FreeDofs())\n",
    "    gfu.vec.data += inv * r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfu.components[0], mesh, \"sigma\")\n",
    "Draw (gfu.components[1], mesh, \"u\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.10.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}