{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "11.2.3 Neumann Laplace Indirect Method\n",
    "======================================================\n",
    "**keys**: homogeneous Neumann bvp, hypersingular operator"
   ]
  },
  {
   "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 solution $u\\in H^1(\\Omega)$ of the above bvp can be written in the following form (representation forumula) \n",
    "\n",
    "$$ u(x) =  \\underbrace{ \\int\\limits_\\Gamma \\displaystyle{\\frac{1}{4\\,\\pi}\\, \\frac{ \\langle n(y), x-y \\rangle }{\\| x-y\\|^3} } \\, m(y)\\, \\mathrm{d}\\sigma_y}_{\\displaystyle{ \\mathrm{DL}(m) }}\\,.$$ \n",
    "\n",
    "Let's carefully apply the Neumann trace $\\gamma_1$ to the representation formula and solve the resulting boundary integral equation for $m \\in H^{\\frac12}(\\Gamma)$ by discretisation of the following variational formulation: \n",
    "\n",
    "$$ \\forall \\, v\\in H^{\\frac12}(\\Gamma): \\left\\langle v, \\gamma_1 \\left(\\mathrm{DL}(m)\\right) \\right\\rangle_{-\\frac12} = \\left\\langle u_1, v\\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 = Sphere( (0,0,0), 1)\n",
    "mesh = Mesh( OCCGeometry(sp).GenerateMesh(maxh=0.3)).Curve(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "Define the finite element space for $H^{\\frac12}(\\Gamma)$ for given geometry :  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fesH1 = H1(mesh, order=2, definedon=mesh.Boundaries(\".*\"))\n",
    "uH1,vH1 = fesH1.TnT()\n",
    "print (\"ndof H1 = \", fesH1.ndof)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "Define the right hand side with given Neumann data $u_1$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "uexa = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )\n",
    "graduexa = CF( (uexa.Diff(x), uexa.Diff(y), uexa.Diff(z)) )\n",
    "\n",
    "n = specialcf.normal(3)\n",
    "u1exa = graduexa*n\n",
    "rhs = LinearForm(u1exa*vH1.Trace()*ds(bonus_intorder=3)).Assemble()"
   ]
  },
  {
   "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{m}= \\mathrm{rhs}\\,, $$ \n",
    "\n",
    "where $\\mathrm{D}$ is the hypersingular operator and the stabilisation $(\\mathrm D + \\mathrm{S})$ is regular and symmetric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "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": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = GridFunction(fesH1)\n",
    "pre = BilinearForm(uH1*vH1*ds).Assemble().mat.Inverse(freedofs=fesH1.FreeDofs()) \n",
    "with TaskManager(): \n",
    "    D=HypersingularOperator(fesH1, intorder=12)\n",
    "    CG(mat = D.mat+S, pre=pre, rhs = rhs.vec, sol=m.vec, tol=1e-8, maxsteps=200, initialize=False, printrates=True)\n",
    "\n",
    "Draw (m, mesh, draw_vol=False, order=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "**Note:** the hypersingular operator is implemented as follows \n",
    "\n",
    "$$ D = \\langle v, \\gamma_1 \\mathrm{DL}(m) \\rangle_{-\\frac12} = \\frac{1}{4\\pi} \\int\\limits_\\Gamma\\int\\limits_\\Gamma \\frac{ \\langle \\mathrm{\\boldsymbol {curl}}_\\Gamma  \\,m(x), \\mathrm{\\boldsymbol{curl}}_\\Gamma\\, v(y) \\rangle}{\\|x-y\\|} \\, \\mathrm{d} \\sigma_{y} \\, \\mathrm{d} \\sigma_x $$\n",
    "\n",
    "Details for instance in [Numerische Näherungsverfahren für elliptische Randwertprobleme](https://link.springer.com/book/10.1007/978-3-322-80054-1), p.127, p.259 (1st edition)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "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
}