{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# 2.13 Interface resisitivity\n",
    "\n",
    "We model a problem with resistivity interface conditions\n",
    "\\begin{eqnarray*} \n",
    "\\lambda_l \\frac{\\partial u_l}{\\partial n_l} + \n",
    "\\lambda_r \\frac{\\partial u_r}{\\partial n_r} & = & 0 \\\\\n",
    "\\alpha (u_l - u_r) & = & \\lambda_l \\frac{\\partial u_l}{\\partial n_l}\n",
    "\\end{eqnarray*}\n",
    "\n",
    "This is easily modeled by an additional Robin-like term on the interface:\n",
    "\n",
    "$$\n",
    "\\int_\\Omega \\lambda \\nabla u \\nabla v \\, dx + \\int_\\gamma \\alpha [u][v] \\, ds,\n",
    "$$\n",
    "where $[u] = u_l - u_r$ denotes the jump across the interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from netgen.occ import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "inner = Rectangle(1,1).Face()\n",
    "inner.edges.name=\"interface\"\n",
    "inner.edges.Max(X).name=\"interface_right\"\n",
    "\n",
    "outer = MoveTo(-1,-1).Rectangle(3,3).Face()\n",
    "outer.edges.name=\"dir\"\n",
    "outer = outer-inner\n",
    "\n",
    "inner.faces.name=\"inner\"\n",
    "outer.faces.name=\"outer\"\n",
    "geo = Glue ([inner,outer])\n",
    "mesh = Mesh(OCCGeometry(geo, dim=2).GenerateMesh(maxh=0.2))\n",
    "\n",
    "Draw (mesh);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print (mesh.GetMaterials(), mesh.GetBoundaries())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "### Interfaces without boundaries\n",
    "\n",
    "In the first case, we assume the interface is the boundary of the inner domain. For this case, we define spaces on both sides separately, and link them via the Robin - term.\n",
    "\n",
    "In our problem we choose zero conductivity across the right edge of the interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fesi = H1(mesh, order=3, definedon=mesh.Materials(\"inner\"))\n",
    "feso = H1(mesh, order=3, definedon=mesh.Materials(\"outer\"), dirichlet=\"dir\")\n",
    "fes = fesi*feso\n",
    "\n",
    "(ui,uo), (vi,vo) = fes.TnT()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 10\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += grad(ui)*grad(vi)*dx(\"inner\")\n",
    "a += grad(uo)*grad(vo)*dx(\"outer\")\n",
    "a += alpha * (ui-uo) * (vi-vo) * ds(\"interface\")\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(1*vi*dx(\"inner\")).Assemble()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "gfu.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs()) * f.vec\n",
    "\n",
    "sol = mesh.MaterialCF({\"inner\":gfu.components[0], \"outer\":gfu.components[1]})\n",
    "Draw (sol, mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "## Interfaces with boundaries\n",
    "\n",
    "Now we model an interface resistivity only at the right edge. For this we need a space discontinuous only at that edge. We obtain that by enriching an H1-space with functions defined by values on the interface, and supported only on one side. The jump is exactly the enrichment function.\n",
    "\n",
    "We define the enrichment space first on one domain, and then restrict degrees of freedom to the interior of the right interface edge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes1 = H1(mesh, order=3, dirichlet=\"dir\")\n",
    "\n",
    "fesext = H1(mesh, order=3, definedon=\"outer\")\n",
    "actdofs = fesext.GetDofs(mesh.Boundaries(\"interface_right\")) & \\\n",
    "    ~fesext.GetDofs(mesh.Boundaries(\"interface\"))\n",
    "fesext = Compress (fesext, active_dofs=actdofs)\n",
    "print (actdofs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "A typical enrichment function looks like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "gfext = GridFunction(fesext)\n",
    "gfext.Set(sin(10*y))\n",
    "Draw (gfext);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = fes1 * fesext\n",
    "(u,uext), (v,vext) = fes.TnT()\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += grad(u)*grad(v)*dx(\"inner\")\n",
    "a += (grad(u)+grad(uext))* (grad(v)+grad(vext))*dx(\"outer\")\n",
    "a += 3 * uext*vext * ds(\"interface_right\")\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(1*v*dx(\"inner\")).Assemble()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "gfu.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs()) * f.vec\n",
    "\n",
    "sol = mesh.MaterialCF({\"inner\":gfu.components[0], \"outer\":gfu.components[0]+gfu.components[1]})\n",
    "Draw (sol, mesh);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "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": 5
}