{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "11.2.4. Neumann Laplace Direct Method \n",
    "======================================================\n",
    "**keys**: homogeneous Neumann bvp, hypersingular operator, unknown Dirichlet 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 Neumann boundary value problem \n",
    "\n",
    "$$ \\left\\{ \\begin{array}{rlc l} \\Delta u &=& 0, \\quad &\\Omega \\subset \\mathbb R^3\\,,\\\\ \\gamma_1 u&=& u_1, \\quad &\\Gamma = \\partial \\Omega\\,.\\end{array} \\right. $$ \n",
    "\n",
    "The unique solution $u\\in H^1(\\Omega)$ can be written in the following form (representation forumula) \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_1( y)\\, \\mathrm{d}\\sigma_y}_{\\displaystyle{ \\mathrm{DL}(u_1) }}\\,.$$ \n",
    "\n",
    "Let's carefully apply the Neumann trace $\\gamma_1$ to the representation formula and solve the resulting boundary integral equation for the Dirichlet trace $u_0$ of $u$ by discretisation of the following variational formulation: \n",
    "\n",
    "$$ \\forall \\, v\\in H^{\\frac12}(\\Gamma):  \\left\\langle v, \\gamma_1 \\left(\\mathrm{DL}(u_0)\\right) \\right\\rangle_{-\\frac12}  = \\left\\langle u_1, v\\right\\rangle_{-\\frac12} - \\left\\langle v, \\gamma_1 \\left(\\mathrm{SL}(u_1)\\right) \\right\\rangle_{-\\frac12}\\,. $$"
   ]
  },
  {
   "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 = Glue(Sphere( (0,0,0), 1).faces)\n",
    "mesh = Mesh( OCCGeometry(sp).GenerateMesh(maxh=0.2)).Curve(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "Define conforming finite element spaces for $H^{\\frac12}(\\Gamma)$ and $H^{-\\frac12}(\\Gamma)$ respectively:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fesL2 = SurfaceL2(mesh, order=2, dual_mapping=False) \n",
    "u,v = fesL2.TnT()\n",
    "fesH1 = H1(mesh, order=2, definedon=mesh.Boundaries(\".*\"))\n",
    "uH1,vH1 = fesH1.TnT()\n",
    "print (\"ndofL2 = \", fesL2.ndof, \"ndof H1 = \", fesH1.ndof)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "Project the given Neumann data $u_1$ in finite element space and have a look at the boundary data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "uexa = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )\n",
    "\n",
    "graduexa = CF( (uexa.Diff(x), uexa.Diff(y), uexa.Diff(z)) )\n",
    "n = specialcf.normal(3)\n",
    "u1 = GridFunction(fesL2)\n",
    "u1.Interpolate(graduexa*n, definedon=mesh.Boundaries(\".*\"))\n",
    "\n",
    "Draw (u1, mesh, draw_vol=False, order=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "The discretisation of the variational formulation leads to a system of linear equations, ie \n",
    "\n",
    "$$ \\left( \\mathrm{D} + \\mathrm{S}\\right) \\, \\mathrm{u}_0 = \\left( \\frac12 \\,\\mathrm{M} - \\mathrm{K}' \\right) \\, \\mathrm{u}_1\\,, $$\n",
    "\n",
    "where the linear operators are as follows\n",
    "- $\\mathrm{D}$ is the hypersingular operator and the stabilisation $(\\mathrm D + \\mathrm{S})$ is regular and symmetric.\n",
    "- $\\mathrm{M}$ is a mass matrix.\n",
    "- $\\mathrm{K}'$ is the adjoint double layer operator. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "The stabilisation matrix $S$ is needed to cope with the kernel of the hypersingular operator. The following stabilisation matrix is used here: \n",
    "\n",
    "$$ S \\in \\mathbb R^{n\\times n}, \\quad S_{ij} = \\int\\limits_{\\Gamma} v_i(x) \\,\\mathrm{d} \\sigma \\cdot \\int\\limits_{\\Gamma} v_j(x) \\,\\mathrm{d} \\sigma\\,,$$ \n",
    "\n",
    "where $v_i, v_j$ being $H^{\\frac12}(\\Gamma)$-conforming shape functions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "vH1m1 = LinearForm(vH1*1*ds(bonus_intorder=3)).Assemble()\n",
    "S = (BaseMatrix(Matrix(vH1m1.vec.Reshape(1))))@(BaseMatrix(Matrix(vH1m1.vec.Reshape(fesH1.ndof))))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "Now we assemble the stabilised system matrix and the right hand side and solve for the Dirichlet data $u_0$ with an iterative solver:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "u0 = GridFunction(fesH1)\n",
    "pre = BilinearForm(uH1*vH1*ds).Assemble().mat.Inverse(freedofs=fesH1.FreeDofs()) \n",
    "with TaskManager():\n",
    "    D = HypersingularOperator(fesH1, intorder=12)\n",
    "    \n",
    "    M = BilinearForm( v.Trace() * uH1 * ds(bonus_intorder=3)).Assemble()\n",
    "    Kt = DoubleLayerPotentialOperator(fesL2, fesH1, intorder=12)    \n",
    "    rhs = ( (0.5 * M.mat.T - Kt.mat) * u1.vec ).Evaluate()\n",
    "\n",
    "    CG(mat=D.mat+S, pre=pre, rhs=rhs, sol=u0.vec, tol=1e-8, maxsteps=200, printrates=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "Let's have a look at the numerical and exact Dirichlet data and compare them quantitatively:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (u0, mesh, draw_vol=False, order=3);\n",
    "u0exa = GridFunction(fesH1)\n",
    "u0exa.Interpolate (uexa, definedon=mesh.Boundaries(\".*\"))\n",
    "Draw (u0exa, mesh, draw_vol=False, order=3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "print (\"L2 error of surface gradients =\", sqrt (Integrate ( (grad(u0exa)-grad(u0))**2, mesh, BND)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}