NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A Slighly More Complex Problem\n",
    "\n",
    "Let's create now a capacitor problem in which we prescribe the potential at the electrodes. The equation that we need to solve is the following:\n",
    "\n",
    "Find $u$ such that\n",
    "\\begin{align*} -\\varepsilon \\Delta u &= 0 \\quad \\text{in } \\Omega, \\\\ u &= 1 \\quad \\text{on } \\text{Electrode 1},\\\\ u &= -1 \\quad \\text{on } \\text{Electrode 2},\\\\ n\\cdot \\nabla u &= 0 \\quad \\text{on } \\text{outer}, \\end{align*}\n",
    "\n",
    "With permettivity given by:\n",
    "\\begin{align*} \\varepsilon = \\begin{cases} 1 & \\text{in air},\\\\ 2 & \\text{in dielectric}. \\end{cases} \\end{align*}\n",
    "And domain as in figure below:\n",
    "<div style=\"text-align:center;\">\n",
    "<svg width=\"400\" height=\"400\" viewBox=\"-2.1 -2.1 4.2 4.2\">\n",
    "    <!-- Circle from (0, 0) with radious 2 and call it \"air\" with black border -->\n",
    "    <circle cx=\"0\" cy=\"0\" r=\"2\" fill=\"white\" stroke=\"black\" stroke-width=\"0.01\"/>\n",
    "    <!-- Electrodes -->\n",
    "    <rect x=\"-0.4\" y=\"0.2\" width=\"0.8\" height=\"0.1\" fill=\"black\"/>\n",
    "    <rect x=\"-0.4\" y=\"-0.2\" width=\"0.8\" height=\"0.1\" fill=\"black\"/>\n",
    "    <!-- Dielectric -->\n",
    "    <rect x=\"-0.9\" y=\"-0.1\" width=\"1.8\" height=\"0.3\" fill=\"gray\"/>\n",
    "    <!--Add the names air Electrode 1 Electrode 2 Dielectric -->\n",
    "    <text x=\"-0.5\" y=\"0.7\" font-size=\"0.09\" fill=\"black\">Air</text>\n",
    "    <text x=\"-0.4\" y=\"0.3\" font-size=\"0.09\" fill=\"white\">Electrode</text>\n",
    "    <text x=\"-0.9\" y=\"0.1\" font-size=\"0.09\" fill=\"black\">Dielectric</text>\n",
    "</svg>\n",
    "</div>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "square = Circle( (0,0), 2).Face() # create a rectangle with boundary conditions\n",
    "square.edges.name = \"outer\"\n",
    "square.faces.name = \"air\"\n",
    "\n",
    "el1 = MoveTo(-0.4, 0.2).Rectangle( 0.8, 0.1).Face() # create a rectangle with boundary conditions\n",
    "el1.edges.name = \"el1\"\n",
    "el1.vertices.name = \"el1\"\n",
    "\n",
    "el2 = MoveTo(-0.4, -0.2).Rectangle( 0.8, 0.1).Face() # create a rectangle with boundary conditions\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",
    "\n",
    "shape.edges[\"el.*\"].maxh = 0.1 # set the mesh size of the edges with the name \"el.*\" to 0.05\n",
    "shape.vertices[\"el.*\"].maxh = 0.05 # set the mesh size of the vertices with the name \"el.*\" to 0.05\n",
    "\n",
    "\n",
    "Draw(shape); # draw the shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "mesh = Mesh(OCCGeometry(shape, dim = 2).GenerateMesh(maxh=0.5)).Curve(4) # create a mesh from the shape\n",
    "\n",
    "eps = CoefficientFunction( [1, 2] ) # define the coefficient function # try sin(5*x) instead of 2 \n",
    "\n",
    "Draw(eps, mesh, \"eps\"); # draw the coefficient function\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bilinear and liner forms now are :\n",
    "\\begin{align*} a(u,v) &= \\int_{\\Omega} \\varepsilon \\nabla u \\cdot \\nabla v \\, dx,\\\\ f(v) &= 0. \\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order = 3\n",
    "fes = H1(mesh, order=order, dirichlet=\"el.*\") # define the finite element space\n",
    "\n",
    "u ,v = fes.TnT() # define the trial and test functions\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += eps * grad(u) * grad(v) * dx\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f.Assemble() # not needed it is zero"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We dont have homogeneous Dirichlet boundary conditions, so we need a homogenization technique.\n",
    "\n",
    "Homogenization technique in a nutshell:\n",
    "Suppose we the solution $u_{\\text{sol}} = u_0 + u_{D}$, where $u_0$ is the solution of the problem with homogeneous boundary conditions and $ u_{D}$ is the solution \n",
    "\n",
    "\\begin{align*}\n",
    "A(u_{\\text{sol}}, v ) = f(v)\\\\\n",
    "A(u_0 + u_D, v ) = f(v)\\\\\n",
    "A(u_0, v )  = f(v)-A(u_D, v )\\\\\n",
    "\\end{align*}\n",
    "\n",
    "After the discretization is equivalent to solve the following linear system:\n",
    "\\begin{align*}\n",
    "u_0 = A^{-1}(f-Au_D)\n",
    "\\end{align*}\n",
    "The solution can be found as\n",
    "\\begin{align*}\n",
    "u_\\text{sol} = A^{-1}(f-Au_{D}) + u_D\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "gfu = GridFunction(fes) # define the grid function\n",
    "electrode = mesh.BoundaryCF({\"el1\":1, \"el2\":-1}, default = 0 ) # define the boundary conditions\n",
    "gfu.Set(electrode, definedon=mesh.Boundaries(\"el.*\")) # set the boundary conditions\n",
    "\n",
    "Draw(gfu, deformation =True); # draw the grid function\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "gfu.vec.data += a.mat.Inverse(freedofs = fes.FreeDofs()) * (f.vec - a.mat * gfu.vec) # solve the system\n",
    "Draw(gfu , deformation =True, scale = 1); # draw the solution\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Electic field\n",
    "Draw(-eps*grad(gfu) , mesh, \"E\", min = 0, max = 4); # draw the solution\n",
    "\n",
    "# if you zoom in you will see the difference in the two materials."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some visualization tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve.webgui import Draw, FieldLines, AddFieldLines\n",
    "fes_flux = HDiv(mesh, order=order-1) # Same as \n",
    "#fes_flux = H1(mesh, order=order-1, dim=mesh.dim) \n",
    "# fes_flux = HDiv(mesh, order=order-1)\n",
    "\n",
    "\n",
    "flux = GridFunction(fes_flux, name=\"flux\")\n",
    "flux.Set( eps*grad(gfu) )\n",
    "\n",
    "N = 10\n",
    "#fac = 0 if mesh.dim == 2 else 1\n",
    "#p = [(-1+4*i/N,-2+4*j/N,fac * 2*k/N) for i in range(1,2*N) for j in range(1,2*N) for k in range(1,N)]\n",
    "p = [(i / N, -1 + (2*j)/ N, 0) for i in range(N+1) for j in range(N+1) ] \n",
    "\n",
    "fieldlines = flux._BuildFieldLines(mesh, p, num_fieldlines=N**3//5, randomized=False, length=0.3)\n",
    "\n",
    "clipping = { \"clipping\" : { \"y\":0, \"z\":-1} }\n",
    "\n",
    "\n",
    "Draw(-eps*grad(gfu), mesh,  \"X\", draw_vol=True, draw_surf=True, objects=[fieldlines], \\\n",
    "     autoscale=True, min = 0, max = 10, settings={\"Objects\": {\"Surface\": False}}, **clipping);\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ngspy",
   "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
}