{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "# 1.5 Spaces and forms on subdomains\n",
    "     \n",
    "\n",
    "In NGSolve, finite element spaces can be defined on subdomains. This is useful for multiphysics problems like fluid-structure interaction.\n",
    "\n",
    "In addition, bilinear or linear forms can be defined as integrals over regions. *Regions* are parts of the domain. They may be subdomains, or  parts of the domain boundary, or parts of subdomain interfaces.\n",
    "\n",
    "In this tutorial, you will learn about \n",
    "\n",
    "- defining finite element spaces on `Region`s,\n",
    "- `Compress`ing such spaces,\n",
    "- integrating on `Regions`s, and \n",
    "- set operations on `Region` objects via `Mask`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *   # Opencascade for geometry modeling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "### Naming subdomains and boundaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "We define a geometry with multiple regions and assign names to these regions below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outer = Rectangle(2, 2).Face()\n",
    "outer.edges.name=\"outer\"\n",
    "outer.edges.Max(X).name = \"r\"\n",
    "outer.edges.Min(X).name = \"l\"\n",
    "outer.edges.Min(Y).name = \"b\"\n",
    "outer.edges.Max(Y).name = \"t\"\n",
    "\n",
    "inner = MoveTo(1, 0.9).Rectangle(0.3, 0.5).Face()\n",
    "inner.edges.name=\"interface\"\n",
    "outer = outer - inner\n",
    "\n",
    "inner.faces.name=\"inner\"\n",
    "inner.faces.col = (1, 0, 0)\n",
    "outer.faces.name=\"outer\"\n",
    "\n",
    "geo = Glue([inner, outer])\n",
    "Draw(geo);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(OCCGeometry(geo, dim=2).GenerateMesh(maxh=0.2))\n",
    "Draw(mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "This mesh has two subdomains that we named \"inner\" and \"outer\" during the geometry construction. They can be identified as `Region` objects.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.Materials(\"inner\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The bottom, right, top and left parts of the outer rectangle's boundaries define boundary segments, which we respectively labeled  \"b\", \"r\", \"t\", \"l\". They are also instances of `Region` objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.Boundaries(\"b\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "### A finite element space on a subdomain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fes1 = H1(mesh, definedon=\"inner\")\n",
    "u1 = GridFunction(fes1, \"u1\")\n",
    "u1.Set(x*y)\n",
    "Draw(u1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Note how $u_1$ is displayed only in the `inner` region."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Do not be confused with `ndof` for spaces on subdomains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fes = H1(mesh)\n",
    "fes1.ndof, fes.ndof"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Although `ndof`s are the same for spaces on one  subdomain and the whole domain, `FreeDofs` show that the real number of degrees of freedom for the subdomain space is much smaller:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(fes1.FreeDofs())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "One can also glean this information by examining `CouplingType` of each degrees of freedom of the space, which reveals that there are many unknowns of type `COUPLING_TYPE.UNUSED_DOF`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "for i in range(fes1.ndof):\n",
    "    print(fes1.CouplingType(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "It is possible to remove these dofs (thus making `GridFunction` vectors smaller) as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fescomp = Compress(fes1)\n",
    "print(fescomp.ndof)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "### Integrating on regions "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "You have already seen how boundary regions are used in setting Dirichlet boundary conditions.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=3, dirichlet=\"b|l|r\")\n",
    "u, v = fes.TnT()\n",
    "gfu = GridFunction(fes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Such boundary regions or subdomains can also serve as domains of integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "term1 = u1 * v * dx(definedon=mesh.Materials(\"inner\")) \n",
    "term2 = 0.1 * v * ds(definedon=mesh.Boundaries(\"t\"))\n",
    "f = LinearForm(term1 + term2).Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Here the functional $f$ is defined as a sum of two integrals, one over the inner subdomain, and another over the top boundary:\n",
    "$$\n",
    "f(v) = \\int_{\\Omega_{inner}} u_1 v \\; dx + \\int_{\\Gamma_{top}} \\frac{v}{10}\\; ds\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The python objects `dx` and `ds` represent volume and boundary integration, respectively. They admit a keyword argument `definedon` that is **either** a `Region` **or** a `str`, as can be seen from the documentation: look for \n",
    "`definedon: Optional[Union[ngsolve.comp.Region, str]]` in the help output below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "help(dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Thus \n",
    "- `dx(definedon=mesh.Materials(\"inner\"))` may be replaced by \n",
    "`dx('inner')`,\n",
    "\n",
    "and similarly, \n",
    "\n",
    "- `ds(definedon=mesh.Boundaries(\"t\")` may be replaced by  `ds('t')`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm(u1*v*dx('inner') + 0.1*v*ds('t')).Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a = BilinearForm(fes)\n",
    "a += grad(u)*grad(v)*dx\n",
    "a.Assemble()\n",
    "gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec\n",
    "Draw(gfu);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "### Operations with `Region`s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Query mesh for all regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mesh.GetBoundaries()   # list boundary/interface regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mesh.GetMaterials()    # list all subdomains"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Print region information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(mesh.Materials(\"inner\").Mask())\n",
    "print(mesh.Materials(\"[a-z]*\").Mask())  # can use regexp\n",
    "print(mesh.Boundaries('t|l').Mask())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Add regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "io = mesh.Materials(\"inner\") + mesh.Materials(\"outer\")\n",
    "print(io.Mask())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Take complement of a region:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "c = ~mesh.Materials(\"inner\")\n",
    "print(c.Mask())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Subtract regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "diff = mesh.Materials(\"inner|outer\") - mesh.Materials(\"outer\")\n",
    "print(diff.Mask())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Set a piecewise constant `CoefficientFunction` using the subdomains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "domain_values = {'inner': 3.7,  'outer': 1}\n",
    "cf = mesh.MaterialCF(domain_values)\n",
    "Draw(cf, mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Use boundary regions to define a coefficient function: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Make a linear function equaling 2 at the bottom right vertex\n",
    "# of (0,2) x (0,2) and equals zero at the remaining vertices:\n",
    "bdry_values = {'b': x, 'r': 2-y}\n",
    "cf = mesh.BoundaryCF(bdry_values, default=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Draw`ing such a `BoundaryCF` does not show anything useful. However, we can use a `GridFunction` and `Set`, which then lets us view an extension of these boundary values into the domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "g = GridFunction(H1(mesh), name='bdry')\n",
    "g.Set(cf, definedon=mesh.Boundaries('b|r'))\n",
    "Draw(g);"
   ]
  }
 ],
 "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.11.3"
  },
  "name": "subdomains.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}