{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "11.2.2. Dirichlet Laplace Direct Method \n",
    "==========================\n",
    "**keys**: homogeneous Dirichlet bvp, double and single layer potential, unknown Neumann data"
   ]
  },
  {
   "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.solvers import CG"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "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. $$ "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "Let us choose the following ansatz for the solution $u\\in H^1(\\Omega)$ (direct ansatz) \n",
    "\n",
    "$$ u(x) = \\underbrace{ \\int\\limits_\\Gamma \\displaystyle{\\frac{1}{4\\,\\pi}\\, \\frac{1}{\\| x-y\\|} } \\, u_1(y)\\, \\mathrm{d}\\sigma_y}_{\\displaystyle{ \\mathrm{ SL}(u_1) }} - \\underbrace{ \\int\\limits_\\Gamma \\displaystyle{\\frac{1}{4\\,\\pi}\\, \\frac{\\langle n(y) , x-y\\rangle }{\\| x-y\\|^3} } \\, u_0(y)\\, \\mathrm{d}\\sigma_y}_{\\displaystyle{ \\mathrm{DL}(u_0) }}$$ \n",
    "\n",
    "and solve for the Neumann data $u_1 \\in H^{-\\frac12}(\\Gamma)$ by the boundary element method, i.e., \n",
    "\n",
    "$$ \\forall \\, v\\in H^{-\\frac12}(\\Gamma): \\quad \\left\\langle \\gamma_0 \\left(\\mathrm{SL}(u_1)\\right), v \\right\\rangle_{-\\frac12}= \\left\\langle u_0, v\\right\\rangle_{-\\frac12} + \\left\\langle \\gamma_0 \\left(\\mathrm{DL}(u_0)\\right), v\\right\\rangle_{-\\frac12}\\,. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "Define the geometry $\\Omega \\subset \\mathbb R^3$ and create a mesh:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "sp = Sphere( (0,0,0), 1)\n",
    "mesh = Mesh( OCCGeometry(sp).GenerateMesh(maxh=0.2)).Curve(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "Create the finite element spaces for $H^{-\\frac12}(\\Gamma)$ and $H^{\\frac12}(\\Gamma)$ according to the given mesh:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "fesL2 = SurfaceL2(mesh, order=3, dual_mapping=True)\n",
    "u,v = fesL2.TnT()\n",
    "fesH1 = H1(mesh, order=4)\n",
    "uH1,vH1 = fesH1.TnT()\n",
    "print (\"ndofL2 = \", fesL2.ndof, \"ndof H1 = \", fesH1.ndof)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "Compute the interpolation of exact Dirichlet data $u_0$ in finite element space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "uexa = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )\n",
    "u0 = GridFunction(fesH1)\n",
    "u0.Interpolate (uexa)\n",
    "Draw (u0, mesh, draw_vol=False, order=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "The discretisation of the above variational formulation leads to a system of linear equations, ie \n",
    "\n",
    "$$ \\mathrm{V} \\, \\mathrm{u}_1 = \\left( \\frac12 \\,\\mathrm{M} + \\mathrm{K} \\right) \\, \\mathrm{u}_0\\,, $$ \n",
    "where the linear operators are as follows\n",
    "- $\\mathrm{V}$ is the single layer operator. $\\mathrm V$ is regular and symmetric.\n",
    "- $\\mathrm{M}$ is a mass matrix.\n",
    "- $\\mathrm{K}$ is the double layer operator. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "We approximate the linear operators as hmatrices and solve for the Neumann data $u_1$ with an iterative solver:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "u1 = GridFunction(fesL2)\n",
    "pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()\n",
    "\n",
    "with TaskManager():\n",
    "    # V = SingleLayerPotentialOperator(fesL2, intorder=12)\n",
    "    V = LaplaceSL(u*ds)*v*ds\n",
    "    \n",
    "    M = BilinearForm( uH1 * v.Trace() * ds(bonus_intorder=3)).Assemble()\n",
    "    # K = DoubleLayerPotentialOperator(fesH1, fesL2, intorder=12)\n",
    "    K = LaplaceDL(uH1*ds)*v*ds\n",
    "    rhs = ( (0.5 * M.mat + K.mat)*u0.vec).Evaluate()\n",
    "    CG(mat = V.mat, pre=pre, rhs = rhs, sol=u1.vec, tol=1e-8, maxsteps=50, initialize=False, printrates=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (u1, mesh, draw_vol=False, order=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "Let's have a look at the exact Neumann data and compute the error of the numerical solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "graduexa = CF( (uexa.Diff(x), uexa.Diff(y), uexa.Diff(z)) )\n",
    "n = specialcf.normal(3)\n",
    "u1exa = graduexa*n\n",
    "Draw (u1exa, mesh, draw_vol=False, order=3);\n",
    "print (\"L2-error =\", sqrt (Integrate ( (u1exa-u1)**2, mesh, BND)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "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
}