{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "11.2.1 Dirichlet Laplace Indirect Method\n",
    "=============================\n",
    "**keys**: homogeneous Dirichlet bvp, single layer potential, ACA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from ngsolve.bem import *\n",
    "from ngsolve import Projector, Preconditioner\n",
    "from ngsolve.krylovspace import CG"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "We consider the Dirichlet boundary value problem \n",
    "\n",
    "$$ \\left\\{ \\begin{array}{rcl l} \\Delta u &=& 0, \\quad &\\Omega \\subset \\mathbb R^3\\,,\\\\ \\gamma_0 u&=& u_0, \\quad &\\Gamma = \\partial \\Omega\\,.\\end{array} \\right. $$ \n",
    "\n",
    "Let us choose the following ansatz for the solution $u\\in H^1(\\Omega)$ (indirect ansatz) \n",
    "\n",
    "$$ u(x) = \\underbrace{ \\int\\limits_\\Gamma \\displaystyle{\\frac{1}{4\\,\\pi}\\, \\frac{1}{\\| x-y\\|} } \\, j(y)\\, \\mathrm{d}\\sigma_y }_{\\displaystyle{ \\mathrm{SL}(j) } }$$ \n",
    "\n",
    "and solve for the density $j\\in H^{-\\frac12}(\\Gamma)$ by the boundary element method, i.e. the numerical solution of the variational formulation \n",
    "\n",
    "$$ \\forall \\, v\\in H^{-\\frac12}(\\Gamma): \\quad \\left\\langle \\gamma_0 \\left(\\mathrm{SL}(j)\\right), v \\right\\rangle_{-\\frac12} = \\left\\langle u_0, v\\right\\rangle_{-\\frac12} \\,. $$\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "Define the geometry $\\Omega \\subset \\mathbb R^3$ and create a mesh:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "sp = Sphere( (0,0,0), 1)\n",
    "mesh = Mesh( OCCGeometry(sp).GenerateMesh(maxh=0.3)).Curve(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "Create test and trial function finite element spaces for $H^{-\\frac12}(\\Gamma)$ according to the given mesh:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fesL2 = SurfaceL2(mesh, order=4, dual_mapping=True)\n",
    "u,v = fesL2.TnT()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "Define Dirichlet data $u_0$ and compute the right hand side vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "u0 = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )\n",
    "rhs = LinearForm (u0*v.Trace()*ds(bonus_intorder=3)).Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "The discretisation of the above variational formulation leads to a system of linear equations, ie \n",
    "\n",
    "$$ \\mathrm{V} \\, \\mathrm{j} =  \\mathrm{rhs} \\,, $$ \n",
    "where $\\mathrm{V}$ is the single layer potential operator. $\\mathrm V$ is regular and symmetric."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "**Demo 1**: Assemble the single layer operator $V$ as dense matrix and solve for unknwon density $j$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "j = GridFunction(fesL2)\n",
    "pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()\n",
    "with TaskManager():\n",
    "    # V = SingleLayerPotentialOperator(fesL2, intorder=10)\n",
    "    V = LaplaceSL(u*ds)*v*ds\n",
    "    CG(mat = V.mat, pre=pre, rhs = rhs.vec, sol=j.vec, tol=1e-8, maxsteps=200, initialize=False, printrates=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (j, order=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "**Notes:**\n",
    "- For details on the analysis of boundary integral equations derived from elliptic partial differential equations, see for instance [Strongly Elliptic Systems and Boundary Integral Equations](https://www.cambridge.org/de/universitypress/subjects/mathematics/differential-and-integral-equations-dynamical-systems-and-co/strongly-elliptic-systems-and-boundary-integral-equations?format=HB&isbn=9780521663328).\n",
    "- The integration of singular pairings is done as proposed in [Randelementmethoden](https://link.springer.com/book/9783519003687)\n",
    "- The adaptive cross approximation is done as proposed in [Hierarchical Matrices](https://link.springer.com/book/10.1007/978-3-540-77147-0).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}