{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adaptivity\n",
    "==="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from netgen.occ import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the geometry by 2D Netgen-OpenCascade modeling: Define rectangles and polygones, and glue them together to one shape object. OCC maintains the full geometry topology of vertices, edges, faces and solids."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MakeGeometry():\n",
    "    base = Rectangle(1, 0.6).Face()\n",
    "    chip = MoveTo(0.5,0.15).Line(0.15,0.15).Line(-0.15,0.15).Line(-0.15,-0.15).Close().Face()\n",
    "    top = MoveTo(0.2,0.6).Rectangle(0.6,0.2).Face()\n",
    "    base -= chip\n",
    "\n",
    "    base.faces.name=\"base\"\n",
    "    chip.faces.name=\"chip\"\n",
    "    chip.faces.col=(1,0,0)\n",
    "    top.faces.name=\"top\"\n",
    "    geo = Glue([base,chip,top])\n",
    "    geo.edges.name=\"default\"\n",
    "    geo.edges.Min(Y).name=\"bot\"\n",
    "    Draw(geo)\n",
    "    return geo\n",
    "\n",
    "geo = MakeGeometry()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Piece-wise constant coefficients in sub-domains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = geo.GenerateMesh(maxh=0.2, dim=2)\n",
    "print (mesh.GetMaterials())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=3, dirichlet=\"bot\", autoupdate=True)\n",
    "u, v = fes.TnT()\n",
    "\n",
    "lam = mesh.MaterialCF( { \"base\" : 1, \"chip\" : 1000, \"top\" : 20 } )\n",
    "a = BilinearForm(lam*grad(u)*grad(v)*dx)\n",
    "\n",
    "# heat-source in inner subdomain\n",
    "f = LinearForm(1*v*dx(definedon=\"chip\"))\n",
    "\n",
    "c = preconditioners.MultiGrid(a, inverse=\"sparsecholesky\")\n",
    "\n",
    "gfu = GridFunction(fes) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assemble and solve problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveBVP():\n",
    "    a.Assemble()\n",
    "    f.Assemble()\n",
    "    inv = CGSolver(a.mat, c.mat)\n",
    "    gfu.vec.data = inv * f.vec\n",
    "    \n",
    "SolveBVP()\n",
    "Draw (gfu, mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gradient recovery error estimator: Interpolate finite element flux \n",
    "\n",
    "$$\n",
    "q_h := I_h (\\lambda \\nabla u_h)\n",
    "$$\n",
    "\n",
    "and take difference as element error indicator:\n",
    "\n",
    "$$\n",
    "\\eta_T := \\tfrac{1}{\\lambda} \\| q_h - \\lambda \\nabla u_h \\|_{L_2(T)}^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = []    # l = list of estimated total error\n",
    "space_flux = HDiv(mesh, order=2, autoupdate=True)\n",
    "gf_flux = GridFunction(space_flux, \"flux\", autoupdate=True)\n",
    "\n",
    "def CalcError():\n",
    "    \n",
    "    flux = lam * grad(gfu)   # the FEM-flux      \n",
    "    gf_flux.Set(flux)        # interpolate into H(div)\n",
    "    \n",
    "    # compute estimator:\n",
    "    err = 1/lam*(flux-gf_flux)*(flux-gf_flux)\n",
    "    eta2 = Integrate(err, mesh, VOL, element_wise=True)\n",
    "    l.append ((fes.ndof, sqrt(sum(eta2))))\n",
    "    print(\"ndof =\", fes.ndof, \" toterr =\", sqrt(sum(eta2)))\n",
    "    \n",
    "    # mark for refinement:\n",
    "    maxerr = max(eta2) \n",
    "    # marking with Python loop:\n",
    "    # for el in mesh.Elements():\n",
    "    #    mesh.SetRefinementFlag(el, eta2[el.nr] > 0.25*maxerr)\n",
    "    \n",
    "    # marking using numpy vectorization:\n",
    "    mesh.ngmesh.Elements2D().NumPy()[\"refine\"] = eta2.NumPy() > 0.25*maxerr\n",
    "  \n",
    "CalcError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adaptive loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "level = 0\n",
    "while fes.ndof < 50000:  \n",
    "    mesh.Refine()\n",
    "    SolveBVP()\n",
    "    CalcError()\n",
    "    level = level+1\n",
    "    if level%5 == 0:\n",
    "        Draw (gfu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfu);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.xlabel(\"ndof\")\n",
    "plt.ylabel(\"H1 error-estimate\")\n",
    "ndof,err = zip(*l)\n",
    "plt.plot(ndof,err, \"-*\")\n",
    "\n",
    "plt.ion()\n",
    "plt.show();"
   ]
  },
  {
   "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.13.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}